Vecchia Approximations of Gaussian-Process Predictions

Abstract

Gaussian processes (GPs) are highly flexible function estimators used for geospatial analysis, nonparametric regression, and machine learning, but they are computationally infeasible for large datasets. Vecchia approximations of GPs have been used to enable fast evaluation of the likelihood for parameter inference. Here, we study Vecchia approximations of spatial predictions at observed and unobserved locations, including obtaining joint predictive distributions at large sets of locations. We consider a general Vecchia framework for GP predictions, which contains some novel and some existing special cases. We study the accuracy and computational properties of these approaches theoretically and numerically, proving that our new methods exhibit linear computational complexity in the total number of spatial locations. We show that certain choices within the framework can have a strong effect on uncertainty quantification and computational cost, which leads to specific recommendations on which methods are most suitable for various settings. We also apply our methods to a satellite dataset of chlorophyll fluorescence, showing that the new methods are faster or more accurate than existing methods and reduce unrealistic artifacts in prediction maps. Supplementary materials accompanying this paper appear on-line.

Appendices

Vector and Indexing Notation

As an example, define $\mathbf{y} = (\mathbf{y} _1,\mathbf{y} _2,\mathbf{y} _3,\mathbf{y} _4,\mathbf{y} _5)$ as a vector of vectors. Vectors of indices are used for defining subvectors. For example, if $o = (4,1,2)$ is a vector of indices, then $\mathbf{y} _o = (\mathbf{y} _4,\mathbf{y} _1,\mathbf{y} _2)$. Unions of vectors are vectors and are defined when the two vectors have the same type and when the ordering is specified. For example, if $\mathbf{z} = (\mathbf{z} _1,\mathbf{z} _2)$, then $\mathbf{y} _o \cup \mathbf{z} = (\mathbf{y} _4,\mathbf{y} _1, \mathbf{z} _1, \mathbf{y} _2, \mathbf{z} _2)$ defines the union of $\mathbf{y} _o$ and $\mathbf{z}$. Intersection is defined similarly and uses the $\cap$ notation. The ordering of the elements of the union or intersection can be defined alternatively via an index function $\#$ taking in an element and a vector and returning the index occupied by the element in the vector. Continuing the example above, $\#(\mathbf{y} _4, \mathbf{y} ) = 4$, whereas $\#(\mathbf{y} _4, \mathbf{y} _o \cup \mathbf{z} ) = 1$. A full description of the vector and indexing notation can be found in Katzfuss and Guinness (2019, App. A).

Computing $\mathbf{U }$

We recapitulate here the formulas for computing $\mathbf{U }$ from Katzfuss and Guinness (2019). Let g(i) denote the vector of indices of the elements in $\mathbf{x}$ on which $x_i$ conditions. Also define $C(x_i,x_j)$ as the covariance between $x_i$ and $x_j$ implied by the true model in Sect. 2; that is, $C(y_i,y_j) = C(z_i,y_j) = K(\mathbf{s} _i,\mathbf{s} _j)$ and $C(z_i,z_j) = K(\mathbf{s} _i,\mathbf{s} _j) + \mathbb {1}_{i=j} \tau ^2_i$. Then, the (j, i)th element of $\mathbf{U }$ can be calculated as

