# Random samples from a normal distribution without explicitly constructing a covariance matrix

I have a set of samples in MxN matrix X that I want to use to approximate a covariance matrix:

$$C = X^\top X + eI$$

Where e is a small scalar and I is the NxN identity matrix.

I need draw a random sample from the normal distribution with covariance matrix $C^{-1}$. I think this is normally done by drawing a random sample z from the standard normal, finding a factorization of $C^{-1} = M M^\top$ and then calculating $y = M z$.

The difficulty is that N will in general be large so I need to avoid explicitly calculating C or other NxN matricies. This seems difficult to me, but I know that some L-BFGS methods do something a little like this so I have some hope. Any good ways of doing this?

One potentially convenient fact is that the $C^{-1} = {L^{-1}}^\top L^{-1}$ where $C = L L^\top$ is the Cholesky decomposition of C. If I could solve $L^\top y = z$ for y without explicitly constructing L, that would solve my problem.

Try the transformation $z=X^T a+\sqrt{\epsilon}b$ where a and b are standerd normal of dim M and N