fast multivariate-Gaussian evaluations

In our paper on factoring products of Gaussiasn, we recommend using the matrix inversion lemma (and we should have recommended the matrix determinant lemma too) to speed calculations. This weekend I wrote the most efficient numpy code I could to implement the log-Gaussian formulae in that paper. Here were my rules: Never construct a zero matrix or a sparse matrix using diag() or anything else that makes a container of mainly zeros. Never use inv(), only solve(). Always store and re-use elements of the calculation, especially those shared between the inversion and determinant lemmas. In the end I got a factor of thousands in speed over naive implementations, which themselves were hundreds of times faster than just constructing the matrices and operating on them.

That was fun, and a rare moment of pure coding for me.