# Perturbation theory without power series:
iterative construction of non-analytic operator spectra

Matteo Smerlak  Max Planck Institute for Mathematics in the Sciences, Leipzig,
Germany

###### Abstract

It is well known that quantum-mechanical perturbation theory often give rise
to divergent series that require proper resummation. Here I discuss simple
ways in which these divergences can be avoided in the first place. Using the
elementary technique of relaxed fixed-point iteration, I obtain convergent
expressions for various challenging ground states wavefunctions, including
quartic, sextic and octic anharmonic oscillators, the hydrogenic Zeeman
problem, and the Herbst-Simon Hamiltonian (with finite energy but vanishing
Rayleigh-Schrödinger coefficients), all at arbitarily strong coupling. These
results challenge the notion that non-analytic functions of coupling constants
are intrinsically “non-perturbative”. A possible application to exact
diagonalization is briefly discussed.

perturbation theory, eigenvalue problem

## Introduction.—

Dyson was the first to observe that, since quantum electrodynamics with a
negative fine structure constant is unstable, the perturbative series in
powers of $\alpha$ must diverge in the physical sector with $\alpha>0$, no
matter how small its numerical value [1]. Unfortunately, such asymptotic
series with zero radius of convergence are the rule rather than the exception:
simple problems such as the one-dimensional quartic anharmonic oscillator or
the hydrogen atom in an external field already give rise to perturbative
series that diverge for any coupling strength [2]. Similar issues arise in
molecular [3] and nuclear [4] physics, quantum chemistry [5], and field theory
[6].

The usual strategy to deal with a divergent perturbative series is to attempt
to _resum_ it, _i.e._ to use its coefficients to construct a convergent
expression in ways other than simple summation [7, 8, 9]. Various procedures
have been developed for this purpose, including sequence transformations [10],
Padé approximants [11], Borel-Laplace resummation [12, 13, 14], and order-
dependent mappings [15]. Applying these procedures in practice can be
difficult. For starters, since the perturbative coefficients grow factorially
or super-factorially, high-precision (or exact) arithmetics is usually
required. Second, methods based on analytic continuations (_e.g._ Borel
resummation or order-dependent mappings) require _a priori_ knowledge of the
singularity structure of the eigenvalues, which must obtained through
semiclassical estimates [7] or some other technique [16], and can give rise to
ambiguities [14].

Worse still, resummation is not always possible. Herbst and Simon exhibited an
anharmonic oscillator with finite ground state energy $E(g)\sim e^{-d/g^{2}}$
(with $d>0$) but vanishing perturbative coefficients at all orders [17],
meaning that perturbation theory contains _no_ information about the ground
state. To make progress, Jentschura and Zinn-Justin leveraged semiclassical
estimates to conjecture a modified Born-Sommerfeld quantization condition [18,
19] which gives the energy as a resurgent expansion [20]. Combined with Borel-
Padé resummation, this approach allows to compute $E(g)$ at small coupling
values $g$ [21]. However—like other resummation procedures—the method does not
by itself allow to construct the corresponding eigenvector.

The basic problem underlying the divergence of perturbative expansions is that
a function $E(g)$ of a variable $g$ with a singularity at
$g_{*}\in\mathbb{C}$— _e.g._ “intruder state” [22]—cannot be expressed as a
convergent powers series in $g$ outside a disk of convergence of radius
$|g_{*}|$ (Abel’s lemma); when such a singularity lies at $g=0$, a power
series expansions can never converge. Resummation theory is an effort to
bypass this obstruction by interpreting the perturbation series as the Taylor
expansion of some member of a more general, yet sufficiently rigid, class of
function of $g$: a rational function for Padé approximants, a resurgent
function for Borel-Ecalle resummation, etc. A common view is that non-analytic
contributions to the energy function $E(g)$ (_e.g._ instantons $\sim
e^{-d/g^{2}}$), being intrinsically “non-perturbative”, can only be captured
from the large order behavior of perturbation theory, if at all.

But there is more to perturbation theory than Taylor expansions in the
theory’s coupling constant $g$. First, given a Hamiltonian $H=F+gI$ (where $F$
stands for “free” and $I$ for “interaction”), nothing dictates that we treat
$F$ as the unperturbed Hamiltonian and $gI$ as the perturbation. When we use
any another partitioning into $H=H_{0}+\lambda H_{1}$ (with $H_{0}$ a spectral
function of $F$, with the same eigenvectors), the Rayleigh-Schrödinger (RS)
algorithm yields a solution which is a power series in $\lambda$, but not in
$g$. Second, whether we use $F$ or some other $H_{0}$ as unperturbed
Hamiltonian, assuming a power series ansatz is not necessarily the most
straightforward approach to solving the perturbation equation. I discuss below
an alternative formulation of perturbation theory, called “iterative
perturbation theory” (IPT), that does not rely on this ansatz and is more
computationally efficient.

These observations are elementary, but also empowering. In the following I
show that convergent approximations of the ground state energy—and
eigenvector—of even anharmonic oscillators, the Zeeman problem and the Herbst-
Simon Hamiltonian can be constructed perturbatively, at any coupling strength
and without the need for educated guesses, extraneous semiclassical estimates,
ingenious resummation techniques, or knowledge of the nature and location of
singular points in the complex $g$ plane.

Many authors have noted that the convergence of perturbation theory can be
improved by using the freedom to repartition the Hamiltonian [23, 24, 15, 25,
26, 27]. However, the large-order behavior of these optimized schemes has
rarely been analyzed. To my knowledge, the resummation-free construction of a
non-trivial quantum-mechanical ground state at arbitrarily high coupling
strength with these methods has not been reported (with the exception of a
zero-dimensional model [28]). More importantly—as I hope to convince the
reader—none of these optimized perturbation methods is as elementary as the
relaxed iterative scheme proposed here.

Figure 1: Ground state energy of the quartic oscillator computed with IPT
($\alpha=1/2$) under EN partitioning (7) (continuous lines) compared to exact
diagonalization in a basis set with $5\cdot 10^{4}$ harmonic oscillator states
(dashed lines), for increasingly large coupling constants $g$. Eigenvectors
and eigenvalues both converge (albeit slowly) without the need for any
resummation procedure.

## Eigenvalue perturbation theory.—

Consider a Hamiltonian of the form $H=H_{0}+\lambda H_{1}$ where $H_{0}$ has a
known spectrum and $H_{1}$ is a perturbation. (The method also applies to non-
Hermitian problems, but we restrict to Hamiltonian operators here for
simplicity). Here $\lambda$ is a perturbation parameter which we will
eventually set to $1$. We focus on an isolated, non-degenerate, stable
111Stability here means that the eigenvalues remains isolated and non-
degenerate for small values of $\lambda$. eigenvalue $E_{0}$ with
corresponding eigenvector $\psi_{0}$ and eigenprojection $P_{0}$. We seek an
eigenvector $\psi$ of $H$ with eigenvalue $E$ and
$\langle\psi_{0}|\psi\rangle=1$. Under this normalization condition, the
eigenvalue equation $H\psi=E\psi$ is equivalent to

$\quad\psi=\psi_{0}+R_{0}(\mathcal{E})(\lambda H_{1}+\mathcal{E}-E)\psi,$ (1)

where $R_{0}(\mathcal{E})=(I-P_{0})(\mathcal{E}-H_{0})^{-1}(I-P_{0})$ is the
reduced resolvent. In the following we focus on the choice
$\mathcal{E}=E_{0}=E-\langle\psi_{0}|\lambda H_{1}\psi\rangle$, which turns
the perturbation equation into the quadratic equation 222Choosing instead
$\mathcal{E}=E$ is the starting point of Brillouin-Wigner perturbation theory,
in which $E$ is computed implicitly rather than explicitly.

$\psi=\psi_{0}+\lambda
R_{0}(E_{0})(H_{1}-\langle\psi_{0}|H_{1}\psi\rangle)\psi.$ (2)

In a broad sense, perturbation theory consists in finding solution to (2)
iteratively. There is more than one way to approach this problem.

## RS perturbation theory.—

Conventional RS theory seeks solutions of (2) as a power series
$\psi_{\textrm{RS}}=\sum_{l\geq 0}\lambda^{l}a_{l}$. Inserting this ansatz
into (2) and collecting terms of the same order in $\lambda$ gives the well-
known recursion relation $a_{l+1}=Q_{\textrm{RS}}(a_{0},\cdots,a_{l})$ with

$Q_{\textrm{RS}}(a_{0},\cdots,a_{l})\equiv\lambda
R_{0}(E_{0})[H_{1}a_{l}-\sum_{s=0}^{l}\langle\psi_{0}|H_{1}a_{s}\rangle
a_{l-s}]$ (3)

with initial condition $a_{0}=\psi_{0}$. In many cases of physical interest
(including all examples below, and others [31]), the perturbation $H_{1}$ is
block-diagonal in the basis of unperturbed eigenstates, hence (3) is
effectively a final-dimensional problem which can be solved at any finite
order $k$ using computer algebra 333In particular, to any finite perturbative
order $k$, all sums over unperturbed eigenstates are all finite sums.. (When
this condition is not met, the perturbation equations must usually be solved
numerically in a finite basis set, _i.e._ perturbation theory becomes a
variational approximation.) From these eigenvector coefficients the energy can
be computed at any order $k$, either directly as
$E^{(k)}_{\textrm{RS}}=\sum_{l\geq
0}^{k-1}\lambda^{l}\langle\psi_{0}|Ha_{l}\rangle$, or more efficiently through
Wigner’s $2n+1$ theorem.

Figure 2: Convergence of IPT energy iterates $E^{(k)}$ for the ground states
of the quartic, sextic and octic anharmonic oscillators, all with $g=100$.
Here we used $\alpha=1/2$ and EN partitioning; $E$ is obtained by exact
diagonalization. Note the doubly-logarithmic axes: convergence is slow for
these examples.

## Iterative perturbation theory.—

An alternative approach, termed here “iterative perturbation theory” (IPT),
starts from the observation that (2) is a fixed-point equation
$\psi=Q_{\textrm{IPT}}(\psi)$ for the non-linear operator

$Q_{\textrm{IPT}}(\psi)\equiv\psi_{0}+\lambda
R_{0}(E_{0})(H_{1}-\langle\psi_{0}|H_{1}\psi\rangle)\psi.$ (4)

