How to convert the half-spaces that constitute a convex hull to a set of extreme points?

Question

I have a convex set in a Euclidean space (3D, but would like answers for nD) that is characterized by a finite set of half-spaces (normal vector + point).

Is there a better algorithm to find the extreme points of the convex set other than compute brute force all points that are intersections of 3 (or, n) half-spaces and eliminate those that are not extreme points?

Acklin

Do you want to find all the extreme points, or just some subset of them? – Queenstown

If I got the theory right, to define the convex set I need all extreme points. Depends on the exact definition of extreme points. I'm thinking of an extreme point as a point that can not be obtained by p= p0 * t + p1*(1-t) for 0<= t <=1 and p0 !=p1, both within the convex shape. Put in other words, I want the minimal set of points that generate the convex set. – Acklin

I see, there could be degenerated cases ... . Edit: thinking twice, I don't see clearly, not immediately. – Acklin

It sounds like you want the convex hull of the polygon, except that instead of given the points, you're given the half-planes. Is that correct? – Queenstown

Exactly. Instead of the points that I need, I have the other representation using half-spaces. (Would be half-planes in 2D). – Acklin

link.springer.com/article/10.1007%2FBF02293050 seems relevant – Huntsman

The con2vert Matlab script implements this using the primal-dual polytope method; the script isn't very long. – Chichihaerh

@Normal thank you, from the description it looks exactly like what I was looking for. Alas, when in click at "Watch this File" I'm asked to register. Won't leave my DNA everywhere easily. Is this patented or licensed thinking, or is it free? – Acklin

"Watch this file" means to be notified of new releases. The download button is to the right. However, there is no license. I see that this file has "inspired" a later submission, which is BSD-licensed and appears to supersede the functionality of the former. – Chichihaerh

@Normal, thank you, it looks like I've been just a little bit over-paranoid about all those clowns trying to charge everyone for about anything. If you could give a little help (link) on how to read this (specially "c=A\b"), I would accept this as an answer. – Acklin

Answer 1

The key term is vertex enumeration of a polytope P. The idea of the algorithm described below is to consider the dual polytope P*. Then the vertices of P correspond to the facets of P*. The facets of P* are efficiently computed with Qhull, and then it remains to find the vertices by solving the corresponding sub-systems of linear equations.

The algorithm is implemented in BSD-licensed toolset Analyze N-dimensional Polyhedra in terms of Vertices or (In)Equalities for Matlab, authored by Matt J, specifically, its component lcon2vert. However, for the purpose of reading the algorithm and re-implementing it another language, it is easier to work with the older and simpler con2vert file by Michael Kleder, which Matt J's project builds on.

I'll explain what it does step by step. The individual Matlab commands (e.g., convhulln) are documented on MathWorks site, with references to underlying algorithms.

---

The input consists of a set of linear inequalities of the form Ax<=b, where A is a matrix and b is a column vector.

Step 1. Attempt to locate an interior point of the polytope

First try is c = A\b, which is the least-squares solution of the overdetermined linear system Ax=b. If A*c<b holds componentwise, this is an interior point. otherwise, multivariable minimization is attempted with the objective function being the maximum of 0 and all numbers A*c-b. If this fails to find a point where A*c-b<0 holds, the program exits with "unable to find an interior point".

Step 2. Translate the polytope so that the origin is its interior point

This is done by b = b - A*c in Matlab. Since 0 is now an interior point, all entries of b are positive.

Step 3. Normalize so that the right hand side is 1

This is just the division of ith row of A by b(i), done by D = A ./ repmat(b,[1 size(A,2)]); in Matlab. From now on, only the matrix D is used. Note that the rows of D are the vertices of the dual polytope P* mentioned at the beginning.

Step 4. Check that the polytope P is bounded

The polytope P is unbounded if the vertices of its dual P* lie on the same side of some hyperplane through the origin. This is detected by using the built-in function convhulln that computes the volume of the convex hull of given points. The author checks whether appending zero row to matrix D increases the volume of the convex hull; if it does, the program exits with "Non-bounding constraints detected".

Step 5. Computation of vertices

This is the loop

for ix = 1:size(k,1)
    F = D(k(ix,:),:);
    G(ix,:)=F\ones(size(F,1),1);
end

Here, the matrix k encodes the facets of the dual polytope P*, with each row listing the vertices of the facet. The matrix F is the submatrix of D consisting of the vertices of a facet of P*. Backslash invokes the linear solver, and finds a vertex of P.

Step 6: Clean-up

Since the polytope was translated at Step 2, this translation is undone with V = G + repmat(c',[size(G,1),1]);. The remaining two lines attempt to eliminate repeated vertices (not always successfully).

Chichihaerh

Answer 2

I am the author of polco, a tool which implements the "double description method". The double description method is known to work well for many degenerate problems. It has been used to compute tens of millions of generators mostly for computational systems biology problems.

The tool is written in Java, runs in parallel on multicore CPUs and supports various input and output formats including text and Matlab files. You will find more information and publications about the software and the double description method via given link to a university department of ETH Zurich.

Walkerwalkietalkie