key: cord-0047231-uvpc245s
authors: Chalkis, Apostolos; Emiris, Ioannis Z.; Fisikopoulos, Vissarion
title: Practical Volume Estimation of Zonotopes by a New Annealing Schedule for Cooling Convex Bodies
date: 2020-06-06
journal: Mathematical Software - ICMS 2020
DOI: 10.1007/978-3-030-52200-1_21
sha: b102ceeb7ecf458b8cfead395fca11be875f791f
doc_id: 47231
cord_uid: uvpc245s

We 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 algorithm is based on Multiphase Monte Carlo (MMC) methods, and our main contributions include: (i) a new uniform sampler employing Billiard Walk for the first time in volume computation, (ii) a new simulated annealing generalizing existing MMC by making use of adaptive convex bodies which fit to the input, thus drastically reducing the number of phases. Extensive experiments on zonotopes show our algorithm requires sub-linear number of oracle calls in the dimension, while the best theoretical bound is cubic. Moreover, our algorithm can be easily generalized to any convex body. We offer an open-source, optimized C++ implementation, and analyze its performance. Our code tackles problems intractable so far, offering the first efficient algorithm for zonotopes which scales to high dimensions (e.g. one hundred dimensions in less than 1 h).

Volume computation is a fundamental problem with many applications. It is #P-hard for explicit polytopes [7, 11] , and APX-hard [9] for convex bodies in the oracle model. Therefore, a significant effort has been devoted to randomized approximation algorithms, starting with the celebrated result in [8] with complexity O * (d 23 ) oracle calls, where O * (·) suppresses polylog factors and dependence on error parameters, and d is the dimension. Improved algorithms reduced the exponent to 5 [13] and further results [5, 14] reduced the exponent to 3. Current theoretical results consider either the general oracle model or polytopes given as an intersection of halfspaces (i.e. H-polytopes). Regarding implementations, the approach of [13] led to the first practical implementation in [10] for high dimensions, followed by another practical implementation [6] based on [5, 14] . However, both implementations can handle only H-polytopes.

An important class of convex polytopes are zonotopes [15] . A zonotope is the Minkowski sum of k d-dimensional segments. Equivalently, given a matrix G ∈ R d×k a zonotope can be seen as the affine projection of the hypercube [−1, 1] k to R d using the matrix G, while the columns of G are the corresponding segments (or generators). Zonotopes are centrally symmetric and each of their faces are again zonotopes. We call the order of a zonotope P the ratio between the number of generators of P over the dimension. For a nice introduction to zonotopes we refer to [20] .

Volume approximation for zonotopes is of special interest in several applications in smart grids [1] , in autonomous driving [2] or human-robot collaboration [16] . The complexity of algorithms that work on zonotopes strongly depends on their order. Thus, to achieve efficient computations, a solution that is common in practice is to over-approximate P , as tight as possible, with a second zonotope P red of smaller order, while vol(P red ) is given by, an easy to compute, closed formula. A good measure for the quality of the approximation is the ratio of fitness, ρ = (vol(P red )/vol(P )) 1/d , which involves a volume computation problem [3] . Existing work (e.g. in [12] ) uses exact -deterministic volume computation [11] , and thus ρ can not be computed for d > 10 in certain applications.

A typical randomized algorithm uses a Multiphase Monte Carlo (MMC) technique, which reduces volume approximation of convex P to computing a telescoping product of ratios of integrals. Then each ratio is estimated by means of random walks sampling from a proper multivariate distribution. In this paper we rely on MMC of [13] which specifies a sequence of convex bodies P m ⊆ · · · ⊆ P 0 = P , assuming P is well-rounded, i.e. B d ⊆ P ⊆ C √ dB d , where C is constant and B d is the unit ball. We define a sequence of scaled copies of B d , and let P i = (2 (m−i)/d B d ) ∩ P, i = 0, . . . , m. One computes vol(P m ) and applies:

