eig() shows positive eigenvalues while chol() failed

Hi guys,

I was LU-decomposing a supposed-to-be positive-definite matrix using "chol" in a loop. However, occasionally it failed with the warning

 " Error using chol
   Matrix must be positive definite. "

eig() shows that the smallest eigenvalue is 1.8*10^(-4) which I think should be well above zero. BTW, the largest one is 7.7*10^9.

How can I get around this dilemma? Your suggestions are greatly appreciated. Many thanks.

[John D'Errico]

If you were using the LU decomposition, I can firmly state that CHOL would not have produced an error. So you were NOT "LU-decomposing" as you state. Instead, you must have been doing a Cholesky decomposition, thus using CHOL. LU will not fail if the matrix is not spd.

[Allen Lau]

Hi John,

Thanks for the reply. Yes, sorry for misusing the word of "LU-decomposing". Actually, I mean Cholesky decomposition. Your codes for generating SPD matrix is great. But I don't think it can be simply applied to my case, because I can't see a reason to twist my matrix which is based on a draw from its posterior distribution. Once it is twisted, it's not a draw from the posterior any more.

I am just confused about a fundamental question i.e. why a symmetric matrix fails "chol" given its eigenvalues are positive (of course, if 10^(-4) is treated as positive in Matlab)? Can you or anyone clarify this?

Thank you.

[John D'Errico]

First of all, it is very often the case that a matrix that LOOKS symmetric is actually not so. You can test that. Try this:

norm(A - A')

or

max(abs(A - A'))

Do these tests return exact zeros? Very often the computations that are done don't create an exactly symmetric matrix. This happens because of floating point trash, that when you change the order of operations, you don't end up with exactly the same result. And in a matrix multiply, it is the BLAS that can do things like that to you.

As far as "twisting" the matrix, the difference in this case will be in the least significant bits of those numbers. If it was failing the chol test, that means your matrix is not SPD as far as MATLAB is concerned. Wanting that to be different is not enough.

Your matrix, if the spread of the eigenvalues is only 1e9 down to 1e-4 (so only 13 powers of 10) would not seem to be numerically singular. That dynamic range would normally need to be on the order of 1e16 for that to be a problem. So my guess is that your matrix is not indeed symmetric, even though you think it is so.

In any case, I have a hard time seeing why you feel it to be a problem to tweak the matrix to be usable. It was not SPD as you created it, so it was NOT truly a valid matrix for your purposes.

[Allen Lau]

Thanks for letting me know the numerical issues in Matlab. I have tested the created matrix and you are right. It's not symmetric.

Correcting the matrix by nearestSPD or other methods will change the distribution of my target. That is the last thing I would do. Thank you very much anyway.

[Ricardo Almeida]

I´m having the same problem. Without going into peculiarities of decomposition methods, I think it might be some technical issue. I'm running chol function in two different computers, both Windows 7 64bits and matlab 2015a. One flags a positive definite matrix and other don't (Maybe it's a coincidence but always return the number of columns).

[John D'Errico]

Read my comment to the other question where you made the same comment.

Answer 1

[John D'Errico]

I wrote nearestSPD (download it from the file exchange) for this purpose. It will tweak the values of your matrix slightly so that chol should always succeed. Of course, if the required tweak is necessarily large, it will still make a minimal change so that the result is as close as possible, yet will be positive definite.