Department of Informatics & Telecommunications
National & Kapodistrian University of Athens<EMAIL_ADDRESS>Department of Informatics & Telecommunications
National & Kapodistrian University of Athens<EMAIL_ADDRESS>Department of Informatics & Telecommunications
National & Kapodistrian University of Athens, Greece<EMAIL_ADDRESS>Apostolos Chalkis, Ioannis Z. Emiris and Vissarion Fisikopoulos

###### Acknowledgements.

# Practical volume estimation by a new annealing schedule for cooling convex
bodies

Apostolos Chalkis  ,  Ioannis Z. Emiris  and  Vissarion Fisikopoulos

###### Abstract.

We study the problem of estimating the volume of convex polytopes, focusing on
H- and V-polytopes, as well as zonotopes. Although a lot of effort is devoted
to practical algorithms for H-polytopes there is no such method for the latter
two representations. We propose a new, practical algorithm for all
representations, which is faster than existing methods. It relies on Hit-and-
Run sampling, and combines a new simulated annealing method with the
Multiphase Monte Carlo (MMC) approach.

Our method introduces the following key features to make it adaptive: (a) It
defines a sequence of convex bodies in MMC by introducing a new annealing
schedule, whose length is shorter than in previous methods with high
probability, and the need of computing an enclosing and an inscribed ball is
removed; (b) It exploits statistical properties in rejection-sampling and
proposes a better empirical convergence criterion for specifying each step;
(c) For zonotopes, it may use a sequence of convex bodies for MMC different
than balls, where the chosen body adapts to the input. We offer an open-
source, optimized C++ implementation, and analyze its performance to show that
it outperforms state-of-the-art software for H-polytopes by Cousins-Vempala
(2016) and Emiris-Fisikopoulos (2018), while it undertakes volume computations
that were intractable until now, as it is the first polynomial-time, practical
method for V-polytopes and zonotopes that scales to high dimensions (currently
$100$).

We further focus on zonotopes, and characterize them by their order (number of
generators over dimension), because this largely determines sampling
complexity. The number of phases in MMC tends to constant as order increases,
while the bodies in (c) are balls. For low orders, the generators’ matrix is
used to define a centrally symmetric convex polytope in (c). We analyze a
related application, where we evaluate methods of zonotope approximation in
engineering.

###### Key words and phrases:

Polytope volume, zonotope, sampling, simulated annealing, V-polytopes

###### 1991 Mathematics Subject Classification:

Design and analysis of algorithms:
Computational geometry, Random walks and Markov chains

##  1\. Introduction

As a special case of integration, volume computation is a fundamental problem
with many applications in science and engineering. From a computational
complexity point of view it is hard even if we restrict to convex sets. In
particular, it is $\\#$P-hard for H- and V-polytopes, including zonotopes
[13]. It is even hard to approximate, namely, APX-hard [10]. Therefore, a
great effort has been devoted to randomized approximation algorithms, starting
with the celebrated result in [9] with complexity $O^{*}(d^{23})$, where
$O^{*}(\cdot)$ indicates a soft-big-Oh hiding polylog factors. It introduced
the Multiphase Monte Carlo (MMC) technique, which reduced volume approximation
to computing a telescoping product of volumes, estimated by uniformly sampling
a sequence of convex bodies (by means of random walks). The following years,
improved algorithms reduced the exponent of dimension $d$ down to $5$ [17].
The latter led to the first practical implementation [11] for high dimensions,
which highlighted the importance of Coordinate Direction Hit-and-Run (HnR).
Further results in [6, 18] reduced the exponent to $3$, followed by another
practical method [7].

This paper proposes a new practical volume estimation for convex polytopes,
improving upon state-of-the-art methods for H-polytopes, while yielding the
first method capable to scale in the case of V-polytopes and zonotopes. We use
an adaptive sequence of convex bodies, simulated annealing, and the
statistical properties of the telescoping ratios in order to drastically
reduce the number of phases in MMC as well as the sample size required to
estimate these ratios. Our aim is to fully exploit probabilistic methods
within the current paradigm of the MMC approach, and optimize the resulting
software by careful algorithmic engineering.

Notation. $P$ is a full-dimensional convex polytope lying in $d$-dimensional
space. An H-polytope (in H-representation) is $P=\\{x\ |\ Ax\leq b,\
A\in{\mathbb{R}}^{q\times d},\ b\in{\mathbb{R}}^{q}\\}$. A V-polytope is the
convex hull of a pointset in ${\mathbb{R}}^{d}$. A zonotope (Z-polytope) is
the Minkowski sum of $k$ $d$-dimensional segments or equivalently given by
matrix $G\in{\mathbb{R}}^{d\times k}$ and seen as a linear map of hypercube
$[-1,1]^{k}$ to ${\mathbb{R}}^{d}$: we call it a Z-representation. The order
of a zonotope is the ratio $k/d$. We study sequences of bodies $C_{i}$
intersecting $P$ where the corresponding convex bodies $P_{i}=C_{i}\cap P$
define the telescoping product.

Membership oracles. A point $x_{0}\in P$ if and only if $Ax_{0}\leq b$ when
$P$ is a H-polytope. If $P$ is a V-polytope we follow [12] and solve the
linear program (LP) below,

max $\displaystyle f=z^{T}x_{0}-z_{0}$ subject to $\displaystyle
z^{T}v_{i}-z_{0}\leq 0,\ \text{for all }i=1,\dots,\varv$ $\displaystyle
z^{T}x_{0}-z_{0}\leq 1$

The $x_{0}\in P$ holds if and only if the optimal value $f$ of the LP is
strictly positive and the LP computes a separating hyperplane. If $P$ is a
zonotope we solve the following linear feasibility problem,

find $\displaystyle\lambda$ satisfying $\displaystyle G\lambda=x_{0}$
$\displaystyle-1\leq\lambda_{i}\leq 1,\ \text{for all }i=1,\dots,k$

The $x_{0}\in P$ holds if and only if the answer is positive.

Paper structure. The rest of the section presents previous work as well as our
contributions. Section 2 discusses our method, while Section 3 presents our
implementation, evaluates its practical complexity, compares to existing
software and offers a concrete application.

###  1.1. Previous work

The prevalent paradigm in volume approximation relies on MMC and sampling with
random walks. We build on the approach which defines a sequence of convex
bodies $P_{0}\subseteq\dots\subseteq P_{m}=P$ such that rejection sampling
would efficiently estimate $\mbox{vol}(P_{i+1})/\mbox{vol}(P_{i})$. Assuming
$P$ is well-rounded, i.e. $B_{d}\subseteq P\subseteq C\sqrt{d}B_{d}$, one
defines a sequence of scaled copies of the unit ball $B_{d}$, and
$P_{i}=(2^{i/d}B_{d})\cap P,\ i=0,\dots,m$. Then, it suffices to compute
$\mbox{vol}(P_{0})$ and apply the following telescopic product from [17]:

$\mbox{vol}(P)=\mbox{vol}(P_{0})\frac{\mbox{vol}(P_{1})}{\mbox{vol}(P_{0})}\cdots\frac{\mbox{vol}(P_{m})}{\mbox{vol}(P_{m-1})},\quad
m=O(d\lg d)$

In practical implementations [11], assuming $rB_{d}\subseteq P\subseteq
RB_{d}$, the construction gives $m=\lceil d\lg(R/r)\rceil$, see Figure 1
(left). The critical complexity issue is to define a sequence that minimizes
$m$ while each ratio remains bounded by a constant. This would permit a larger
approximation error per ratio without compromising overall error, while it
would require a smaller sample.

Each ratio is estimated by sampling uniform points from $P_{i}$, obtained by
random walks. The main approach today being HnR while, more recently, a
convergence rate is given for Hamiltonian walk [16] applicable to H-polytopes
only. The Vaidya walk is even faster [JMLR:v19:18-158] when the number of
facets $\gg d$. A recent Hamiltonian walk with reflections [5] reduces mixing
time: It can enhance our method for H-polytopes and may offer better stability
than Coordinate Directions HnR; it can also be integrated to the two existing
implementations discussed below. But its application to V- and Z-polytopes is
unclear.

In [18] they construct and estimate the volume of a $d+1$ dimensional convex
body, called the "pencil", and then they use rejection sampling to estimate
the volume of the input convex body. Moreover, they generalize the telescopic
product to a sequence of functions by fixing a sequence of exponential
distributions approximating the uniform; the total complexity is
$O^{*}(d^{4})$. In [6], which holds the current record in asymptotic analysis
for volume approximation of general convex bodies, they consider a sequence of
spherical Gaussian distributions on $P$; the total complexity is
$O^{*}(d^{3})$. The sequence of spherical Gaussians is not deterministic, but
approaches the uniform distribution fast. In [16] they improve the asymptotic
complexity to $O^{*}(qd^{\frac{2}{3}})$ for H-polytopes as they combine
Hamiltonian Monte Carlo with a sequence of Gibbs distributions instead of
Gaussians.

Current state-of-the-art software handles H-polytopes based on the above
paradigms and, typically, HnR. The software of [11] scales up to hundreds of
dimensions and uses Coordinate-Direction HnR. We also juxtapose the software
of [7] (for H-polytopes), implementing [6] with an annealing schedule [18] of
a sequence of Gaussians. For the ratio, they use a sliding window and stop
sampling when the maximum and minimum values meet a convergence criterion.

The main reason that these implementations cannot handle efficiently Z- nor
V-polytopes (cannot scale beyond, say, $d\geq 15$) is that the boundary and
membership oracles require to solve linear programs (LPs). Moreover both of
them require inscribed balls. If $P$ is a zonotope, checking whether a ball
$B\subseteq P$ is in co-NP, but it is not known whether it is co-NP-complete.
When $P$ is a V-polytope, given $p\in P$ the computation of the largest
inscribed ball centered at $p$ is NP-hard [19]. Additionally, the software of
[7] requires the number of facets which is typically exponential in the
dimension for both Z- and V-polytopes (Section 3 and Figure 6).

To sum up the discussion on previous work we should mention the rich area of
implementations of deterministic algorithms; notable examples are VINCI and
qHull but the list is quite long. As expected, those implementations do not
scale beyond, say, $d\geq 15$ dimensions for general polytopes.

For zonotopes and V-polytopes, computing the largest inscribed ball is a key
issue for both methods. If $P$ is a zonotope, checking whether a ball
$B\subseteq P$ is in co-NP, but it is not known whether it is co-NP-complete.
When $P$ is a V-polytope given $p\in P$ the computation of the largest
inscribed ball centered at $p$ is NP-hard [19].

Figure 1. Balls computed by MMC in [17] (left) and by the annealing schedule
of Algorithm 2 with $r$=0.25 and $\delta$=0.05 (right).

###  1.2. Our contribution

