solving system of equations using singular value decomposition and fsolve

Hello, i have a system of equations which i want to solve. My simplified system is shown below.

k1=0.6828;
k2=1.3656;
k3=2.0485;
k4=2.7313;
l=0.1920;

M=        [cos(k1*l)^2,          (sin(k1*l)*cos(k1*l))/k1,       sin(k1*l)*cos(k1*l),         (sin(k1*l)^2)/k1;         
          -sin(k1*l)*cos(k1*l),         -(sin(k1*l)^2)/k1,             cos(k1*l)^2,         (sin(k1*l)*cos(k1*l))/k1;
          cos(k2*l)^2,           (sin(k2*l)*cos(k2*l))/k2,       sin(k2*l)*cos(k2*l),         (sin(k2*l)^2)/k2;         
          -sin(k2*l)*cos(k2*l),         -(sin(k2*l)^2)/k2,             cos(k2*l)^2,         (sin(k2*l)*cos(k2*l))/k2;
          cos(k3*l)^2,           (sin(k3*l)*cos(k3*l))/k3,       sin(k3*l)*cos(k3*l),         (sin(k3*l)^2)/k3;         
          -sin(k3*l)*cos(k3*l),         -(sin(k3*l)^2)/k3,             cos(k3*l)^2,         (sin(k3*l)*cos(k3*l))/k3;
          cos(k4*l)^2,           (sin(k4*l)*cos(k4*l))/k4,       sin(k4*l)*cos(k4*l),         (sin(k4*l)^2)/k4;         
          -sin(k4*l)*cos(k4*l),         -(sin(k4*l)^2)/k4,             cos(k4*l)^2,         (sin(k4*l)*cos(k4*l))/k4];

F= [0.0051; -0.0092; 0.0199; -0.0157; 0.0436; -0.0168; 0.0745; -0.0103];
X= [xi; xpi; yi; ypi];
% F=M*X is the equation to solve with elements of X being unknowns

i have a wide range of k values to solve (and more complex M also but for simplicity i stuck to this case). I have tried to solve using Singular value decomposition (SVD), fsolve and but all these functions reveals different results. Please advice what will be the best way to solve such system. SVD was selected because in some cases, matrix m became rank deficit. My solutions are as follows.

[U,S,V] = svd(M,0)
U1=U(:,1:2);
S1=S(1:2,1:2);
V1=V(:,1:2);
c = U1.'*(F);
p1 = inv(S1)*c;
X = V1*p1

I also tried using "solve" or "fsolve" but all of these functions give me different values of X. Please advice on how to solve such a system.

Answer 1

