How to get all the solutions to an underdetermined system of linear equations using lsqnonneg()?

MUKTIKANTA (view profile)

on 13 Oct 2018
Latest activity Commented on by MUKTIKANTA

on 14 Oct 2018

Bruno Luong (view profile)

I am trying find out all the positive solutions to an underdetermined system of linear equations and my code is as follows:
format rat
a=[3 22 119 620 3 1 3 3 2 20 20 22 22 119 119 3 3 3 3 3 3 1 1 1 20 20 20 20 2 1 2 20 20 3 3 3 3 3 3 1 1 1 2 2 1;
3 0 0 0 3 1 3 3 2 0 0 0 0 0 0 3 3 3 3 3 3 1 1 1 0 0 0 0 2 1 2 0 0 3 3 3 3 3 3 1 1 1 2 2 1;
0 22 0 0 3 1 0 0 0 20 20 22 22 0 0 3 3 3 3 0 0 1 1 1 20 20 20 20 2 1 0 20 20 3 3 3 3 3 3 1 1 1 2 2 1;
0 0 119 0 0 0 3 0 0 20 20 0 0 119 119 3 3 0 0 3 3 1 0 0 20 20 20 20 0 0 2 20 20 3 3 3 3 3 3 1 1 1 2 2 1;
0 0 0 620 0 0 0 3 2 0 0 22 22 119 119 0 0 3 3 3 3 0 1 1 20 20 20 20 2 1 2 20 20 3 3 3 3 3 3 1 1 1 2 2 1];
b=[620;3;22;119;620];
w=lsqnonneg(a,b);
Each time I run this code, it gives me only one solution. How can I get all the solution? Is there any method available in MATLAB to achieve this?

Walter Roberson

Walter Roberson (view profile)

on 13 Oct 2018
That is rather the wrong function to use.
An underdetermined system of linear equations describes hyperplane (or possibly a plane or a line.)
You have 45 variables and 5 equations. Given any arbitrary values for the last 40 variables, exact values of the first 5 can be calculated.
The set of all positive solutions could only be described by boundary equations. I would have to work out how many of those you would need. It is late at night here and one avenue of thought is suggesting that you might need 45!/5! so I would want to think about that when I am more awake. Maybe it is 45!/40! ?

Bruno Luong (view profile)

on 13 Oct 2018
Edited by Bruno Luong

Bruno Luong (view profile)

on 13 Oct 2018

The number of solutions is infinity. Explanation:
>> size(a)
ans =
5 45
In your case you have 5 eqt and 45 unknowns, and
>> rank(a)
ans =
5
confirms your matrix is full row-rank. meaning all the 5 rows of A are independent vectors in R^45.
One particular solution can be obtained with backshlash (or lsqnonneg if you prefer)
>> x0=a\b;
>> fprintf('x0 = %s\n', mat2str(x0));
x0 = [0;0;5.93708057585062e-17;0.77741935483871;0;0;0;0;0;0;0;0.863636363636363;0;0.974789915966387;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0.999999999999998;0;0;0;0;0;0;0;0;0;0;0]
The kernel of a, meaning a basis of equation A*X = 0 can be obtained with NULL command.
>> X=null(a);
>> size(X)
ans =
45 40
It tells that the kernel has dimension 40 and a maths theorem tells it's equals to size(a,2)-rank(a).
To confirm X is the kernel, you can check with frobenus norm of max values (they are 0 at the numerical precision)
>> y=a*X;
>> norm(y,'fro')
ans =
1.5079e-13
>> max(max(abs(y)))
ans =
7.5533e-14
Now all solutions of a*x=b can be written as
x = x0 + X*r
where r is any arbitrary vector of size 40 x 1. The set of solutions is an affine space of dimension 40 and there are of course an infinity vectors in this huge space. To check it try:
r = rand(40,1); % any random vector
x = x0 + X*r; % general form of the solution a*x=b
check if a*x = b
>> max(abs(a*x-b))
ans =
1.1369e-13
If you add the constraints (x >= 0), meaning
x(i) >= 0 for all i = 1,...45
Then the general solution is the intersection of 45 half-hyperspace and the affine space of dimension 40.
It will be a polytope of dimension 40 living in R^45. It is still huge. This polytope is actually a an interior of the convex combination of n vertices, n can be as high as:
>> nchoosek(45,40)
ans =
1221759
each of them is one particular solution of
a*x = b
x >= 0
There is no way you can enumerate them all and make on top of that an arbitrary sub-convex combination if that is your intention.

Walter Roberson

Walter Roberson (view profile)

on 13 Oct 2018
Thanks, Bruno. I would have found it tough to pull all those steps together.
Bruno Luong

Bruno Luong (view profile)

on 13 Oct 2018
Note: what I wrote is not entirely correct, sometime the polytope is unbounded if no other constraint on x is introduced and it can not express as the convex combination of finite set of vertexes.
However I hop it discourage people trying to beat the market with such brute force method.
MUKTIKANTA

MUKTIKANTA (view profile)

on 14 Oct 2018
Thanks, Bruno for such a well-explained answer.