Error using chol Matrix must be positive definite.

I have a positive definite matrix C for which R=chol(C) works well. I want to apply the chol function to a new matrix A = U*C*U' where U is a unitary matrix obtained as output from SVD, i.e. from [V,S,U] = dvd(T); but I get an error telling me that A is not positive definite. I checked that det(U) = 1.0 so I don't understand why the symmetric matrix A is not positive definite.

Given that C is positive definite then y'*C*y>0 and if I let y = U'*x then x'*U*C*U'*x>0 which implies that U*C*U'is also positive definite.

Is this problem due to round off or am I missing some important linear algebra concept. If not is there a way around this problem?

[Matt J]

Attach a .mat file with C and U. We'll need to play with the data.

[John D'Errico]

That det(A)==1 is NOT any assurance that the matrix is not numerically singular. In fact, it is trivial to create a matrix that has a determinant equal to ANY value, yet it still be singular in double precision.

A = diag([1e10,1e-10]);

Clearly, the determinant is 1. As clearly, it is also effectively a numerically singular matrix in double precision.

A = diag([1e10,1e-10]);
det(A)
ans =
     1

rank(A)
ans =
     1

cond(A)
ans =
                     1e+20

NEVER use the determinant as a measure of singularity. NEVER. NEVER. NEVER. That you may have read it in a book is irrelevant. That you may have seen it in some text that is 40 years old is irrelevant. Sadly, the authors of books today are still referring back to those texts they learned from 40+ years ago, still teaching their own students the wrong things about numerical methods.

[Youssef Khmou]

that is correct, what about the condition number : lambda_max/lambda_min ?

[John D'Errico]

cond returns that value. see my example. It is a good predictor of numerical singularity, certainly far better than det. Of course, a random number generator can be as good as det in that respect.

[Francois]

Thank you all for your answers and suggestions. The conditioning of my matrix was indeed the problem. It turned out that my matrix U was well conditioned (condition number of 1) but my matrix C was not. It had a condition number on the order of 2*10^24. I guess the fact that chol(C) worked ok was just a fluke. I tried the nearestSPD and it worked well.

Thanks for the quick and most useful advice.

[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]

The most common reason for this is NOT the difference in code, which should not be, but how you pass the array between. Too often people think they can pass an ascii file between the two machines, that this is sufficient.

Unless the array is passed EXACTLY between machines as a .mat file, you are NOT making a proper comparison. Without use of a .mat file, there will be tiny errors in the least significant bits.

[Zhiyong Niu]

The diagnal of a positive definite matrix is real. However, if you obtain A by A = U*C*U' ,the diagnal of A may have imagenary parts, even though they are extremely tiny, on the order of 1e-17i. if so, the chol() may give you an error when the elements of diagnal was checked.

Answer 1

[John D'Errico]

One solution is to use my nearestSPD code as found on the file exchange. It handles the semi-definite matrix, finding the smallest perturbation into a positive definite matrix, one that will be ASSUREDLY factorizable using chol.

[Junho Kweon]

I totally agree with this. Due to the limitation of the computation in MATLAB, sometimes minor computational error cause sight asymmetric or slightly-not-positive-definite matrix. I had same issue with @Francois, but the code as @John D'Errico showed solved this problem very well.

Answer 2

[Youssef Khmou]

this an interesting problem,

Generally, the matrix C must contain some negative and positive eigenvalues ( eig(C)) according the description, in the other hand, the matrix A is positive semi definite only if C is diagonal matrix with the diagonal elements being the eigenvalues corresponding the eigenvectors U(:,1),....U(:,N).

In this case you multiply C whether diagonal or not with non corresponding eigenvectors, so A can not be positive semi definite .

[Matt J]

A is positive semi definite only if C is diagonal matrix with the diagonal elements being the eigenvalues corresponding the eigenvectors U(:,1),....U(:,N).

Not true. Suppose U=eye(N). Then A=C and both are positive (semi) definite simultaneously, regardless of whether C is diagonal.