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.