Take the 2-minute tour ×
MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

Finding square root of matrices using Cholesky decomposition is limited to positive definite matrices. Any other method to find square root of matrix which has some diagonal values approximately zero (0.00000001) ?

share|improve this question
    
What is the relation between the title and the body of your question? –  Mariano Suárez-Alvarez Oct 28 '10 at 5:52
    
If you have diagonal values close to zero or zero, the matrix can be non-positive-definite (semi, negative definite). –  Anbu Oct 28 '10 at 6:04
    
By square root of $A$, one usually means a (symmetric) matrix $R$ such that $R^2=A$. On the other hand, Cholesky gives a (nonsymmetric) matrix so that $L^TL=A$. The two are not related, as far as I know. –  Federico Poloni Oct 28 '10 at 7:13
add comment

2 Answers

Actually Cholesky can be shown to extend to semidefinite matrices:

http://en.wikipedia.org/wiki/Cholesky_decomposition#Proof_for_positive_semi-definite_matrices

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.139.8064

Evidently implemented in R.

http://tolstoy.newcastle.edu.au/R/e6/help/09/04/9980.html

Meanwhile, covariance matrices are semidefinite, read some background at

http://en.wikipedia.org/wiki/Covariance_matrix

share|improve this answer
    
Note however that for stability purposes, symmetric pivoting might be required in the positive semidefinite case. –  J. M. Oct 28 '10 at 6:02
    
J.M., you should make a full answer and talk about symmetric pivoting. –  Will Jagy Oct 28 '10 at 6:11
    
Thanks Will. I'm trying to find square root of a covariance matrix. I think covariance matrix is positive definite for some probability distributions (eg. Normal distribution). –  Anbu Oct 28 '10 at 6:11
add comment

Personally, I'm still a bit torn while giving this prescription. I am of the opinion that for truly reliable computation, one should use the singular value decomposition whenever one can tolerate the performance penalty; it truly is a very stable algorithm that can also be used to compute a bunch of important diagnostic quantities. (As an aside, if it is truly the square root you want, you could exploit the fact that having a singular value decomposition of $\mathbf{A}$ is effectively equivalent to having an eigendecomposition of either of $\mathbf{A}^T \mathbf{A}$ or $\mathbf{A}\mathbf{A}^T$)

Having said this, if you're still dead-set on using Cholesky on a positive semidefinite matrix, while in exact arithmetic you are supposed to encounter a zero, what might actually happen with a (large enough) matrix with inexact entries is that your Cholesky routine encounters a tiny quantity not detected as zero, and the routine happily carries over this tiny quantity to divide other entries with. Disaster!

The cure is that one does a symmetric pivoting $\mathbf{A}\to\mathbf{P}\mathbf{A}\mathbf{P}^T$ (where $\mathbf{P}$ is a permutation matrix), which reorders the diagonal entries of $\mathbf{A}$ (no off-diagonal entries are moved into the diagonal). This has the effect that the (near-)zero quantities are not encountered until one already has proceeded through the $r$-th iteration of the main loop in the Cholesky decomposition, where $r$ is the (perceived) rank of the matrix (this depends on what tolerance you have set, or the value such that anything whose absolute value is lower than it is effectively treated as zero).

I have only given a brief sketch, since what I really should be doing is to point out this paper by Nick Higham, which discusses the nuances of Cholesky decomposition with symmetric pivoting for symmetric positive semidefinite matrices.

share|improve this answer
    
This is nice, J.M. My sense is that the OP needs a canned routine, and evidently LINPACK has that, as Higham says. Good to know there is a simple way to push all the near-zero quantities away, while still in the middle of the algorithm. –  Will Jagy Oct 28 '10 at 16:18
    
@Will: yep. Here: netlib.org/linpack/dchdc.f in fact is the LINPACK FORTRAN routine for pivoted Cholesky. However, there is now LAPACK, the successor to LINPACK, and it has a FORTRAN routine as well: netlib.org/lapack/double/dpstrf.f . –  J. M. Oct 28 '10 at 23:09
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.