Finding nonnegative solutions to an underdetermined linear system

Question

Here's the environment of my problem: I have a linear system of 4 equations in 8 unknowns (i.e. $Ax = b$, where $A$ is $4 \times 8$, $x$ is $8 \times 1$, and $b$ is $4 \times 1$, with $A$ given and $b$ varying with some exogenous parameters, i.e. essentially given). So, after reducing the system, the solution space can be described by a linear function $f \, : \, \mathbb{R}^4 \rightarrow \mathbb{R}^4$, which I can find explicitly.

Because of the nature of the problem, I am interested only in nonnegative values for each of the 8 unknowns. So, I can restrict the function $f$ to the domain $\mathbb{R}^4_{\geq 0}$ (i.e. examine only nonnegative values for each of the 4 free variables). What I'm interested in knowing (and classifying for different parameter values, which will change $b$) is whether solutions exist for which all 8 unknowns are nonnegative--or equivalently, whether the function $f$ takes nonnegative values in the nonnegative domain $\mathbb{R}^4_{\geq 0}$.

If there's a way to figure this out analytically, that would be great; alternatively, if there's a good way to do this numerically (I use MATLAB), that would be only slightly less great. Thanks in advance!

---

Thank you all so much for your answers--this has been a huge help. Working out the solution by hand using the LL Dines paper may prove to be too cumbersome, as it requires knowing the sign of the coefficients, and we examine thousands of different sets of coefficients (as we vary several parameters); however, if nothing else, it's really nice to know that there is an algorithm for (analytically) determining the existence of nonnegative solutions. Moreover, all of the software advice given has been extremely helpful. I've looked into all of your suggestions, and for now I think we'll try using lsqnonneg. Thanks again for your help!

Answer 1

Matlab has a function lsqnonneg which can minimize $\| Ax - b\|_2$ subject to the constraint that $x \geq 0$.

The minimum value will be $0$ if and only if there exists a nonnegative solution to $Ax = b$.

Answer 2

This is a standard linear optimization problem. To summarize, the set of all nonnegative solutions can be completely described by either a finite system of linear inequalities, or a finite "generating set" such that every solution is a positive linear combination of elements in the generating set.

A $8\times 4$ matrix is very small and reasonable in this context, so if you have exact coefficients for your matrix you will be able to determine the set of nonnegative solutions exactly.

http://en.wikipedia.org/wiki/Linear_programming#Solvers_and_scripting_.28programming.29_languages

Answer 3

Fourier-Motzkin elimination may be of some help. You can write your equations as opposed pair of inequalities [$a'x = b$ becomes $a'x \ge b$ and $a'x \le b$].

Computer algebra software may be more useful than MATLAB. Here is the Mathematica documentation. The free Wolfram Alpha site can also solve systems of inequalities.

Answer 4

An interesting article on the positivity of the solutions of a linear system can be found here: http://www.jstor.org/stable/1968384?seq=1 (if I remember correctly, if you register, you have the right to read 3 articles for free, so you don't need to pay to see the article) An algorithm that can be done by hand is presented.

A Matlab toolbox which can solve optimization problems with constraints is Yalmip. It is free, easy to install (just copy the folder in the Matlab path) and looking a bit in the tutorial will show you that handling the constraints is really easy. You do not have to worry how you could write them in the matrix form. The software does that for you.

Taking your example, you have a $4 \times 8$ matrix $A$ which is given and a $8\times 1$ vector $b$ also given. The unknown is the vector $x$ with positive components.

A possible Matlab code (look up the syntax yourself on the documentation site):

A = [...]; b = [...]; %variable x = sdpvar(8,1); %constraints F = [A*x == b, x(:)>=0]; %find a feasible x solvesdp(F); %display double(x)

Answer 5

Perhaps not useful to the original poster at this point, but as I came across this question, the following was useful to me. The set $\{x \mid Ax = b, x \geq 0 \}$ is a polyhedron (as mentioned in another answer, it can also be represented by a linear combination of generators, if needed). There are libraries that can handle polyhedra efficiently, e.g., PPL, Apron. Some parts of PPL are accessible from Sage. These libraries provide API to build polyhedra from constraints (or generators), to check for emptiness, to project (will be useful for the OP's question "for which values of $b$ the polyhedron is non-empty") and so on.