# Sampling from Gaussian with very large covariance matrix in block form

I'm interested in sampling from a Gaussian with zero-mean and covariance given by: $$\Sigma = \begin{bmatrix} \Sigma_{11} & \Sigma_{12} & \cdots &\Sigma_{1,100}\\ \Sigma_{21} & \Sigma_{22} & \cdots &\Sigma_{2,100}\\ \vdots & \vdots & \cdots &\vdots \\ \Sigma_{100,1} & \Sigma_{100,2} & \cdots &\Sigma_{100,100} \end{bmatrix}$$ where $$\Sigma_{ij}$$ is a square matrix of dimension $$p \times p$$. In words, the covariance matrix is massive and I am not able to load the entire array into memory. Are there any approaches that allow me to sample from such a Gaussian but only access a subset of the blocks at a time?

• There is the standard method but that just begs the question: how to do the Cholesky decomposition. Will you allow yourself to resort to approximations such as PCA? Commented Mar 13 at 18:57
• @BenjaminWang yes, i didn't mention it in my post but if there is a way to compute the cholesky here that would be sufficient to solve my problem. I think PCA is worth a shot, but also not really sure how to do PCA here either.. Commented Mar 13 at 19:08
• If the covariance is low rank, you can use a randomized svd or similar algorithm to find a low rank decomposition.
– p.s.
Commented Mar 26 at 0:59

You can use the Cholesky–Banachiewicz and Cholesky–Crout algorithms. They allow you to calculate each row of the Cholesky decomposition while using only one row from $$\Sigma$$ at a time.