Find a nonnegative basis of a matrix nullspace / kernel

Question

I have a matrix $S$ and need to find a set of basis vectors $\{\mathbf{x_i}\}$ such that $S\mathbf{x_i}=0$ and $\mathbf{x_i} \ge \mathbf{0}$ (component-wise, i.e. $x_i^k \ge 0$).

This problem comes up in networks of chemical reactions where $S$ is the stoichiometric matrix and $\mathbf{x_i}$ are extreme vectors that define a proper, polyhedral cone at steady state (I don't understand entirely what all of this means - have a picture in my head though).

More detail: Suppose you're looking at a concentration network involving $n$ chemical species , summarized in a vector $\textbf{u}$, and $m$ reactions with rates $\textbf{w}(\textbf{u})$. Then an ODE that describes the behaviour of this dynamical system is $$\frac{d \textbf{u}}{dt}=S\textbf{w}(\textbf{u}).$$ The stoichiometric matrix $S \in \mathbb{Z}^{n \times m}$ describes how much of the species involved in a reaction is consumed / produced relative to one other.

At steady state, $d\textbf{u}/dt|_{\textbf{u}=\textbf{u}^*}=0$, we now look for the kernel / null space of $S$, i.e. reaction rates $\textbf{w}$ such that $S\mathbf{w}=0$.

Apparently (however I don't fully understand this), the intersection of $\text{ker}(S)$ and $\mathbb{R}_+^m$ forms a proper, polyhedral cone. This cone can be represented by a nonnegative combination of the finite set of extreme vectors that are unique up to scaling by a positive constant - the $\textbf{x}_i$ above are these extreme vectors.

(I quoted the latter bit mostly verbatim from http://www.springerlink.com/content/g5093r283160p518/)

Possible route of solving this:

So there appears to be a way to define this as a linear programming (LP) problem, but I don't quite see this. This seems to be suggested here: http://r.789695.n4.nabble.com/convex-nonnegative-basis-vectors-in-nullspace-of-matrix-td4548822.html

Any elaboration on the LP approach or any other way of solving this is greatly appreciated.

Answer 1

First note that $ker(S)$ is a linear subspace of $\mathbb{R}^m$ (I am assuming $S$ is an $n \times m$ matrix). Next observe that $\mathbb{R}_+^m$ is a polyhedral cone (a cone is polyhedral if it can be expressed as $Ax \leq 0$). The intersection of a linear subspace and the cone $\mathbb{R}_+^m$ will be another polyhedral cone. It is worth trying to visualize this. Consider $\mathbb{R}^3$. Then the non-negative orthant $\{(x,y,z) \mid x\geq 0, y\geq0, z\geq0\}$ is a polyhedral cone. Another description of this cone is that it is the convex hull of the following three extreme rays $\{t(1,0,0) \mid t \geq 0\}$, $\{t(0,1,0) \mid t \geq 0\}$, and $\{t(0,0,1) \mid t \geq 0\}$. Now you should be able to convince yourself that the intersection of the non-negative orthant with any hyperplane passing through the origin will result in another polyhedral cone.

Polyhedral cones have a finite number of extreme rays. There are various methods to enumerate these extreme rays. You can truncate the cone and enumerate the extreme points of the resulting polytope; this seems to be the approach in the link you included. There is also very good software that will explicitly enumerate the extreme rays: PORTA and cdd/cdd+ are two that come to mind.

Answer 2

It is a difficult problem, however such spaces are characterized by duality theorems of alternative, in particular you could use the Motzkin theorem of the alternative to check if a given set of vectors form a convex basis.

Answer 3

There is not much information in your question to go on, at present, but certainly there is always the possibility of using any one of numerous methods to solve the linear system $Sx=0$ with the additional constraint imposed by your nonnegativity condition.

As Godot mentioned, this is not possible in general. Maybe there are some additional properties about stoichiometric matrices that would be good to know?

Answer 4

Hard way: Let the stoichiometry matrix (Michaelis-Menten reaction) be, $$A =\begin{bmatrix} -1& 1 & 0& 0\\ -1& 1 & 1& -1\\ 1 & -1 & -1 & 1\\ 0 & 0& 1 & -1 \end{bmatrix}$$ Step 1: find null space of the given matrix A ($null(A)$), $$N=\begin{bmatrix} 0 & 1\\ 1 & -1\\ 1 & 0\\ 0 & 1\end{bmatrix}$$ Step 2: add 2 or more columns of (null(A) and get a positive vector. For N, adding (in general, linear combination) two columns gives a vector $v=[1, 0, 1, 1]^T$. Then the vector v and the first column of N, will define the positive null space. i.e., $$N=\begin{bmatrix} 1 & 0\\ 0 & 1\\ 1 & 1\\ 1 & 0\end{bmatrix}$$

But sometimes it will be hard to find which columns to be added.

Easiest way: For a matrix 'A' which has nullspace N (in MATLAB), $$N=rref(null(A','r')')'$$ This will give you the null space after rref . Sometimes leads to a positive null space. If you are getting negative values, you can find the positive null space by doing step 2. But here finding columns to be added will be much easier.

$$N=rref(null(A','r')')'=\begin{bmatrix} 1 & 0\\ 0 & 1\\ 1 & 1\\ 1 & 0\end{bmatrix}$$

This is not the exact way. But an easy way.

Answer 5

You could make use of a projection onto two convex sets (POCS). Since $S$ is known, the pseudo-inverse $S^\dagger$ can be found. Now the projection onto the kernel of $S$ is defined as $P_{\text{ker}(S)}(x) := x-A^\dagger Ax$. The projection onto the positive set is very easy, you just replace all negative values with zero. By iterating, you converge to the true solution.