There is a closed-form expression to compute vol(P m ) = vol(B d ) . Each ratio r i = vol(P i+1 )/vol(P i ) in Eq. (1) can be estimated within arbitrary small error i by sampling uniformly distributed points in P i and accept/reject points in P i+1 so vol(P ) can be derived after m multiplications. The estimation of r i shows how sampling comes into the picture. In [10] , assuming rB d ⊆ P ⊆ RB d for r < R, they get m = d lg(R/r) . The issue is to minimize m while each ratio remains bounded by a constant, and to use a random walk that converges, after a minimum number of steps, to the uniform distribution. The first would permit a larger approximation error per ratio without compromising overall error, while it would require a smaller uniform sample to estimate each ratio. The second would reduce the cost per sample point. Total complexity is determined by the number of ratios, or phases, multiplied by the number of points, or steps, to estimate each ratio, multiplied by the cost to generate a point. The first two factors are determined by the MMC and the third by the random walk.

Previous Work. Exact volume computation for zonotopes can be reduced to a sum of absolute values of determinants, with an exponential number of summands in d [11] . Practical algorithms for volume computation of zonotopes are limited to low dimensions (typically ≤ 10 in [6] ). This is due to two main reasons: current algorithms create a long sequence of phases in MMC for zonotopes, and the boundary and membership oracles are costlier than for H-polytope, as they both reduce to Linear Programs (LP). In [6] , they consider low dimensional zonotopes: in R 10 with k = 20 generators, the algorithm performs 1.92 × 10 5 Boundary Oracle Calls (BOC), whereas our algorithm requires only 8.50 × 10 3 BOCs, and for d = 100, k = 200 it performs 6.51 × 10 4 BOCs (see Table 1 ).

In [19] they present an implementation of an efficient algorithm that computes Minkowski sums of polytopes (generalization of zonotopes). In [18] they propose a randomized algorithm for enumerating the vertices of a zonotope.

Our Contribution. We focus on zonotopes and introduce crucial algorithmic innovations to overcome the existing barriers, by reducing significantly the number of oracle calls. Thus, our method scales to high dimensions (d = 100 in ≤1 h), performing computations which were intractable till now. We use a new simulated annealing method in order to define a sequence of appropriate convex bodies, instead of balls, in MMC, and we exploit the fast convergence of Billiard Walk (BW) [17] to the uniform distribution. We experimentally analyze complexity by counting the number of BOCs, since BW uses boundary reflections.

The new simulated annealing specifies the P i 's by exploiting the statistical properties of the telescoping ratios to drastically reduce the number of phases. In particular, we bound each ratio r i = vol(P i+1 )/vol(P i ) to a given interval [r, r + δ] with high probability, for some real r. Moreover, our MMC generalizes balls, used in [13] and previous papers, by taking as input any convex body C and constructing the sequence by only scaling C. It does not need an enclosing body of P nor an inscribed ball (or body), unlike [10, 13] .

Most of the previous algorithms use a rounding step before volume computation, as preprocessing, to reduce the number of phases in MMC. However, rounding requires uniform sampling from P which makes it costly for zonotopes because of the expensive oracle calls. Our approach is to exploit the fact that the schedule uses any body C and skip rounding by letting C be an H-polytope that fits well to P . The idea is to construct C fast and reduce the number of phases and the total runtime more than a rounding preprocessing would do in practice.

We prove that the number of bodies defined in MMC is, with high probability, m = O(lg(vol(P )/vol(P m ))), where P m = qC ∩ P , for some q ∈ R, is the body with minimum volume, and vol(qC∩P ) vol(qC) ∈ [r, r + δ]. The bound on m is not surprising, as it does not improve worst-case complexity [5] , if C is a ball, but offers crucial advantages in practice. First, the hidden constant is small. More importantly, if C is a good fit to P , vol(P m ) increases and m decreases (Fig. 1) . We also show that, for constant d, and k (number of generators) increasing, m decreases to 1, when we use ball in MMC, since the schedule constructs an enclosing ball of P . Intuitively, while order increases for constant d, a random zonotope approximates the hypersphere. The latter can be approximated up to in the Hausdorff metric by a zonotope with k ≤ c(d) [4] . This does not directly prove our claim on m but strengthens it intuitively. So, in our experiments, the number of phases is m ≤ 3 for any order, without rounding for d ≤ 100.

