Re: calculating matrix inverses (was Re: modifying matrix allocationfunctions for use with R)

[Edwin Robert Tisdale <E.Robert.Tisdale@jpl.nasa.gov>]

Faheem Mitha wrote:

> I want to calculate the inverse of (HH^t)^{-1} H
> where H is a given matrix
> as part of simulating from a multivariate normal distribution
> (the given expression is the mean).

For a nonsingular square matrix

	A = (HH^T)

you need to solve

	AX = I

where I is the identity matrix and

	X = A^{-1}

is the inverse of A.

The problem with matrix inversion is that
there is no one right way to do it.
Perhaps, in this case, a Cholesky decomposition
might be the preferred solution method.
Computing the inverse in two steps --
decomposition followed by back substitution --
has the added advantage of allowing you
to evaluate the condition of the matrix
after the decomposition and
before you attempt to divide by zero
during the back substitution.
to divide by zero