[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