Our contribution is a new volume approximation method for H-, V-, and
Z-polytopes, also applicable to general convex bodies, though in this paper we
focus on convex polytopes. For V- and Z-polytopes, our algorithm requires
solving two (related) LPs per step of the random walk (Section 2.3). However,
we drastically reduce the number of such steps, hence offering the first
practical algorithm for such bodies in high dimensions. Regarding LPs, each
step solves two LPs with a common basic feasible solution. In Section 3 we
experimentally analyze our method to show it scales up to 100 dimensions for
V-polytopes and low order zonotopes. Hence, it outperforms both
implementations in [11, 7] on V- and Z-polytopes (Figure 6). In fact, it
performs volume computations which were intractable until now (Table 2). On
H-polytopes our method is faster than [11] for every $d$ and faster for $d\leq
100$ than [7] (Table 3 and Figure 6). The main algorithmic features follow:

  * •

We design a new simulated annealing method for cooling convex bodies in MMC
(Section 2.1). We exploit the fact that rejection sampling works efficiently
for bounded ratios, which are actually smaller than those in [17] (see Figure
1). To reduce the number of phases, we employ statistical tests to bound each
$r_{i}=\mbox{vol}(P_{i+1})/\mbox{vol}(P_{i})$ with high probability. Our
annealing schedule does not need an enclosing body of $P$ as they do in [17,
11]: it suffices to set $P_{0}=P$. In addition, our method does not require
computing an inscribed ball, as older methods did: the ball (or any body we
use in MMC) with minimum volume is computed by the annealing schedule, thus
further reducing the number of phases in practice. Finally, we prove that the
annealing schedule terminates successfully with constant probability (Section
2.1). This adaptive MMC sequence reduces significantly the number of phases:
we prove that this number in [17, 11] upper bounds the number of phases in our
method with high probability (Section 2.4).

  * •

When we sample $N$ uniform points from $P_{i}$, the number of points in
$P_{i+1}$ follows the binomial distribution. The main task here is to estimate
the ratio $\mbox{vol}(P_{i+1})/\mbox{vol}(P_{i})$ with minimum $N$. We exploit
the binomial proportion confidence interval and modify it by using the
standard deviation of a sliding window, in order to specify a new empirical
convergence criterion for the ratio, as $N$ increases (Section 2.2); this
drastically reduces sample size (Figure 6). An analogous technique was used in
[7], but the window here is of half the length.

  * •

We allow other convex bodies besides balls in MMC, the choice being a function
of the input, aiming at reducing the number of phases. For H- and V-polytopes
we actually use balls as in classical MMC. We leave it as an open question
whether there are more suitable bodies such as those defined in [2]. For
zonotopes we show that, as order grows, most suitable is the ball: our method
requires only a small constant number of them (Figure 3 and Table 2). However,
for low order, e.g. $\leq 4$, we use the matrix of generators to define a
centrally-symmetric H-polytope that accelerates the algorithm so as to scale
to, say, $d=100$, which used to be intractable (Section 2.5). For instance,
for a random zonotope with $k=2d$ generators, $d=100$, our software takes
$<10$ hr (Table 2).

We prove that, in our method, the number of phases is
$m=O(\lg(\mbox{vol}(P)/\mbox{vol}(B_{d})))$ with high probability (Section
2.4), when we use balls in MMC. This yields $m=O(d)$ for some well-rounded
polytopes such as the cross polytopes, a clear improvement over the general
bound of $O(d\lg d)$ [17]. Specifically, for a cross V-polytope in $d=100$ our
method requires just two balls and takes $406$ sec (Figure 5 and Table 1),
while the problem is intractable under the H-representation. Moreover, when
applying a rounding step to random V-polytopes our method requires very small
$m$, in practice $\leq 4$, even for $d=100$ (Section 3, and Table 1). If $P$
is a zonotope, we experimentally show that, for constant $d$, and $k$ (number
of generators) increasing, the number of phases $m$ decreases to 1 (Table 2
and Figure 3). An intuitive property of zonotopes is that, while order
increases for constant $d$, a random zonotope approximates the hypersphere. In
[3] they prove that for $d\geq 2$ the unit ball $B_{d}$ can be approximated up
to $\epsilon$ in the Hausdorff distance by a zonotope defined by $k$ segments
of equal length, $k\leq c(d)(\epsilon^{2}|\lg\epsilon|)^{(d-1)/(d+2)}$, where
$c(d)$ is a constant. This result cannot be used straightforwardly to prove
our claim but strengthens it intuitively. For some instances, when order is
large, our method creates just one ball and the method reduces to just one or
two rejection-sampling steps (Table 2).

Last but not least, we offer an open source efficient implementation of the
new method in
C++111https://github.com/GeomScale/volume_approximation/tree/CoolingBodies.

##  2\. Volume algorithm

Our volume algorithm relies on Multiphase Monte Carlo (MMC) method and samples
from uniform target distribution with Hit-and-Run (HnR) random walk. Hence,
the first part of the algorithm (Algorithm 1) is to construct a sequence of
convex bodies $C_{1}\supseteq\cdots\supseteq C_{m}$ intersecting the given
polytope $P$ using simulated annealing (Algorithm 2). Then it estimates each
ratio in the telescopic product of Equation (1) using HnR and a new empirical
criterion convergence (Algorithm 3). A typical choice for the $C_{i}$’s is a
sequence of co-centric balls but other convex bodies may be used (see, e.g.,
Section 2.5). Then,

$\mbox{vol}(P)=\frac{\frac{\mbox{vol}(P_{m})}{\mbox{vol}(C_{m})}}{\frac{\mbox{vol}(P_{1})}{\mbox{vol}(P_{0})}\frac{\mbox{vol}(P_{2})}{\mbox{vol}(P_{1})}\cdots\frac{\mbox{vol}(P_{m})}{\mbox{vol}(P_{m-1})}}\mbox{vol}(C_{m}),\
\mbox{ where }P_{0}=P,P_{i}=C_{i}\cap P.$ (1)

In the sequel, we write $r_{i}=\mbox{vol}(P_{i+1})/\mbox{vol}(P_{i}),\
i=0,\dots,m-1$ and $r_{m}=\mbox{vol}(P_{m})/\mbox{vol}(C_{m})$.

Algorithm 1 VolumeAlgorithm ($P$, $\epsilon$, $r$, $\delta$, $\alpha$, $\nu$,
$N$, $k$)

Construct convex body $C\subseteq{\mathbb{R}}^{d}$ s.t. $C\cap P\neq\emptyset$
and set interval $[q_{\min},q_{\max}]$

$\\{P_{0},\dots,P_{m},C_{m}\\}=$ AnnealingSchedule
$(P,C,r,\delta,\alpha,\nu,N,q_{\min},q_{\max})$

Set $\epsilon_{i},\ i=0,\dots,m$ s.t.
${\sum_{i=0}^{m}\epsilon_{i}^{2}}=\epsilon^{2}$

for $i=0,\dots m$ do

if $i<m$ then

$r_{i}=$ EstimateRatio ($P_{i},P_{i+1},\epsilon_{i},m,k$)

else

$r_{m}=$ EstimateRatio ($C_{m},P_{m},\epsilon_{m},m,k$)

end if

end for

return $\mbox{vol}(C_{m})/r_{0}/\cdots/r_{m-1}\cdot r_{m}$

The behavior of Algorithm 1 is parameterized by: the error of approximation
$\epsilon$, cooling parameters $0<r+\delta<1$, $r,\delta>0$, which are used in
the schedule, significance level (s.l.) $\alpha>0$ of the statistical tests,
$\nu$ the degrees of freedom for the t-student used in t-tests (all in Section
2.1), and parameter $N$ that controls the number of points $\nu N$ generated
in $P_{i}$. We generate uniform samples in $P_{i}$ using HnR (Section 2.3).

Following the telescopic product of Equation (1), it is clear that in
practical estimations $C_{m}$ has to be a convex body whose volume is obtained
much faster than $\mbox{vol}(P)$ (ideally by a closed formula) and easy to
sample. The maximum number of constructed convex bodies by Algorithm 1 can be
bounded by upper and lower probabilistic bounds (Section 2.4). When the input
is a zonotope we study other choices for the sequence of bodies, i.e.
sequences of balls and H-polytopes (Section 2.5).

###  2.1. Annealing schedule for convex bodies

Given a convex polytope $P$, an error parameter $\epsilon$ and
$r,\delta,\alpha$ s.t. $0<r+\delta<1$, the annealing schedule (Algorithm 2)
generates the sequence of convex bodies $C_{1}\supseteq\cdots\supseteq C_{m}$
defining $P_{i}=C_{i}\cap P,\ i=1,\dots,m$ and $P_{0}=P$. The main goal is to
restrict each ratio $r_{i}$ in the interval $[r,r+\delta]$ with high
probability. We call $r,\delta$ cooling parameters and $\alpha$ the
significance level (s.l.) of the schedule.

We introduce some notions from statistics needed to define two tests and refer
to [8] for details. Given $\nu$ observations from a r.v.
$X\sim\mathcal{N}(\mu,\sigma^{2})$ with unknown variance $\sigma^{2}$, the
(one tailed) t-test checks the null hypothesis that the population mean
exceeds a specified value $\mu_{0}$ using the statistic
$t=\frac{\bar{x}-\mu_{0}}{s/\sqrt{\nu}}\sim t_{\nu-1}$, where $\bar{x}$ is the
sample mean, $s$ the sample standard deviation and $t_{\nu-1}$ is the
t-student distribution with $\nu-1$ degrees of freedom. Given a s.l.
$\alpha>0$ we test the null hypothesis for the mean value of the population,
$H_{0}:\mu\leq\mu_{0}$ against $H_{1}:\mu>\mu_{0}$. We reject $H_{0}$ if,

$t\geq
t_{\nu-1,\alpha}\Rightarrow\bar{x}\geq\mu_{0}+t_{\nu-1,\alpha}{s}/{\sqrt{\nu}},$

which implies $\Pr(\text{reject }H_{0}\ |\ H_{0}\ \text{true})=\alpha$.
Otherwise we fail to reject $H_{0}$. The schedule algorithm uses the following
two statistical tests:  testL($P_{1},P_{2},r,\delta,\alpha,\nu,N$):
testR($P_{1},P_{2},r,\alpha,\nu,N$):
$H_{0}:\ \mbox{vol}(P_{2})/\mbox{vol}(P_{1})\geq r+\delta\hskip 122.34692pt\
H_{0}:\ \mbox{vol}(P_{2})/\mbox{vol}(P_{1})\leq r$
$H_{1}:\ \mbox{vol}(P_{2})/\mbox{vol}(P_{1})\leq r+\delta\hskip 122.34692pt\
H_{1}:\ \mbox{vol}(P_{2})/\mbox{vol}(P_{1})\geq r$
Successful if $H_{0}$ is rejected Successful if $H_{0}$ is rejected They are
used to restrict each $r_{i}$ to $[r,r+\delta]$. In the sequel we write
testL($P_{1},P_{2}$) and testR($P_{1},P_{2}$).

