Radial Neighbors for Provably Accurate Scalable Approximations of Gaussian Processes

Let the spatial domain $\Omega$ be a connected subset of $\mathbb{R}^{d}$. Let $Z=(Z_{s}:s\in\Omega)$ be a real valued Gaussian process on $\Omega\subset\mathbb{R}^{d}$ with zero mean and covariance function $K:\Omega\times\Omega\to\mathbb{R}$ such that for all $s_{1},s_{2}\in\Omega$, $\mathbb{E}(Z_{s_{1}})=0$ and $\mathrm{Cov}(Z_{s_{1}},Z_{s_{2}})=K(s_{1},s_{2})$. The objective of this paper is to approximate the process $Z$ with another process $\hat{Z}$ such that the distance between $\hat{Z}$ and $Z$ is small in some metric and finite marginal densities of $\hat{Z}$ are easy to compute.

Let $\mathcal{D}=\{w_{1},w_{2},\ldots,w_{n}\}$ be a finite collection of locations in $\Omega$. A $\rho$-radial neighbors graph on $\mathcal{D}$ is defined to be a directed acyclic graph such that for any $w_{i},w_{j}\in\mathcal{D}$ with $\|w_{i}-w_{j}\|_{2}<\rho$, there is a directed edge between $w_{i}$ and $w_{j}$. Such radial neighbors graphs are not unique. Section \arabicsection.\arabicsubsection provides one practical way of constructing a radial neighbors graph that will be used in our experimental studies, while Section \arabicsection.\arabicsubsection presents how to build a Gaussian process from any radial neighbors graph. Our theory in Section \arabicsection applies to any Gaussian process built from a radial neighbors graph, not just the specific graph we constructed in Section \arabicsection.\arabicsubsection.

We present a procedure to construct the radial neighbors graph that will be used for experiments in Section \arabicsection. First, we specify an order of all locations in $\mathcal{D}$. This is done by choosing a “center” location of $\mathcal{D}$ and sorting all the locations in $\mathcal{D}$ according to their distances to this “center” in a non-decreasing order. If the locations in $\mathcal{D}$ are roughly evenly distributed in the domain $\Omega$, then the coordinates of this center are set as the arithmetic mean of the coordinates of all locations in $\mathcal{D}$.

Next, once an order of locations in $\mathcal{D}$ is obtained, for each location $w_{i}\in\mathcal{D}$, we set the parent set of $w_{i}$ to be all the locations preceeding $w_{i}$ in this order and within $\rho$ distance to $w_{i}$. As a consequence, for all locations within $\rho$ radius of $w_{i}$, the ones that precede it are its parents, whereas the ones that follow it are its children. This construction satisfies the definition of radial neighbors graphs in Section \arabicsection.\arabicsubsection that any pair of locations $(w_{i},w_{j})$ with $\|w_{i}-w_{j}\|_{2}<\rho$ have a directed edge connecting them.

Depending on the spatial dispersion of locations in $\mathcal{D}$, some locations may not have neighbors within $\rho$ distance. We thus populate the parent set of a location $w_{i}\in\mathcal{D}$ without $\rho$-neighbors by drawing an edge from $w_{i}$’s nearest neighbor among all locations that precede $w_{i}$ in the order. This step ensures that the resulting directed acyclic graph is connected.

If the locations in $\mathcal{D}$ are roughly evenly distributed in the domain $\Omega$, once the first $i$ locations are chosen and the directed graph on them is constructed, our procedure roughly chooses the $(i+1)$th location to be as close to the first $i$ locations as possible. Since spatial covariance function $K(s_{1},s_{2})$ is often relatively large when $s_{1}$ and $s_{2}$ are close to each other, our procedure aims to minimize the conditional variance of the $(i+1)$th location given the first $i$ locations, or in other words, maximize the spatial information that the $(i+1)$th location can borrow from the first $i$ locations. Figure \arabicfigure illustrates the process of building this radial neighbors graph when $\mathcal{D}$ is a $15\times 15$ grid on $[0,1]^{2}$.

