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.)
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.
Any elaboration on the LP approach or any other way of solving this is greatly appreciated.
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.