Symmetric, positive-semidefinite matrix has (big) negative eigenvalues

Hi there,

when I run the following code:

whos J;
JProd = (J.'*J);
whos JProd;
eigenVals = eig(JProd);
min(eigenVals)

I get the output: Name Size Bytes Class Attributes

    J         65x2963            1540760  double              

    Name          Size                 Bytes  Class     Attributes

    JProd      2963x2963            70234952  double              

ans =

    -1.0126e+03

How can the smallest eigenvalue be -1012? The matrix JProd should only have non-negative eigenvalues, right? Is there anything wrong with my syntax above when multiplying the matrix J with its transposed?

Thanks in advance, Joerg

Answer 1

[Matt J]

There's no immediate reason to think that -1012 is big. We haven't seen how big it is relative to the other eigenvalues. If max(eigenvals) = 1e23 or something huge like that, then -1012 is obviously just round-off error in the eig() computation.

[Joerg]

Thanks a lot for your answer, Matt. Okay, fair enough, "big" is relative. I will check max(eigenVals) tomorrow at work.

Nevertheless, that would mean that I would have to live with those round-off errors? My original problem is that I want to invert the matrix

JProd + alpha*eye(2963)

for alpha a double value (smaller than 10). If the eigenvalues are all poisitive that should result in the above matrix being regular. However, MATLAB keeps telling me that this matrix is almost singular. If the rounding errors are so big, what can I do? I mean, I can't change my mathematical model or the theoretical algorithm, just because MATLAB produces rounding errors, can I?

Do you have any advice for me? I would be so grateful to you.

[Matt J]

If the eigenvalues are all poisitive that should result in the above matrix being regular.

Well, no, even that's not true. Consider the following matrix which definitely has strictly positive eigenvalues. Still, it has problems when you want to invert it numerically,

    >> A=diag([1e20,1]); A\[1;1]
  Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
  RCOND =  1.000000e-20.

You shouldn't be composing the matrix

    (J.'*J)+alpha*eye(2963)

with the intention of doing

    inv((J.'*J)+alpha*eye(2963))*something

You should be doing this instead

    [J;sqrt(alpha)*eye(2963)]\something

Note that rcond(J.'*J) is the square of rcond(J), so the former form makes your equations less stably invertible.

I mean, I can't change my mathematical model or the theoretical algorithm, just because MATLAB produces rounding errors, can I?

Yes, that's exactly what you need to do. If it were a good model/algorithm, it wouldn't lead you to a situation where you need to invert a poorly conditioned matrix.