Preconditioned Krylov Subspace Methods for Sampling Multivariate Gaussian Distributions

A common problem in statistics is to compute sample vectors from a multivariate Gaussian distribution with zero mean and a given covariance matrix $A$. A canonical approach to the problem is to compute vectors of the form $y = S z$, where $S$ is the Cholesky factor or square root of $A$, and $z$ is a standard normal vector. When $A$ is large, such an approach becomes computationally expensive. This paper considers preconditioned Krylov subspace methods to perform this task. The Lanczos process provides a means to approximate $A^{1/2} z$ for any vector $z$ from an $m$-dimensional Krylov subspace. The main contribution of this paper is to show how to enhance the convergence of the process via preconditioning. Both incomplete Cholesky preconditioners and approximate inverse preconditioners are discussed. It is argued that the latter class of preconditioners has an advantage in the context of sampling. Numerical tests, performed with stationary covariance matrices used to model Gaussian processes, illustrate ...

[1]  R. A. Silverman,et al.  Special functions and their applications , 1966 .

[2]  M. Shinozuka,et al.  Digital simulation of random processes and its applications , 1972 .

[3]  J. Pasciak,et al.  Computer solution of large sparse positive definite systems , 1982 .

[4]  Donald Geman,et al.  Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images , 1984, IEEE Transactions on Pattern Analysis and Machine Intelligence.

[5]  William F. Moss,et al.  Decay rates for inverses of band matrices , 1984 .

[6]  O. Axelsson,et al.  On the rate of convergence of the preconditioned conjugate gradient method , 1986 .

[7]  M. Fixman.  Construction of Langevin forces in the simulation of hydrodynamic interaction , 1986 .

[8]  H. V. D. Vorst,et al.  The rate of convergence of Conjugate Gradients , 1986 .

[9]  Goodman,et al.  Multigrid Monte Carlo method. Conceptual foundations. , 1989, Physical review. D, Particles and fields.

[10]  Y. Saad,et al.  On acceleration methods for coupled nonlinear elliptic systems , 1991 .

[11]  Leslie Greengard,et al.  The Fast Gauss Transform , 1991, SIAM J. Sci. Comput..

[12]  L. Knizhnerman,et al.  Two polynomial methods of calculating functions of symmetric matrices , 1991 .

[13]  Yousef Saad,et al.  Efficient Solution of Parabolic Equations by Krylov Approximation Methods , 1992, SIAM J. Sci. Comput..

[14]  Y. Saad.  Analysis of some Krylov subspace approximations to the matrix exponential operator , 1992 .

[15]  L. Kolotilina,et al.  Factorized Sparse Approximate Inverse Preconditionings I. Theory , 1993, SIAM J. Matrix Anal. Appl..

[16]  Y. Saad.  Theoretical Error Bounds and General Analysis of afew , 2022 .

[17]  Vladimir Druskin,et al.  Krylov subspace approximation of eigenpairs and matrix functions in exact and computer arithmetic , 1995, Numer. Linear Algebra Appl..

[18]  C. R. Dietrich,et al.  Fast and Exact Simulation of Stationary Gaussian Processes through Circulant Embedding of the Covariance Matrix , 1997, SIAM J. Sci. Comput..

[19]  Michele Benzi,et al.  Robust Approximate Inverse Preconditioning for the Conjugate Gradient Method , 2000, SIAM J. Sci. Comput..

[20]  A. A. Nikishin,et al.  Factorized sparse approximate inverse preconditionings. III. Iterative construction of preconditioners , 2000 .

[21]  Michele Benzi,et al.  Orderings for Factorized Sparse Approximate Inverse Preconditioners , 1999, SIAM J. Sci. Comput..

[22]  Edmond Chow,et al.  A Priori Sparsity Patterns for Parallel Sparse Approximate Inverse Preconditioners , 1999, SIAM J. Sci. Comput..

[23]  H. Rue.  Fast sampling of Gaussian Markov random fields , 2000 .

[24]  H. V. D. Vorst,et al.  Numerical methods for the QCDd overlap operator. I. Sign-function and error bounds , 2002, hep-lat/0202025.

[25]  Alan S. Willsky,et al.  A Krylov Subspace Method for Covariance Approximation and Simulation of Random Processes and Fields , 2003, Multidimens. Syst. Signal Process..

[26]  Yousef Saad,et al.  Iterative methods for sparse linear systems , 2003 .

[27]  Christopher K. I. Williams,et al.  Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) , 2005 .

[28]  Oliver G. Ernst,et al.  A Restarted Krylov Subspace Method for the Evaluation of Matrix Functions , 2006, SIAM J. Numer. Anal..

[29]  Nicholas J. Higham,et al.  Functions of matrices - theory and computation , 2008 .

[30]  A. Pettitt,et al.  Fast sampling from a Gaussian Markov random field using Krylov subspace approaches , 2008 .

[31]  Valeria Simoncini,et al.  Acceleration Techniques for Approximating the Matrix Exponential Operator , 2008, SIAM J. Matrix Anal. Appl..

[32]  Carl E. Rasmussen,et al.  Gaussian processes for machine learning , 2005, Adaptive computation and machine learning.

[33]  I. Turner,et al.  A restarted Lanczos approximation to functions of a symmetric matrix , 2010 .

[34]  Lei Wang,et al.  Fast Evaluation of Multiquadric RBF Sums by a Cartesian Treecode , 2011, SIAM J. Sci. Comput..

[35]  Jie Chen.  A DEFLATED VERSION OF THE BLOCK CONJUGATE GRADIENT ALGORITHM WITH AN APPLICATION TO GAUSSIAN PROCESS MAXIMUM LIKELIHOOD ESTIMATION , 2011 .

[36]  Mihai Anitescu,et al.  Computing f(A)b via Least Squares Polynomial Approximations , 2011, SIAM J. Sci. Comput..

[37]  Colin Fox,et al.  Sampling Gaussian Distributions in Krylov Spaces with Conjugate Gradients , 2012, SIAM J. Sci. Comput..

[38]  J. Skolnick,et al.  Krylov subspace methods for computing hydrodynamic interactions in brownian dynamics simulations. , 2012, The Journal of chemical physics.

[39]  Jo Eidsvik,et al.  Norges Teknisk-naturvitenskapelige Universitet Iterative Numerical Methods for Sampling from High Dimensional Gaussian Distributions Iterative Numerical Methods for Sampling from High Dimensional Gaussian Distributions , 2022 .

[40]  Stochastic Relaxation , 2014, Computer Vision, A Reference Guide.

[41]  Owe Axelsson,et al.  Reaching the superlinear convergence phase of the CG method , 2014, J. Comput. Appl. Math..