[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Help-gsl] Inverse of a matrix
From: |
Gongyi Liao |
Subject: |
Re: [Help-gsl] Inverse of a matrix |
Date: |
Sun, 30 Apr 2006 21:34:59 +0800 |
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!