Cholesky decomposition NR and MATLAB
vivekanand
07-02-2008, 02:07 AM
Hi,
To check the positivedefiniteness of a matrix, I used the cholesky decomposition of NR and MATLAB. In the NR, I am getting that matrix is not positive definite while in MATLAB, it is showing that matrix is positive definite.
To check my code, I used some other matrices too, in those the result was consistent i.e., both were giving the same result. In this particular (attached with the post) matrix, I am getting different result. The matrix and the implemented code is attahced with this post.
Please help me out in this and suggest me the other ways to check for the positivedefiniteness of a matrix.
davekw7x
07-02-2008, 10:21 AM
...In the NR, I am getting that matrix is not positive definite I get that also. the Cholesky constructor in Numerical Recipes Version 3 throws a "Cholesky failed" exception. I get the same results with Numerical Recipes version 2 choldc function. (FYI, the failure in the "sum" loop is with a value of -0.000124042. This occurred with i == 7 and j == 7.)
while in MATLAB, it is showing that matrix is positive definite.
I don't have Matlab, but I used GNU octave, which should give the same results.
After reading in your matrix, here's what I see:
octave:2> M
M =
1.00000 0.93790 0.97230 0.93240 0.91940 0.62610 0.92160 0.91430 0.79730 0.89360
0.93790 1.00000 0.98620 0.86100 0.94930 0.70050 0.97870 0.92820 0.79650 0.95960
0.97230 0.98620 1.00000 0.87840 0.95540 0.72150 0.97300 0.96280 0.80050 0.95880
0.93240 0.86100 0.87840 1.00000 0.88940 0.40870 0.79690 0.76290 0.77950 0.74850
0.91940 0.94930 0.95540 0.88940 1.00000 0.66400 0.94070 0.92330 0.79720 0.93240
0.62610 0.70050 0.72150 0.40870 0.66400 1.00000 0.77180 0.85410 0.78000 0.82780
0.92160 0.97870 0.97300 0.79690 0.94070 0.77180 1.00000 0.95100 0.79500 0.99380
0.91430 0.92820 0.96280 0.76290 0.92330 0.85410 0.95100 1.00000 0.80170 0.96640
0.79730 0.79650 0.80050 0.77950 0.79720 0.78000 0.79500 0.80170 1.00000 0.79640
0.89360 0.95960 0.95880 0.74850 0.93240 0.82780 0.99380 0.96640 0.79640 1.00000
octave:3> chol(M)
error: chol: matrix not positive definite
I also tried
octave:4> for i=1:10
> det(M(1:i,1:i))
> end
ans = 1
ans = 0.12034
ans = 0.0010572
ans = 8.9126e-05
ans = 3.8326e-06
ans = 9.0578e-07
ans = 4.0255e-10
ans = -4.9933e-14
ans = -1.0147e-17
ans = 1.3198e-21
Since not all determinants are greater than zero, we can conclude that the matrix is not positive-definite.
Note that the failure shows up on the eighth row, consistent with my observation from the Numerical Recipes loop.
In this particular (attached with the post) matrix, I am getting different result. The matrix and the implemented code is attahced with this post.
The code in your attachment is not Numerical Recipes code. Since I don't have your Matrix class, I can't test it. If it gives different results than your Numerical Recipes, compare the code with any version of Numerical Recipes code.
A final note: The error message in your code indicates that it is testing for positive-semidefiniteness, whereas (if it actually works) it is testing for positive-definiteness.
Regards,
Dave
vivekanand
07-03-2008, 01:41 AM
Hi,
I agree with you. For the completeness I am also attaching the result from the MATLAB. It shows the matrix, then it has the result after cholesky decomposition and then I have also included the eigen values of the matrix. It shows that each eigen values are positive, hence should be Positive Definite. I have also calculated the determinants of the matrix as done by you in Octave. The result is
>> for i=1:10
det(C2(1:i,1:i))
end
ans =
1
0.1203
0.0011
8.8421e-005
3.7836e-006
8.9367e-007
5.5906e-013
3.7671e-019
2.5650e-034
5.6235e-049
It shows that all the determinants of the matrix is positive, hence positive definite.
Is GNU octave the same product as MATLAB?
I just dont know what should I do now.
davekw7x
07-03-2008, 09:07 AM
...It shows the matrixIt shows the matrix values printed to four decimal place accuracy. What were the input values of the matrix? (That is: how was the matrix created? By doing some calculations inside Matlab or what?)
Here's the reason I ask:
When I do a Singular Value Decomposition on the matrix consisting of your printed values, I get a condition number of something like 2.5e+06 (similar numbers from Octave and Numerical Recipes code). That means that computations on the matrix (including the elementary row operations used in the Cholesky decomposition) are very susceptible to roundoff errors. In other words a roundoff error of one-half of the least significant digit of a few of the coefficients can lead to enough loss of significance that makes positive-definite matrices out of non-positive-definite matrices (and vice-versa). See Footnote.
Run the following through Matlab and tell me what you get:
M = [
1.0000 0.9379 0.9723 0.9324 0.9194 0.6261 0.9216 0.9143 0.7973 0.8936;
0.9379 1.0000 0.9862 0.8610 0.9493 0.7005 0.9787 0.9282 0.7965 0.9596;
0.9723 0.9862 1.0000 0.8784 0.9554 0.7215 0.9730 0.9628 0.8005 0.9588;
0.9324 0.8610 0.8784 1.0000 0.8894 0.4087 0.7969 0.7629 0.7795 0.7485;
0.9194 0.9493 0.9554 0.8894 1.0000 0.6640 0.9407 0.9233 0.7972 0.9324;
0.6261 0.7005 0.7215 0.4087 0.6640 1.0000 0.7718 0.8541 0.7800 0.8278;
0.9216 0.9787 0.9730 0.7969 0.9407 0.7718 1.0000 0.9510 0.7950 0.9938;
0.9143 0.9282 0.9628 0.7629 0.9233 0.8541 0.9510 1.0000 0.8017 0.9664;
0.7973 0.7965 0.8005 0.7795 0.7972 0.7800 0.7950 0.8017 1.0000 0.7964;
0.8936 0.9596 0.9588 0.7485 0.9324 0.8278 0.9938 0.9664 0.7964 1.0000
]
for i=1:10
det(M(1:i,1:i))
end
chol(M)
Here's what I get from octave:
M =
1.00000 0.93790 0.97230 0.93240 0.91940 0.62610 0.92160 0.91430 0.79730 0.89360
0.93790 1.00000 0.98620 0.86100 0.94930 0.70050 0.97870 0.92820 0.79650 0.95960
0.97230 0.98620 1.00000 0.87840 0.95540 0.72150 0.97300 0.96280 0.80050 0.95880
0.93240 0.86100 0.87840 1.00000 0.88940 0.40870 0.79690 0.76290 0.77950 0.74850
0.91940 0.94930 0.95540 0.88940 1.00000 0.66400 0.94070 0.92330 0.79720 0.93240
0.62610 0.70050 0.72150 0.40870 0.66400 1.00000 0.77180 0.85410 0.78000 0.82780
0.92160 0.97870 0.97300 0.79690 0.94070 0.77180 1.00000 0.95100 0.79500 0.99380
0.91430 0.92820 0.96280 0.76290 0.92330 0.85410 0.95100 1.00000 0.80170 0.96640
0.79730 0.79650 0.80050 0.77950 0.79720 0.78000 0.79500 0.80170 1.00000 0.79640
0.89360 0.95960 0.95880 0.74850 0.93240 0.82780 0.99380 0.96640 0.79640 1.00000
ans = 1
ans = 0.12034
ans = 0.0010572
ans = 8.9126e-05
ans = 3.8326e-06
ans = 9.0578e-07
ans = 4.0255e-10
ans = -4.9933e-14
ans = -1.0147e-17
ans = 1.3198e-21
error: chol: matrix not positive definite
Is GNU octave the same product as MATLAB?
No, it is an open-source (free) Matlab look-alike with a similar user command-line interface and lots of the functionality of Matlab.
I seriously doubt that the actual matrix functions in Octave are "better than" Matlab. (In fact, I would almost guarantee that they are not better.) I just have a feeling that the matrix that you printed out is not the matrix that your Matlab program is dealing with internally.
Regards,
Dave
Footnote:
Notice the extremely small magnitudes (relative to the magnitude of the largest eigenvalue) of the last three eigenvalues that you posted. That's a really big clue about the ill-conditioning of the matrix.
"The purpose of computing is insight, not numbers."
---Richard W. Hamming
vivekanand
07-04-2008, 07:35 AM
Hi
yes you are correct. This difference is due to some roundoff errors. Actually I have generated that Matrix from another matrix and internally MATLAB is using some other matrix.
Thanks alot for your help.