Gaussian process approximations

In statistical modeling, it is often convenient to assume that $y\; \backslash in\; \backslash mathcal\{Y\}$, the phenomenon under investigation is a Gaussian process indexed by $X\; \backslash in\; \backslash mathcal\{X\}\; =\; \backslash mathcal\{X\}\_1\; \backslash times\; \backslash mathcal\{X\}\_2\; \backslash dots\; \backslash mathcal\{X\}\_d$ which has mean function $\backslash mu:\; \backslash mathcal\{X\}\; \backslash rightarrow\; \backslash mathcal\{Y\}$ and covariance function $K:\; \backslash mathcal\{X\}\; \backslash times\; \backslash mathcal\{X\}\; \backslash rightarrow\; \backslash mathbb\{R\}$.

One can also assume that data $\backslash mathbf\{y\}\; =\; (y\_1,\; \backslash dots,\; y\_n)$ are values of a particular realization of this process for indices $\backslash mathbf\{X\}\; =\; X\_1,\; \backslash dots,\; X\_n$.

Consequently, the joint distribution of the data can be expressed as

:$$

\mathbf{y} \sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma}),

where $\backslash mathbf\{\backslash Sigma\}\; =\; \backslash left[\; K(X\_i,\; X\_j)\; \backslash right]\_\{i,j=1\}^n$ and $\backslash mathbf\{\backslash mu\}\; =\; \backslash left(\backslash mu(X\_1),\; \backslash mu(X\_2),\; \backslash dots,\; \backslash mu(X\_d)\backslash right)^\{\backslash top\}$, i.e. respectively a matrix with the covariance function values and a vector with the mean function values at corresponding (pairs of) indices.

The negative log-likelihood of the data then takes the form

:$$

-\log\ell(\mathbf{y}) = \frac{d}{2\pi} + \frac{1}{2}\log\det(\mathbf{\Sigma}) + \left(\mathbf{y}-\mathbf{\mu}\right)^\top\mathbf{\Sigma}^{-1}\left(\mathbf{y}-\mathbf{\mu}\right)

Similarly, the best predictor of $\backslash mathbf\{y\}^*$, the values of $y$ for indices $\backslash mathbf\{X\}^*\; =\; \backslash left(X\_1^*,\; X\_2^*,\; \backslash dots,\; X\_d^*\backslash right)$, given data $\backslash mathbf\{y\}$ has the form

:$$

\mathbf{\mu}^*_\mathbf{y} = \mathbb{E}\left[\mathbf{y}^*|\mathbf{y}\right] = \mathbf{\mu}^* - \mathbf{\Sigma}_{\mathbf{y}^*\mathbf{y}} \mathbf{\Sigma}^{-1}\left(\mathbf{y} - \mathbf{\mu}\right)

In the context of Gaussian models, especially in geostatistics, prediction using the best predictor, i.e. mean conditional on the data, is also known as kriging.

The most computationally expensive component of the best predictor formula is inverting the covariance matrix $\backslash mathbf\{\backslash Sigma\}$, which has cubic complexity $\backslash mathcal\{O\}(n^3)$. Similarly, evaluating likelihood involves both calculating $\backslash mathbf\{\backslash Sigma\}^\{-1\}$ and the determinant $\backslash det(\backslash mathbf\{\backslash Sigma\})$ which has the same cubic complexity.

Gaussian process approximations can often be expressed in terms of assumptions on $y$ under which $\backslash log\backslash ell(\backslash mathbf\{y\})$ and $\backslash mathbf\{\backslash mu\}^*\_\backslash mathbf\{y\}$ can be calculated with much lower complexity. Since these assumptions are generally not believed to reflect reality, the likelihood and the best predictor obtained in this way are not exact, but they are meant to be close to their original values.

This class of approximations is expressed through a set of assumptions which are imposed on the original process and which, typically, imply some special structure of the covariance matrix. Although most of these methods were developed independently, most of them can be expressed as special cases of the sparse general Vecchia approximation.

