Gaussian process approximations

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 ${\displaystyle d}$ over ${\displaystyle {\mathcal {X}}}$ and assumes that for ${\displaystyle X,{\tilde {X}}\in {\mathcal {X}}}$ we have ${\displaystyle Cov(y(X),y({\tilde {X}}))\neq 0}$ only if ${\displaystyle d(X,{\tilde {X}})<r}$ for some radius ${\displaystyle r}$. The second method assumes that there exist ${\displaystyle {\mathcal {X}}^{(1)},\dots ,{\mathcal {X}}^{(K)}}$ such that ${\displaystyle \bigcup _{k=1}^{K}{\mathcal {X}}^{(k)}}$. Then with appropriate distribution of indices among partition elements and ordering of elements of ${\displaystyle X}$ the covariance matrix is block diagonal.

This family of methods assumes that the precision matrix ${\displaystyle \mathbf {\Lambda } =\mathbf {\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 ${\displaystyle d=2}$ and when ${\displaystyle {\mathcal {X}}}$ has a defined metric and takes advantage of the fact, that the Markov property holds which makes ${\displaystyle \mathbf {\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 ${\displaystyle {\mathcal {X}}}$ and takes advantage of the so-called screening effect assuming that ${\displaystyle \mathbf {\Lambda } _{i,j}\neq 0}$ only if ${\displaystyle d(x_{i},x_{j})<r}$, for some ${\displaystyle r>0}$.

In many practical applications, calculating ${\displaystyle \mathbf {\Lambda } }$ is replaced with computing first ${\displaystyle \mathbf {L} }$, the Cholesky factor of ${\displaystyle \mathbf {\Sigma } }$, and second its inverse ${\displaystyle \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 ${\displaystyle \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 ${\displaystyle y}$, the Gaussian process of interest, is effectively low-rank. More precisely, it is assumed, that there exists a set of indices ${\displaystyle {\bar {X}}=\{{\bar {x}}_{1},\dots ,{\bar {x}}_{p}\}}$ such that every other set of indices ${\displaystyle X=\{x_{1},\dots ,x_{n}\}}$

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

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

More generally, on top of selecting ${\displaystyle {\bar {X}}}$, one may also find an ${\displaystyle n\times p}$ matrix ${\displaystyle \mathbf {A} }$ and assume that ${\displaystyle X=\mathbf {A} \mathbf {\eta } }$, where ${\displaystyle \mathbf {\eta } }$ are ${\displaystyle p}$ values of a Gaussian process possibly independent of ${\displaystyle 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 ${\displaystyle \mathbf {A} }$ and ${\displaystyle \mathbf {\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 ${\displaystyle 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 (${\displaystyle {\mathcal {O}}(n\log n)}$) complexity.

Probabilistic graphical models provide a convenient framework for comparing model-based approximations. In this context, value of the process at index ${\displaystyle x_{k}\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 ${\displaystyle 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.