Sampling Gaussian Distributions in Krylov Spaces with Conjugate Gradients

This paper introduces a conjugate gradient sampler that is a simple extension of the method of conjugate gradients (CG) for solving linear systems. The CG sampler iteratively generates samples from a Gaussian probability density, using either a symmetric positive definite covariance or precision matrix, whichever is more convenient to model. Similar to how the Lanczos method solves an eigenvalue problem, the CG sampler approximates the covariance or precision matrix in a small dimensional Krylov space. As with any iterative method, the CG sampler is efficient for high dimensional problems where forming the covariance or precision matrix is impractical, but operating by the matrix is feasible. In exact arithmetic, the sampler generates Gaussian samples with a realized covariance that converges to the covariance of interest. In finite precision, the sampler produces a Gaussian sample with a realized covariance that is the best approximation to the desired covariance in the smaller dimensional Krylov space. ...