Simulating Gaussian processes on a lattice

How can I simulate a Gaussian Process 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 points $\mathrm{\Omega }$. $\mathrm{\Omega }=(x_0,x_1,\dots ,x_m)$.

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

The problem then reduces to generating a vector $\text{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 $\epsilon \sim \mathcal{N}(0,I)$ is an $m+1$-dimensional normal random variable, and $A{A}^T=R$, then $\text{vv}y=\text{mm}A\text{vv}\epsilon$ has the required distribution.

1 The circulant embedding trick

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

Suppose further that our points form a grid, $\mathrm{\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 $\text{vv}{x}^t\text{mm}R\text{vv}x\ge 0\mathrm{\forall }\text{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 Loper et al. (2021), as well as Kronecker and Toeplitz matrices when working with regularly spaced grids (Dietrich and Newsam 1997; Grace Chan and Wood 1997).