Simulating Gaussian processes on a lattice

Assumed audience:

ML people

How can I simulate a Gaussian Processes on a lattice with a given covariance?

The general (non-lattice) case is given in historical overview in Liu et al. (2019), but in this notebook we are interested in specialising a little. Following the introduction in Dietrich and Newsam (1993), let’s say we wish to generate a stationary Gaussian process \(Y(x)\) on a points \(\Omega\). \(\Omega=(x_0, x_1,\dots, x_m)\).

Stationary in this context means that the covariance function \(r\) is translation-invariance and depend only on distance, so that it may be given \(r(|x|)\). Without loss of generality, we assume that \(\bb E[Y(x)]=0\) and \(\var[Y(x)]=1\).

The problem then reduces to generating a vector \(\vv y=(Y(x_0), Y(x_1), \dots, Y(x_m) )\sim \mathcal{N}(0, R)\) where \(R\) has entries \(R[p,q]=r(|x_p-x_q|).\)

Note that if \(\bb \varepsilon\sim\mathcal{N}(0, I)\) is an \(m+1\)-dimensional normal random variable, and \(AA^T=R\), then \(\vv y=\mm A \vv \varepsilon\) has the required distribution.

The circulant embedding trick

If we have additional structure, we can work more efficiently.

Suppose further that our points form a grid, \(\Omega=(x_0, x_0+h,\dots, x_0+mh)\); specifically, equally-spaced-points on a line.

We know that \(R\) has a Toeplitz structure. Moreover it is non-negative definite, with \(\vv x^t\mm R \vv x \geq 0\forall \vv x.\) (Why?) 🏗

Wilson et al. (2021) credits the following authors:

Well-known examples of this trend include banded and sparse matrices in the context of one-dimensional Gaussian processes and Gauss–Markov random fields [;Durrande et al. (2019)], as well as Kronecker and Toeplitz matrices when working with regularly-spaced grids (Dietrich and Newsam 1997; Grace Chan and Wood 1997).


