A Practical Algorithm for Volume Estimation based on Billiard Trajectories and Simulated Annealing
Article No.: 1.3, Pages 1 - 34
Abstract
We tackle the problem of efficiently approximating the volume of convex polytopes, when these are given in three different representations: H-polytopes, which have been studied extensively, V-polytopes, and zonotopes (Z-polytopes). We design a novel practical Multiphase Monte Carlo algorithm that leverages random walks based on billiard trajectories, as well as a new empirical convergence test and a simulated annealing schedule of adaptive convex bodies. After tuning several parameters of our proposed method, we present a detailed experimental evaluation of our tuned algorithm using a rich dataset containing Birkhoff polytopes and polytopes from structural biology. Our open-source implementation tackles problems that have been intractable so far, offering the first software to scale up in thousands of dimensions for H-polytopes and in the hundreds for V- and Z-polytopes on moderate hardware. Last, we illustrate our software in evaluating Z-polytope approximations.
1 Introduction
Volume computation is a fundamental problem with many applications [10, 33, 58, 62, 63]. It is computationally hard: \(\#\)P-hard for explicit polytopes [21, 28], and APX-hard [23] for convex bodies in the oracle model. Therefore, a great effort has been devoted to randomized approximation algorithms, based on Monte Carlo methods and Markov chains for sampling. The first celebrated result given in [22] with complexity \({\mathcal {O}^*} (d^{23})\), where \({\mathcal {O}^*} (\cdot)\) suppresses polylog factors and dependence on error parameters, and d is the dimension. Improved algorithms reduced the exponent to 5 [46]. The latter led to the first practical implementation [24] for high dimensional polytopes. Further results [16, 48] reduced the exponent to 3, followed by another practical implementation [17]. There are two limitations in current algorithm implementations. First, they only support polytopes given by a set of inequalities (H-polytopes) and, second, they scale typically only up to a few hundreds of dimensions. Our aim is to provide the first practical method and its implementation that scales up to thousands dimensions for H-polytopes and to a few hundreds for other polytope representations.
Interestingly, Monte Carlo methods based on Markov chains have gained a lot of interest in the literature for several important problems. For example, in maximum likelihood estimation [27], machine learning [8], finance [10], optimal control [31, 32], Bayesian inference [25], optimization [5, 19, 37, 50, 59], and sampling [20, 47, 51, 53, 54].
Considering volume approximation, in the most general setting, a convex body is given by a (membership or separation) oracle. Thus, the complexity of a randomized volume algorithm is given by an upper bound on the number of oracle calls. Generally, a randomized algorithm uses a Multiphase Monte Carlo (MMC) technique that reduces volume approximation of a convex body P to computing a product of ratios of integrals. The MMC we employ defines a sequence of convex bodies \(P_k\subseteq \dots \subseteq P_0=P\). Then, the volume of P is given by the telescoping product of volume ratios:Clearly, the number of volume ratios is equal to the number of bodies in MMC.
\begin{equation} \mbox{vol}(P) = \mbox{vol}(P_k) \frac{\mbox{vol}(P_{k-1})}{\mbox{vol}(P_k)} \cdots \frac{\mbox{vol}(P_0)}{\mbox{vol}(P_1)}. \end{equation}
(1)
Geometric random walks1 come into the picture to estimate each ratio by sampling from an appropriate distribution. A random walk in polytope P starts at some interior point and, at each step, moves to a “neighbouring” point, that we choose according to some distribution depending only on the current point. The complexity of a random walk is determined by two factors: its mixing time, which equals the number of steps (or walk length) required to make the distance between the current and the target distribution close to null, and the complexity of the basic operations that we perform at each step of the walk, called per step cost.
A crucial observation is that, if the ratios are bounded by some constant, then acceptance-rejection sampling would suffice to accurately estimate them, i.e., from Chebyshev inequality we obtain that for some volume ratio \(i\lt k\), about \(k\cdot (\mbox{vol}(P_{i-1})/\mbox{vol}(P_i))\) uniformly distributed points in \(P_i\) suffice to estimate that volume ratio. All in all, the total number of operations performed by an MMC volume algorithm is
\begin{equation} k \times \#\text{points per volume ratio}\, \times \#\text{operations per point}. \end{equation}
(2)
The relations \(rB_d\subseteq P\subseteq RB_d\), for \(r\lt R\) and \(B_d\) the d-dimensional unit ball, imply that \(k=\lceil d\lg (R/r)\rceil\). We call \(R/r\) the sandwiching ratio of P. For general convex bodies, the sandwiching ratio is minimized when the body is well-rounded, i.e., \(B_d\subseteq P\subseteq \mathcal {O} (\sqrt {d})B_d\). In [46] they put P in isotropic position, which implies well-roundness, and use scaled copies of the d-dimensional unit ball \(B_d\) to define the MMC, i.e., \(P_i=(2^{(k-i)/d}B_d)\cap P\), for \(i=0, \ldots , k=\Theta (d\lg d)\). These techniques combined with the so called speedy walk–a special case of the ball wall–delivers an \({\mathcal {O}^*} (d^5)\) algorithm.
In [48] they fix a sequence of functions \(f_0,\dots , f_k\) and \(\mbox{vol}(P)\) is given by a telescopic product of ratios of integrals generalizing Equation (1). They use a sequence of exponential functions; after putting P in well-rounded position in \({\mathcal {O}^*} (d^4)\) oracle calls, the number of integral ratios in MMC \(k={\mathcal {O}^*} (\sqrt {d})\) and the total number of oracle calls is \({\mathcal {O}^*} (d^4)\) by using Hit-and-Run (HR). In [16] they employ a sequence of spherical Gaussian distributions and an improved mixing time for the ball walk. They assume that the input convex body is well-rounded and the total number of oracle calls is \({\mathcal {O}^*} (d^3)\). In [34] they provide an improved rounding algorithm which puts a convex body in well-rounded position in \(O(d^{3.5})\) oracle calls (that becomes \(O(d^3)\) after recent progress in KLS conjecture [15]). Thus, the latter complexity also determines the bound for volume approximation of general convex bodies. For H-polytopes, with m facets, in [45] they extend the simulated annealing of [48] to the manifold setting and together with the Riemannian Hamiltonian Monte Carlo sampling, they obtain a \(\mathcal {O} ^*(m^2d^{\omega -1/3})\) bound, where \(\omega\) stands for the matrix multiplication exponent. In [49], they show a \(\mathcal {O} ^*(md^{4.5}+md^{4})\) bound via a sub-linear ball walk and Gaussian cooling.
However, it is impractical2 to use an implementation of the above volume approximation algorithms due to extremely large constants in the complexity and pessimistic upper bounds on the mixing time of the random walks they employ. Considering practical volume algorithms, the objective is to perform efficient computations to obtain the volume of a high dimensional polytope. Typically, to achieve this goal, a practical algorithm on one hand, is based on a theoretical algorithmic scheme and on the other hand, provides practical adjustments, heuristics and efficient tuning. Moreover, practical algorithms instead of theoretical guarantees, provide strong empirical evidences for the run-time and the accuracy of the method through extended experimental results. For example, the method in [24] is based on the ball sequence in [22] and performs \(\sim md^4\) operations to achieve a relative error at most 0.1 in practice. The method in [17] is based on the Gaussian cooling in [16], and performs \(\sim md^3\) operations to achieve the same relative error.
Considering sampling in practice, the dominant paradigm for random walks is CDHR [38]. It is a version of HR that uses directions parallel to the axes. HR mixes after \({\mathcal {O}^*} (d^3)\) [47] steps for isotropic log-concave distributions. Very recently, the mixing time of CDHR has been bounded by \({\mathcal {O}^*} (d^9(R/r)^2)\) [44, 53]. However, experiments [24, 30] indicate that both CDHR and HR have similar convergence rates in practice. Thus, the faster step of CDHR compared to HR—\(O(m)\) vs. \(O(md)\)—is the main reason why CDHR overshadowed, until recently, all other random walks in practical computations on H-polytopes.
A certainly interesting and overlooked problem is to compute in practice the volume of a polytope given by other representations, namely, V-polytopes (given by the set of extreme points or vertices) and Z-polytopes or zonotopes, given as a Minkowski sum of segments. For those representations all known practical algorithms [17, 24, 26] fail to compute the volume for dimensions beyond, say, \(d\gt 15\). There are two main reasons for that, besides the costly oracles (which are equivalent to a linear program). First, it is not possible to compute efficiently the largest inscribed ball needed as an initialization by all algorithms. In particular, for a Z-polytope P and a ball B, checking whether \(B\subseteq P\) is in co-NP while for P being a V-polytope, given \(p\in P\), the computation of the largest inscribed ball centered at p is NP-hard [52]. Second, the number of facets–which is requested by the implementation of [17]–is typically exponential in d for both Z- and V-polytopes.
A powerful approach to obtain well roundness is to put P in near isotropic position [16, 48]. In [1], they prove that \(\mathcal {O} ^*(d)\) uniformly distributed points in P suffice to achieve 2-isotropic position after some proper linear transformations. However, this could be quite computationally expensive in practice due to extensive uniform sampling from skinny polytopes. In [30], they provide an alternative method which achieves efficiency comparable in practice to that in [17]. They compute the maximum volume ellipsoid in P, they map it to the unit ball and apply to P the same transformation. They experimentally show that a few iterations suffice to put P in John’s position [35]. However, when P is a V-polytope there is no known algorithm to compute the largest inscribed ellipsoid of P. For a Z-polytope, the John ellipsoid E is tighter, as for any centrally symmetric convex body compared to non-symmetric convex bodies. However, there is no algorithm for the computation of the John ellipsoid of a Z-polytope. The algorithm in [12] computes an ellipsoid \(E: \ E\subseteq P\subseteq dE\) for a given Z-polytope; the complexity is polynomial but the degree is not determined. Thus, it is unclear whether it may lead to an efficient rounding and whether it can be practical for Z-polytopes.
Interestingly, for well-rounded convex polytopes, using any sequence of bodies to fix the telescopic product leads to \(k = \Theta ^*(d)\) phases in the telescopic product [48], while using a sequence of any log-concave functions leads to \(k = \Theta ^*(\sqrt {d})\) phases [16]. Therefore, the fact that our practical algorithm is faster than that in [17] apparently contradicts theory [16, 46]. The main reason behind this “contradiction” is the efficiency of the improved version of Billiard Walk we use. In particular, extended experiments have shown that—for rounded convex polytopes—Billiard Walk mixes after \(O^*(1)\) steps in practice [14], while its per step cost is \(O(md)\). On the other hand, both HR and CDHR mix after \(O^*(d^2)\) steps in practice [17, 30].
1.1 Our Contributions
We propose the first practical volume algorithm that scales to thousands of dimensions for H-polytopes outperforming all the existing practical algorithms. Additionally, it efficiently handles two more representations of polytopes, namely V- and Z-polytopes, up to a few hundred dimensions for the first time. We provide an extended experimental analysis that provides empirical evidences on the accuracy as well as the performance of our practical algorithm. Moreover, our experiments show that our proposed practical algorithm outperforms the state-of-the-art, namely [17] (e.g. for unit cubes, 20 times faster in \(d=50\) and 95 times faster in \(d=500\); for \(d=1,000\) our implementation is expected to be \(\sim 130\) times faster; see Section 5.2 and Figure 6 for more details).
The algorithm (Algorithm 3) is based on the MMC framework and combines the sequence of convex bodies with a new simulated annealing schedule, while leveraging uniform sampling and the superior practical performance of Billiard walk [56]. In brief, the algorithm fixes a convex body \(C\in {\mathbb {R}}^d\) and a sequence \(C_1\supseteq \cdots \supseteq C_k\) that intersects P, where \(C_i = q_iC\), \(q_i\in {\mathbb {R}}_+\), using a cooling or annealing schedule (Algorithm 5). The algorithm uses a telescopic product (Equation (7)) of ratios of volumes, each bounded by a constant with high probability by the annealing schedule. With high probability, \(k=\mathcal {O} (\lg (\mbox{vol}(P)/\mbox{vol}(P_k)))\), where P is the input polytope and \(P_k = C_k\cap P\) the body with minimum volume in MMC (Section 4.2.1), which is not surprising, as—if C is a ball—this results to \({\mathcal {O}^*} (d)\) bodies in MMC with \({\mathcal {O}^*} (d)\) points required per volume ratio. However, our scheme offers crucial advantages in practice compared to the state-of-the-art by (a) optimizing the hidden constants in the complexity (Propositions 3 and 4), (b) taking advantage of fast practical uniform samplers, and (c) exploiting good choices of C in the Z-polytope case.
To fix the sequence of \(C_i\) we introduce a new simulated annealing schedule. The schedule minimizes the number of phases k given a body C with high probability—assuming perfect uniform sampling. To achieve this we define a new statistical test that restricts each ratio in the telescopic product in Equation (1) in the interval \([r, r+\delta ]\), where r and \(\delta\) are predefined constants.
For sampling, we leverage the Billiard walk for the first time in volume computation (Section 2). We show that, with the right selection of parameters, this walk outperforms CDHR even if the walk length is 1. Using this walk length, our algorithm only needs \({\mathcal {O}^*} (1)\) points per volume ratio (Section 5 and Figures 8 and 9), to ensure the desired accuracy (Figure 2). We provide efficient implementations for the step of billiard_walk for all polytope representations. For Z- and V-polytopes, the steps are reduced to solving linear programs while for H-polytopes, with m facets, we exploit the fast implementation of billiard_walk in [14] to extend it to the case of the intersection with a ball. Our efficient implementation requires \(O(dm)\) operations per billiard_walk point.
For rounding, we propose methods for all types of polytopes (Section 3) that are faster than the ones currently available [17, 24, 30] as reported by our experiments (Section 5). For H-polytopes, we simply enhance the method from [17] with the billiard sampler. For V-polytopes, our method brings P to approximate John’s position by computing the smallest enclosing ellipsoid using P’s vertices and the algorithm in [61]. For Z-polytopes, we construct an H-polytope that approximates P well, and set it as C in MMC. We experimentally show that this choice, with no rounding, is more efficient than estimating the volume after applying various rounding methods.
We offer an optimized implementation of the proposed algorithm that scales in thousands of dimensions for H-polytopes and hundreds of dimensions for V- and Z-polytopes on moderate hardware (Section 5) for the first time. Our software performs computations which were intractable until now (Tables A.3 and A.4), outperforms current state-of-the-art implementations [17, 24] in all representations, and is illustrated in evaluating Z-polytope approximations in the context of mechanical engineering [41]. Our implementation3 builds upon and enhances volesti, a C++ open source library for high-dimensional sampling and volume computation with an R interface [13].
Notation. We denote the dimension that the polytope P lies with d. For H-polytopes m denotes the number of facets. For V- and Z-polytopes we denote the number of vertices and the number of generators respectively with n. Considering algorithms, with volume we refer to our practical algorithm and its implementation, and with CoolingGaussian we refer to the implementation of the practical algorithm in [17].
2 Uniform Sampling and Polytope Oracles
In this section we introduce the proposed uniform sampler that the volume algorithm will use and discuss polytope operations.
We start with some notation. P is a full-dimensional convex polytope lying in d-dimensional space and \(\partial P\) is its boundary. An H-polytope (in H-representation) iswith m facets in d dimensions. A V-polytope is the convex hull of a pointset in \({\mathbb {R}}^d\); equivalently it can be seen as a linear projection, given by a matrix \(V\in {\mathbb {R}}^{d\times n}\), of the canonical simplex, i.e.,while the columns of V corresponds to the n vertices of P. A zonotope (Z-polytope) is the Minkowski sum of n d-dimensional segments or equivalently given by matrix \(G\in {\mathbb {R}}^{d\times n}\) and seen as a linear map of hypercube \([-1,1]^n\) to \({\mathbb {R}}^d\): we call it a Z-representation. The order of a Z-polytope is the ratio \(n/d\). We denote by \([n]\) the set \(\lbrace 1,2,\dots ,n\rbrace\) for some natural number n.
\begin{equation} P=\lbrace x\ |\ Ax\le b,\ A\in {\mathbb {R}}^{m\times d},\ b\in {\mathbb {R}}^{m} \rbrace \end{equation}
(3)
\begin{equation*} \left\lbrace x\in {\mathbb {R}}^n\ |\ \sum _{i=1}^nx_i=1,\ x_i\ge 0 \right\rbrace , \end{equation*}
Random sampling is fundamental for efficient volume approximation. Our algorithm employs Billiard walk (Algorithm 1, billiard_walk) to sample approximately uniform points from a polytope P. billiard_walk starts from a given point \(p_0 \in P\), selects uniformly at random a direction, say \(v_0\), and it moves along the direction of \(v_0\) for length L; it reflects on the boundary if necessary. This results in a new point \(p_1\) inside P. We repeat the procedure from \(x_1\). Asymptotically it converges to the uniform distribution over P. The length is \(L=-\tau \ln \eta\), where \(\eta\) is a uniform number in \((0,1)\), that is \(\eta \sim \mathcal {U}(0,1)\), and \(\tau\) is a predefined constant. It is useful to set a bound, say \(\rho\), on the number of reflections to avoid computationally hard cases where the trajectory may be stuck in corners. In [56] they set \(\tau \approx \mathop {diam}(P)\) and \(\rho =10d\).
billiard_walk requires an access to a boundary oracle, that is, to compute the intersection of a ray \(\ell :=\lbrace p+tv,\ t\in {\mathbb {R}}_+ \rbrace\) with \(\partial P\), where \(p\in P\). Additionally, the walk needs the tangent plane on the intersection point to compute the reflection of the billiard trajectory.
When the input is an H-polytope, see Equation (3), the inner vector of the boundary point’s facet is given by a row of A. To compute the boundary point \(\partial P\cap \ell\) for the first reflection, we solve \(a_j^T(p_0 + tv_0) = b_j,\ j\in [m]\), for t, and keep the smallest positive root, where \(a_j\) is the j-th row of A. This takes dm arithmetic operations in total. A straightforward approach for billiard_walk would consider that each reflection costs \(\mathcal {O} (md)\) and thus each step costs \(\mathcal {O} (\rho md)\). However, we exploit the fast version of billiard_walk in [14]. It improves both point and direction updates (see related comments in Algorithm 1) by storing computations from the previous iteration combined with preprocessing. The preprocessing computes the inner products of all possible pairs of normal vectors of the facets, that takes \(m^2 d\) operations. Then, the first reflection of each billiard_walk step costs \(\mathcal {O} (md)\) operations and the rest ones cost \(\mathcal {O} (m)\) operations each. Thus, the amortized per-step complexity to sample from a H-polytope with billiard_walk becomes \(\mathcal {O} ((\rho + d)m)\).
The amortized per-step complexity of the billiard_walk (Algorithm 1) for a polytope \(P\in {\mathbb {R}}^d\) with m facets, given in H-representation, is \(\mathcal {O} ((\rho + d)m)\) operations after preprocessing that takes \(\mathcal {O} (m^2d)\) operations, where \(\rho\) is the maximum number of reflections per step.
In the context of the volume algorithm (Section 4) when the body used in MMC is a ball we have to sample from the intersection of an H-polytope with a ball using billiard_walk. Assuming that the center of the ball is the origin with radius R, to compute the intersection of the ray \(\ell :=\lbrace p+tv,\ t\in {\mathbb {R}}_+ \rbrace\) with the boundary of the ball we compute the positive root of the equation,which takes \(\mathcal {O} (d)\) operations. The normal to the tangent plane of a ball at point p is \(p/ \Vert p\Vert _2\). Thus, to compute the reflection of the ray when it hits the boundary of the ball takes \(\mathcal {O} (d)\) operations. Since \(m\gt d\) the amortized per-step complexity of billiard_walk is \(\mathcal {O} ((\rho + d)m)\).
\begin{equation} t^2 + 2(p^Tv)t + \Vert p\Vert _2^2 - R^2 = 0, \end{equation}
(4)
For both Z- and V-polytopes, given a ray \(\ell\), computing \(\partial P\cap \ell\) reduces to a linear program (LP). One may consider the pre-image of \(\ell\), hence computation of \(\partial P\cap \ell\) becomes an optimization problem in the region defined by the intersection of a linear subspace with the hypercube or the simplex, respectively. Let \(g_i\) be the generators of a Z-polytope, \(i\in [n]\). The intersection of \(\partial P\cap \ell\) is given by the following LP:If \(v_i\) are the vertices of a V-polytope, \(i\in [n]\), then the intersection of \(\partial P\cap \ell\) is given as follows:
\begin{equation} \begin{split}\max \ t,\hspace{8.5359pt}s.t.\ \hspace{5.69046pt} & p+tv = \sum _{i=1}^{n} x_ig_i,\\ & -1\le x_i\le 1,\ i\in [n]. \end{split} \end{equation}
(5)
\begin{equation} \begin{split}\max \ t,\hspace{8.5359pt}s.t.\ \hspace{5.69046pt} & p+tv = \sum _{i=1}^{n} x_iv_i, \\ & x_i\ge 0,\ \sum _{i=1}^nx_i=1,\ i\in [n]. \end{split} \end{equation}
(6)
To compute the reflection of the ray we exploit the solutions of the above LPs. For V-polytopes we keep the vertices that correspond to the positive \(x_i\), which generate the facet that \(\ell\) hits. Then, we compute the normal to the hyperplane they define. Similarly, for Z-polytopes we keep the generators that correspond to \(x_i\ne \lbrace -1,1\rbrace\). Computing the normal in both cases requires to solve a \(d\times d\) linear system.
Volume estimation also requires a membership oracle which decides whether point p lies in P (see estimation of \(r_{k+1}\) in Equation (7)). For an H-polytope one has to check at most m inequalities, thus it costs \(\mathcal {O} (dm)\) operations. For both Z- and V-polytopes, membership oracles are feasibility problems reducing to LP. For a V-polytope, \(p\in P\) if and only if the following region is feasible:For a Z-polytope, \(v_i\) corresponds to the generators of P and the rest of the constraints become \(-1\le x_i\le 1\).
\begin{equation*} \sum _{i=1}^n x_i v_i = p,\ x_i\ge 0,\ \sum _{i=1}^n x_i =1. \end{equation*}
3 Preprocessing: to Round or Not to Round
We introduce two different preprocessing methods. First, a rounding method that transforms a polytope of any representation to a so-called near isotropic position. Second, we propose a preprocessing method for Z-polytopes without rounding.
We employ the rounding in [17] to introduce a faster method that brings P to near-isotropic position, using billiard_walk for uniform sampling instead of Hit-an-Run. The method samples uniformly from P, then computes a linear transformation that puts the sample to isotropic position and applies to P the same transformation. To map the pointset to isotropic position, we compute the SVD decomposition of the matrix that contains the generated points in rows [3]. We iterate this procedure until the ratio of the maximum over the minimum singular value falls below a threshold, e.g. 4.
Next, we propose a different preprocessing with no rounding that can be more efficient than rounding for some specific polytopes such as Z-polytopes (Section 5). The key idea is the choice of C at MMC so that it is both a good approximation of P and we can efficiently sample from it. The first property affects the number of bodies in MMC while the second the runtime per volume ratio in the telescopic product in Equation (7).
When P is a Z-polytope with generator matrix \(G\in {\mathbb {R}}^{d\times n}\), we compute a centrally symmetric H-polytope \(C\subset P\) (Algorithm 2). Since P is the projection of the hypercube \([-1,1]^n\) onto \({\mathbb {R}}^d\) by G, our method projects a linear subspace of \([-1,1]^n\) onto \({\mathbb {R}}^d\) such that (a) it is easy to compute the H-representation of the underlying polytope C and (b) we minimize the information loss for P implied by the restriction of \([-1,1]^n\) to that subspace. We exploit the basic idea of Principal Component Analysis (PCA) [36] and obtain C by projecting with G the linear subspace spanned by the kernel of \(G^TG\) restricted in \([-1,1]^n\). The matrix \(G^TG\) has \(n-d\) zero eigenvalues. We compute an orthonormal basis of its kernel matrix \(E\in {\mathbb {R}}^{n\times (n-d)}\), whose SVD is given below, and use the orthogonal complement \(W_{\perp }\) defined as follows:
\begin{equation*} E = 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. \end{equation*}
In short, we use the linear map induced by \(W_{\perp }\) to compute an H-representation of C. By construction, \(C\subset P\) since the domain of C is a subset of the domain of P (subspace of the hypercube) under the same mapping G.
For an illustration of an example of C and P see Figure 1 (right).
Fig. 1.
4 Volume Algorithm
In this section, we detail the volume algorithm (Algorithm 3) based on uniform sampling and the preprocessing of Sections 2 and 3. The algorithm uses the following telescopic product:where \(P_0=P,\ P_i=C_i\cap P,\, i=1, \ldots , k\). The annealing schedule bounds each ratio \(r_i\) by a constant, with high probability. The algorithm uses a ratio estimation (Algorithm 6) to estimateusing billiard_walk (Algorithm 1) to sample uniformly from \(P_{i}\), and acceptance/rejection in \(P_{i+1}\). The algorithm also estimates the ratiousing uniform sampling from \(C_k\) and acceptance/rejection in \(P_{k}\). Lastly, estimate_ratio uses a sliding window of the last l ratios and an empirical criterion to declare convergence in a statistical sense.
\begin{equation} \mbox{vol}(P) = \frac{\frac{\mbox{vol}(P_k)}{\mbox{vol}(C_k)}}{\frac{\mbox{vol}(P_1)}{\mbox{vol}(P_0)}\frac{\mbox{vol}(P_2)}{\mbox{vol}(P_1)}\cdots \frac{\mbox{vol}(P_{k})}{\mbox{vol}(P_{k-1})}} \, \mbox{vol}(C_k), \end{equation}
(7)
\begin{equation*} r_i = \frac{\mbox{vol}(P_i)}{\mbox{vol}(P_{i-1})},\; 1\le i\le k, \end{equation*}
\begin{equation*} r_{k+1} = \frac{\mbox{vol}(P_k)}{\mbox{vol}(C_k)} \end{equation*}
Algorithm volume is parameterized by: the error of approximation \(\epsilon\), cooling parameters \(r,\delta \gt 0\), \(0\lt r+\delta \lt 1\), significance level (s.l.) \(\alpha \gt 0\) of the statistical tests, the degrees of freedom \(\nu\) for the t-student used in t-tests, and N that controls the number of points \(\nu N\) generated per \(P_i\) all used in the annealing schedule.
4.1 Annealing Schedule for Convex Bodies
Given P and \(q_{\min }C \subseteq P \subseteq q_{\max }C\), annealing_schedule generates the sequence of convex bodies \(C_1\supseteq \cdots \supseteq C_k\) defining \(P_i=C_i\cap P\) and \(P_0=P\). The computation of \(q_{\min },q_{\max }\) is not trivial but in Section 5 we give practical choices depending on C which are very efficient in practice. The main goal is to restrict each ratio \(r_i\) to the interval \([r, r+\delta ]\) with high probability.
We introduce some notions from statistics needed to define two tests and refer to [18] 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 st.d. and \(t_{\nu -1}\) is the t-student distribution with \(\nu -1\) degrees of freedom. Given a significance level \(\alpha \gt 0\) we test the null hypothesis for the mean value of the population, \(H_0:\mu \le \mu _0\) against \(H_1:\mu \gt \mu _0\). We reject \(H_0\) if,which implies \(\Pr [\text{reject }H_0\ |\ H_0\ \text{true} ] =\alpha\). Otherwise we fail to reject \(H_0\). In the sequel, we define two statistical tests for a volume ratio, which can be reduced to t-tests:
\begin{align*} t\ge t_{\nu -1,\alpha } \iff \bar{x}\ge \mu _0+t_{\nu -1,\alpha }{s}/{\sqrt {\nu }}, \end{align*}
\begin{align} \text{U-test}(P_1,P_2)\quad & H_0: \mbox{vol}(P_2)/\mbox{vol}(P_1)\ge r+\delta \end{align}
(8)
\begin{align} \text{L-test}(P_1,P_2)\quad & H_0: \mbox{vol}(P_2)/\mbox{vol}(P_1)\le r \end{align}
(9)
The U-test and L-test are successful if and only if null hypothesis \(H_0\) is rejected, namely \(r_i\) is upper bounded by \(r+\delta\) or lower bounded by r, with high probability, respectively.
If we sample N uniform points from a body \(P_{i-1}\) then r.v. X that counts points in \(P_i\), 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)\) follows the Normal distribution. The Algorithm 4 provides all the steps to perform both statistical tests.
Remark. The normal approximation of r.v. \(Y=X/N\) 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)\gt 10\).
Then each sample proportion that counts successes in \(P_i\) over N is an unbiased estimator for \(\mathbb {E}[Y]\), which is \(r_i\). So if we sample \(\nu N\) points from \(P_{i-1}\) 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 U-test and L-test. Let \(\hat{\mu }\) be the sample mean of the ratios. Ifthen both U-test and L-test are successful and \(r_i\) is restricted to \([r,r+\delta ]\) with high probability.
\begin{equation} r+\delta -t_{\nu -1,\alpha }\frac{s}{\sqrt {\nu }}\ge \hat{\mu }\ge r +t_{\nu -1, \alpha }\frac{s}{\sqrt {\nu }} \end{equation}
(10)
Let us now describe annealing_schedule in more detail: Given \(q_{\min }C \subseteq P \subseteq q_{\max }C\), annealing_schedule computes the sequence \(C_1\supseteq \dots \supseteq C_k\). Each body \(C_i\) 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. The initialization step of annealing_schedule computes the body with minimum volume, denoted by \(C^{\prime }\) which is set as \(C_k\) at termination of annealing_schedule. To be precise, the last body in the sequence shall be \(C_k\), s.t. \(r_{k+1} = \mbox{vol}(P_k)/\mbox{vol}(C_k) \in [r,r+\delta ]\) with high probability. In particular, at initialization, annealing_schedule binary searches in the interval \([q_{\min },q_{\max }]\) to compute a scalar \(q^{\prime }\) s.t. \(C^{\prime }=q^{\prime }C\) and both U-test(\(C^{\prime },C^{\prime }\cap P\)) and L-test(\(C^{\prime },C^{\prime }\cap P\)) are successful.
To compute \(C_i = q_iC\), annealing_schedule computes the scalar \(q_i\) using binary search in a proper interval, such that both U-test(\(P_{i-1},q_iC\cap P\)) and L-test(\(P_{i-1},q_iC\cap P\)) are successful, where \(P_0=P,\ P_i = C_i\cap P\) and \(i=1, \ldots , k\). As \(C^{\prime }=q^{\prime }C\) is the body of smallest volume in the sequence of MMC the interval to binary search for \(q_i\) has to be \([q^{\prime },q_{i-1}],\ i=1,\dots , k\), where \(q_0 = q_{\max }\). At termination the annealing_schedule sets \(q_k = q^{\prime }\). Notice that the binary search in the interval \([q^{\prime },q_{i-1}]\) implies that \(\mbox{vol}(C^{\prime }\cap P) \lt \mbox{vol}(P_i) \lt \mbox{vol}(P_{i-1})\) and if both L-test and U-test are successful then \(r_i\in [r,r+\delta ]\) with high probability. Also recall that in each step of binary search which computes \(q_i\), annealing_schedule samples \(\nu N\) points from \(q_{i-1}C\cap P\) to check both U-test and L-test.
annealing_schedule employs \(C^{\prime }\) to decide stopping at the i-th step after computing \(C_i\); if the criterion fails, the algorithm computes \(C_{i+1}\). In particular, it checks whether \(\mbox{vol}(P_i)/ \mbox{vol}(C^{\prime }\cap P) \ge r\) with high probability, using L-test. Formally,
Stop in step i if L-test(\(P_i,C^{\prime }\cap P\)) holds. Then, set \(k=i+1\), \(q_k=q^{\prime }\) and \(C_k = q_kC,\ P_k=C_k\cap P\).
To perform L-test, \(\nu N\) uniformly distributed points are sampled from \(P_i\).
4.2 Termination
This section proves that annealing_schedule terminates with constant probability.
Let J be the minimum number of steps by annealing_schedule, corresponding to no errors occurring in the t-tests. Let also annealing_schedule actually perform \(M\ge J\) steps and \(\beta _{\max }\), \(\beta _{\min }\) be the maximum and minimum among all \(\beta\)’s in the M pairs of t-tests in the U-test and L-test, respectively. Then, annealing_schedule terminates with constant probability, namely:
\begin{align*} \Pr [{{\rm\small ANNEALING\_SCHEDULE}\ terminates}] \ge \ & 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}. \end{align*}
We demonstrate halting of annealing_schedule. In the t-tests, errors of different types may occur, thus, binary search may enter intervals that do not contain ratios in \([r,r+\delta ]\), which implies that there is a (smaller) probability that annealing_schedule fails to terminate. We bound the probability that annealing_schedule enters inappropriate intervals in a certain step–and thus we bound the probability of failing to terminate–by a constant. Let \(\beta\) capture the power of a t-test: \(pow = \Pr [\text{reject}\ H_0\ |\ H_0\ \text{false}] = 1-\beta\).
Let \(M\ge J\). If annealing_schedule fails to terminate after M pairs of U-test and L-test, 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 quantile function of t-student. For the latter probability we write \(\beta _L,\ \beta _R\) for U-test and L-test, respectively. If, for a pair of tests, both null hypotheses are false then an error occurs with probabilitySimilarly,Then,□
\begin{align*} p_1 &= \Pr [\text{error occurs}\ |\ \text{both }H_0\text{ false}] =\beta _L(1-\beta _R) + \beta _R(1-\beta _L) + \beta _L\beta _R \\ &= \beta _L + \beta _R - \beta _L\beta _R \le 2\beta _{\max } - \beta ^2_{\min } \end{align*}
\begin{align*} p_2 &= \Pr [\text{error occurs}\ |\ \text{U-test }H_0 \text{ false } \text{and L-test }H_0\text{ true}]\\ &= \alpha (1-\beta _L) + \alpha \beta _L + \beta _L(1-\alpha) \\ & = \alpha + \beta _L - \alpha \beta _L \le \alpha + \beta _{\max } - \alpha \beta _{\min } \end{align*}
\begin{align*} p_3 &= \Pr [\text{error occurs}\ |\ \text{L-test }H_0 \text{ false } \text{and U-test }H_0\text{ true}]\\ &= \alpha (1-\beta _R) + \alpha \beta _R + \beta _R(1-\alpha) \\ & = \alpha + \beta _R - \alpha \beta _R \le \alpha + \beta _{\max } - \alpha \beta _{\min } \end{align*}
\begin{align*} \Pr &[{{\rm{\small {ANNEALING\_SCHEDULE}}}\ \text{fails to terminate}}] \le \sum _{i=1}^M p_1^i + \sum _{i=1}^M p_2^i + \sum _{i=1}^M p_3^i\\ &\le \sum _{i=1}^{\infty } p_1^i + \sum _{i=1}^{\infty } p_2^i + \sum _{i=1}^{\infty } p_3^i = \frac{1}{1-p_1} + \frac{1}{1-p_2} + \frac{1}{1-p_3} - 3 \\ &\le 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}. \end{align*}
4.2.1 Number of Bodies in MMC.
We give probabilistic bounds on the number of bodies in MMC. Let us assume that (i) we sample perfect uniform points in each step of annealing_schedule and (ii) that annealing_schedule terminates successfully. We offer a probabilistic upper bound on the number of bodies in MMC k and then a probabilistic interval where k lies.
Given a convex polytope \(P\subset {\mathbb {R}}^d\) with sandwiching ratio \(R/r\), i.e., \(rB_d\subseteq P \subseteq RB_d\), and cooling parameters \(r,\delta\) such that \(r+\delta \lt 1/2\) and parameters \(\alpha\), N, \(\nu\), \(q_{\min }\), \(q_{\max }\) let k be the number of convex bodies in MMC returned by Algorithm annealing_schedule, when C is the unit ball. Then \(\Pr [k\le \lceil d\lg (R/r)\rceil ] \ge 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 L-test.
If \(k \gt \lceil d\lg (R/r)\rceil\) holds then \(k\ge 1\) type I errors of U-test occurred in annealing_schedule i.e., \(H_0\) holds but the test rejects it, while L-test 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 L-test is \(1-\beta _R\). Let \(\beta _{\min }\) be the minimum among the values of the quantile function appearing in all instances of L-test. Then,□
\begin{align*} \Pr [k\gt \lceil d\lg (R/r)\rceil ] &\le \sum _{i=1}^k (a(1-\beta _i))^i \le \sum _{i=1}^{\infty } (a (1-\beta _{\min }))^i\\ &= \frac{1}{1-\alpha (1-\beta _{\min })}-1\Rightarrow \\ \Pr [k\le \lceil d\lg (R/r)\rceil ] &\ge 2-\frac{1}{1-\alpha (1-\beta _{\min })}\\ &= 1 - \frac{\gamma }{1-\gamma }, \mbox{ where } \gamma =\alpha (1-\beta _{\min }). \end{align*}
Let a convex polytope \(P\subset {\mathbb {R}}^d\), cooling parameters \(r,\delta\) such that \(r+\delta \lt 1/2\), parameters \(\alpha\), N, \(\nu\), \(q_{\min }\), \(q_{\max }\), and k be the number of bodies in MMC by annealing_schedule Ṫhen,where \(\gamma _L = \alpha (1-\beta ^L_{\min })\), \(\gamma _R = \alpha (1-\beta ^R_{\min })\) and \(\beta ^R_{\min },\ \beta ^L_{\min }\) are minimum among all values of the quantile function appearing in U-test and L-test, respectively.
\begin{equation*} \Pr \bigg [\bigg \lfloor \log _{\frac{1}{r+\delta }}\bigg (\frac{{vol}(P)}{{vol}(P_k)}\bigg)\bigg \rfloor \\ \le k\le \bigg \lceil \log _{\frac{1}{r}}\bigg (\frac{{vol}(P)}{{vol}(P_k)}\bigg)\bigg \rceil \bigg ]\ge 1- \frac{\gamma _L}{1-\gamma _L} - \frac{\gamma _R}{1-\gamma _R}, \end{equation*}
Let,
\begin{align*} & \Pr \bigg [\bigg \lfloor \log _{\frac{1}{r+\delta }}\bigg (\frac{\mbox{vol}(P)}{\mbox{vol}(P_k)}\bigg)\bigg \rfloor \le k \le \bigg \lceil \log _{\frac{1}{r}}\bigg (\frac{\mbox{vol}(P)}{\mbox{vol}(P_k)}\bigg)\bigg \rceil \bigg ] \\ & = 1 - \Pr \bigg [k \gt \bigg \lceil \log _{\frac{1}{r}}\bigg (\frac{\mbox{vol}(P)}{\mbox{vol}(P_k)}\bigg)\bigg \rceil \ \text{or}\ k \lt \bigg \lfloor \log _{\frac{1}{r+\delta }}\bigg (\frac{\mbox{vol}(P)}{\mbox{vol}(P_k)}\bigg)\bigg \rfloor \bigg ] \\ & = 1 - \Pr \bigg [k \gt \bigg \lceil \log _{\frac{1}{r}}\bigg (\frac{\mbox{vol}(P)}{\mbox{vol}(P_k)}\bigg)\bigg \rceil \bigg ] - \Pr \bigg [k \lt \bigg \lfloor \log _{\frac{1}{r+\delta }}\bigg (\frac{\mbox{vol}(P)}{\mbox{vol}(P_k)}\bigg)\bigg \rfloor \bigg ]. \end{align*}
Similarly to the proof of Proposition 3, we have \(\Pr [ k \gt \lceil \log _{\tfrac{1}{r}}(\tfrac{\mbox{vol}(P)}{\mbox{vol}(P_k)})\rceil ] \le \tfrac{\gamma _R}{1 - \gamma _R}\), where \(\gamma _R = \alpha (1-\beta _{\min }^R)\), and \(\beta _{\min }^R\) is the minimum among all values of the quantile function appearing in L-test.
If \(k\lt \lfloor \log _{\tfrac{1}{r+\delta }}(\tfrac{\mbox{vol}(P)}{\mbox{vol}(P_k)})\rfloor\), then k type-I errors of L-test occurred, where \(k\ge 1\), while U-test was successful with probability \(1-\beta\). This implies that \(\Pr [ k \lt \lfloor \log _{\tfrac{1}{r+\delta }}(\tfrac{\mbox{vol}(P)}{\mbox{vol}(P_k)})\rfloor ] \le \tfrac{\gamma _L}{1 - \gamma _L}\), where \(\gamma _L = \alpha (1-\beta _{\min }^L)\), and \(\beta _{\min }^L\) is the minimum among all values of the quantile function appearing in U-test.
Putting everything together, it follows that:□
\begin{equation*} \Pr \bigg [\bigg \lfloor \log _{\frac{1}{r+\delta }}\bigg (\frac{\mbox{vol}(P)}{\mbox{vol}(P_k)}\bigg)\bigg \rfloor \le k \le \bigg \lceil \log _{\frac{1}{r}}\bigg (\frac{\mbox{vol}(P)}{\mbox{vol}(P_k)}\bigg)\bigg \rceil \bigg ]\ge 1 - \frac{\gamma _L}{1-\gamma _L} - \frac{\gamma _R}{1-\gamma _R} \end{equation*}
Therefore, the number of bodies is \(k=O(\lg (\mbox{vol}(P)/\mbox{vol}(P_k)))\) with high probability. This implies that k decreases when C is a good fit to P. Clearly, the body that minimizes the number of bodies is the one that maximizes \(\mbox{vol}(C_k\cap P)\).
4.3 Ratio Estimation
volume requires that we estimate \(k+1\) ratios, which have been formed by annealing_schedule. This section describes how to perform these estimations (Algorithm 6). First, we bound the error in each ratio estimation in order to use it for the definition of the stopping criterion. From standard error propagation analysis, for each volume ratio \(r_i = \mbox{vol}(P_i) / \mbox{vol}(P_{i-1})\) in Equation 7, we have to bound the corresponding error by \(\epsilon _i\) such thatThen, the telescopic product in Equation (7) approximates \(\mbox{vol}(P)\) with error \(\le \epsilon\). In Section 5, we further discuss efficient error splitting. To estimate \(r_i\) we generate n billiard_walk points in \(P_{i-1}\) and we compute the proportion of the points that lie in \(P_i\). For the ratio \(r_{k+1} = \mbox{vol}(P_k)/\mbox{vol}(C_k)\) we follow the same procedure, but we sample from \(C_k\) and we count the number of points in \(P_k\). Computationally, we would like to determine an as small as possible integer n where we are within our target accuracy \(\epsilon _i\).
\begin{equation} {\sum _{i=1}^{k+1} \epsilon _i^2}=\epsilon ^2, \end{equation}
(11)
The main issue that one should address, is that volume generates correlated samples/points in the polytope P using billiard_walk. Therefore, the best known bounds for the number of billiard_walk points required are far too large for practical computations. Thus, to estimate \(r_i\), we use a sliding window to keep the last l estimation values of the ratio. That is a queue with length l; each time a new sample point is generated by inserting the new ratio value of \(\hat{r}_i\) and by popping out the oldest ratio value.
To determine an empirical stopping criterion, first, we assume that we have n i.i.d. uniformly distributed samples in \(P_{i-1}\). Then, and we employ the binomial confidence interval for the estimator of \(r_i\) to bound n. The number of points in \(P_i\) follows the binomial distribution \(b(n,r_i)\). A binomial proportion confidence interval for \(\hat{r}_i\)—the current approximation to \(r_i\)—is given by \(\hat{r}_i \pm z_{\alpha /2} \sqrt {\tfrac{\hat{r}_i(1-\hat{r}_i)}{n}}\), where \(z_{\alpha /2}\) is the \(1-\alpha /2\) quantile of the Gaussian distribution. As n increases, the interval tightens around \(\hat{r}_i\), thus, stopping for that value of n whereobtains an estimation of \(r_i\) within error \(\epsilon _i\) with probability \((1-\alpha)\), and Equation (7) would estimate \(\mbox{vol}(P)\) up to at most \(\epsilon\) with probability \((1 - \alpha)^{k+1}\).
\begin{equation*} \frac{z_{\alpha /2} \sqrt {\hat{r}_i(1-\hat{r}_i)/n}}{\hat{r}_i - z_{\alpha /2} \sqrt {\hat{r}_i(1-\hat{r}_i)/n}} \le \epsilon _i, \end{equation*}
Notice that \(\sqrt {\hat{r}_i(1-\hat{r}_i)/n}\) is an estimator of the st.d. of all sample fractions of size n. Thus, to define an empirical criterion we replace that quantity with the st.d., s, of the last l ratio values, i.e., the ratios stored in the sliding window. Moreover, we consider the average, \(\bar{r}\), of the ratios in the sliding window. Finally, we stop sampling when \(\bar{r}\) and the st.d. s meet the criterion of Equation (12): then we say they meet convergence. Clearly, for the first l points sampled, we do not check for convergence. For a probability p sufficiently close to 1, the empirical criterion for declaring convergence is as follows:
\begin{align} \frac{(\bar{r}+z_{p/2}s) - (\bar{r}-z_{p/2}s)}{\bar{r}-z_{p/2}s} = \frac{2z_{p/2}s}{\bar{r}-z_{p/2}s}\le \frac{\epsilon _i}{2}.\; \end{align}
(12)
In Section 5 we experimentally show that this empirical rule suffice to achieve error \(\epsilon\) for \(p=1-\sqrt [k+1]{3/4}\), and for \(l = O(1)\), i.e., when the length of the sliding window is constant—independent from the dimension.
5 Implementation and Experiments
In this section we discuss our implementation. We perform extended experiments analysing various aspects such as the tuning of the algorithm’s parameters and its practical complexity. Finally, we apply our software to compute the volume of high-dimensional Birkhoff polytopes, and we test the quality of approximation of various methods for low-order reduction of Z-polytopes.
We use the eigen library [29] for linear algebra and lpsolve [4] for solving linear programs. All experiments were performed on a PC with Intel® Core™ i7-6700 3.40GHz 8 CPU and 32GB RAM running Ubuntu 18. Runtimes reported in the plots and tables are averaged over 10 runs unless otherwise stated. We denote by CoolingGaussian the implementation of [17]. Brief instructions for how to run the implementation and reproduce the computational results of this article are described in the Appendix (Section B). All experiments are run without multi-threading.
It is of special interest to create and maintain a database of convex polytopes to evaluate the performance and the accuracy of various algorithms for sampling and volume approximation. Our polytope database, presented in Table 1, is constructed by merging and extending the databases used in [9, 17, 24].
Table 1.
| Polytope | Definition | H-rep. | V-rep. |
|---|---|---|---|
| cube-d | \(\lbrace x=(x_1,\dots ,x_d)\, |\, -1\le x_i\le 1, x_i\in {\mathbb {R}},\, i=1,\dots ,d \rbrace\) | \(\checkmark\) | \(\checkmark\) |
| cross-d | cross polytope, the dual of cube, i.e., \(\mbox{conv}(\lbrace -e_i,e_i : \, i=1,\dots ,d\rbrace)\) | \(\checkmark\) | \(\checkmark\) |
| \(\Delta\)-d | d-dimensional simplex \(\mbox{conv}(\lbrace e_i : \, i=1,\dots ,d\rbrace)\) | \(\checkmark\) | \(\checkmark\) |
| \(\Delta\)-d-d | the product of two simplices, i.e., \(\lbrace (p,p^{\prime })\in {\mathbb {R}}^{2d}\ |\ p\in \Delta -d,\ p^{\prime }\in \Delta -d \rbrace\), | \(\checkmark\) | |
| \(\mathcal {B}_n\) Birkhoff polytope | \(\lbrace x\in {\mathbb {R}}^d\ |\ x_{ij}\ge 0,\ \sum _i x_{ij}=1,\ \sum _j x_{ij}=1,\ 1\le i,j\le n \rbrace\) | \(\checkmark\) | |
| E-n-s Everest polytope | see [39] | \(\checkmark\) | |
| K-d dual Knapsack polytope | \(\mbox{conv}(\pm e_1, \ldots , \pm e_d, \mathbf {\alpha })\), where \(\mathbf {\alpha }\in \mathbb {Z}^d_{\ge 0}\) is generated randomly | \(\checkmark\) | |
| cc-8-k | 8d product of two 4d cyclic polyhedra | \(\checkmark\) | |
| ccp-k | complete cut polytope on k vertices | \(\checkmark\) | |
| rvc-d-n | polytope with n vertices in \([-1,1]^d\): generate uniform points in \([-1,1]^d\), stop when n lie in convex position | \(\checkmark\) | |
| rvs-d-n | polytope with n vertices uniformly distributed on the unit sphere | \(\checkmark\) | |
| rhs-d-n | polytope with n facets, with normals uniformly distributed on the unit sphere | \(\checkmark\) | |
| metabolic polytopes | polytopes that correspond to the flux space of a metabolic network [14, 30], see Table A.1 | \(\checkmark\) | |
| \(Z_{\mathcal {U}}\)-d-n | n generators, length of each generator selected uniformly from \([0,100]\) | ||
| \(Z_{\mathcal {N}}\)-d-n | n generators, length of each generator from \(\mathcal {N}(50,(50/3)^2)\) truncated to \([0,100]\) | ||
| \(Z_{Exp}\)-d-n | n generators, length of each generator from \(Exp(1/30)\) truncated to \([0,100]\) |
Table 1. The Polytope Database, Where d is the Dimension, \(e_1, \ldots , e_d\) are the Standard Basis Vectors in \({\mathbb {R}}^d\)
The direction of each generator for Z-polytopes: \(Z_{\mathcal {U}}\)-d-n, \(Z_{\mathcal {N}}\)-d-n, \(Z_{Exp}\)-d-n, is chosen uniformly from the \((d-1)\)-dimensional unit hypersphere.
For the special case of polytope that comes from a metabolic network, the polytope is given as a set of equalities and inequalities; that is a low-dimensional polytope. To obtain the full dimensional polytope \(P\subset {\mathbb {R}}^d\) we follow the preprocessing in [14] which we compute with library cobra [57].
The estimation of the telescopic product in Equation (7) might lead to numerical overflows and underflows in floating-point arithmetic. To avoid that our implementation estimates the following expression,
\begin{equation} \log (\mbox{vol}(P)) = \log (\mbox{vol}(C_k)) - \sum _{i=1}^k \log ({r_i}) +\log (r_{k+1}). \end{equation}
(13)
To evaluate the efficiency of our implementation, when the input polytope is a V- or a Z-polytope, we count the number of reflections (or boundary oracle calls) that billiard_walk performs to achieve a relative error \(\le \epsilon\). This is the number of Linear Programs (LP) in Equations (5, 6) that our implementation solves. When the input polytope is a H-polytope, reflections have different complexities i.e., the first takes \(\mathcal {O} (md)\) while the rest \(\mathcal {O} (m)\) operations. Therefore, we count the number of points that billiard_walk generates to achieve a relative error \(\le \epsilon\).
5.1 Parameter Tuning
Let us explain how we fine-tune the volume algorithm and the sampling procedures presented in Sections 2 and 4.
To start sampling from each body \(P_i=C_i\cap P\) in MMC of Equation (7) we use a central point of P. For H-polytopes we use the center of the Chebychev ball, i.e., the largest inscribed ball in P, which requires to solve a linear program [6]. For V-polytopes, we compute an approximation of the minimum-volume enclosing ellipsoid of the vertices [61] and use its center. For Z-polytopes we use the center of symmetry as they are centrally symmetric convex bodies.
5.1.1 Billiard Walk’s Parameters.
For billiard_walk first we have to set the parameter \(\tau\) that controls the length of the linear trajectory of each step of the random walk. In [56] they suggest to set \(\tau\) equal to the diameter of the body we are about to sample. Since computing the diameter could be hard, we are suggesting the following heuristics depending on the representation of the input convex body.
(1)
For the intersection of a polytope P, given in any representation, with a ball B, we set \(\tau\) equal to the diameter of B.
(2)
For H-polytopes computing the diameter is hard even by randomized algorithms [7]. Thus, we set \(\tau =4*\sqrt {d}r\), where r is the radius of the largest inscribed (Chebychev) ball in P.
(3)
For V-polytopes, the computation of the diameter is straightforward and takes \(O(dn^2)\) operations, where n is the number of vertices.
(4)
For Z-polytopes the diameter can be computed by a non-convex optimization problem:This optimization is NP-hard. However, we use the following heuristic: (a) compute the covariance matrix of the generators of P, (b) compute its eigenvector w that corresponds to the maximum eigenvalue, (c) compute \(\max \langle w,x \rangle ,\ x\in P\), in \(O(dn)\), where n is the number of generators.
\begin{equation} \min \ -\Vert x\Vert _2^2 = -\Vert Gy\Vert _2^2 = -y^TG^TGy, \text{ subject to: } y\in [-1,1]^n. \end{equation}
(14)
(5)
To compute the diameter of the intersection of a Z-polytope with a scaled copy of the polytope computed by Algorithm 2, i.e., \(Z-approx: =\lbrace x\in {\mathbb {R}}^d\ |\ Mx\le b \rbrace\), the program in Equation (14) becomes:In practice, we sample a set S of \(O(d^2)\) points from \([-1,1]^d\cap \Phi\), with CDHR, keeping the boundary points, and return \(2\max \limits _{x\in S} ||x||^2_2\).
\begin{equation} \begin{split}\min \ -\Vert x\Vert _2^2 = -\Vert Gy\Vert _2^2 = -y^TG^TGy, \text{ subject to: } y \in [-1,1]^d\cap \Phi ,\\ \text{ where }\Phi :=\lbrace y\in {\mathbb {R}}^d\ |\ MTy\le b\rbrace . \end{split} \end{equation}
(15)
Both heuristics for Z-polytopes and \(P\cap Z-approx\) return better solutions in our experiments than solvers like Sequential QP [42] from package NLopt for the optimization problems of Equations (14) and (15).
The second parameter of billiard_walk is the upper bound on the number of reflections \(\rho\). In our experiments, the average number of reflections in a billiard_walk’s step increases linearly or sub-linearly with the dimension (see Figures 8, 5, and 9). Thus, we tested upper bounds given by a linear function \(\sim \gamma d\) for some \(\gamma \in \mathbb {N}\). We choose to set \(\gamma = 20\) by following the experimental evidences in Figures 9 and 8 where \(\gamma \lt 16\) up to dimension 100.
Last, we have to set the walk length, which is the number of billiard_walk’s steps we burn until we use a point for ratio estimation. We aggressively set this to 1 following the observation in [17] that the fastest convergence of the empirical distribution happens for this particular value of walk length.
Interestingly, in our experiments, the number of points that billiard_walk with walk length equals to 1 generates to estimate each volume ratio, depends only on the number of bodies k in MMC and not on the dimension d. For example, in Figures 8 and 9 (down left plots), when volume fixes the same number of bodies in MMC, for different dimensions, billiard_walk generates the same number of points to achieve a relative error \(\le 0.1\). Consequently, the mixing rate of billiard_walk does not depend on the dimension on those instances.
5.1.2 Annealing Schedule’s Parameters.
For the cooling parameters in Section 4.1, we set \(r=0.1\) and \(\delta =0.05\) and thus each volume ratio would be restricted in \([r,r+\delta ]=[0.1,0.15]\) with high probability. We choose the significance level to be \(\alpha =0.10\). A smaller \(\alpha\) can be chosen for a tighter test around \(r+\delta\) which will result to more iterations in each step of annealing schedule. For larger \(\alpha\) values the number of iterations could be reduced but the method becomes unstable, as the probability that annealing_schedule fails to terminate increases while \(\alpha\) increases as well (see Theorem 2).
To set the parameters \(\nu , N\) we follow the empirical rule of Remark 4.1. Thus, assuming perfect uniform sampling and \(r_i=0.1\), one has to set \(N\ge 112\). At the initialization step of the annealing schedule we have to sample from a ball we set \(N=120\) as we employ perfect uniform sampling from a ball. Otherwise, when billiard_walk is used in all the other cases of uniform sampling from a convex body we set \(N=125\). We also set \(\nu =10\) to keep the product \(\nu N\) as small as possible. Smaller \(\nu\) values would result to more unstable iterations in annealing_schedule. Hence, following well known empirical rules we use t-student distribution in the t-tests. We apply the following optimization, for the i-th volume ratio, we sample from \(P_i\) and if the stopping criterion fails, we employ binary search by reusing the same sample. Hence, we sample only once the \(\nu N\) points per step of annealing_schedule.
annealing_schedule also has to compute \(q_{\min },q_{\max }\) such that \(q_{\min }C \subseteq P\subseteq q_{\max }C\). Ideally, the first would correspond to the largest inscribed scaled copy of the body C in the polytope P; the second to the minimally scaled copy of C that encloses P. In our implementation we set \(q_{\min }=0\) and to set \(q_{\max }\) we sample \(N\nu\) points and compute \(q_{\max }\) such that all points belong to \(q_{\max }C\). In the special case of C being the H-polytope of Algorithm 2 we compute the scaling factor \(q_{\max }\) to compute the polytope \(q_{\max }C\) that tightly encloses the Z-polytope P and we set \(q_{\min }=0\).
5.1.3 Error Splitting and Sliding Window.
Recall that when we estimate each ratio in Multiphase Monte Carlo (MMC), we have to do an error splitting for the volume ratios that appear in the telescopic product in Equation (7)—i.e., given a requested error \(\epsilon\) for the algorithm compute an error \(\epsilon _i\) for each volume ratio \(r_i\) (see Section 4.3), such that,
\begin{equation*} \sum _{i=1}^{k+1}\epsilon _i^2 = \epsilon ^2. \end{equation*}
We do not split the error \(\epsilon\) equally to all ratios. In particular, we set the requested error \(\epsilon _{k+1}=\epsilon /2\sqrt {k+1}\) for the last ratio \(r_{k+1} = \mbox{vol}(P_k)/\mbox{vol}(C_k)\) to be the smallest ratio in the telescopic product, where k is the number of bodies in MMC and the rest errors \(\epsilon _i,\ i=1,\dots , k\) are set to be equally weighted. The reason is that \(\mbox{vol}(P_k)/\mbox{vol}(C_k)\) converges faster than the other ratios in practice. The latter occurs because sampling from \(C_k\) is usually faster and more accurate than sampling from any other \(P_i\), e.g. when the body C used in MMC is either a ball while \(P_i\) is the intersection of P with a ball, or an H-polytope while P is a Z-polytope.
Then we have to split \(\epsilon ^{\prime }=\epsilon \sqrt {4(k+1)-1}/2\sqrt {k+1}\) to the remaining ratios, since \(\epsilon _{k+1}^2 + \epsilon ^{\prime 2} = \epsilon ^2\). Thus, we set \(\epsilon _i=\epsilon ^{\prime }/\sqrt {k},\ i=1, \ldots , k\) so that \({\sum _{i=0}^{k+1} \epsilon _i^2}=\epsilon ^2\) holds.
In the special case of P being a Z-polytope, the body \(C_k\) is an H-polytope. We estimate \(\mbox{vol}(C_k)\) by calling volume using balls in MMC, while the computational time of \(\mbox{vol}(C_k)\) is a small portion of the total runtime. For the computation of \(\mbox{vol}(C_k)\) we set the requested error \(e^{\prime \prime }=\epsilon /2\sqrt {k+1}\) and then we equally split \(\epsilon ^{\prime } =\epsilon \sqrt {2k+1}/\sqrt {2k+2}\) to the \(k+1\) ratios respecting \({\sum _{i=0}^{k+1} \epsilon _i^2}=\epsilon ^2\).
For the length of the sliding window, we set \(l=250\). Our experiments on the error show this choice offers stability, e.g. in Figure 2 volume always compute a smaller error than the requested error 0.1 for unit cubes, unit simplices, and cross polytopes never exceeds the requested value for \(d\le 1000\).
Fig. 2.
For each new generated point, we update the average mean and the variance of the sliding window in \(O(1)\) as follows. Let \(\hat{\mu }\) be the average mean of the ratios in the sliding window, then we write the variance asWe store the sum of the window’s ratios \(\sum _{i=1}^{l}\hat{r}_i\) and the sum of the squared ratios \(\sum _{i=1}^{l}\hat{r}_i^2\). For each new generated point we obtain an updated ratio and the oldest ratio is popped out. We use both the updated and the popped out ratios to update both \(\sum _{i=1}^{l}\hat{r}_i^2\) and \(\sum _{i=1}^{l}\hat{r}_i\) and to compute the updated average mean value and st.d. of the current ratios in the sliding window.
\begin{equation*} \frac{1}{l}\sum _{i=1}^{l}(\hat{r}_i-\hat{\mu })^2=\frac{1}{l}\left(\sum _{i=1}^{l}\hat{r}_i^2-2\hat{\mu }\sum _{i=1}^{l}\hat{r}_i+l\hat{\mu }^2\right). \end{equation*}
In the sequel, we benchmark the proposed algorithm volume tuned as described in the previous sections. We present the result obtained from our experiments for each polytope representation separately.
5.2 Volume Experiments for H-polytopes
The implementation of our algorithm, volume, scales up to thousands of dimensions within a few hours. In Tables A.1 and A.2 we report the runtime of volume for cube-d, \(\Delta\)-d, \(\Delta\)-d-d, rhs-m-d, metabolic polytopes and \(\mathcal {B}_n\), while volume achieves relative error equal at most 0.1. Except of cube-d, \(\Delta\)-d, \(\Delta\)-d-d we apply the rounding preprocess on Section 3 before volume estimation. Our implementation takes for \(d = 100\) a few seconds and for \(d = 500\) a few minutes except for iSDY_1059 that takes almost an hour as it has a larger number of facets (i.e., \(m\approx 3000\)). When \(d=1,000\) our implementation takes a few hours, while the runtime increases with the number of facets. However, the runtime for cube-1000 is smaller than the \(\Delta\)-1000 as the unit simplex has a larger isotropic constant than the unit hypercube; consequently the number of bodies in MMC is larger for the unit simplex. Our implementation is the first one that estimates the volume of very high dimensional Birkhoff polytopes extending the computational results of [17, 24]. In particular, the \(\mathcal {B}_n\) of order \(16\le n\le 33\) with dimension \(225\le d\le 1024\) are computed in \(\le 4\)hr.
To evaluate the performance of our implementation of volume we count the number of points that billiard_walk generates for the cube-d, \(\Delta\)-d, rhs-d-m and \(\mathcal {B}_n\). First, for the number of bodies in MMC we notice in Figure 3 that it grows as \(O(d\lg d)\) which agrees with the analysis in Proposition 4. Then, due to the error splitting of Section 5.1.3 the generated points per volume ratio grows as \(O^*(d)\). Thus, the total number of billiard_walk points grows as \(O^*(d^2)\) for those polytopes (left plot in Figure 3). Thus, the total run-time of our implementation grows as \(O(d^3m)\), because the per step cost of billiard_walk is \(O(md)\).
Fig. 3.
For the rounding method of Section 3 we experimentally evaluate its quality by counting the number of generated bodies in MMC. In particular, given a skinny polytope P we round it and then we compare the number of bodies in MMC with that of isotropic polytopes (e.g. cube-d, \(\Delta\)-d-d). In Figure 3 (right) for rhs-d-m and \(\mathcal {B}_n\) we notice that for both rounded and isotropic bodies the number of bodies grows as \(O(d\lg d)\) which is a strong evidence that our rounding method transforms the polytope to a near isotropic position.
Next, we experimentally study how the complexity of our implementation depends on the input error parameter. In Figure 4 we estimate the volumes of cube-d for \(d = 50,100,\dots , 500\) for two different values for the error parameter \(\epsilon = 0.1,\ 0.2\). In the left plot we notice that for both values the number of generated points by volume increase as \(O(d^2)\), while a smaller value of \(\epsilon\) just increases the slope of the line, i.e., increases the run-time by a constant. The right plot also confirms this observation; it shows that the run-time of volume increases as \(O(d^4)\) for cube-d. This is expected as the number of points increases as \(O(d^2)\) and the per step cost of billiard_walk is \(O(md)\), while \(m=O(d)\) for cube-d.
Fig. 4.
Finally, we examine the time spent in different steps in the algorithm i.e., rounding, constructing the sequence, and volume computation. We conclude that those are highly depended on the convex body. For example for a 100-cube the times for the sequence construction and volume computation is 2.06836 and 4.29843 secs, respectively. While for a 100-simplex those numbers are 2.85144 and 17.2575. Intuitively, this effect is related to the isotropy of the convex set. Regarding preprocessing/rounding this also depends on the geometry of the body. For example, the rounding of the unit 100-cube takes 0.01427 secs while for a skinny 100-cube it takes 6.92041 secs.
5.2.1 Comparison with Cooling Gaussian.
We compare our implementation, i.e., volume, against the state-of-the-art MATLAB implementation of CoolingGaussian in [17] for cube-d, simplices \(\Delta\)-d and random polytopes rhs-d-m with \(m=5d\). In Figure 5, we report the total number of generated points for both algorithms to achieve a relative error \(\le 0.1\). CoolingGaussian performs computations up to \(d=500\) while volume compute up to \(d=1000\) in a shorter time frame. More interestingly, the ratio between the number of points generated by CoolingGaussian over the number of points generated by volume increases linearly with the dimension. Then, the comparison of the run-times is determined by the per step cost of the random walks used in each implementation.
Fig. 5.
The plot in Figure 6 compares the run-times of both implementations for cube-d, \(\Delta\)-d in H-representation (left plot), and rhs-d-m with \(m=5d\) and \(d=50,100, \ldots , 1000\). It confirms that volume is faster than CoolingGaussian in [17]. Moreover, notice that the gap on the run-time increases with the dimension as expected. In particular, for \(d=50\) our implementation is 20, 2, and 3 times faster for cube-d, \(\Delta\)-d and rhs-d-\(m,\) respectively. For \(d=500\) our implementation is 95, 8, and 10 times faster for the same classes of polytopes. For \(d=1000\) we use polynomial interpolation to estimate the runtime of CoolingGaussian, and our implementation is expected to be around 130, 13, and 20 times faster for cube-d, \(\Delta\)-d and rhs-d-\(m,\) respectively.
Fig. 6.
5.3 Volume Experiments for Z-polytopes
Our implementation is the first one that scales up to hundreds of dimensions in order of hours for Z-polytopes. In Table A.3 we report the runtimes for both high and low order Z-polytopes. Our implementation, for a 100-dimensional Z-polytope of order equal to 2 performs \(6.37\cdot 10^4\) reflections (equivalent to boundary oracle calls in [17]) in less than an hour. Considering high order Z-polytopes our implementation scales up to order 150 for \(d\le 15\) and to order 100 and 20 for \(d=20\) and \(d=30,\) respectively, in at most half an hour by performing \(5.56\cdot 10^3\) reflections.
The only alternative implementation for Z-polytopes is that in [17], however they estimate volumes only for low dimensional Z-polytopes, while our implementation is undoubtedly superior. Since, we were unable to reproduce the results reported on [17] for Z-polytopes we compare against the data given in [17]. In particular, for a 2-order, 10-dimensional Z-polytope our implementation performs \(9.79\cdot 10^3\) reflections (boundary oracle calls) in 2.7 seconds to achieve a relative error smaller than 0.1; the implementation of CoolingGaussian in [17] generates \(9.58\cdot 10^4\) CDHR points which corresponds to \(1.92\cdot 10^5\) boundary oracle calls in 3010 seconds to achieve a relative error higher than 0.1 and smaller than 0.2. For a 100-order, 10-dimensional Z-polytope our implementation performs \(2.95\cdot 10^3\) reflections in 180 seconds while the CoolingGaussian in [17] performs \(2.50\cdot 10^5\) boundary oracle calls in 4,900 seconds for the same values of relative errors. Clearly, for those instances, our implementation achieves a better accuracy, while the run-time is at least two order of magnitudes smaller than that of CoolingGaussian in [17].
Exact volume computation of Z-polytopes [28] require an exponential to the dimension d number of operations. Thus, it cannot scale beyond \(d = 10\) for low order Z-polytopes.
In the second column of Table A.3 we report the type of body we use in MMC. For low order Z-polytopes we use the H-polytope from Section 3 while for high order Z-polytopes we use the unit ball. Our experiments suggest to set the threshold for the order \(n/d\le 5\) so that the Z-polytope is considered as a low order one. In particular, Figure 7 shows that for order \(\le 4\), if we use the H-polytope in MMC, we get a bound on the the number of bodies in MMC \(k\le 3\) for \(d\le 100\). Note that when we use the H-polytope the number of bodies k is smaller for all pairs \((d,n)\) compared to using the unit ball without applying rounding to the Z-polytope. When we use balls in MMC, n decreases for constant d as n increases. For orders equal to 5 the number of balls in MMC, without rounding, is equal or smaller than the number of H-polytopes in MMC when we use the H-polytope of Section 3. Table A.3 shows that, for high-order Z-polytopes, \(k=1\), which implies one or two acceptance-rejection steps. Moreover, Table A.5 illustrates that the rounding method combined with balls in MMC results to a larger number of bodies in MMC and runtime than the case of using the H-polytope in MMC, while for order equal to 5 the best runtime occurs when using a ball in MMC without rounding. Thus, we use the H-polytope in MMC if the order \(n/d \lt 5\). We conclude that for a Z-polytope of any order, no rounding preprocessing is needed, if we make the right choice of C depending on the order; the maximum number of bodies is \(k\le 3\) for \(d\le 100\).
Fig. 7.
To evaluate the efficiency of our implementation for Z-polytopes we wish to count the average number of reflections. We run our implementation of volume for \(Z_{\mathcal {U}}\)-d-n, \(Z_{\mathcal {N}}\)-d-n and \(Z_{Exp}\)-d-n Z-polytopes. An example for Z-polytopes of order \(n/d=2\) is illustrated in Figure 8. First, as mentioned above our experiments imply that the number of bodies in MMC is \(\le 3\) for any order for \(d\le 100\). Moreover, the average number of reflections per point grows sub-linearly in d (down-right plot in Figure 8). Considering that billiard_walk generates \(O^*(1)\) points per volume ratio, the total number of reflections grows sub-linearly in d.
Fig. 8.
Application: Evaluate Z-polytope approximation.
We propose an efficient algorithm for evaluating an over-approximation of a given Z-polytope P. Z-polytopes are critical in applications such as autonomous driving [2] or human-robot collaboration [55]. Complexity strongly depends on the order of the encountered Z-polytopes. Thus, a practical solution is to over-approximate P, as tight as possible, with another Z-polytope \(P_{red}\supseteq P\) of smaller order. A good measure of the approximation quality (fitness) isThus, to compute the approximation quality we need to compute volumes of Z-polytopes. In [41] they compute volumes exactly and deterministically, therefore it is impossible to compute the quality of approximation for \(d\gt 10\).
\begin{equation} R= (\mbox{vol}(P_{red})/\mbox{vol}(P))^{1/d} \end{equation}
(16)
Here, we exploit our software to test the quality of such approximations. The two methods that are able to scale for \(d\ge 20\) are, primarily, the Principal Component Analysis (PCA) and the BOX method from [41]. Both adopt similar approximations and are of comparable reliability; here we focus on PCA: Let \(X= [G | -G]^T\) and \(\text{SVD}(X^T X)=U\Sigma V^T\), then \(G_{red}=U\cdot IH(U^T G)\) is square and generates \(P_{red}\), where \(IH(\cdot)\) is the “interval hull” from [43]. Over-approximation can be seen as a reduction problem, so that the covariance among the d generators of \(P_{red}\) must be null.
Table A.3 shows experimental results for Z-polytopes up to \(d=100\) and various orders. \(\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. To our knowledge this is the first time practical volume estimation is used to test approximation methods in that high dimensional spaces.
5.4 Volume Experiments for V-polytopes
In Table A.4 we report the runtimes of our implementation for several V-polytopes. It estimates volume within a relative error at most 0.1 in less than an hour. To round a V-polytope, we introduce a new method to control the sandwiching ratio \(R/r\). Since the vertices are known, we put P to John’s position [35]. First, we compute the smallest enclosing ellipsoid E of P’s vertices by applying the approximation algorithm in [61]. Notice that E is also the smallest enclosing ellipsoid of P and, if E is a ball, then P is in John’s position. To round P, we apply to it the transformation that maps E to the unit ball and repeat until the ratio between E’s longest and shortest axes falls under a given threshold. Our experiments suggest this is much faster than the method in Section 3, as sampling is costly for this representation due to expensive boundary oracle calls. In the left plot of Figure A.1, notice that, for rounded V-polytopes, the number of bodies in MMC k increases as a linear function of d. Moreover, as the number of vertices increases for constant d, then k decreases. All rounding steps in our experiments for \(d\le 100\) took \(\le 3\) seconds.
To evaluate the performance of our algorithm we wish to count the average number of reflections (or boundary oracle calls). We run volume for rvc-d-n and rvs-d-n polytopes after a rounding step. An example of V-polytopes with number of vertices \(n=1.5d,\ 2d,\ 3d\) is illustrated in Figure 9: The average number of reflections grows sublinearly in d (down-right plot in Figure 9). The total number of oracle calls grows as \(O^*(d)\) (up-left plot in Figure 9). volume performs similarly for both rvc-d-n and rvs-d-n polytopes. However, the number of boundary oracle calls decreases as the number of vertices n increases for constant d. This holds because after a rounding step the more the vertices of a random V-polytope the smaller the sandwiching ratio.
Fig. 9.
6 Conclusions and Future Work
We propose a new practical algorithm that computes for the first time volumes of very high-dimensional H-polytopes i.e., in the order of thousands of dimensions and high-dimensional V-, and Z- polytopes i.e., in the order of hundreds of dimensions. We provide strong empirical evidences on the efficiency of the algorithm and its accuracy through extended experiments on benchmark polytopes in computational geometry and bioinformatics.
We expect our software to address hard counting problems, e.g. counting the number of linear extensions of a partially ordered set [60]. Our MMC scheme can be extended for other hard computational problems as Markov Chain Monte Carlo integration for multivariate integrals over a convex polytope. In particular, one could estimate the volume and generate the uniform samples required by one pass. Last but not least, we could exploit parallelism and faster software libraries for linear programs to scale to higher dimensions for V- and Z-polytopes.
Footnotes
1
i.e., Markov Chain Monte Carlo (MCMC) algorithms to sample from multivariate distributions constrained in (non-)convex bodies.
2
A recent exception for sampling is the implementation of [40].
A Figures and Tables
Table A.1.
| H-polytope | volume | std/mean | error | k | points | time |
|---|---|---|---|---|---|---|
| cross-10 | 2.81e-04 | 0.06 | 0.01 | 1.0 | 1.50e+03 | 0.11 |
| cross-13 | 1.30e-06 | 0.07 | 0.01 | 1.0 | 1.50e+03 | 1.17 |
| cross-15 | 2.49e-08 | 0.04 | 0.01 | 1.0 | 1.50e+03 | 12.8 |
| cube-100 | 1.18e+30 | 0.15 | 0.07 | 9.2 | 2.80e+04 | 1.5 |
| cube-250 | 1.96e+75 | 0.16 | 0.08 | 34.5 | 1.96e+05 | 36.2 |
| cube-500 | 3.38e+150 | 0.25 | 0.03 | 84.7 | 7.09e+05 | 425 |
| cube-750 | 5.75e+225 | 0.22 | 0.03 | 139.9 | 1.52e+06 | 1501 |
| cube-1000 | 2.11e+301 | 0.29 | 0.03 | 201.0 | 2.70e+06 | 6180 |
| \(\Delta\)-100 | 1.02e-158 | 0.16 | 0.05 | 21.0 | 1.01e+05 | 4.8 |
| \(\Delta\)-250 | 3.20e-493 | 0.15 | 0.03 | 66.3 | 5.05e+05 | 66.4 |
| \(\Delta\)-500 | 7.55e-1135 | 0.20 | 0.08 | 149.1 | 1.65e+06 | 847.8 |
| \(\Delta\)-750 | 3.70e-1833 | 0.17 | 0.04 | 237.9 | 3.28e+06 | 4088 |
| \(\Delta\)-1000 | 2.67e-2568 | 0.24 | 0.07 | 328.4 | 5.34e+06 | 10711 |
| \(\Delta\)-50-50 | 1.07e-129 | 0.21 | 0.01 | 22.3 | 9.90e+04 | 3.3 |
| \(\Delta\)-125-125 | 2.91e-419 | 0.18 | 0.03 | 67.0 | 5.07e+05 | 107.2 |
| \(\Delta\)-250-250 | 9.10e-986 | 0.23 | 0.05 | 151.3 | 1.69e+06 | 958.5 |
| \(\Delta\)-375-375 | 6.26e-1609 | 0.21 | 0.06 | 239.3 | 3.35e+06 | 4344 |
| \(\Delta\)-500-500 | 6.22e-2269 | 0.24 | 0.07 | 332.5 | 5.39e+06 | 9432 |
| rhs-100-1000 | 4.19e-46 | 0.11 | ?? | 10.6 | 3.33e+04 | 5.25 |
| rhs-300-2000 | 2.75e-205 | 0.20 | ?? | 50.0 | 3.22e+05 | 253 |
| rhs-1000-6000 | 4.18e+275 | 0.23 | ?? | 202.1 | 2.61e+06 | 19431 |
| rhs-1000-10000 | 1.47e+216 | 0.21 | ?? | 172.3 | 1.90e+06 | 22833 |
| e_coli_core [174-24] | 3.37e-47 | 0.08 | ?? | 2.0 | 3.52e+03 | 0.16 |
| iAB_RBC_283 [906-130] | 1.73e+43 | 0.28 | ?? | 23.9 | 1.11e+05 | 31.3 |
| iJR904 [1132-227] | 6.00e+53 | 0.25 | ?? | 50.8 | 3.38e+05 | 169.3 |
| iAT_PLT_636 [2016-289] | 3.65e-854 | 0.25 | ?? | 66.3 | 4.99e+05 | 512.7 |
| iSDY_1059 [2966-509] | 2.25e-350 | 0.27 | ?? | 136.0 | 1.43e+06 | 3486 |
| Recon1 [4934-931] | 8.01e-5321 | 0.17 | ?? | 266.1 | 3.84e+06 | 18225 |
Table A.1. Volume Estimation for H-polytopes by volume
For each polytope volume performs 10 runs; volume: mean volume computed; std/mean: ratio of std and mean of volumes; error: the relative error; k: average number of generated bodies in MMC; points: average number of generated points; time: average time in seconds; \(\epsilon = 0.1\) used in all cases; ??: the exact volume is unknown. The last six are metabolic polytopes with m facets in d dimensions denoted as \([m-d]\). For metabolic polytopes we follow the preprocessing of cobra package [57] also used in [14, 30].
Table A.2.
| \(\mathcal {B}_n\) | d | exact | asymptotic | volume | error | std/mean | k | points | time |
|---|---|---|---|---|---|---|---|---|---|
| \(\mathcal {B}_3\) | 4 | 1.12 | 1.410849 | 1.12e+00 | 0.00 | 0.02 | 1 | 1.50e+03 | 0.01 |
| \(\mathcal {B}_4\) | 9 | 6.21e-02 | 7.61e-02 | 6.17e-02 | 0.01 | 0.02 | 1 | 1.50e+03 | 0.01 |
| \(\mathcal {B}_5\) | 16 | 1.41e-04 | 1.69e-04 | 1.38e-04 | 0.02 | 0.01 | 1 | 1.52e+03 | 0.03 |
| \(\mathcal {B}_6\) | 25 | 7.35e-09 | 8.62e-09 | 7.35e-09 | 0.00 | 0.05 | 2.3 | 3.63e+03 | 0.06 |
| \(\mathcal {B}_7\) | 36 | 5.64e-15 | 6.51e-15 | 5.49e-15 | 0.03 | 0.02 | 4.0 | 7.98e+03 | 0.16 |
| \(\mathcal {B}_8\) | 49 | 4.42e-23 | 5.03e-23 | 4.24e-23 | 0.04 | 0.05 | 6.4 | 1.62e+04 | 0.37 |
| \(\mathcal {B}_9\) | 64 | 2.60e-33 | 2.93e-33 | 2.54e-33 | 0.02 | 0.03 | 9.2 | 2.82e+04 | 0.86 |
| \(\mathcal {B}_{10}\) | 81 | 8.78e-46 | 9.81e-46 | 8.46e-46 | 0.04 | 0.05 | 13.1 | 4.69e+04 | 1.9 |
| \(\mathcal {B}_{11}\) | 100 | ?? | 1.49e-60 | 1.24e-60 | ?? | 0.05 | 17.6 | 7.17e+04 | 3.5 |
| \(\mathcal {B}_{12}\) | 121 | ?? | 8.38e-78 | 7.28e-78 | ?? | 0.10 | 23.1 | 1.07e+05 | 6.7 |
| \(\mathcal {B}_{13}\) | 144 | ?? | 1.43e-97 | 1.205e-97 | ?? | 0.19 | 29.2 | 1.49e+05 | 11.9 |
| \(\mathcal {B}_{14}\) | 169 | ?? | 6.24e-120 | 5.13e-120 | ?? | 0.08 | 35.9 | 2.04e+05 | 20.7 |
| \(\mathcal {B}_{15}\) | 196 | ?? | 5.94e-145 | 4.92e-145 | ?? | 0.05 | 43.4 | 2.68e+05 | 34.5 |
| \(\mathcal {B}_{16}\) | 225 | ?? | 1.06e-172 | 9.62e-173 | ?? | 0.12 | 51.9 | 3.49e+05 | 56.3 |
| \(\mathcal {B}_{17}\) | 256 | ?? | 3.10e-203 | 2.50e-203 | ?? | 0.10 | 60.9 | 4.41e+05 | 84.2 |
| \(\mathcal {B}_{18}\) | 289 | ?? | 1.30e-236 | 1.03e-237 | ?? | 0.10 | 72.0 | 5.65e+05 | 129 |
| \(\mathcal {B}_{19}\) | 324 | ?? | 6.96e-273 | 5.41e-273 | ?? | 0.09 | 82.0 | 6.88e+05 | 198 |
| \(\mathcal {B}_{20}\) | 361 | ?? | 4.23e-312 | 3.31e-312 | ?? | 0.11 | 94.9 | 8.46e+05 | 297 |
| \(\mathcal {B}_{21}\) | 400 | ?? | 2.62e-354 | 1.77e-354 | ?? | 0.27 | 101.6 | 9.40e+05 | 325 |
| \(\mathcal {B}_{22}\) | 441 | ?? | 1.49e-399 | 1.49e-399 | ?? | 0.17 | 115.0 | 1.12e+06 | 453 |
| \(\mathcal {B}_{23}\) | 484 | ?? | 7.10e-448 | 4.33e-448 | ?? | 0.30 | 128.8 | 1.33e+06 | 624 |
| \(\mathcal {B}_{24}\) | 529 | ?? | 2.57e-499 | 2.10e-499 | ?? | 0.27 | 143.7 | 1.55e+06 | 845 |
| \(\mathcal {B}_{25}\) | 576 | ?? | 6.46e-554 | 5.56e-554 | ?? | 0.28 | 158.8 | 1.82e+06 | 1060 |
| \(\mathcal {B}_{26}\) | 625 | ?? | 1.04e-611 | 7.80e-612 | ?? | 0.22 | 176.3 | 2.11e+06 | 1405 |
| \(\mathcal {B}_{27}\) | 676 | ?? | 9.81e-673 | 1.00e-672 | ?? | 0.30 | 194.0 | 2.41e+06 | 1918 |
| \(\mathcal {B}_{28}\) | 729 | ?? | 5.05e-737 | 6.13e-737 | ?? | 0.31 | 211.1 | 2.75e+06 | 2573 |
| \(\mathcal {B}_{29}\) | 784 | ?? | 1.31e-804 | 1.17e-804 | ?? | 0.31 | 231.5 | 3.15e+06 | 3539 |
| \(\mathcal {B}_{30}\) | 841 | ?? | 1.60e-875 | 1.35e-875 | ?? | 0.35 | 249.1 | 3.53e+06 | 4982 |
| \(\mathcal {B}_{31}\) | 900 | ?? | 8.55e-950 | 8.43e-950 | ?? | 0.40 | 272.9 | 4.02e+06 | 6967 |
| \(\mathcal {B}_{32}\) | 961 | ?? | 1.86e-1027 | 1.58e-1027 | ?? | 0.39 | 293.8 | 4.47e+06 | 8905 |
| \(\mathcal {B}_{33}\) | 1024 | ?? | 1.56e-1108 | 1.37e-1108 | ?? | 0.42 | 314.1 | 4.97e+06 | 11105 |
Table A.3.
| Z-polytope | Body | order | volume | std/mean | k | Refl | time | \(\mbox{vol}(P_{red})\) | R |
|---|---|---|---|---|---|---|---|---|---|
| \(Z_{\mathcal {U}}\)-10-200 | Ball | 20 | 5.47e+34 | 0.06 | 1 | 3.15e+03 | 13 | 2.24e+37 | 1.82 |
| \(Z_{\mathcal {U}}\)-10-500 | Ball | 50 | 3.19e+38 | 0.07 | 1 | 3.00e+03 | 50 | 1.42e+41 | 1.84 |
| \(Z_{\mathcal {U}}\)-10-1000 | Ball | 100 | 2.57e+41 | 0.04 | 1 | 2.95e+03 | 180 | 1.08e+44 | 1.83 |
| \(Z_{\mathcal {U}}\)-10-1500 | Ball | 150 | 1.85e+43 | 0.05 | 1 | 2.88e+03 | 320 | 7.30e+45 | 1.82 |
| \(Z_{\mathcal {U}}\)-15-300 | Ball | 20 | 3.79e+51 | 0.07 | 1 | 4.34e+03 | 50 | 3.63e+56 | 2.15 |
| \(Z_{\mathcal {U}}\)-15-750 | Ball | 50 | 8.14e+57 | 0.05 | 1 | 3.39e+03 | 183 | 7.84e+62 | 2.15 |
| \(Z_{\mathcal {U}}\)-15-1500 | Ball | 100 | 3.45e+62 | 0.07 | 1 | 3.22e+03 | 503 | 3.20e+62 | 2.14 |
| \(Z_{\mathcal {U}}\)-15-2250 | Ball | 150 | 1.47e+65 | 0.05 | 1 | 3.14e+03 | 1071 | 1.31e+70 | 2.14 |
| \(Z_{\mathcal {U}}\)-20-300 | Ball | 15 | 1.31e+67 | 0.09 | 1 | 4.53e+03 | 72 | 7.91e+74 | 2.45 |
| \(Z_{\mathcal {U}}\)-20-400 | Ball | 20 | 5.58e+69 | 0.05 | 1 | 4.46e+03 | 110 | 2.90e+77 | 2.43 |
| \(Z_{\mathcal {U}}\)-20-1000 | Ball | 50 | 1.12e+77 | 0.05 | 1 | 3.76e+03 | 455 | 4.91e+84 | 2.41 |
| \(Z_{\mathcal {U}}\)-20-2000 | Ball | 100 | 1.93e+83 | 0.06 | 1 | 3.51e+03 | 1736 | 7.56e+90 | 2.40 |
| \(Z_{\mathcal {U}}\)-30-450 | Ball | 15 | 8.15e+99 | 0.10 | 1 | 5.45e+03 | 280 | 5.73e+113 | 2.89 |
| \(Z_{\mathcal {U}}\)-30-600 | Ball | 20 | 1.49e+104 | 0.06 | 1 | 5.46e+03 | 503 | 1.04e+118 | 2.89 |
| \(Z_{Exp}\)-60-78 | Hpoly | 1.3 | 1.47e+16 | 0.07 | 1 | 7.26e+04 | 357 | 9.25e+54 | 4.43 |
| \(Z_{Exp}\)-80-104 | Hpoly | 1.3 | 6.39e+08 | 0.08 | 3.0 | 1.37e+05 | 1604 | 5.03e+71 | 6.11 |
| \(Z_{Exp}\)-100-130 | Hpoly | 1.3 | 4.73e+17 | 0.05 | 3.0 | 1.97e+05 | 2851 | 1.01e+101 | 6.81 |
| \(Z_{Exp}\)-60-90 | Hpoly | 1.5 | 9.24e+19 | 0.07 | 1 | 5.68e+04 | 332 | 1.04e+58 | 4.31 |
| \(Z_{Exp}\)-80-120 | Hpoly | 1.5 | 3.46e+21 | 0.05 | 2.7 | 7.28e+04 | 901 | 9.70e+75 | 4.79 |
| \(Z_{Exp}\)-100-150 | Hpoly | 1.5 | 5.30e+37 | 0.09 | 3 | 1.27e+05 | 2647 | 2.66e+108 | 5.09 |
| \(Z_{\mathcal {N}}\)-60-102 | Hpoly | 1.7 | 2.06e+132 | 0.07 | 2 | 4.18e+04 | 329 | 2.65e+175 | 5.23 |
| \(Z_{\mathcal {N}}\)-80-136 | Hpoly | 1.7 | 1.74e+176 | 0.10 | 2.8 | 5.17e+04 | 857 | 2.63e+237 | 5.81 |
| \(Z_{\mathcal {N}}\)-100-170 | Hpoly | 1.7 | 3.74e+222 | 0.14 | 3 | 8.06e+04 | 2402 | 1.98e+303 | 6.42 |
| \(Z_{\mathcal {U}}\)-60-120 | Hpoly | 2 | 1.85e+135 | 0.12 | 2 | 3.38e+04 | 337 | 1.89e+177 | 5.01 |
| \(Z_{\mathcal {U}}\)-80-160 | Hpoly | 2 | 2.27e+186 | 0.15 | 2.6 | 5.72+04 | 1227 | 2.56e+245 | 5.47 |
| \(Z_{\mathcal {U}}\)-100-200 | Hpoly | 2 | 5.86e+231 | 0.17 | 3 | 6.37e+04 | 2652 | 1.85e+310 | 6.10 |
| \(Z_{\mathcal {U}}\)-60-180 | Hpoly | 3 | 2.95e+155 | 0.21 | 2 | 1.98e+04 | 417 | 6.61e+191 | 4.03 |
| \(Z_{\mathcal {U}}\)-80-240 | Hpoly | 3 | 5.59e+206 | 0.19 | 2.8 | 2.97e+04 | 2497 | 9.47e+260 | 4.76 |
Table A.3. Volume Estimation for Z-polytopes by volume
For each polytope volume performs 10 runs. Body: type of body in MMC; order: the ratio between the number of generators over the dimension (\(n/d\)); volume: the mean computed volume; error: relative error; k: average number of bodies in MMC; Refl: average number of reflections (or boundary oracle calls) performed by billiard_walk; time: average time in seconds; \(\mbox{vol}(P_{red})\): volume of the over-approximation; R: ratio of fitness by Equation (16); \(\epsilon = 0.1\) in all cases.
Table A.4.
| V-polytope | volume | std/mean | k | Refl | error | time | exact time |
|---|---|---|---|---|---|---|---|
| cross-50 | 4.546e-13 | 0.05 | 1.0 | 10.9e+03 | 0.05 | 3 | \(--\) |
| cross-100 | 1.43e-128 | 0.05 | 2.0 | 5.80e+04 | 0.08 | 112 | \(--\) |
| cross-150 | 1.50e-64 | 0.07 | 2.0 | 2.99e+04 | 0.08 | 24 | \(--\) |
| cross-200 | 1.65e-95 | 0.10 | 3.0 | 3.92e+04 | 0.02 | 53 | \(--\) |
| \(\Delta\)-40 | 1.17e-48 | 0.15 | 6.1 | 2.42e+05 | 0.05 | 230 | 0.008 |
| \(\Delta\)-60 | 1.12e-82 | 0.14 | 11.2 | 5.13e+05 | 0.07 | 989 | 0.02 |
| \(\Delta\)-80 | 1.35e-119 | 0.21 | 15.6 | 1.21e+06 | 0.03 | 4140 | 0.07 |
| cube-10 | 1052.4 | 0.07 | 1 | 1.50e+03 | 0.03 | 54 | \(--\) |
| cube-11 | 1930.2 | 0.06 | 1 | 1.5e+03 | 0.06 | 155 | \(--\) |
| cube-12 | 4240.6 | 0.08 | 1 | 1.81e+03 | 0.04 | 567 | \(--\) |
| cube-13 | 7538.2 | 0.06 | 1 | 1.87e+03 | 0.08 | 2937 | \(--\) |
| K-60 | 2.43e-62 | 0.10 | 1.0 | 4.08e+04 | ?? | 178 | \(--\) |
| K-70 | 2.05e-77 | 0.12 | 1.0 | 4.78e+04 | ?? | 298 | \(--\) |
| K-80 | 3.62e-93 | 0.14 | 1.0 | 5.43e+04 | ?? | 466 | \(--\) |
| K-90 | 1.72e-109 | 0.08 | 1.9 | 7.32e+04 | ?? | 756 | \(--\) |
| K-100 | 3.16e-126 | 0.12 | 2.0 | 7.79e+04 | ?? | 1066 | \(--\) |
| E-3-4 | 1.35e-01 | 0.07 | 1.0 | 1.22e+04 | 0.02 | 35.4 | \(--\) |
| E-3-5 | 8.98e-03 | 0.05 | 1.0 | 1.49e+04 | 0.00 | 358 | \(--\) |
| E-5-3 | 1.04e-01 | 0.08 | 1.0 | 1.63e+04 | 0.00 | 403 | \(--\) |
| E-4-4 | 1.43e-02 | 0.07 | 1.0 | 1.74e+04 | 0.02 | 321 | \(--\) |
| ccp-5 | 2.31 | 0.07 | 1.0 | 9.24e+03 | 0.00 | 3.1 | 0.02 |
| ccp-6 | 1.40 | 0.06 | 1.0 | 1.28e+04 | 0.04 | 8.8 | 44.6 |
| ccp-7 | 5.21e-01 | 0.06 | 1.0 | 1.67e+04 | ?? | 29.5 | \(--\) |
| cc-8-5 | 1.77e-03 | 0.04 | 1.0 | 8.79e+03 | 0.02 | 6.9 | 0.03 |
| cc-8-7 | 2.75e+01 | 0.06 | 1.0 | 7.60e+03 | 0.00 | 7.2 | 0.07 |
| cc-8-9 | 1.36e+04 | 0.05 | 1.0 | 7.38e+03 | 0.02 | 9.3 | 0.5 |
| cc-8-11 | 1.40e+06 | 0.07 | 1.0 | 7.22e+03 | 0.01 | 14.1 | 2.4 |
| rvs-10-20 | 1.51e-05 | 0.07 | 1 | 1.11e+04 | 0.06 | 3.2 | 0.01 |
| rvc-10-30 | 5.30e-02 | 0.05 | 1 | 9.65e+03 | 0.05 | 4.5 | 0.09 |
| rvs-10-40 | 3.85e-04 | 0.09 | 1 | 7.19e+03 | 0.05 | 4.8 | 0.4 |
| rvc-10-50 | 0.440 | 0.05 | 1 | 7.27e+03 | 0.02 | 3.7 | 0.9 |
| rvs-15-45 | 1.01e-03 | 0.07 | 1 | 9.34e+03 | 0.04 | 6.2 | 87 |
| rvc-15-60 | 8.72e-04 | 0.05 | 1 | 8.98e+03 | 0.02 | 9.1 | 549 |
| rvs-20-40 | 2.86e-10 | 0.06 | 1 | 1.23e+04 | ?? | 7.7 | \(--\) |
| rvc-50-200 | 1.04e-18 | 0.08 | 1 | 2.25e+03 | ?? | 215 | \(--\) |
| rvc-60-90 | 1.48e-39 | 0.10 | 3.1 | 5.18e+04 | ?? | 267 | \(--\) |
| rvc-80-120 | 5.94e-57 | 0.08 | 4.5 | 7.88e+04 | ?? | 763 | \(--\) |
| rvc-100-150 | 1.31e-75 | 0.07 | 5.8 | 1.00e+05 | ?? | 2237 | \(--\) |
| rvs-60-120 | 2.04e-72 | 0.11 | 2.2 | 5.39e+04 | ?? | 422 | \(--\) |
| rvs-80-160 | 5.17e-106 | 0.19 | 3.0 | 7.50e+04 | ?? | 1069 | \(--\) |
| rvs-100-200 | 9.13e-142 | 0.13 | 4.1 | 1.05e+05 | ?? | 2520 | \(--\) |
| rvc-60-120 | 2.72e-33 | 0.14 | 2.1 | 5.25e+04 | ?? | 354 | \(--\) |
| rvc-80-160 | 5.52e-49 | 0.15 | 3.0 | 7.93e+04 | ?? | 1052 | \(--\) |
| rvc-100-200 | 1.04e-65 | 0.16 | 4.0 | 1.04e+05 | ?? | 2278 | \(--\) |
Table A.4. Volume Estimation for V-polytopes by volume
For each polytope volume performs 10 runs; volume: average estimated volume; k: average number of bodies in MMC; Refl: average number of reflections (or boundary oracle calls) performed by billiard_walk; error: average relative error; time: average time in seconds; exact time: time in seconds for exact volume computation by qhull (using R package geometry); \(--\) implies the execution failed due to memory issues or exceeded 1hr; ??: the exact volume is unknown; \(\epsilon = 0.1\) in all cases.
Table A.5.
| no rounding (C is a ball) | rounding | H-polytope approx. | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Z-polytope | k | Refl | time | k | Refl | time | k | Refl | time |
| \(Z_{\mathcal {U}}\)-30-60 | 3 | 3.09e+04 | 54 | 1 | 2.80e+04 | 50 | 1 | 1.57e+04 | 35.9 |
| \(Z_{\mathcal {U}}\)-40-80 | 4 | 4.23e+04 | 126 | 1 | 3.89e+04 | 115 | 1 | 1.71e+04 | 67.1 |
| \(Z_{\mathcal {U}}\)-50-100 | 5 | 5.53e+04 | 282 | 1 | 5.10e+04 | 270 | 1 | 1.84e+04 | 133 |
| \(Z_{\mathcal {U}}\)-60-120 | 7 | 9.04e+04 | 825 | 1 | 6.61e+04 | 575 | 2 | 3.52e+04 | 369 |
| \(Z_{\mathcal {U}}\)-30-150 | 1 | 7.01e+03 | 57 | 1 | 1.61e+04 | 111 | 1 | 7.24e+03 | 64 |
| \(Z_{\mathcal {U}}\)-40-200 | 1 | 7.79e+03 | 126 | 1 | 2.07e+04 | 323 | 1 | 7.81e+03 | 163 |
| \(Z_{\mathcal {U}}\)-50-250 | 1 | 8.33e+03 | 319 | 1 | 2.67e+04 | 858 | 2 | 1.24e+04 | 414 |
| \(Z_{\mathcal {U}}\)-60-300 | 1 | 9.53e+03 | 721 | 1 | 3.35e+04 | 2121 | 2 | 1.37e+04 | 1168 |
Table A.5. Comparisons between Rounding and the H-polytope Approximation used as Body in MMC of Section 3
For each Z-polytope volume performs 10 runs. k: average number of bodies in MMC; Refl: average number of reflections (or boundary oracle calls) performed by billiard_walk; time: average runtime of volume in seconds. We set \(\epsilon = 0.1\) in all cases. Bold marks best runtimes. The volumes for each Z-polytope agree up to at most \(\epsilon = 0.05\) and are thus omitted.
Figure A.1.
B Software Tutorial
This section aims in providing a brief outline on how to reproduce the computational results of this article. First, to run the C++ implementation of volume we suggest to use the R interface. To install the R interface follow the instructions given at the repository:
After installation the user should open an R terminal or alternatively use Rstudio and load the volesti package running:
library(volesti)
Then we start using the package by illustrating polytope generators that could be used to generate polytopes in various representations.
Note that when the value of the input flag of generator in the gen_rand_zonotope() is ’uniform’ then the function generates a random \(Z_{\mathcal {U}}\)-d-n. Other options are ’gaussian’ for random \(Z_{\mathcal {N}}\)-d-n and ’exponential’ for random \(Z_{Exp}\)-d-n.
To estimate the volume of a polytope with our implementation of volume the user should use the function volume(). The output contains both the estimated volume and its logarithm (useful in overflow or underflow cases).
To request polytope rounding before volume estimation using the methods we introduced in Section 3 run:
To estimate the volume with algorithm CoolingGaussian as described in [17],
To compute the PCA over-approximation of a Z-polytope,
References
[1]
Radosław Adamczak, Alexander Litvak, Alain Pajor, and Nicole Tomczak-Jaegermann. 2010. Quantitative estimates of the convergence of the empirical covariance matrix in log-concave ensembles. J. of the American Mathematical Society 23, 2 (2010), 535–561.
[2]
M. Althoff and J. M. Dolan. 2014. Online verification of automated road vehicles using reachability analysis. IEEE Trans. Robotics 30, 4 (2014), 903–918.
[3]
Shiri Artstein-Avidan, Haim Kaplan, and Micha Sharir. 2020. On Radial Isotropic Position: Theory and Algorithms. (2020). arxiv:cs.CG/2005.04918.
[4]
M. Berkelaar, K. Eikland, and P. Notebaert. 2004. lp_solve 5.5, Open Source (Mixed-Integer) Linear Programming System. http://lpsolve.sourceforge.net/5.5/.
[5]
Dimitris Bertsimas and Santosh Vempala. 2004. Solving convex programs by random walks. J. ACM 51, 4 (July2004), 540–556.
[6]
Stephen Boyd and Lieven Vandenberghe. 2004. Convex Optimization. Cambridge University Press, New York, NY, USA.
[7]
A. Brieden, P. Gritzmann, R. Kannan, V. Klee, L. Lovász, and M. Simonovits. 1998. Approximation of diameters: Randomization doesn’t help. In 39th Symp. Foundations of Computer Science (FOCS). 244–251.
[8]
Andrew Brock, Jeff Donahue, and Karen Simonyan. 2019. Large scale GAN training for high fidelity natural image synthesis. In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019. OpenReview.net.
[9]
B. Büeler and A. Enge. 2000. VINCI. (2000). http://www.math.u-bordeaux1.fr/aenge/index.php?category=software&page=vinci.
[10]
L. Calès, A. Chalkis, I. Z. Emiris, and V. Fisikopoulos. 2018. Practical volume computation of structured convex bodies, and an application to modeling portfolio dependencies and financial crises. In 34th International Symposium on Computational Geometry (SoCG 2018) (LIPIcs), Bettina Speckmann and Csaba D. Tóth (Eds.), Vol. 99. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Dagstuhl, Germany, 19:1–19:15.
[11]
E. Rodney Canfield and Brendan D. McKay. 2007. The asymptotic volume of the Birkhoff polytope. (2007). arxiv:math.CO/0705.2422.
[12]
Michal Cerný. 2012. Goffin’s algorithm for zonotopes. Kybernetika 48 (2012), 890–906.
[13]
Apostolos Chalkis and Vissarion Fisikopoulos. 2021. Volesti: Volume approximation and sampling for convex polytopes in R. The R Journal 13, 2 (2021), 642–660.
[14]
Apostolos Chalkis, Vissarion Fisikopoulos, Elias Tsigaridas, and Haris Zafeiropoulos. 2021. Geometric algorithms for sampling the flux space of metabolic networks. In International Symposium on Computational Geometry, SoCG (LIPIcs). Schloss Dagstuhl - Leibniz-Zentrum für Informatik. https://arxiv.org/pdf/2012.05503. to appear.
[15]
Yuansi Chen. 2021. An almost constant lower bound of the isoperimetric coefficient in the KLS conjecture. Geom. Funct. Anal. 31 (2021), 34–61.
[16]
B. Cousins and S. Vempala. 2015. Bypassing KLS: Gaussian cooling and an \({O}^*(n^3)\) volume algorithm. In Proc. ACM STOC. 539–548.
[17]
B. Cousins and S. Vempala. 2016. A practical volume algorithm. Mathematical Programming Computation 8 (2016). Issue 2.
[18]
H. Cramer. 1946. Mathematical Methods of Statistics. Princeton University Press.
[19]
F. Dabbene, P. S. Shcherbakov, and B. T. Polyak. 2010. A randomized cutting plane method with probabilistic geometric convergence. SIAM Journal on Optimization 20, 6 (2010), 3185–3207. arXiv:
[20]
L. Devroye. 1984. Random variate generation for unimodal and monotone densities. Computing 32, 1 (April1984), 43–68.
[21]
M. Dyer and A. Frieze. 1988. On the complexity of computing the volume of a polyhedron. SIAM J. Comput. 17, 5 (1988), 967–974.
[22]
M. Dyer, A. Frieze, and R. Kannan. 1991. A random polynomial-time algorithm for approximating the volume of convex bodies. J. ACM 38, 1 (1991), 1–17.
[23]
G. Elekes. 1986. A geometric inequality and the complexity of computing volume. Discr. Comput. Geom. 1, 4 (1986), 289–292.
[24]
I. Z. Emiris and V. Fisikopoulos. 2018. Practical polytope volume approximation. ACM Trans. Math. Soft. 44, 4 (2018), 38:1–38:21. Prelim. version: Proc. SoCG 2014.
[25]
Dani Gamerman and Hedibert F. Lopes. 2006. Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference. CRC Press.
[26]
C. Ge and F. Ma. 2015. A fast and practical method to estimate volumes of convex polytopes. In Frontiers in Algorithmics, J. Wang and C. Yap (Eds.). Springer, 52–65.
[27]
Charles J. Geyer. 1991. Markov Chain Monte Carlo maximum likelihood. Interface Foundation of North America. http://conservancy.umn.edu/handle/11299/58440. Accepted: 2010-02-24T20:38:06Z.
[28]
E. Gover and N. Krikorian. 2010. Determinants and the volumes of parallelotopes and zonotopes. Linear Algebra Appl. 413 (2010), 28–40.
[29]
Gaël Guennebaud, Benoît Jacob, et al. 2010. Eigen v3. http://eigen.tuxfamily.org.
[30]
S. Haraldsdóttir, B. Cousins, I. Thiele, R. M. T. Fleming, and S. Vempala. 2017. CHRR: Coordinate Hit-and-Run with rounding for uniform sampling of constraint-based models. Bioinformatics 33, 11 (2017), 1741–1743. arXiv:http://oup.prod.sis.lan/bioinformatics/article-pdf/33/11/1741/25153723/btx052.pdf.
[31]
Runxin He and Humberto Gonzalez. 2017. Numerical synthesis of Pontryagin optimal control minimizers using sampling-based methods. In 2017 IEEE 56th Annual Conference on Decision and Control (CDC). IEEE, 733–738.
[32]
Vu Anh Huynh, Sertac Karaman, and Emilio Frazzoli. 2012. An incremental sampling-based algorithm for stochastic optimal control. In 2012 IEEE International Conference on Robotics and Automation. IEEE, 2865–2872.
[33]
S. Iyengar. 1988. Evaluation of normal probabilities of symmetric regions. SIAM J. Sci. Statist. Comput. 9, 3 (1988), 418–423.
[34]
H. Jia, A. Laddha, Y. T. Lee, and S. S. Vempala. 2020. Reducing Isotropy and Volume to KLS: An \(O(n^3\psi ^2)\) Volume Algorithm. (2020). arxiv:cs.DS/2008.02146
[35]
Fritz John. 2014. Extremum problems with inequalities as subsidiary conditions. In Traces and Emergence of Nonlinear Programming, Giorgio Giorgi and Tinne Hoff Kjeldsen (Eds.). Springer, Basel, 197–215.
[36]
Ian T. Jolliffe and Jorge Cadima. 2016. Principal component analysis: A review and recent developments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374, 2065 (2016), 20150202. arXiv:
[37]
Adam Tauman Kalai and Santosh Vempala. 2006. Simulated annealing for convex optimization. Mathematics of Operations Research 31, 2 (2006), 253–266. https://www.jstor.org/stable/25151723. Publisher: INFORMS.
[38]
D. E. Kaufman and R. L. Smith. 1998. Direction choice for accelerated convergence in Hit-and-Run sampling. Operations Research 46, 1 (1998), 84–95. https://www.jstor.org/stable/223065.
[39]
M. Kerber, R. Tichy, and M. F. Weitzer. 2017. Constrained triangulations, volumes of polytopes, and unit equations. In 33rd Intern. Symp. Computational Geometry (SoCG 2017). Schloss Dagstuhl - Leibniz-Zentrum für Informatik, Germany, 46:1–46:15.
[40]
Yunbum Kook, Yin Tat Lee, Ruoqi Shen, and Santosh S. Vempala. 2022. Sampling with Riemannian Hamiltonian Monte Carlo in a constrained space. CoRR abs/2202.01908 (2022). arXiv:2202.01908. https://arxiv.org/abs/2202.01908.
[41]
A. K. Kopetzki, B. Schürmann, and M. Althoff. 2017. Methods for order reduction of zonotopes. In Proc. IEEE Annual Conf. Decision & Control (CDC). 5626–5633.
[42]
D. Kraft. 1988. A software package for sequential quadratic programming. Technical Report, DFVLR-FB 88-28, Institut für Dynamik der Flugsysteme (1988).
[43]
W. Kühn. 1998. Rigorously computed orbits of dynamical systems without the wrapping effect. Computing 61 (1998), 47–67.
[44]
Aditi Laddha and Santosh Vempala. 2020. Convergence of Gibbs Sampling: Coordinate Hit-and-Run Mixes Fast. (2020). arxiv:cs.DS/2009.11338.
[45]
Y. T. Lee and S. S. Vempala. 2018. 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). Association for Computing Machinery, New York, NY, USA, 1115–1121.
[46]
L. Lovász, R. Kannan, and M. Simonovits. 1997. Random walks and an \({O}^*(n^5)\) volume algorithm for convex bodies. Random Structures and Algorithms 11 (1997), 1–50.
[47]
L. Lovász and S. Vempala. 2006. Hit-and-Run from a corner. SIAM J. Comp. 35, 4 (2006), 985–1005.
[48]
L. Lovász and S. Vempala. 2006. Simulated annealing in convex bodies and an O\(^*(n^4)\) volume algorithms. J. Computer & System Sciences 72 (2006), 392–417.
[49]
O. Mangoubi and N. K. Vishnoi. 2019. Faster polytope rounding, sampling, and volume computation via a sub-linear ball walk. In 2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS). 1338–1357.
[50]
L. Martino, V. Elvira, D. Luengo, J. Corander, and F. Louzada. 2016. Orthogonal parallel MCMC methods for sampling and optimization. Digital Signal Processing 58 (Nov.2016), 64–84.
[51]
Luca Martino, David Luengo, and Joaquín Míguez. 2018. Independent sampling for multivariate densities. In Independent Random Sampling Methods, Luca Martino, David Luengo, and Joaquín Míguez (Eds.). Springer International Publishing, Cham, 197–247.
[52]
K. G. Murty. 2009. Ball Centers of Special Polytopes. Technical Report. Department of Industrial & Operations Engineering, University of Michigan.
[53]
Hariharan Narayanan and Piyush Srivastava. 2020. On the mixing time of coordinate Hit-and-Run. (2020). arxiv:cs.DS/2009.14004.
[54]
Radford M. Neal et al. 2011. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo 2, 11 (2011), 2.
[55]
A. Pereira and M. Althoff. 2015. Safety control of robots under computed torque control using reachable sets. In Proc. IEEE Inter. Conf. Robotics Automation (ICRA). 331–338.
[56]
B. T. Polyak and E. N. Gryazina. 2014. Billiard walk - a new sampling algorithm for control and optimization. IFAC Proceedings Volumes 47, 3 (2014), 6123–6128. 19th IFAC World Congress.
[57]
Chaitra Sarathy, Martina Kutmon, Michael Lenz, Michiel E. Adriaens, Chris T. Evelo, and Ilja C. W. Arts. 2020. EFMviz: A COBRA toolbox extension to visualize elementary flux modes in genome-scale metabolic models. Metabolites 10, 2 (2020).
[58]
J. Schellenberger and B. O. Palsson. 2009. Use of randomized sampling for analysis of metabolic networks. J. Biological Chemistry 284, 9 (2009), 5457–5461.
[59]
M. Schumer and K. Steiglitz. 1968. Adaptive step size random search. IEEE Trans. Automat. Control 13, 3 (June1968), 270–276. Conference Name: IEEE Transactions on Automatic Control.
[60]
T. Talvitie, K. Kangas, T. Niinimäki, and M. Koivisto. 2018. Counting linear extensions in practice: MCMC versus exponential Monte Carlo. In AAAI Conference on Artificial Intelligence. https://aaai.org/ocs/index.php/AAAI/AAAI18/paper/view/16957.
[61]
Michael J. Todd and E. Alper Yildirim. 2007. On Khachiyan’s algorithm for the computation of minimum-volume enclosing ellipsoids. Discrete Applied Mathematics 155, 13 (2007), 1731–1744.
[62]
A. Venzke, D. K. Molzahn, and S. Chatzivasileiadis. 2021. Efficient creation of datasets for data-driven power system applications. Electric Power Systems Research 190 (2021), 106614.
[63]
Frederick Wong, Christopher K. Carter, and Robert Kohn. 2003. Efficient estimation of covariance selection models. Biometrika 90, 4 (2003), 809–830. http://www.jstor.org/stable/30042090.
Index Terms
- A Practical Algorithm for Volume Estimation based on Billiard Trajectories and Simulated Annealing
Recommendations
Practical Volume Estimation of Zonotopes by a New Annealing Schedule for Cooling Convex Bodies
Mathematical Software – ICMS 2020AbstractWe study the problem of estimating the volume of convex polytopes, focusing on zonotopes. Although a lot of effort is devoted to practical algorithms for polytopes given as an intersection of halfspaces, there is no such method for zonotopes. Our ...
Practical Polytope Volume Approximation
We experimentally study the fundamental problem of computing the volume of a convex polytope given as an intersection of linear halfspaces. We implement and evaluate randomized polynomial-time algorithms for accurately approximating the polytope’s ...
Efficient Random-Walk Methods for Approximating Polytope Volume
SOCG'14: Proceedings of the thirtieth annual symposium on Computational geometryWe experimentally study the fundamental problem of computing the volume of a convex polytope given as an intersection of linear inequalities. We implement and evaluate practical randomized algorithms for accurately approximating the polytope's volume in ...
Comments
Information & Contributors
Information
Published In
Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s).
Publisher
Association for Computing Machinery
New York, NY, United States
Publication History
Published: 11 May 2023
Online AM: 03 March 2023
Accepted: 01 February 2023
Revised: 15 December 2022
Received: 08 April 2022
Published in JEA Volume 28
Author Tags
Qualifiers
- Research-article
Contributors
Other Metrics
Bibliometrics & Citations
Bibliometrics
Article Metrics
- 0Total Citations
- 1,249Total Downloads
- Downloads (Last 12 months)811
- Downloads (Last 6 weeks)93
Reflects downloads up to 10 Nov 2024
Other Metrics
Citations
View Options
View options
View or Download as a PDF file.
PDFeReader
View online with eReader.
eReaderHTML Format
View this article in HTML Format.
HTML FormatGet Access
Login options
Check if you have access through your login credentials or your institution to get full access on this article.
Sign in