# Faces & Vertices of N Linear Constraints

Here’s a fun one, so instead of doing what I normally do(hacking away at code), I figured I’d ask if people had tools to recommend for solving a problem. I know you people working in optimization and geometry live in this world, so my guess is there’s a bunch of tools to make life easier.

Problem statement. I have M linear inequalities pertaining to a finite volume of an N dimensional space. How do I find the vertices of the shape, and the barycenters of the faces?

For example:
0 < X < 1.
0 < Y < 1.
Makes a square. (4 vertices, 4 faces)

But we could
0 < X < Y
Now we have a traingle (3 vertices, 3 faces)

But now what if I have 20 dimensions and 30 constraints, where are my vertices?

I can imagine a naive approach to doing this… But, it would likely be pretty slow, and I’d like to hear how others could approach the problem (as easily as possible) using the julia ecosystem. Turns out this is a useful thing for a missing component to our ecosystem, but also - a fun algorithm problem.

To get the vertices and rays of the polytope from the inequality representation, you can use Polyhedra.jl. For the barycenters, please let me know if you figure something out. I know that the Chebyshev center of the polytope can be obtained via convex optimization for example. I would first retrieve the facets of the polytope, and then try to formulate the barycenter as an optimization.

If you are interested in computational geometry, take a look at Meshes.jl. I am actively working there.

2 Likes

Cool! Never seen polyhedra.jl. I’ll take a crack at it using it - might ask some questions along the way. Also will likely lead to a PR to another package.

I just found this paper, didn’t get a chance to read it yet: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.579.9132&rep=rep1&type=pdf “Analysis of backtrack algorithms for listing all vertices and all faces of a convex polyhedron”. Computational Geometry Volume 8, Issue 1, June 1997, Pages 1-12.

Yes I think that the naive thing is best for the barycenter calculations. Quickly determining the facets is a neat problem too.

I am very interested in computational geometry, my biggest issue is - finding time to formally learn a lot of it. I’ll star Meshes.jl though and keep track of it.

1 Like

Thanks for sharing the papers. My research is requiring a lot of computational geometry lately, so I have no excuse to find the time. 1 Like

Hey the good thing about comp. geom. is it get’s used everywhere. Graphics, physics, data stuff, stat’s, pure math/opt, etc. So having a backbone in it will only benefit you.

So I got a MWE for how to get the vertices, really happy with how this comes together:

``````using Polyhedra, CDDLib, JuMP

model = Model()

@variables(model, begin
x
y
end)

con = @constraints(model, begin
0. <= x <= 1.0
0. <= y <= 1.0
x <= y
end)

poly = polyhedron(model, CDDLib.Library(:exact))
points(vrep(poly))
``````
2 Likes

any idea why I can’t get rays/lines for the example above? I believe it’s due to the use of `<=` and `>=` but if you swap them for just the `<`/`>` oprators everything crashes.

I am not using Polyhedra.jl, I think you will have a better chance opening an issue there.

1 Like

To get each facet, you can do

``````for hs in halfspaces(poly)
hp = Polyhedra.hyperplane(hs)
facet = intersect(poly, hp)
end
``````

You can get the Chebyshev center with `chebyshevcenter(facet)` but the barycenter is not implemented.
The backtracking algorithm you are mentioning is implemented in https://github.com/JuliaPolyhedra/LRSLib.jl

4 Likes

I just tried the chebryshevcenter and it gave me an error but - that’s okay. I just whipped up a little script:

``````plane_point_map = Dict()
for (h,halfspace) in enumerate(halfspaces(hrep(poly))), (p,point) in enumerate(points(vrep(poly) ) )
if ( in(point,Polyhedra.hyperplane(halfspace)) )
push!(plane_point_map[h], p)
else
plane_point_map[h] = [p]
end
end
end
``````

gives the points on each plane/facet. So basically you could take the column means of the coordinates(with > 1 point) and get the barycenter of each face? Not a “fast” algorithm but - speed isn’t a huge concern yet and it could be put into embarassingly parallel no problem.

