[EM] Exact volume of a simplex

[Kristofer Munsterhjelm]

On 08.02.2022 01:17, Daniel Carrera wrote:
> 
> 
> On Mon, Feb 7, 2022 at 12:16 PM Forest Simmons
> <forest.simmons21 at gmail.com <mailto:forest.simmons21 at gmail.com>> wrote:
>>     Recursively triangulate the boundary of C into a union beta of (N-1)
>>     dim'l simplices.
>>     Then let p be any point of the interior of C, and cone every member
>>     of beta over p. Then C is the union of the resulting N dimensional
>>     simplices.
> 
> 
> I'm pretty new to the field, but that sounds like it should work. I
> can't imagine how it wouldn't. I can already get a list of vertices from
> the Polytope library and I can give those to QHull (via scipy.spatial)
> and QHull can give me a list of simplices that make up the faces.
> Lastly, to get p we can grab any two vertices that do not share a face
> and pick the midpoint. It's hard to see how that could go wrong, but
> then again, what do I know?
> Kristofer pointed out that my previous algorithm wouldn't work
> efficiently in cases where most sets of N+1 vertices share a face (e.g.
> a square base pyramid in 3D). But your idea has no such problem as far
> as I can tell. I'm sure that the general problem of dividing up a
> polytope into simplices is very hard, but the specific case of a Voronoi
> cell honestly looks easy.
> 
> If there is a reason why it is better to divide the Voronoi cell into as
> few simplices as possible, a variation of my algorithm could work well.
> Instead of grabbing vertices at random you'd have to do something more
> clever to avoid grabbing vertices that don't all share the same face,
> but that doesn't sound hard.

Computational geometry isn't my field either, so I can't tell why a
particular method will or won't work. But from what I understand, the
problem of determining the volume of a general convex polytope if you
only have one of its dual representations (the half-space representation
or the vertex representation) is #P-hard, which is the counting analog
of NP-hard.[1]

So my heuristic is basically "unless you've just coincidentally solved a
Millennium Prize problem, your algorithm probably won't work, or it
won't be quick".

But this shouldn't matter in practice. To quote another paper:

> It is important to note that the algorithms differ in the input they 
> require; because the transformation V<->H can be prohibitively costly
> for some polytopes, this can influence decisively the choice of
> algorithms. Though highly depending on the concrete structure of the
> polytope under consideration, one can expect tractability on 
> workstations for polytopes with up to 10^2 hyperplanes and 10^4 
> vertices, or, 10^4 hyperplanes and 10^2 vertices, in dimension up to 
> 15, requiring memory up to the range of 10^2 - 10^3 MB ...
So although the problem is asymptotically nasty, it starts off gentle
and we shouldn't have a problem for our particular domain... though
something more optimized than Python might be required for 7-10 dimensions.

There might also be the case that some algorithm works nicely for
Voronoi polytopes (or intersections of them), but not for general
polytopes - that's one way your algorithm could possibly evade the
"pretty much sure this is impossible" bound.

(I think part of the problem lies in going from V-representation to
H-representation and back; this can take a lot of time and space.)

-km