This suggests that $\psi$ can also be computed dynamically, as the limit of
the sequence of iterates
$\psi^{(k+1)}_{\textrm{IPT}}=Q_{\textrm{IPT}}(\psi^{(k)})$ starting from
$\psi^{(0)}=\psi_{0}$. The corresponding energy given by
$E^{(k)}_{\textrm{IPT}}=\langle\psi_{0}|H\psi^{(k-1)}_{\textrm{IPT}}\rangle$.
This approach is _a priori_ more efficient than the RS recursion (3): IPT is a
first-order recurrence relation, whereas (3) involves all previous $a_{l}$
coefficients. Moreover, IPT makes it easy to prove convergence results using
Banach’s fixed point theorem. For instance, we can show that, if $H_{1}$ is
bounded, $\psi^{(k)}$ converges as $k\to\infty$ for any
$|\lambda|<(3-2\sqrt{2})/\|H_{1}\|\delta$ where $\delta$ the distance between
$E_{0}$ and the rest of $H_{0}$’s spectrum [33]. The corresponding result in
RS theory is “not at all trivial” _dixit_ Kato [34]. Finally, IPT is directly
related to RS as follows: after $k$ iterations, $\psi^{(k)}_{\textrm{IPT}}$
contains powers of $\lambda$ with exponents up to $2^{k}-1$, but low-order
terms coincide with RS, _viz._
$\psi^{(k)}_{\textrm{IPT}}=\psi^{(k)}_{\textrm{RS}}+\mathcal{O}(\lambda^{k+1})$,
as can be easily seen by recursion.

## Relaxation.—

Both RS and IPT involve recurrence relations of the form $X_{n+1}=Q(X_{n})$.
The main purpose of this paper is to direct attention to the fact that the
convergence of perturbative calculations is much improved if we use instead
the _relaxed_ iteration procedure

$X_{n+1}=\alpha Q(X_{n})+(1-\alpha)X_{n}\ \textrm{for some}\ 0<\alpha<1.$ (5)

Relaxation is a well-known technique in numerical analysis: it simply amounts
to choosing a smaller time step in the Euler discretization of a dynamical
system—a natural strategy to improve convergence. From the perspective of
perturbation theory, (5) turns out to be equivalent to a repartitioning of the
Hamiltonian into

$H=[\alpha^{-1}H_{0}]+[H_{1}+(1-\alpha^{-1})H_{0}].$ (6)

Eq. (6) with a special choise of $\alpha$ is known in quantum chemistry as
“Feenberg scaling” [35, 36, 37]. Such a repartitioning is also the starting
point of order-dependent mappings [15] and of self-consistent expansions [38].
They key point is that, upon such a repartitioning, neither RS nor IPT give
rise to power series of $\lambda$ 444IPT never does: even with $\alpha=1$,
$\psi^{(k)}_{\textrm{IPT}}$ is a polynomial in $\lambda$ with $k$-dependent
coefficients., but instead to expansions of the form
$\psi^{(k)}=\sum_{l=0}^{c(k)}b_{k,l}(\alpha)\lambda^{l}$ with $c(k)=k$ for RS
theory and $c(k)=2^{k}-1$ for IPT. Abel’s lemma does not apply to such
expansions: $\psi=\lim_{k\to\infty}\psi^{(k)}$ can admit a singularity at some
$\lambda=\lambda_{*}$, and nevertheless converge outside the disk with radius
$|\lambda_{*}|$. For an illustration of this fact using a simple two-
dimensional example see Supplementary Material I.

## Epstein-Nesbet partitioning.—

We have noted before that, given $H=F+gI$ where $F$ is a free theory (say a
harmonic oscillator or a Fockian) and $g$ a physical coupling constant,
setting $H_{0}=F$ and $H_{1}=gI$ is not necessarily the most judicious choice.
Indeed another natural choice, Epstein-Nesbet (EN) partitioning [40], is often
preferable. In the basis of eigenvectors of $F$, this is defined by

$H_{0}=\textrm{diag}(H_{nn})_{n\geq 0}=F+g\,\textrm{diag}(I_{nn})_{n\geq 0}$
(7)

and $H_{1}=H-H_{0}$. Perhaps because it is incompatible with the notion that
perturbation theory should provide an expansion in powers of $g$, this choice
is rarely used in the physical literature.

We now proceed to illustrate the virtues of combining EN partitioning with
relaxed iteration through several examples. Although both RS and IPT give
comparable results, in the following I focus on the latter (omitting
subscripts) owing to its simpler structure and greater computational
efficiency.

## Even anharmonic oscillators.—

Even-order one-dimensional anharmonic oscillor $H=p^{2}+x^{2}+gx^{2s}$ are
commonly cited to illustrate the difficulties encountered in RS theory [2, 9].
It was noted in [41] that, using (7) but no relaxation, the RS series for the
quartic oscillator $s=2$ may converge for very small values of $g$. We checked
that this is no longer true for $s=3,4$.

Relaxed perturbation theory with EN partitioning, on the other hand, always
converges. Fig. 1 shows IPT energy iterates for the quartic oscillator at
increasingly large coupling $g$; Fig. 2 in turn shows the relative errors for
the quartic, sextic and octic oscillators with $g=100$. Without restrictions,
we find that $E^{(k)}$ converges to the exact eigenvalue $E$ as $k\to\infty$,
albeit increasingly slowly as $s$ or $g$ gets larger. We emphasize that the
convergence of $E^{(k)}$ is underpinned by the convergence of the eigenvector
$\psi^{(k)}$ itself, as can be checked by computing the residuals
$\|H\psi^{(k)}-\langle\psi_{0}|H\psi^{(k)}\rangle\psi^{(k)}\|$ (not shown).

## Herbst-Simon Hamiltonian.—

The Herbst-Simon anharmonic oscillator with potential
$V(x;g)=2gx-2gx^{3}+g^{2}x^{4}$ has a purely instantonic positive ground state
energy $E\sim e^{-d/g^{2}}$ for some $d>0$. This unusual non-analytic
behaviour precludes a meaningful expansion in power of $g$, and shows that the
standard RS expansion (identically zero in this case) can sometimes converge
to the wrong value. As already noted, successful perturbative calculations of
$E$ for g up to $\sqrt{0.3}$ have so far involved educated guesses and
sophisticated resummation procedures [18, 21].

Figure 3: The Herbst-Simon ground state energy vanishes to all orders in $g$,
but can be accurately computed using relatex perturbation theory, here IPT
with $\alpha=1/2$ and EN partitioning (7).

Using IPT with $\alpha=1/2$ (again in EN partitioning) we easily compute $E$
for any $g$, to arbitrary accuracy. For instance, to obtain
$E(\sqrt{0.3})\simeq 1.11$ to two significant figures (the state of art using
the methods of Ref. [21]), we only need a dozen iterations of IPT. With
$k=220$ iterations, we obtain all $13$ figures of the exact result
$[E(\sqrt{0.3})-1]/2=5.318357438655\cdot 10^{-2}$ cited in [21, Table III].
Fig. 3 illustrates the convergence of $E^{(k)}$ to the exact value
(indistinguishable from the red curve on this figure), with its non-polynomial
behavior near $g=0$.

## Zeeman problem.—

We close with an example of historical as well as experimental significance,
the hydrogen atom in a magnetic field, $H=-\Delta/2-1/r+g(x^{2}+y^{2})$ with
$g=B^{2}/8$ (in atomic units). To compute Zeeman shifts it is useful to make
use of the $SO(4,2)$ dynamical symmetry of the problem [42, 31]. Thus, for the
ground state $\psi$ with energy $E=-1/2+\Delta E$, the Schrödinger equation
may be reformulated as the generalized eigenvalue problem

$(T_{3}-1+gW)\psi=(\Delta E)S\psi$ (8)

where $T_{3}$ is a generators of a $\mathfrak{so}(2,1)$ sub-algebra (with
discrete spectrum), $W$ and $S$ can be expressed in terms of the other
generators of $\mathfrak{so}(4,2)$ [31, Chap. 9]. Using this representation
with the partitioning $H_{0}=T_{3}-I$, it is possible to compute the RS
coefficients of $\psi$ to large orders [3]. The resulting series diverges for
any $g\neq 0$, but Padé approximants provide accurate estimates of $E$ up to
$g\simeq 1$. In stronger fields, Borel resummation combined with order
dependent mappings [43] or sequence transformations [44] give better results.
These methods are all rather sophisticated and generally use additional
information about the $g\to\infty$ limit, _e.g._ via a scaling of space which
maps $g$ to the unit interval 555The RS series for the Zeeman has been called
“one of the most difficult summable divergent series encountered in physics”
[44]..

A straightforward generalization of IPT to account for the $S$ matrix on the
RHS of (8) reproduces these results without using _ad hoc_ information , and
without any limitation on $g$ (other than the slow convergence in very strong
fields). For $B=1$ (a strong field of $2\cdot 10^{5}$ T) we obtain the correct
result $E=-0.3312$ with just $k=100$ iterations of IPT ($\alpha=0.3$). For
$B=10$, $k=1000$ iterations give $E^{(k)}=3.72$, to be compared with $E=3.25$
obtained in Ref. [44] using Weniger-type sequence transformations. Using
elementary Aitken extrapolation this improves to $E^{(k)^{\prime}}=3.35$
666Aitken extrapolation of a sequence $s_{n}$ is
$s_{n}^{\prime}=(s_{n}s_{n+2}-s_{n+1}^{2})/(s_{n}+s_{n+2}-2s_{n+1})$..

Disorder $h$ | IPR | Iterations / CPU time ($s$)
---|---|---
|  | IPT-AA | KS | LOBPCG
1 | 0.004 | 396/19.4 | 8/2.1 | 85/4.4
2 | 0.11 | 64/3.2 | 7/1.9 | 65/3.4
5 | 0.68 | 21/1.0 | 7/1.8 | 62/3.3
10 | 0.92 | 14/0.64 | 5/1.3 | 45/2.4
50 | 0.99 | 7/0.30 | 4/1.1 | 32/1.7
100 | 0.999 | 5/0.22 | 3/0.93 | 27/1.4
Table 1: Exact (numerical) diagonalization of the ground state of the random-
field Heisenberg spin chain. For disorder strengths $h\geq 3.7$ this system is
many-body localized, implying that eigenstates are close to states in the
computional basis [serbyn2013local] and display high inverse participation
ratio (IPR). We used $L=20$ spins, periodic boundary conditions and set the
tolerance to $10^{-10}$. The Hamiltonian matrix was built using quimb [47],
and Krylov-Schur (KS) and LOBPCG calculations were done with SLEPc [48]. For
LOBPCG we used a trivial preconditioner, which proved more efficient than a
diagonal preconditioner in this case. IPT was applied with Anderson
acceleration (IPT-AA) with memory $M=10$ [49]. Timings are with a $3.6$ GHz
10-core Intel i9 processor.

## IPT as eigensolver.—

Like conventional RS theory, IPT is a scheme for computing analytical
approximations of perturbed eigenvectors. However, in cases where $H_{1}$ is
not block diagonal or can only be obtained in a finite basis set, as is common
in quantum chemistry, the latter can also be used as an efficient numerical
eigensolver. Given a matrix $H$, the simplest numerical implementation of IPT
consists in replicating the relaxed perturbative calculations above: set
$H_{0}=\textrm{diag}(H)$, start from a basis vector, say
$\psi_{0}=(1,0,\cdots,0)$, and iterate $\psi^{(k+1)}=\alpha
Q(\psi^{(k)})+(1-\alpha)\psi^{(k)}$ until a convergence condition is met.
However, instead of using a fixed $\alpha$, it is more efficient to use an
optimal convex combination of $Q_{\textrm{IPT}}(\psi_{k-1})$ and previous
iterates. This method, known as Anderson acceleration [50], constructs the new
iterate as

