Sampling multivariate normal with low-rank covariance

Question

In the context of sparse Gaussian Processes we get an approximation for an $N$ by $N$ positive definite covariance matrix that is of rank $M$:

$\Sigma = \Lambda + VV^T$

where $\Lambda$ is diagonal and $V$ is an $N$ by $M$ matrix.

This is very useful for finding the inverse (via the Woodbury formula), but what I need is to sample from the resulting multivariate normal. This usually requires a Cholesky or an eigendecomposition, but I cannot find a way to exploit the low-rank structure of $\Sigma$ in those cases.

So, my question is whether there is a way to compute Cholesky, Eigen/SVD, or to sample from the normal distribution without explicitly forming the $N$ by $N$ matrix?

Answer 1

No need to make Cholesky decomposition to sample from a multivariate Gaussian with covariance $\Sigma = \Lambda + V^\top V$.

Let the shapes of $\Lambda = \text{diag}(\lambda)$ be (N, N) and $V$ be $(N, K)$. First, we need samples from i.i.d. zero-mean one-variance normal with size $(N, )$ and $(K, )$. Let these samples be $\epsilon_\lambda$ and $\epsilon_V$.

The sample we want can be computed from

$$\Lambda^{1/2} \epsilon_\lambda + V \epsilon_V.$$

The variance of the first term is $\Lambda$ and that for the second is $V^T V$.