Finally, if we have a new data set $\tilde{\mathcal{D}}$, the radial neighbors graph constructed on $\mathcal{D}$ can be extended to a radial neighbors graph on $\mathcal{D}\cup\tilde{\mathcal{D}}$. Specifically, once an ordering of locations in $\tilde{\mathcal{D}}$ is obtained, we concatenate the ordering of locations in $\mathcal{D}$ and $\tilde{\mathcal{D}}$ to obtain an ordering of locations $\mathcal{D}\cup\tilde{\mathcal{D}}$. We then construct a radial neighbors graph on $\mathcal{D}\cup\tilde{\mathcal{D}}$ via the same procedure as described above. The final graph will contain the graph on $\mathcal{D}$ as a subgraph. This extension framework will be useful when extending a graph from training data to test data.

For some dataset, the spatial dispersion of both training and testing locations can be highly uneven. While a radial neighbor graph ensures all locations are connected in the same directed graph, it can still result in poor predictions on some test locations. We can mitigate this issue by artificially adding more test locations, so that each test locations receives a decent number of parents.

Given a $\rho$-radial neighbors graph on $\mathcal{D}=\{w_{1},w_{2},\ldots\,w_{n}\}$ defined in Section \arabicsection.\arabicsubsection, we can build a RadGP on $\mathcal{D}$ in which each $w_{i}$ is conditionally only dependent on its parent set. Let the symbol $\overset{d}{=}$ denote equality in distribution and define $\Sigma_{A,B}$ as the covariance matrix between two spatial location sets $A$ and $B$ under the covariance function $K(\cdot,\cdot)$. The conditional distribution of $\hat{Z}_{w_{i}}$ is defined as

where for $i=1$ the parent set $\mathrm{pa}(w_{1})$ is empty.

Equation (\arabicsection.\arabicsubsection) defines the distribution of $\hat{Z}$ on all finite subsets of $\mathcal{D}$, which essentially says the conditional distribution of $\hat{Z}_{w_{i}}$ given the process $\hat{Z}$ on all locations preceding $w_{i}$ has the same law as the conditional distribution of $Z_{w_{i}}$ given the process $Z$ on the parent set $\mathrm{pa}(w_{i})$. We then extend this process $\hat{Z}$ to arbitrary finite subsets of the whole space $\Omega$. For all $s\in\Omega\backslash\mathcal{D}$, we define the parents of $s$ as $\mathrm{pa}(s)=\{s^{\prime}\in\mathcal{D}:\|s^{\prime}-s\|_{2}<\rho\}$. For any finite set $U\subset\Omega\backslash\mathcal{D}$, we define the conditional distribution of $\hat{Z}_{U}$ given $\hat{Z}_{\mathcal{T}_{1}}$ as

Equations (\arabicsection.\arabicsubsection) and (\arabicsection.\arabicsubsection) together give the joint distribution of $\hat{Z}$ on any finite subset of $\Omega$. Specifically, for a generic finite set $A=U\cup V$ with $U\subset\Omega\backslash\mathcal{D}$ and $V\subset\mathcal{D}$, the joint density of $\hat{Z}$ on $A$ implied by Equations (\arabicsection.\arabicsubsection) and (\arabicsection.\arabicsubsection) is

Equations (\arabicsection.\arabicsubsection)-(\arabicequation) complete the definition of our RadGP.

The proof of Lemma \arabicsection.\arabictheorem relies on the Kolmogorov extension theorem and is available at Section S1 of the Supplementary Material.

The objective of this section is to establish theoretical upper bounds for the Wasserstein distance between the original Gaussian process and any RadGP on $\mathcal{D}$ described above. For two generic random vectors $Z_{1},Z_{2}$ defined on the same space $\mathcal{Z}$ following the probability distributions $\mu_{1},\mu_{2}$, their Wasserstein-2 ($W_{2}$) distance is

where $\mathcal{M}$ is the set of all probability distributions on the product space $\mathcal{Z}\times\mathcal{Z}$ whose marginals are $\mu_{1}$ and $\mu_{2}$. Let $Z_{\mathcal{D}}$ and $\hat{Z}_{\mathcal{D}}$ be the random vectors from the original Gaussian process and the RadGP on the set $\mathcal{D}$. Denote the covariance matrix of $Z_{\mathcal{D}}$ and $\hat{Z}_{\mathcal{D}}$ as $\Sigma_{\mathcal{D}\mathcal{D}}$ and $\hat{\Sigma}_{\mathcal{D}\mathcal{D}}$, respectively. Then the Wasserstein-2 distance between $Z_{\mathcal{D}}$ and $\hat{Z}_{\mathcal{D}}$ has a closed form (Gelbrich, 1990):

