chol() returns error even the matrix is positive definite.

I have a matrix , where I is a identity matrix and B is a matrix. It is clear that M is a positive definite matrix. However, I am facing a situation in which the Matlab Cholesky decomposition function

chol(M);

returns an error stating that M must be positive definite. I examine matrix B and realise that it's elements are pretty small, resulting in the elements in matrix are extreme large in magnitude (the magnitudes are about e+15 to e+17).

Any advice/suggestion how to handle this issue is much appreciated !

Update on 9 Dec 22:

To help explaning the issue more cleary, I attached the data containing the matrix. Please load the data 'example.mat' and run the code below

[p, r] = size(B);    
    B1 = d.\B; 
    z1 = z./d;  
    z2 = B1'*z1; 
    R = chol(eye(r) + B1'*B1);

[Walter Roberson]

Taking B'*B cannot result in large values if the entries are all small... not unless there are a huge number of entries. Each element of B'*B is the dot product between a row and another row of B. The maximum possible value would be max(abs(B))^2 * number of columns

Perhaps you are looking at the pseudoinverse of the matrix instead of the transpose

[Hung Dao]

Thank you for your comment. My apologies for my late reply as I am travelling. You are right. Matrix B is quite a big matrix (it has the size of ). Many of the elements are small (23,036 entries are smaller than 1e-1) but I do realise that B doeve has some very large elements (26 elements are greater than 1e+6 with the maximum is about 1e+9).

I attached the data and the code. Hopefully it helps explain my issue better.

Thank you for answering my question. I appreciate it !

Answer 1

[Bruno Luong]

This statement

B1 = d.\B; 

looks very suspecious to me (the .\ operator), if you remove the dot it will work.

Here you do element-wise of d that oscillate to 0, create a huge number in B1 (it can goes to infinity indeed s such noisy data)

[p, r] = size(B);    
B1 = d\B;  
R = chol(eye(r) + B1'*B1);

Now wht you trying to do and implement only you know, but I suspect your code is incorrect for whatever the goal you want to achieve.

[Hung Dao]

Thank you for your answer. I appreciate your support. However, I respectfully disagree with you.

I am using

B1 = d.\B

as I want to calculate matrix , where , a diagonal matrix with the main diagonal is vector d. Hence, replacing my formula using the .\ operator by

B1 = d\B

is incorrect.

[Bruno Luong]

Only you know what is correct, but your d looks like this

It looks like it contains noise close to 0 and it cross randomly 0 mny time. 1./ d.\1 is then randomy big.

It that make sense for you then OK, but don't expect getting anything meaningful from that.

B1'*B1 have condition number and large element of the order of probably 1e20. Add eye(r) to it is meaning less, doing cholesky doesn't make any sense numerically. Your code won't produce anything usuful but gardebage with such data.

Answer 2

[Askic V]

Without actual code it is very hard to know for sure, but using chol() function in Matlab is actually a preferred way to test if the matrix is positive definite.

For me, it is very hard to understand how values of B'*B can be extremly large if values of B are small. Something is wrong. Please have a look at the example:

B = randn(3,3)*1e-6;
C = B'*B

C = 3×3

1.0e-11 * 0.3916 0.0891 -0.0314 0.0891 0.0303 -0.0417 -0.0314 -0.0417 0.1399

max(C(:))

ans = 3.9156e-12

min(C(:))

ans = -4.1702e-13

Regarding your problem, I suggest you to read these advices:

https://www.mathworks.com/matlabcentral/answers/225177-chol-gives-error-for-a-barely-positive-definite-matrix?s_tid=answers_rc1-1_p1_MLT

[Hung Dao]

Thank you very much for your prompt answer. My apologies for my late reply as I am travelling. I cannot provide the actual code because this calculation is a small part of my algorithm. In addition, matrix B is quite a big matrix (it has the size of ), hence it is hard to tell what results in an issue. Many of the elements are small (23,036 entries are smaller than 1e-1) but I do realise that B doeve has some very large elements (26 elements are greater than 1e+6 with the maximum is about 1e+9).

Attached is the Matlad data.

If you load the data and run the code below, you will see the issue I am facing better :)

[p, r] = size(B);    
    B1 = d.\B; 
    z1 = z./d;  
    z2 = B1'*z1; 
    R = chol(eye(r) + B1'*B1);

[Askic V]

Well, by definition, Cholesky decomposition or Cholesky factorization is a decomposition of a Hermitian, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose.

In this particular example, you do have symmetric but not positive definite matrix, so you cannot apply chol decomposition.

C = B1'*B1;
[V,D] = eig(C);

min(D(:))

This shows that you have negative eigenvalues, so C is not positive definite. The minimal value is -2773.5, so adding identity matrix will not change anything.

I'm not sure about the meaning of your data in the columns, but it could be possible to first normlaize data (for example using z score) and then apply cholesky decomposition if that will not change the meaning and purpose of the application.

For example,

[Bn, Bc, Bs] = normalize(B1);

R = chol(eye(r) + Bn'*Bn);

[Bruno Luong]

@Askic V "In this particular example, you do have symmetric but not positive definite matrix"

I claim is wrong

C = eye(r) + B1'*B1

is always spd theorically since

x'*C*x is equal to |x|^2+|B1*x|^ 2 >= |x|^2 > 0 for any vector x ~= 0. Which is the requirement of spd (symmetric is easy to see).

Even B1'*B1 is non negative, eig fails to estimate the smallest eigen value.

OP simply has numerical error due to bad scaling issue.

[Askic V]

I stand corrected. Thank you Bruno.