Gaussian Process simulation and circulant embeddings

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

TODO

Discuss in the context of

• deterministic covariance
• expected autocorrelation in covariance

Simulating Gaussian RVs with the desired covariance structure

[Following the introduction in (DiNe93)]

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

Refs

• DiNe93: C. R. Dietrich, G. N. Newsam (1993) A fast and exact method for multidimensional gaussian stochastic simulations. Water Resources Research, 29(8), 2861–2869. DOI
• GKNS17a: Ivan G. Graham, Frances Y. Kuo, Dirk Nuyens, Rob Scheichl, Ian H. Sloan (2017a) Analysis of circulant embedding methods for sampling stationary random fields. ArXiv:1710.00751 [Math].
• StSL17: Jonathan R. Stroud, Michael L. Stein, Shaun Lysen (2017) Bayesian and Maximum Likelihood Estimation for Gaussian Processes on an Incomplete Lattice. Journal of Computational and Graphical Statistics, 26(1), 108–120. DOI
• GuFu16: Joseph Guinness, 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. DOI
• GKNS17b: Ivan G. Graham, Frances Y. Kuo, Dirk Nuyens, Rob Scheichl, Ian H. Sloan (2017b) Circulant embedding with QMC -- analysis for elliptic PDE with lognormal coefficients. ArXiv:1710.09254 [Math].
• Whit53a: P. Whittle (1953a) Estimation and information in stationary time series. Arkiv För Matematik, 2(5), 423–434. DOI
• YeLi16: Ke Ye, Lek-Heng Lim (2016) Every Matrix is a Product of Toeplitz Matrices. Foundations of Computational Mathematics, 16(3), 577–598. DOI
• ElLa92: Robert L. Ellis, David C. Lay (1992) Factorization of finite rank Hankel and Toeplitz matrices. Linear Algebra and Its Applications, 173, 19–38. DOI
• HeRo11: Georg Heinig, Karla Rost (2011) Fast algorithms for Toeplitz and Hankel matrices. Linear Algebra and Its Applications, 435(1), 1–59. DOI
• Powe14: Catherine E. Powell (2014) Generating realisations of stationary gaussian random fields by circulant embedding. Matrix, 2(2), 1.
• DaBr13: Tilman M. Davies, David Bryant (2013) On Circulant Embedding for Gaussian Random Fields in R. Journal of Statistical Software, 55(9). DOI
• ChSi16: Krzysztof Choromanski, Vikas Sindhwani (2016) Recycling Randomness with Structure for Sublinear time Kernel Expansions. ArXiv:1605.09049 [Cs, Stat].
• ChWo99: G. Chan, A. T. A. Wood (1999) Simulation of stationary Gaussian vector fields. Statistics and Computing, 9(4), 265–268. DOI
• Whit52: Peter Whittle (1952) Some results in time series analysis. Scandinavian Actuarial Journal, 1952(1–2), 48–60. DOI
• Whit53b: P. Whittle (1953b) The Analysis of Multiple Stationary Time Series. Journal of the Royal Statistical Society. Series B (Methodological) , 15(1), 125–139.
• Gray06: Robert M. Gray (2006) Toeplitz and Circulant Matrices: A Review. Foundations and Trends® in Communications and Information Theory, 2(3), 155–239. DOI