How can i solve this linear system ?

Hello,

I have a linear system of equation composed of 14 equations and 49 unknowns. I know that my system can be writed in the Ax=b form and then performing the division A\b to find the unknown vector x. I'm new to Matlab so sorry if I'm asking trivial questions.

My questions are the following:

1) My system should have multiple solutions since it is underdetermined. A\b only gives 1 solution. Is there any way to get more than 1 solution ?

2) I already know that all my 49 unknown are within the range [0;1]. How can I set this constrains to all my unkown ?

3) I may have some unknown that I already know before solving the system. Is it possible to restrict some value of the x vector without having to rewrite the whole system like it is said in the following thread: http://www.mathworks.fr/matlabcentral/answers/17795-basic-matrix-equation-ax-b-with-restrictions

4) Is it possible to import a matrix from a text file, because my A matrix is quite big ?

Regards, Debzz

Answer 1

[John D'Errico]

Your problem is such that there are infinitely many solutions. backslash will return ONE solution, but you want more. How many more? 2? 37? 454233531211 of them?

The general solution can be written in the form

X = X0 + null(A)*k

where k is any general vector, here of size [49-14,1], and X0 is given by

X0 = A\b

If you further require the solution to fall in the 49-dimensional unit box [0,1]^49, it is entirely possible that no solution exists at all. BUT, lsqlin can solve for A solution if any exists. Then you use the same trick that I showed above to find a new solution.

If you don't have the optimization toolbox, so no lsqlin, then you can use lsqnonneg, employing slack variables to transform the problem into one that it can solve. This will require some minor coding of course.

If some of the parameters are already known, then you could simply eliminate them from the problem. Subtract those terms from the right hand side of the linear system. Or, again, if you have lsqlin, then simply set up equality constraints. (This is a far better solution than setting the lower and upper bounds to the same value, which will probably introduce numerical problem.)

And as far as reading in a text matrix, of course MATLAB can read in a file. 14x49 is not even remotely what anyone would call big though.

Answer 2

[debzz]

Thanks for your help.

After some readings about lsqlin I came up with this code:

MATLAB code
C=[1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;

    0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1;
    9957.81 0 0 0 0 0 0 1168.06 0 0 0 0 0 0 4817.60 0 0 0 0 0 0 5417.43 0 0 0 0 0 0 520.87 0 0 0 0 0 0 86.75 0 0 0 0 0 0 210.55 0 0 0 0 0 0;
    0 9957.81 0 0 0 0 0 0 1168.06 0 0 0 0 0 0 4817.60 0 0 0 0 0 0 5417.43 0 0 0 0 0 0 520.87 0 0 0 0 0 0 86.75 0 0 0 0 0 0 210.55 0 0 0 0 0;
    0 0 9957.81 0 0 0 0 0 0 1168.06 0 0 0 0 0 0 4817.60 0 0 0 0 0 0 5417.43 0 0 0 0 0 0 520.87 0 0 0 0 0 0 86.75 0 0 0 0 0 0 210.55 0 0 0 0;
    0 0 0 9957.81 0 0 0 0 0 0 1168.06 0 0 0 0 0 0 4817.60 0 0 0 0 0 0 5417.43 0 0 0 0 0 0 520.87 0 0 0 0 0 0 86.75 0 0 0 0 0 0 210.55 0 0 0;
    0 0 0 0 9957.81 0 0 0 0 0 0 1168.06 0 0 0 0 0 0 4817.60 0 0 0 0 0 0 5417.43 0 0 0 0 0 0 520.87 0 0 0 0 0 0 86.75 0 0 0 0 0 0 210.55 0 0;
    0 0 0 0 0 9957.81 0 0 0 0 0 0 1168.06 0 0 0 0 0 0 4817.60 0 0 0 0 0 0 5417.43 0 0 0 0 0 0 520.87 0 0 0 0 0 0 86.75 0 0 0 0 0 0 210.55 0;
    0 0 0 0 0 0 9957.81 0 0 0 0 0 0 1168.06 0 0 0 0 0 0 4817.60 0 0 0 0 0 0 5417.43 0 0 0 0 0 0 520.87 0 0 0 0 0 0 86.75 0 0 0 0 0 0 210.55];
d=[1 1 1 1 1 1 1 11613.45 399.11 3475.97 5311.40 1139.73 87.63 151.79];
lb=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
ub=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
A=[];
b=[];
Aeq=[];
beq=[];
options = optimoptions('lsqlin','Algorithm','active-set','MaxIter',1000000);
x=lsqlin(C,d,A,b,Aeq,beq,lb,ub,x,options);

Matlab keep telling me that the number of iterations exceeded. I increased MaxIter parameter to 1 million iterations but the number of iterations keep being exceeded. I'm I doing something wrong ?

[John D'Errico]

By the way, you don't need to add answers. Use comments instead. An answer should be exactly that, an answer.

[John D'Errico]

I looked at your problem as posed here. First, the fact that C and d do not conform properly for linear algebra, in the sense that d should be 14x1, not 1x14. So...

d = d.';

Next, your problem is probably poorly scaled. Large disparities in the magnitude of those matrix elements will potentially cause numerical problems with many optimization schemes. So...

r = max(abs(C),[],2);
Chat = bsxfun(@times,C,1./r);
dhat = d'./r;

This may have helped a bit. See that Chat is a bit better behaved than was C. But as importantly, see that C was singular! It has a zero singular value there, so the rank of your matrix was only 13, not 14. RANK would have told us the same thing, but SVD tells us a bit more. Note that a singular matrix tells us that an exact solution need not exist.

[svd(C),svd(Chat)]
ans =
      12386       2.9235
      12386       2.6458
      12386       2.6458
      12386       2.6458
      12386       2.6458
      12386       2.6458
      12386       2.6458
     2.6458       1.2438
     2.6458       1.2438
     2.6458       1.2438
     2.6458       1.2438
     2.6458       1.2438
     2.6458       1.2438
 1.4763e-13   2.7105e-16

Does any exact solution exist? If we augment Chat with dhat, then test the rank, if that rank were to change, then in fact there is no exact solution, even disregarding the bound constraints.

svd([Chat,dhat])
ans =
     4.0219
     2.7745
     2.6458
     2.6458
     2.6458
     2.6458
     2.6458
      1.415
     1.2438
     1.2438
     1.2438
     1.2438
     1.2438
  2.093e-07

Ah, so the augmented matrix is no longer rank deficient. So no exact solution exists.

Next, before I ever try lsqlin, use backslash. What does it do? (Besides tell us that C was singular, which we already knew.)

Chat\dhat
Warning: Rank deficient, rank = 13, tol =  1.538691e-14. 
ans =
     1.1663
    0.04008
    0.34907
  -0.010648
    0.11446
   -0.67447
   0.015243
          0
          0
          0
          0
          0
          1
          0
          0
          0
          0
          0
          0
          1
          0
          0
          0
          0
          1
          0
          0
          0
          0
          0
          0
          0
          0
          1
          0
          0
          0
          0
          0
          0
          1
          0
          0
          0
          0
          0
          0
          1
          0

BACKSLASH gave a valid answer, although one that fails your constraints.

Remember, that since we already knew there was no exact solution, we cannot expect perfection. Try lsqlin now. I'll be lazy with the options, using the defaults.

xhat = lsqlin(Chat,dhat,[],[],[],[],lb,ub)
Warning: The trust-region-reflective algorithm requires at least as many equations as variables
    in C matrix; using active-set algorithm instead. 
> In lsqlin at 301 
Maximum number of iterations exceeded; increase options.MaxIter.
To continue solving the problem with the current solution as the
starting point, set x0 = x before calling lsqlin.
xhat =
      0.12429
            0
      0.40281
  -2.1684e-19
      0.19805
     0.092399
     0.098842
            1
  -5.0388e-17
            0
            0
            0
  -5.1065e-12
   -1.548e-16
      0.95955
            0
            0
  -1.1297e-16
            0
            0
            0
            1
  -3.4694e-18
            0
            1
  -1.7347e-18
            0
            0
            0
            1
            0
            0
            0
   1.1102e-16
            0
            0
   1.7347e-18
            1
            0
            0
            0
            0
            0
            0
            1
            0
            0
            0
            0

Again, it found a solution that satisfies the bounds. But NO exact solution exists, so I doubt that we can do much better here. Allowing it to use an infinite number of iterations will not help. My guess is the singularity of your system is causing LSQLIN to keep looking.

Answer 3

[debzz]

Ok, so just to explain better my problem I did some test. I took the problem backwards. I had a similar linear system (Cx=d) to the one I posted above. In this one I already know a solution ( S ) so I wanted to see if lsqlin could solve x. To see if the solution was indeed a solution, I computed A*S and compare the result with d.

Matlab code
A*S=[1 1 1 1 1 1 1 9939.3 1169.5 4820.8 5396.5 521.0 84.4 209.6]
d=[1 1 1 1 1 1 1 9939.44 1169.60 4821.32 5396.32 520.68 84.52 210.32]

As you can see A*S and d are pretty close so the solution S is a good solution and now we are pretty sure that a solution exist.

Then I tried to use lsqlin to find a solution. I used the same script as I posted above (but with a different C and d matrix of course).

lsqlin can't find any solutions. And no matter how many iterations I set (I increased MaxIter parameter up to 1 million) the number of iteration is always exceeded and the x vector doesn't seem to improve at all.

In addition, I set the x0 initial solution vector as the vector S (the solution that I already knew). I was expecting the solver to do only 1 or 2 iterations since I fed it with the solution. Unfortunatly it didn't work. In fact the solver still can't find a solution eventhough I gave it to him.

Finally I swithed the algorithm to 'active-set' because when I use the default one I get the following warning:

MATLAB Code
Warning: Large-scale algorithm requires at least as many equations as variables
  in C matrix; using medium-scale algorithm instead.

One more thing: I try to find a solution without the use of lower and upper bound and it worked. The problem is that the solution is useless because some of the solution are not within the range of [0;1].

There must be something that I am missing but I can't see it. Do you have a hint or something ?

[John D'Errico]

You say that you have a problem that DOES have a known solution, but you don't show the problem. How can I add anything?

Read my response to your LAST answer, as a rank deficient system probably causes problems.

And please STOP ADDING an infinite number of answers that are not answers! This only confuses things. USE COMMENTS!

[debzz]

Sorry, my bad. Thank you very much for your answer. It is very helpful and clarifying.

Here is the C and d Matrix (Cx=d)

C=[
  1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1;
  8724.80 0 0 0 0 0 0 1549.60 0 0 0 0 0 0 5453.64 0 0 0 0 0 0 5643.12 0 0 0 0 0 0 394.52 0 0 0 0 0 0 82.88 0 0 0 0 0 0 293.64 0 0 0 0 0 0;
  0 8724.80 0 0 0 0 0 0 1549.60 0 0 0 0 0 0 5453.64 0 0 0 0 0 0 5643.12 0 0 0 0 0 0 394.52 0 0 0 0 0 0 82.88 0 0 0 0 0 0 293.64 0 0 0 0 0;
  0 0 8724.80 0 0 0 0 0 0 1549.60 0 0 0 0 0 0 5453.64 0 0 0 0 0 0 5643.12 0 0 0 0 0 0 394.52 0 0 0 0 0 0 82.88 0 0 0 0 0 0 293.64 0 0 0 0;
  0 0 0 8724.80 0 0 0 0 0 0 1549.60 0 0 0 0 0 0 5453.64 0 0 0 0 0 0 5643.12 0 0 0 0 0 0 394.52 0 0 0 0 0 0 82.88 0 0 0 0 0 0 293.64 0 0 0;
  0 0 0 0 8724.80 0 0 0 0 0 0 1549.60 0 0 0 0 0 0 5453.64 0 0 0 0 0 0 5643.12 0 0 0 0 0 0 394.52 0 0 0 0 0 0 82.88 0 0 0 0 0 0 293.64 0 0;
  0 0 0 0 0 8724.80 0 0 0 0 0 0 1549.60 0 0 0 0 0 0 5453.64 0 0 0 0 0 0 5643.12 0 0 0 0 0 0 394.52 0 0 0 0 0 0 82.88 0 0 0 0 0 0 293.64 0;
  0 0 0 0 0 0 8724.80 0 0 0 0 0 0 1549.60 0 0 0 0 0 0 5453.64 0 0 0 0 0 0 5643.12 0 0 0 0 0 0 394.52 0 0 0 0 0 0 82.88 0 0 0 0 0 0 293.64];

d=[1 1 1 1 1 1 1 9939.44 1169.60 4821.32 5396.32 520.68 84.52 210.32];

And here is a known solution to the problem. Eventhough it's not exact, it's pretty close as I have shown above.

Sol=[0.8837; 0.0112; 0.0494; 0.0411; 0.0095; 0.0015; 0.0036; 0.1334; 0.6379; 0.1679; 0.0022; 0.0567; 0.0016; 0.0003; 0.2357; 0.0095; 0.6938; 0.0538; 0.0056; 0.0005; 0.0010; 0.1027; 0.0006; 0.0566; 0.8376; 0.0007; 0.0000; 0.0017; 0.1206; 0.0688; 0.0496; 0.0027; 0.7462; 0.0108; 0.0012; 0.2085; 0.0072; 0.0106; 0.0188; 0.0376; 0.4778; 0.2394; 0.3154; 0.0011; 0.0206; 0.0403; 0.0620; 0.0757; 0.4849];

I have used your code to try to find a solution on both linear systems and like you said the results are quite bad because after running lsqlin I computed C*x and it is different from d. Do you think that restricting a little bit more the solution with lb and ub vector could improve things.

In fact the x vector represent quantity of change of land cover. For exemple the first value represent the quantity of forest that is going to persist and remain forest in the future. The second value of the x vector represent the quantity of grasslands that is going to transform into forest and so on.

At first I just wanted to find a solution without putting restriction (except [0;1] because probabilities have to be within that range).

Since it doesn't work quite well I'm probably going to try restricting the solution with equalities and inequalities. Do you think putting some equality restriction on some values could improve the result given by lsqlin ?

[John D'Errico]

Again, C is a rank deficient matrix, so for a solution to exist, we need to have d be representable as a linear combination of the columns of C. We can test if that is true by a simple test:

rank([C,d'])
ans =
    13

rank(C)
ans =
    13

Both results must yield the same result, else an exact solution will not exist. The solution you provided is not really the exact solution of course, probably because of rounding issues.

[C*Sol,d',C*Sol - d']
ans =
          1            1            0
          1            1            0
     0.9999            1      -0.0001
     0.9999            1      -0.0001
     0.9999            1      -0.0001
     0.9999            1      -0.0001
          1            1            0
     9939.3       9939.4     -0.17258
     1169.5       1169.6     -0.13423
     4820.8       4821.3     -0.50531
     5396.5       5396.3      0.21858
     520.95       520.68      0.27128
     84.383        84.52     -0.13719
     209.62       210.32     -0.69796

It did not do too badly, but lsqlin will be trying to chase down an exact solution, and we have a rank deficient problem.

Note that if I drop the last row of C, it does not change the rank of C. There are ways of finding which row might be the best one to drop.

rank(C(1:13,:))
ans =
    13

A simple one is...

[Q,R,E] = qr(C',0);
E
E =
   8     9    10    11    12    13    14     6     7     5     2     3     4     1

That last element of E is the row we could have dropped. Here E(end) was 1, which tells me that the first row is mathematically the logical choice to drop.

How does LSQLIN do here?

options = optimset('lsqlin');
options.Algorithm = 'active-set';

x = lsqlin(C(2:14,:),d(2:14)',[],[],[],[],lb,ub,[],options)
Optimization terminated.

x =
  -5.5511e-17
            0
       0.3815
       0.6185
            0
            0
            0
      0.52827
            0
            0
            0
      0.33601
            0
      0.13573
      0.49631
      0.21446
      0.27373
            0
            0
     0.015498
            0
            1
            0
            0
            0
            0
            0
            0
            1
            0
            0
  -6.6174e-24
            0
            0
            0
            1
            0
            0
            0
            0
            0
            0
            1
            0
            0
            0
            0
            0
   4.2352e-22

[C*x,d',C*x - d']
ans =
            1            1  -3.8636e-14
            1            1    2.176e-13
            1            1            0
            1            1            0
            1            1            0
            1            1            0
            1            1            0
       9939.4       9939.4            0
       1169.6       1169.6            0
       4821.3       4821.3            0
       5396.3       5396.3            0
       520.68       520.68            0
        84.52        84.52            0
       210.32       210.32            0

The result is now pretty much exact, with no warning messages generated, and no problem with the number of iterations needed. That the answer is not the same as Sol is irrelevant, since there are infinitely many valid solutions, and Sol was not truly an exact solution.

So clearly the problem LSQLIN (specifically the active-set algorithm) is having is with your rank deficiency. The fact that the system is underdetermined is not a problem.

Of course, we were able to happily drop one row of the matrix because we knew that d was representable as an exact linear combination of the columns of C. And valid solution to the last 13 rows must be as equally valid for the first row.

[debzz]

I have tried you code with the other linear system and it worked very well to. I'm pretty satisfied with this result right now but I might play and tweak the equalities, inequalities and bounds to see if I can have more control over the solution.

Thank you very much for your help and responsiveness M. John D'Errico.

[John D'Errico]

You have many dimensions to play with. A minor problem is that this gives you little control over which solution is chosen, but unless you know the answer, or you have more information than you supply to the system, then LSQLIN will be pretty arbitrary about what it gives.

Since I don't know what is good or bad here about the solution, I can't offer advice. However an alternative option you might try is to use either LSQLIN or QUADPROG to find that solution which is closest to 0.5*ones(49,1), whiie still satisfying the linear system and the bound constraints.