[Numpy-discussion] inversion of large matrices

[Dan Elliott]

Thanks for the reply.

David Warde-Farley <dwf <at> cs.toronto.edu> writes:
> On 2010-08-30, at 11:28 AM, Daniel Elliott wrote:
> > Large matrices (e.g. 10K x 10K)
> 
> > Is there a function for performing the inverse or even the pdf of a
> > multinomial normal in these situations as well?
> 
> There's a function for the inverse, but you almost never want to use it, 
especially if your goal is the
> multivariate normal density. A basic explanation of why is available here:
> http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/ 
> 
> In the case of the multivariate normal density the covariance is assumed to be 
positive definite, and thus a
> Cholesky decomposition is appropriate. scipy.linalg.solve() (NOT 
numpy.linalg.solve()) with the
> sym_pos=True argument will do this for you.

You don't think this will choke on a large (e.g. 10K x 10K) covariance matrix?  

Given what you know about how it computes the log determinant and how the 
Cholesky decomposition, do you suppose I would be better off using eigen-
decomposition to do this since I will also get the determinant from the sum of 
the logs of the eigenvalues?

> What do you mean by a "multinomial normal"? Do you mean multivariate normal? 
Offhand I can't remember if it
> exists in scipy.stats, but I know there's code for it in PyMC.

Oops, I meant multivariate normal.  I found code in the pymix library and so far 
I like what I see there.  I am sure it works well for most cases.  However, with 
the exception of one R library, I have never found a packaged implementation 
that works on covariance matrices as large as what I am talking about above.