Sampling multivariate normal with low-rank covariance

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?

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$$.
$$\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$$.