Accelerated Gibbs sampling of normal distributions using matrix splittings and polynomials

Standard Gibbs sampling applied to a multivariate normal distribution with a specified precision matrix is equivalent in fundamental ways to the Gauss-Seidel iterative solution of linear equations in the precision matrix. Specifically, the iteration operators, the conditions under which convergence occurs, and geometric convergence factors (and rates) are identical. These results hold for arbitrary matrix splittings from classical iterative methods in numerical linear algebra giving easy access to mature results in that field, including existing convergence results for antithetic-variable Gibbs sampling, REGS sampling, and generalizations. Hence, efficient deterministic stationary relaxation schemes lead to efficient generalizations of Gibbs sampling. The technique of polynomial acceleration that significantly improves the convergence rate of an iterative solver derived from a \emph{symmetric} matrix splitting may be applied to accelerate the equivalent generalized Gibbs sampler. Identicality of error polynomials guarantees convergence of the inhomogeneous Markov chain, while equality of convergence factors ensures that the optimal solver leads to the optimal sampler. Numerical examples are presented, including a Chebyshev accelerated SSOR Gibbs sampler applied to a stylized demonstration of low-level Bayesian image reconstruction in a large 3-dimensional linear inverse problem.

[1]  William Feller,et al.  An Introduction to Probability Theory and Its Applications , 1951 .

[2]  P. Whittle.  ON STATIONARY PROCESSES IN THE PLANE , 1954 .

[3]  J. Gillis,et al.  Matrix Iterative Analysis , 1961 .

[4]  William Feller,et al.  An Introduction to Probability Theory and Its Applications , 1967 .

[5]  L. Fox,et al.  Chebyshev polynomials in numerical analysis , 1970 .

[6]  V. Turchin.  On the Computation of Multidimensional Integrals by the Monte-Carlo Method , 1971 .

[7]  Louis A. Hageman,et al.  Iterative Solution of Large Linear Systems. , 1971 .

[8]  S. Adler.  Over-relaxation method for the Monte Carlo evaluation of the partition function for multiquadratic actions , 1981 .

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

[10]  Charles R. Johnson,et al.  Matrix analysis , 1985, Statistical Inference for Engineers and Data Scientists.

[11]  J. Besag.  On the Statistical Analysis of Dirty Pictures , 1986 .

[12]  E. Loh.  Multigrid Monte Carlo methods , 1988 .

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

[14]  Adrian F. M. Smith,et al.  Sampling-Based Approaches to Calculating Marginal Densities , 1990 .

[15]  P. Barone,et al.  Improving Stochastic Relaxation for Gussian Random Fields , 1990, Probability in the Engineering and Informational Sciences.

[16]  Y. Amit.  On rates of convergence of stochastic relaxation for Gaussian and non-Gaussian distributions , 1991 .

[17]  U. Grenander,et al.  Comparing sweep strategies for stochastic relaxation , 1991 .

[18]  P. Green,et al.  Metropolis Methods, Gaussian Proposals and Antithetic Variables , 1992 .

[19]  J. Besag,et al.  Spatial Statistics and Bayesian Computation , 1993 .

[20]  O. Nevanlinna.  Convergence of Iterations for Linear Equations , 1993 .

[21]  J Besag,et al.  DISCUSSION ON THE MEETING ON THE GIBBS SAMPLER AND OTHER MARKOV CHAIN-MONTE CARLO METHODS , 1993 .

[22]  Jun S. Liu,et al.  Covariance Structure and Convergence Rate of the Gibbs Sampler with Various Scans , 1995 .

[23]  David B. Dunson,et al.  Bayesian Data Analysis , 2010 .

[24]  David J. Spiegelhalter,et al.  Introducing Markov chain Monte Carlo , 1995 .

[25]  Coordinate Selection Rules for Gibbs Sampling , 1996 .

[26]  Stephen S. Wilson,et al.  Random iterative models , 1996 .

[27]  G. Roberts,et al.  Updating Schemes, Correlation Structure, Blocking and Parameterization for the Gibbs Sampler , 1997 .

[28]  Charles M. Grinstead,et al.  Introduction to probability , 1999, Statistics for the Behavioural Sciences.