Considering uniform sampling, BW defines a linear trajectory starting at the current point, using boundary reflections [17] . No theoretical mixing time exists. We show that with the right selection of parameters, BW behaves like an almost perfect uniform sampler even if the walk length is 1. In particular, for this walk length, it generates just O * (1) points per phase, with sub-linear number of reflections per point, and provides the desired accuracy. To stop sampling when estimating ratio r i we modify the binomial proportion confidence interval. We use the standard deviation of a sliding window of the last l ratios, thus defining a new empirical convergence criterion; l = O(1) suffices with BW.

Our software contributions build upon and enhance volesti 1 a C++ open source library for high dimensional sampling and volume computation with an R interface. We experimentally show that the total number of oracle calls grows as o * (d) for random zonotopes; the best available theoretical bound is O * (d 3 ) [5] .

The algorithm first constructs a sequence of convex bodies C 1 ⊇ · · · ⊇ C m intersecting the zonotope P ; the C i 's are determined by simulated annealing.

A typical choice of C i 's in this paper is co-centric balls, or centrally symmetric H-polytopes. C m is chosen for its volume to be computed faster than vol(P ) and easily sampled. Then,

where P 0 = P . Let r i = vol(Pi+1) vol(Pi) , i = 0, . . . , m − 1, r m = vol(Pm) vol(Cm) .

We use BW to sample approximately uniform points in P i at each phase i. BW picks a uniformly distributed line through the current point. It walks on a linear trajectory of length L = −τ ln η, η ∼ U(0, 1), reflecting at the boundary. BW can be used to sample only uniform points; in [17] they experimentally show that BW converges fast to the uniform distribution when τ ≈ diam(P ). The membership oracle is a feasibility problem. A point p ∈ P iff the following region is feasible: 

The second extreme point which corresponds to a negative value of λ is not used by BW. For the BW we need the normal of the facet that intersects to compute the reflection of the trajectory if needed. We keep the generators that corresponds to x i = −1, 1 and then the normal vector is computed straightforwardly.

Given P , the annealing schedule generates the sequence of convex bodies C 1 ⊇ · · · ⊇ C m defining P i = C i ∩ P and P 0 = P . The main goal is to restrict each ratio r i in the interval [r, r + δ] with high probability. We define the following two statistical tests, which can be reduced to t-tests:

The U-test and L-test are successful iff null hypothesis H 0 is rejected, namely r i is upper bounded by r+δ or lower bounded by r, with high probability, respectively. If we sample N uniform points from P i then r.v. X that counts points in P i+1 , follows X ∼ b (N, r i ) , the binomial distribution, and Y = X/N ∼ N (r i , r i (1 − r i )/N ). Then each sample proportion that counts successes in P i+1 over N is an unbiased estimator for the mean of Y , which is r i . Perform L-test and U-test Input : convex bodies P 1 , P 2 , cooling parameters r, δ, s.l. α, ν, N ∈ N Sample νN uniform points from P 1 Partition νN points to lists S 1 , . . . , S ν , each of length N Compute ratiosr i = |{p ∈ P 2 : p ∈ S i }|/ N, i = 1, . . . , ν Compute the mean,μ, and st.d., s, of the ν ratios

Let us now describe the annealing schedule: Each C i in C 1 ⊇ · · · ⊇ C m is a scalar multiple of a given body C. Since our algorithm does not use an inscribed body, initialization computes the body with minimum volume, denoted by C or C m . This is the last body in the sequence. The algorithm sets P 0 = P and employs C to decide stopping at the i-th phase.

Initialization. Given C, and interval [q min , q max ], one employs binary search to compute q ∈ [q min , q max ] s.t. both U-test(qC, qC ∩ P ) and L-test(qC, qC ∩ P ) are successful. Let q = (q min + q max )/2. If U-test(qC, qC ∩ P ) succeeds and L-test(qC, qC ∩P ) fails, we continue to the left-half of the interval. With inverse outcomes, we continue to the right-half of the interval. If both succeed, stop and set C = qC. The output is C , denoted by C m at termination.

At iteration i, the algorithm determines P i+1 s.t. volume ratio r i ∈ [r, r + δ] with high probability. The schedule samples νN points from P i and binary searches for a q i+1 in an updated interval [q min , q max ] s.t. both U-test(P i , q i+1 C ∩P ) and L-test(P i , q i+1 C ∩P ) are successful. Then set P i+1 = q i+1 C ∩ P .

