Multivariate normal covariance matrix

Dear Stan gurus,

I have a simple question and hope it makes sense.

Let Y_{(s)} be a response variable, being

Y_{(s)} \sim \mathcal{N}(\mu_{(s)}, \Sigma_{(s)} )

where \Sigma_{(s)} is a spatial covariance matrix.
If the inverse of \Sigma_{(s)} (e.g. Sigma_inverse in the codes) is available and it is a sparse matrix.

Would be (1) better and (2) faster to use

target += multi_normal_prec_lpdf(y| mu, Sigma_inverse );

rather than

target += multi_normal_cholesky_lpdf(y| mu, cholesky_decompose(Sigma) );

Thank you

Without looking at how exactly Stan implemented the likelihood from either matrix I’d say using the inverse (\Sigma^{-1}) is fastest (if you already have it), since it involves simply computing the matrix operation (\mathbf{x} - \mathbf{\mu})^T \Sigma^{-1} (\mathbf{x} - \mathbf{\mu}) (the determinant |\Sigma| is just 1/|\Sigma^{-1}|).

Computing a Cholesky decomposition is usually preferred to computing the inverse, but there you’d be computing the inverse and then the decomposition, right? And whether in practice it’s better to do it separately and then computing the likelihood will depend on whether there’s some clever mathematical or computational optimization to the function. Like I said with regard to that I don’t know what is preferred.

1 Like

It’s a good question, and unfortunately the answer is ‘it depends’. The multi_normal_cholesky distribution is more optimised and has analytic gradients, whereas the multi_normal_precision distribution relies on autodiff. Additionally, if your matrix is large enough, the cholesky decomposition can be GPU-accelerated.

So it’s really a case of try both and see how it works


Hi @caesoma and @andrjohns Thank you for your prompt reply. It is very helpful.

@andrjohns is the GPU-based option that you mentioned what it was introduced in the paper?
OR is there another option/tutorial?

That was the initial introduction yes. To use the GPU-accelerated functions you’ll need to use the cmdstanr package instead of rstan (or you can use the cmdstanPy package for python), and setup the appropriate OpenCL drivers. You can follow this guide to get that all setup:

And then you can run your models as normal, and any functions (like cholesky_decompose) with GPU implementations will automatically be used

1 Like

Thank you @andrjohns. That’s great!

1 Like