$\psi^{(k+1)}=\sum_{m=0}^{M}\beta_{m}Q(\psi^{(k-m)}),$ (9)

where $M$ is a memory parameter and the positive coefficients $\beta_{m}$’s
are recomputed at each step by minimizing the residual norm.

With this modification, IPT provides an eigensolver that can be remarkably
efficient: a python implementation given in Supplementary Material II proves
up to $5$x faster that state-of-the-art implementations of the Implicitly
Restarted Lanczos and Locally Optimal Block Preconditioned Conjugate Gradient
(LOBPCG) algorithms 777IPT can also be applied to non-Hermitian problems, in
which case speed-up factors reach $100$ or more [33].. (Other recent methods,
such as Generalized Davidson or Jacobi-Davidson, were slower in this case.)
Table 1 reports the number of iterations and CPU time required to compute the
ground state of a random-field Heisenberg spin chain at various disorder
strengths $h$, showing significant speed-ups in the many-body localized phase
$h\geq 3.7$ (for notation and motivation see e.g. [52]). It is important to
emphasize, however, that unlike Lanczos and other general-purpose eigenvalue
algorithms, IPT has a limited applicability: due to its perturbative nature,
the performance of IPT degrades rapidly when off-diagonal elements become
large compared to the separation between its diagonal elements. For more
details on the numerical aspects of IPT, I refer the reader to Ref. [33].

## Conclusion.—

“Relaxed perturbation theory” is the idea of applying relaxed iteration to
compute the eigenvectors of perturbed operators, either using the well-known
RS scheme, or, more conveniently, through the application of the non-linear
operator $Q_{\textrm{IPT}}$. Combined with the Epstein-Nesbet partitioning
prescription, this technique provides convergent approximations for the ground
state wavefunctions (not just energies) of challenging Hamiltonians, including
the purely instantonic Herbst-Simon Hamiltonian and the hydrogenic Zeeman
problem, up to arbitrarily high field coupling strengths. My interpretation of
these results is that “non-perturbative physics” is, in fact, squarely within
the scope of perturbation theory.

Future work should aim to put these results on stronger mathematical footing.
We noted earlier that sufficient conditions for the convergence of IPT are
easily derived in the special case where $H_{1}$ is bounded, a restrictive
assumption which does not cover the physical examples discussed in this paper.
Numerical evidence shows that these conditions are far too strong, but
mathematical progress towards relaxing them has so far remained beyond my
reach.

###### Acknowledgements.

It is a pleasure to thank Maxim Kenmoe for stimulating my interest in
perturbation theory, Anton Zadorin for a productive collaboration on the
mathematical aspects of IPT, and an anonymous referee for constructive
comments. Funding for this work was provided by the Alexander von Humboldt
Foundation in the framework of the Sofja Kovalevskaja Award endowed by the
German Federal Ministry of Education and Research.

## References

  * Dyson [1952] F. J. Dyson, Phys. Rev. 85, 631 (1952).
  * Simon [1982] B. Simon, Int. J. Quantum Chem. 21, 3 (1982).
  * Čížek and Vrscay [1982] J. Čížek and E. R. Vrscay, International Journal of Quantum Chemistry 21, 27 (1982).
  * Roth and Langhammer [2010] R. Roth and J. Langhammer, Phys. Lett. B 683, 272 (2010).
  * Olsen  _et al._ [1996] J. Olsen, O. Christiansen, H. Koch, and P. Jo/rgensen, J. Chem. Phys. 105, 5082 (1996).
  * Rivasseau [1991] V. Rivasseau,  _From Perturbative to Constructive Renormalization_ , 2nd ed. (Princeton University Press, Princeton, N.J, 1991).
  * Guillou and Zinn-Justin [1990] J. C. L. Guillou and J. Zinn-Justin, eds.,  _Large-Order Behaviour of Perturbation Theory_ (North Holland, Amsterdam ; New York, 1990).
  * Janke and Kleinert [2007] W. Janke and H. Kleinert,  _Resummation of Divergent Perturbation Series: Introduction to Theory and Guide to Practical Applications_ (World Scientific Pub Co Inc, s.l., 2007).
  * Caliceti  _et al._ [2007] E. Caliceti, M. Meyer-Hermann, P. Ribeca, A. Surzhykov, and U. D. Jentschura, Physics Reports 446, 1 (2007).
  * Weniger [1989] E. J. Weniger, Comput. Phys. Rep. 10, 189 (1989).
  * Loeffel  _et al._ [1969] J. J. Loeffel, A. Martin, B. Simon, and A. S. Wightman, Phys. Lett. B 30, 656 (1969).
  * Graffi  _et al._ [1970] S. Graffi, V. Grecchi, and B. Simon, Phys. Lett. B 32, 631 (1970).
  * Mera  _et al._ [2015] H. Mera, T. G. Pedersen, and B. K. Nikolić, Phys. Rev. Lett. 115, 143001 (2015).
  * Delabaere  _et al._ [2016] E. Delabaere, M. Loday-Richaud, C. Mitschi, and D. Sauzin,  _Divergent Series, Summability and Resurgence I-III_ , Lecture Notes in Mathematics (Springer International Publishing, 2016).
  * Seznec and Zinn-Justin [1979] R. Seznec and J. Zinn-Justin, J. Math. Phys. 20, 1398 (1979).
  * Kvaal  _et al._ [2011] S. Kvaal, E. Jarlebring, and W. Michiels, Phys. Rev. A 83, 032505 (2011).
  * Herbst and Simon [1978] I. W. Herbst and B. Simon, Phys. Lett. B 78, 304 (1978).
  * Zinn-Justin and Jentschura [2004] J. Zinn-Justin and U. D. Jentschura, Ann. Phys. (N. Y.) 313, 197 (2004).
  * Jentschura and Zinn-Justin [2004] U. D. Jentschura and J. Zinn-Justin, Phys. Lett. B 596, 138 (2004).
  * Aniceto  _et al._ [2019] I. Aniceto, G. Başar, and R. Schiappa, Phys. Rep. 809, 1 (2019).
  * Surzhykov  _et al._ [2006] A. Surzhykov, M. Lubasch, J. Zinn-Justin, and U. D. Jentschura, Phys. Rev. B 74, 205317 (2006).
  * Choe  _et al._ [2001] Y.-K. Choe, H. A. Witek, J. P. Finley, and K. Hirao, J. Chem. Phys. 114, 3913 (2001).
  * Hirschfelder [1963] J. O. Hirschfelder, J. Chem. Phys. 39, 2099 (1963).
  * Caswell [1979] W. E. Caswell, Annals of Physics 123, 153 (1979).
  * Shirokov and Blinov [1995] I. Shirokov and N. Blinov, J. Exp. Theor. Phys. 80, 3 (1995).
  * Scherer [1995] W. Scherer, Phys. Rev. Lett. 74, 1495 (1995).
  * Szabados and Surján [1999] Á. Szabados and P. R. Surján, Chem. Phys. Lett. 308, 303 (1999).
  * Remez and Goldstein [2018] B. Remez and M. Goldstein, Phys. Rev. D 98, 056017 (2018).
  * Note [1] Stability here means that the eigenvalues remains isolated and non-degenerate for small values of $\lambda$.
  * Note [2] Choosing instead $\mathcal{E}=E$ is the starting point of Brillouin-Wigner perturbation theory, in which $E$ is computed implicitly rather than explicitly.
  * Adams [1994] B. G. Adams,  _Algebraic Approach to Simple Quantum Systems: With Applications to Perturbation Theory_ (Springer-Verlag, Berlin Heidelberg, 1994).
  * Note [3] In particular, to any finite perturbative order $k$, all sums over unperturbed eigenstates are all finite sums.
  * Kenmoe  _et al._ [2021] M. Kenmoe, R. Kriemann, M. Smerlak, and A. S. Zadorin, arXiv:2012.14702 [math.NA] (2021), arXiv:2012.14702 [math.NA] .
  * Kato [2013] T. Kato,  _Perturbation Theory for Linear Operators_ (Springer, 2013).
  * Feenberg [1958] E. Feenberg, Annals of Physics 3, 292 (1958).
  * Dietz  _et al._ [1993] K. Dietz, C. Schmidt, M. Warken, and B. A. Hess, J. Phys. B: At. Mol. Opt. Phys. 26, 1885 (1993).
  * Schmidt  _et al._ [1993] C. Schmidt, M. Warken, and N. C. Handy, Chemical Physics Letters 211, 272 (1993).
  * Schwartz and Edwards [1992] M. Schwartz and S. F. Edwards, EPL 20, 301 (1992).
  * Note [4] IPT never does: even with $\alpha=1$, $\psi^{(k)}_{\textrm{IPT}}$ is a polynomial in $\lambda$ with $k$-dependent coefficients.
  * Surján and Szabados [2004a] P. R. Surján and Á. Szabados, in  _Fundamental World of Quantum Chemistry: A Tribute to the Memory of Per-Olov Löwdin Volume III_, edited by E. J. Brändas and E. S. Kryachko (Springer Netherlands, Dordrecht, 2004) pp. 129–185.
  * Surján and Szabados [2004b] P. R. Surján and Á. Szabados, Collect. Czech. Chem. Commun. 69, 105 (2004b).
  * Čížek and Paldus [1977] J. Čížek and J. Paldus, International Journal of Quantum Chemistry 12, 875 (1977).
  * Le guillou and Zinn-justin [1990] J. C. Le guillou and J. Zinn-justin, in  _Current Physics–Sources and Comments_, Large-Order Behaviour of Perturbation Theory, Vol. 7, edited by J. C. Le guillou and J. Zinn-justin (Elsevier, 1990) pp. 259–286.
  * Čížek  _et al._ [2003] J. Čížek, J. Zamastil, and L. Skála, J. Math. Phys. 44, 962 (2003).
  * Note [5] The RS series for the Zeeman has been called “one of the most difficult summable divergent series encountered in physics” [44].
  * Note [6] Aitken extrapolation of a sequence $s_{n}$ is $s_{n}^{\prime}=(s_{n}s_{n+2}-s_{n+1}^{2})/(s_{n}+s_{n+2}-2s_{n+1})$.
  * Gray [2018] J. Gray, Journal of Open Source Software 3, 819 (2018).
  * Hernandez  _et al._ [2005] V. Hernandez, J. E. Roman, and V. Vidal, ACM Trans. Math. Softw. 31, 351 (2005).
  * Zhang  _et al._ [2020] J. Zhang, B. O’Donoghue, and S. Boyd, SIAM J. Optim. 30, 3170 (2020).
  * Walker and Ni [2011] H. F. Walker and P. Ni, SIAM J. Numer. Anal. 49, 1715 (2011).
  * Note [7] IPT can also be applied to non-Hermitian problems, in which case speed-up factors reach $100$ or more [33].
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 91, 081103 (2015).


