# Generating all extreme rays

I am trying to understand a problem and would like to generate all extreme rays for a given set of linear constraints. With the Python interface of CPLEX, I was able to generate a single ray (not sure if it is guaranteed to be extreme) but is there a way to get all extreme rays?

On a theoretical level, I know that in my problem, 0 is always a feasible solution and if $$x$$ is feasible and $$\lambda \ge 0$$ then $$\lambda x$$ is feasible (I think this means that the constraints describe a pointed polyhedral cone). So as far as I understand that means that there is a ray through every feasible point. But how do I find the extreme ones out of this infinite set?

Edit Thanks to the helpful comments below, I learned about the double description algorithm and found implementations in sagemath.org and polymake. Unfortunately, I had problems with both implementations. The first could not handle a matrix with a rank smaller than the number of columns. The second returns a set of rays but I think it must be missing some because there are solutions to the original constraints that cannot be expressed as a non-negative combination of the rays. I'm not certain about using polymake correctly, so let me give an example:

Example (in polymake)

$$inequalities=new Matrix([ [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0,-1, 0, 0, 0, 0], [0, 1, 0, 0, 1,-1, 0, 0, 0], [0, 0, 1, 0, 0, 1,-1, 0, 0], [0, 0, 0, 1, 0, 1, 0,-1, 0], [0, 1, 0, 0, 0, 0, 1,-1, 0], [0, 0, 0, 0, 0, 0, 0, 1,-1]]);$$p=new Polytope<Rational>(INEQUALITIES=>$$inequalities); print_constraints($$p->INEQUALITIES);
0: x4 >= 0
1: -x4 >= 0
2: x1 + x4 - x5 >= 0
3: x2 + x5 - x6 >= 0
4: x3 + x5 - x7 >= 0
5: x1 + x6 - x7 >= 0
6: x7 - x8 >= 0
7: 0 >= -1
# All variables are general, i.e., not restricted to non-negative values.

print $p->VERTICES; 1 0 0 0 0 0 0 0 0 0 1 -1 -1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 -1 0 0 0 0 1 -1 0 0 1 0 1 0  Now, I think the first coordinate indicates whether the following row vector is an extreme point (1) or ray (0). Since this is a pointed cone, it makes sense that the only extreme point is 0 and all other rows are rays. However, I thought that every solution to the constraints should be expressible as a non-negative combination of these rows. But there are solutions with $$x_8 = 1$$ even though all rays have $$x_8 = 0$$. For example, $$(1, 0, 0, 0, 1, 1, 1, 1)$$ should be a solution. Can someone help me understand this? Second Edit The trick with polymake was that rays $$r$$ where $$r$$ and $$-r$$ are extreme rays are not part of RAYS. Instead they are given as LINEALITY_SPACE. Including them and their negative values gave me the answer I was looking for. print$p->LINEALITY_SPACE;
0 -1 2 1 0 -1 1 0 0
0 1/2 0 1/2 0 1/2 1/2 1 1

• What is the dimension of your set? If it is not "too big" then you should be googling "double description algorithm". A list of codes that do polyhedral computation is at: cs.mcgill.ca/~fukuda/soft/polyfaq/node41.html. If the dimension is large, then you probably can't enumerate all extreme rays, there are likely to be too many unless your set is specially-structured. Aug 15, 2019 at 19:34
• @JeffLinderoth Thanks for the pointer. I found a description and even an implementation on sagemath.org. My example is hopefully small enough to try it: the matrix A 7x8 but doesn't have full rank. Is this a problem for the algorithm? (Here is my attempt on sagemath: tinyurl.com/y6jdv4hv) Aug 15, 2019 at 23:55
• Maybe you could try to use polymake polymake.org - they also have various functions to compute different representations. See forum.polymake.org/viewtopic.php?t=24 for some example. Aug 16, 2019 at 7:39

