Re: [Help-gsl] Inverse of a matrix

[Gongyi Liao]

On Sat, 2006-04-29 at 13:39 -0400, James Bergstra wrote:
> One reason why these guides might be confusing, is that actually you don't 
> need
> (and shouldn't want) the inverse of your covariance matrix.
> 
> Assuming your covariance matrix S is square and positive definite (which it 
> must
> be in order to be invertible) and symmetric (which it is if it is a covariance
> matrix) then what you can do is factorize the matrix into "Cholesky factors" 
> L, L'
> such that L L' = S, and L is a triangular matrix.  The gsl call which does 
> this
> is called "gsl_linalg_cholesky_decomp".  The LAPACK call which does this is
> called "dpotrf".
> 
> Your equation is actually easier to compute in the log-domain:
> 
> log(1/sqrt(det(2 pi SIGMA)) * exp(- 0.5 (Y - MU) inv(SIGMA) (Y-MU)))
> = -N/2 log(2 pi) - log(det(L)) - 0.5 ||inv(L)(Y-MU)||^2
> 
> log(det(L)) is the sum of the logarithms of the terms on the diagonal of L.
> 
> inv(L)(Y-MU) is a vector that you compute using the function
> gsl_linalg_cholesky_svx.  You first calculate Y-MU, and then the function 
> uses L to
> transform it directly to inv(L)(Y-MU), without ever inverting L.
> 
> If you are calculating the PDF of many points at once, you can get more speed 
> by
> replacing the level2 blas calls with level3 ones in the source of
> gsl_linalg_cholesky_svx.
> 
> 
> Some of this work is actually implemented in the libmsl/Layer extension in my 
> libmsl
> package.  If you want, check out the implementation of the layer_gaussian 
> routine:
> http://www-etud.iro.umontreal.ca/~bergstrj/msl/html/a00004.html#a14
> 
> The catch is that layer_gaussian expects that your covariance matrix is 
> already
> factored on input, eg. it expects L as input, not SIGMA.
> 
> 
Thanks for your comprehensive answer,  you have given a good approach
for programming, I am now trying using GSL to do MCMC simulation and
calculation, I wonder if there is some "simple and intuitive" way in GSL
for computing the inverse of a matrix without exam it is symmetric,
positive-definite or singular(if it is , we need generalized inverse)
it seems that many "ad hoc" solutions are there, thanx for your
suggestion!