These methods approximate the true model in a way the covariance matrix is sparse. Typically, each method proposes its own algorithm that takes the full advantage of the sparsity pattern in the covariance matrix. Two prominent members of this class of approaches are covariance tapering and domain partitioning. The first method generally requires a metric $d$ over $\backslash mathcal\{X\}$ and assumes that for $X,\; \backslash tilde\{X\}\; \backslash in\; \backslash mathcal\{X\}$ we have $Cov(y(X),\; y(\backslash tilde\{X\}))\backslash neq\; 0$ only if

This family of methods assumes that the precision matrix $\backslash mathbf\{\backslash Lambda\}\; =\; \backslash mathbf\{\backslash Sigma\}^\{-1\}$ is sparse and generally specifies which of its elements are non-zero. This leads to fast inversion because only those elements need to be calculated. Some of the prominent approximations in this category include the approach based on the equivalence between Gaussian processes with Matern covariance function and stochastic PDEs, periodic embedding, and Nearest Neighbour Gaussian processes. The first method applies to the case of $d=2$ and when $\backslash mathcal\{X\}$ has a defined metric and takes advantage of the fact, that the Markov property holds which makes $\backslash mathbf\{\backslash Lambda\}$ very sparse. The second extends the domain and uses Discrete Fourier Transform to decorrelate the data, which results in a diagonal precision matrix. The third one requires a metric on $\backslash mathcal\{X\}$ and takes advantage of the so-called screening effect assuming that $\backslash mathbf\{\backslash Lambda\}\_\{i,j\}\; \backslash neq\; 0$ only if , for some $r>0$.

In many practical applications, calculating $\backslash mathbf\{\backslash Lambda\}$ is replaced with computing first $\backslash mathbf\{L\}$, the Cholesky factor of $\backslash mathbf\{\backslash Sigma\}$, and second its inverse $\backslash mathbf\{L\}^\{-1\}$. This is known to be more stable than a plain inversion. For this reason, some authors focus on constructing a sparse approximation of the Cholesky factor of the precision or covariance matrices. One of the most established methods in this class is the Vecchia approximation and its generalization. These approaches determine the optimal ordering of indices and, consequently, the elements of $\backslash mathbf\{x\}$ and then assume a dependency structure which minimizes in-fill in the Cholesky factor. Several other methods can be expressed in this framework, the Multi-resolution Approximation (MRA), Nearest Neighbour Gaussian Process, Modified Predictive Process and Full-scale approximation.

While this approach encompasses many methods, the common assumption underlying them all is the assumption, that $y$, the Gaussian process of interest, is effectively low-rank. More precisely, it is assumed, that there exists a set of indices $\backslash bar\{X\}\; =\; \backslash \{\backslash bar\{x\}\_1,\; \backslash dots,\; \backslash bar\{x\}\_p\backslash \}$ such that every other set of indices $X\; =\; \backslash \{x\_1,\; \backslash dots,\; x\_n\backslash \}$

$y(X)\; \backslash sim\; \backslash mathcal\{N\}\backslash left(\backslash mathbf\{A\}\_X\backslash bar\{\backslash mathbf\{\backslash mu\}\},\; \backslash mathbf\{A\}\_X^\{\backslash top\}\backslash bar\{\backslash mathbf\{\backslash Sigma\}\}\backslash mathbf\{A\}\_X\; +\; \backslash mathbf\{D\}\backslash right)$

where $\backslash mathbf\{A\}\_X$ is an $p\; \backslash times\; k$ matrix, $\backslash bar\{\backslash mathbf\{\backslash mu\}\}\; =\; \backslash mu\backslash left(y\backslash left(\backslash bar\{X\}\backslash right)\backslash right)$ and $\backslash bar\{\backslash mathbf\{\backslash Sigma\}\}\; =\; K\backslash left(\backslash bar\{X\},\; \backslash bar\{X\}\backslash right)$ and $\backslash mathbf\{D\}$ is a diagonal matrix. Depending on the method and application various ways of selecting $\backslash bar\{X\}$ have been proposed. Typically, $p$ is selected to be much smaller than $n$ which means that the computational cost of inverting $\backslash bar\{\backslash mathbf\{\backslash Sigma\}\}$ is manageable ($\backslash mathcal\{O\}(p^3)$ instead of $\backslash mathcal\{O\}(n^3)$).