where $\mathrm{tr}(A)$ denotes the trace of a square matrix $A$. The right-hand side of equation (\arabicequation) involves various matrix powers that are difficult to analyze. Fortunately, for Gaussian measures, the squared Wasserstein-2 distance can be upper bounded by the trace norm of the difference between their covariance matrices, which greatly simplifies our analysis.

We consider the Gaussian process $Z=(Z_{s}:s\in\Omega)$ with zero mean and an isotropic nonnegative covariance function $K(\cdot,\cdot)$. Thus we can reformulate $K$ as $K(s_{1},s_{2})=K_{0}(\|s_{1}-s_{2}\|_{2})$ for all $s_{1},s_{2}\in\Omega$, where $K_{0}:[0,+\infty)\to[0,+\infty)$. For all $r>0$, define the function $v_{r}(\cdot):\mathbb{R}\to\mathbb{R}_{+}$ as $v_{r}(x)=\sum_{k=0}^{+\infty}|x|^{k}/(k!)^{r}$. Then $1/v_{r}(x)$ is a monotone decreasing function of $x$, where a larger $r$ results in a slower decay rate in $1/v_{r}(x)$. Specifically, $1/v_{1}(x)=\exp(-x)$ and $1/v_{r}(x)$ with $r>1$ has a decay rate faster than all polynomials but slower than the exponential function. We also define a series of polynomials $c_{r}(x)=1+x^{r}$ for $r>0$. For each $v_{r}$ with $r>1$ and each $c_{r}$ with $r>0$, we define the following families of Gaussian processes:

Intuitively, $\mathscr{Z}_{v_{r}}$ is the family of isotropic Gaussian processes whose covariance function decays no slower than some subexponential rate of the spatial distance, while $\mathscr{Z}_{c_{r}}$ is the family with covariance functions decaying no slower than an $r$th order polynomial of the spatial distance.

The objective of our theoretical study is to bound the Wasserstein-2 distance between the marginal distribution of RadGP and original Gaussian process on $\mathcal{D}$ when the original process belongs to the two spatial decaying families in (\arabicequation) and (\arabicequation):

For the finite set $\mathcal{D}$ described above, we define the minimal separation distance among all spatial locations in $\mathcal{D}$ as $q=\min_{1\leq i<j\leq|\mathcal{D}|}\|w_{i}-w_{j}\|_{2}$. This minimal separation distance is a parsimonious statistic for describing the spatial dispersion of $\mathcal{D}$, and is useful in bounding various quantities such as the maximal eigenvalue and condition number of $\Sigma_{\mathcal{D}\mathcal{D}}$; see Lemma S2 of the Supplementary Material for details. Let $\hat{K}_{0}(w)=(2\pi)^{-d/2}\int K_{0}(x)\exp(-\imath w^{{\mathrm{\scriptscriptstyle T}}}x)dx$ be the Fourier transform of $K_{0}$ where $\imath^{2}=-1$, and let $\phi_{0}(x)=\inf_{\|w\|_{2}\leq 2x}\hat{K}_{0}(w)$ for all $x>0$. For two positive sequences $a_{n}$ and $b_{n}$, we use $a_{n}\lesssim b_{n}$ to denote the relation that $a_{n}\leq Cb_{n}$ for all $n$ and some finite constant $C$. We have the following theorems on the Wasserstein-2 distance between the RadGP and the original Gaussian process. The constants $c_{1},c_{2}$ are from Lemma S2 while $c_{3}$ is from the proof of Lemma S5 in the Supplementary Material, all of which only depend on the dimension $d$ of the spatial domain $\Omega$.

The conditions for the above theorems, namely the equations regarding $c_{6},c^{\prime}_{6},c_{7},c^{\prime}_{7}$, are necessary but not sufficient conditions for the upper bounds of Wasserstein-2 distance to go to zero. Thus, they can be safely ignored if we aim for accurate approximations. We now give a detailed analysis of these upper bounds for Wasserstein-2 distance involving three variables: the sample size $n$, the minimal separation distance $q$, and the approximation radius $\rho$ in the RadGP. We will discuss from the perspectives of both finite samples and asymptotics. From a finite sample perspective, we fix two out of three variables and discuss how and why the upper bounds change as the remaining variable changes; from an asymptotic perspective, we describe the relations between $n,q$ and $\rho$ such that as $n$ goes to infinity, the Wasserstein-2 distance goes to zero.

