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?
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).