Generating all extreme rays

Question

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<Rational>([
    [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

Answer 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: 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.

Answer 2

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