# Polytope volume in Normaliz

Winfried Bruns  Universität Osnabrück, Institut für Mathematik, 49069
Osnabrück, Germany<EMAIL_ADDRESS>To the memory of Wolmer Vasconcelos
and to Rafael Villareal on his 70th birthday

###### Abstract.

We survey the computation of polytope volumes by the algorithms of Normaliz to
which the Lawrence algorithm has recently been added. It has enabled us to
master volume computations for polytopes from social choice in dimension
$119$. This challenge required a sophisticated implementation of the Lawrence
algorithm.

###### Key words and phrases:

polytope, volume, Lawrence algorithm

###### 2010 Mathematics Subject Classification:

52B55, 52A38

##  1\. Introduction

About 20 years ago Amelia Taylor asked the author whether Normaliz [10] could
compute polytope volumes. It was easy to include this computation goal into
the triangulation based “primal algorithm”. Since then, polytope volumes have
played an important role in the development of Normaliz, and in recent years
specific algorithms have been added.

Polytope volumes can be interpreted as degrees of projective toric varieties
and multiplicities of monomial algebras. In 2011, Bogdan Ichim pointed out
their applications in social choice. Since then they have been a driving
challenge for the volume algorithms in Normaliz whose history we sketch
briefly. Before going on, let us emphasize that Normaliz computes lattice
normalized volumes that for rational polytopes are rational numbers. Euclidean
volumes, if asked for by the user, are derived from them.

In social choice, polytope volumes are interpreted as probabilities of certain
paradoxa and quality measures of voting schemes. See the books by Gehrlein and
Lepelley [16], [17]. These applications become rapidly very difficult since
the relevant polytopes explode in dimension: for $n$ candidates they are
cross-sections of cones of dimension $n!$. In their paper [19, p. 382] of 2008
Lepelley, Louichi and Smaoui state:

> Consequently, it is not possible to analyze four candidate elections, where
> the total number of variables (possible preference rankings) is $24$. We
> hope that further developments of these algorithms will enable the
> overcoming of this difficulty.

With the efficient parallelization of Normaliz in 2012 and the addition of
Schürmann’s symmetrization method [23], Normaliz could compute a wide variety
of social choice polytopes for $4$ candidates; see Bruns, Ichim and Söger [9].

When the limitations of the available algorithms became visible in 2017, the
author implemented an algorithm for polytope volumes by descent in the face
lattice. It is implicitly based on a reverse lexicographic triangulation of
the polytope, but does not compute the triangulation explicitly. It brought a
significant improvement in computation times for polytopes defined by
inequalities, and made more computations for elections with $4$ candidates
possible, as shown in Bruns and Ichim [6].

But the case $n=5$, $n!=120$, remained elusive. The breakthrough came with the
Normaliz implementation of Lawrence’s algorithm [18]. It is based on a duality
between “generic” triangulations of the dual cone and signed decompositions of
the “primal” polytope into simplices. In principle, signed decompositions are
as good as ordinary ones for volume computation, but they present hard
numerical problems. The rational arithmetic of Normaliz can cope with them,
but it must pay by computation time. The applications to $5$ candidates
elections have been documented by Bruns and Ichim [7]. They would have been
unreachable without the sophisticated implementation that we explain in
Section 4.

Simultaneously with the Lawrence algorithm, we introduced a refinement of the
descent algorithm that identifies isomorphic faces in the descent. Isomorphism
classes, as explained in [4], are computed by McKay and Piperno’s package
nauty [21]. Even for nauty, isomorphism classes are expensive, but they help
in the volume computation of some classical polytopes.

Section 2 gives a very brief introduction to the terminology of this note. It
explains basic results that are used in the computation of lattice normalized
volume. Section 3 contains an overview of the volume algorithms and explains
them, in particular Lawrence’s algorithm. The final section 5 lists
computation times, with emphasis on the new algorithms, signed decomposition
and descent with the exploitation of isomorphism types. They are not only
applied to polytopes from social choice, but also to classical polytopes, for
example cubes, Birkhoff polytopes and linear ordering polytopes.

The package vinci [12] contains algorithms for polytope volumes. See [6] for a
comparative study of computation times and memory usages. Because of its
floating point arithmetic, vinci is often faster than Normaliz, but its
results come without an error bound. Because of the extreme numerical
difficulty, its implementation of the Lawrence algorithm fails reliable
results already for polytopes coming from $4$ candidates elections; see Remark
6.

This note is dedicated to my friends Wolmer Vasconcelos and Rafael Villarreal.
Their constant support has been very encouraging in the $25$ years of the
Normaliz project. One of the first third party publications citing Normaliz is
their paper [14] with Delfino, Taylor and Weininger. The example collection of
Normaliz still contains input files supplied by Rafael a quarter of a century
ago, and his book [24] documents numerous applications.

_Acknowledgement._ The author was partially supported by the DFG grant BR
688/26-1. He thanks Ulrich von der Ohe for fruitful discussions.

##  2\. Preliminaries

We refer the reader to [5] for discrete convex geometry. Here we content
ourselves to a very brief overview.

###  2.1. Cones and polytopes

A _cone_ $C$ in the real space ${\mathbb{R}}^{n}$ is the intersection of
finitely many linear halfspaces:

$C=\bigcap_{i=1}^{s}H_{i}^{+},$ (2.1)

and for each $i$ the halfspace $H_{i}^{+}$ is the set
$\\{x\in{\mathbb{R}}^{n}:\sigma_{i}(x)\geq 0\\}$ for a linear form
$\sigma_{i}$ in the dual space $({\mathbb{R}}^{n})^{*}$. By the theorem of
Minkowski-Weyl, one can equivalently describe cones as the conical set
generated by finitely many vectors $v_{j}\in{\mathbb{R}}^{n}$,

$C=\\{q_{1}v_{1}+\dots+q_{n}v_{m}:q_{1},\dots,q_{m}\geq 0\\}.$ (2.2)

Since we want to deal only with polytopes and cones derived from them, we can
restrict our cones to a subclass: $C$ is a _pointed_ cone: if $-x\in C$ for
$x\in C$, then $x=0$. If $C$ is pointed, then the elements in a minimal set of
generators as in (2.2) are uniquely determined up to positive scalars, and the
sets ${\mathbb{R}}_{+}v_{i}$ are the _extreme rays_ of $C$.

The cone $C$ is _rational_ if the vectors $v_{i}$ can be chosen in
${\mathbb{Q}}^{n}$, and therefore in ${\mathbb{Z}}^{n}$. Then each extreme ray
contains exactly one _primitive_ integral vector, namely one with coprime
coordinates. It is called an _extreme integral generator_.

The dimension of $C$ is the dimension of the vector subspace
${\mathbb{R}}_{+}C$. If $\dim C<n$, then the halfspaces $H_{i}^{+}$ in (2.1)
are not uniquely determined, but the halfspaces $H_{i}^{+}\cap{\mathbb{R}}C$
of ${\mathbb{R}}C$ in an irredundant representation
$C=\bigcap_{i}(H_{i}^{+}\cap{\mathbb{R}}C)$ are. They intersect $C$ in its
_facets_. More generally, a _face_ of $C$ is the intersection of $C$ with a
hyperplane that has $C$ inside one of the two closed halfspaces it defines. A
face of $C$ is again a cone.

In Sections 3.5 and 4 the dual cone $C^{*}$ will play a central role. Its
definition does not only depend on the intrinsic structure of $C$, but also on
the ambient space. Therefore we will then assume that $C$ is a full
dimensional pointed cone: $\dim C=n$. Consequently the halfspaces $H_{i}^{+}$
in an irredundant representation (2.1) are uniquely determined, and there is a
unique primitive choice for $\sigma_{i}$. These linear forms $\sigma_{i}$ are
called the _support forms_ of $C$. In this note the hyperplanes $H_{i}$ are
the _support hyperplanes_ of $C$. The _dual cone_

$C^{*}=\\{\lambda\in({\mathbb{R}}^{n})^{*}:\lambda(x)\geq 0\text{ for all
}x\in C\\}$

is full dimensional and pointed as well. Under the natural identification
${\mathbb{R}}^{n}=({\mathbb{R}}^{n})^{**}$ the bidual cone $C^{**}$ is
identified with $C$: the extreme rays of $C^{*}$ are the linear forms defining
the facets of $C$, and vice versa. In the rational case the extreme integral
generators of $C^{*}$ are the support forms of $C$, and vice versa.

A _polytope_ $P$ is the convex hull of finitely many points in a real space
${\mathbb{R}}^{n}$. Our polytopes will be rational: such polytopes have
vertices in ${\mathbb{Q}}^{n}$. Computationally, polytopes are treated as
compact intersections of pointed cones and hyperplanes. The hyperplane is
defined by a linear form with integral coefficients, called _degree_ , such
that

$P=\\{x\in C:\deg x=1\\}.$ (2.3)

The intersection $P$ is compact (and nonempty) if and only if $C\neq 0$ and
$\deg x>0$ for $x\in C$, $x\neq 0$. This is not a restriction of generality:
if $P\subset{\mathbb{R}}^{n}$ is not given as in (2.3), then we can easily re-
embed it suitably: we identify $P$ with
$P^{\prime}=P\times\\{1\\}\subset{\mathbb{R}}^{n+1}$, and choose
$C={\mathbb{R}}_{+}P^{\prime}$.

###  2.2. Lattice normalized volume

Normaliz computes lattice normalized volume. We review this notion with
emphasis of its computation. The reader can find more details in [6, Sect. 3].
Let $P\subset{\mathbb{R}}^{n}$ be a rational polytope. The affine hull
$A=\operatorname{aff}P$ is a rational affine subspace of ${\mathbb{R}}^{n}$.
First assume that $0\in A$. Then $L=(\operatorname{aff}P)\cap{\mathbb{Z}}^{n}$
is a subgroup of ${\mathbb{Z}}^{n}$ of rank $d=\dim P$ (and
${\mathbb{Z}}^{n}/L$ is torsionfree). Choose a ${\mathbb{Z}}$-basis
$v_{1},\dots,v_{d}$ of $L$. The _lattice (normalized) volume_
$\operatorname{Vol}$ on $A$ is the Lebesgue measure on $A$ scaled in such a
way that the simplex $\operatorname{conv}(0,v_{1},\dots,v_{d})$ has measure
$1$. The definition is independent of the choice of $v_{1},\dots,v_{d}$ since
all invertible $d\times d$ matrices over ${\mathbb{Z}}$ have determinant $\pm
1$. If $0\notin A$, then we replace $A$ by a translate $A_{0}=A-w$, $w\in A$,
and set $\operatorname{Vol}X=\operatorname{Vol}(X-w)$ for $X\subset A$. This
definition is independent of the choice of $w$ since $\operatorname{Vol}$ is
translation invariant on $A_{0}$. Note that the polytope containing a single
point $x\in{\mathbb{Q}}^{n}$ has lattice volume $1$. If desired, the
definition of lattice volume can be extended to arbitrary measurable subsets
of $A$, and Normaliz does it for algebraic polytopes.