[John D'Errico]

Solving a singular (or nearly so) linear system of equations will have infinitely many solutions, all equally good (or bad).

Fsolve is not a good way to solve the problem, using an iterative solver to find a solution will at best only provide a different solution for every set of starting values, all equally bad.

Use pinv, or lsqminnorm.

X = pinv(M)*F;

or (in a recent MATLAB release)

X = lsqminnorm(M,F);

If you have additional information about the variables, then there are ways you can bring that into the problem. But without that, there is nothing you can simply do that is better than the above solutions.

[Sumera Yamin]

thank you for your detail reply. So now lets say i modify my system, and make it a determined one, with 4 equations and 4 unknowns, Then can i still use pinv(M)*F safely? For a determined system, what should be the best way to solve it?

k1=0.6828;
k2=1.3656;
l=0.1920;
M=        [cos(k1*l)^2,          (sin(k1*l)*cos(k1*l))/k1,       sin(k1*l)*cos(k1*l),         (sin(k1*l)^2)/k1;         
          -sin(k1*l)*cos(k1*l),         -(sin(k1*l)^2)/k1,             cos(k1*l)^2,         (sin(k1*l)*cos(k1*l))/k1;
          cos(k2*l)^2,           (sin(k2*l)*cos(k2*l))/k2,       sin(k2*l)*cos(k2*l),         (sin(k2*l)^2)/k2;         
          -sin(k2*l)*cos(k2*l),         -(sin(k2*l)^2)/k2,             cos(k2*l)^2,         (sin(k2*l)*cos(k2*l))/k2];
F= [0.0051; -0.0092; 0.0199; -0.0157];
X= [xi; xpi; yi; ypi];
% F=M*X is the equation to solve with elements of X being unknowns

[John D'Errico]

Yes. The beauty of pinv as a solution, is it works as well as possible if the system is singular or nearly so, or if it is perfectly well conditioned.

If your matrix is non-singular, then X=pinv(M)*F will yield the same solution as does M\F.

You can absolutely be safe and just use pinv always there. The only downside to using pinv is it will be slightly less efficient. Thus pinv takes just a bit more work. But that difference is TINY for a 4x4 problem. For example:

M = rand(4);
F = rand(4);

timeit(@() M\F)
ans =
    6.497e-06

timeit(@() pinv(M)*F)
ans =
   6.4058e-05

So, it took roughly 10 times as long. But in terms of absolute numbers, who cares if you take an extra tiny fraction of a millisecond? :)

If you are doing it many millions or billions of times, well, yes, it may be significant. But even for a few thousand times, or tens or even hundreds of thousands, it simply is not worth worrying about the time differential there.

Answer 2

[Anton Semechko]

This is an overdetermined system (i.e., number of equations > number of variables). You can obtain least squares solution as X=(M'*M)\(M'*F). Note that (M'*M) is the matrix that is being inverted here, not M (because M is not a square matrix).

[John D'Errico]

No. Don't tell people to use that form to solve a system of equations. It is a poor solution. MATLAB provides this:

X = M\F;

which works even for over-determined systems.

If M is nearly singular, then use of the equation form you suggest merely makes the system MORE singular, squaring the condition number.

Use pinv instead, which uses the SVD.

X = pinv(M)*F

Or use lsqminnorm.

[Anton Semechko]

If you want to use SVD to solve (M'*M)\(M'*F), you can do so as follows:

[U,D]=svd(M'*M); % note M'*M=U*D*U'
X=(U*diag(1./diag(D))*U')*M'*F; % here, I made the assumption that all diagonal entries of D are non-zero

This will yield the same answer as M\F, which automatically recognizes that the system is overdetermined.

[Anton Semechko]

In case, M'*M is singular (i.e., doesn't have an inverse), you can obtained regularized solution as follows:

[U,D]=svd(M'*M);
d=diag(D);
d_inv=d;
d_inv(d>eps)=1./d_inv(d>eps);
X=(U*diag(d_inv)*U')*M'*F;

[Anton Semechko]

John D'Erico, in an overdetermined system, the matrix M will have more rows than columns, that is why I said the system must be solved using least squares. In this specific situation, the matrix that must be inverted is M'*M not M^2. When M has more rows than columns, it is not even possible to compute M^2, so I am not sure how your comment is relevant here.

[John D'Errico]

Because you don't understand the linear algebra. Nor do you apparently understand my comment at all. You are throwing around ideas like the SVD, but you don't seem to understand the linear algebra.

In an overdetermined system, I NEVER said you compute M^2. However, M'*M does indeed square the condition number of the problem.

Apparently it is time for a (sadly, long) lesson on some linear algebra.

Consider a non-square, but nearly singular matrix M. Do you know about the condition number of a matrix?

x = rand(20,1);
M = x.^(0:10);
b = rand(20,1);

cond(M)
ans =
   1.6808e+08

cond(M'*M)
ans =
   2.9607e+16

cond(M).^2
ans =
   2.8251e+16

What you don't want to do is compute M'*M. That basically throws away information, leaving much of the information in that matrix now lost for double precision computations. You don't want to tell people to compute M'*M, then compute the svd of that product. Again, that is a bad thing!

The point is, you are telling people to use the normal equations, thus (M'*M)\(M'*b), when there is no need to do so, and when that can be an actively bad thing to do.

MATLAB knows how to compute the solution of a NON-square poorly conditioned system. Those methods are already built into backslash, which applies at least for non-singular matrices M. backslash uses a pivoted QR decomposition. We can get that as:

[Q,R,E] = qr(M);

This is a decomposition with the form

M*E = Q*R

So if M*x=b, then we can write it as

Q*R*E'*x = b

E is just a permutation matrix, so effectively permuting the elements of x.

We can get x most simply as

x = E*(R\(Q'*b));

As you see here, we never needed to square the condition number of M, because we never needed to compute M'*M. The only equations that needed to be solved at all were in the back substitution with the matrix R. And R has the same condition number as M.

cond(R)
ans =
   1.6808e+08

So, how did that work?

x = E*(R\(Q'*b))
x =
       48.647
        -1705
        24325
  -1.8643e+05
   8.5977e+05
  -2.5191e+06
   4.7907e+06
  -5.8833e+06
   4.4945e+06
  -1.9398e+06
   3.6096e+05

M\b
ans =
       48.647
        -1705
        24325
  -1.8643e+05
   8.5977e+05
  -2.5191e+06
   4.7907e+06
  -5.8833e+06
   4.4945e+06
  -1.9398e+06
   3.6096e+05

(M'*M)\(M'*b)
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.371326e-17. 

ans =
         50.8
        -1788
        25636
  -1.9763e+05
   9.1758e+05
  -2.7088e+06
   5.1948e+06
  -6.4393e+06
   4.9698e+06
   -2.169e+06
   4.0855e+05

As you can see, the QR computation gave the same results here as backslash, with M\b. Whereas the normal equations form you seem to be wedded to introduces junk into the solution. In fact, does the form you seem to love produce a minimum norm solution at all?

x = M\b;
xn = (M'*M)\(M'*b);
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.371326e-17. 

norm(M*x - b)
ans =
       1.0155

norm(M*xn - b)
ans =
       1.0159

And M is not actually that obscenely poorly conditioned, with a condition number on the order of 1e8.

You can do similar things using the svd. But what you DON'T want to do is compute the SVD of M'*M. Again, that is a bad thing to do. You compute the svd of M. Use it to compute a pseudo-inverse of M, as pinv does. Then the solution is trivial:

x = pinv(M)*b;

This is a nice, numerically stable solution to a nearly singular system.

[Anton Semechko]

Whooaa, your hostility is very palpable and unnecessary. If you feel compelled to offer corrections and point to more numerically stable approaches, you are welcome do to so, but that can be done without condescension.

[John D'Errico]

That is called irritation, induced by someone who is determined to lecture me about linear algebra that they don't understand. Need I quote you?

" in an overdetermined system, the matrix M will have more rows than columns, that is why I said the system must be solved using least squares. In this specific situation, the matrix that must be inverted is M'*M not M^2. When M has more rows than columns, it is not even possible to compute M^2, so I am not sure how your comment is relevant here."

I'm sorry, but that does not even make sense.

% A random matrix
v = rand(20,1);
M = v.^(0:10);

% Remember that the condition number of M is fairly big.
cond(M)
ans =
   1.4195e+08

% With a known solution vector x0
x0 = (1:11)';

% So the right hand side is given as
rhs = M*x0
rhs
rhs =
     1.2844
     9.0625
     1.0041
     4.4738
     1.3807
     1.1894
     2.1065
     13.996
      2.617
     1.6089
      1.847
     1.3431
     1.4598
     3.1627
     6.7219
     1.3192
         40
     9.2536
     3.6106
     1.3765

Just some garbage numbers.

But now, suppose we perturb the right hand side, by a tiny amount? Just junk, down in the least significant bits of rhs.

noise = randn(size(rhs))*1e-15;

Really tiny junk too. Now solve a perturbed problem, using backslash, as well as the normal equations form.

xb = M\(rhs + noise)
xb =
          1
          2
          3
          4
          5
          6
          7
          8
          9
         10
         11

As you can see, backslash recovered the original vector nicely.

xn = (M'*M)\(M'*(rhs + noise))

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  2.348908e-17. 

xn =
            1
            2
       3.0001
       3.9988
       5.0083
       5.9668
       7.0829
       7.8693
       9.1261
        9.932
       11.016

As you can see however, just that tiny amount of crap in the least significant bits was sufficient to corrupt the solution.

How do these vectors compare to the base vector?

norm(xb - x0)
ans =
   1.7191e-07

norm(xn - x0)
ans =
      0.21427

You can think of the condition number as an amplification factor. It tells you how much any junk in the data may potentially perturb the solution.

So when we solved the system using backslash, it amplified the norm of the noise I added in by a factor of roughly 1e8. Th=at just happens to be the condition number of M. But when we created M'*M, we squared that condition number, which squares the amplification factor. Now our noise is starting to show up as highly significant.

Can you form the svd from M'*M? That too is dangerous. Once you compute M'*M, the svd is no longer a stable way to solve this problem. Once you created M'*M, that blew you out of the water.

And that is why you don't solve it as you have been trying to solve. When you form the product M'*M, you create a matrix with the square of the condition number of M.

Why do I see this as a problem? Because you are teaching others to use a poor solution. Where did you learn it? Form someone else, who also did not know any better.

If nobody ever tells you why you are wrong, then you will never learn that there are both good ways and bad ways to solve the problem. Worse, you will then continue to teach others the wrong solution. And THEY will teach those who they have influence upon. This is a bad meme that propagates forever, the blind leading the blind.

[Anton Semechko]

I am familiar with the notion of a condition number, John. And so I concede to your point that the LS solution to M*F=X can be obtained without squaring the non-singular eigenvalues of M. Using a SVD-based approach, the matrix (M'*M)\M' can be rewritten as

((U*D*V')'*(U*D*V'))\(U*D*V')' =

(V*D*U'*U*D*V')\(U*D*V')' =

(V*D^2*V')\(V*D*U') =

V*inv(D^2)*V'*V*D*U' =

V*inv(D^2)*D*U' =

V*inv(D)*U' <== Moore-Penrose pseudo-inverse of M

where [U,D,V]=svd(M,'econ');

And so resulting LS solution would be

d=diag(D); d(abs(d)<eps)=eps; X=V*diag(1./diag(d))*U'*F;

[John D'Errico]

Yet, while you were able to parrot the equations from a pseudo-inverse, you were also willing to tell someone to use a poor solution that fails on singular problems, and makes nearly singular problems worse. After all, you did claim that what you proposed is THE way to solve a least squares problem. And while you claim to know what a condition number is, you seemed to be clueless about it at first. So it seems you did some recent reading, but are trying to make it seem like you understood the numerical linear algebra from the beginning.

All this on a problem that was asked because the person had rank deficient systems.

And that is why I've now had to spend an hour or so of my time to explain why what you were trying to teach was a bad thing to spread - because someone might read what you said and believe it to be good, merely because you said so.

[Anton Semechko]

John, I admit that the answer I offered at the beginning was a bad way to obtain LS solutions to overdermined systems. Thanks for educating me!