[Numpy-discussion] inversion of large matrices

[David Warde-Farley]

On 2010-08-30, at 11:28 AM, Daniel Elliott wrote:

> Hello,
> 
> I am new to Python (coming from R and Matlab/Octave).  I was preparing
> to write my usual compute pdf of a really high dimensional (e.g. 10000
> dimensions) Gaussian code in Python but I noticed that numpy had a
> function for computing the log determinant in these situations.

Yep. Keep in mind this is a fairly recent addition, in 1.5 I think, so if you ship code make sure to list this dependency.

> 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.

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.

David