[29]  Radford M. Neal.  Monte Carlo Implementation of Gaussian Process Models for Bayesian Regression and Classification , 1997, physics/9701026.

[30]  Radford M. Neal,et al.  Suppressing Random Walks in Markov Chain Monte Carlo Using Ordered Overrelaxation , 1995, Learning in Graphical Models.

[31]  S. Paddock,et al.  Confocal laser scanning microscopy. , 1999, BioTechniques.

[32]  Persi Diaconis,et al.  Iterated Random Functions , 1999, SIAM Rev..

[33]  Y. Saad,et al.  Iterative solution of linear systems in the 20th century , 2000 .

[34]  Alain Galli,et al.  Rate of Convergence of the Gibbs Sampler in the Gaussian Case , 2001 .

[35]  H. Haario,et al.  An adaptive Metropolis algorithm , 2001 .

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

[37]  Julian Stander,et al.  Over-relaxation methods and coupled Markov chains for Monte Carlo simulation , 2002, Stat. Comput..

[38]  Sw. Banerjee,et al.  Hierarchical Modeling and Analysis for Spatial Data , 2003 .

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

[40]  David J. C. MacKay,et al.  Information Theory, Inference, and Learning Algorithms , 2004, IEEE Transactions on Information Theory.

[41]  Peter Guttorp,et al.  On the Whittle-Matérn correlation family , 2005 .

[42]  Leonhard Held,et al.  Gaussian Markov Random Fields: Theory and Applications , 2005 .

[43]  David Higdon,et al.  A Primer on Space-Time Modeling from a Bayesian Perspective , 2006 .

[44]  T. Gneiting,et al.  Fast and Exact Simulation of Large Gaussian Lattice Systems in ℝ2: Exploring the Limits , 2006 .

[45]  G. Meurant.  The Lanczos and Conjugate Gradient Algorithms: From Theory to Finite Precision Computations , 2006 .

[46]  I. Klapper,et al.  A Multidimensional Multispecies Continuum Model for Heterogeneous Biofilm Development , 2007, Bulletin of mathematical biology.

[47]  Haluk Beyenal,et al.  Fundamentals of biofilm research , 2007 .

[48]  J. Rosenthal,et al.  Coupling and Ergodicity of Adaptive Markov Chain Monte Carlo Algorithms , 2007, Journal of Applied Probability.

[49]  Jeffrey S. Rosenthal,et al.  Coupling and Ergodicity of Adaptive MCMC , 2007 .

[50]  Colin Fox,et al.  TECHNICAL REPORTS from the ELECTRONICS GROUP at the UNIVERSITY of OTAGO A Conjugate Direction Sampler for Normal Distributions, with a Few Computed Examples , 2008 .

[51]  G. Meurant.  The Lanczos and conjugate gradient algorithms , 2008 .

[52]  Nicolai Schipper Jespersen,et al.  An Introduction to Markov Chain Monte Carlo , 2010 .

[53]  George Casella,et al.  A Short History of Markov Chain Monte Carlo: Subjective Recollections from Incomplete Data , 2008, 0808.2902.

[54]  H. Rue,et al.  An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach , 2011 .

[55]  Haavard Rue,et al.  Think continuous: Markovian Gaussian models in spatial statistics , 2011, 1110.6796.

[56]  Tiangang Cui,et al.  Bayesian calibration of a large‐scale geothermal reservoir model by a new adaptive delayed acceptance Metropolis Hastings algorithm , 2011 .

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

[58]  Johnathan M. Bardsley,et al.  MCMC-Based Image Reconstruction with Uncertainty Quantification , 2012, SIAM J. Sci. Comput..

[59]  Colin Fox,et al.  Convergence in Variance of Chebyshev Accelerated Gibbs Samplers , 2014, SIAM J. Sci. Comput..

[60]  Barbara Kaltenbacher,et al.  Iterative Solution Methods , 2015, Handbook of Mathematical Methods in Imaging.

[61]  Colin Fox,et al.  Efficiency and computability of MCMC with Langevin, Hamiltonian, and other matrix-splitting proposals , 2015, 1501.03150.