{'fraction_non_alphanumeric': 0.033336820448443, 'fraction_numerical': 0.011472608710813544, 'mean_word_length': 4.490714148956538, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient computation of the volume of a polytope in high-dimensions using Piecewise Deterministic Markov Processes

Computing the volume of a polytope in high dimensions is computationally challenging but has wide applications. Current state-of-the-art algorithms to compute such volumes rely on efficient sampling of a Gaussian distribution restricted to the polytope, using e.g. Hamiltonian Monte Carlo. We present a new sampling strategy that uses a Piecewise Deterministic Markov Process. Like Hamiltonian Monte Carlo, this new method involves simulating trajectories of a non-reversible process and inherits similar good mixing properties. However, importantly, the process can be simulated more easily due to its piecewise linear trajectories - and this leads to a reduction of the computational cost by a factor of the dimension of the space. Our experiments indicate that our method is numerically robust and is one order of magnitude faster (or better) than existing methods using Hamiltonian Monte Carlo. On a single core processor, we report computational time of a few minutes up to dimension 500.


Volume of polytopes
High dimensional integration and sampling is a pervasive challenge in modern science. A subproblem of prime importance consists of computing the volume of a polytope, a bounded convex set of R d defined as the intersection of a fixed set of half spaces (H-polytope) or the convex hull of a finite set of vertices (Vpolytope). This challenge, of estimating the volume of a set, is closely related to problems in Statistics of computing marginal likelihoods [1] or Bayes factors for comparing competing models [2,3]. It is also related to problems in Physics of estimating partition functions and free energies [4].
Other applications appear across a range of disciplines. For example, in systems biology, the study of genome wide metabolic models based on systems of ODEs requires studying polytopes defined as the intersection between a hyperplane and the null space of a stoichiometry matrix [5,6]. In robotics, the computation of reachable sets for time-varying linear systems is based upon special polytopes called zonotopes [7]. In finance, the cross sectional score of a portfolio is defined from the intersection between a simplex (representing assets) and hyperplanes or ellipsoids [8]. In artificial intelligence, ReLU networks can be characterized by the conjunction of a set of linear inequalities which define a polytope in the input domain known as the activation condition [9]. In a related vein, the regularization of neural networks using piecewise affine functions involves the evaluation of the volume of polytopes defined by intersections of these hyperplanes [10].
Complexity-wise, computing the volume is #-P hard irrespective of the representation of the polytope (H-polytope or V-polytope) [11]. Intuitively, any deterministic algorithm using a polynomial number of points and computing the corresponding convex hull omits an exponentially large fraction of the volume. This observation prompted the development of approximation algorithms delivering (ε, δ) approximations, that is volume approximations within a factor 1+ε with a probability at least 1−δ [12,13]. The complexity of such algorithms is measured by the number of calls to an oracle stating whether a point is inside the polytope, or alternatively computing the intersection between a line and the polytope boundary. Remarkably, over the years, the complexity has been lowered from O (d 23 ) [14] to O (d 5 ) [15] using a sequence of balls intersecting the convex body, then from O (d 4 ) [16] to O (d 3 ) [17] using a sequence of smooth probability densities restricted to the convex body. The reader is referred to [18] for the full history.
Importantly, it should be stressed that the backbone of such algorithms is a telescoping product (Section 2) whose individual terms are the ratio of the integrals (i) of Gaussians on the polytope, or (ii) of the identity function integrated on the intersection between the polytope and a convex body. The former strategy is theoretically faster but the latter puts less burden on the samplers-see next section, since uniform distributions are used.
While the previous works are remarkable from the theoretical standpoint, it is difficult to turn these algorithms into effective implementations. This task is indeed challenging due to possibly huge constants in the complexities and un-realistic worst cases. This state of affairs recently motivated the development of strategies relaxing the theoretical guarantees, based on novel algorithmic and statistical techniques [19,20,21,22]. While these methods are not provably correct in general, their performances in terms of accuracy and running time have proven satisfactory.

Samplers
The most recent volume computation algorithms mentioned above rely on a procedure that samples a Gaussian distribution restricted to the polytope. The standard way to sample a probability distribution π is to build a Markov chain with invariant distribution π.
For sampling from a Gaussian distribution restricted to a compact region, a number of Markov chains have been proposed, including the Ball Walk sampler [23], the Hit and Run and coordinate Hit and Run samplers [24,25,26,5], and Hamiltonian Monte Carlo [27,21].
The Ball Walk is a simple Markov chain where the next point is proposed in a ball centered around the current point. If the proposed point is outside of the polytope, or if it rejected by Metropolis-Hasting, the chain stays at the current point.
By comparison, Hit and Run is a rejection free algorithm: a line passing through the current point is chosen at random, and the next point is sampled from the intersection of this line and the inside of the polytope. While this Hit and Run mixes faster than the Ball Walk sampler, it can get stuck in corners [26,21]. The computational cost of Hit and Run can be reduced by restricting the lines chosen to axis-coordinate Hit and Run.
Hamiltonian Monte Carlo can also be viewed as a rejection free method. It produces trajectories that are based on Hamiltonian dynamics for a physical system whose potential energy is defined by the Gaussian distribution within the polytope, and reflects the trajectory if it hits the boundary of the polytope. This sampler does not get stuck in corners like Hit and Run and has a very good mixing. However, unlike Hit and Run, the computation of the intersection of the trajectories and the boundary requires inverse trigonometric functions, and the computational cost cannot be reduced like Coordinate Hit and Run.
Another sampler is based on the billiard walk [28], which can be seen as a special case of HMC for uniform distributions, and is amenable to a computational complexity reduction similar to the one of coordinate Hit and Run [20]. This sampler can only target uniform distributions and not Gaussians.
Finally, a new class of samplers based on Piecewise Detereministic Markov Processes (PDMPs) has emerged in the computational statistics community [29]. Recall that a PDMP is a deterministic process inbetween random jump events, which occur at a certain rate. Hamiltonian Monte Carlo can actually be seen as a PDMP. However, Hamiltonian Monte Carlo relies entirely on the Hamiltonian flow to preserve the target measure, while most of the PDMP samplers developed to date have straight paths. The target measure is then preserved by changing direction at opportune random times. These methods are non reversible, which prevents the diffusive behavior of reversible chains, and have good properties in high dimension [30].

Contributions
In this work we introduce the use of a PDMP, the Bouncy Particle Sampler [31], to sample from a Gaussian restricted to a polytope within algorithms for calculating the volume of the polytope. Like Hamiltonian Monte Carlo, the Bouncy Particle Sampler has good mixing properties. Moreover, we show how we can re-use many calculations so that it has a much lower computational overhead. The key idea is that the main computation in sampling from a polytope is checking whether and when a trajectory hits the boundary of the polytope. For the Bouncy Particle Sampler, the trajectories are straight lines. We can use properties of how the trajectories change, for example after reflecting off a boundary, to reduce the cost of recalculating when the trajectory will next hit a given face of the boundary of a d-dimensional polytope from O(d) to O(1).
We provide detailed experiments up to dimension 500, and a comparison with Hamiltonian Monte Carlo up to dimension 100. These show that our method can be an order of magnitude faster, for the same level of precision, as the current best algorithms using Hamiltonian Monte Carlo.

Volume estimation algorithms
Without loss of generality, we consider a polytope H in d dimensions with the origin 0 strictly inside H. We define the polytope H through a d × d matrix A and a d-vector b: We wish to estimate the volume of H, Our approach is based on the algorithm of [19]. Consider a sequence of isotropic Gaussians with marginal variances σ 2 0 < σ 2 1 < ... < σ 2 m . Let f i denote the density of the Gaussian with variance σ i . We discuss the choice of the variances below, but [19] assumes that σ m is sufficiently large so that the Gaussian of variance σ m is nearly flat on H. These assumptions mean that Thus the volume of H can be rewritten as: We can rewrite each ratio in the product as: which is just the expectation of fi(X) fi−1(x) where X is distributed as a Gaussian with variance σ 2 i−1 restricted to H. Finally, [19] assumes σ 0 is sufficiently small that almost all the mass of the Gaussian with variance σ 2 0 lies within in H. This means that H f 0 (x)dx ≈ 1. This leads to the following approach to estimate V ol(H): (1) Choose m, the variances σ 2 0 < · · · < σ 0 m and Monte Carlo sizes N 1 , . . . , N m .
(2) For each i = 1, . . . , m, use an MCMC algorithm, or other, to get N i draws x Ni from a Gaussian with variance σ 2 i−1 restricted to H.
(3) For each i = 1, . . . , m construct the estimator of (1) aŝ . (2) In the original algorithm [32], the number of phases, the variances σ i and the number of steps per phase were chosen deterministically, leading to a O (d 3 ) complexity. In the practical implementation [19], a heuristic was added to find a better sequence of Gaussians, reducing the number of steps. In addition, to reduce the number of steps N i , a convergence diagnosis based on a heuristic was added. It uses a sliding window of size W (tied to the mixing time of the random walk) to estimate the ratio of Eq. 2. The number of samples N i is therefore not a fixed number, and tuning W is a non trivial issue [21]. Currently the best implementations of this approach [21], uses Hamiltonian Monte Carlo [33] to sample from the restricted Gaussians in step (2).
In this work, we improve this algorithm in three respects. First, in using the first Gaussian, σ 0 is chosen so that the Gaussian with variance σ 2 0 has probability mass of between c min and c max within H. Practically, we use c min = 0.1 and c max = 0.2. This introduces an additional term to the estimator -we need to add an estimate of the log of the probability mass that f 0 places with H to (3). To estimate this probability mass, we sample points from the Gaussian and use rejection sampling. The advantage of this adaptation is that it allows us to take a much larger variance, σ 2 0 , for the initial Gaussian, and thus a smaller value of m.
Second, we remove the stopping criterion on the window size used as a proxy for the number of samples N i . The size of the sliding windows W has to be tuned, which proved complicated to do. Furthermore, the target precision was not reliably obtained. We use a different strategy where a global given computational budget N is used (outside of the tuning phase of the random walks). For each Gaussian of variance σ i , we compute ess i the Effective Sample Size (ESS) per iteration of the random walk using the method described in [34,Section 16.4.2]. Then we choose the target number of samples (N i ) i≤m such that each ratio has about the same number of independent samples, i.e. N i × ess i ≈ N j × ess j for all i and j, and Third, as explained in the next section, we introduce a novel strategy to sample the Gaussians.

Bouncy Particle Sampler on Unbounded Space
We describe the PDMPs, and proceed with the Bouncy Particle Sampler on an unbounded space.

PDMP.
A PDMP, z t , is a continuous time Markov process defined by: 2. a jump rate λ(z), and 3. a jump kernel q(·|z).
The process z t follows the deterministic flow until a jump event happens. Jumps happen with probability If a jump event happens at time t, the process z t jumps using the jump kernel q: Bouncy Particle Sampler. For a given target density π on R d , the state space of the Bouncy Particle Sampler is extended by adding a velocity vector, yielding the state space: Thus the state of our PDMP is z = (x, v) and the target density becomes µ( is the density of a multivariate normal of mean 0 and covariance matrix I d . The Bouncy Particle Sampler is defined as follow: Note that the latter formula corresponds to a reflection of vector v with respect to the gradient of the potential. In practice, computing the jump times requires simulating a Poisson process of intensity λ(z t ). This is done by sampling u from an exponential law of rate 1, then solving for t the equation For a Gaussian target of the form π(x) ∝ exp(−a x 2 ), finding the event times amounts to solving: which gives a quadratic equation.

Bouncy Particle Sampler Restricted to a Polytope
The Bouncy Particle Sampler just introduced requires the target density to be continuous and almost everywhere differentiable. Therefore restriction of a Gaussian to a polytope requires adding jumps whenever the process reaches the boundary [35]. We will write q b (· mod z) the jump kernel at the boundary. Our target density is of the form for some constant a that depends on the variance of the Gaussian, and where 1 H is the indicator function of H. At a point x of the boundary ∂H, we write n(x) the outward normal, , v < 0} the set of in-going velocities. Then the process will target the correct distribution if the jump kernel at the boundary q b satisfies the following condition [35]: where x is on the boundary, n(x) the outward normal at x, and v ∈ V + an outgoing velocity. The simplest dynamics that satisfy (6) are to reflect the trajectory off the boundary, so the new velocity becomes with n the normal to the boundary that was hit. Finally, the Bouncy Particle Sampler is not always ergodic [31]. To ensure ergodicity, we add a refresh event with constant rate λ ref resh [31]. At refresh events, the velocity is resampled from its marginal invariant distribution, the multivariate normal distribution.
The full algorithm is described in Algorithm 1, and an example trajectory can be found in Figure 1.

Algorithm 1 Bouncy Particle Sampler (BPS).
while t ≤ t max do compute τ H the intersection time with H compute τ evt the next event time by solving (5).

Efficient Implementation
An important feature of the Bouncy Particle Sampler is that the velocity is constant between events which leads to deterministic trajectories that are straight-lines. This makes it possible to substantially reduce the computational cost of calculating the times at which the trajectory will hit a boundary. By pre-calculating suitable matrix vector products we can reduce the cost of calculating when the process next hits a boundary from O(dk), where k is the number of hyperplanes that defines H, to O(d). This efficient implementation is not possible, for Hamiltonian Monte Carlo samplers due to its non-constant velocity.
Starting from a point x with velocity v, the intersection time τ i for hyperplane i is the solution of Finding the intersections for all hyperplanes requires computing the products Ax and Av, which leads to a complexity of dk. However, in our cases, we keep track of the values of Ax and Av and update them without having to fully recompute the matrix vector product. To update the value of Ax for each event, we use By precomputing An for each hyperplane of the polytope, we can avoid the matrix multiplication. Case 2: PDMP jump event. The new velocity v is the reflection of v with respect to the gradient of the potential. In the case of Gaussians, the gradient is colinear with x. Hence v = v − 2 x,v x x, and we can write Since Ax is known, we can also avoid the matrix-vector product. Finally, note that for a refresh event, we have no choice but to recompute the matrix-vector product. In practice the proportion of events that are of this type is small. Remark 1. In [6], a similar strategy is used for reflection on the boundary, in a setting where there is no PDMP jump event.

Tuning parameters of the Bouncy Particle Sampler
The Bouncy Particle Sampler has two parameters that requires to be tuned. We describe here the automatic tuning strategy we implemented.
The Bouncy Particle Sampler is currently described as a continuous time process. To extract points from the trajectory z t , we add another Poisson rate λ out . Each time an event happens with respect to this Poisson process, the current point z t is passed to the volume computation algorithm. To tune λ out , we use the following heuristic: to get a sample independent from the starting point, we require d events (see Cases 1 and 2 above) to happen. Hence, we automatically tune λ out at runtime so that there is on average d events between points passed to the volume algorithm.
Furthermore, as noted before, we have a refresh rate λ ref resh , whose tuning is important for the mixing time. To tune this parameter, we use an optimization procedure described in more details in Section 5.2.

Code overview
Our code is written in C++, using the number types discussed below. ESS calculations were carried out using the Autoppl library https://github.com/JamesYang007/autoppl.

Multiprecision
Trajectories escaping the polytope. As noted in [21], Hamiltonian Monte Carlo trajectories can escape the polytope due to numerical issues -a phenomenon becoming prevalent in high dimension. The solution introduced in [21] is to recompute a trajectory leaving the convex by using interval arithmetic with increased precision.
Since the trajectory escaping the polytope can lead to a completely wrong estimate at best, and a crash of the program at worse, a similar strategy is used in this work. Assuming the T k are the output times of the trajectories that are passed to the volume algorithm (i.e. the sequence used by the volume algorithm is (x T k ) k∈N ), the strategy is as follow: It should also be noticed that precision cannot be increased indefinitely. As a fallback, if the required precision reaches a certain threshold, we go back to time T i and sample a new velocity. In the Experiments below, we report on these outcomes using two statistics: • #M : the number of times the else above is entered.
• #R: the number of times the precision limit is reached.
Number types. To handle precision refinements, we use the number type boost::multiprecision::mpfr_float. Numerical values for volumes. Multiprecision is also required to compute the volume, in two guises. On the one hand, the exact volume of the polytopes used in our tests cannot be represented using double precision, when increasing the dimension. On the other hand, high precision floating points must be used when computing the exponential of the sum of the log of the ratios m i=1 log(Î i ). To handle these difficulties, we use the number type boost::multiprecision::cpp_dec_float, which makes it possible to specify the number of decimal digits used.

Setup and statistics of interest
Polytopes. We study our algorithm for three polytopes where the exact volume is known: • The cube, −1 ≤ x i ≤ 1, for i = 1, . . . , d.
Targeting a given error. Assessing the complexity of our algorithm as a function of the dimension-for a given polytope, requires finding the number of samples for which a target error on the volume estimate is obtained. To this end, consider a minibatch of repeats (24 in our case to exploit brute force parallelism on our computer), each using N samples. We run a binary search on N to obtain the smallest value for which the median error of the minibatch lies in the interval 4% ± 1%. (Nb: to speed up calculations, the binary search is also exited if the number of samples between two consecutive runs varies less than 5%.) In this experiment, the dimensions d = 50, 70, 100, 140, 175, 250 are used for each polytope. For a given statistic of interest (the median running time or the final number of samples N ), we then perform a linear regression in log log scale, to assess the polynomial scaling of the statistic.
The output of each computation contains • the final estimation of the volume, • the run-time (in seconds), and • the number of times the precision had to be increased to stay in the polytope.
Using a fixed number of samples. To compare PDMP against state-of-the-art volume algorithms based on Hamiltonian Monte Carlo [21], we resort to experiments at a fixed number of samples. The volume computation is launched for dimension d = 100, 500, using sample sizes N = 10 5 , 10 6 , 10 7 -the latter for d = 500 only.
Computer used. Calculations were run on a desktop DELL Precision 7920 Tower (Intel Xeon Silver 4214 CPU at 2.20GHz, 64 Gb of RAM), under Linux Fedora core 34.

Tuning
Based on the results in [36], heuristically, a small refresh rate will lead to better mixing for individual coordinates of the samples, whereas a higher refresh rate will mean better mixing for functionals such as the norm of the position. In our case, bouncing on the boundary of the polytope can play a similar role to refreshing the velocity. Hence we expect the optimal refresh rate to depend on the variance of the sampled restricted Gaussian. Thus we compute the ESS for the projection on each coordinate, and take the minimum, which we call ess min . Further we compute ess norm , the ESS associated to the norm of our samples. If ess min < ess norm , we decrease the refresh rate λ ref resh , and if ess min > ess norm , we increase the refresh rate.
However, refreshing the velocity is expensive, since it requires recomputing the matrix-vector product Av, with a cost of O(d 2 ) instead of O(d) for other events. Thus the new refresh rate is only accepted if it leads to a higher ESS per second, with the overall ESS being min(ess min , ess norm ). If not, the tuning is stopped (Fig. 2).
The number of consecutive samples used to evaluate the ESS results from a tradeof between the quality of the tuning, which increases with number of samples, and the time spent tuning the walk with In our experiments, the number of consecutive samples is restricted by runs in dimension 500, which leads us to use 100 consecutive samples to evaluate the ESS. is advised.

Results
PDMP: complexity. The experiments at prescribed error rate show a clear polynomial time complexity (Fig. 3), with a scaling around O(d 3.5 ) for the three test polytopes (Table 1). In particular, we can see that the time complexity is close to the number of samples times d 2 . This is consistent with intuition, since each sample requires on average O(d) events, and each event (except for refresh events) requires O(d) computations.  PDMP versus HMC: complexity. Hamiltonian Monte Carlo based methods require the specification of a so-called window size used to determine convergence of the ratio of (1) [21]. In the sequel, we use two window sizes: • Low profile comparison: the window size W = 250 is used to compare Hamiltonian Monte Carlo against PDMP with 10 5 samples -short-hand: HMC1; • High profile comparison: the window size W = 250+d √ d is used to compare Hamiltonian Monte Carlo against PDMP with 10 6 samples -short-hand: HMC2.
Two facts stand out (Tables S3 and S4). On the one hand, for the three polytopes, the two algorithms compared yield equivalent error estimates. On the other hand, PDMP yields a speedup between one and two orders of magnitude over HMC. To be precise, in dimension d = 100 e.g.: between ∼ 56 and 135 for the  Table 1. low profile comparison, and between ∼ 14 and 44 for the high profile comparison. This comparison shows the superiority of linear trajectories over curved ones in HMC.
As a final analysis, we challenge PDMP with calculations up to dimension d = 500, using a fixed number of samples (Tables S5 and S6). Obtaining a relative error below 10% now requires using a number of samples N one order of magnitude larger. Still, with a running time of the order of hours, such large dimensions remain tractable.
PDMP versus HMC: numerics. Another important issue is that of numerical robustness faced by geometric algorithms in general [37]. Along with the volumes, we monitor the number of multi-precision refinements - Table S2. In sharp contrast with the method based on Hamiltonian Monte Carlo [21], such refinements are not triggered by the PDMP sampler -the statistic #M is null and therefore so is #R. We conjecture that multiprecision is required with HMC due to intersections between curved trajectories and the polytope boundary, requiring the (numerically challenging) computation of an arctan.
Cooling schedules: using Gaussians versus balls. As mentioned in the Introduction, an alternative algorithm uses balls intersected with the polytope instead of Gaussians restricted to the polytope [20]. Since this algorithm only requires uniform sampling, the simpler billiard walk sampler can be used [20]. It should be noted than our sampling strategies are essentially the same than the one provided in [20], the main difference being that our random walk handle Gaussians while still having a low computational complexity.
Ball cooling essentially amounts to changing the functions f i in Eq. (2) from Gaussian densities to indicator functions of balls. In that case, the ratio found in Eq. (2) is exactly 0 or 1. In principle, this should lead to a higher variance ofÎ i of Eq. (2) than for Gaussians, which would be less efficient. This is consistent with known theoretical results [15,17] which seems to indeed indicate a greater efficiency of the Gaussian cooling compared to balls.
Interestingly, the implementation provided in the preprint [20] yields better results, which contradicts theory [15,17]. This calls for further analysis in two directions. The first is the analysis of the cooling schedule based on a sequence of balls [20] versus a sequence of Gaussians [19]. The second one is the tuning of the random walk generating samples (in our case: refresh rate λ ref resh , output rate λ output ).
We believe that two things might be happening. First, the tuning strategy used for balls might outperform ours, which is expensive since it relies on calculating the ESS. Second, we rely on the sequence of Gaussians provided by [19], while a new one, possibly more optimal for balls, is used in [20].

Conclusion
Recent developments for polytope volume calculation algorithms have exploited two strategies, namely using a sequence of Gaussians or a sequence of balls intersecting the polytope. Our work presents improvements for the former, yielding a substantial speedup (more than one order of magnitude on running times) over state-of-the-art HMC based methods. On the other hand, recent work in the latter vein has resulted in faster and more accurate results [20], an unexpected observation given the theoretical bounds known to date for both classes of methods. An interesting avenue for future work will therefore consist of hybridising both approaches.
For convex bodies with piecewise C 1 boundaries, our algorithm could also be applied. While the improvement to the oracle complexity might not apply, the intersection of linear trajectories with the boundary would improve on HMC.
On the practical aspect, there are a few improvements that could be made in future work. The most obvious one would be to re-use the points sampled at previous phases using importance sampling. The second would be to introduce a way of refreshing velocities that do not require recomputing the matrixvector product.
Polytope volume calculations algorithms first underwent major improvements in the theoretical realm (with bounds on the mixing times but algorithms lacking efficiency in practice), and more recently in the practical realm (with efficient algorithms lacking theoretical bounds). Combining both is clearly an outstanding challenge ahead.