Exploring underdetermined linear system with non-negative solution I haven't had much luck searching for this specific problems. Any pointers would be greatly appreciated.
I have an underdetermined system where $ A $ and $ b $ are known. $ x $ is a real vector with $x_i > 0 $. 
$$ Ax = b $$
I'm not looking to optimize a function but to explore the solution space.
I'd like to be able to incrementally move around the solution space, using the "current" solution to evaluate a completely different problem (that would be absolutely horrible and infeasible to integrate into this). How would I be able to define this space?
Could anyone give an idea of where to start? Thanks!
 A: The solution space is an affine space and can be described as the vector space of the solutions of the corresponding homogeneous system $Ax=0$ shifted by a specific solution $x_0$ of the inhomogeneous system $Ax=b$. Thus, the general solution takes the form $x=x_0+Sv$, where $S$ is a matrix whose columns form a basis of the solution space of the homogeneous equation and $v$ is a column vector of coefficients. You can move around the solution space by varying $v$.
A: Say you have $x_1,\ldots,x_{n+1}$ with $x_k > 0$ and $Ax_k = b$ for all $1 \leq k \leq n+1$. Then every convex combination of these vectors, i.e. every sum $$
  y = \sum_{k=1}^{n+1} \lambda_kx_k
  \quad\text{where}\quad
  \sum_{i=1}^{n+1} \lambda_k = 1
  \quad\text{and for all $k$}\quad
  \lambda_k \geq 0
$$
also obeys $Ay=b$ and $y > 0$.
Now, the question becomes whether you can find "enough" such $x_i$ to guarantee that every solution $y$ can be obtained this way. The strictness in $x > 0$ causes trouble here, I think. Take for example the system $$
    \left(\begin{matrix}1&1\\0&0\end{matrix}\right)x
  = \left(\begin{matrix}1\\0\end{matrix}\right) \text{,}
$$ i.e. the system whose solution is a line through $(1,0)$ and $(0,1)$. Due to the strictness in your requirement, you have to $x_1=(1-\epsilon_1,\epsilon_1)$, $(\epsilon_2,1-\epsilon_2)$ for $\epsilon_1,\epsilon_2 > 0$ as your convex basis. But that means you'll miss e.g. the solution $(1-\frac{\epsilon_1}{2},\frac{\epsilon_1}{2})$.
So fix this, you'll have to replace $x > 0$ with $x \geq 0$ in your requirement, and then pick only convex combinations which satisfy the stricter condition. Usually simply requirement $\lambda_i > 0$ instead of $\lambda_i \geq 0$ should be enough, but there might cases where it's not. I have the feeling that in these cases, the stricter system has no solution at all, but I haven't proved this so I might be wrong.
You'll also have to figure out what "enough" means. $n$ can quite obviously be at most $\dim \ker A$, but there are cases where you won't find that many affinely independent vectors with only positive coordinates. Take for example the system $$
    \left(\begin{matrix}1&1&0&0\\0&0&1&1\\0&0&0&0\\0&0&0&0\end{matrix}\right)x
  = \left(\begin{matrix}1\\0\\0\\0\end{matrix}\right) \text{.}
$$
If the third or forth component of a solution is strictly greather than zero, the other one needs to be strictly less than zero, and you thus get only the two vectors $x_1=(1,0,0,0)$ and $x_2=(0,1,0,0)$.
I'd take a look at the simples algorithm, it should help to find suitable $x_k$. 
