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
4 Likes
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? https://arxiv.org/pdf/1907.01063.pdf
OR is there another option/tutorial?
Regards
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: https://github.com/bstatcomp/stan_gpu_install_docs/blob/master/install_docs/INSTALL_CMDSTANR.md
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