Getting Extremal Rays of Cone

Question

So I have a set of linear homogeneous equations $A\vec{x}=0$. I would like to solve this for non-negative solutions. I can solve the system in general and I get the two vectors that span the solution space, but one of the vectors has negative entries. I would like to restrict the solution space to $\mathbb{R}_{\geq0}^n$. I should add that the vector $\vec{x}$ is a vector of symbolic variables (Sympy), not sure if this is important.

I have been told that one way to do this is to use the simplex method of linear programming and set the objective function to be the zero function. When I do this in SciPy I only get one solution vector and not the two I want (the output seems to be the sum of the two I want).

Questions

• How could I go tweak the algorithm so that it spits out the vectors (extremal rays) that span the solution space (cone)?

• Are there any other libraries that I can use that will do this for me in Python?

Example Let's say I have the system:

1. $\lambda_0 - \lambda_3 = 0$
2. $\lambda_1 + \lambda_4 - \lambda_0 = 0$
3. $\lambda_2 - \lambda_1 = 0$
4. $\lambda_3 - \lambda_2 - \lambda_4 = 0$

Solving this system is equivalent to finding the nullspace of,

\begin{bmatrix} 1&0&0&-1&0\\ -1&1&0&0&1\\ 0&-1&1&0&0\\ 0&0&-1&1&-1\\ \end{bmatrix}

Using a SciPy method to find the nullspace of the matrix gives me the solution space spanned by the vectors, $$ [1, 1, 1, 1, 0] , [0, -1, -1, 0, 1]$$ I would like to restrict this to non-negative solutions. I believe the two vectors that span the solution space restricted to the positive reals is, $$ [1, 1, 1, 1, 0] , [1, 0, 0, 1, 1] $$ So far I have tried to use the linear programming method with the simplex method by setting the objective function to the zero function. The output there is, $$ [2, 1, 1, 2, 2]$$ This is the sum of the two I want. The least-squares method gives the same result under bounds (0.5,2).

Answer 1

You can use a randomized approach to solve this, as I mentioned in the comments. Ultimately, we have the equality constraint for our solution vector $\lambda$ being $A \lambda = \boldsymbol{0}$, for some given $A$. Since we are considering elements of the nullspace, their scale does not really matter but what we do care about is that the elements of $\boldsymbol{\lambda}$ are greater than or equal to $0$. Thus, let us set the box constraints $0 \leq \lambda_i \leq 1$, for all $i$.

Now, we assume there are potentially multiple solutions that satisfy these constraints, so how can we potentially identify them? Let us seek to minimize an objective of the form $c^T \lambda$ where we ensure that our choice for $c$ is such that $c_i < 0$ for all $i$ to ensure we do not get the trivial solution $\lambda = \boldsymbol{0}$ under our constraints.

Suppose then we use a randomized approach, perhaps choosing for all $i$ that $c_i \sim \text{unif}(-1,0)$ and then generating a set of $k$ random vectors $\lbrace c^{(1)}, \cdots, c^{(k)}\rbrace$. If we solve the $k$ linear programs minimizing $\left(c^{(j)}\right)^T \lambda$ under our constraints, for sufficiently large $k$ we should obtain all the solutions that you desire with high probability. I'm sure one can work out some bounds using Hoeffding's inequality or something along those lines and in turn figure out the sample complexity needed to ensure you find all the solutions with high probability, but ultimately this should work.

I coded up a sample Matlab script to illustrate this algorithm

% test randomized algorithm idea for finding nullspace with constraints

% setup the linear program values
A = [1, 0, 0, -1, 0;
    -1, 1, 0, 0, 1;
    0, -1, 1, 0, 0;
    0, 0, -1, 1, -1];
b = zeros(4,1);
lb= zeros(5,1);
ub= ones(5,1);

% solve k linear programs with randomized approach
k = 15;
X = zeros(5,k);
for i = 1:k
X(:,i) = linprog(-rand(5,1), A, b, [], [], lb, ub);
end

% find the unique solutions
U = unique(X','rows')'

For $k = 15$, I have run the script many times and it always produces the correct two solutions, namely producing the output

U =

     1     1
     0     1
     0     1
     1     1
     1     0

If I use $k = 2$, sometimes it will produce both, sometimes it does not. You can play around with your approach, perhaps sampling until you have as many unique solutions as you expect, but the point is this algorithm will give you the results you are looking for with enough samples.