$$\begin{aligned} \mathbf{U }_{ji} = {\left\{ \begin{array}{ll} d_i^{-1/2}, &{} i=j,\\ -b_{i}^{(j)} d_i^{-1/2}, &{} j \in g(i), \\ 0, &{}\text {otherwise}, \end{array}\right. } \end{aligned}$$

(9)

where $\mathbf{b} _i'= C(x_i,\mathbf{x} _{g(i)}) C(\mathbf{x} _{g(i)},\mathbf{x} _{g(i)})^{-1}$, $d_i = C(x_i,x_i) - \mathbf{b} _i' C(\mathbf{x} _{g(i)},x_i)$, and $b_i^{(j)}$ denotes the kth element of $\mathbf{b} _i$ if j is the kth element in g(i) (i.e., $b_i^{(j)}$ is the element of $\mathbf{b} _i$ corresponding to $x_j$).

Prediction with Unknown Parameters

In practice, most GP models depend on an unknown parameter vector $\varvec{\theta }$, which we will make explicit here. Likelihood approximation for parameter inference is discussed in detail in Katzfuss and Guinness (2019), but we will review it briefly here. Integration of $\widehat{f}(\mathbf{x} |\varvec{\theta })$ in (3) with respect to $\mathbf{y}$ results in the following Vecchia likelihood (Katzfuss and Guinness 2019, Prop. 2):

$$\begin{aligned} \textstyle -2 \log \widehat{f}(\mathbf{z} _o|\varvec{\theta })= & {} -2 \sum _{i} \log \mathbf{U }_{ii} + 2\sum _{i} \log \mathbf{V} _{ii} + {\tilde{\mathbf{z }}}'{\tilde{\mathbf{z }}} \nonumber \\&- (\mathbf{V} ^{-1}\mathbf{U }_{\ell ,\bullet }{\tilde{\mathbf{z }}})'(\mathbf{V} ^{-1}\mathbf{U }_{\ell ,\bullet }{\tilde{\mathbf{z }}}) + n \log (2\pi ), \end{aligned}$$

(10)

where $\mathbf{U }$ and $\mathbf{V}$ implicitly depend on $\varvec{\theta }$, and ${\tilde{\mathbf{z }}} = \mathbf{U }_{r,\bullet }'\mathbf{z} _o$. The computational cost for evaluating this Vecchia likelihood is often low, and Katzfuss and Guinness (2019) provide conditions on the g(i) under which the cost is guaranteed to be linear in n.

This allows for various forms of likelihood-based parameter inference. In a frequentist setting, we can compute ${\hat{\varvec{\theta }}} = \hbox {arg max}_{\varvec{\theta }} \log \widehat{f}(\mathbf{z} _o|\varvec{\theta })$ and then compute summaries of the posterior predictive distribution $\widehat{f}(\mathbf{y} |\mathbf{z} _o) = \mathcal {N}_n(\varvec{\mu }({\hat{\varvec{\theta }}}),\varvec{\Sigma }({\hat{\varvec{\theta }}}))$ as described in Sect. 3.3. This often has low computational cost but ignores uncertainty in ${\hat{\varvec{\theta }}}$. An example is given in Sect. 6.

In a Bayesian setting, given a prior distribution $f(\varvec{\theta })$, a Metropolis–Hastings sampler can be used for parameters whose posterior or full-conditional distribution are not available in closed form. At the $(l+1)$th step of the algorithm, one would propose a new value $\varvec{\theta }^{(P)} \sim q(\varvec{\theta }|\varvec{\theta }^{(l)})$ and accept it with probability $\min (1,h(\varvec{\theta }^{(P)},\varvec{\theta }^{(l)})/h(\varvec{\theta }^{(l)},\varvec{\theta }^{(P)}))$, where $h(\varvec{\theta },{\tilde{\varvec{\theta }}}) = f(\varvec{\theta })\widehat{f}(\mathbf{z} _o|\varvec{\theta }) q({\tilde{\varvec{\theta }}}|\varvec{\theta })$. After burn-in and thinning, this results in a sample, say, $\varvec{\theta }^{(1)},\ldots ,\varvec{\theta }^{(L)}$, leading to a Gaussian-mixture prediction: $\widehat{f}(\mathbf{y} |\mathbf{z} _o) = (1/L) \sum _{l=1}^L \mathcal {N}_n(\varvec{\mu }(\varvec{\theta }^{(l)}),\varvec{\Sigma }(\varvec{\theta }^{(l)}))$. See Katzfuss and Guinness (2019, App. E) for an example for the use of RF-full in this setting. For more complicated Bayesian hierarchical models, inference can be carried out using a Gibbs sampler in which $\mathbf{y}$ is sampled from its full-conditional distribution as described in item 4. in Sect. 3.3.

Numerical Nonzeros in $\mathbf{V}$

For simplicity, we focus here on RF-full, although numerical nonzeros can similarly occur for RF-stand, LF-full, and LF-ind (see Figure S1). For RF-full, the upper triangle of $\mathbf{W}$ is at least as dense as the upper triangle of $\mathbf{V}$. Specifically, for $j<i$, $\mathbf{V} _{ji}=\mathbf{U }_{\ell _j\ell _i}=0$ unless $j \in q_y(i)$. From Katzfuss and Guinness (Katzfuss and Guinness (2019), Prop. 3.2), we have that $\mathbf{W} _{ji} =0$ unless $j \in q_y(i)$ or $\exists k > i$ such that $i,j \in q_y(k)$. Thus, for any pair $j<i$ such that $j \notin q_y(i)$ but $i,j \in q_y(k)$ for some $k>i$, we generally have $\mathbf{V} _{ji} = 0$ and $\mathbf{W} _{ji} \ne 0$.

From the standard Cholesky algorithm, we can derive that the algorithm for $\mathbf{V} = {{\,\mathrm{rchol}\,}}(\mathbf{W} )$ computes $\mathbf{V} _{ji}$ as

$$\begin{aligned} \textstyle \mathbf{V} _{ji} = (\mathbf{W} _{ji} - \sum _{k=i+1}^n\mathbf{V} _{ik}\mathbf{V} _{jk})/\mathbf{V} _{ii}. \end{aligned}$$

Thus, for any pair $j<i$ such that $\mathbf{V} _{ji}=0$ but $\mathbf{W} _{ji} \ne 0$, we know that $\mathbf{W} _{ji} = \sum _{k=i+1}^n\mathbf{V} _{ik}\mathbf{V} _{jk}$ theoretically, but due to potential numerical error it is not guaranteed that this equation holds exactly. A numerical nonzero is introduced in $\mathbf{V}$ whenever a rounding error occurs in $\sum _{k=i+1}^n\mathbf{V} _{ik}\mathbf{V} _{jk}$, which relies on (potentially many) previous calculations in the Cholesky algorithm. Such numerical nonzeros are avoided by extracting $\mathbf{V} =\mathbf{U }_{\ell \ell }$ as a submatrix of $\mathbf{U }$ (as proposed in Sect. 4.1), instead of explicitly carrying out the Cholesky factorization $\mathbf{V} = {{\,\mathrm{rchol}\,}}(\mathbf{W} )$.

1.1 Implications for Selected Inverse

When $\mathbf{V}$ is computed by copying a submatrix from $\mathbf{U }$ to avoid numerical nonzeros, the selected inverse of this $\mathbf{V}$ is not guaranteed to return the exact posterior variances of $\mathbf{y}$, unless $\mathbf{V}$ is “padded” with zeros, which results in additional costs. This is because the selected inverse operates on the symbolic nonzero elements; that is, it operates on all elements that have to be computed in the Cholesky, even if they cancel to zero numerically (which is the case for many entries in our case). Denoting by $\mathbf{S}$ the selected inverse of $\mathbf{W}$ based on $\mathbf{V}$, a close look at the Takahashi recursions reveals that for all j, k with $j,k \in q_y(i)$, we need $\mathbf{S} _{ji}$ and $\mathbf{S} _{kj}$. The latter element is only calculated if $j \in q_y(k)$. However, if $j \notin q_y(k)$, $\mathbf{S} _{kj} ={{\,\mathrm{cov}\,}}(y_j,y_k|\mathbf{z} )$ will typically be very small (if m is reasonably large), because their corresponding locations will likely be far away from each other and data can be observed in between. In our experiments, the additional approximation error introduced by the selected inverse was negligible relative to the error introduced by the Vecchia approximation itself. When m becomes large enough for the Vecchia approximation to be accurate, the additional approximation error introduced by SelInv goes to zero as well.

If the exact variances implied by the Vecchia approximation $\widehat{f}(\mathbf{y} |\mathbf{z} _o)$ are desired, they can be computed as ${{\,\mathrm{diag}\,}}({{\,\mathrm{var}\,}}(\mathbf{y} _p|\mathbf{z} _o)) = ((\mathbf{V} ^{-1}\mathbf{I} _{\bullet p})\circ (\mathbf{V} ^{-1}\mathbf{I} _{\bullet p}))'\mathbf 1 _{n}$. For RF-full and RF-stand, this requires $\mathcal {O}(nmn_P)$ time, as described in Sect. 4.3.3, and so the overall computational complexity would not be increased if $n_P = \mathcal {O}(m^2)$.

Proofs

In this section, we provide proofs for the propositions stated throughout the article.

Proof of Proposition 1

Part 1 of this proposition is equivalent to Thm. 1 in Guinness (2018). We prove the statement here again in a different way, which can be easily extended to prove Part 2.

First, consider a generic Vecchia approximation of a vector $\mathbf{x}$ of length n: $\widehat{f}(\mathbf{x} ) = \prod _{i=1}^n f(x_i|\mathbf{x} _{g(i)}) = \prod _{i=1}^n \mathcal {N}(x_i|\mu _{i|g(i)},\sigma ^2_{i|g(i)})$. Then,

$$\begin{aligned} \textstyle (-2) {{\,\mathrm{\mathbb {E}}\,}}^\mathbf{x} \log \widehat{f}(\mathbf{x} ) = \sum _{i=1}^n \log {{\,\mathrm{var}\,}}(x_i|\mathbf{x} _{g(i)}) + \sum _{i=1}^n {{\,\mathrm{\mathbb {E}}\,}}^\mathbf{x} w_i^2 + n \log (2\pi ), \end{aligned}$$

where $w_i = (x_i - \mu _{i|g(i)})/\sigma _{i|g(i)}$. We have ${{\,\mathrm{\mathbb {E}}\,}}(w_i) = 0$, because ${{\,\mathrm{\mathbb {E}}\,}}(\mu _{i|g(i)}) = {{\,\mathrm{\mathbb {E}}\,}}{{\,\mathrm{\mathbb {E}}\,}}(x_i|\mathbf{x} _{g(i)}) = {{\,\mathrm{\mathbb {E}}\,}}(x_i)$. We also have ${{\,\mathrm{var}\,}}(w_i)=1$, because

$$\begin{aligned} {{\,\mathrm{var}\,}}(x_i-\mu _{i|g(i)})= & {} {{\,\mathrm{var}\,}}{{\,\mathrm{\mathbb {E}}\,}}\big (x_i-{{\,\mathrm{\mathbb {E}}\,}}(x_i|\mathbf{x} _{g(i)})|\mathbf{x} _{g(i)}\big ) + {{\,\mathrm{var}\,}}\big (x_i-{{\,\mathrm{\mathbb {E}}\,}}(x_i|\mathbf{x} _{g(i)})|\mathbf{x} _{g(i)}\big ) \\= & {} 0 + \sigma ^2_{i|g(i)}. \end{aligned}$$

Hence, $w_i \sim \mathcal {N}(0,1)$, and so ${{\,\mathrm{\mathbb {E}}\,}}^\mathbf{x} w_i^2 = 1$ and

$$\begin{aligned} \textstyle (-2) {{\,\mathrm{\mathbb {E}}\,}}^\mathbf{x} \log \widehat{f}(\mathbf{x} ) = \sum _{i=1}^n \log {{\,\mathrm{var}\,}}(x_i|\mathbf{x} _{g(i)}) + n + n \log (2\pi ). \end{aligned}$$

Because the exact density $f(\mathbf{x} )$ is a special case of Vecchia with $g(i) = (1,\ldots ,i-1)$, we have

$$\begin{aligned} \textstyle \text {KL}\big (f(\mathbf{x} )\Vert \widehat{f}(\mathbf{x} )\big ) = {{\,\mathrm{\mathbb {E}}\,}}^\mathbf{x} \log f(\mathbf{x} ) - {{\,\mathrm{\mathbb {E}}\,}}^\mathbf{x} \log \widehat{f}(\mathbf{x} ) = \frac{1}{2} \sum _{i=1}^n \log \frac{{{\,\mathrm{var}\,}}(x_i|\mathbf{x} _{g(i)})}{{{\,\mathrm{var}\,}}(x_i|\mathbf{x} _{1:i-1})}. \end{aligned}$$

(11)

For $g_1(i) \subset g_2(i)$, we can write, say, $g_2(i) = g_1(i) \cup c(i)$. Using the law of total variance,

$$\begin{aligned} {{\,\mathrm{var}\,}}(x_i|\mathbf{x} _{g_1(i)}) = {{\,\mathrm{var}\,}}(x_i|\mathbf{x} _{g_1(i)}, \mathbf{x} _{c(i)}) + {{\,\mathrm{var}\,}}({{\,\mathrm{\mathbb {E}}\,}}(x_i|\mathbf{x} _{g_1(i)}, \mathbf{x} _{c(i)})|\mathbf{x} _{c(i)}) \ge {{\,\mathrm{var}\,}}(x_i|\mathbf{x} _{g_2(i)}).\nonumber \\ \end{aligned}$$

(12)

Now, Part 1 follows by combining (11) and (12) with the assumption that $g_1(i) \subset g_2(i)$ for all i.

For Part 2, we consider $\mathbf{x} =(\mathbf{x} ^{(1)},\mathbf{x} ^{(2)})$ with $\mathbf{x} ^{(1)} = (x_1,\ldots ,x_{n_1})$ and $\mathbf{x} ^{(2)} = (x_{n_1+1},\ldots ,x_{n_1+n_2})$. Then,

$$\begin{aligned} \textstyle&\text {CKL}\big (f(\mathbf{x} ^{(2)}|\mathbf{x} ^{(1)})\Vert \widehat{f}(\mathbf{x} ^{(2)}|\mathbf{x} ^{(1)})\big )\\&\quad = \int f(\mathbf{x} ^{(1)}) \int f(\mathbf{x} ^{(2)}|\mathbf{x} ^{(1)}) \log \big (f(\mathbf{x} ^{(2)}|\mathbf{x} ^{(1)})/\widehat{f}(\mathbf{x} ^{(2)}|\mathbf{x} ^{(1)})\big ) d\mathbf{x} ^{(2)} d\mathbf{x} ^{(1)} \\&\quad = {{\,\mathrm{\mathbb {E}}\,}}^\mathbf{x} \log f(\mathbf{x} ^{(2)}|\mathbf{x} ^{(1)}) - {{\,\mathrm{\mathbb {E}}\,}}^\mathbf{x} \log \widehat{f}(\mathbf{x} ^{(2)}|\mathbf{x} ^{(1)}) \\&\quad = \frac{1}{2} \sum _{i=n_1+1}^{n_1+n_2} \log \frac{{{\,\mathrm{var}\,}}(x_i|\mathbf{x} _{g(i)})}{{{\,\mathrm{var}\,}}(x_i|\mathbf{x} _{1:i-1})}, \end{aligned}$$

where the last equality can be shown almost identically to (11) above, noting that $\widehat{f}(\mathbf{x} ^{(2)}|\mathbf{x} ^{(1)}) = \widehat{f}(\mathbf{x} )/\widehat{f}(\mathbf{x} ^{(1)}) = \prod _{i=n_1+1}^{n_1+n_2} f(x_i|\mathbf{x} _{g(i)})$. Part 2 follows by combining this result and (12) with the assumption that $g_1(i) \subset g_2(i)$ for all $i=n_1+1,\ldots ,n_1+n_2$. $\square$

Proof of Proposition 2

For all parts of this proposition, note that all variables in the model are conditionally independent of $z_j$ given $y_j$, and so conditioning on $y_j$ is equivalent to conditioning on $y_j$ and $z_j$. For Part 1, we can thus verify easily that g(i) is equivalent to $(1,\ldots ,i-1)$, for all i and all methods under consideration, and so Part 1 follows from (11). Using Proposition 1, the proof for all other parts simply consists of showing that certain conditioning vectors contain certain other conditioning vectors (similar to Katzfuss and Guinness 2019, Prop. 4). For example, for Part 3, all response-first methods are based on the ordering $\mathbf{x} = (\mathbf{z} _o,\mathbf{y} _o,\mathbf{y} _p)$ with nearest neighbor conditioning (under some restrictions for RF-stand and RF-ind), and so it is easy to see that $g_m(i) \subset g_{m+1}(i)$ for all $i \in \ell _p=(n+1,\ldots ,n+n_P)$. LF-full and LF-ind can equivalently be defined based on the ordering $(\mathbf{y} _o,\mathbf{z} _o,\mathbf{y} _p)$, in which case we also have $g_m(i) \subset g_{m+1}(i)$ for all $i \in \ell _p=(n+1,\ldots ,n+n_P)$. For Part 5, note that in RF-stand and RF-ind, $\mathbf{y} _p$ does not condition on $\mathbf{y} _o$. Hence, the distribution $\widehat{f}(\mathbf{y} _p|\mathbf{z} _o)$ is equivalent to the distribution obtained under a simplified Vecchia approximation based on $\mathbf{x} =(\mathbf{z} _o,\mathbf{y} _p)$ (i.e., with $\mathbf{y} _o$ removed completely). It is straightforward to show that Part 5 holds for this simplified approximation. For Part 6, note that g(i) is the same for RF-full and RF-stand for all $i=1,\ldots ,n$. For $i \in p$, letting $a(i) = q_y^{\text {RF-full}}(i) \cap o$, we have $\mathbf{x} _{g^{\text {RF-full}}(i)}=(\mathbf{y} _{q_y^{\text {RF-stand}}},\mathbf{y} _{a(i)})=(\mathbf{y} _{q_y^{\text {RF-stand}}},\mathbf{y} _{a(i)},\mathbf{z} _{a(i)})$ and $\mathbf{x} _{g^{\text {RF-stand}}(i)}=(\mathbf{y} _{q_y^{\text {RF-stand}}},\mathbf{z} _{a(i)})$, and so $g^{\text {RF-stand}}(i) \subset g^{\text {RF-full}}(i)$ for all $i \in p$. For Part 7, note that a GP with exponential covariance in 1-D is a Markov process, and so in (6) with $q_y(i)=i-1$ and left-to-right ordering, we have $f(y_i|\mathbf{y} _{q_y(i)})=f(y_i|y_{1},\ldots ,y_{i-1})$ and hence $\widehat{f}(\mathbf{x} ) = f(\mathbf{x} )$. $\square$