[Help-gsl] inverting a symmetric positive definite matrix

[Chenyang Wang]

Dear Sir or mam,

I was using the gsl library to do some symmetric positive definite matrix
inversion task. I tried a few methods:

1. find the eigensystem of the symmetric positive definite matrix, then
invert the eigenvalue matrix.
result: takes 25 seconds to invert a 800x800 symmetric positive definite
matrix

2. use the gsl_linalg_LU_decomp and then gsl_linalg_LU_invert.
result: takes 6 seconds to invert a 800x800 symmetric positive definite
matrix

3. use gsl_linalg_cholesky_decomp and then loop all the columns and do
gsl_linalg_cholesky_solve to solve Ax=b. here, A is the matrix that needs to
be inverted, and b is vector in forms like: (1,0,0,0)' , (0,0,1,0)' etc,
being each column of an identity matrix.
result: takes 2 seconds to invert a 800x800 symmetric positive definite
matrix

Looks like using cholesky method is the fastest.

My question is: is there any faster way to invert a symmetric positive
definite matrix? the cholesky method looks good, but it did not take
advantage of the fact that the inverse of a symmetric matrix is also a
symmetric matrix. I think since many of the matrix inversions in real life
is done on symmetric positive definite matrix, it definitly helps to make a
function that inverts a symmetric positive definite matrix.

Thanks a lot!
Chenyang