If we sample $N$ uniform points from a body $P_{i}$ then r.v. $X$ that counts
points in $P_{i+1}$, follows $X\sim b(N,r_{i})$, the binomial distribution and
$Y=X/N\sim\mathcal{N}(r_{i},r_{i}(1-r_{i})/N)$ is Gaussian.

######  Remark 2.1.

This normal approximation suffices when $N$ is large enough and we adopt the
well known rule of thumb to use it only if $Nr_{i}(1-r_{i})>10$.

Then each sample proportion that counts successes in $P_{i+1}$ over $N$ is an
unbiased estimator for the mean of $Y$, which is $r_{i}$. So if we sample $\nu
N$ points from $P_{i}$ and split the sample into $\nu$ sublists of length $N$,
the corresponding $\nu$ ratios are experimental values that follow
$\mathcal{N}(r_{i},r_{i}(1-r_{i})/N)$ and can be used to check both null
hypotheses in testL and testR. So, using the mean $\hat{\mu}$ of the ratios,
assuming

$r+\delta-t_{\nu-1,\alpha}\frac{s}{\sqrt{\nu}}\geq\hat{\mu}\geq
r+t_{\nu-1,\alpha}\frac{s}{\sqrt{\nu}},$

then $r_{i}$ is restricted to $[r,r+\delta]$ with high probability. Details
and bounds on the number of phases are given in Section 2.4.  Perform testR
and testL
Input: convex bodies $P_{1},\ P_{2}$, cooling parameters $r,\delta$,
significance level $\alpha$ and $\nu,N\in\mathbb{N}$
————————————————
Sample $\nu N$ uniform points from $P_{1}$
Partition $\nu N$ points to lists $S_{1},\dots,S_{\nu}$, each of length $N$
Compute ratios $\hat{r}_{i}=|\\{q\in P_{2}:\ q\in S_{i}\\}|/\,N$,
$i=1,\dots,\nu$
Compute the mean, $\hat{\mu}$, and standard deviation, $s$, of the $\nu$
ratios
if $\hat{\mu}\geq r+t_{\nu-1,\alpha}\frac{s}{\sqrt{\nu}}$ then testR holds,
otherwise testR fails
if $\hat{\mu}\leq r+\delta-t_{\nu-1,\alpha}\frac{s}{\sqrt{\nu}}$ then testL
holds, otherwise testL fails

Let us now describe the annealing schedule: Each body $C_{i}$ in
$C_{1}\supseteq\dots\supseteq C_{m}$ is a scalar multiple of a given body $C$.
When $C$ is the unit ball, the body used in each step is determined by a
radius. Since our algorithm does not use an inscribed ball, the initialization
step computes the body with minimum volume, denoted by $C^{\prime}$ or
$C_{m}$, s.t. $r_{m}\in[r,r+\delta]$ with high probability. The algorithm
employs $C^{\prime}$ to decide stopping at the $i$-th step; if the criterion
fails, the algorithm computes $P_{i+1}$ by a regular step.

Initialization step. The schedule is given convex body $C$, and an interval
$[q_{\min},q_{\max}]$. At initialization, one computes
$q\in[q_{\min},q_{\max}]$ s.t. both testL($qC,qC\cap P$) and testR($qC,qC\cap
P$) are successful. Now $q_{\min}C$ corresponds to a body that fails in testR
with unit probability, e.g. when $q_{\min}C\subseteq P$, while $q_{\max}C$ to
one for which testL fails with probability arbitrarily close to 1. The
computation of $q_{\min},q_{\max}$ is not trivial but in Section 3 we give
practical selections depending on body $C$ which are very efficient in
practice.

The algorithm performs binary search in $[q_{\min},q_{\max}]$. Let
$q=(q_{\min}+q_{\max})/2$ then:

  1. (1)

If testL($qC,qC\cap P$) succeeds and testR($qC,qC\cap P$) fails, continue to
the left-half of the interval.

  2. (2)

If testL($qC,qC\cap P$) fails and testR($qC,qC\cap P$) succeeds, continue to
the right-half of the interval.

  3. (3)

If both testL($qC,qC\cap P$) and testR($qC,qC\cap P$) succeed, stop and set
$C^{\prime}=qC$.

  4. (4)

If both testL($qC,qC\cap P$) and testR($qC,qC\cap P$) fail (contradiction)
then sample a new set of $\nu N$ uniform points from $qC$ and repeat both
tests.

Note that in each step of binary search, the schedule samples $\nu N$ points
from $qC$ to check both testL and testR. The output is $C^{\prime}$ which
shall be denoted by $C_{m}$ at termination.

Algorithm 2 AnnealingSchedule ($P$, body $C$, $r$, $\delta$, $\alpha$, $\nu$,
$N$, $q_{\min}$, $q_{\max}$)

$q_{1}=q_{\min}$, $q_{2}=q_{\max}$, $q=(q_{1}+q_{2})/2$

loop of the binary search for the initialization:

$C^{\prime}=qC$

if testR($C^{\prime}$,$C^{\prime}\cap P$) and
testL($C^{\prime}$,$C^{\prime}\cap P$) succeed stop loop

if testR($C^{\prime}$,$C^{\prime}\cap P$) succeeds and
testL($C^{\prime}$,$C^{\prime}\cap P$) fails set $q_{1}$=$q$,
$q$=$(q_{1}$+$q_{2})/2$

if testR($C^{\prime}$,$C^{\prime}\cap P$) fails and
testL($C^{\prime}$,$C^{\prime}\cap P$) succeeds set $q_{2}$=$q$,
$q$=$(q_{1}+q_{2})/2$

if testR($C^{\prime}$,$C^{\prime}\cap P$) and
testL($C^{\prime}$,$C^{\prime}\cap P$) fail sample $\nu N$ new uniform points
from $C^{\prime}$

end loop

set $P_{0}=P$, $i=0$, $q_{\min}=q$

loop for the definition of the sequence

$q_{1}=q_{\min}$, $q_{2}=q_{\max}$, $q=(q_{1}+q_{2})/2$

if testR($P_{i}$,$C^{\prime}\cap P$) succeeds set $m=i+1$,
$P_{m}=C^{\prime}\cap P$ and stop loop

loop of the binary search to define $P_{i+1}$

$P^{\prime}=qC\cap P$

if both testR($P_{i}$,$P^{\prime}$) and testL($P_{i}$,$P^{\prime}$) succeed
set $P_{i+1}=P^{\prime}$, $i=i+1$ stop loop

if testR($P_{i}$,$P^{\prime}$) succeeds and testL($P_{i}$,$P^{\prime}$) fails
set $q_{2}=q$, $q$=$(q_{1}$+$q_{2})/2$

if testR($P_{i}$,$P^{\prime}$) fails and testL($P_{i}$,$P^{\prime}$) succeeds
set $q_{1}=q$, $q$=$(q_{1}$+$q_{2})/2$

if testR($C^{\prime}$,$C^{\prime}\cap P$) and
testL($C^{\prime}$,$C^{\prime}\cap P$) fail sample $\nu N$ new uniform points
from $P_{i}$

end loop

$q_{\max}=q$

end loop

return $\\{P_{0},\dots,P_{m},C^{\prime}\\}$

Regular step. At step $i$, the algorithm determines $P_{i+1}$ by computing a
scaling factor of $C$ s.t. volume ratio $r_{i}\in[r,r+\delta]$ with high
probability. The schedule samples $\nu N$ points from $P_{i}$ and binary
searches for a $q_{i+1}$ in an updated interval $[q_{\min},q_{\max}]$ s.t.
both testL($P_{i},q_{i+1}C\cap P$) and testR($P_{i},q_{i+1}C\cap P$) are
successful. Then set $P_{i+1}=q_{i+1}C\cap P$. To update the interval,
Algorithm 2 uses the $q$ value of $C^{\prime}$ computed in the initialization
step as $q_{\min}$ and the $q$ value of $P_{i}$ computed in the previous step
as $q_{\max}$. The updated interval implies that $\mbox{vol}(P_{i+1})$ has to
lie between $\mbox{vol}(P_{i})$ and $\mbox{vol}(C^{\prime}\cap P)$.

Stopping criterion. The algorithm uses $C^{\prime}\cap P$ in the $i$-th step
for checking whether $\mbox{vol}(P_{i})/\mbox{vol}(C^{\prime}\cap P)\geq r$
with high probability, using only testR. Formally, $C^{\prime}\cap P$ is the
body with minimum volume in the sequence:

Stop in step $i$ if testR($P_{i},C^{\prime}\cap P$) holds. Then, set $m=i+1$,
and $P_{m}=C^{\prime}\cap P$.

To perform testR, $\nu N$ points are sampled from $P_{i}$. The schedule stops
when $\mbox{vol}(C^{\prime}\cap P)/\mbox{vol}(P_{i})$ is large enough
according to testR. It is clear that, at termination, $C_{m}=C^{\prime}$.
Otherwise, the algorithm determines the next convex body $P_{i+1}$ in a
regular step.

Termination. We demonstrate halting of Algorithm 2 for a given input polytope
$P$ and a set of parameters. Before stating the theorem, let us introduce the
notion of the power of a t-test: $pow=\Pr[\text{reject}\ H_{0}\ |\ H_{0}\
\text{false}]=1-\beta$. The power of a t-test cannot be usually calculated in
practice. It is well known that it depends on s.l. $\alpha$, sample size
$\nu$, and the magnitude of the effect on the mean value of the population.
For example, for testR, assuming $r_{i}=\theta>r$, we have

$pow(\theta)=\Pr\bigg{[}\frac{\hat{\mu}-r}{s/\sqrt{\nu}}>t_{\nu-1,\alpha}\
\bigg{|}\
r_{i}=\theta\bigg{]}=1-F^{-1}_{t_{\nu-1}}\bigg{(}t_{\nu-1,\alpha}-\frac{\theta-r}{s/\sqrt{\nu}}\bigg{)}=1-\beta,$

where $F^{-1}_{t_{\nu-1}}$ is the quantile function of t-student with $\nu-1$
degrees of freedom. A similar analysis for the power of testL is
straightforward.

In the t-tests of Algorithm 2 we might have some errors of type I or II and,
thus, binary search in intervals that do not contain values corresponding to
ratios in $[r,r+\delta]$. Therefore, there is a probability that Algorithm 2
fails to terminate. The following theorem states that this probability is
bounded by a constant when Algorithm 2 performs at least as many steps as the
minimum number required for it to terminate: this number is denoted by $Q$,
depends on the inputs, and occurs when there are no errors in the performed
t-tests.