Viewing Theorems \arabicsection.\arabictheorem and \arabicsection.\arabictheorem as finite sample bounds, the upper bounds are linearly dependent on the sample size $n$. This directly follows from the definition of Wasserstein distance which measures the difference between two $n$-dimensional Gaussian distributions. Provided both the minimal separation distance $q$ and approximation radius $\rho$ are fixed, as the sample size $n$ increases, the local structure for existing locations will not change and the approximation accuracy for them will not improve. Thus, adding more locations will lead to a linearly increasing approximation error in the squared $W_{2}$ distance.

Second, the upper bounds on the squared $W_{2}$ distance decay with respect to the approximation radius $\rho d^{-1/2}$ at the same rate as the spatial covariance function with respect to the spatial distance, up to a $(d+1)$th order polynomial. For Theorem \arabicsection.\arabictheorem, the covariance function is required to decay faster than $1/v_{r}(x)$ by a polynomial $1+x^{d+1}$ as defined in equation (\arabicequation); for Theorem \arabicsection.\arabictheorem, the covariance function decays faster than an $r$th order polynomial while the upper bound of $W_{2}$ distance decays at an $(r-d-1)$th order polynomial. The difference between the decay rate of the covariance function and the $W_{2}$ distance in Theorems \arabicsection.\arabictheorem and \arabicsection.\arabictheorem comes from the fact that high dimension reduces the effective decay rate of covariance function. For example, under a fixed minimal separation distance, when $d=1$, as long as the covariance function decays faster than $1/x$ by a polynomial term, the matrix $\Sigma_{\mathcal{D}\mathcal{D}}$ has a bounded $l_{1}$ norm regardless of sample size; however, for a general $d$ dimensional space, the covariance function must decay faster than $1/x^{d}$ for $\Sigma_{\mathcal{D}\mathcal{D}}$ to have a bounded $l_{1}$ norm. These upper bounds show that the approximation radius $\rho$ indeed controls the approximation accuracy. Specifically, for fixed $q$ and $n$, we can reach an arbitrarily small $W_{2}$ error by increasing $\rho$ over some threshold value; see our detailed conditions on $\rho$ in Corollary \arabicsection.\arabictheorem for various types of covariance functions.

Third, for fixed $\rho$ and $n$, both upper bounds in the two theorems increase as the minimal separation distance $q$ decreases. This is because as the minimal separation distance decreases, spatial locations get closer to each other and the covariance matrix gets close to singular, leading to a larger approximation error.

From the asymptotic perspective, we can increase the sample size $n$ while viewing the minimal separation distance $q$ and the approximation radius $\rho$ as some functions of $n$. The objective is to choose an approximation radius $\rho$ dependent on $n$ and $q$ such that the squared Wasserstein distance goes to zero as $n$ goes to infinity. We first present a corollary showing how to choose such $\rho$ for three commmonly used covariance functions, where the convergence is understood in the sense that $q$ is also a function of $n$.

We then analyze our main theorems from the asymptotic perspective and verify the analysis using the three examples in Corollary \arabicsection.\arabictheorem. Depending on whether and how the minimal separation distance $q$ changes with respect to the sample size $n$, there are three commonly used asymptotic frameworks: fixed domain, mixed domain and increasing domain asymptotics.

Under the increasing domain setting, as $n$ goes to infinity, the radius of the domain $\Omega$ increases while the minimal separation distance $q$ is fixed. Therefore, as long as $n\gtrsim d^{1/2}v_{r}^{-1}(n)$ under conditions of Theorem \arabicsection.\arabictheorem or $n\gtrsim d^{1/2}n^{1/(r-d-1)}$ under conditions of Theorem \arabicsection.\arabictheorem, we have $o(1)$ approximation error. When applied to Matérn covariance in case (1) and Gaussian covariance in case (2) of Corollary \arabicsection.\arabictheorem, this yields a slowly increasing rate of $O(\ln^{3}n)$ for the approximation radius $\rho$, which is much smaller than the sample size $n$. On the other hand, when the covariance function decays polynomially such as the generalized Cauchy in case (3) of Corollary \arabicsection.\arabictheorem, the approximation radius $\rho$ needs to increase at a polynomial rate of sample size $n$. This is because for covariance functions that decrease slowly in the spatial distance, one needs more spatial locations and a larger approximation radius to capture the dependence information.

