# Algorithm for finding the volume of a convex polytope

It's easy to find the area of a convex polygon by division into triangles, but what is the optimal way of finding the volume of higher-dimensional convex bodies? I tried a few methods for dividing them into simplices, but gave it up and went with a Monte Carlo estimation scheme instead. Bonus question: How to find the surface area of those same convex bodies?

EDIT: To answer David's question: the data set is a Voronoi tessellation of an n-dimensional volume (n usually 4) with a periodic boundary (like a torus). So I have the coordinates of the vertices of the convex bodies as well as the connectivity of all the facets, faces, etc. For the Monte Carlo I mentioned, I did convert everything to half-spaces, so I think that was not very difficult.

-

As a follow-up to Barton's response: for hardness results, see I. Bárány & Z. Füredi, Computing the volume is difficult, Discrete and Computational Geometry, 1987. But there are polynomial time approximation schemes for volume of convex bodies independent of dimension, based on random walks within the body: see e.g. M. Dyer, A. Frieze, & R. Kannan, A random polynomial-time algorithm for approximating the volume of convex bodies, J. ACM 1991 and R. Kannan, L. Lovász, & M. Simonovits, Random walks and an O*(n^5) volume algorithm for convex bodies, Random structures and algorithms, 1997.

Before getting to these answers, though, it's important to ask: how is your input represented? Is it a convex hull of a set of points, and all you know is the points? Is it an intersection of halfspaces? Is it given to you as an entire face lattice? Also, what range of dimensions do you care about The upper and lower bounds above are quite general but the dimension and input representation may still make a difference to the answer.

-

I think this problem is hard--the known algorithms are both slow and nontrivial to implement. See Exact Volume Computation for Polytopes for a survey. An interesting feature is that there are various algorithms which are well suited for different kinds of polytopes.

As a practical answer, Qhull can compute volumes and surface areas.

-

Matthias Beck and Dennis Pixton used almost 17 GigaHertz years to compute the value of the 10-dimensional Birkhoff polytope. Read about it, and find the answer, here.

If you can get the same result quicker, I'm sure they'd be delighted to know how you did it.

-

Here's a fairly straightforward solution for polyhedra (3 dimensions), with running time O(v+ve), where v is the number of vertices and e is the number of edges. I suppose it could be extended to higher dimensions, but it would probably have much worse running time (I fear roughly exponential as in O(vn), where n is the number of dimensions).

Let our polyhedron have n vertices, defined by their x,y,z coordinates: v1, v2, ..., vn and let the lowest point be v1 and the origin (modify the values for the others accordingly), and let it have e edges, defined by the vertices which they connect. Then, since we have coordinates for the vertices (as that is how we defined them), there must be a "ground level" plane p0 running through the x and z axes (the y-axis being height, and the ground never having an elevation). Then, let v2 be the point closest to the ground plane (shortest line perpendicular to the plane), and let v3 be the next closest, etc, through vn.

Through each of the points v2 through vn, draw a plane perpendicular to the ground, and let them be numbered pm, where m is the subscript of the vertex through which it was drawn. Then, the volume of our polyhedron is equal to the sum of the volumes of the figures between the planes. We should have something resembling this:

Let the heights between the segments be h1 through hn-1, where height hj is the height between planes pj and pj+1.

Now, through each plane, we have a polygon (or more, if the figure is concave), whose vertices' coordinates can be calculated easily as follows:
Let the edge that runs through the plane pj have endpoints va and vb. Then, the displacement vector is vb - va (assuming the coordinates of v are in vector-form), and the percentage travelled up is $\frac{h_j-h_a}{h_{b-1}-h_a}$. Multiply this by vb - va and add to va to calculate the new point of intersection for that edge:
Intersection point = $(v_b-v_a)\frac{h_j-h_a}{h_{b-1}-h_a}+v_a$
The area of these polygons can be determined using triangles, or a simplification of this very process in just 2 dimensions.

PlanetMath says that the volume of a prismatoid (which is the type of figure contained between sequential planes) is $h\frac{B_1 + B_2 + 4M}{6}$, where the Bs are the areas of the parallel polygons and M is the area of the midway polygon, which is exactly halfway between them (and parallel to them). Since we already know the area of each of the end polygons, and we can easily calculate the vertices of the midway polygon (using the previous paragraph's method), we can calculate the volume of the resulting prismatoids. Adding them up yields the total volume of the polyhedron.

I suppose that the only real issue in this case, then, is, via code, determining which edges run through any particular plane, but if we were to actually look at it, we could tell very easily.

A simpler version of this can be used to figure out the area of any polygon; simply draw lines through the vertices parallel to the x-axis and calculate the area of the resulting trapezoids as (b1+b2)/2

-

There is a randomized algorithm with O(n^4) running time based on the simulated annealing using the hit-and-run convex body sampling algorithm along with isotropy approximation by Vempala and Lovasz: here. The running time is the number of oracle calls for membership in the convex body and has an implicit polylogarithmic term.

-

As documented above the paper by Bueler, Enge and Fukuda is a good place to start and give many techniques for computing volumes of polytopes. In my experience using lrs to build the triangulation is a very good method. If your polytope has symmetries, then it is possible to use them to reduce the size of the computation. See http://www.liga.ens.fr/~dutour/Polyhedral/index.html for the relevant programs.

-

An implementation of randomized algorithms using hit-and-run and comparison with exact software presented in Bueler, Enge and Fukuda is presented here: http://homepages.ulb.ac.be/~vfisikop/publications/EF_socg14.pdf It works efficiently in high dimensions (say 100) but it only computes an estimation, not the exact volume.

-

There is a result by J. Borwein et al on volumes of convex polytopes constructed by intersecting symmetric hyperplanes. These volumes are obtained from multivariate integrals involving the function sin(x)/x. These integrals have very nice closed forms. Here are two references on this work:

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.40.8186

and

http://algo.inria.fr/seminars/sem01-02/borwein1.pdf

-

Since you're given so much information about the faces of the polytope, I think that the division into simplices method ought to be straightforward (if tedious) to carry out. The idea is to divide the polytope into orthoschemes, with signed volumes, as described below. This has the advantage that the formula for the volume of an orthoscheme is very simple, and computing the coordinates of the orthoschemes is simple linear algebra. The disadvantage is that there will be one orthoscheme for each flag of the polytope, so the number of them could be quite large. Computing the surface area of the boundary could be done in tandem, since one would just take the area of each orthoschemes intersection with the supporting hyperplanes. I suspect this method won't work for you though, if the polytope is complicated (and maybe this is the approach you already tried).

To compute the vertices of the orthoschemes, first choose a point, then take its orthogonal projection to each supporting hyperplane of the polytope. Inductively do this for each facet. The for each flag, take the corresponding points in each face of the flag, and form an orthoscheme. The sign of the volume of the orthoscheme will be determined in each dimension by whether the vertex lies inside or outside the corresponding hyperplane times the sign of the lower dimensional one it is a cone on. Some of the orthoschemes will lie partly outside of the polytope, but the volumes outside will cancel with this sign convention.

-

You could use a finite element method. If you stuff the surface with tetrahedra http://www.cs.berkeley.edu/~jrs/papers/stuffing.pdf‎ (pdf link not allowed), you could then sum the volumes of those tetrahedra.

-
Such a method cannot work efficiently in high dimensions, can it? You need (super-)exponentially many simplices to fill a cube, for example. –  Noah Stein Aug 21 '13 at 17:42