######  Theorem 1.

Let Algorithm 2 perform some number $M\geq Q$ of steps. Let $\beta_{\max}$,
$\beta_{\min}$ be the maximum and the minimum among all the values of the
quantile function in the $M$ pairs of t-tests of testL and testR respectively.
Then Algorithm 2 terminates with constant probability, namely:

$\Pr[\text{Algorithm\leavevmode\nobreak\ \ref{schann} terminates}]\geq
1-2\frac{\alpha(1-\beta_{\min})+\beta_{\max}}{1-\alpha(1-\beta_{\min})+\beta_{\max}}-\frac{2\beta_{\max}-\beta_{\min}^{2}}{1-2\beta_{\max}-\beta_{\min}^{2}}.$

######  Proof 2.2.

Let $M\geq Q$. If Algorithm 2 fails to terminate after $M$ pairs of testL and
testR, then some type I or type II error occurred in the t-tests. An error of
type I occurs when the null Hypothesis is true and the test rejects it, while
type II occurs when the null Hypothesis is false and the test fails to reject
it. The respective probabilities are $\Pr[\text{reject }H_{0}\ |\ H_{0}\text{
true}]=\alpha$ and $\Pr[\text{fail to reject }H_{0}\ |\ H_{0}\text{
false}]=\beta$, which is a value of the quintile function of t-student. For
the latter probability we write $\beta_{L},\ \beta_{R}$ for testL and testR
respectively. If, for a pair of tests, both null hypotheses are false then an
error occurs with probability

$\displaystyle p_{1}=\Pr[\text{error occurs}\ |\ \text{both }H_{0}\text{
false}]$
$\displaystyle=\beta_{L}(1-\beta_{R})+\beta_{R}(1-\beta_{L})+\beta_{L}\beta_{R}$
$\displaystyle=\beta_{L}+\beta_{R}-\beta_{L}\beta_{R}\leq
2\beta_{\max}-\beta^{2}_{\min}$

Similarly,

$\displaystyle p_{2}=\Pr[\text{error occurs}\ |\ \text{{\tt testL}\
}H_{0}\text{ false }$ $\displaystyle\text{and {\tt testR}\ }H_{0}\text{
true}]=$
$\displaystyle=\alpha(1-\beta_{L})+\alpha\beta_{L}+\beta_{L}(1-\alpha)$
$\displaystyle=\alpha+\beta_{L}-\alpha\beta_{L}\leq\alpha+\beta_{\max}-\alpha\beta_{\min}$
$\displaystyle p_{3}=\Pr[\text{error occurs}\ |\ \text{{\tt testR}\
}H_{0}\text{ false }$ $\displaystyle\text{and {\tt testL}\ }H_{0}\text{
true}]=$
$\displaystyle=\alpha(1-\beta_{R})+\alpha\beta_{R}+\beta_{R}(1-\alpha)$
$\displaystyle=\alpha+\beta_{R}-\alpha\beta_{R}\leq\alpha+\beta_{\max}-\alpha\beta_{\min}$

Then,

$\displaystyle\Pr[\text{Algorithm\leavevmode\nobreak\ \ref{schann} fails to
terminate}]$
$\displaystyle\leq\sum_{i=1}^{M}p_{1}^{i}+\sum_{i=1}^{M}p_{2}^{i}+\sum_{i=1}^{M}p_{3}^{i}\leq\sum_{i=1}^{\infty}p_{1}^{i}+\sum_{i=1}^{\infty}p_{2}^{i}+\sum_{i=1}^{\infty}p_{3}^{i}$
$\displaystyle=\frac{1}{1-p_{1}}+\frac{1}{1-p_{2}}+\frac{1}{1-p_{3}}-3$
$\displaystyle\leq
2\frac{\alpha(1-\beta_{\min})+\beta_{\max}}{1-\alpha(1-\beta_{\min})+\beta_{\max}}+\frac{2\beta_{\max}-\beta_{\min}^{2}}{1-2\beta_{\max}-\beta_{\min}^{2}}.$

###  2.2. Empirical ratio estimation

As described in the previous section, annealing schedule returns $m$ bodies
intersecting $P$ that is we estimate $m+1$ ratios in total. In this section we
describe how this estimation is performed. First, we bound the error in the
estimation of each ratio in order to use it for the definition of the stopping
criterion. For each ratio $r_{i}$, we bound error by $\epsilon_{i}$ s.t.

${\sum_{i=0}^{m}\epsilon_{i}^{2}}=\epsilon^{2}$ (2)

then, from standard error propagation analysis, Equation (1) estimates
$\mbox{vol}(P)$ with error at most $\epsilon$. In section 3 we discuss
efficient error splitting in practical implementations.

Algorithm 3 EstimateRatio($P_{1},P_{2}$, $e$, $m$, $k$)

set $p=1-\sqrt[m+1]{3/4}$, $convergence=$false, $j=0$, $count\\_in=0$, sliding
window $W$ of length $k$

while $convergence==$ false do

$j=j+1$

Get a uniform point $q_{j}$ from $P_{1}$

if $q_{j}\in P_{2}$ then

$count\\_in=count\\_in+1$

end if

$\hat{r}=\frac{count\\_in}{j}$

if $W$ is full then remove $W_{1}$

store $\hat{r}$ to $W$

if $j>k$ then

set $s=std(W)$, $a=\hat{r}-z_{p/2}s$, $b=\hat{r}+z_{p/2}s$

if $(b-a)/a\leq e/2$ then

$convergence=$ true

end if

end if

end while

return $\hat{r}$

For fixed step $i$ of the schedule, and for each new sample point generated in
$P_{i}$, we update and keep the value of the $i$-th ratio. If we assume
uniform sampling then the number of points in $P_{i+1}$ follows the binomial
distribution $b(n,r_{i})$, where $n$ is the number of points we have generated
in $P_{i}$. Then a confidence interval of $r_{i}$ is given by

$\hat{r}\pm z_{\alpha/2}\sqrt{\frac{\hat{r}(1-\hat{r})}{n}}$ (3)

where $\hat{r}$ is the proportion of the number of points in $P_{i+1}$ over
$n$ and $z_{\alpha/2}$ is the $1-\alpha/2$ quintile of the Gaussian
distribution. Notice that while $n$ increases the interval is becoming tighter
around $\hat{r}$ so that a natural choice is to stop when
$z_{\alpha/2}\sqrt{\frac{\hat{r}(1-\hat{r})}{n}}\bigg{/}\bigg{(}\hat{r}-z_{\alpha/2}\sqrt{\frac{\hat{r}(1-\hat{r})}{n}}\bigg{)}\leq\epsilon_{i}$
and to sample $O(1/\epsilon_{i}^{2})$ points from $P_{i}$. Then Equation (1)
would estimate $\mbox{vol}(P)$ up to at most $\epsilon$ error with probability
$(1-\alpha)^{m+1}$. In practice we generate approximate uniform samples which
makes that criterion completely useless, because for small random walk steps
the number of points we obtain in $P_{i+1}$ do not follow the binomial
distribution and the criterion usually results to false positives.

Recall that the quantity $\sqrt{\hat{r}(1-\hat{r})/n}$ is an estimator of the
standard deviation of all the sample proportions of size $n$. In our empirical
criterion we replace that quantity with the standard deviation of a sliding
window, which stores consecutive estimators of $r_{i}$. In particular, we
store the last $k$ ratios in a queue called sliding window denoted by $W$ and
we say that $W$ has length $k$. We update $W$ each time a new sample point is
generated by inserting the new value of the $i$-th ratio and by popping out
the oldest ratio value in $W$. We stop sampling when $\hat{r}$ and the st.d.
of the values in $W$, called $s$, meet the criterion of Equation (4): they we
say they meet convergence. Clearly, for the first $k$ points sampled in
$P_{i}$, we do not check for convergence. Let $p=1-\sqrt[m+1]{3/4}$ and the
criterion is as follows:

$\text{If
}\frac{(\hat{r}+z_{p/2}s)-(\hat{r}-z_{p/2}s)}{\hat{r}-z_{p/2}s}=\frac{2z_{p/2}s}{\hat{r}-z_{p/2}s}\leq\frac{\epsilon_{i}}{2},\text{
then declare convergence.}$ (4)

The size $k$ of the sliding window is determined experimentally (Section 3).

###  2.3. Sampling

We use HnR with uniform target distribution for sampling from $P_{i}$ at step
$i$ of the annealing schedule or ratio estimation. With $\ell$ directed along
a random vector on the boundary of the $d$-dimensional unit hypersphere, we
have Random Direction HnR (RDHR). If $\ell$ is defined by a random vector
parallel to one of the axes, we have Coordinate Direction HnR (CDHR). If $P$
is given as a set of $m$ inequalities, RDHR costs $O(md)$ and CDHR $O(m)$ per
step. For zonotopes each step in both CDHR and RDHR solves the following LP to
compute one extreme point on $\ell\cap P$:

$\begin{split}\min\ \alpha,\ s.t.\quad p+\alpha
v=\sum_{i=1}^{k}\lambda_{i}g_{i}\quad-1\leq\lambda_{i}\leq 1.\end{split}$ (5)

For the second extreme point, keep the same constraints and minimize
$-\alpha$. This LP uses the basic feasible solution of the first one. For
V-polytopes we use the same LP with constraint $0\leq\lambda_{i}\leq 1$ and
$\sum_{i=1}^{v}\lambda_{i}=1$ while $g_{i}$ equals the $i$-th vertex. In
Section 3 we discuss practical choices of the walk step.

###  2.4. Number of phases

In this section we give probabilistic bounds for the number of phases. To do
this, we assume that i) polytope $P$ is sandwiched, $rB_{d}\subseteq
P\subseteq RB_{d}$, ii) we sample uniform points in each step of Algorithm 2
and iii) that Algorithm 2 terminates successfully. The construction in [17]
defines a sequence of convex bodies of length
$m=\lceil\lg(\mbox{vol}(RB_{d})/\mbox{vol}(rB_{d}))\rceil=\lceil
d\lg(R/r)\rceil$ and each ratio is restricted
$\frac{1}{2}\leq\mbox{vol}(K_{i})/\mbox{vol}(K_{i+1})\leq 1$. If $P$ is well-
rounded then $m=O(d\lg d)$. In our method we can use any convex body for the
sandwiching. Moreover we give a corollary which states which convex body
minimizes the number of phases for a given input $P$.

######  Proposition 2.

Given $P\subset{\mathbb{R}}^{d}$ and cooling parameters $r,\delta$ s.t.
$r+\delta<1/2$ and parameters $\alpha$, $N$, $\nu$, $q_{\min}$, $q_{\max}$ let
$m$ be the number of convex bodies in MMC returned by Algorithm 2. Then
$\Pr[m\leq\lceil d\lg(R/r)\rceil]\geq
2-\frac{1}{1-\gamma}=1-\frac{\gamma}{1-\gamma}$, where
$\gamma=\alpha(1-\beta_{\min})$ and $\beta_{\min}$ is the minimum among all
the values of the quantile function appearing in testR.