Under the mixed domain setting, as $n$ goes to infinity, the radius of the domain $\Omega$ increases with $n$ while the minimal separation distance $q$ also decreases with $n$. By solving the inequalities on the right hand side of equations (\arabicequation) and (\arabicequation) being no greater than $o(1)$, we obtain the lower bound of $\rho$ that guarantees the $o(1)$ approximation error. Corollary \arabicsection.\arabictheorem directly shows such lower bounds for Matérn, Gaussian and generalize Cauchy covariance functions. We can still achieve the $o(1)$ approximation error as long as these lower bounds of $\rho$ in the three cases of Corollary \arabicsection.\arabictheorem increase slower than the radius of the domain $\Omega$. Otherwise, if $\rho$ increases exactly at the same rate as the radius of $\Omega$, RadGP is identical to the original Gaussian process and does not make any approximation at all. Our theory is not applicable when the lower bounds of $\rho$ exceed the radius of domain $\Omega$. In particular, for the fixed domain setting, the radius of $\Omega$ is fixed as a constant and the minimal separation distance $q$ decreases to zero at a rate no slower than $O(n^{-1/d})$. In such case, the lower bounds of $\rho$ in Corollary \arabicsection.\arabictheorem will diverge with $n$, and therefore Theorems \arabicsection.\arabictheorem and \arabicsection.\arabictheorem do not provide a vanishing bound for Wasserstein-2 distances.

Our results fill a major gap in the literature by providing the previously lacking theoretical support for popular local or neighborhood-based Gaussian process approximations. Our theory only requires knowledge of the spatial decaying pattern of the covariance functions (e.g., exponential and polynomial in $\mathscr{Z}_{v_{r}}$ and $\mathscr{Z}_{c_{r}}$, respectively). These conditions are easily verifiable, with derivations for the Matérn, Gaussian and the generalized Cauchy covariance functions provided in Corollary \arabicsection.\arabictheorem. Similar techniques can be used for a wider range of stationary covariance functions. Furthermore, our novel theory explicitly links the approximation quality from choosing radius $\rho$ with the minimal separation distance $q$ and the sample size $n$, both of which are readily available from observed data.

Throughout Section 5, we consider the specific RadGP described in Section 2. We compare the RadGP prior with the nearest-neighbor Gaussian process (NNGP, Datta et al. 2016) implemented on a maximin ordering of the data (Guinness, 2018), in terms of their approximation accuracy. The true Gaussian process has mean zero and an isotropic Matérn covariance function

with the true parameters $(\phi,\tau^{2},\nu)=(31.63,1,1.5)$. The value $\phi=31.63$ is obtained by setting $K_{0}(0.15)=0.05$, a choice suggested by Katzfuss et al. (2020). The experiment is carried out on a $40\times 40$ grid in $[0,1]^{2}$. According to the analysis of computational complexity in Section \arabicsection, we use the average number of nonzero elements per column in the precision matrix as a measure of complexity for both methods. The quality of approximation is measured by squared Wasserstein-2 distance between each approximate process and the true Gaussian process, which can be computed efficiently using the R package transport available from CRAN (Schuhmacher et al. 2024). As shown in Figure \arabicfigure, the results are generally in favor of RadGP, as it achieves smaller approximation error under various complexity of the directed graphs.

We study the performance of RadGP for the spatial regression model (\arabicequation) with no covariates, where the true covariance function for spatial effects $Z(s)$ is the same as that in Section \arabicsection.\arabicsubsection. The white noise $\epsilon(s)$ follows $N(0,\sigma^{2})$ with $\sigma^{2}=0.01$. The parameters $\phi,\tau^{2},\sigma^{2}$ are unknown and need to be estimated from training data, while the Matérn smoothness parameter $\nu=3/2$ is assumed to be fixed and known. The domain is set as $\Omega=[0,1]^{2}$, with training data consisting of $1600$ grid locations in $\Omega$. We replicate the analysis on 50 datasets. For each dataset, the test set is generated as $1000$ random samples from the uniform distribution on $\Omega$.

