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?