If $P$ is a lattice polytope, i.e., a polytope with vertices in
${\mathbb{Z}}^{n}$, then $\operatorname{Vol}P$ is an integer. For an arbitrary
rational polytope we have $\operatorname{Vol}P\in{\mathbb{Q}}$. As a
consequence, $\operatorname{Vol}P$ can be computed precisely by rational
arithmetic. This is not true for Euclidean volume in general: the diagonal of
the unit square has length $\sqrt{2}$.

A second invariant we need is the lattice height of a _rational_ point $x$
over a rational subspace $H\neq\emptyset$. More generally, one can consider
points $x$ such that $\operatorname{aff}(H,x)$ is again rational; for example,
this is the case if $H$ is a hyperplane in ${\mathbb{R}}^{n}$. If $x\in H$, we
set $\operatorname{Ht}_{H}(x)=0$. Otherwise let $A=\operatorname{aff}(H,x)$ so
that $H$ is a hyperplane in $A$.

Assume first that $0\in H$. Then $H$ is cut out from $A$ by an equation
$\lambda(y)=0$ with a primitive ${\mathbb{Z}}$-linear form $\lambda$ on
$L=A\cap{\mathbb{Z}}^{n}$. With this choice of $\lambda$,
$\operatorname{Ht}_{H}(x)=|\lambda(x)|$ is called the _lattice height_ of $x$
over $H$. (There are exactly two choices for $\lambda$, differing by the
factor $-1$.) If $0\notin H$, then we choose an auxiliary point $v\in H$,
replace $H$ by $H-v$, $A$ by $A-v$ and $x$ by $x-v$. In the algorithms we will
only have to deal with the case $0\in H$. If $P$ is a rational polytope and
$F$ is a facet or, more generally, a face of $P$, then we set
$\operatorname{Ht}_{F}(x)=\operatorname{Ht}_{H}(x)$ where
$H=\operatorname{aff}F$.

The following proposition relates lattice volume and lattice height.

######  Proposition 1.

Let $P$ be a rational polytope and $v\in P$ a vertex of $P$ such that there is
a single facet $F$ of $P$ with $v\notin F$. Then

$\operatorname{Vol}P=\operatorname{Ht}_{F}(v)\operatorname{Vol}F.$

This is part of [6, Prop. 1], to which we refer for the proof. The next basic
result tells us how to compute the volume of a _simplex_ , which is a polytope
of dimension $d$ with $d+1$ vertices.

######  Proposition 2.

Let $S\subset{\mathbb{R}}^{n}$ be a rational simplex with vertices
$v_{0},\dots,v_{d}$. Choose a basis $u_{1},\dots,u_{d}$ of the lattice
$\operatorname{aff}(S-v_{0})\cap{\mathbb{Z}}^{n}$. Define the $d\times d$
matrix $T=(t_{ij})$ by the representations
$v_{i}-v_{0}=\sum_{j=1}^{d}t_{ij}u_{j}$, $i=1,\dots,d$. Then

$\operatorname{Vol}S=|\det T|.$

This follows immediately from the transformation formula of Lebesgue measure.
See [5, 2.C] for an algebraic proof.

As mentioned already, we present rational polytopes $P$ in the form $P=C\cap
H$ where $C$ is a pointed cone and $H$ is defined by the condition $\deg x=1$
with a ${\mathbb{Z}}$-linear form $\deg$. This brings a second polytope into
play, namely $\overline{P}=\operatorname{conv}(0,P)$ as in Figure 1.

$P$$\overline{P}$$\deg=1$ Figure 1. $P$ and $\overline{P}$

All algorithms of Normaliz compute $\operatorname{Vol}\overline{P}$, and then
derive $\operatorname{Vol}P$ from it:

######  Proposition 3.

With the notation introduced, let $L={\mathbb{R}}C\cap{\mathbb{Z}}^{n}$ and
$\deg|L=k\deg^{\prime}$ with a primitive linear form $\deg^{\prime}$ on $L$
and $k>0$. Then

$\operatorname{Vol}P=k\operatorname{Vol}\overline{P}.$

###### Proof.

Let $F=P$ be the unique facet of $\overline{P}$ opposite to $0$. We can use
$\deg^{\prime}$ to measure lattice height over $F$. Since
$\deg^{\prime}x=(1/k)\deg x=1/k$ for $x\in F$, one has
$\operatorname{Ht}_{H}(0)=1/k$, and the claim follows from Proposition 1. ∎

The number $k$ in Proposition 3 is called the _grading denominator_ in
Normaliz. The reason is that $\deg^{\prime}=\deg/k$ is considered as the
“true” grading on the cone $C$. The user can choose between the given grading
$\deg$ or the divided one, $\deg^{\prime}$.

As our final tool we formulate a homogeneous version of Proposition 2:

######  Proposition 4.

Let the simplex $S$ be given in the form $S=C\cap H$ where $H$ is the
hyperplane of degree $1$ points and $C={\mathbb{R}}_{+}S$. Let
$v_{1},\dots,v_{d}$, $d=\dim S+1$, be nonzero points in the $d$ extreme rays
of $C$, for example the extreme integral generators. Then

$\operatorname{Vol}\overline{S}=\frac{1}{g_{1}\cdots g_{d}}|\det T|,\qquad
g_{i}=\deg v_{i},\ i=1,\dots,d,$

where $T=(t_{ij})$ is the $d\times d$ matrix with $v_{i}=\sum_{j}t_{ij}u_{j}$
for a basis $u_{1},\dots,u_{d}$ of the lattice
$L={\mathbb{Z}}^{n}\cap{\mathbb{R}}S$.

This follows immediately from Proposition 2 if we set $v_{0}=0$, observing
that $v_{1}/g_{1},\dots,\allowbreak v_{d}/g_{d}$ are the remaining vertices of
$\overline{S}$.

##  3\. Volume algorithms in Normaliz

There are three basic algorithms:

  1. (1)

the _primal_ volume algorithm: Normaliz computes a lexicographic
triangulation, and finds the volume as the sum of the volumes of the simplices
in the triangulation;

  2. (2)

volume by _descent in the face lattice_ : there is a reverse lexicographic
triangulation in the background, but it is not computed explicitly;

  3. (3)

volume by _signed decomposition_ , the _Lawrence algorithm_ : Normaliz
computes a triangulation of the dual cone and converts it into a signed
decomposition of the polytope.

Normaliz also computes the exact volume of full dimensional polytopes defined
over real algebraic number fields. For them only (1) is implemented at
present. One could extend (3) to them, whereas (2) is not suitable. The
algorithms (1) and (3) are also used in the computations of integrals of
rational polynomials over polytopes.

By rule of thumb one can say that the best choice is

  1. (1)

if the polytope has few vertices, but potentially many facets;

  2. (2)

if the number of vertices and the number of facets are of the same order of
magnitude;

  3. (3)

if there are _very_ few facets and many vertices.

This recommendation will be confirmed by the computational data in Section 5.
There are variants:

  1. (a)

exploitation of isomorphism types of faces in the descent algorithm;

  2. (b)

symmetrization as explained below.

Normaliz checks the default conditions of the algorithms in the order

signed decomposition $\to$ descent $\to$ symmetrization.

If the default conditions are not satisfied for any of them, the primal
triangulation algorithm is used. These decisions must often be made on the
basis of partial information. Therefore it can be useful to choose a certain
variant explicitly or to exclude others. The exploitation of isomorphism types
must always be asked for by the user.

Normaliz uses OpenMP for parallelization. Unless the user insists on
computations with GMP integers, Normaliz tries $64$ bit arithmetic first, and
restarts the computation with GMP integers if it recognizes an overflow.

###  3.1. The primal volume algorithm

Mathematically there is not much to say: if a polytope $P$ is decomposed into
simplices with non-overlapping interiors, then its volume is the sum of the
volumes of the simplices forming the decomposition.

Figure 2. A triangulation

Since the computation of Hilbert bases and Hilbert series is based on
(lexicographic) triangulations as well, Normaliz has a sophisticated algorithm
for them, using pyramid decomposition; see [8]. Normaliz tries to avoid
determinant computations by the “exploitation of unimodularity”; see [8, Prop.
7].

###  3.2. Volume by descent in the face lattice

The idea is to exploit the following proposition:

######  Proposition 5.

Let $P\subset{\mathbb{R}}^{n}$ be a rational polytope, and $v\in P$. Then

$\operatorname{Vol}P=\sum_{F\text{ facet of
}P}\operatorname{Ht}_{F}(v)\operatorname{Vol}F.$ (3.1)

Proposition 5 follows immediately from Proposition 1 since the polytopes
$\operatorname{conv}(v,F)$ constitute a polyhedral decomposition of $P$.
Usually $v$ is a vertex of the polytope $P$ with as few opposite facets
$F_{i}$ as possible, as illustrated by Figure 3.

$v$$F_{1}$$F_{2}$$F_{3}$$F_{4}$ Figure 3. Pyramid decomposition of a polytope

The recursive application results in building a _descent system_ , i.e., a
subset $\mathcal{F}$ of the face lattice so that for each face
$F\in\mathcal{F}$, to which (3.1) is applied, all facets of $F$ that are
opposite to the selected vertex are contained in $\mathcal{F}$. However, if a
face is simplicial, its multiplicity is computed by the standard determinant
formula. The algorithm is implemented in such a way that all data are
collected in the descent and no backtracking is necessary. The RAM usage is
essentially determined by the two largest layers. For a detailed discussion we
refer the reader to [6].

###  3.3. Exploitation of isomorphism classes

If the integral automorphism group of the cone $C$ over the polytope $P$ is
large enough, one can expect that each face in the descent system
${\mathcal{F}}$ has many isomorphic copies in ${\mathcal{F}}$. These can be
detected and identified so that only one representative of every isomorphism
class must be kept in ${\mathcal{F}}$. This reduces ${\mathcal{F}}$ in size
and can significantly speed up the volume computation. It must be used with
care since the computation of isomorphism classes is rather expensive. See [4]
for a discussion of their computation.

If the polytope is specified by generators and the number of facets is large,
then the first step in the descent system is built differently. Normaliz
computes the automorphism group of the polytope and selects one representative
in each orbit of facets. The vertex $v$ above is replaced by a fix point of
the automorphism group, and the first step in the volume computation is the
formula

