Does the function chol correctly indicates that a Matrix is positive definite ?

According to the MATLAB documentation for the function chol: "[R,p] = chol(A) for positive definite A, produces an upper triangular matrix R from the diagonal and upper triangle of matrix A, satisfying the equation R'*R=A and p is zero. If A is not positive definite, then p is a positive integer and MATLAB® does not generate an error"

I have found an example contradicting the behavior described so I do not understand if I can trust this function or whether there is an error in the documentation.

Let's have the following Matrix.

A = 
[[0.957166948  0.421761283  0.655740699  0.655740699];
[0.485375649  0.915735525  0.035711679  0.035711679];
[0.800280469  0.79220733  0.849129306  0.849129306];
[0.141886339  0.959492426  0.933993248  0.933993248]];

Then the Cholesky factorization gives the following result:

[R,p] = chol(A)

R =
      0.9783    0.4311    0.6703    0.6703
           0    0.8543   -0.2964   -0.2964
           0         0    0.5586    0.5586
           0         0         0    0.2913
p =
     0

So the p = 0 indicates according to the interpretation I have of the documentation that A is positive definite, which also indicates that the matrix A is nonsingular.

However this matrix is clearly singular given that the last two columns are identical. When computing the eigenvalues I get a zero eigenvalue which confirms this singularity. Also the determinant gives a numerical zero value.

> eig(A)

ans =

     2.5108 + 0.0000i
    -0.0000 + 0.0000i
     0.5726 + 0.3574i
     0.5726 - 0.3574i

>> abs(det(A)) < eps

ans =

    logical

     1

Conclusion: I can not conclude from P=0 in chol that the matrix is Positive definite.

Is there an error in this function, in the documentation or I am missing/misunderstanding something ?

I appreciate your help with this, dear MATLAB users.

Answer 1