Abrahamsen, P., V. Kvernelv, and D. Barker. 2018. Simulation Of Gaussian Random Fields Using The Fast Fourier Transform (Fft).” In, 2018:1–14. European Association of Geoscientists & Engineers.
Alexanderian, Alen. 2015. A Brief Note on the Karhunen-Loève Expansion.” arXiv:1509.07526 [Math], October.
Bolin, David. 2016. Models and Methods for Random Fields in Spatial Statistics with Computational Efficiency from Markov Properties.
Chan, Grace, and Andrew T.A. Wood. 1997. Algorithm AS 312: An Algorithm for Simulating Stationary Gaussian Random Fields.” Journal of the Royal Statistical Society: Series C (Applied Statistics) 46 (1): 171–81.
Chan, G., and A. T. A. Wood. 1999. Simulation of Stationary Gaussian Vector Fields.” Statistics and Computing 9 (4): 265–68.
Choromanski, Krzysztof, and Vikas Sindhwani. 2016. Recycling Randomness with Structure for Sublinear Time Kernel Expansions.” arXiv:1605.09049 [Cs, Stat], May.
Cotter, S. L., G. O. Roberts, A. M. Stuart, and D. White. 2013. MCMC Methods for Functions: Modifying Old Algorithms to Make Them Faster.” Statistical Science 28 (3): 424–46.
Davies, Tilman M., and David Bryant. 2013. On Circulant Embedding for Gaussian Random Fields in R.” Journal of Statistical Software 55 (9).
Dietrich, C. R., and G. N. Newsam. 1993. A Fast and Exact Method for Multidimensional Gaussian Stochastic Simulations.” Water Resources Research 29 (8): 2861–69.
———. 1997. Fast and Exact Simulation of Stationary Gaussian Processes Through Circulant Embedding of the Covariance Matrix.” SIAM Journal on Scientific Computing 18 (4): 1088–1107.
Durrande, Nicolas, Vincent Adam, Lucas Bordeaux, Stefanos Eleftheriadis, and James Hensman. 2019. Banded Matrix Operators for Gaussian Markov Models in the Automatic Differentiation Era.” In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, 2780–89. PMLR.
Ellis, Robert L., and David C. Lay. 1992. Factorization of Finite Rank Hankel and Toeplitz Matrices.” Linear Algebra and Its Applications 173 (August): 19–38.
Erhel, Jocelyne, Mestapha Oumouni, Géraldine Pichot, and Franck Schoefs. n.d. “Analysis of Continuous Spectral Method for Sampling Stationary Gaussian Random Fields,” 26.
Galy-Fajou, Théo, Valerio Perrone, and Manfred Opper. 2021. Flexible and Efficient Inference with Particles for the Variational Gaussian Approximation.” Entropy 23 (8): 990.
Gilboa, E., Y. Saatçi, and J. P. Cunningham. 2015. Scaling Multidimensional Inference for Structured Gaussian Processes.” IEEE Transactions on Pattern Analysis and Machine Intelligence 37 (2): 424–36.
Graham, Ivan G., Frances Y. Kuo, Dirk Nuyens, Rob Scheichl, and Ian H. Sloan. 2017a. Analysis of Circulant Embedding Methods for Sampling Stationary Random Fields.” arXiv:1710.00751 [Math], October.
———. 2017b. Circulant Embedding with QMC – Analysis for Elliptic PDE with Lognormal Coefficients.” arXiv:1710.09254 [Math], October.
Gray, Robert M. 2006. Toeplitz and Circulant Matrices: A Review. Vol. 2.
Guinness, Joseph, and Montserrat Fuentes. 2016. Circulant Embedding of Approximate Covariances for Inference From Gaussian Data on Large Lattices.” Journal of Computational and Graphical Statistics 26 (1): 88–97.
Haran, Murali. 2011. Gaussian Random Field Models for Spatial Data.” In Handbook of Markov Chain Monte Carlo, edited by Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. Vol. 20116022. Chapman and Hall/CRC.
Heinig, Georg, and Karla Rost. 2011. Fast Algorithms for Toeplitz and Hankel Matrices.” Linear Algebra and Its Applications 435 (1): 1–59.
Lang, Annika, and Jürgen Potthoff. 2011. Fast Simulation of Gaussian Random Fields.” Monte Carlo Methods and Applications 17 (3).
Latz, Jonas, Marvin Eisenberger, and Elisabeth Ullmann. 2019. Fast Sampling of Parameterised Gaussian Random Fields.” Computer Methods in Applied Mechanics and Engineering 348 (May): 978–1012.
Liu, Yang, Jingfa Li, Shuyu Sun, and Bo Yu. 2019. Advances in Gaussian Random Field Generation: A Review.” Computational Geosciences 23 (5): 1011–47.
Powell, Catherine E. 2014. “Generating Realisations of Stationary Gaussian Random Fields by Circulant Embedding.” Matrix 2 (2): 1.
Rue, Havard. 2001. Fast Sampling of Gaussian Markov Random Fields.” Journal of the Royal Statistical Society. Series B (Statistical Methodology) 63 (2): 325–38.
Rue, Håvard, and Leonhard Held. 2005. Gaussian Markov Random Fields: Theory and Applications. Monographs on Statistics and Applied Probability 104. Boca Raton: Chapman & Hall/CRC.
Sidén, Per. 2020. Scalable Bayesian Spatial Analysis with Gaussian Markov Random Fields. Vol. 15. Linköping Studies in Statistics. Linköping: Linköping University Electronic Press.
Stroud, Jonathan R., Michael L. Stein, and Shaun Lysen. 2017. Bayesian and Maximum Likelihood Estimation for Gaussian Processes on an Incomplete Lattice.” Journal of Computational and Graphical Statistics 26 (1): 108–20.
Teichmann, Jakob, and Karl-Gerald van den Boogaart. 2016. Efficient Simulation of Stationary Multivariate Gaussian Random Fields with Given Cross-Covariance.” Applied Mathematics 7 (17): 2183–94.
Whittle, P. 1954. “On Stationary Processes in the Plane.” Biometrika 41 (3/4): 434–49.
Whittle, P. 1953a. The Analysis of Multiple Stationary Time Series.” Journal of the Royal Statistical Society: Series B (Methodological) 15 (1): 125–39.
———. 1953b. Estimation and Information in Stationary Time Series.” Arkiv För Matematik 2 (5): 423–34.
Whittle, Peter. 1952. Some Results in Time Series Analysis.” Scandinavian Actuarial Journal 1952 (1-2): 48–60.
Wilson, James T, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter Deisenroth. 2021. Pathwise Conditioning of Gaussian Processes.” Journal of Machine Learning Research 22 (105): 1–47.
Yang, Linxiao, Jun Fang, Huiping Duan, Hongbin Li, and Bing Zeng. 2018. Fast Low-Rank Bayesian Matrix Completion with Hierarchical Gaussian Prior Models.” IEEE Transactions on Signal Processing 66 (11): 2804–17.
Ye, Ke, and Lek-Heng Lim. 2016. Every Matrix Is a Product of Toeplitz Matrices.” Foundations of Computational Mathematics 16 (3): 577–98.

No comments yet. Why not leave one?

GitHub-flavored Markdown & a sane subset of HTML is supported.