$\operatorname{Vol}P=\sum_{i=1}^{c}O(F_{i})\operatorname{Ht}_{F_{i}}(v)\operatorname{Vol}F_{i}$

where $F_{i}$ represents one of the $c$ orbits and $O(F_{i})$ is the number of
facets in the orbit of $F_{i}$. Then $F_{1},\dots,F_{c}$ form the first layer
in the descent system. This allows the application of descent in cases where
the number of facets is too large for a successful computation without
exploitation of isomorphism classes. If the user does not prohibit it,
Normaliz encodes isomorphism classes by their SHA256 checksums.

###  3.4. Symmetrization

To understand the computation of volumes through symmetrization one must take
a detour through Ehrhart series. As usual, assume that our polytope $P$ is
given as the intersection $P=C\cap H$ where $C\subset{\mathbb{R}}^{d}$ is a
pointed rational cone and $H=\\{x\in{\mathbb{R}}^{d}:\deg x=1\\}$ is the
hyperplane of degree $1$ points. For symmetrization we assume that $\deg$ is
primitive.

Under certain conditions one can count lattice points of degree $k$,
$k\in{\mathbb{N}}$, in $C$ by mapping $C$ to a cone $C^{\prime}$ of lower
dimension and then counting each degree $k$ lattice point $y$ in $C^{\prime}$
with the number of its lattice preimages. This approach works well if the
number of preimages is given by a polynomial in the coordinates of $y$. Since
$C^{\prime}$ has lower dimension, one can hope that its combinatorial
structure is much simpler than that of $C$. One must of course pay a price:
instead of counting each lattice point with the weight $1$, one must count it
with a polynomial weight.

The availability of this approach depends on symmetries in the coordinates of
$C$, and therefore we call it _symmetrization_. Normaliz tries symmetrization
under the following condition: $C$ and the relevant lattice are given by
constraints (inequalities, equations, congruences) and the inequalities
contain the sign conditions $x_{i}\geq 0$ for all coordinates $x_{i}$ of $C$.
Then Normaliz groups coordinates that appear in all constraints and the
grading (!) with the same coefficients, and, roughly speaking, replaces them
by their sum. The number of preimages that one must count for the vector $y$
of sums is then a product of binomial coefficients—a polynomial as desired.
More precisely, if $y_{j}$, $j=1,\dots,m$, is the sum of $u_{j}$ variables
$x_{i}$ then

$f(y)=\binom{u_{1}+y_{1}-1}{u_{1}-1}\cdots\binom{u_{m}+y_{m}-1}{u_{m}-1}.$

is the number of preimages of $(y_{1},\dots,y_{m})$.

Since the Lebesgue measure can be approximated by scaled counting measures,
one obtains

$\operatorname{Vol}P=\int_{P}h\,d\negthinspace\lambda$

where $h$ is the highest homogeneous component of $f$ with respect to total
degree, and $\lambda$ is the suitably scaled Lebesgue measure. We learnt this
approach from Schürmann [23]. The Normaliz algorithm for integrals is
described in [11]. This note contains a complete elementary treatment and
several references to advanced aspects.

Symmetrization can have stunning effects. Nevertheless we do not include it in
the computations of Section 5 since it does not help for any of them, at least
not in the present implementation. Plenty of examples are contained in [9],
where it is often very useful in the computation of Hilbert series.

###  3.5. Volume by signed decomposition

This algorithm uses that a generic triangulation of the dual cone induces a
signed decomposition of the primal polytope, as we will now explain.

Let $C\subset{\mathbb{R}}^{d}$ be a pointed cone of dimension $d$ (it is
important that $C$ is full dimensional). The polytope $P$ is the intersection
of $C$ with the hyperplane $H$ defined by a grading $\deg$:
$H=\\{x:\deg(x)=1\\}$. The grading is an interior element of the dual cone
$C^{*}=\\{\lambda\in({\mathbb{R}}^{d})^{*}:\lambda(x)\geq 0\text{ for all
}x\in C\\}$. In order to visualize the situation we take an auxiliary
(irrelevant) cross-section $Q$ of the dual cone as in Figure 4.

$P$

$Q$$\deg$

Figure 4. A square and a cross-section of the dual cone

Now suppose that we have a _generic_ triangulation $\Delta$ of the dual cone
where genericity is defined as follows: $\deg$ is not contained in any
hyperplane through a facet of any $\delta\in\Delta$. Let $\delta\in\Delta$ be
given, and denote the linear forms on $({\mathbb{R}}^{d})^{*}$ defining its
facets by $\ell_{1},\dots\ell_{d}\in({\mathbb{R}}^{d})^{**}={\mathbb{R}}^{d}$.
( $\ell_{1},\dots\ell_{d}$ are the extreme rays of the dual of $\delta$.) The
hyperplanes defined by the vanishing of $\ell_{1},\dots\ell_{d}$ decompose
$({\mathbb{R}}^{d})^{*}$ into “orthants” that can be labeled by a sign vector
$\sigma=(s_{1},\dots,s_{d})\in\\{\pm 1\\}^{d}$:

$D(\delta,\sigma)=\\{\alpha:(-1)^{s_{i}}\ell_{i}(\alpha)\geq 0\\}.$

By the assumption on $\deg$, there is _exactly on_ e sign vector $\sigma$ such
that $\deg$ lies in the interior of $D(\delta,\sigma)$. Consequently the
hyperplane $H$ intersects the dual $D(\delta,\sigma)^{*}$ in a polytope
$R_{\delta}$. Set $e(\delta)=|\\{i:s_{i}=-1\\}|$.

Let $\iota_{X}$ denote the indicator function of a subset
$X\subset{\mathbb{R}}^{d}$. Then

$\iota_{P}(x)=\sum_{\delta\in\Delta}(-1)^{e(\delta)}\iota_{R_{\delta}}(x)$
(3.2)

for all $x\in{\mathbb{R}}^{d}$ outside a union of finitely many hyperplanes.
Since volume (lattice normalized or Euclidean) is additive on indicator
functions, this formula can be used for the computation $\operatorname{Vol}P$,
and more generally for the computation of integrals over $P$.

In order to find a generic triangulation, Normaliz first computes a
triangulation $\Delta$ of $C^{*}$ and saves the induced _hollow triangulation_
$\Gamma$ that $\Delta$ induces on the boundary of $C^{*}$. Then it finds a
_generic_ element $\omega\in C^{*}$ such that the _star triangulation_
$\Sigma$ of $C^{*}$, in which every simplicial cone is generated by the
_center_ $\omega$ and a facet of the hollow triangulation, is generic. Figure
5 illustrates the signed decomposition of a square into $4$ simplices.

$P$$+$$-$$-$$+$

$Q$$\deg$$\omega$$+$$+$$-$$-$

Figure 5. Generic triangulation of the dual and signed decomposition

The algorithm has been developed by Lawrence [18] in the language of linear
programming. We have learnt it from Filliman’s paper [15], which contains a
proof of equation (3.2). See Beck, Haase and Sottile [1] for an elementary
approach and the relationship to other duality theorems.

##  4\. The implementation of the Lawrence algorithm

The complete course of the volume computation consists of $4$ steps that can
be clearly delimited from each other:

  1. (1)

computation of a triangulation $\Delta$ of $C^{*}$;

  2. (2)

computation of the induced hollow triangulation $\Gamma$;

  3. (3)

choice of the generic element $\omega$;

  4. (4)

evaluation of the star triangulation $\Sigma$.

Each of the $4$ steps is highly critical if one wants to reach the
applications in social choice that were our driving challenge. For (1) we
could essentially rely on the standard triangulation algorithm of Normaliz.
Step (2) and the arithmetic for (3) and (4) are described in the following.
Both (3) and (4) are iterations over the hollow triangulation and star
triangulations derived from it.

###  4.1. The hollow triangulation

Suppose the triangulation $\Delta$ of $C^{*}$ has been computed. For each
simplicial cone $\delta\in\Delta$ we must now find the facets of $\delta$ that
lie in the boundary of $C^{*}$. There are various solutions for this task. The
first that comes to mind is to compute the facets of $C^{*}$ and match the
facets of $\delta$ with it. But $C^{*}$ can have an enormous number of facets
that one does better not compute since they can easily exhaust RAM. The facets
of $C^{*}$ are of course extreme rays of the cone $C$ over $P$, but for signed
decomposition Normaliz only computes them if asked for by the user. A second
approach that is much better in terms of RAM is to compute the facets of
$\delta$ and select those that have all extreme rays of $C^{*}$ on the same
side as $\delta$. However, this requires an enormous number of scalar products
that in high dimension are expensive.

Instead of these geometric approaches, Normaliz uses only the combinatorics of
triangulations of cones. It forms the set of all pairs $(F,\delta)$ where $F$
is a facet of $\delta$, $\delta\in\Delta$. Both components are encoded as
$0$-$1$-vectors that indicate the extreme rays of $C^{*}$ spanning $F$ and
$\delta$, respectively. From this set one must discard all pairs $(F,\delta)$
for which there exists a pair $(F,\delta^{\prime})$,
$\delta^{\prime}\in\Delta$, $\delta^{\prime}\neq\delta$. In principle one
could eliminate all $F$ that appear a second time without remembering the
“mother” $\delta$, but $\delta$ helps in several ways. The first is that one
can store the hollow triangulation as a set of pairs $(\delta,\phi(\delta))$
where $\phi(\delta)$ is a second $0$-$1$-vector indicating those extreme rays
of $\delta$ whose omission yields a facet of the hollow triangulation.

In order not to blow up memory for large $\Delta$, the pairs
$(\delta,\phi(\delta))$ are computed in small portions controlled by
“patterns”. Each pattern is an increasing sequence $(p_{1},\dots,p_{r})$ of
indices, and it is required that the facet $F$ satisfies the following
condition: if $q_{1},\dots,q_{d-1}$, $d=\dim C^{*}$, are the indices of the
extreme rays of $F$ in ascending order, then $q_{i}=p_{i}$ for $i=1,\dots,r$.

###  4.2. Piggyback simplices

After the purely combinatorial computation of the hollow triangulation,
arithmetic must be used in steps (3) and (4) above, namely in finding a
generic linear form $\omega\in C^{*}$ and then in the volume computation. Both
tasks are accelerated significantly if one takes advantage of the fact that
simplices $G$ and $G^{\prime}$ of the star triangulation are in “piggyback”
position to each other, if the facets $F$ and $F^{\prime}$ of the hollow
triangulation that define them belong to the same simplex $\delta\in\Delta$.
By “piggyback” position we mean that the simplices share a facet and lie on
different sides of it, as indicated in Figure 6.

$F$$F^{\prime}$$G$$G^{\prime}$$\lambda_{1}$$\lambda_{3}$$\lambda_{2}$$\mu$
Figure 6. Piggyback simplices