######  Proof 2.3.

If $m>\lceil d\lg(R/r)\rceil$ holds then $k\geq 1$ type I errors of testL
occurred in Algorithm 2, i.e. $H_{0}$ holds but the test rejects it, while
testR was successful, i.e. $H_{0}$ is false and the test rejects it. Type I
error occurs with probability $\alpha$ and the probability of the success of
testR is $1-\beta_{R}$. Let $\beta_{\min}$ be the minimum among the values of
the quantile function appearing in all instances of testR. Then,

$\displaystyle\Pr[m>\lceil
d\lg(R/r)\rceil]\leq\sum_{i=1}^{k}(a(1-\beta_{i}))^{i}\leq\sum_{i=1}^{\infty}(a(1-\beta_{\min}))^{i}=\frac{1}{1-\alpha(1-\beta_{\min})}-1\Rightarrow$
$\displaystyle\Pr[m\leq\lceil d\lg(R/r)\rceil]\geq
2-\frac{1}{1-\alpha(1-\beta_{\min})}=1-\frac{\gamma}{1-\gamma},\quad\mbox{where
}\gamma=\alpha(1-\beta_{\min}).$

Proposition 2 implies that the number of phases in [17] upper bounds the
number of phases in Algorithm 1 with high probability.

######  Proposition 3.

Given $P\subset{\mathbb{R}}^{d}$ and cooling parameters $r,\delta$ s.t.
$r+\delta<1/2$ and parameters $\alpha$, $N$, $\nu$, $q_{\min}$, $q_{\max}$ let
$m$ be the number of convex bodies in MMC returned by Algorithm 2. Then
$\Pr\bigg{[}\bigg{\lfloor}\log_{\frac{1}{r+\delta}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rfloor}\leq
m\leq\bigg{\lceil}\log_{\frac{1}{r}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rceil}\bigg{]}\geq
1-\frac{\gamma_{L}}{1-\gamma_{L}}-\frac{\gamma_{R}}{1-\gamma_{R}}$, where
$\gamma_{L}=\alpha(1-\beta^{L}_{\min})$,
$\gamma_{R}=\alpha(1-\beta^{R}_{\min})$ and $\beta^{R}_{\min},\
\beta^{L}_{\min}$ are the minimum among all the values of quantile function
appearing in testL and testR respectively.

######  Proof 2.4.

$\displaystyle\Pr\bigg{[}\bigg{\lfloor}\log_{\frac{1}{r+\delta}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rfloor}\leq
m\leq\bigg{\lceil}\log_{\frac{1}{r}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rceil}\bigg{]}$
$\displaystyle=1-\Pr\bigg{[}m>\bigg{\lceil}\log_{\frac{1}{r}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rceil}\
\text{or}\
m<\bigg{\lfloor}\log_{\frac{1}{r+\delta}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rfloor}\bigg{]}$
$\displaystyle=1-\Pr\bigg{[}m>\bigg{\lceil}\log_{\frac{1}{r}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rceil}\bigg{]}-\Pr\bigg{[}m<\bigg{\lfloor}\log_{\frac{1}{r+\delta}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rfloor}\bigg{]}$

Similar to the proof of Proposition 2,
$\Pr\bigg{[}m>\bigg{\lceil}\log_{\frac{1}{r}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rceil}\bigg{]}\leq\frac{\gamma_{R}}{1-\gamma_{R}}$,
where $\gamma_{R}=\alpha(1-\beta_{\min}^{R})$ and $\beta_{\min}^{R}$ is the
minimum among all the values of quantile function appearing in testR. If
$m<\bigg{\lfloor}\log_{\frac{1}{r+\delta}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rfloor}$
holds then $k\geq 1$ type I errors of testR occurred, while testL was
successful with probability $1-\beta$. Then,
$\Pr\bigg{[}m<\bigg{\lfloor}\log_{\frac{1}{r+\delta}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rfloor}\bigg{]}\leq\frac{\gamma_{L}}{1-\gamma_{L}}$,
where $\gamma_{L}=\alpha(1-\beta_{\min}^{L})$ and $\beta_{\min}^{L}$ is the
minimum among all the values of quantile function appearing in testL. Then

$\Pr\bigg{[}\bigg{\lfloor}\log_{\frac{1}{r+\delta}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rfloor}\leq
m\leq\bigg{\lceil}\log_{\frac{1}{r}}\bigg{(}\frac{\mbox{vol}(P)}{\mbox{vol}(P_{m})}\bigg{)}\bigg{\rceil}\bigg{]}\geq
1-\frac{\gamma_{L}}{1-\gamma_{L}}-\frac{\gamma_{R}}{1-\gamma_{R}}$

From Proposition 3 it is easily derived that the number of phases in Algorithm
1 is $m=O(\lg(\mbox{vol}(P)/\mbox{vol}(P_{m})))$ with high probability. If we
use balls in MMC and assume that $rB_{d}\subseteq P$ then
$m=O(\lg(\mbox{vol}(P)/\mbox{vol}(rB_{d})))$ which is always smaller, due to a
constant, than the number of phases in [17]. As we mentioned, there are well-
rounded convex bodies that our method improves this bound by $\lg d$, e.g. for
the cross polytope $m=O(d)$, as $\mbox{vol}(P)=2^{d}/d!$. We leave it as an
open question if this result can be extended to other convex bodies.
Proposition 3 also shows the importance of the type of body $C$.

######  Corollary 4.

Given $P$, the body $C$ that minimizes the number of phases of Algorithm 1,
with probability $p\geq
1-\frac{\gamma_{L}}{1-\gamma_{L}}-\frac{\gamma_{R}}{1-\gamma_{R}}$ for given
cooling parameters $r,\delta$ and s.l. $\alpha$, is the one that maximizes
$\mbox{vol}(C_{m}\cap P)$ in the annealing schedule of Algorithm 2.

###  2.5. Multiphase Monte Carlo for zonotopes

In this section we study different types of convex bodies used in the MMC
sequence to approximate the volume of a zonotope, i.e. the Minkowski sum of
segments. Let $k>d$ be the number of generators of zonotope
$P\subset{\mathbb{R}}^{d}$, and $G$ the $d\times k$ matrix they define. Note
$G^{T}G$ has $k-d$ zero eigenvalues; the corresponding eigenvectors form
$Q\in{\mathbb{R}}^{k\times(k-d)}$. The intersection of the hypercube
$[-1,1]^{k}$ with the $d$-dimensional affine subspace defined by $Q^{T}=0$
equals a $d$-dimensional polytope $C^{\prime}$ in ${\mathbb{R}}^{k}$. SVD
yields an orthonormal basis for the linear constraints, and its orthogonal
complement $W_{\perp}$:

$Q=USV^{T}=\left[\begin{array}[]{c}W\\\ W_{\perp}\\\
\end{array}\right]^{T}\left[\begin{array}[]{cc}S_{1}&0\\\ 0&0\\\
\end{array}\right]V^{T}$

Let $Ay\leq b_{0},\ A\in{\mathbb{R}}^{2k\times k}$ be a H-representation of
$[-1,1]^{k}$, then $AW^{T}_{\perp}x\leq b_{0},\
AW^{T}_{\perp}\in{\mathbb{R}}^{2k\times d}$ is an H-representation of a
$d$-dimensional polytope and $Mx\leq b_{0},\
M=AW^{T}_{\perp}(GW^{T}_{\perp})^{-1}\in{\mathbb{R}}^{2k\times d}$ is a
H-representation of the full-dimensional, centrally symmetric polytope
$C\subset P$ with $\leq 2k$ facets. Each $C_{i}$ in MMC arises from a parallel
shifting of the facets of $C$. This type of $C$ improves the schedule when
order is low, e.g. $\leq 4$.

Let $Mx\leq b_{0},\ x\in{\mathbb{R}}^{d}$ the H-representation of $C$. We
introduce a second improvement in the schedule: In each step $i-1$ we do not
compute a $q_{i}$ in an interval $[q_{\min},q_{\max}]$ as suggested in Section
2.1 but we consider two vectors
$b^{i}_{\min},b^{i}_{\max}\in{\mathbb{R}}^{2k}$ and compute with binary search
a $t_{i}\in[0,1]$ s.t. $b^{i}=b^{i}_{\min}+t_{i}(b^{i}_{\max}-b^{i}_{\min})$
and $C_{i}:=\\{x\ |\ Mx\leq b^{i}\\}$ results to successions in both testL and
testR. For a $\hat{t_{i}}$ if the testL is succeeded and testR is failed we
continue to the right-half of the interval and if testL is failed and testR is
succeeded we continue to the left-half; if both fail (contradiction) we sample
a new set of $\nu N$ points and repeat both tests.

######  Proposition 5.

Let $A\in{\mathbb{R}}^{q\times d}$ the matrix that contains row-wise the
normals of the facets of a convex polytope $C\subset P$, where $P$ is a
$d$-dimensional zonotope. Let $G\in{\mathbb{R}}^{d\times k}$ the generators’
matrix and $r_{l}=AG_{l}$, the $l_{th}$ row of $AG$ and
$b_{\max}\in{\mathbb{R}}^{m}$ s.t. $b_{\max}^{l}=\sum_{j=1}^{d}|r_{lj}|,\
l=1,\dots q$, where $b_{\max}^{l}$ is the $l_{th}$ coefficient of $b_{\max}$.
If $C_{\max}:=\\{x\ |\ Ax\leq b_{\max}\\}$ then $C_{\max}\supseteq P$ holds.
Moreover for every facet $f_{j}$ of $C_{\max}$, $f_{j}\cap P\neq\emptyset,\
j=1,\dots,q$ holds.

######  Proof 2.5.

Let $G\in{\mathbb{R}}^{d\times k}$ the matrix of zonotope $P$,
$H:=\\{c^{T}x\leq z_{0}\\},\ c\in{\mathbb{R}}^{d}$ a halfspace intersecting
$P$. Let $z_{\max}=\max\\{c^{T}y\ |\ y\in P\\}$. For all $y_{0}\in P$ there is
a $\lambda\in[-1,1]^{k}$ s.t. $y_{0}=G\lambda$, so
$z_{\max}=\max\\{c^{T}G\lambda\ |\ \lambda\in[-1,1]^{k}\\}$. The $\lambda$
that gives the maximum inner product with a vector $v\in{\mathbb{R}}^{k}$ is
the following,