What is the dimension of your set? If it is not "too big" then you should be googling "double description algorithm". A list of codes that do polyhedral computation is at: https://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node41.html. If the dimension is large, then you probably can't enumerate all extreme rays, there are likely to be too many unless your set is specially-structured.

• Including laying off beer for a month?
– prubin
Aug 17, 2019 at 18:33
• Thank you. Using Polymake, I was able to find the rays I was looking for. There is a subtle issue that RAYS do not contain all of the extreme rays but some of them are in LINEALITY_SPACE. I'll update my example in the question above. Aug 19, 2019 at 7:53

I'm assuming that your variables ($$x$$) are nonnegative. If you take a cross-section of the cone by adding a constraints such as $$\sum_i x_i = 1$$, you get a polytope, and I believe that there is a 1-1 correspondence between extreme rays of the original cone and extreme points of the polytope. IIRC, there are programs for computing all extreme points of a polytope. You can also try a Monte Carlo approach, which will not provably find all of them but will likely find most if not all.

I tried the Monte Carlo approach on your algebraic representation (constraints 0 through 6) above. What I did was minimize a randomly weighted combination of the variables, then maximize the same weighted sum, then repeat with new weights, until a certain number of consecutive LPs had been solved without finding a new solution. (I set the limit as 1,000 consecutive failures. CPLEX hit the limit after solving 1,370 LPs in under a second. Using this method, I found what I believe are 16 extreme rays:

[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0]
[0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0]
[0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0]
[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
[1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0]
[1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0]
[1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0]
[1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0]
[1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]
[1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0]
[1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0]
[1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 2.0, 0.0]
[1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 2.0, 2.0]


Cranking the failure limit to 10,000 did not alter the results (and, somewhat scarily, still only needed about a second, with output turned off).

• Than you for the pointers. My variables are unfortunately all general, so I don't think the additional constraint would be correct. The Monte Carlo approach sound like it would also work with general variables. What I didn't understand is how to prevent the LP from giving the same answer every time. Is there a way to exclude a discovered extreme ray? Aug 19, 2019 at 7:49
• The added constraint will work with an polar cone with vertex 0, and it is totally arbitrary. You just want to take a "slice" across the cone. You could set the right side to any positive value, and you could use a different set of nonnegative coefficients (but all ones is easiest). What keeps you from getting the same solution every time is changing the objective coefficient vector each time (to a different random vector). You do get repetitions ... lots and lots of them ... but you will not get just a single solution over and over (unless your random number generator really hates you).
– prubin
Aug 19, 2019 at 22:07
• The dimensions in your example are low enough that the rays (more precisely, the extreme points of the cross section I created) can be enumerated explicitly. I cobbled together an R notebook to do so, and wound up with the same 16 rays I got from the Monte Carlo method. The workload for the brute force approach is a combinatorial function of the dimension (8 here) and number of constraints (7 here), so it won't work for much higher dimensions.
– prubin
Aug 19, 2019 at 22:47
• Thank you again for spending your time on this. What still confuses me is that the results of the MC method and the ones from polymake (see my edited question) don't match. For example, polymake finds a ray $(0, -1/2, 0, -1/2, 0, -1/2, -1/2, -1, -1)$ and these negative values cannot be a non-negative combination of the 16 rays found by the MC method. I also don't see how the constraint $\sum_i x_i = 1$ can generate a slice if the variable are general. Shouldn't it be possible to non-negatively scale every solution until it "hits the slice"? I don't see how to do this with a negative solution. Aug 21, 2019 at 9:14
• I've been assuming your variables are nonnegative. In that case, $(0,−1/2,0,−1/2,0,−1/2,−1/2,−1,−1)$ cannot be a ray but its negation $(0, 1/2, 0, 1/2, 0, 1/2, 1/2, 1, 1)$ could be ... and is (double it and you get the last ray I listed above). As far as hitting the slice, multiply by 1/4 and you're there.
– prubin
Aug 22, 2019 at 14:59