Suppose that $G$ has extreme rays
$\lambda_{1},\dots,\lambda_{d}\in({\mathbb{R}}^{d})^{*}$ and its facets are
given by linear forms
$\ell_{1},\dots,\ell_{d}\in{\mathbb{R}}^{d}=({\mathbb{R}}^{d})^{**}$ where
$\ell_{d}$ defines the common facet, $\lambda_{1},\dots,\lambda_{d-1}$ are the
extreme rays shared by $G$ and $G^{\prime}$ and $\mu$ is the remaining extreme
ray of $G^{\prime}$. Then the facets of $G^{\prime}$ are determined by
$-\ell_{d}$ and

$k_{i}=-\ell_{d}(\mu)\ell_{i}+\ell_{i}(\mu)\ell_{d},\qquad i=1,\dots,d-1.$
(4.1)

Since the computation of the $k_{i}$ and $-\ell_{d}$ from
$\lambda_{1},\dots,\lambda_{d-1},\mu$ alone amounts to the inversion of a
matrix, it is clear that the use of (4.1) is a significant advantage, even if
the computation of the values $\ell_{i}(\mu)$ needs $d^{2}$
multiplications.There is actually no need to compute the $k_{i}$ completely.
We will only need their values on elements in the dual space, for which the
values of the $\ell_{i}$ are known, for example the degree.

In the primal space, $\ell_{1},\dots,\ell_{d}$ and
$k_{1},\dots,k_{d-1},-\ell_{d}$ are extreme rays of the simplices dual to $G$
and $G^{\prime}$. For the volume of the corresponding simplices we need the
determinants. By standard rules

$|\det(k_{1},\dots,k_{d-1},-\ell_{d})|=|\ell_{d}(\mu)^{d-1}\det(\ell_{1},\dots,\ell_{d})|.$

So the piggyback relation between $G$ and $G^{\prime}$ pays off a second time.

In dealing with the simplices of the star triangulation that belong to the
same simplicial cone $\delta$ of $\Delta$, we pick one of them, say $G_{1}$,
and take all others piggyback. For $G_{1}$ we must indeed invert the matrix
$M$ with rows $\lambda_{1},\dots,\lambda_{d}$ over ${\mathbb{Q}}$, using
$D=|\det M|$ as the denominator: $M^{-1}=(1/D)N$ with a matrix
$N\in{\mathbb{Z}}^{d\times d}$. After extraction of their greatest common
divisors, the columns $\ell_{1},\dots,\ell_{d}$ of $N$ are the support forms
of $G_{1}$, equivalently, the extreme integral generators of the dual cone of
$G_{1}$. For the volume computation we need the determinant of the matrix
$N^{\prime}$ with columns $\ell_{1},\dots,\ell_{d}$. There is no need to
compute it directly: since $MN^{\prime}$ is the diagonal matrix with entries
$\ell_{i}(\lambda_{i})$, $i=1,\dots,d$, one has

$\det N^{\prime}=\frac{\prod_{i=1}^{d}\ell_{i}(\lambda_{i})}{\det M},$

and $\det M$ has already been computed.

###  4.3. Finding a generic element

The generic element $\omega$ of $C^{*}$ must satisfy the following condition:
the grading $\deg$ does not lie on any hyperplane through a facet of a
simplicial cone in the hollow triangulation and $\omega$. But this condition
is symmetric in $\omega$ and $\deg$! It is much better to first take the star
triangulation with center $\deg$ and check that $\omega$ is not on any of the
critical hyperplanes. The main difference between $\deg$ and $\omega$ is the
size of the coordinates: those of $\deg$ are usually very small and those of
$\omega$ very large. So, in working with $\deg$ as the center, there is a very
good chance to get away with $64$ bit arithmetic. The computation based on
$\omega$ which is necessary for the volume, must very often be done with GMP
integers.

Instead of choosing one vector $\omega$ at random and verifying that it is
generic, Normaliz takes two vectors $\omega_{1}$ and $\omega_{2}$ and checks
that not both of them lie together on a critical hyperplane. If this condition
is satisfied, then a suitable linear combination
$\omega=a_{1}\omega_{1}+a_{2}\omega_{2}$, $a_{1},a_{2}\in{\mathbb{Z}}$,
$a_{1},a_{2}\geq 0$, is generic. For the check we apply the piggyback trick of
Section 4.2, and especially (4.1). If $\omega_{1}$ and $\omega_{2}$ do not
work, then their coordinates are increased.

###  4.4. The addition of fractions

There is one more extremely critical aspect, namely the addition of fractions.
If one accumulates the volume as an alternating sum of simplex volumes by
successive addition of fractions, one can easily spend $99\%$ of the
computation time on this addition, or, in extreme situations, not finish at
all, creating fractions whose numerators and denominators fill gigabytes. The
reason is that the extreme integral generators of the duals of cones involving
the generic element can be very large: we must divide by them.

In extreme cases there is no other choice but to work with fixed precision. If
the user asks for it, then the simplex volumes are still computed exactly, but
for addition they are truncated to a fixed number of decimal digits. The
default choice is truncation to an integral multiple of $10^{-100}$. With this
choice the final volume is computed up to an error $\leq|\Gamma|\cdot
10^{-100}$. The user can set a higher or lower precision.

For computations with full precision, Normaliz uses an _addition pyramid_ , an
extensible vector $(q_{0},\dots,q_{n})$ of rational numbers. It serves as an
accumulator. A new summand $s$ is added to $q_{0}$, provided the number of
summands that have already been accumulated in $q_{0}$ is smaller than the
_capacity_. If the capacity has been reached, then $q_{0}$ is added to $q_{1}$
and set to $s$. The addition of $q_{0}$ to $q_{1}$ is handled in the same way
as that of $s$ to $q_{0}$ etc. At present the capacity is $8$. This scheme has
proved to be very efficient. Of course, at the very end, all entries of the
pyramid must be added.

######  Remark 6.

vinci [12] contains the Lawrence algorithm in a floating point implementation.
As the authors state in [13], it is numerically unstable, and our preceding
discussion should also indicate this problem.

If one has a look at the volumes whose alternating sum must be formed, then
they easily reach absolute values of $10^{100}$. For any precision of their
alternating sum, which may be of order $10^{-6}$, one therefore needs a very
high number of significant digits that a standard floating point format does
not offer.

For example, this becomes visible already in the comparison of $4$ voting
schemes for $4$ candidates [6, Sect. 6.1]. Let $P$ be the corresponding
polytope. With its algorithm HOT (based on the same principle as descent)
vinci correctly computes the Euclidean volume of $1.260510232743\cdot
10^{-25}$, for $\overline{P}$, whereas the vinci Lawrence algorithm yields
$9.287423132835\cdot 10^{-8}$. (We are grateful to Bogdan Ichim for these
computations.) For this reason we are not comparing the Lawrence algorithms in
Normaliz and vinci.

##  5\. Computational data

All computations have been done on the Dell R640 server of the Institute of
Mathematics at Osnabrück. It is equipped with two Intel Xeon Gold 6152 cards
(a total of $44$ cores) and $1$ TB of RAM. The computations use $32$ parallel
threads (of the maximum of $88$). The listed times are “wall clock” times. In
order to avoid overloading the tables, information about RAM usage has been
inserted into the text.

That we allow $32$ threads for a computation does of course not mean that they
can be used. The percentage of CPU that the computations got varies from $\sim
400\%$ to almost $3200\%$.

###  5.1. Polytopes defined by vertices

The computation times for several polytopes defined by vertices are listed in
Table 1. In the tables, $\dim$ is the dimension of the cone $C$ over the
polytope, #ext the number of extreme rays of $C$ and $\sup$ the number of its
support hyperplanes. The size of the triangulation computed by the primal
triangulation is to be found in the column #tri. The number of determinants
computed is usually smaller.

