Vecchia approximation

Let ${\displaystyle x}$ be a Gaussian process indexed by ${\displaystyle {\mathcal {S}}}$ with mean function ${\displaystyle \mu }$ and covariance function ${\displaystyle K}$. Assume that ${\displaystyle S=\{s_{1},\dots ,s_{n}\}\subset {\mathcal {S}}}$ is a finite subset of ${\displaystyle {\mathcal {S}}}$ and ${\displaystyle \mathbf {x} =(x_{1},\dots ,x_{n})}$ is a vector of values of ${\displaystyle x}$ evaluated at ${\displaystyle S}$, i.e. ${\displaystyle x_{i}=x(s_{i})}$ for ${\displaystyle i=1,\dots ,n}$. Assume further, that one observes ${\displaystyle \mathbf {y} =(y_{1},\dots ,y_{n})}$ where ${\displaystyle y_{i}=x_{i}+\varepsilon _{i}}$ with ${\displaystyle \varepsilon _{i}{\overset {\text{i.i.d.}}{\sim }}{\mathcal {N}}(0,\sigma ^{2})}$. In this context the two most common inference tasks include evaluating the likelihood

${\displaystyle {\mathcal {L}}(\mathbf {y} )=\int f(\mathbf {y} ,\mathbf {x} )\,d\mathbf {x} ,}$

or making predictions of values of ${\displaystyle x}$ for ${\displaystyle s^{*}\in {\mathcal {S}}}$ and ${\displaystyle s\not \in S}$, i.e. calculating

${\displaystyle f(x(s^{*})\mid y_{1},\dots ,y_{n}).}$

The original Vecchia method starts with the observation that the joint density of observations ${\displaystyle f(\mathbf {y} )=\left(y_{1},\dots ,y_{n}\right)}$ can be written as a product of conditional distributions

${\displaystyle f(\mathbf {y} )=f(y_{1})\prod _{i=2}^{n}f(y_{i}\mid y_{i-1},\dots ,y_{1}).}$

Vecchia approximation assumes instead that for some ${\displaystyle k\ll n}$

${\displaystyle {\hat {f}}(\mathbf {y} )=f(y_{1})\prod _{i=2}^{n}f(y_{i}\mid y_{i-1},\dots ,y_{\max(i-k,1)}).}$

Vecchia also suggested that the above approximation be applied to observations that are reordered lexicographically using their spatial coordinates. While his simple method has many weaknesses, it reduced the computational complexity to ${\displaystyle {\mathcal {O}}(nk^{3})}$. Many of its deficiencies were addressed by the subsequent generalizations.

To describe extensions of the Vecchia method in its most general form, define ${\displaystyle z_{i}=(x_{i},y_{i})}$ and notice that for ${\displaystyle \mathbf {z} =(z_{1},\dots ,z_{n})}$ it holds that like in the previous section

${\displaystyle f(\mathbf {z} )=f(x_{1},y_{1})\left(\prod _{i=2}^{n}f(x_{i}\mid z_{1:i-1})\right)\left(\prod _{i=2}^{n}f(y_{i}\mid x_{i})\right)}$

because given ${\displaystyle x_{i}}$ all other variables are independent of ${\displaystyle y_{i}}$.

It has been widely noted that the original lexicographic ordering based on coordinates when ${\displaystyle {\mathcal {S}}}$ is two-dimensional produces poor results.[6] More recently another orderings have been proposed, some of which ensure that points are ordered in a quasi-random fashion. Highly scalable, they have been shown to also drastically improve accuracy.[4]

Similar to the basic version described above, for a given ordering a general Vecchia approximation can be defined as

${\displaystyle {\hat {f}}(\mathbf {z} )=f(x_{1},y_{1})\left(\prod _{i=2}^{n}f(x_{i}\mid z_{q(i)})\right)\left(\prod _{i=2}^{n}f(y_{i}\mid x_{i})\right),}$

where ${\displaystyle q(i)\subset \left\{1,\dots ,i-1\right\}}$. Since ${\displaystyle y_{i}\perp x_{-i},y_{-i}\mid x_{i}}$ it follows that ${\displaystyle f(x_{i}\mid z_{q(i)})=f(x_{i}\mid x_{q}(i),y_{q}(i))=f(x_{i}\mid x_{q}(i))}$ since suggesting that the terms ${\displaystyle f(x_{i}\mid z_{q(i)})}$ be replaced with ${\displaystyle f(x_{i}\mid x_{q(i)})}$. It turns out, however, that sometimes conditioning on some of the observations ${\displaystyle z_{i}}$ increases sparsity of the Cholesky factor of the precision matrix of ${\displaystyle (\mathbf {x} ,\mathbf {y} )}$. Therefore, one might instead consider sets ${\displaystyle q_{y}(i)}$ and ${\displaystyle q_{x}(i)}$ such that ${\displaystyle q(i)=q_{y}(i)\cup q_{x}(i)}$ and express ${\displaystyle {\hat {f}}}$ as

${\displaystyle {\hat {f}}(\mathbf {z} )=f(x_{1},y_{1})\left(\prod _{i=2}^{n}f(x_{i}\mid x_{q_{x}(i)},y_{q_{y}(i)})\right)\left(\prod _{i=2}^{n}f(y_{i}\mid x_{i})\right).}$

Multiple methods of choosing ${\displaystyle q_{y}(i)}$ and ${\displaystyle q_{x}(i)}$ have been proposed, most notably the nearest-neighbour Gaussian process (NNGP),[7] meshed Gaussian process[8] and multi-resolution approximation (MRA) approaches using ${\displaystyle q(i)=q_{x}(i)}$, standard Vecchia using ${\displaystyle q(i)=q_{y}(i)}$ and Sparse General Vecchia where both ${\displaystyle q_{y}(i)}$ and ${\displaystyle q_{x}(i)}$ are non-empty.[5]