# How to find all vertices of a polyhedron

I have a convex polyhedron given by a set of linear inequalities, for example:

$$x_1 \geq 0,~~ x_2 \geq 0, ~~x_3\geq 0 \\ x_1+x_2\leq 1,~~ x_2+x_3\leq 1,~~ x_3+x_1\leq 1$$ I want to list all the extreme points of the polyhedron. In this case, these points would be: $$(0,0,0),~~(1,0,0),~~(0,1,0),~~(0,0,1),~~(1/2,1/2,1/2)$$

In python, there are several linear programming libraries, such as scipy.linprog or cvxpy, that can return one such extreme point using the Simplex method. But I want to list all of them. How can I do this?

• Take a look at the double description method (link.springer.com/chapter/10.1007/3-540-61576-8_77). – batwing Jul 17 at 15:25
• AFAIK, SCIP solver has a facility to enumerate all feasible solutions. please, see this link. Also, modern solvers like CPLEX or Gurobi has a solution pool facility to do that. – A.Omidi Jul 18 at 9:41
• @A.Omidi, I think all three of those solvers enumerate only the equivalence classes of mixed integer solutions, where two solutions are considered equivalent if they agree on the integer variables. In particular, they would not return $(1/2,1/2,1/2)$ if you define any of the $x_i$ as integer, and they would return only a single solution otherwise. – RobPratt Jul 18 at 18:24
• @RobPratt, thanks. You are right. What I mentioned is especially about MILPs. About the LPs, it is a bit complicated. For a system of equations polymake would be interesting. – A.Omidi Jul 20 at 5:12
• This is not a proper answer as it is not guaranteed to enumerate all possible vertices. However, in a pinch I often use a “pea-shooter experiment”. Sample a vector $c$ from the unit sphere in $n$ dimensions, and solve the resulting LP with objective $c$. Record the optimal solution. Repeat. After a while this stops producing new solutions. If you warm-start each LP solve with the previous solution this can work very quickly. A fun programming exercise, at least. – David M. Jul 25 at 0:28

The problem of enumerating all vertices of a polytope has been studied, see for example Generating All Vertices of a Polyhedron Is Hard by Khachiyan, Boros, Borys, Elbassioni & Gurvich (available free online at Springer's website) and A Survey and Comparison of Methods for Finding All Vertices of Convex Polyhedral Sets by T. H. Matheiss and D. S. Rubin. That's a pretty old survey though (1980), so newer methods may be available.

A naive brute force approach can be deduced from the definition of vertex/extreme point. Let's call the polytope $$P$$. Pseudocode can be as follows:

1. Select a subset of $$n$$ inequalities (in you example $$n = 3$$), getting a smaller linear system of inequalities with submatrix $$A'$$ and vector $$b'$$.

2. Solve the linear system $$A'x = b'$$. There are three cases here:

a. The system has no solution: Then, return to (1) and select another subset (not chosen before).

b. The system has no unique solution: Then, $$A'$$ is linearly dependent. Return to (1) and select a new subset.

c. The system has a unique solution: If that solution is feasible for $$P$$, then it's a vertex. Go back to (1).

The algorithm ends when no new subsets can be chosen. Note that different subsets of rows could yield the same vertex.

A second alternative can be treating the polyhedron's vertices and edges as a graph (might work faster than the brute force solution above):

1. Start at any vertex $$x$$ of the polytope. For example, the one you found using the Simplex, Interior Point or Ellipsoid method with some cost function.
2. Find all $$P$$'s edges incident to $$x$$. That is, all 1-dimensional faces of $$P$$. This can be done similar to pivoting on nonbasic variables (with respect to the current vertex). Note that vertices are 0-dimensional faces of $$P$$.
3. Explore this graph (with the analogy of vertices and edges) using breadth-first search or depth-first search.

As @batwing mentioned, another alternative is using the Double Description Method by Motzkin et al. to generate all extreme points and extreme rays of a general convex polyhedron represented as a system of linear inequalities $$Ax \leq b$$. An implementation called cdd can be found at Komei Fukuda's website here, while this GitHub repo contains pycddlib, a Python wrapper to interact with that library. Finally, at this repo the package pypoman is developed to interact with the Python wrapper to get the extreme points for $$Ax \leq b$$ starting from $$A$$ and $$b$$.

• Wonderful answer, thanks! Do you know of any Python packages that implement one of these more advanced algorithms? – Erel Segal-Halevi Jul 18 at 18:16
• I'm not aware of python packages that implement those solutions. But I edited adding a third procedure, the double description method, with python wrappers to implementations available by the authors. – dhasson Jul 18 at 22:24
• @dhasson If it's not too much, could you give an explanation as to why vertices of polyhedra in $\mathbb{R}^n$ correspond to solutions of non-singular subsystems of $n$ equations? – Blue Oct 13 at 13:45
• @Blue A proof of this equivalence can be found in section 2.2 of Introduction to Linear Optimization by Bertsimas and Tsitsiklis. They show that the definitions of (a) vertex, (b) extreme point and (c) basic feasible solution are equivalent. a) => b) is simple; b) => c) is deduced from showing that if $x$ is not a basic feasible solution, then it can be written as the convex combination of 2 points in the polyhedron; some algebraic manipulation yields c) => a) – dhasson Oct 13 at 14:36
• @dhasson Thank you so much! – Blue Oct 13 at 15:00

You obtain all vertices of a polytope using polymake.

You can directly try the online version.

New contributor
Graph4Me Consultant is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.

It seems to me that cdd libraries can be useful to solve this problem. Description is available at cdd. There is an implementation of this function in R: rcdd. You can use the following instruction to solve this problem:

install.packages("rcdd")
require(rcdd)
scdd(makeH(rbind(-diag(3),c(1,1,0),c(0,1,1),c(1,0,1)),c(rep(0,3),rep(1,3))))