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.

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.

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