I’ll try to switch to the LRSLib backend and see if I notice anything different. The CDDLib backend is working great so far, but is a little slow to load on the first go.

Looks like I was able to get what I wanted. Code quality isn’t top notch but I’m just playing with an idea for now.

``````using Polyhedra, CDDLib, JuMP, Statistics

model = Model()

@variables(model, begin
x
y
end)

@constraints(model, begin
0 <= x <= 1.0
0 <= y <= 1.0
x <= y
end)

print(model)
poly = polyhedron(model, CDDLib.Library())

plane_point_map = Dict()
for (h,halfspace) in enumerate(halfspaces(hrep(poly))), (p,point) in enumerate(points(vrep(poly) ) )
if ( in(point,Polyhedra.hyperplane(halfspace)) )
push!(plane_point_map[h], p)
else
plane_point_map[h] = [p]
end
end
end
plane_point_map

pnts = hcat(points( vrep( poly ) )...)'
vars = last(size(pnts))
barycenters = []
for (_,v) in plane_point_map
if length( v ) > 1
push!( barycenters, mean(pnts[v,:], dims = 1) )
end
end

barycenters = vcat( barycenters... )
pnts

using Plots
scatter(pnts[:,1], pnts[:,2], legend = false)
scatter!(barycenters[:,1], barycenters[:,2], legend = false)
``````

You can use `incidentpoints` to simplify the code.
CDD computes the incidence information and `incidentpoints` just retrieve this information. In your code sample, you recompute it.

``````using Polyhedra, CDDLib, JuMP, Statistics

model = Model()

@variables(model, begin
x
y
end)

@constraints(model, begin
0 <= x <= 1.0
0 <= y <= 1.0
x <= y
end)

print(model)
poly = polyhedron(model, CDDLib.Library())

removehredundancy!(poly)
barycenters = [mean(incidentpoints(poly, hidx)) for hidx in eachindex(halfspaces(poly))]
``````

I get the output

``````Feasibility
Subject to
x - y ≤ 0.0
x ∈ [0.0, 1.0]
y ∈ [0.0, 1.0]
3-element Array{Array{Float64,1},1}:
[0.5, 0.5]
[0.0, 0.5]
[0.5, 1.0]
``````

Also: try to avoid things like `hcat(points( vrep( poly ) )...)` or `vcat( barycenters... )`. If you have millions of points, it will compile a `hcat` or `vcat` with millions of arguments which will be quite slow to compile !

1 Like

I did try to dereference/unsafe_load the pointers from the representation to get a 2 array, but there’s some issues with that. I think I can use them just fine as a list of 1-arrays, I do like what you did for the barycenters.

Thank you for sharing this snippit this is a really nice example of how to use the outputs from the types in the package. I think I’m going to use this code to build the skeleton of a package today. If you replied in this thread, please message me if you want a special form of credit otherwise I will give a shout out in the readme. I’ll let you all know when the git is made and all of that.

Here’s what this has helped me do. https://github.com/caseykneale/ExtremeVertexDesigns.jl

Not a released package - yet, but, it’ll be a workspace for some design of experiments things to accompany: https://github.com/phrb/ExperimentalDesign.jl

Nice package, good to see the `incidentpoints` features in action! Thanks for the credits in the README 1 Like

A curiosity: I have never used CDDLib before and I’m wondering how far I could push it as a teaching aid. For example could it handle 50 or 100 or a couple of hundred constraints and demonstrate to students the task that LP solvers must confront (in terms of number of vertices of the polyhedron)?

Míle buíochas.

1 Like

It really depends on the dimension, you can handle billions of constraints in 2D as there is an n log n algorithm in that case. For higher dimension, it depends also on the number of vertices you expect to have. See http://cgm.cs.mcgill.ca/~avis/doc/avis/ABS96a.ps for a discussion on this

1 Like