More generally, on top of selecting $\backslash bar\{X\}$, one may also find an $n\; \backslash times\; p$ matrix $\backslash mathbf\{A\}$ and assume that $X\; =\; \backslash mathbf\{A\}\backslash mathbf\{\backslash eta\}$, where $\backslash mathbf\{\backslash eta\}$ are $p$ values of a Gaussian process possibly independent of $x$. Many machine learning methods fall into this category, such as subset-of-regressors (SoR), relevance vector machine, sparse spectrum Gaussian Process and others and they generally differ in the way they derive $\backslash mathbf\{A\}$ and $\backslash mathbf\{\backslash eta\}$.

The general principle of hierarchical approximations consists of a repeated application of some other method, such that each consecutive application refines the quality of the approximation. Even though they can be expressed as a set of statistical assumptions, they are often described in terms of a hierarchical matrix approximation (HODLR) or basis function expansion (LatticeKrig, MRA, wavelets). The hierarchical matrix approach can often be represented as a repeated application of a low-rank approximation to successively smaller subsets of the index set $X$. Basis function expansion relies on using functions with compact support. These features can then be exploited by an algorithm who steps through consecutive layers of the approximation. In the most favourable settings some of these methods can achieve quasi-linear ($\backslash mathcal\{O\}(n\backslash log\; n)$) complexity.

Probabilistic graphical models provide a convenient framework for comparing model-based approximations. In this context, value of the process at index $x\_k\; \backslash in\; X$ can then be represented by a vertex in a directed graph and edges correspond to the terms in the factorization of the joint density of $y(X)$. In general, when no independent relations are assumed, the joint probability distribution can be represented by an arbitrary directed acyclic graph. Using a particular approximation can then be expressed as a certain way of ordering the vertices and adding or removing specific edges.

This class of methods does not specify a statistical model or impose assumptions on an existing one. Three major members of this group are the meta-kriging algorithm, the gapfill algorithm and Local Approximate Gaussian Process approach. The first one partitions the set of indices into $K$ components $\backslash mathcal\{X\}^\{(1)\},\; \backslash dots,\; \backslash mathcal\{X\}^\{(k)\}$, calculates the conditional distribution for each those components separately and then uses geometric median of the conditional PDFs to combine them. The second is based on quantile regression using values of the process which are close to the value one is trying to predict, where distance is measured in terms of a metric on the set of indices. Local Approximate Gaussian Process uses a similar logic but constructs a valid stochastic process based on these neighboring values.

• {{cite journal|arxiv=1807.01065|last1=Liu|first1=Haitao|last2=Ong|first2=Yew-Soon|last3=Shen|first3=Xiaobo|last4=Cai|first4=Jianfei|title=When Gaussian Process Meets Big Data: A Review of Scalable GPS|journal=IEEE Transactions on Neural Networks and Learning Systems|year=2020|volume=PP|issue=11 |pages=1–19|doi=10.1109/TNNLS.2019.2957109|pmid=31944966}}
• {{cite journal|last1=Heaton|first1=Matthew J.|last2=Datta|first2=Abhirup|last3=Finley|first3=Andrew O.|last4=Furrer|first4=Reinhard|last5=Guinness|first5=Joseph|last6=Guhaniyogi|first6=Rajarshi|last7=Gerber|first7=Florian|last8=Gramacy|first8=Robert B.|last9=Hammerling|first9=Dorit|last10=Katzfuss|first10=Matthias|last11=Lindgren|first11=Finn|last12=Nychka|first12=Douglas W.|last13=Sun|first13=Furong|last14=Zammit-Mangion|first14=Andrew|title=A Case Study Competition Among Methods for Analyzing Large Spatial Data|journal=Journal of Agricultural, Biological and Environmental Statistics|volume=24|issue=3|year=2018|pages=398–425|issn=1085-7117|doi=10.1007/s13253-018-00348-w|doi-access=free|pmid=31496633 |pmc=6709111}}
• {{cite journal|last1=Banerjee|first1=Sudipto|title=High-Dimensional Bayesian Geostatistics|journal=Bayesian Analysis|volume=12|issue=2|year=2017|pages=583–614|doi=10.1214/17-BA1056R|pmid=29391920|doi-access=free|pmc=5790125}}

Category:Geostatistics

Category:Computational science

Category:Computational statistics