Stopping and Termination. The algorithm uses C ∩ P in the i-th iteration for checking whether vol(C ∩ P )/vol(P i ) > r with high probability, using only L-test, and then stops if L-test(P i , C ∩ P ) holds. Then, set m = i + 1, and P m = C ∩ P .

In the t-tests, errors of different types may occur, thus, binary search may enter intervals that do not contain ratios in [r, r +δ]. Hence, there is a probability that annealing schedule fails to terminate. Let β capture the power of a t-test: Fig. 2 . Number of bodies in MMC. For each dimension we generate 10 random zonotopes and we compute the number of bodies, m, in MMC when C is ball. We keep the zonotope with the larger m and then, for that one, we compute m when C is the P-approx.

The annealing schedule allows as to use any C which must (a) be a good fit to P , (b) allow for more efficient sampling than in P , and (c) for faster volume calculation than of vol(P ). For low order ones C shall be an enclosing H-polytope that fits well to P . Indeed it is possible that with certain choices for C rounding is not needed. We define a centrally symmetric H-polytope with ≤ 2k facets:

Construct P-approx Input : The generator matrix G ∈ R d×k of zonotope P Output: An H-polytope C ⊃ P compute the eigenvectors of G T G (has k − d zero eigenvalues) let the eigenvectors of k − d zero eigenvalues of G T G form E ∈ R k×(k−d) . compute an orthonormal basis for E, and the orthogonal complement

We perform extended experiments analyzing practical complexity. We use eigen 2 for linear algebra and lpSolve 3 for LPs. All experiments were performed on a PC with Intel R Core TM i7-6700 3.40 GHz 8 CPU and 32 GB RAM. We use three zonotope generators. All of them pick uniformly a direction for each one of the k segments. Then, (a) Z U -d-k: the length of each segment is uniformly sampled from [0, 100], (b) Z N -d-k: the length of each segment is random from In particular, when we use P-approx, m is smaller for order ≤ 4 compared to using balls without rounding. For order equal to 5 the number of balls in MMC is smaller compared to the number of bodies when the choice is the Papprox. Notice than when we use balls in MMC, m decreases for constant d as k increases. Table 1 shows that, for high-order zonotopes, m = 1, which implies one or two rejection steps, while the run-time is smaller when we use ball in MMC. It also reports the difference in the run-time for random zonotopes of order = 2 between the cases of using ball and the P-approx in MMC. In all our experiments, BW performs only O * (1) steps per phase with just a factor of i hidden in the complexity. The plot that counts the BW reflections per point in Fig. 3 imply this number grows sub-linearly in d. Hence, the total number of BOCs grows sub-linearly in d. Table 1 . Volume estimation for zonotopes. For each Z-d-k we approximate its volume using ball and the P-approx in MMC. Body stands for the type of body in MMC; order = k/d, Vol the average of volumes over 10 runs; m the average number of bodies in MMC; OracleCalls is the average number of BOCs; time is average time in seconds. We set the error parameter = 0.1 in all cases.

Formal and compositional analysis of power systems using reachable sets

Online verification of automated road vehicles using reachability analysis

Reachability Analysis and Its Application to the Safety Assessment of Autonomous Cars

Approximating the ball by a Minkowski sum of segments with equal length

Bypassing KLS: Gaussian cooling and an O * (n 3 ) volume algorithm

A practical volume algorithm

On the complexity of computing the volume of a polyhedron

A random polynomial-time algorithm for approximating the volume of convex bodies

A geometric inequality and the complexity of computing volume

Practical polytope volume approximation

Determinants and the volumes of parallelotopes and zonotopes

Methods for order reduction of zonotopes

Random walks and an O * (n 5 ) volume algorithm for convex bodies

Simulated annealing in convex bodies and an O * (n 4 ) volume algorithms

On zonotopes

Safety control of robots under computed torque control using reachable sets

Billiard walk -a new sampling algorithm for control and optimization

A randomized algorithm for enumerating zonotope vertices

Implementation and parallelization of a reverse-search algorithm for Minkowski sums

Lectures on Polytopes