Incorrect volume calculation when H-rep contains redundant inequalities

[maxkapur]

I have a function that generates a bunch of polyhedra and computes their volume. All of the polyhedra are bounded by [0,1] in all variables (hence, their volume is finite and less than one), and I add additional constraints according to a certain rule. Sometimes, one of these constraints is completely redundant, and in such cases, volume() returns an overestimate of the volume.

MWE:

using Polyhedra

poly = polyhedron(HalfSpace([-1.0, 0.0], 0.0) ∩
                  HalfSpace([0.0, -1.0], 0.0) ∩
                  HalfSpace([1.0, 0.0], 1.0) ∩
                  HalfSpace([0.0, 1.0], 1.0) ∩
                  HalfSpace([-0.2, -0.8], -0.0) ∩
                  HalfSpace([0.6, 0.4], 0.6))

volume(poly)
# 1.1666666666666665

The vertices of this polyhedron, as correctly given by vrep(poly), are

 [ 1/3 , 1 ]
 [ 1   , 0 ]
 [ 0   , 0 ]
 [ 0   , 1 ]

and a simple sketch shows that the true area is 2/3.

• Polyhedra version: 0.6.13
• Julia version: 1.5.4

[schillic]

Confirmed in Julia v1.6.0 (and Polyhedra v0.6.13):

julia> using Polyhedra

julia> poly = polyhedron(HalfSpace([-1.0, 0.0], 0.0) ∩
                         HalfSpace([0.0, -1.0], 0.0) ∩
                         HalfSpace([1.0, 0.0], 1.0) ∩
                         HalfSpace([0.0, 1.0], 1.0) ∩
                         HalfSpace([-0.2, -0.8], -0.0) ∩
                         HalfSpace([0.6, 0.4], 0.6));

julia> volume(poly)
1.1666666666666665

# --- trying with LazySets

julia> using LazySets

julia> H = HPolytope(poly);

julia> area(H)
0.6666666666666666

# --- plotting the polygon

julia> using Plots

julia> plot(poly)  # plot(H) gives the same result

[maxkapur]

Here is another one whose area is miscalculated. The issue seems to arise when one of the redundant constraints just "kisses" a vertex. In this case, commenting out a redundant constraint that nowhere holds with inequality in the polyhedron still yields the incorrect volume (shown below), but commenting out the second constraint, which holds with equality at (1, 0), makes the area compute correctly.

using Polyhedra

poly = polyhedron(# HalfSpace([-1.0, 0.0], 0.0) ∩
                  HalfSpace([0.0, -1.0], 0.0) ∩ # Comment here to get correct area
                  HalfSpace([1.0, 0.0], 1.0) ∩
                  HalfSpace([0.0, 1.0], 1.0) ∩
                  HalfSpace([-0.6, -0.4], -0.6) ∩
                  HalfSpace([0.2, 0.8], 0.95))
@show volume(poly)
# volume(poly) = 0.6510416666666655
plot(poly, xlim=(0,1), ylim=(0,1))

[inechita]

I confirm this issue in Julia 1.9.

Here's another situation where volume returns an incorrect value:

using Printf
using Polyhedra
using Combinatorics
using GLPK; solver = GLPK.Optimizer
lib = DefaultLibrary{BigInt}(solver)
L3 = polyhedron(vrep([
    0 0 0 0 0 0 
    0 0 1 0 0 0 
    0 1 0 0 0 0 
    0 1 1 0 0 1 
    1 0 0 0 0 0 
    1 0 1 0 1 0 
    1 1 0 1 0 0 
    1 1 1 1 1 1 
    
]), lib)
[dim(L3) fulldim(L3) volume(L3)]

returns volume = 41//240, while the correct answer is 1//180, see this paper.

In both cases, volume seems to overestimate the correct value

[My-Kingdom-for-a-Cat]

I get the same issue with Julia 1.9.4 and GLPK or other libraries (e.g. CDDLib):

using Polyhedra
using CDDLib
lib=CDDLib.Library(:exact)
L3 = polyhedron(vrep([
    0 0 0 0 0 0 
    0 0 1 0 0 0 
    0 1 0 0 0 0 
    0 1 1 0 0 1 
    1 0 0 0 0 0 
    1 0 1 0 1 0 
    1 1 0 1 0 0 
    1 1 1 1 1 1 
    ]),lib)
volume(L3)

returns 41/240 instead of 1/180. I think this is due to a bug in the function triangulation_indices() (source code) used in volume() (source code), that finds a Delaunay simplicial decomposition of a polyhedron. In the example above, the function triangulation_indices() returns 123 simplices (defined by their indices), many of them being the same! Their indices are obtained for instance through

delaunay = triangulation_indices(L3)
toomanysimplices=[[index.value for index in simplexindices] for simplexindices in delaunay]

which returns a 123-element Vector{Vector{Int64}} containing (truncating the output)

 [1, 2, 3, 4, 5, 6, 7]
 [1, 2, 3, 4, 5, 6, 7]
 [1, 2, 3, 4, 6, 7, 8]
 [1, 2, 3, 4, 6, 7, 8]
 [1, 2, 3, 4, 6, 7, 8]
 [1, 2, 3, 4, 6, 7, 8]
 ⋮
 [1, 3, 4, 5, 6, 7, 8]
 [1, 3, 4, 5, 6, 7, 8]
 [1, 3, 4, 5, 6, 7, 8]
 [1, 3, 4, 5, 6, 7, 8]

while in this list there are only four distinct list of indices. These are obtain through union(toomanysimplices) which returns

4-element Vector{Vector{Int64}}:
 [1, 2, 3, 4, 5, 6, 7]
 [1, 2, 3, 4, 6, 7, 8]
 [1, 2, 4, 5, 6, 7, 8]
 [1, 3, 4, 5, 6, 7, 8]

To compute the correct volume, one find the (non-redundant) list of simplices through

simplices=map(Δ -> vrep(get.(L3, Δ)), union(triangulation_indices(L3)))

and the corresponding volume is computed through

reduce(+, unscaled_volume_simplex(Δ) for Δ in simplices; init=0) / factorial(fulldim(L3))

which returns 1/180, as announced by inechita.

The function volume() yields wrong results on other examples: for instance, if P is a bounded convex polytope in high dimension, the volume of P is found to be smaller than the sum of the volumes of P∩H⁺ and P∩H⁻ where H⁺ and H⁻ are complementary halfspaces. A quick workaround is to replace

    return map(Δ -> vrep(get.(p, Δ)), triangulation_indices(p))

by

    return map(Δ -> vrep(get.(p, Δ)), union(triangulation_indices(p)))

on line 50 of triangulation.jl. This solve the inconsistencies above, but I presume there is better to do, i.e. to ensure that triangulation_indices() does not generate redundant simplices indices.