| combinatorial data |  | computation times
---|---|---|---
| dim | #ext | #supp | #tri |  | primal | descent | isotypes
( lo-6 | 16 | 720 | 910 | $5.8\cdot 10^{9}$ |  | 19:20.80 m | 3:17.48 m | 0:04.57 m
( lo-7 | 22 | 6040 | 87,472 |  |  |  |  | 21:39:51 h
( cr-20 | 21 | 40 | $2^{20}$ | $2^{19}$ |  | 0:08.50 m | 0:08.94 m | 0:15.02 m
( cr-24 | 25 | 48 | $2^{24}$ | $2^{23}$ |  | 2:11.91 m | 3:29.67 m | 5:42.24 m
( cr-28 | 29 | 48 | $2^{28}$ | $2^{27}$ |  | 42:12.11 m | 1:39:37 h | 2:21:09 h
( A543 | 36 | 60 | 29,387 | $103\cdot 10^{6}$ |  | 0:24.09 m | 36:56.56 m | 0:18.59 m
( A553 | 43 | 75 | 306,955 | $9.2\cdot 10^{9}$ |  | 44:53.26 m |  | 7:10.36 m
( cy-60 | 17 | 60 | 656,100 |  |  | 0:46.45 m |  | 0:44.27 m
Table 1. Polytopes defined by vertices

####  5.1.1. Linear ordering polytopes

`Lo<n>` is the linear ordering polytope for a set of $n$ elements. These
polytopes have been investigated in combinatorial optimization; see [20]. The
maximum $n$ reachable is $7$. For $n=8$ not even the number of facets is
known. It is however $>800\cdot 10^{6}$. The computation of the volumes is
surprisingly fast if one exploits the isomorphism classes of faces. Note that
for $n=7$ the computation of the support hyperplanes alone takes $>20$ h so
that the computation of the volume needs $\sim 1$ h.

The maximum RAM usage of `lo-6` is $2.1$ GB for descent, the other two
algorithms need $<1.5$ GB. The computation for `lo-7` takes $14.7$ GB.

####  5.1.2. Cross polytopes

`cr-<n>` is the unit cross polytope of dimension $n$. We have computed their
(known) volumes for $n=20,24,28$. They have only $2n$ vertices, but $2^{n}$
facets. But all facets are simplices, and therefore the descent algorithm and
its variant exploiting isomorphism classes are applicable. However, the primal
algorithm behaves better for two reasons: (i) it avoids the administrative
overhead of the descent algorithm, and (ii) the formation of the single orbit
of facets takes rather long—it cannot be parallelized. So the saving in the
computations of determinants is overcompensated. Since all facets of a cross
polytope are simplicial, there is only one descent step, namely from the full
polytope to the facets opposite to the chosen vertex.

For $n=20$ the primal algorithm gets away with $713$ MB, whereas the two
descent algorithms need about $1$ GB. For $n=24$ the numbers are $4.2$ GB and
$21$ GB. For $n=28$ they rise to $194$ GB and $\sim 300$ GB. It takes a lot of
space to accommodate the $2^{28}$ extreme rays.

####  5.1.3. Other polytopes

`A543` and `A553` are taken from the Ohsugi-Hibi classification [22] of
polytopes related to contingency tables. `A553` shows that descent with
isomorphism types can be favorable if the automorphism group of the polytope
is sufficiently large. For `A543` this effect is already visible, but still
small. This applies to `cy-60` as well, the cyclotomic polytope of order $60$
defined by Beck and Hoşten [2]. That the pure descent algorithm is not
suitable for this type of polytope is shown by `A543`.

The RAM usage of `A543` is about $1.6$ GB for the primal algorithm and $1$ GB
for descent with isomorphism classes. For `A553 ` the corresponding numbers
are $4.3$ GB and $101$ GB. For `cy-60` they are $316$ MB and $1.4$ GB.

###  5.2. Polytopes defined by inequalities

We now turn to polytopes defined by inequalities and equations. Among them we
have chosen Birkhoff polytopes, cubes and polytopes from social choice—as said
already, the latter were the driving challenge for our implementation of the
Lawrence algorithm. In the tables, #tri dual is the size of the triangulation
of the dual cone, and #hollow that of the associated hollow triangulation.

####  5.2.1. Birkhoff polytopes

The Birkhoff polytope of order $n$ is the set of doubly stochastic $n\times n$
matrices. Its vertices are the $n\times n$ permutation matrices, and their
number $n!$ is rapidly growing. Their volumes have been computed for $n\leq
10$ by Beck and Pixton [3] with residue methods that are not (yet) available
in Normaliz.

For $n\leq 5$ any of the Normaliz algorithm does the job very quickly, but for
$n=6$ the primal algorithm must already give up since the lexicographic
triangulation becomes too large. In the table we start with this case. The
Lawrence algorithm reaches $n=8$. As one can see, even the triangulations of
the dual cone grow too quickly for the next step. The bulk of the computation
time for $n=8$ is taken by the computation of the hollow triangulation, namely
$\sim 13$ h. The coordinates of the generic element $\omega$ are small enough
to allow 64 bit arithmetic for the volume computation ($<3$ h).

| combinatorial data |  | computation times
---|---|---|---
| dim | #ext | #supp | #tri dual | #hollow |  | isotypes | signed dec
( `bi-6` | 26 | 720 | 36 | 142,755 | 933,120 |  | 0:03.64 m | 0:03.10 m
( `bi-7` | 37 | 5040 | 49 | $11\cdot 10^{6}$ | $85\cdot 10^{6}$ |  | 38:05.85 m | 5:36.92 m
( `bi-8` | 50 | 40,320 | 64 | $1.2\cdot 10^{9}$ | $11\cdot 10^{9}$ |  |  | 17:59:20 h
( `cu-20` | 21 | $2^{20}$ | 40 | $2^{19}$ | $2^{20}$ |  | 0:08.73 m | 0:07.00 m
( `cu-24` | 25 | $2^{24}$ | 48 | $2^{23}$ | $2^{24}$ |  | 4:56.11 m | 2:12.64 m
( `cu-28` | 29 | $2^{28}$ | 56 | $2^{27}$ | $2^{28}$ |  | 1:53:27 h | 1:01:37 h
Table 2. Birkhoff polytopes and cubes

RAM usage for `bi-6` is $44$ MB and $473$ MB. The computations for `bi-7` take
$7.9$ GB and $12.5$ GB. That for `bi-8` needs $216$ GB.

####  5.2.2. Cubes

`cu-<n>` is the unit cube of dimension $n$. It is a good test object since its
volume is known. Since faces of the same codimension are isomorphic, descent
with isomorphism types is expected to be fast, and it is indeed. However, it
must use the huge number of vertices explicitly, and for this reason signed
decomposition is even faster. The Normaliz binary in the distribution does
never reach any of these algorithms since Normaliz recognizes parallelotopes
$P$, computes the volume of a “corner simplex” and multiplies it by $n!$,
$n=\dim P$. This takes $\sim 0.01$ s, even for $n=28$. It would certainly be
possible to go to $n=32$ with the Lawrence algorithm or descent with
isomorphism types.

In [6] the reader can find performance data for the descent algorithm applied
to `cu-20` and `cu-24`. In addition, more general parallelotopes of the same
dimensions are computed there. The computation times show that the arithmetic
is secondary and the times are dominated by the combinatorial complexity.

RAM usage for `cu-20` is $\sim 1.1$ GB for both algorithms, and for `cu-24` we
need $23$ GB and $2$ GB. `cu-28` takes $428$ GB and $127$ GB. That signed
decomposition is so much better, is due to the fact that it does not store the
extreme rays.

####  5.2.3. Polytopes from social choice

Computational data for polytopes from social choice are contained in [9] for
the primal algorithm and symmetrization and in [6] for descent in the face
lattice. Voting schemes with $5$ candidates are essentially inaccessible to
them. The Lawrence algorithm has now reached them, and [7] contains data for
them. Tables 3 4 are imported from there. The names of the polytopes are
explained in [7].

For two polytopes we have included the number of extreme rays to show the
order of magnitude. For the Lawrence algorithm they are not needed explicitly,
and in particular they need not be stored.

| dim $C$ | #ext |  $\\#$ supp |  $\\#$ tri dual |  $\\#$ hollow
---|---|---|---|---|---
( strictBorda 4cand | 24 |  | 33 | 100,738 | 324,862
( CondEffAppr 4cand | 74 |  | 80 | 1,620,052 | 30,564,920
( Condorcet | 120 | 290,064 | 124 | 137,105 | 6,572,904
( PlurVsRunoff | 120 | 80,912,472 | 125 | 4,912,369 | 93,749,784
( CWand2nd | 120 |  | 126 | 15,529,730 | 608,572,514
( CondEffPlurRunoff | 120 |  | 127 | 246,310,369 | 5,456,573,880
( CondEffPlur | 120 |  | 128 | 2,388,564,481 | 39,390,184,920
Table 3. Combinatorial data

The stages (1)–(3) of all computations could be done by $64$ bit arithmetic,
and this holds even for the volume computations of the first and third
polytope. The volume computations of the two largest had to be done with fixed
precision.

For the two largest examples it was necessary to use distributed computation
on a high performance cluster (indicated by HPC). For this reason we have
split the computation times. For `CondEffPlur` on the HPC of the University of
Osnabrück the time was $<9$ h. We refer the reader to [7] for more
information.

| RAM | time
---|---|---
| in GB | stages (1) – (3) | stage (4) | total
( strictBorda 4cand | 0.35 | 1.278 s | 0.464 s | 1.742 s
( CondEffAppr 4cand | 7.4 | 97.8 s | 14:31 m | 16:09 m
( Condorcet | 1.67 | 18.0 s | 52.493 s | 1:10 m
( PlurVsRunoff | 26.2 | 12:40 m | 1:29:21 h | 1:42:01 s
( CWand2nd | 56.4 | 49:55 m | 10:21:36 h | 11:11:31 h
( CondEffPlurRunoff | 113 | 13:30:22 h | HPC | —
( CondEffPlur | 646 | 125:27:20 h | HPC | —
Table 4. Memory usage and times for parallelized volume computations

## References

  * [1] M. Beck, Ch. Haase and F. Sottile, _Formulas of Brion, Lawrence, and Varchenko on rational generating functions for cones._ Math. Intelligencer 31 (2009), 9–17.
  * [2] M. Beck and S. Hoşten, Cyclotomic polytopes and growth series of cyclotomic lattices. Math. Res. Lett. 13 (2006), 607–622.
  * [3] M. Beck and D. Pixton, _The Ehrhart polynomial of the Birkhoff polytope._ Discrete Comput. Geom. 30 (2003), 623–637.
  * [4] W. Bruns, _Automorphism groups and normal forms in Normaliz._ Prepr4int arXiv:2112.08145.
  * [5] W. Bruns and J. Gubeladze, _Polytopes, rings and K-theory_. Springer, 2009.
  * [6] W. Bruns and B. Ichim, _Polytope volume by descent in the face lattice and applications in social choice._ Math. Program. Comput. 13 (2021), 415–442.
  * [7] W. Bruns and B. Ichim, _Computations of volumes in five candidates elections._ Preprint arXiv:2109.00473.
  * [8] W. Bruns, B. Ichim and C. Söger, The power of pyramid decomposition in Normaliz. Journal of Symbolic Computation 74 (2016), 513–536.
  * [9] W. Bruns, B. Ichim and C. Söger, Computations of volumes and Ehrhart series in four candidates elections. Annals of Operations Research 280 (2019), 241–265.
  * [10] W. Bruns, B. Ichim, C. Söger and U. von der Ohe, _Normaliz. Algorithms for rational cones and affine monoids._ Available at https://normaliz.uos.de.
  * [11] W. Bruns and C. Söger, Generalized Ehrhart series and Integration in Normaliz. J. Symb. Comp. 68 (2015), 75–86.
  * [12] B. Büeler and A. Enge, _Vinci._ Package available from https://www.math.u-bordeaux.fr/~aenge/
  * [13] B. Büeler, A. Enge and K. Fukuda, _Exact volume computation for polytopes: a practical study._ In: Polytopes - combinatorics and computation (Oberwolfach, 1997), pp. 131–154, DMV Sem., 29, Birkhäuser, Basel, 2000.
  * [14] D. Delfino, A. Taylor, W. V. Vasconcelos, N. Weininger, and R. H. Villarreal, _Monomial ideals and the computation of multiplicities_. In Commutative ring theory and applications (Fez, 2001), 87–106, M. Dekker, 2003.
  * [15] P. Filliman, _The volume of duals and sections of polytopes._ Mathematika 39 (1992), 67–80.
  * [16] W.V. Gehrlein and D. Lepelley, _Voting Paradoxes and Group Coherence._ Springer, 2011.
  * [17] W.V. Gehrlein and D. Lepelley, _Elections, Voting Rules and Paradoxical Outcomes._ Springer, 2017.
  * [18] J. Lawrence, _Polytope volume computation._ Mathematics of Computation 57 (1991), 259–271.
  * [19] D. Lepelley, A. Louichi and H. Smaoui, _On Ehrhart polynomials and probability calculations in voting theory._ Social Choice and Welfare 30 (2008), 363–383.
  * [20] R. Martí and G. Reinelt, _The linear ordering problem. Exact and heuristic methods in combinatorial optimization_. Springer, 2011.
  * [21] B. D. McKay and A. Piperno, _Practical graph isomorphism, II_. J. Symbolic Comput. 60 (2014), 94–112.
  * [22] H. Ohsugi and T. Hibi, Toric ideals arising from contingency tables. In: Commutative Algebra and Combinatorics. Ramanujan Mathematical Society Lecture Note Series 4 (2006), 87–111.
  * [23] A. Schürmann, _Exploiting polyhedral symmetries in social choice._ Social Choice and Welfare 40 (2013), 1097–1110.
  * [24] R. H. Villarreal, _Monomial algebras_. Second edition. CRC Press 2015.