How one weird trick helps you evaluate correlated Normal distributions quickly


Sorry, couldn’t resist the opportunity to buzzfeed this research blog :)

I’ve been trying to get really efficient at writing samplers for various Bayesian spatial models. And, typically, this involves clever numerical tricks, trying to avoid computing either log determinants or matrix inverses by keeping around matrix factorizations or finding derived products that you can persist between iterations. Some of this has come to fruition in my sparse log determinant work, but I’m always looking for computation speed gains, especially as some of the targets I’ve optimized are drying up.

One I’ve just identified might be in sampling kernels that look like multivariate normal kernels. Naively, when you write the logp of a normal distribution, you’d expect to do the following kinds of computations. For respone Y, linear predictor XBeta, and a dense correlated covariance matrix Sigma, I’d naively compute a normal kernel using the following code:

e = Y - XBeta
Sigma_i = np.linalg.inv(Sigma)
kernel = -.5*np.linalg.multi_dot(e.T, Sigma_i, e.T)

But, there’s actually a siginficant gain if you use the solve method instead of inverting Sigma by itself.

e = Y - XBeta
kernel = np.linalg.solve(Sigma, e)

Who knew? How much faster? Timings are attached for covariance matrices from 5^2 to 20^2 in dimension. Left is the solve method, right is the inverse method. While I’m scaling these up for much larger spatially-correlated error models as I write, it looks to me like a promising speed gain!

edit: scaling quite nicely up to 50^2. I’ll have to see where the sparse math goes with this trick. 

imported from: yetanothergeographer