$\lambda_{i}=\left\\{\begin{array}[]{ll}1&v_{i}\geq 0\\\ -1&v_{i}<0\\\
\end{array}\right.,\ i=1,\dots,k$

So $z_{\max}=\sum_{i=1}^{k}|c^{T}G_{i}|$. Then for the halfspace
$H^{\prime}:=\\{c^{T}x\leq z_{\max}\\}$, $H^{\prime}\supseteq P$ holds.
Moreover $H^{\prime}\cap P\neq\emptyset$. We apply the same for every facet of
$C\subset P$.

Following this proposition, we use $b_{0},b_{\max}$ to compute $C^{\prime}$ in
the initialization step of the schedule obtaining a vector $b^{\prime}$. Then
in each step, $b_{\min}=b^{\prime}$ and $b_{\max}=b_{i-1}$. When $C^{\prime}$
is an H-polytope we estimate $\mbox{vol}(C_{m})$ using Algorithm 1 by using
balls in MMC and, moreover, we sample from $C_{m}$ with HnR to estimate ratio
$r_{m}$.

##  3\. Implementation and experiments

In this section we discuss our implementation and setup the parameters of the
algorithm. We perform extended experiments analyzing various aspects of our
method such as practical complexity and how it is affected by the bodies used
in MMC. Finally, we apply our software to test the quality of approximation of
various methods for low order reduction of zonotopes.

We use the eigen library (http://eigen.tuxfamily.org) for linear algebra and
lpSolve (http://lpsolve.sourceforge.net) for solving LPs. All experiments have
been performed on a personal computer with Intel® Core™ i7-6700 3.40GHz × 8
CPU and 32GB RAM. All runtimes reported in the plots and tables are averaged
over 10 runs unless otherwise stated. When we use balls in MMC we denote by
CoolingBall the implementation of Algorithm 1, and by CoolingHpoly when we use
the H-polytope from Section 2.5. We denote by CoolingGaussian the
implementation of [7] and by SeqOfBalls that of [11].

Polytope database. To perform our experiments is vital to adopt a polytope
database with polytopes that explore average cases as well as corners of our
method. We use

  * •

cube-$d$: $\\{x=(x_{1},\dots,x_{d})\,|\,x_{i}\leq
1,x_{i}\geq-1,x_{i}\in{\mathbb{R}}\text{ for all }i=1,\dots,d\\}$,

  * •

cross-$d$: cross polytope, the dual of cube, i.e.
$\mbox{conv}(\\{-e_{i},e_{i},\,i=1,\dots,d\\})$,

  * •

$\Delta$-$d$: the $d$-dimensional simplex $\mbox{conv}(\\{e_{i},\text{ for
}i=1,\dots,d\\})$

  * •

rh-$d$-$m$: polytopes constructed by choosing $m$ random hyperplanes tangent
to the sphere,

  * •

rv-$d$-$n$: dual to rh-$d$-$m$, i.e. polytopes with $n$ vertices randomly
distributed on the sphere,

  * •

z-$d$-$k$: $\sum_{s\in S}s$; choose a random direction for each segment $s\in
S$ and pick a random length in the interval $[0,\sqrt{d}]$

In the experiments of this paper we do not apply to the polytope any rounding
step before the volume computation. However, in the case of random
V-polytopes, we compute the minimum volume enclosing ellipsoid of the vertices
and then we apply to the polytope a linear map that maps the ellipsoid to the
unit ball, to test our method for random polytopes that are near to well-
rounded position (Table 1).

Annealing Schedule. For the cooling parameters we set $r=0.1$ and
$\delta=0.05$ in order to define the next convex body in MMC with about $10\%$
of the volume of the previous body, and we choose the significance level to be
$\alpha=0.10$. A smaller $\alpha$ can be chosen for a tighter test around
$r+\delta$. We set the number of points that are generated from $P_{i}$ in
each step to be $\nu N=1200+2d^{2}$ and $\nu=10$. The value of $\nu N$ was
experimentally determined; HnR generates approximate samples, so we set $\nu
N=O(d^{2})$ to obtain, in practice, the normal approximation that Remark 2.1
implies. The choice of these parameters might be improved in the future. In
the $i$-th step of the schedule we sample from $P_{i}$ and then check if the
stopping criterion holds. If it fails we binary search reusing the sample.
Hence we sample only once a set of $\nu N$ points from $P_{i}$ in each step of
Algorithm 2.

When $C$ is a ball we can sample almost perfectly uniform points from it (no
need of random walks). Following Remark 2.1, the normal approximation suffices
when $r_{i}=0.1$ and $N=112$, thus we set $\nu N=1200=O(1)$, $\nu=10$. When
$C$ is the polytope of Section 2.5 we set a large step for HnR, namely
$10+2d$, to get more accurate samples, and set $\nu N=1200$ as we can assume
that the sample is close to uniform. For the rest of the steps in the schedule
we use HnR to sample from $P_{i}$ where the value $1200+2d^{2}$ was picked
experimentally so that we have stable cooling as the dimension grows.

When using balls, let $C$ be the unit ball, then $q_{\min},q_{\max}$ are
minimum and maximum radii in computing $C^{\prime}$. We set $q_{\min}=0$ and
sample $1200+2d^{2}$ points to set $q_{\max}$ s.t. all points belong to

Experimental results for V-polytopes. $P$ Vol m steps error time ex. Vol ex.
time cross-$20$ 4.54e-13 1 3.42e+03 0.05 2.81 $--$ $--$ cross-$40$ 1.22e-36 1
8.23e+03 0.09 14.5 $--$ $--$ cross-$60$ 1.27e-64 1 20.6e+03 0.08 60.1 $--$
$--$ cross-$80$ 1.62e-95 1 29.3e+03 0.04 122 $--$ $--$ cross-$100$ 1.43e-128 2
94.2e+03 0.05 406 $--$ $--$ $\Delta$-$20$ 4.24e-19 2 1.49e+04 0.03 6.970
4.11e-19 0.07 $\Delta$-$40$ 1.32e-48 5 21.2e+04 0.08 210.3 1.22e-48 0.008
$\Delta$-$60$ 1.08e-82 10 77.4e+04 0.1 1442 1.203-82 0.02 $\Delta$-$80$
1.30e-119 13 187e+04 0.07 6731 1.39e-119 0.07 cube-$10$ 1052.4 1 1851 0.03
54.39 $--$ $--$ cube-$11$ 1930.2 1 1935 0.06 155.9 $--$ $--$ cube-$12$ 4240.6
1 2017 0.04 567.0 $--$ $--$ cube-$13$ 7538.2 1 2127 0.08 2937 $--$ $--$
rv-$10$-$20$ 9.11e-06 1 0.185e+04 0.08 1.57 9.87e-06 0.015 rv-$10$-$40$
5.04e-04 1 0.185e+04 0.09 1.92 4.64e-04 0.46 rv-$10$-$80$ 3.74e-03 1 0.185e+04
0.08 3.35 3.46e-03 6.8 rv-$10$-$160$ 1.59e-02 1 0.140e+04 0.06 5.91 1.50e-03
59 rv-$15$-$30$ 2.73e-10 1 0.235e+04 0.02 3.46 2.79e-10 2.1 rv-$15$-$60$
4.41e-08 1 0.235e+04 $??$ 6.46 $--$ $--$ rv-$20$-$40$ 1.27e-15 1 0.305e+04
$??$ 6.24 $--$ $--$ rv-$30$-$60$ 6.06e-28 1 0.914e+04 $??$ 25.7 $--$ $--$
rv-$40$-$80$ 9.37e-42 2 3.02e+04 $??$ 148 $--$ $--$ rv-$50$-$100$ 8.34e-57 2
4.29e+04 $??$ 325 $--$ $--$ rv-$20$-$200$ 2.57e-10 1 0.305e+04 $??$ 35.3 $--$
$--$ rv-$20$-$2000$ 2.89e-07 1 0.305e+04 $??$ 457 $--$ $--$ rv-$80$-$160$
5.84e-106 3 11.3e+04 $??$ 325e+01 $--$ $--$ rv-$100$-$200$ 1.08e-141 4
24.5e+04 $??$ 133e+02 $--$ $--$

Table 1.  Vol: the average estimated volume; m: the maximum number of phases
in MMC; steps: the average number of steps of CoolingBall; $error$ the error
of each computation; time: the average time in seconds for CoolingBall; ex.
Vol: the exact volume; ex. time: the time in seconds for the exact volume
computation i.e. qhull in R (package geometry). The $--$ implies that the
execution failed due to memory issues or exceeded $1$ hr. We set the requested
error $\epsilon=0.1$ for all the above.

$q_{\max}C$. For low order zonotopes we use the polytope in Section 2.5 and
follow those steps.

Error spliting. We do not split the error equally to all ratios, but set
$\epsilon_{m}=\epsilon/2\sqrt{m+1}$ in order to be the smallest one as the
ratio $\mbox{vol}(P_{m})/\mbox{vol}(C_{m})$ converges faster than the other
ratios in practice. The latter occurs because sampling from $C_{m}$ is usually
faster and more accurate, e.g. when $C_{m}$ is a ball. Then we split
$\epsilon^{\prime}=\epsilon\sqrt{4(m+1)-1}/2\sqrt{m+1}$ equally to the
remaining ratios, i.e. $\epsilon_{i}=\epsilon^{\prime}/\sqrt{m},\
i=0,\dots,m-1$ so that Equation (2) holds. As we mentioned in section 2.5, if
$C_{m}$ is an H-polytope and $P$ a zonotope we estimate $\mbox{vol}(C_{m})$
calling Algorithm 1 and by using balls in MMC. For the computation of
$\mbox{vol}(C_{m})$ we set $e^{\prime\prime}=\epsilon/2\sqrt{m+1}$ and then we
equally split $\epsilon^{\prime}=\epsilon\sqrt{2m+1}/\sqrt{2m+2}$ to the $m+1$
ratios respecting Equation (2).

Sliding window. For the ratio estimation we reuse the points from annealing
schedule but we use only the value of the ratio between the number of
successes over $\nu N$. For the sliding window length, set $k=2d^{2}+250$
since our experiments on the error show this choice offers stability. In
Figure 2 we see that, for unit cubes, the error of Algorithm 1 when
$\epsilon=0.1$ is around that value and converges to $0.057$ while we increase
the number of experiments and take the average error for the $20$-dimensional
unit cube. If we set $k=O(d)$ the error of our method seems to exceed
$\epsilon$ while the dimension increases, but we believe further improvements
can be made. For each new generated point we update the mean and the variance
of the sliding window in $O(1)$ instead of $O(d)$: Let $\hat{\mu}$ be the mean
of the sliding window, then we write the sample variance,

$\frac{1}{k}\sum_{i=1}^{k}(\hat{r}_{i}-\hat{\mu})^{2}=\frac{1}{k}(\sum_{i=1}^{k}\hat{r}_{i}^{2}-2\hat{\mu}\sum_{i=1}^{k}\hat{r}_{i}+k\hat{\mu}^{2})$

We store the sum of the window’s ratios and the sum of the squared ratios. For
each new generated point we obtain an updated ratio and the oldest ratio is
popped out. We use both values to update both $\sum_{i=1}^{k}\hat{r}_{i}^{2}$
and $\sum_{i=1}^{k}\hat{r}_{i}$ and to compute the updated mean value and
st.d. of the window. In practice the length of the window is always small
enough to achieve numerical robustness.

Figure 2. Left: the errors of CoolingBall using RDHR and CDHR for unit cubes
in H-representation, $d=5,10,\dots,100$. Right: We run CoolingBall multiple
times. The plot shows the errors of CoolingBall averaged in each iteration
using RDHR and CDHR for the $20$-dimensional unit cube in H-representation as
well; RDHR converges to $0.057$ and CDHR converges to $0.129$.

Sampling by HnR. Theoretical bounds for the mixing time are too pessimistic,
thus our approach is more aggressive and we set the step equals to one.
Moreover we do not compute a warm start for the random walk; we always start
HnR from the feasible point we compute in $P$ when we estimate ratio $r_{i}$.
This approach is faster because the sliding window can handle efficiently
these disadvantages. For zonotopes and V-polytopes we use RDHR since
coordinate directions do not have any advantage. When $P$ is a H-polytope we
use CDHR but our experiments show that RDHR is statistically more stable as
you can see in Figure 2.

Body selection for MMC. Algorithm 1 needs a point in $P$ in order to start the
sampling. If $P$ is an H-polytope we compute the Chebychev center [4]. When
$P$ is a zonotope we use the origin. For V-polytopes, we compute an
approximation of the minimum volume enclosing ellipsoid and use the center of
the ellipsoid as a feasible point. In general, for zonotopes, one selects
convex body $C$ as follows: Estimate $\mbox{vol}(C^{\prime}\cap P)$ for $C$
being a ball and the H-polytope in Section 2.5; according to Corollary 4 pick
the one that maximizes this volume.

Figure 3.  Left: Number of bodies in MMC using balls. Right: Number of bodies
in MMC using centrally symmetric H-polytopes. We set the number of generators
$k=d+1,\ \lfloor 1.1d\rfloor,\ \lfloor 1.2d\rfloor,\ \lfloor 1.3d\rfloor,\
\lfloor 1.4d\rfloor,\ \lfloor 1.5d\rfloor,\ 2d,\ 3d,\ 4d$ and we generate a
random zonotope for each experiment.

Figure 3 shows that for random zonotopes of order $\leq 4$, if we use the
H-polytope, the number of bodies used in MMC is smaller for all pairs $(d,k)$
compare to the case of using balls. Notice that, when we use balls in MMC, the
number of phases is decreasing for a constant $d$ while $k$ increases. Table 2
shows that the number of phases for high-order zonotopes is $m=1$, whereas for
low-order zonotopes the H-polytope we defined reduces significantly the number
of phases and runtime. The maximum number of phases for zonotopes (up to what
our software computes in $<10$hr) is $m\leq 3$.

Experimental results for zonotopes. z-$d$-$k$ $Body$ $order$ $Vol$ $m$ $steps$
$time$ z-$5$-$500$ Ball 100 4.63e+13 1 0.1250e+04 22.26 z-$10$-$1000$ Ball 100
2.62e+29 1 0.1400e+04 130.1 z-$15$-$1500$ Ball 100 5.00e+45 1 0.1650e+04 506.1
z-$20$-$2000$ Ball 100 2.79e+62 1 0.2000e+04 1428 z-$50$-$65$ Hpoly 1.3
1.42e+62 1 1.487e+04 173.9 z-$60$-$78$ Hpoly 1.3 2.88e+75 2 3.800e+04 529.5
z-$70$-$91$ Hpoly 1.3 7.64e+90 2 5.470e+04 1067 z-$80$-$104$ Hpoly 1.3
1.10e+107 2 8.685e+04 2277 z-$90$-$117$ Hpoly 1.3 1.76e+122 3 10.17e+04 2810
z-$100$-$130$ Hpoly 1.3 1.37e+138 3 17.19e+04 6073 z-$40$-$60$ Hpoly 1.5
4.23e+51 1 0.7851e+04 104.5 z-$45$-$67$ Hpoly 1.5 4.17e+58 1 0.9551e+04 135.0
z-$50$-$75$ Hpoly 1.5 2.96e+66 2 1.615e+04 253.6 z-$55$-$82$ Hpoly 1.5
9.83e+74 2 3.469e+04 541.2 z-$60$-$90$ Hpoly 1.5 5.81e+82 2 5.355e+04 943.9
z-$65$-$97$ Hpoly 1.5 1.40e+90 2 4.869e+04 1296 z-$70$-$105$ Hpoly 1.5
8.66e+98 2 5.495e+04 1352 z-$75$-$112$ Hpoly 1.5 1.45e+107 3 4.915e+04 1632
z-$80$-$120$ Hpoly 1.5 8.48e+114 3 12.35e+04 4180 z-$85$-$127$ Hpoly 1.5
1.04e+123 3 12.09e+04 4507 z-$90$-$135$ Hpoly 1.5 5.11e+131 3 14.44e+04 6710
z-$95$-$142$ Hpoly 1.5 2.21e+141 3 15.55e+04 8524 z-$100$-$150$ Hpoly 1.5
2.32+149 3 15.43e+04 10060 z-$70$-$140$ Hpoly 2 8.71e+111 2 5.059e+04 2695
z-$75$-$150$ Hpoly 2 6.18e+121 3 11.12e+04 4036 z-$80$-$160$ Hpoly 2 2.01e+131
3 11.31e+04 5356 z-$100$-$200$ Hpoly 2 5.27e+167 3 15.25e+04 34110

Table 2.  $Body$ stands for the type of body in MMC; $order=k/d$, $Vol$ the
average of volumes by CoolingHpoly over $10$ runs; $m$ the maximum number of
phases; $steps$ is the number of steps; $time$ is average time in seconds;
$\epsilon=0.1$ except for z-$80$-$160$ and z-$100$-$200$ where $\epsilon=0.2$.

Figure 4. $\sqrt{\\#steps}$ as a function of $d$ for unit cubes (left) and
Birkhoff polytopes (right). Left: the dimension $d=5,10,\dots,200$. Right:
$n=7,4,\dots 15$ for Birkhoff polytopes $B_{n}$. We use balls in MMC.

Figure 5. Left: Number of balls in MMC for cross polytopes. Right: The
$\sqrt{\\#steps}$ for cross polytopes. In both plots, $d=5,10,\dots,100$.

Experimental complexity. To study practical complexity we experimentally
correlate the total number of HnR steps $n$ with $d$. Then, runtime grows
accordingly to the cost of each step of HnR. We test unit H-cubes, Birkhoff
H-polytopes, and cross V-polytopes, and found a linear relation between
$\sqrt{n}$ and $d$ (Figures 4 and 5). The Pearson correlation coefficients
between $\sqrt{n}$ and $d$ are $p_{1}=0.991$, $p_{2}=0.994$ and $p_{3}=0.957$
respectively. These coefficients imply almost perfect positive linear
correlation for all cases. We conclude that our method needs $O(d^{2})$ steps
for these polytopes. In Figure 5 (left) notice that for $d\leq 85$ our method
defines only one ball in MMC for the cross polytope while for $d=100$ it
defines two balls, thus verifying the observation in Section 2.4 for the
particular polytope. In Table 1 notice that the maximum number of phases for
random V-polytopes after rounding is $\leq 4$, showing that our method defines
just a few rejection sampling steps for random convex polytopes that are near
to well-rounded position for $d\leq 100$.

Comparison with other implementations. We compare against the (only) two
available implementations for high dimensional volume estimation namely, the
matlab implementation of CoolingGaussian and the C++ implementation of
SeqOfBalls. In Table 3 we use the matlab implementation CoolingGaussian for
the number of steps. For more fair time comparisons we implement in C++ the
method of CoolingGaussian. Interestingly, we found our implementation to be
around 10 times faster; then we compare the performance of the two C++
implementations in Table 3. The experimental results in Table 3 show that
CoolingBall is faster than both CoolingGaussian and SeqOfBalls for $d\leq 100$
and for the $120$ dimensional hypercube. Figure 6 (left plot) shows the number
of steps of CoolingBall, CoolingGaussian and of SeqOfBalls on unit H-cubes. In
[7] they show that the number of steps CoolingGaussian needs is $O(d^{2})$
even for unit cubes. Our experiments confirm this complexity and show that
CoolingBall is faster than CoolingGaussian for $d\leq 100$ and asymptotically
faster than SeqOfBalls.

Figure 6. Left: the number of steps for for unit cubes in H-representation,
$d=5,10,\dots,100$. Right: the number of steps for random zonotopes of order
$2$, $d=5,10,\dots,80$ . In both plots we use $log_{10}$ scale for the
$y$-axis

To compare CoolingGaussian and SeqOfBalls with Algorithm 1 for zonotopes we
define an inscribed ball and for CoolingGaussian we give the known upper bound
on the number of facets. We compute $r=\max\\{\gamma\ |\ \gamma e_{i}\in
P\\}$, then the ball $B(\mathbf{0},r/\sqrt{d})\subseteq P$. In Figure 6 (right
plot) we compare the number of steps for random zonotopes of order 2, between
CoolingHpoly and both

Experimental results for H-polytopes. $P$ cB Vol m cB steps cB error cB time
cG Vol. cG steps cG error cG time SoB Vol. SoB steps SoB error SoB time
cube-$20$ 1.06e+05 1 0.0419e+05 0.01 0.10 9.94e+05 0.997e+05 0.05 0.030
1.05e+06 0.363e+07 0.002 1.32 cube-$40$ 9.91e+11 2 0.347e+05 0.09 0.20
1.13e+12 5.31e+05 0.03 0.24 1.15e+12 2.84e+07 0.05 16.5 cube-$60$ 1.23e+18 5
2.05e+05 0.07 0.63 1.24e+18 12.8e+05 0.07 0.70 1.19e+18 9.42e+07 0.03 58.6
cube-$80$ 1.27e+24 7 5.37e+05 0.05 1.6 1.16e+24 26.7e+05 0.04 4.1 1.22e+24
22.2e+07 0.02 156 cube-$100$ 1.23e+30 9 10.8e+05 0.03 3.7 1.34e+30 45.5e+05
0.05 9.4 1.17e+30 43.3e+07 0.08 301 cube-$120$ 1.25e+36 12 20.1e+05 0.06 7.1
1.19e+36 83.8e+05 0.1 14 1.37e+36 75.1e+07 0.03 897 $\Delta$-$60$ 1.07e-82 11
10.7e+05 0.02 0.99 1.26e-82 22.1e+05 0.04 2.9 1.15e-82 21.6e+07 0.04 70.2
$\Delta$-$80$ 5.46e-113 16 42.3e+05 $>1$ 2.9 1.39e-119 45.9e+05 0.003 6.4
1.46e-119 0.876e+07 0.04 173 ${}^{*}\Delta$-$80$ 1.31e-119 16 29.0e+05 0.06 35
1.39e-119 45.9e+05 0.003 6.4 1.46e-119 0.876e+07 0.04 173 $\Delta$-$100$
2.65e-149 21 92.4e+05 $>1$ 7.4 1.13e-158 79.9e+05 0.05 16 1.06e-158 1.35+07
0.01 528 ${}^{*}\Delta$-$100$ 9.80e-159 21 55.2e+05 0.06 90 1.13e-158 79.9e+05
0.05 16 1.06e-158 1.35+07 0.01 528 rh-$20$-$60$ 2.66e+27 3 0.259e+05 ?? 0.074
2.35e+27 1.83e+05 ?? 0.16 2.44e+27 0.449e+07 ?? 3.24 rh-$40$-$160$ 3.82e+53 5
1.82e+05 ?? 0.47 2.66e+53 7.51e+05 ?? 1.8 2.52e+53 3.14e+07 ?? 5.54
∗rh-$40$-$160$ 2.26e+53 5 2.07e+05 ?? 2.5 2.66e+53 7.51e+05 ?? 1.8 2.52e+53
3.14e+07 ?? 5.54 rh-$60$-$300$ 5.68e+77 7 5.14e+05 ?? 2.1 5.68e+77 16.8e+05 ??
6.8 6.10e+77 9.48e+07 ?? 305 rh-$100$-$600$ 7.68e+126 1 19.0e+05 ?? 13
7.40e+126 59.0e+05 ?? 51 – – – $>1$ hr

Table 3.  $P$ the type of polytope; $m$ for the maximum number of phases in
MMC among $10$ experiments for CoolingBall. cb denotes CoolingBall, cg denotes
CoolingGaussian and SoB denotes SeqOfBalls. $Vol$ the averaged estimated
volume, $steps$ the averaged number of steps, $error$ the error of each
method, $time$ the averaged time in seconds. The ∗ denotes that we use RDHR
where CDHR for CoolingBall fails. We set the requested error $\epsilon=0.1$
for all the above.

CoolingGaussian and SeqOfBalls. CoolingGaussian fails to estimate volumes for
$d>15$ as the upper bound for the number of facets is the bottleneck for this
implementation while SeqOfBalls takes $>1$ hr for $d>15$. Moreover
CoolingHpoly needs a smaller number of steps for $d=80$ than CoolingGaussian
for a random zonotope in $d=15$. In Table 1 we give volume estimations for
V-polytopes and exact computation using qhull. Notice CoolingBall scales
efficiently up to $d=100$ for random V-polytopes after the rounding step we
described. Moreover, it is very efficient for the case of cross polytope while
qhull fails for $d\geq 20$. In addition, CoolingBall estimates the volume of
the hypercube for $d\leq 13$ and takes $\geq 1$ hr for larger dimensions as
the number of vertices grows exponentially in $d$. Finally, notice that qhull
is faster for simplices as the computation of the convex hull is very fast.

Summarizing our implementations of Algorithm 1 for V-polytopes and zonotopes
outperforms both CoolingGaussian and SeqOfBalls and, additionally, perform
volume computations which were intractable until now. CoolingBall is
asymptotically better than SeqOfBalls and faster for $d\leq 100$ than
CoolingGaussian for H-polytopes.

###  3.1. Application: test zonotope approximations

We propose an efficient algorithm for evaluating over-approximation of a
zonotope $P$. Zonotopes are critical in applications such as autonomous
driving [1] or human-robot collaboration [20]. Algorithm complexity strongly
depends on the order of the encountered zonotopes. Thus, a practical solution
is to over-approximate $P$, as tight as possible, with another zonotope
$P_{red}\supseteq P$ of smaller order. A good measure of the approximation’s
quality (fitness) is

$R=(\mbox{vol}(P_{red})/\mbox{vol}(P))^{1/d}.$ (6)

This reduces to volume computation. In [14] they compute volumes exactly and
deterministically, therefore they cannot compute the quality of approximation
for $d>10$.

Here, we employ our software to test the quality of such approximations.
Methods that are able to scale for $d\geq 20$ are, primarily, Principal
Component Analysis (PCA) and the BOX method [14]. Both adopt similar
approximations and are of comparable reliability; here we focus on PCA.  PCA
(zonotope $P$ with generators’ matrix $G\in{\mathbb{R}}^{d\times k}$)
$X=[G|-G]^{T}$
$USV^{T}=\text{SVD}(X^{T}X)$
Return: $G_{red}=U\cdot IH(U^{T}G)$ Notice that $G_{red}=U\cdot IH(U^{T}G)$ is
square and generates $P_{red}$; the $IH(\cdot)$ is the “interval hull” from
[15]. Over-approximation can be seen as a reduction problem, so that the
covariance among the $d$ generators of $P_{red}$ must be null.

Table 4 shows experimental results for zonotopes up to $d=30$ of order up to
$15$. We use balls in MMC (CoolingBall); $\mbox{vol}(P_{red})$ is obtained
exactly by computing one determinant. For PCA over-approximations we show that
$R$ increases as $d$ grows but the same does not occur for fixed $d$ as order
increases. This is probably the first time practical volume estimation is used
to test approximation methods.

Testing ratios of PCA approximations. z-$d$-$k$ Order $\mbox{vol}(P)$ Time
$\mbox{vol}(P_{red})$ $R$ z-$10$-$50$ $5$ 1.19e+16 2.10 7.09e+18 1.90
z-$10$-$100$ $10$ 1.93e+19 3.11 6.51e+21 1.79 z-$10$-$150$ $15$ 1.24e+21 6.84
5.41e+23 1.84 z-$10$-$500$ $50$ 1.90+26 37.5 9.87e+28 1.86 z-$15$-$65$ $5$
4.21e+25 7.11 4.72e+30 2.17 z-$15$-$150$ $10$ 8.29e+30 14.8 3.11e+35 2.02
z-$15$-$225$ $15$ 3.46e+33 26.3 2.19e+38 2.09 z-$15$-$750$ $50$ 1.36e+41 113
1.39e+46 2.16 z-$20$-$200$ $10$ 5.25e+42 36.7 8.58e+49 2.29 z-$20$-$300$ $15$
1.78e+46 66.4 2.55e+53 2.28 z-$30$-$300$ $10$ 1.46e+66 180 1.15e+80 2.91
z-$30$-$450$ $15$ 6.57e+71 359 2.62e+85 2.84

Table 4.  $\mbox{vol}(P)$ is the estimated volume of zonotope z-$d$-$k$ with
CoolingBall, in Time (sec), with requested error $\epsilon=0.1$;
$\mbox{vol}(P_{red})$ is the volume of the over-approximation; $R$ the ratio
of fitness in Equation (6).

##  4\. Acknowledgements

I.Z.E. and A.C. are partially supported by the European Union’s Horizon 2020
research and innovation programme under grant agreement No 734242 (Project
LAMBDA). We would like to thank Prof. Althoff for discussions on zonotope
approximation methods (Section 3.1).

## References

  * [1] M. Althoff and J. M. Dolan.  Online verification of automated road vehicles using reachability analysis.  IEEE Transactions on Robotics, 30(4):903–918, Aug 2014.
  * [2] U. Betke and M. Henk.  Approximating the volume of convex bodies.  Discrete & Computational Geometry, 10:15–21, 1993.
  * [3] J. Bourgain and J. Lindenstrauss.  Approximationg the ball by a Minkowski sum of segments with equal length.  Discrete Computational Geometry, 9:131–144, 1993.
  * [4] Stephen Boyd and Lieven Vandenberghe.  Convex Optimization.  Cambridge University Press, New York, NY, USA, 2004.
  * [5] A. Chevallier, S. Pion, and F. Cazals.  Hamiltonian Monte Carlo with boundary reflections, and application to polytope volume calculations.  Research Report RR-9222, INRIA Sophia Antipolis, France, 2018.
  * [6] B. Cousins and S. Vempala.  Bypassing KLS: Gaussian cooling and an ${O}^{*}(n^{3})$ volume algorithm.  In Proc. ACM STOC, pages 539–548, 2015.
  * [7] B. Cousins and S. Vempala.  A practical volume algorithm.  Mathematical Programming Computation, 8, June 2016.
  * [8] H. Cramer.  Mathematical methods of statistics.  Princeton University Press, 1946.
  * [9] M. Dyer, A. Frieze, and R. Kannan.  A random polynomial-time algorithm for approximating the volume of convex bodies.  J. ACM, 38(1):1–17, 1991.
  * [10] G. Elekes.  A geometric inequality and the complexity of computing volume.  Discr. Comput. Geom., 1(4):289–292, 1986.
  * [11] I.Z. Emiris and V. Fisikopoulos.  Practical polytope volume approximation.  ACM Trans. Math. Soft., 44(4):38:1–38:21, 2018.  Prelim. version: Proc. Symp. Comp. Geometry, 2014.
  * [12] Komei Fukuda.  Frequently asked questions in polyhedral computation.  ftp://ftp.math.ethz.ch/users/fukudak/reports/polyfaq040618.pdf, 2004\.
  * [13] E. Gover and N. Krikorian.  Determinants and the volumes of parallelotopes and zonotopes.  Linear Algebra and its Applications, 413:28–40, 2010.
  * [14] A.K. Kopetzki, B. Schürmann, and M. Althoff.  Methods for order reduction of zonotopes.  pages 5626–5633, 2017.
  * [15] W. Kühn.  Rigorously computed orbits of dynamical systems without the wrapping effect.  Computing, 61:47–67, 1998.
  * [16] Y.T. Lee and S. Vempala.  Convergence rate of Riemannian Hamiltonian Monte Carlo and faster polytope volume computation.  In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2018, pages 1115–1121, 2018.
  * [17] L. Lovász, R. Kannan, and M. Simonovits.  Random walks and an ${O}^{*}(n^{5})$ volume algorithm for convex bodies.  Random Structures and Algorithms, 11:1–50, 1997.
  * [18] L. Lovász and S. Vempala.  Simulated annealing in convex bodies and an O${}^{*}(n^{4})$ volume algorithms.  J. Computer and System Sciences, 72:392–417, 2006.
  * [19] K.G. Murty.  Ball centers of special polytopes.  Department of Industrial ad Operations Engineering, University of Michigan, 2009.
  * [20] A. Pereira and M. Althoff.  Safety control of robots under computed torque control using reachable sets.  In 2015 IEEE International Conference on Robotics and Automation (ICRA), pages 331–338, May 2015.