Gaussian Process simulation and circulant embeddings

I might shoehorn Whittle likelihoods in here too

An converse problem to covariance estimation. Related: phase retrieval. Gaussian process regression.


Discuss in the context of

  • deterministic covariance
  • expected autocorrelation in covariance

Simulating Gaussian fields with the desired covariance structure

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\bb \varepsilon\) has the required distribution.

The circulant embedding trick

But what if that’s too slow? If we have additional structure, we can do better.

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?) 🏗

Simulating point processes with the desired covariance structure

For now, see spatial point processes. 🏗


Alexanderian, Alen. 2015. “A Brief Note on the Karhunen-Loève Expansion.” October 26, 2015.
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.” May 29, 2016.
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.
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.
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.” October 2, 2017.
———. 2017b. “Circulant Embedding with QMC – Analysis for Elliptic PDE with Lognormal Coefficients.” October 25, 2017.
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.
Heinig, Georg, and Karla Rost. 2011. “Fast Algorithms for Toeplitz and Hankel Matrices.” Linear Algebra and Its Applications 435 (1): 1–59.
Powell, Catherine E. 2014. “Generating Realisations of Stationary Gaussian Random Fields by Circulant Embedding.” Matrix 2 (2): 1.
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.
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.
Ye, Ke, and Lek-Heng Lim. 2016. “Every Matrix Is a Product of Toeplitz Matrices.” Foundations of Computational Mathematics 16 (3): 577–98.