Three methods are compared: RadGP proposed in this paper, nearest-neighbor Gaussian process with maximin ordering, and Vecchia Gaussian process predictions (V-Pred) (Katzfuss et al., 2020) with maximin ordering. The last two methods share the same training procedure but differ in how they make predictions: NNGP assumes that testing locations are independent given the training locations whereas V-Pred allows dependencies among testing locations. All methods are fitted via the MCMC algorithm described in Section \arabicsection and thus only differ in terms of their assumed graphs.

To measure the accuracy of approximations for prediction in the spatial regression model (\arabicequation), we compute the Wasserstein-2 distance between the posterior predictive distributions from each of the three approximation methods and the true Gaussian process regression model. The Wasserstein-2 distance is computed using MCMC samples. The results are shown in Figure \arabicfigure. When the complexity of the directed graph is sufficiently large, NNGP and RadGP tend to have the same approximation error while V-Pred is slightly worse. However, RadGP is able to achieve the limiting approximation accuracy with a much smaller complexity compared to the other two methods, indicating a faster approximation by RadGP. Unlike prior approximation in Section \arabicsection.\arabicsubsection, the posterior approximation error does not go to zero when the complexity of directed graphs increases. This nonzero error comes from two parts: first, the Gaussian process regression model (\arabicequation) has unknown parameters that need to be estimated from a finite amount of data; second, the Wasserstein-2 distance computed based on a finite number of MCMC samples contains Monte Carlo error.

We consider the land surface temperature data on a $499\mathrm{km}\times 581\mathrm{km}$ region southeast of Addis Ababa in Ethiopia at April 1st, 2020. The data are a product of the MODIS (Moderate Resolution Imaging Spectroradiometer) instrument on the Aqua satellite. The original data set comes as a $499\times 581$ image with each pixel corresponding to a $1\mathrm{km}\times 1\mathrm{km}$ spatial region, as shown in Fig. \arabicfigure left. However, atmosphere conditions like clouds can block earth surfaces, resulting in missing data for many pixels. There are $211,683$ spatial locations with observed temperature values and $78,236$ locations with missing values. We randomly split the training data into $3$ folds, using two folds for training and one fold for out-of-sample validation. This results in a total of $141,122$ training locations and $148,797$ testing locations.

We consider the same Gaussian process regression model as Section 5.1 with no covariates, zero mean function and exponential covariance function $K_{0}(x)=\tau^{2}\exp(-\phi x)$. The observed temperature values are centered to compensate for the lack of intercept term. We evaluate the performances of RadGP with radius $4.01\mathrm{km}$ and nearest neighbor Gaussian process with neighbor size $15$. The priors for both methods are set as $p(\phi)\propto 1$, $\tau^{2}\sim\mathrm{IG}(2,1)$ and $\sigma^{2}\sim\mathrm{IG}(1,0.1)$. Estimates for covariance parameters, out-of-sample mean squared errors and coverage for $95\%$ credible intervals are shown in Table \arabictable. The results are quite similar between the two methods.

The observed temperature and predicted temperature by RadGP are visualized in Fig. \arabicfigure. In general, the land surface temperature changes smoothly with respect to geological locations. The east and south part of the figure, which is covered by less vegetation, features slightly higher surface temperature. Small variations of temperature also occur in many local regions, indicating the geographical complexity of Ethiopia. There is a small blue plate in the west border of the figure that has significantly lower surface temperature than its surroundings, which corresponds to several lakes in the Great Rift Valley.

We further compare the performance of the two methods in terms of joint predictions at multiple held out locations. As we do not know the ground truth, we cannot calculate Wasserstein-2 error and instead focus on the joint likelihood ratio for held out data. Specifically, we select $567$ pairs of locations in the set reserved for out-of-sample validation. Each pair consists of two locations with $1$km distance from each other. Any two locations from different pairs are at least $10$km apart. For each pair of locations, we compute the logarithm of the likelihood ratio between RadGP and NNGP. The distribution of the log likelihood ratio can provide important information about joint inference. Specifically, the mean of the log likelihood ratio is $2.152$ and there are $68\%$ pairs with positive log likelihood ratios. Since the distribution of the log likelihood ratio has heavy tails, we remove the bottom and top $5\%$ of the values and visualize the rest in Fig. \arabicfigure, where the red vertical line is $\mathrm{ratio}=0$. We can see in the majority of cases, RadGP performs better than NNGP in terms of characterizing dependencies between temperature at two close locations.