[Numpy-discussion] inversion of large matrices

[Daniel Elliott]

Sturla Molden <sturla <at> molden.no> writes:

> > Yes, this is what I am computing.  I am computing the pdf of a very high-
> > dimensional multivariate normal.  Is there a specialized method to compute
> > this?
> 
> If you use cho_solve and cho_factor from scipy.linalg, you can proceed
> like this:
> 
>   cx = X - m
>   sqmahal = (cx*cho_solve(cho_factor(S),cx.T).T).sum(axis=1)
> 
> where X is the data (n x p), m is mean (1 x p) and S is the covariance (p
> x p).
> 
> We could do this twice as fast using a hypothetical method tri_solve that
> calls LAPACK subroutine DTRTRS.

Sorry for the laziness of this question...

As a reminder, I am working with high-dimensional data (10K dimensions).  I was 
computing the log of the MVN pdf because the probabilities would almost all go 
to zero.

Do you suppose the method you have shown me will be numerically stable (the 
probabilities will be small but they stay above zero)?

By the way, I would be happy to implement the method you desire in a couple 
months if you are willing to do a little hand-holding.

Thanks for the excellent suggestion.

- dan