Sampling from Gaussian with very large covariance matrix in block form

Question

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?

Answer 1

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.