Re: [Help-gsl] gsl_linalg_LU_invert

[Patrick Alken]

Hi,

  If you're solving generalized eigensystems Ax = sBx where
A and B are symmetric, and B is positive definite, then you
want to be using the Cholesky decomposition instead of LU.
The reason is if you compute A B^{-1}, then you have a
nonsymmetric eigenvalue problem to solve, which requires more
work than you need for this. If you use the Cholesky decomposition
on B, you will end up with a symmetric eigenvalue problem to solve:

A x = s B x
    = s L L^t x
A L^{-t} L^t x = s L L^t x
[ L^{-1} A L^{-t} ] L^t x = s x

so now you must solve:

C y = s y
C = [L^{-1} A L^{-t} ]
y = L^t x

note that C is symmetric, so you can use the symmetric eigensolver
which is n^2 instead of n^3. The main problem here is to compute
the matrix C efficiently. The lapack routine dsygst will show you
how to do this. Once you have C, just use gsl_eigen_symmv() to
get the eigenvalues and eigenvectors, then use one of the triangular
BLAS solvers to get the 'x' eigenvectors from the 'y' eigenvectors.

If performance isn't a factor for you, then it would be simpler
just to compute A B^{-1} and use gsl_eigen_nonsymmv(). The
symmetric-definite generalized eigenproblem is pretty straightforward
(compared to the nonsymmetric case), so it will probably be
implemented in gsl eventually.

Patrick Alken