[John D'Errico]

For the 10 millionth time, det is a bad way to check is a matrix is singular. (Not your fault for asking though.) In fact, it is NEVER a good way to do so, except in homework, where teachers love to tell you to use det, but not explain why det is bad so often. Sadly, that propagates, because their students will do the same thing. And if nobody ever tells you, then your potential students would learn the same fallacies.

Next, giving us a matrix of elements that are only written in 9 decimal digits of precision for a matrix of doubles tells us nothing.

You are confusing the use of chol to test for a positive definite matrix, with testing for singularity. As well, the matrix you have shown is not even symmetric. That tells me it will usually have complex eigenvalues.

I am a bit surprised that chol does not test to see if the metrix is symmetric. But it looks as if chol only uses the upper triangle of the input array. We can verify that by the following test:

[R,P] = chol(A);
A - R'*R
ans =
 1.1102e-16   5.5511e-17   1.1102e-16   1.1102e-16
   0.063614            0   6.9389e-18   6.9389e-18
    0.14454       0.7565            0            0
   -0.51385      0.92378     0.084864            0

So chol is indeed only looking at the upper triangle. (I think that chol used to test for symmetry too. So that may have changed in some recent release.)

If you want to test to see if A is SINGULAR, NEVER use det. Chol is completely inappropriate too.

rank(A)
ans =
     3

A non-singular 4x4 matrix would have a rank of 4. So A is rank deficient.

Understanding the SVD would be even better for some people, but rank is perfectly adequate in this respect.

[Juan Diego Quesada]

I would like to thank you for your explanation. You mentioned that "You are confusing the use of chol to test for a positive definite matrix, with testing for singularity", but there is a link between positive definite and singularity. A positive definite matrix is non singular. I was already aware that the determinant is not a good method to test for singularity, that is way I also displayed the eigenvalues. I am aware of the rank method but it is extremely slow. Is there a fast way (I found the chol as a good candidate) to determine that a covariance matrix is "well defined" ? I wanted first to check for positive definite at first and only if it is not then check for PSD for instance. The eigenvalues function is also quite slow compared with chol.

[Juan Diego Quesada]

I think an important point you made is the fact that I used a non symmetric matrix as example. I think if the matrix is symmetric the chol function will always correctly indicate the positive definite nature. I will test this. Thanks.

[John D'Errico]

NO. I'm sorry, but you are mistaken. There is a link between singularity and positive definite, but that link goes only ONE way. And the link is a terribly poor one when implemented with floating point numbers.

It is true that a positive definite matrix will be non-singular, at least in terms of mathematics.

By the way, if you will go about telling people that a matrix is singular or not, and trying to prove that based on the determinant computed using floating point arithmetic, then be prepared to have people tell you that you don't know what you are talking about. And they would be right. Unless of course, you can prove that

A = 0.1*eye(1000);
det(A)
ans =
     0

is actually singular as det seems to think, or that it is in any way materially different from

A = 10*eye(1000);

As far as negative eigenvalues go, it is very easy to create a matrix that has negative eigenvalues, yet is not singular.

A = ones(2) - eye(2)/2
A =
          0.5            1
            1          0.5

eig(A)
ans =
         -0.5
          1.5

A is even symmetric, is not singular. It is obviously not positive definite. Although chol will correctly identify A as not positive here.

chol(A)
Error using chol
Matrix must be positive definite.

So chol cannot be used to identify if a matrix is singular.

In fact, a matrix can even be numerically semi-definite (which is really all that chol can test for), and still be numerically singular. CHOL can be easily confused, with a matrix that is effectively numerically singular.

As an example, I'll use a simple function I wrote for a recent totally different answer.

Cfun = @(k) (1-k)*eye(4) + k*ones(4);

A = Cfun(-1/3 + eps)
A =
            1     -0.33333     -0.33333     -0.33333
     -0.33333            1     -0.33333     -0.33333
     -0.33333     -0.33333            1     -0.33333
     -0.33333     -0.33333     -0.33333            1

So A is a simple matrix. (Be careful, as I wrote out only 4 digits there in those numbers.) Chol actually succeeds here.

chol(A)
ans =
            1     -0.33333     -0.33333     -0.33333
            0      0.94281      -0.4714      -0.4714
            0            0       0.8165      -0.8165
            0            0            0   5.6742e-08

But, in fact, A is numerically singular, certainly by any rational test for singularity.

rank(A)
ans =
     3

svd(A)
ans =
       1.3333
       1.3333
       1.3333
   8.0491e-16

That the matrix is symmetric or not, chol does not seem to test for, nor does it even use the lower triangle of the input. I showed that for your matrix.

[John D'Errico]

Your last question is how best to test if the matrix is positive definite. chol is the accepted test in MATLAB, because even if the matrix is semi-definite and chol succeeds, then essentially anything you will do with a covariance matrix, the Cholesky factor is all you will need. So chol is essentially sufficient as a test.

Rank is more accurate as a test for singularity, but as you say, rank is slower.

Perhaps what you really need is a tool that will ensure that your matrix is SPD. A few years ago, I wrote nearestSPD, based on an algorithm from Nick Higham.

A = rand(3)
A =
      0.25754      0.73656      0.55376
      0.21045      0.99486     0.089536
      0.91036      0.86264      0.89538

chol(A)
Error using chol
Matrix must be positive definite. 

Ahat = nearestSPD(A)
Ahat =
      0.43488      0.44497      0.62977
      0.44497      0.99946      0.49255
      0.62977      0.49255      0.95438

chol(Ahat)
ans =
      0.65945      0.67475      0.95499
            0      0.73768     -0.20583
            0            0   1.5805e-08

nearestSPD is on the File Exchange.

A = Cfun(-1/3 - eps);
chol(A)
Error using chol
Matrix must be positive definite. 

Ahat = nearestSPD(A);

chol(Ahat)
ans =
            1     -0.33333     -0.33333     -0.33333
            0      0.94281      -0.4714      -0.4714
            0            0       0.8165      -0.8165
            0            0            0   1.0537e-08

norm(A-Ahat)
ans =
   7.0009e-16

Most of the time, the perturbation is a tiny one, since most of the time, one would only apply that code to a matrix where you think it is at least semi-definite.

[Steven Lord]

John wrote "I am a bit surprised that chol does not test to see if the metrix is symmetric. But it looks as if chol only uses the upper triangle of the input array." That's almost correct. From the first line in the Description section on the documentation page for chol:

" R = chol(A) produces an upper triangular matrix R from the diagonal and upper triangle of matrix A, satisfying the equation R'*R=A."

The third paragraph in the Description section of that page covers the two output syntax:

" [R,p] = chol(A) for positive definite A, produces an upper triangular matrix R from the diagonal and upper triangle of matrix A, satisfying the equation R'*R=A and p is zero."

chol uses the upper triangle and the main diagonal. So while the matrix the asker passed into chol was not positive definite, the matrix whose Cholesky factorization chol computed was.

[Juan Diego Quesada]

Very clear explanation. I appreciate you took the time to help me with these questions. I have accepted your answer.

Answer 2

[Zheng]

Thanks John!

I have also encountered this problem when studying "chol()" and luckily found someone else has realized this. What's more, I found that Chol will behave quite differently towards two similar matrices even if we do not consider decimal digits. I have tried some examples and wrote down here. My current version is Matlab r2022a.

The base matrix is just the first example from the Chol help documents.

A = [1 0 1; 0 2 0; 1 0 3]

A = 3×3

1 0 1 0 2 0 1 0 3

• 1. If we change the element in A(3,1) to any other number except "1", for instance "1000", Chol will always regard this matrix as positive definite and return "flag = 0"

A1 = [1 0 1; 0 2 0; 1000 0 3]

A1 = 3×3

1 0 1 0 2 0 1000 0 3

[R1,flag1] = chol(A1)

R1 = 3×3

1.0000 0 1.0000 0 1.4142 0 0 0 1.4142

flag1 = 0

• 2. However, if we change the element in A(1,3) to any other number except "1", for instance "2", Chol will always correctly detect non-positive definite matrix and return non-zero values.

A2 = [1 0 2; 0 2 0; 1 0 3]

A2 = 3×3

1 0 2 0 2 0 1 0 3

[R2,flag2] = chol(A2)

R2 = 2×2

1.0000 0 0 1.4142

flag2 = 3

[Steven Lord]

This is Not a Bug.

From the first entry in the Description section of the documentation for the chol function: "R = chol(A) factorizes symmetric positive definite matrix A into an upper triangular R that satisfies A = R'*R. If A is nonsymmetric , then chol treats the matrix as symmetric and uses only the diagonal and upper triangle of A." [Emphasis added.]

In your first case you're changing the lower triangular part of A which chol doesn't use. In your second case you're changing the upper triangular part of A which it does.

[Zheng]

Thanks Steven.

Yes, I found the description which indicates how it treats a nonsymmetric matrix. However, the thing is that both A1 and A2 are nonsymmetric in my examples, if I am right, while Chol treat them differently. You can also see that from the value of "flag". I am afraid that both of them should be non-zero numbers.

[Zheng]

Ah I think I found the reason.

The second example, A2, could not be positive definite under this case even if Chol only uses the upper triangular part of A2. Symmetry is just a necessary condition...

[Steven Lord]

if I am right, while Chol treat them differently.

Yes, that is correct and is the correct behavior. The chol function ignores the lower triangular part of the input entirely unless you explicitly specify the second input as 'lower' to tell it to ignore the upper triangular part instead.

From the description of the A input argument in the documentation page to which I linked in my previous comment: "chol assumes that A is symmetric for real matrices or Hermitian for complex matrices. chol uses only the upper or lower triangle of A to perform its computations, depending on the value of triangle." [Emphasis added.]