In my understanding, the Stan function for computing the log density

\log \text{MultiNormal}(y \mid \mu, \Sigma),

i.e., multi_normal_lpdf, requires inversion of the n \times n matrix \Sigma, which is an O(n^3) operation. What about having also a function multi_normal_structured_lpdf or something which could be used if the user knows that

\Sigma = A + U C V

where

A is diagonal n \times n

U is n \times m

C is m \times m

V is m \times n

with m < n? The complexity would then be O(m^2 n) (I think) if the Woodbury formula is used in the underlying implementation to invert \Sigma.

Yes, multi_normal_lpdf requires an inversion (technically, we donâ€™t invert, we use matrix division for more stabilty). If there are multiple y (e.g., an array of vectors), the the solve only happens once.

multi_normal_precision_lpdf does not require an inversion if you can work with precision matrices instead of covariance. Also, multi_normal_cholesky_lpdf works with a Cholesky factor rather than the full matrix, which makes inversion easy. And itâ€™s a much more stable parameterization if the covariance matrix is a parameter. But obviously it wonâ€™t work in something like a GP, where we first build the covariance matrix, then factor it.

That shouldnâ€™t be too challenging. It only requires working out derivatives of inverse(A + U * C * V) to be efficient. As with other additions to Stan, itâ€™s a matter of motivating someone to do it.

When does this structure of covariance matrix come up? Will that structure give us something full rank? As is, we require positive-definite covariance and canâ€™t deal with only positive semi-definite covariance matrices.

I have been doing GPs so that L = cholesky_decompose(Ky) and then y ~ multi_normal_cholesky(mu, L), but I am not sure if it has much advantage over just y ~ multi_normal(mu, Ky).

Yeah, maybe I will give it a shot myself if I have time and no one else has done it by then.

It comes up with GPs if you do a low-rank approximation for the kernel matrix K \approx \tilde{K} = U \Lambda U^T and then you have to compute

\text{MultiNormal}(y \mid \mu, \Sigma)

where \Sigma = \tilde{K} + \sigma^2 I. So \tilde{K} is not full rank but \Sigma is so it is not a problem.