Find positive solutions to underdetermined linear system of equations

Question

I'm a bit new to matlab so sorry if this is horribly simple.

Consider a problem of the following:

Find x_1, x_2, x_3 > 0 such that

67.5 = 60*x_1 +  90*x_2 + 120*x_3 and
60   = 30*x_1 + 120*x_2 +  90*x_3

In this case I want the solution

0 < x_3 < 3/7,
x_2 = 7/20 - 4/10*x_3, and
x_1 = 2/5  -  7/5*x_3

Is there a easy way to make Matlab solve such a problem for me?

Answer 1

The easy answer, since you just need non-negativity constraints on the parameters, is to use lsqnonneg. lsqlin is not needed at all for this problem and since lsqnonneg is part of MATLAB, the optimization TB is not needed.

A = [60 90 120; 30 120 90];
b = [67.5; 60];

x = lsqnonneg(A,b)
x =
                         0
         0.178571428571428
         0.428571428571429

We can test the result to see if it solved the original problem.

A*x - b
ans =
     0
     0

Of course you COULD have used lsqlin, but why bother?

In fact, there is an issue we must consider since an underdetermined system has infinitely many solutions. We can add any amount of the null-space of the array A to our solution. In this case that null-space has rank 1. It is characterized by the vector N here:

N = null(A)
N =
        -0.792593923901216
         -0.22645540682892
         0.566138517072299

An easy way to understand it is to recognize what the nullspace of A means. N is a vector (or a set of basis vectors that span a subspace in case the null-space has dimension higher than 1) such that

A*N = 0

Essentially, N is a vector that is orthogonal to all of the rows of A. If the null-space has dimension higher than 1, then N can be any linear combination of basis vectors of the null-space. Therefore, if X is ANY solution of the under-detemined system

A*x = b

then it must be true that x+c*N is also a solution, for any arbitrary constant c. (Remember that A*N will be zero.)

So for example, I'll pick an arbitrary coefficient for N:

x2 = x + N*(-.1)
x2 =
        0.0792593923901216
          0.20121696925432
         0.371957576864199

Again, x2 is also a solution. It too has entirely positive coefficients. (You can trivially find that range of values for the coefficient on N such that the solution is entirely positive.)

A*x2 - b
ans =
      -1.4210854715202e-14
       -7.105427357601e-15

(Note that these are effectively zeros, to within the double precision trash found in floating point arithmetic.)

So, IF we wanted to do so, it is easy enough to start with the lsqnonneg or backslash or pinv solution, and find the complete set of solutions to your system such that the coefficients are entirely positive. Hint: all you need do is consider the vectors x and N, then look for solutions to the problem

(x + c*N) > 0

where c is some scalar constant. Clearly, C cannot be positive, else the first element of the sum will be negative.

C = -x./N
C =
            0
      0.78855
     -0.75701

x + C(1)*N
ans =
            0
      0.17857
      0.42857

x + C(3)*N
ans =
          0.6
         0.35
            0

As we can see, when c is any value in the closed interval [-.75701,0], we get an entirely positive solution to the problem, in the form of (x+c*N). Go any further than those limits though, and one or more of the elements in the solution will be negative.

On some problems there will be no solution at all that yields an exact solution with all positive elements of the solution vector. This is entirely possible. For example, suppose we changed the original problem into:

A = [60 90 120; 30 120 90];
b = [-67.5; -60];

Now what happens when we apply lsqnonneg?

lsqnonneg(A,b)
ans =
     0
     0
     0

An all-zero solution results. Since that solution clearly can not exactly solve the original problem, it must be true that no such positive solution exists.