Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Obsolete Editions Forum > Methods: Chapters 2, 11, and 18

Reply
 
Thread Tools Display Modes
  #1  
Old 07-02-2008, 03:07 AM
vivekanand vivekanand is offline
Registered User
 
Join Date: Dec 2007
Posts: 5
Cholesky decomposition NR and MATLAB

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.
Attached Files
File Type: txt matrix.txt (1.5 KB, 1408 views)

Last edited by vivekanand; 07-02-2008 at 06:13 AM. Reason: corrected sentences
Reply With Quote
  #2  
Old 07-02-2008, 11:21 AM
davekw7x davekw7x is offline
Registered User
 
Join Date: Jan 2008
Posts: 453
Quote:
Originally Posted by vivekanand View Post
...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.)


Quote:
Originally Posted by vivekanand
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:
Code:
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
Code:
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.

Quote:
Originally Posted by vivekanand View Post
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

Last edited by davekw7x; 07-02-2008 at 01:27 PM.
Reply With Quote
  #3  
Old 07-03-2008, 02:41 AM
vivekanand vivekanand is offline
Registered User
 
Join Date: Dec 2007
Posts: 5
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.
Attached Files
File Type: doc matlab.doc (2.9 KB, 1108 views)

Last edited by vivekanand; 07-03-2008 at 02:46 AM. Reason: added information
Reply With Quote
  #4  
Old 07-03-2008, 10:07 AM
davekw7x davekw7x is offline
Registered User
 
Join Date: Jan 2008
Posts: 453
Quote:
Originally Posted by vivekanand View Post
...It shows the matrix
It 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:

Code:
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:
Code:
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
Quote:
Originally Posted by vivekanand
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

Last edited by davekw7x; 07-04-2008 at 09:42 AM.
Reply With Quote
  #5  
Old 07-04-2008, 08:35 AM
vivekanand vivekanand is offline
Registered User
 
Join Date: Dec 2007
Posts: 5
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.
Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump


All times are GMT -5. The time now is 01:24 PM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.