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

13 views (last 30 days)
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);
  2 Comments
Walter Roberson
Walter Roberson on 6 Dec 2022
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
Hung Dao on 9 Dec 2022
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 !

Sign in to comment.

Answers (2)

Bruno Luong
Bruno Luong on 9 Dec 2022
Edited: Bruno Luong on 9 Dec 2022
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.
  2 Comments
Hung Dao
Hung Dao on 9 Dec 2022
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
Bruno Luong on 9 Dec 2022
Edited: Bruno Luong on 9 Dec 2022
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.

Sign in to comment.


Askic V
Askic V on 6 Dec 2022
Edited: Walter Roberson on 7 Dec 2022
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:
  4 Comments
Bruno Luong
Bruno Luong on 9 Dec 2022
Edited: Bruno Luong on 9 Dec 2022
@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.

Sign in to comment.

Categories

Find more on Linear Algebra in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!