[Numpy-discussion] inversion of large matrices

Daniel Elliott danelliottster@gmail....
Mon Oct 11 20:05:51 CDT 2010


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





More information about the NumPy-Discussion mailing list