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

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

  • $\begingroup$ Thanks a lot! Sometimes I forget about all of the nice properties of Gaussians. $\endgroup$ – Sergiy Protsiv Dec 9 '18 at 19:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.