I have a problem where I can directly compute the (sparse) precision matrix (inverse of the covariance) of a multivariate normal distribution, but the covariance itself is not sparse and I don't want to invert things. I would like to sample points using the precision matrix. What's a fast way to do this?
I do know that the standard procedure is to compute the cholesky decomposition of the covariance, that is find lower triangular matrix L such that $LL' = \Sigma$ and then use a univariate generator to compute a vector of iid normal points $u$ so that finally the vector $z = \mu + Lu$ has the correct covariance. But that would require me computing the covariance AND then taking its Cholesky decomposition. Is there anything faster, using the fact that I have the precision matrix?