help-gsl
[Top][All Lists]

## Re: [Help-gsl] Inverse of a matrix

 From: James Bergstra Subject: Re: [Help-gsl] Inverse of a matrix Date: Sat, 29 Apr 2006 13:39:56 -0400 User-agent: Mutt/1.4.2.1i

```On Sat, Apr 29, 2006 at 03:07:17PM +0800, 廖宮毅 wrote:
> Thanks for attention, I have a question about computing the inverse of a
> matrix, for example, I want to calculate the value of the P.D.F. with a
> set of sample(size N) Y from a multivariate normal distribution
> (p-variate), with mean MU and covariance matrix SIGMA,  the P.D.F is :
>    det(SIGMA)^(0.5*N)*exp(-0.5*(Y - MU)%*%inv(SIGMA)%*%(Y -MU))
> my problem is: how to use GSL calculating the inverse of SIGMA?
> I have found some solutions on the mailing list, but the solutions are
> not easy to understand (due to my lack of knowledge), is there any
> suggestion computing the inverse of a matrix (in a clear and direct
> way)? Any comment is appreciated!

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.

--
james bergstra
http://www-etud.iro.umontreal.ca/~bergstrj

```