A Divide-and-Conquer Bayesian Approach to Large-Scale Kriging
Distributed Kriging

Rajarshi<EMAIL_ADDRESS>Department of Statistics, University of California, Santa Cruz
Cheng<EMAIL_ADDRESS>Department of Statistics and Applied Probability, National University of Singapore
Terrance<EMAIL_ADDRESS>U. S. Bureau of Labor Statistics
Sanvesh SrivastavaT1<EMAIL_ADDRESS>Department of Statistics and Actuarial Science, The University of Iowa

T2Equally contributing authors.
T1Corresponding author.

Guhaniyogi, Li, Savitsky, Srivastava

We propose a three-step divide-and-conquer strategy within the Bayesian paradigm that delivers massive scalability for any spatial process model. We partition the data into a large number of subsets, apply a readily available Bayesian spatial process model on every subset, in parallel, and optimally combine the posterior distributions estimated across all the subsets into a pseudo posterior distribution that conditions on the entire data. The combined pseudo posterior distribution replaces the full data posterior distribution for predicting the responses at arbitrary locations and for inference on the model parameters and spatial surface. Based on distributed Bayesian inference, our approach is called “Distributed Kriging” (DISK) and offers significant advantages in massive data applications where the full data are stored across multiple machines. We show theoretically that the Bayes $L_2$-risk of the DISK posterior distribution achieves the near optimal convergence rate in estimating the true spatial surface with various types of covariance functions, and provide upper bounds for the number of subsets as a function of the full sample size. The model-free feature of DISK is demonstrated by scaling posterior computations in spatial process models with a stationary full-rank and a nonstationary low-rank Gaussian process (GP) prior. A variety of simulations and a geostatistical analysis of the Pacific Ocean sea surface temperature data validate our theoretical results.

Distributed Bayesian inference
Gaussian process
low-rank Gaussian process
modified predictive process
massive spatial data
Wasserstein distance
Wasserstein barycenter

§ INTRODUCTION

A fundamental challenge in geostatistics is the analysis of massive spatially-referenced data. Massive spatial data provide scientists with an unprecedented opportunity to hypothesize and test complex theories, see for example 22, 12, 4. This has led to the development of complex and flexible hierarchical GP-based models that are computationally intractable for a large number of spatial locations, denoted as $n$, due to the $O(n^3)$ computational cost and the $O(n^2)$ storage cost. We develop a three-step general distributed Bayesian approach, called Distributed Kriging (DISK), for boosting the scalability of any state-of-the-art spatial process model based on GP prior or its variants to multiple folds using the divide-and-conquer technique.

There is an extensive literature on scalable Gaussian process (GP)-based modeling of massive spatial data due to its great practical importance [31]. We provide a brief overview of basic ideas, deferring detailed comparisons of the existing literature with DISK to Section <ref>. A common idea in GP-based modeling is to seek dimension-reduction by endowing the spatial covariance matrix either with a low-rank or a sparse structure. Low-rank structures represent a spatial surface using a small number of  a priori chosen basis functions such that the posterior computations scale in the cubic order to the number of chosen basis functions (rather than the number of spatial locations), resulting in reduced storage and computational costs. Sparse structured models assume that the spatial correlation between two distantly located observations is nearly zero. If the assumption is true, then little information is lost by assuming independence between data at distant locations. Another approach introduces sparsity in the inverse covariance matrix using conditional independence assumptions or composite likelihoods. Some variants of dimension-reduction methods partition the spatial domain into sub-regions containing fewer spatial locations. Each of these sub-regions is modeled using a GP which are then hierarchically combined by borrowing information across the sub-regions.

The proposed DISK framework does not belong to any of these classes of methods, but it enhances the scalability of any of these methods by embedding each within the three-step DISK framework.  The outline of the DISK framework is as follows. First, the $n$ spatial locations are divided into $k$ subsets such that each subset has representative data samples from all regions of the spatial domain with the $j$th subset containing $m_j$ data samples. Second, posterior computations are implemented in parallel on the $k$ subsets using any chosen spatial process model after raising the model likelihood to a power of $n/m_j$ in the $j$th subset. The pseudo posterior distribution obtained using the modified likelihood is called the “subset posterior distribution.” Since $j$th subset posterior distribution conditions on $(m_j/n)$-fraction of the full data, the modification of the likelihood by raising it to the power of $n/m_j$ ensures that variance of each subset posterior is of the same order (as a function of $n$) as that of the full data posterior distribution. Third, the $k$ subset posterior distributions are combined into a single pseudo probability distribution, called the DISK pseudo posterior (henceforth,  DISK posterior), that conditions on the full data and replaces the computationally expensive full data posterior distribution for prediction and inference.

Our novel contributions to the growing literature on distributed Bayesian inference are two-fold. Computationally, the main innovations are in the second and third steps because the literature on general sampling and combination schemes is sparse in process-based modeling of spatial data using the divide-and-conquer technique. No restrictive data- or model-specific assumptions, such as the independence between data subsets or independence between blocks of parameters, are adopted and the DISK framework still allows principled Bayesian inference with parameter estimation, surface interpolation, and prediction. Theoretically, we provide guarantees on the accuracy of performance in estimating the true spatial surface using the DISK posterior as a function of $n$, $k$, and analytic properties of the true spatial surface.

We show that when $k$ is controlled to increase in some proper order of $n$ as $n$ tends to infinity, the Bayes $L_2$-risk of the DISK posterior achieves near minimax optimal convergence rates under different types of covariance functions.
There are some theoretical results in this direction [52, 10, 60], but DISK is the first general Bayesian framework addressing these theoretical problems with a focus on computationally efficient posterior computations in massive data applications with complex nonparametric models, while avoiding restrictive assumptions that limit wide applicability.

We illustrate the application of DISK for enhancing the scalability of two representative GP priors. One is the usual full-rank GP prior with a stationary covariance kernel and the other is a low-rank GP prior with a nonstationary covariance kernel called the modified predictive process (MPP) prior [19]. The latter prior is commonly used for estimating nonstationary surfaces in large spatial data. MPP constructs a low-rank approximation of covariance matrix for the generating distribution of the spatial surface to reduce computation time, but if the rank is moderately large, then MPP struggles to provide accurate inference in a manageable time even for $10^4$ observations. Our numerical results presented later establish that DISK with MPP prior scales to $10^6$ observations without compromising on either computational efficiency or accuracy in inference and prediction. An interesting empirical observation is that under a fixed computation budget the accuracy of MPP prior in detecting local surface features is enhanced by embedding it within the DISK framework in the sense that we are able to increase the spatial resolution.  We expect this conclusion to hold for all of the popular structured GP priors.

§.§ DISK and Existing Methods for GP-Based Modeling of Massive Spatial Data

The DISK framework does not compete with existing methods for analyzing massive spatial data, but aims to boost their scalability using the divide-and-conquer technique. With this in mind, we compare DISK with existing approaches for GP-based spatial modeling based on variants of dimension-reduction technique and refer to 31 for a more comprehensive review. Low-rank structures on the spatial covariance matrix are the most widely used tool for computationally efficient spatial computation. They represent the spatial surface using $r$  apriori chosen basis functions with associated computational complexity of $O(nr^2+r^3)$ [11, 6, 19, 26, 5, 49, 66]; however, practical considerations entail that $r$ grows roughly as $O(\sqrt{n})$ for accurate estimation, implying that $O(nr^2)$ flops are also expensive in low-rank structures. In fact, with a small $(r/n)$-ratio, scientists have observed shortcomings in many of the above methods for approximating GPs such as the propensity to oversmooth the data [58, 53]. DISK offers a solution to this problem. If $m_j \ll n$, then $(r/m_j)$-ratio is relatively large on the subsets, yielding accurate and computationally efficient inference using subset posteriors. Our theoretical results guarantee that the DISK posterior has better accuracy than any subset posterior, which can potentially outperform the full data posterior estimated using the same prior. Our simulations empirically confirm  this claim for the MPP prior.

A specific form of sparse structure uses compactly supported covariance functions to create sparse spatial covariance matrices that approximate the full covariance matrix [34, 20]. They are useful for parameter estimation and interpolation of the response (“kriging”), but not for more general inference on the latent processes due to an expensive determinant computation of the massive covariance matrix. An alternative approach is to introduce sparsity in the inverse covariance (precision) matrix of the GP likelihoods using products of lower dimensional conditional distributions [65, 48, 59], or via composite likelihoods [16, 3].
There are recent approaches, extending these ideas, that can introduce sparsity in the inverse covariance (precision) matrix of process realizations and hence enable “kriging” at arbitrary locations  [14, 28]. In related literature on computer experiments, localized approximations of GP models are proposed, see, for example,  [23]. DISK relaxes the trade-off between computation time and the accuracy in modeling a spatial surface.  In current practice, approximation methods are used with the intent to make the computations feasible at the expense of accuracy.  Reduced rank simplifications of the covariance matrix may produce over-smoothing that limits the ability to detect local features, while sparse covariance structures may underestimate correlations.  Yet, both reduced rank and sparse covariance structures may be easily embedded in our DISK framework to dramatically scale the computations such that the degree of approximation required may be notably reduced, which we demonstrate in the sequel.

The remaining variants of dimension-reduction methods combine the benefits of low-rank and sparse structure covariance functions. Examples include non-stationary models [4] and multi-level and multi-resolution models [21, 43, 32, 27]. These models usually achieve scalability by assuming block-independence at some level of the hierarchy, usually across sub-regions, but may lose scalability when they borrow information across sub-regions. Multi-resolution models are in general difficult to implement,  do not generally come with desirable theoretical guarantees concerning large sample behavior, and may become less amenable to various modification to suit different applications. In contrast, DISK makes no independence assumptions across subregions to accomplish predictions at new locations on a spatial surface and can fit a multiresolution model in each subset for enhancing its scalability. 38 proposed an approximation based on viewing a GP with Matérn
covariance as the solution to the corresponding stochastic partial differential equation, but
this approach is only applicable to covariance functions of Matérn type and may not be applicable in scaling GP with low-rank kernels.

There is a class of methods, of which DISK is a member, that divide the data into a large number of subsets, draw inference in parallel on the subsets, and combine the inferences by some mechanism that approximates the inference conditional on the full data. 7 propose combining point estimates of spatial parameters obtained from different subsamples, but they do not provide combined inference on the spatial processes or predictions. Similarly, 30 partition the spatial domain and assume independence between the data in different partitions. Although computationally attractive, assuming independence across subdomains may trigger loss in predictive uncertainty as demonstrated in [31]. In a similar effort to the DISK posterior, [25] propose drawing subset inferences and combine the posterior distributions in subsets using the idea of “meta-posterior”. This approach has an added advantage over that of 30 in that it does not assume independence across data blocks and enables prediction with accurate characterization of uncertainty [31]; however, it produces desirable inference only when a stationary GP model is fitted in each subset and is not accurate in estimation of the spatial surface when nonstationary low-rank models (e.g. MPP) are fitted in each subset. This limits the applicability of the meta-posterior. Also, 25 do not offer any theoretical guidance on choosing the number of subsets for optimal inference on the  spatial surface. The proposed DISK framework fills both these gaps. Our experiments also demonstrate that the DISK posterior provides accurate uncertainty quantification unlike
some of the divide-and-conquer approaches popularly used in the machine learning literature such as Consensus Monte Carlo [51].

Sampling algorithms are computationally inefficient in massive data settings, so this has motivated significant interest in developing general approaches to scalable Bayesian inference using the divide-and-conquer technique. The DISK framework builds on the recent works that combine the subset posterior distributions through their geometric centers, such as the mean or the median, and guarantee wide applicability under general assumptions [39, 54, 37, 40, 50, 55]. A major limitation of the current distributed approaches is that the theory and practice is limited to parametric models. By contrast, the DISK framework is tuned for accurate and computationally efficient posterior inference in nonparametric Bayesian models based on GP priors. In particular, we develop (a) a new approach to modify the likelihood for computing the subset posterior distribution of an unknown function, an infinite-dimensional parameter, (b) generalizations of existing algorithms for a full-rank and a low-rank GP prior to general MCMC samples from a subset distribution with modified likelihood, and (c) theoretical guarantees on the convergence rate of the DISK posterior to the true function, and guidance on choosing $k$ depending on the covariance function and $n$, such that the DISK posterior maintains near minimax optimal performance as $n$ tends to infinity.

The remainder of the manuscript evolves as follows. In Section 2 we outline a Bayesian hierarchical mixed model framework that incorporates models based on both the full-rank and the low-rank GP priors. Our DISK approach will work with posterior MCMC samples from such models. Section 3 develops the framework for DISK, discusses how to compute the DISK posterior distribution, and offers theoretical insights into the DISK for general GPs and their approximations. A detailed simulation study followed by an analysis of the Pacific ocean sea surface temperature data are illustrated in Section 4 to justify the use of DISK for real data. Finally, Section 5 discusses what DISK achieves, and proposes a number of future directions to explore. Proofs of the theoretical results in Section 3 are offered in the supplementary material. It also offers additional theoretical results concerning convergence rate of the DISK posterior.

§ BAYESIAN INFERENCE IN GP-BASED SPATIAL MODELS

Consider the univariate spatial regression model for the data observed at location $\sbb$ in a compact domain $\Dcal$,

\begin{align}\label{parent_proc}
y(\sbb)=\xb(\sbb)^T \betab + w(\sbb) + \epsilon(\sbb),
\end{align}

where $y(\sbb)$ and $\xb(\sbb)$ are the response and a $p\times 1$ predictor vector respectively at $\sbb$, $\betab$ is a $p\times 1$ predictor coefficient, $w(\sbb)$ is an unknown spatial function $w(\cdot)$ at $\sbb$, and $\epsilon(\sbb)$ is the white-noise process $\epsilon(\cdot)$ at $\sbb$, which is independent of $w(\cdot)$. The Bayesian implementation of the model in  (<ref>) customarily assumes  (a) that $\betab$ apriori follows N($\mub_{\betab}$, ${\boldsymbol \Sigma}_{\betab}$) and (b) that $w(\cdot)$ and $\epsilon(\cdot)$ apriori follow mean 0 GPs with covariance functions $C_{\alphab}(\sbb_1,\sbb_2)$ and $D_{\alphab}(\sbb_1,\sbb_2)$ that model  $\cov\{w(\sbb_1),w(\sbb_2)\}$ and $\cov\{\epsilon(\sbb_1), \epsilon(\sbb_2)\}$, respectively, where
$\alphab$ are the process parameters indexing the two families of covariance functions and $\sbb_1, \sbb_2 \in \Dcal$; therefore, the model parameters are $\Omegab = \{\alphab, \betab\}$.

The training data consists of predictors and responses observed at $n$ spatial locations, denoted as $\Scal = \{\sbb_1, \ldots, \sbb_n\}$.

Standard Markov chain Monte Carlo (MCMC) algorithms exist for performing posterior inference on $\Omegab$ and the values of $w(\cdot)$ at a given set of locations $\Scal^* = \{\sbb_1^*, \ldots, \sbb_l^*\}$, where $\Scal^* \cap \Scal = \emptyset$, and for predicting $y(\sbb^*)$ for any $\sbb^* \in \Scal^*$ [4]. Given $\Scal$, the prior assumptions on $w(\cdot)$ and $\epsilon(\cdot)$ imply that $\wb^T = \{w(\sbb_1), \ldots, w(\sbb_n)\}$ and $\epsilonb^T = \{\epsilon(\sbb_1), \ldots, \epsilon(\sbb_n)\}$   are independent and follow $N \left\{ \zero, \Cb({\alphab}) \right\}$ and $N \left\{ \zero, \Db({\alphab}) \right\}$, respectively, with the $(i,j)$th entries of $\Cb(\alphab)$ and $\Db(\alphab)$ are $C_{\alphab}(\sbb_i,\sbb_j)$ and $D_{\alphab}(\sbb_i,\sbb_j)$, respectively.
The hierarchy in (<ref>) is completed by assuming that $ \alphab $ apriori follows a distribution with density $\pi(\alphab)$. The MCMC algorithm for sampling $\Omegab$, $\wb^{*T} = \{w(\sbb_1^*), \ldots, w(\sbb_l^*)\}$, and $\yb^{*T} = \{y(\sbb_1^*), \ldots, y(\sbb_l^*)\}$ cycle through the following three steps until sufficient MCMC samples are drawn post convergence:

* Integrate over $\wb$ in (<ref>) and

* sample $\betab$ given  $\yb,\Xb,\alphab$ from $N(\mb_{\betab}, \Vb_{\betab})$, where

\begin{align}
\label{eq:samp-beta}
\Vb_{\betab} = \left\{ \Xb^T \Vb(\alphab)^{-1} \Xb + \Sigmab_{\betab}^{-1}   \right\}^{-1}, \quad \mb_{\betab} = \Vb_{\betab} \left\{ \Xb^{T} \Vb(\alphab)^{-1} \yb +  \Sigmab_{\betab}^{-1} \mub_{\betab}   \right\},
\end{align}

where $\Xb=[\xb(\sbb_1): \cdots :\xb(\sbb_n)]^T$ is the $n\times p$ matrix of predictors, with $p < n$, and $ \Vb(\alphab) = \Cb(\alphab) + \Db(\alphab)$; and

* sample $\alphab$ given $\yb,\Xb,\betab$ using the Metropolis-Hastings algorithm with a normal random walk proposal.

* Sample $\wb^*$ given $\yb,\Xb,\alphab,\betab$ from $N(\mb_{*}, \Vb_*)$, where

\begin{align}
\label{eq:samp-w}
\Vb_* = \Cb_{*,*}(\alphab) - \Cb_{*}(\alphab) \Vb(\alphab)^{-1} \Cb_{*}(\alphab)^T, \quad \mb_* = \Cb_{*}(\alphab) \Vb(\alphab)^{-1} (\yb - \Xb \betab),
\end{align}

$\Cb_*(\alphab)$  and $\Cb_{*,*}(\alphab)$ are $l \times n$ and $l \times l$ matrices, respectively, and the  $(i,j)$th entries of $\Cb_{*,*}(\alphab)$ and $\Cb_*(\alphab)$ are $C_{\alphab}(\sbb_i^*,\sbb_j^*)$ and $C_{\alphab}(\sbb_i^*,\sbb_j)$, respectively.

* Sample $\yb^*$ given $\alphab,\betab,\wb^*$ from $N\left\{ \Xb^* \betab + \wb^*, \Db(\alphab) \right\}$, where
$\Xb^{*T} = [\xb(\sbb_1^*) : \cdots : \xb(\sbb_l^*)]$.

Many Bayesian spatial models can be formulated in terms of (<ref>) by assuming different forms of $C_{\alphab}(\sbb_1,\sbb_2)$ and $D_{\alphab}(\sbb_1,\sbb_2)$; see [4] and supplementary material for details on the MCMC algorithm. Irrespective of the form of $\Db(\alphab)$, if no additional assumptions are made on the structure of $\Cb(\alphab)$, then the three steps require $O(n^3)$ flops in computation and $O(n^2)$ memory units in storage in every MCMC iteration. Spatial models with this form of posterior computations are based on a full-rank GP prior. In practice, if $n \geq 10^4$, then posterior computations in a model based on a full-rank GP prior are infeasible due to numerical issues in matrix inversions involving an unstructured $\Cb(\alphab)$.

There are methods which either impose a low-rank structure or a sparse structure on $\Cb(\alphab)$ to address this computational issue [4]. Methods with a low-rank structure on $\Cb(\alphab)$ expresses $\Cb(\alphab)$ in terms of $r \ll n$ basis functions (with $r=O(\sqrt{n})$ is desirable for accurate inference), in turn inducing a low-rank GP prior. Again,
a class of sparse structure uses compactly supported covariance functions to create $\Cb(\thetab)$ with overwhelming zero entries [34, 20], where as another variety of sparse structure imposes a Markov random field model on the joint distribution of $\yb$ [65, 48, 59] or $\wb$ [14, 28]. We use the MPP prior as a representative example of this broad class of computationally efficient methods. Let $\mathcal{S}^{(0)}=\{\sbb_{1}^{(0)},...,\sbb_{r}^{(0)}\}$ be a set of $r$ locations, known as the “knots,” which may or may not intersect with $\Scal$. Let $\cb(\sbb, \mathcal{S}^{(0)})= \{ C_{\alphab}(\sbb,\sbb_{1}^{(0)}), \ldots,C_{\alphab}(\sbb,\sbb_{r}^{(0)}) \}^T$ be an $r\times 1$ vector and $\Cb(\mathcal{S}^{(0)})$ be an $r\times r$ matrix whose $(i,j)$th entry is $C_{\alphab}(\sbb_{i}^{(0)},\sbb_{j}^{(0)})$. Using $\cb(\sbb_{1}, \mathcal{S}^{(0)}), \ldots, \cb(\sbb_{n}, \mathcal{S}^{(0)})$ and $\Cb(\mathcal{S}^{(0)})$, define the diagonal matrix $\deltab =\diag\{\delta(\sbb_1), \ldots, \delta(\sbb_n)\}$ with

$\delta(\sbb_i) =  C_{\alphab}(\sbb_i,\sbb_i) - \cb^T(\sbb_i, \mathcal{S}^{(0)}) \Cb(\mathcal{S}^{(0)})^{-1}\cb(\sbb_i, \Scal^{(0)})$,  $i = 1, \ldots, n$.

Let $\one({\ab=\bb}) = 1$ if $\ab=\bb$ and 0 otherwise. Then, MPP is a GP with  covariance function

\begin{align}
\label{eq:mpp-mdl-cov}
\tilde C_{\alphab}(\sbb_1,\sbb_2) = \cb^T(\sbb_1, \mathcal{S}^{(0)}) \Cb(\mathcal{S}^{(0)})^{-1}\cb(\sbb_2, \mathcal{S}^{(0)}) + \delta(\sbb_1) \one({\sbb_1=\sbb_2}), \quad \sbb_1, \sbb_2 \in \Dcal,
\end{align}

where $\tilde C_{\alphab}(\sbb_1,\sbb_2)$ depends on the covariance function of the parent GP and the selected $r$  knots, which define $\Cb(\mathcal{S}^{(0)})$, $\cb^T(\sbb_1, \mathcal{S}^{(0)})$, and $\cb^T(\sbb_2, \mathcal{S}^{(0)})$. We have used a  $\, \tilde {   } \,$ in (<ref>) to distinguish the covariance function of a low-rank GP prior from that of its parent full-rank GP. If $\tilde \Cb(\alphab)$ is a matrix with $(i,j)$th entry $\tilde C_{\alphab}(\sbb_i,\sbb_j)$, then the posterior computations using MPP, a low-rank GP prior, replace $\Cb(\alphab)$ by $\tilde \Cb(\alphab)$ in the steps 1(a), 1(b), and 2. The (low) rank $r$ structure imposed by $\Cb(\mathcal{S}^{(0)})$ implies that $\tilde \Cb(\alphab)^{-1}$ computation requires $O(nr^2)$ flops using the Woodbury formula [29]; however, massive spatial data require that $r = O(\sqrt{n})$, leading to the computational inefficiency of low-rank methods.
The next section develops our DISK framework, which uses the divide-and-conquer technique to scale the posterior computations using full-rank and low-rank GP priors.

§ DISTRIBUTED KRIGING

§.§ First step: partitioning of spatial locations

We partition the $n$ spatial locations into $k$ subsets. The value of $k$ depends on the chosen covariance function used in the spatial model, and it is set to be large enough to ensure computationally efficient posterior computations on any subset. The default partitioning scheme is to randomly allocate the locations into $k$ possibly overlapping subsets (referred to as the random partitioning scheme hereon) to ensure that each subset has representative data samples from all subregions of the domain.

Let $\Scal_j= \{\sbb_{j1}, \ldots, \sbb_{jm_j}\}$ denote the set of $m_j$ spatial locations in subset $j$ ($j=1, \ldots, k$). A spatial location can belong to multiple subsets so that  $\sum_{j=1}^{k} m_j \geq n$ but $\cup_{j=1}^k \mathcal{S}_j = \mathcal{S}$, where $\sbb_{ji} = \sbb_{i'}$ for some $\sbb_{i'} \in \Scal$ and for every $i=1, \ldots, m_j$ and $j=1, \ldots, k$. Denote the data in the $j$th partition as  $\{\yb_j, \Xb_j\}$  ($j = 1,\ldots, k$), where $\yb_j = \{y(\sbb_{j1}), \ldots, y(\sbb_{jm_j})\}^T$ is a $ m_j \times 1$ vector and $\Xb_j =
[\xb(\sbb_{j1}): \cdots :\xb(\sbb_{jm_j})]^T$ is a $m_j \times p$ matrix of predictors corresponding to the spatial locations in $\Scal_j$ with $p < m_j$.  In modern grid or cluster computing environments, all the machines in the network have similar computational power, so the performance of DISK is optimized  by choosing similar values of $m_1, \ldots, m_k$.

One can choose more sophisticated partitioning schemes than random partitioning. For example, it is possible to cluster the data based on centroid clustering [35] or hierarchical clustering based on spatial
gradients [2, 30], and then construct subsets such that each subsets contains representative data samples from each cluster. Detailed exploration later shows that even random partitioning leads to desirable inference in the various simulation settings and in the sea surface data example, hence inferential improvement with any other sophisticated partitioning should be marginal in these examples. Perhaps more sophisticated blocking methods may provide further improvement in the cases where spatial locations are drawn based on specific designs; for example, sophisticated partitioning schemes have inferential benefits when a sub-domain shows substantial local behavior compared to the others [27], or sampled locations are chosen based on a specific survey design. Since they are atypical examples in the spatial context, we will pursue them elsewhere in greater detail.

The univariate spatial regression models using either a full-rank or a low-rank GP prior for the data observed at any location $\sbb_{ji} \in \Scal_j \subset \Dcal$ is given by

\begin{align}\label{eq:mdl-j}
y(\sbb_{ji}) = \xb(\sbb_{ji})^T \betab + w(\sbb_{ji}) + \epsilon(\sbb_{ji}), \quad i=1, \ldots, m_j.
\end{align}

Let $\wb_j^T = \{w(\sbb_{j1}), \ldots, w(\sbb_{jm_j})\}$ and $\epsilonb_j^T = \{\epsilon(\sbb_{j1}), \ldots, \epsilon(\sbb_{jm_j})\}$ be the realizations of GP $w(\cdot)$ and white-noise process $\epsilon(\cdot)$, respectively, in the $j$th subset.
After marginalizing over $\wb_j$ in the GP-based model for the $j$th subset, the likelihood of $\Omegab = \{\alphab, \betab\}$ is given by

$\ell_{j}(\Omegab) = N \{\yb_j \mid \Xb_j \betab, \Vb_j(\alphab)\}$,

where $\Vb_j(\alphab) = \Cb_j(\alphab) + \Db_j(\alphab)$ and  $\Vb_j(\alphab) = \tilde \Cb_j(\alphab) + \Db_j(\alphab)$ for full-rank and low-rank GP priors, respectively, and $\Cb_j(\alphab), \tilde \Cb_j(\alphab), \Db_j(\alphab)$ are obtained by extending the definitions of $\Cb(\alphab), \tilde \Cb(\alphab), \Db(\alphab)$ to the $j$th subset. In a model based on full-rank or low-rank GP prior, the likelihood of $\wb_j$ given $\yb_j$, $\Xb_j$, and $\Omegab$  is

$\ell_{j }(\wb_j) = N \{\yb_j - \Xb_j \betab \mid \wb_j, \Db_j(\alphab) \}.$

The likelihoods in $\ell_{j}(\Omegab)$ and $\ell_{j}(\wb_j)$ are used to define the
posterior distributions for $\betab, \alphab, \wb^*$,  $\yb^*$ based on a full-rank or a low-rank GP prior in subset $j$ and are called $j$th subset posterior distributions.

§.§ Second step: sampling from subset posterior distributions

We define subset posterior distributions by modifying the likelihoods in $\ell_{j}(\Omegab)$ and $\ell_{j}(\wb_j)$. More precisely, the density of the $j$th subset posterior distribution of $\Omegab$ is given by

\begin{align}\label{eq:fullsub-j}
\pi_{m_j} (\Omegab \mid \yb_j) &= \frac{\{\ell_{j }(\Omegab) \}^{n/m_j}\pi(\Omegab) } {\int \{\ell_{j }(\Omegab) \}^{n/m_j}\pi(\Omegab) d \Omegab},
\end{align}

where we assume that $\int \{\ell_{j }(\Omegab) \}^{n/m_j}\pi(\Omegab) d \Omegab < \infty$,
and the subscript `$m_j$' denotes that the density conditions on $m_j$ data samples in the $j$th subset. The modification of likelihood to yield the subset posterior density in (<ref>) is called stochastic approximation [39]. Raising the likelihood to the power of $n/m_j$ is equivalent to replicating every $y(\sbb_{ji})$ $n/m_j$ times ($i=1, \ldots, m_j$), so stochastic approximation accounts for the fact that the $j$th subset posterior distribution conditions on a $(m_j/n)$-fraction of the full data and ensures that its variance is of the same order (as a function of $n$) as that of the full data posterior distribution. Unlike parametric models, stochastic approximation in spatial regression models has not been studied previously in the literature. We address this gap next.

With the proposed stochastic approximation in (<ref>), the full conditional densities of $j$th subset posterior distributions for prediction and inference follow from their full data counterparts. The $j$th full conditional densities of $\betab$ and $\alphab$ in the GP-based models are

\begin{align}
\label{eq:sub-j}
\pi_{m_j} (\betab \mid \yb_j, \alphab) &= \frac{\{\ell_{j }(\Omegab) \}^{n/m_j}\pi(\betab) } {\int \{\ell_{j }(\Omegab) \}^{n/m_j}\pi(\betab) d \betab}, \quad
\pi_{m_j} (\alphab \mid \yb_j, \betab) = \frac{\{\ell_{j}(\Omegab) \}^{n/m_j}\pi(\alphab) } {\int \{\ell_{j }(\Omegab) \}^{n/m_j}\pi(\alphab) d \alphab},
\end{align}

where $\pi(\betab) = N(\mub_{\betab}, {\boldsymbol \Sigma}_{\betab})$, $\pi(\alphab)$ is the prior density of $\alphab$, and we assume that $\int \{\ell_{j}(\Omegab) \}^{n/m_j}\pi(\betab) d \betab$
and $\int \{\ell_{j }(\Omegab) \}^{n/m_j}\pi(\alphab) d \alphab$ respectively are finite.
The $j$th full conditional densities of $\yb^*$ and $\wb^*$ are calculated after modifying the likelihood of $\wb_j$ using stochastic approximation. Given $\yb_j$, $\Xb_j$, and $\Omegab$, straightforward calculation yields that the $j$th subset posterior predictive density of $\wb^*$ is $\pi_{m_j} (\wb^* \mid \yb_j, \Omegab) = N(\wb^* \mid \mb_{j*}, \Vb_{j*})$, with

\begin{align}
\label{eq:sub-w-j}
&\Vb_{j*} = \Cb_{*,*}(\alphab) - \Cb_{*j}(\alphab) \Vb_j(\alphab)^{-1} \Cb_{*j}(\alphab)^T, \quad \mb_{j*} = \Cb_{*j}(\alphab) \Vb_j(\alphab)^{-1} (\yb_j - \Xb_j \betab),%\nonumber\\
%\tilde \Vb_{j*} &= \tilde \Cb_{*,*} - \tilde \Cb_{*j} \tilde \Vb_j(\alphab)^{-1} \tilde \Cb_{*j}^T, \quad \tilde \mb_{j*} = \tilde \Cb_{*j} \tilde \Vb_j(\alphab)^{-1} (\yb_j -\Xb_j \betab),  \nonumber\\
%&\qquad\qquad\qquad\qquad\qquad \Vb_j(\alphab) = \Cb_j(\alphab) + k^{-1} \Db_j(\alphab), %\quad \tilde \Vb_j(\alphab) = \tilde \Cb_{j} + k^{-1} \Db_j(\alphab),
\end{align}

where $\Vb_j(\alphab) = \Cb_j(\alphab) + (n/m_j)^{-1} \Db_j(\alphab)$ and $\Vb_j(\alphab) = \tilde \Cb_j(\alphab) + (n/m_j)^{-1} \Db_j(\alphab)$ for full-rank and low-rank GP priors, respectively,  and
$\Cb_{*,*}(\alphab), \Cb_{*j}(\alphab)$ are $l \times l$, $l \times m_j$ matrices obtained by extending the definition in (<ref>) to subset $j$ for full-rank and low-rank GP priors with covariance functions $C_{\alphab}(\cdot, \cdot)$ and $\tilde C_{\alphab}(\cdot, \cdot)$, respectively. We note that the stochastic approximation exponent, $n/m_j$, scales $\Db_j(\alphab)$ in $\Vb_j(\alphab)$
so that the uncertainty in subset and full data posterior distributions are of the same order (as a function of $n$). The $j$th subset posterior predictive density of $\yb^*$ given the MCMC samples of $\wb^*$  and $\Omegab$ in the $j$th subset is
$N\{\yb^* \mid \Xb^* \betab + \wb^*, \Db_j(\alphab)\}$.
We employ the same three-step sampling algorithm, as earlier introduced, specialized to subset $j$  ($j = 1,\ldots,k$), sampling $\{\betab, \alphab, \yb^*, \wb^*\}$ in each subset across multiple MCMC iterations; see supplementary material for detailed derivations of subset posterior sampling algorithms in the full-rank and low-rank GP priors. The computational complexity of $j$th subset posterior computations follows from their full data counterparts if we replace $n$ by $m_j$. Specifically, the computational complexities
for sampling a subset posterior distribution are $O(m^3)$ and $O(mr^2)$ flops per iteration if the model in (<ref>) uses a full-rank or a low-rank GP prior, respectively, where $m = \max_{j} m_j$. Performing subset posterior computations in parallel across $k$ servers also alleviates the need to store large covariance matrices.

The combination of subset posteriors outlined below is more widely applicable compared to other divide-and-conquer type approaches as it is free of model- or data-specific assumptions, such as independence of samples in training data, except that every subset posterior distribution has a density and has finite second moments.

§.§ Third step: combination of subset posterior distributions

The combination step relies on the notion of Wasserstein barycenter, as used in some related scalable Bayes literature for independent data [37, 55]. We first provide some background on this topic. Let $(\Theta, \rho)$ be a complete separable metric space and $\Pcal(\Theta)$ be the space of all probability measures on $\Theta$. The Wasserstein space of order $2$ is a set of probability distributions defined as
$\Pcal_2(\Theta) = \{\mu \in \Pcal(\Theta): \int_{\Theta} \rho^2(\theta, \theta_0) \mu(d \theta) < \infty \}$,
where $\theta_0 \in \Theta$ is arbitrary and $\Pcal_2(\Theta) $ does not depend on the choice of $\theta_0$. The Wasserstein distance of order 2, denoted as $W_2$, is a metric on $\Pcal_2(\Theta)$. Let $\mu, \nu$ be two probability measures in $\Pcal_2(\Theta)$ and $\Pi (\mu, \nu)$ be the set of all probability measures on $\Theta \times \Theta$ with marginals $\mu$ and $\nu$, then $W_2$ distance between $\mu$ and $\nu$ is defined as
$W_2(\mu, \nu) = \{\underset{\pi \in \Pi (\mu, \nu)} {\mathrm{inf}} \, \int_{\Theta \times \Theta} \rho^2(x, y) \, d \pi(x, y) \}^{1/2}.$
Let $\nu_1, \ldots, \nu_k \in \Pcal_2(\Theta)$, then the Wasserstein barycenter of $\nu_1, \ldots, \nu_k$ is defined as

\begin{align}
\overline \nu = \underset{ \nu \in \Pcal_2(\Theta)} {\argmin} \, \frac{1}{k} \sum_{j=1}^{k} W_2^2 (\nu, \nu_j). \label{eqn:w-bary}
\end{align}

It is known that $\overline \nu$ exists and is unique [1].

In the DISK framework, for any parameter of interest $\theta$, either a scalar or a vector, the DISK posterior is defined to be the Wasserstein barycenter of the $k$ subset posterior distributions of $\theta$. Here, $\theta$ can be taken as $\betab$, $\alphab$, $\wb^*$, $\yb^*$, their individual components, or any functionals of these parameters. In other words, for our DISK approach, $\nu_1, \ldots, \nu_k$ in (<ref>) are taken as the $k$ subset posterior distributions of $\theta$. Hence the DISK posterior, mathematically computed from the Wasserstein barycenter $\overline \nu$ in (<ref>), provides a general notion of obtaining the mean of $k$ possibly dependent subset posterior distributions. For Bayesian inference, the exact subset posteriors of $\theta$ ($\nu_1, \ldots, \nu_k$ in (<ref>)) are analytically intractable in general, but they can be well approximated by the subset posterior MCMC samples of $\theta$, and we can conveniently estimate the empirical version of the Wasserstein barycenter $\overline \nu$ by efficiently solving a sparse linear program as described in [13, 54, 56]. It has been shown that for independent data, the Wasserstein barycenter is a preferable choice to several other combination methods [37, 55]; for example, directly averaging over many subset posterior densities with different means can usually result in an undesirable multimodal pseudo posterior distribution, but the Wasserstein barycenter does not have this problem and can recover a unimodal posterior; see, for example, Figure 1 in 55. Besides, it does not rely on the asymptotic normality of
the subset posterior distributions as in other approches, such as consensus Monte Carlo [51].

If $\theta$ represents a one-dimensional functional of interest (a functional of $\betab$, $\alphab$, $\wb^*$, or $\yb^*$), then the DISK posterior of $\theta$ can be easily obtained by averaging empirical subset posterior quantiles [37]. This is because the $W_2$ distance between two univariate distributions is the same as the $L_2$ distance between their quantile functions (Lemma 8.2 of 9). In particular, let $\nu$ and $\nu_j$ be the full posterior and $j$th subset posterior distribution of $\theta$, and $\overline \nu$ be the Wasserstein barycenter of $\nu_1, \ldots, \nu_k$ as in (<ref>). For any $q\in (0,1)$, let $\hat {\nu}_{j}^q$ be the $q$th empirical quantile of $\nu_j$ based on the MCMC samples from $\nu_j$, and $\hat {\overline \nu}^q$ be the $q$th quantile of the empirical version of $\overline \nu$. Then, $\hat {\overline \nu}^q$ can be computed as

\begin{align}
\label{eq:emp-quant}
\hat {\overline \nu}^q = \frac{1}{k} \sum_{j=1}^k \hat {\nu}_{j}^q, \quad q = \xi, 2\xi, \ldots, 1 - \xi,
\end{align}

where $\xi$ is the grid-size of the quantiles [37]. If the $\xi$-grid is fine enough in (<ref>), then the parameter MCMC samples from the marginal DISK distribution are obtained by inverting the empirical distribution function supported on the quantile estimates.

In practice, the primary interest often lies in the marginal distributions of model parameters and predicted values; that is, the posterior of some one-dimensional functional $\theta$. Therefore, the univariate Wasserstein barycenter obtained by averaging quantiles in (<ref>) accomplishes this with great generality and convenient implementation. For this reason, in the following sections, we only focus on the case where $\theta$ is one-dimensional and use (<ref>) to compute the DISK posterior through its empirical quantiles. Nonetheless, we emphasize that the DISK posterior for a multivariate $\theta$ can still be efficiently computed using the sparse linear program for Wasserstein barycenters as described in 13, 54, 56; however,  these methods are computationally expensive and lead to only a small amount of the improvement over the univariate quantile combination in (<ref>) as revealed by our experiments.

A key feature of the DISK combination scheme is that given the subset posterior MCMC samples, the combination step is agnostic to the choice of a model. Specifically,  given MCMC samples from the $k$ subset posterior distributions,
(<ref>) remains the same for models based on a full-rank GP prior, a low-rank GP prior, such as MPP, or any other model described in Section <ref>. Since the averaging over $k$ subsets takes $O(k)$ flops and $k < n$, the total time for computing the empirical quantile estimates of the DISK posterior in inference or prediction requires $O(k) + O(m^3)$ and $O(k) + O({rm^2})$  flops in models based on full-rank and low-rank GP priors, respectively. Assuming that we have abundant computational resources, $k$ is chosen large enough so that $O(m^3)$ computations are feasible. This would enable applications of the DISK framework in models based on both full-rank and low-rank GP priors in massive $n$ settings.

§.§ Bayes $L_2$-risk of DISK: convergence rates and the choice of $k$

In the divide-and-conquer Bayesian setup, it is already known that when the data are independent and identically distributed (i.i.d.), the combined posterior distribution using the Wasserstein barycenter of subset posteriors approximates the full data posterior distribution at a near optimal parametric rate, under certain conditions as $n, k,  m_1, \ldots, m_k \rightarrow \infty$ [37, 55]; however, in models based on spatial process, data are not i.i.d. and inference on the infinite dimensional true spatial surface is of primary importance. Few formal theoretical results are available in this nonparametric divide-and-conquer Bayes setup. A notable exception is the recent paper [60], which shows that combination using Wasserstein barycenter has optimal Bayes risk and adapts to the smoothness of $w_0(\cdot)$, the true but unknown $w(\cdot)$, in the Gaussian white noise model. The Gaussian white noise model is a special case of (<ref>) with additional smoothness assumptions on $w_0(\cdot)$.

We investigate the theoretical properties of the DISK predictive posterior of the mean surface $\xb(\cdot)^T \betab + w(\cdot)$. For ease of presentation, we assume that $m_1 = \cdots = m_k=m$. Determining the appropriate order for $k$ in terms of $n$ is one of the key issues for all divide-and-conquer statistical methods. If $k$ is too small, then the biases in subset posterior distributions are small due to a large subset size $m$, while the overall variance of the DISK posterior is large due to too few subsets. In contrast, if $k$ is very large, then the biases in subset posterior distributions are large due to a small subset size $m$, while the variance of the DISK posterior can be small due to the large number of subsets.
Ideally, $k$ should be controlled to increase in some order of $n$, such that the bias and the variance can be balanced and the Bayes $L_2$-risk of the DISK posterior can be minimized.

We formally explain the model setup for our theory development. Suppose that the data generation process follows the model (<ref>) with the true parameter value $\Omegab_0=(\alphab_0,\betab_0)$ and the true spatial surface $w_0(\cdot)$. We focus on the Bayes $L_2$-risk of the DISK predictive posterior for the mean function in (<ref>); that is, $\xb(\sbb^*)^T \betab + w(\sbb^*)$ for any testing location $\sbb^* \in \Scal$. To ease the complexity of our theory, we first present two theorems below for the simplified model

\begin{align} \label{theory-mdl}
y(\sbb_i)= w(\sbb_i) + \epsilon(\sbb_i), \quad \epsilon(\sbb_i) \sim N(0, \tau^2), \quad i = 1, \ldots, n, \quad  w(\cdot) \sim \text{GP}\{0, C_{{\alphab}}(\cdot, \cdot) \}.
\end{align}

Compared to the spatial model (<ref>), the model (<ref>) does not contain the regression term $\xb(\sbb)^T \betab$; however, our theory includes this regression term later by modifying the covariance function; see our Corollary <ref> below. Our theoretical setup is a general one that subsumes GP priors with Matern covariance kernel [57] and the wide class of low-rank GP priors. While $\tau^2$ and $\alphab$ are unknown in practice and are assigned prior distributions, our setup in (<ref>) assumes that $\tau^2$ and $\alphab$ are known because this is a common assumption in the asymptotic study of GP-based regression; see 61 and the references therein. Furthermore, it is known that the theoretical results obtained by assuming $\tau^2$ and $\alphab$ as known constants are equivalent to those obtained by assigning priors with bounded supports on these parameters [62, 64]. Generalization to unboundedly supported priors is well known to be a difficult in Bayesian GP theory, especially for $\alphab$, and only partial solutions are available [64, 68].

We introduce some definitions used in stating the results in this section. Let $\PP_{\sbb}$ be a design distribution of $\sbb$ over $\Dcal$, $L_2(\PP_{\sbb})$ be the $L_2$ space under $\PP_{\sbb}$, the inner product in $L_2(\PP_{\sbb})$ is defined as $\langle f, g\rangle_{L_2(\PP_{\sbb})} = \EE_{\PP_{\sbb}} (fg)$ for any $f, g \in L_2(\PP_{\sbb})$. For any $f\in L_2(\PP_{\sbb})$ and $\sbb \in \Dcal$, define the linear operator $(T_{\alphab}f)(\sbb)=\int_{\Dcal} C_{\alphab}(\sbb,\sbb')f(\sbb')d\PP_{\sbb}(\sbb')$. According to the Mercer's theorem, there exists an orthonormal basis $\left\{\phi_i(\sbb)\right\}_{i=1}^\infty$ in $L_2(\PP_{\sbb})$, such that $C_{\alphab}(\sbb, \sbb') = \sum_{i=1}^{\infty} \mu_i \phi_i(\sbb) \phi_i(\sbb')$, where $\mu_1\geq \mu_2 \geq \ldots\geq 0$ are the eigenvalues and $\{\phi_i(\sbb)\}_{i=1}^\infty$ are the eigenfunctions of $T_{\alphab}$.
The trace of the kernel $C_{\alphab}$ is defined as $\tr(C_{\alphab})=\sum_{i=1}^{\infty} \mu_i$.
Any $f\in L_2(\PP_{\sbb})$ has the series expansion $f(\sbb) = \sum_{i=1}^{\infty} \theta_i \phi_i(\sbb)$, where $\theta_i = \langle f,\phi_i\rangle_{L_2(\PP_{\sbb})} $. The reproducing kernel Hilbert space (RKHS) $\HH$ attached to $C_{\alphab}$ is the space of all functions $f \in L_2(\PP_{\sbb})$ such that the $\HH$-norm $\|f\|_{\HH} = \sum_{i=1}^{\infty} \theta_i^2 / \mu_i<\infty$. The RKHS $\HH$ is the completion of the linear space of functions defined as $\sum_{i=1}^I a_i C_{\alphab}(\sbb_i, \cdot)$, where $I$ is a positive integer, $\sbb_i \in \Dcal$, and  $ a_i \in \RR$ ($i=1, \ldots, I$); see [63] for more details on RKHS.

For theory development, we consider a random design scheme with the observed locations $\Scal=\left\{\sbb_1,\ldots,\sbb_n\right\}$ and $\Scal^*=\{\sbb^*\}$. We assume that the locations $\sbb_1, \ldots, \sbb_n, \sbb^*$ are independently drawn from the same sampling distribution $\PP_{\sbb}$. We further impose the following assumptions.

*  (RKHS) The true function $w_0$ is an element of the RKHS $\HH$ attached to the kernel $C_{\alphab}$.

*  (Trace class kernel) $\tr(C_{\alphab})<\infty$.

*  (Moment condition) There are positive constants $\rho$ and, with a slight abuse of notation, $r \geq 2$ such that $\EE_{\PP_{\sbb}} \{\phi_i^{2r}(\sbb)\}  \leq \rho^{2r}$ for every $i$, and $\text{var} \left\{\epsilon(\sbb)\right\} = \tau^2 <\infty$ for any $\sbb\in \Dcal$.

In general, the RKHS $\HH$ can be a much smaller space relative to the support of the GP prior. While we use $w_0\in \HH$ in Assumption <ref> mainly for technical simplicity, this assumption can be possibly relaxed by considering sieves with increasing $\HH$-norms; in the same vein as Assumption B$'$ and Theorem 2 in 72. We expect that similar convergence rate results to our Theorems <ref> and <ref> can be shown for much larger classes of functions than $\HH$; see the additional posterior convergence theory in Section 2 of supplementary material. In Assumption <ref>, $\tr(C_{\alphab})$ measures the size of the covariance kernel and imposes conditions on the regularity of functions that DISK can learn. Assumption <ref> controls the error in approximating $C_{\alphab}(\sbb, \sbb')$ by a finite sum, and the superscript $r$ here should not be confused with the number of knots in low-rank GP priors. Our results are valid for any error distribution with a finite variance, and it holds trivially for the model in (<ref>).

We begin by examining the Bayes $L_2$-risk of the DISK posterior for estimating $w_0$ in (<ref>). Under the setup of (<ref>), let $\EE_{\sbb^*}$, $\EE_{0}$, $\EE_{\Scal}$, and $\EE_{0 \mid \Scal}$  respectively be the expectations with respect to the distributions of $\sbb^*$, $(\Scal, \yb)$, $\Scal$, and $\yb$ given $\Scal$. If $\overline w(\sbb^*)$ is a random variable that follows the DISK posterior for estimating $w_0(\sbb^*)$, then $\overline w(\sbb^*)$ has the density $N(\overline m, \overline v)$, where

\begin{align}
\label{eq:e3}
\overline m = \frac{1}{k} \sum_{j=1}^k \cb_{j,*}^T (\Cb_{j, j} + \tfrac{\tau^2}{k} \Ib)^{-1}\yb_j , \;
\overline v^{1/2} = \frac{1}{k}\sum_{j=1}^k v_j^{1/2}, \;
v_j = c_{*,*} - \cb_{j,*}^T (\Cb_{j,j} + \tfrac{\tau^2}{k} \Ib)^{-1} \cb_{j, *},
\end{align}

$c_{*,*} = \cov\{w(\sbb^*), w(\sbb^*)\}$, and $\cb_{j,*}^T = [\cov\{w(\sbb_{j1}), w(\sbb^*)\}, \ldots, \cov\{w(\sbb_{jm}), w(\sbb^*)\}]$. The Bayes $L_2$-risk of the DISK posterior
in estimating $w_0$ is $\EE_{0} \left[ \EE_{\sbb^*} \{\overline w(\sbb^*) - w_0(\sbb^*)\}^2 \right]$. This risk can be used to quantify how quickly the DISK posterior concentrates around the unknown true surface $w_0(\cdot)$ as the total sample size $n$ increases to infinity. It is straightforward to show (see the proof of Theorem <ref> in the Supplementary Material) that this Bayes $L_2$-risk can be decomposed into the squared bias, the variance of subset posterior means, and the variance of DISK posterior terms as

\begin{align}
&\text{bias}^2 =  \EE_{\sbb^*} \EE_{\Scal} \{  \cb_{*}^T(k \Lb + \tau^2 \Ib)^{-1}\wb_0 - w_0(\sbb^*)\}^2, \quad \nonumber \\
&\text{var}_{\text{mean}} =  \tau^2 \EE_{\sbb^*} \EE_{\Scal} \left\{\cb^T_* (k \Lb + \tau^2 \Ib)^{-2} \cb_* \right\}, \quad \nonumber \\
&\text{var}_{\text{DISK}} =  \EE_{\sbb^*}\EE_{\Scal}  (\overline v(\sbb^*)), \label{asymp-mse}
\end{align}

where  $\overline v(\sbb^*) = \EE_{0 \mid \Scal} \left[ \var \{\overline w(\sbb^*) \mid \yb \} \right]$, $\cb^T_* = (\cb^T_{1,*}, \ldots, \cb^T_{k,*})$,  $\wb_{0j} = \{w_0(\sbb_{j1}), \ldots, w_0(\sbb_{jk})\}$ ($j=1, \ldots, k$), $\wb_0^T = (\wb_{01}, \ldots, \wb_{0k})$, and $\Lb$ is a block-diagonal matrix with $\Cb_{1,1}, \ldots, \Cb_{k,k}$ along the diagonal. The next theorem provides theoretical upper bounds for each of the three terms in (<ref>).

If Assumptions <ref>–<ref> hold, then

\begin{align}
\label{eq:bayesrisk-ineq}
\text{Bayes } L_2\text{ risk} &= \EE_{\Scal} \EE_{0|\Scal} \EE_{\sbb^*} \{\overline w(\sbb_*) - w_0(\sbb_*)\}^2 =  \text{bias}^2 + \text{var}_{\text{mean}} + \text{var}_{\text{DISK}} ,\nonumber \\
\text{bias}^2 & \leq \frac{8\tau^2}{n} \| w_0 \|_{\HH}^2 + \| w_0 \|_{\HH}^2 \; \underset{d \in \NN}{\inf} \, \left[ \frac{8n}{\tau^2} \rho^4 \tr(C_{\alphab}) \tr(C_{\alphab}^{d}) + \tr(C_{\alphab}) R(m,n,d,r) \right],\nonumber \\
\text{var}_{\text{mean}} &\leq \frac{2 n + 4 \| w_0 \|^2_{\HH} }{k} \inf_{d\in \NN} \left[\left\{ \frac{12n}{\tau^2} \rho^4 \tr(C_{\alphab}) +1\right\}\tr(C_{\alphab}^d)  + R(m,n,d,r) \right]  \nonumber \\
&~~~ +
12 \frac{\tau^2}{n} \gamma \left( \frac{\tau^2}{n} \right) +   \frac{12}{k} \frac{\tau^2}{n} \| w_0 \|_{\HH}^2 ,  \nonumber \\
\text{var}_{\text{DISK}} &\leq 3 \frac{\tau^2}{n} \gamma \left( \frac{\tau^2}{n} \right) + \underset{d \in \NN}{\inf} \,\left[ \left\{ \frac{4n}{\tau^2}\tr(C_{\alphab}) + 1 \right\} \tr(C_{\alphab}^d) + \tr(C_{\alphab}) R(m,n,d,r) \right],
\end{align}

where $\NN$ is the set of all positive integers, $A$ is a global positive constant that does not depend on any of the quantities here, and

\begin{align*}
& b(m, d, r) = \max \left( \sqrt{\max(r, \log d)}, \frac{\max(r,\log d)}{m^{1/2 - 1/r}} \right), \\
& R(m,n,d,r)= \left\{  \frac{A\rho^2  b(m, d, r) \gamma(\tfrac{\tau^2}{n})}{\sqrt{m}} \right\}^r,\\
& \gamma(a) = \sum_{i=1}^{\infty} \frac{\mu_i}{\mu_i+a} \text{ for any } a>0, \quad \tr(C^d_{\alphab}) = \sum_{i=d+1}^{\infty} \mu_i.
\end{align*}

These bounds are similar to the bounds obtained in Theorem 1 of 72 for the frequentist divide-and-conquer estimator in kernel ridge regression. Although the upper bounds in (<ref>) appear very complicated and involve many terms, the dominant term among them is $\frac{\tau^2}{n} \gamma \left( \frac{\tau^2}{n} \right)$, where the function $\gamma(\cdot)$ is related to the “effective dimensionality" of the covariance kernel $C_{\alphab}$ [71]. This term determines how fast the Bayes $L_2$-risk converges to zero, as long as $k$ is chosen to be some proper order of $n$ such that all the other terms in the upper bounds of (<ref>) can be made negligible compared to $\frac{\tau^2}{n} \gamma \left( \frac{\tau^2}{n} \right)$.

In contrast to the frequentist literature such as 72, a significant difference in our Theorem <ref> is that our risk bounds involve two different variance terms. While our analysis naturally introduces the variance term $\text{var}_{\text{DISK}}$ that corresponds to the variance of the DISK posterior distribution, any frequentist kernel regression method only finds a point estimate of $w_0$ and thus does not include this variance term.  As a by-product of the proof of Theorem <ref>, the upper bound for $\text{var}_{\text{DISK}}$ can be used to show that the integrated predictive variance of GP decreases to zero as the subset sample size $m\to\infty$ for various types of covariance kernels. A closely related work by 24 studies the asymptotic behavior for the mean squared error of GP, but unrealistically assumes that the error variance increases with the sample size, which prevents their predictive variance of GP from converging to zero.

Each of the three upper bounds in Theorem <ref> can be made close to zero as $n$ increases to $\infty$ and $k$ is chosen to appropriately depend on $n$. The next theorem finds the appropriate order for $k$ in terms of $n$, such that the DISK posterior achieves nearly minimax optimal rates in its Bayes $L_2$-risk (<ref>), for three types of commonly used covariance kernels, (i) degenerate covariance kernels, (ii) covariance kernels with exponentially decaying eigenvalues, and (iii) covariance kernels with polynomially decaying eigenvalues.
The covariance kernel $C_{\alphab}$ is a degenerate kernel of rank $d^*$ if there is some constant positive integer $d^*$ such that $\mu_1\geq \mu_2 \geq \ldots \geq \mu_{d^*}>0$ and $\mu_{d^*+1}=\mu_{d^*+2}=\ldots=\mu_{\infty}=0$. The covariance kernels in subset of regressors approximation [45] and predictive process [6] are degenerate with their ranks equaling the number of “inducing points” and knots, respectively. The squared exponential kernel is very popular in machine learning. Its RKHS belongs to the class of RKHSs of kernels with exponentially decaying eigenvalues. Similarly, the class of RKHSs of kernels with polynomially decaying eigenvalues includes the Sobolev spaces with different orders of smoothness and the RKHS of the Matérn kernel. This kernel is most relevant for spatial applications, but we provide the other two results for a more general audience.

If Assumptions <ref>–<ref> hold and $r>4$ in Assumption <ref>, then, as $n \rightarrow \infty$,

(i) if $C_{\alphab}$ is a degenerate kernel of rank $d^*$ and $k \leq c n^{\frac{r-4}{r-2}}/(\log n)^{\frac{2r}{r - 2}}$ for some constant $c>0$, then the Bayes $L_2$-risk of DISK posterior satisfies $\EE_{\sbb^*} \EE_{\Scal} \EE_{0|\Scal} \{\overline w(\sbb^*) - w_0(\sbb^*)\}^2 = O\left( n^{-1}\right)$;

(ii) if $\mu_i \leq c_{1\mu} \exp\left(-c_{2\mu} i^\kappa \right)$ for some constants $c_{1\mu}>0, c_{2\mu}>0, \kappa>0$ and all $i\in \NN$, and for some constant $c>0$, $k\leq c n^{\frac{r-4}{r-2}}/(\log n)^{\frac{2(r\kappa+r-1)}{\kappa(r-2)}}$, then the Bayes $L_2$-risk of DISK posterior satisfies $\EE_{\sbb^*} \EE_{\Scal} \EE_{0|\Scal} \{\overline w(\sbb^*) - w_0(\sbb^*)\}^2 = O\left\{ (\log n)^{1/\kappa} /n \right\}$; and

(iii) if $\mu_i \leq c_{\mu} i^{-2\nu}$ for some constants $c_{\mu}>0,\nu>\tfrac{r-1}{r-4}$ and all $i\in \NN$, and for some constant $c>0$, $k\leq c n^{\frac{(r-4)\nu-(r-1)}{(r-2)\nu}} / (\log n)^{\frac{2r}{r-2}}$, then the Bayes $L_2$-risk of DISK posterior satisfies $\EE_{\sbb^*} \EE_{\Scal} \EE_{0|\Scal} \{\overline w(\sbb^*) - w_0(\sbb^*)\}^2= O\left( n^{-\frac{2\nu-1}{2\nu}}\right)$.

The rate of decay of the $L_2$-risks in (i) and (ii) with $\kappa=2$ are known to be minimax optimal [47, 69], whereas the rate of decay of the $L_2$-risk in (iii) is slightly larger than the minimax optimal rate by a factor of $n^{\frac{1}{2 \nu (2 \nu + 1)}}$. Since $\Dcal$ is compact in all spatial applications, $r$ in Assumption <ref> can be chosen as large as possible. Letting $r\rightarrow \infty$, the upper bounds on $k$ in (i), (ii), and (iii) reduce to $k = O\{n / (\log n)^{2} \} $, $k = O\{n / (\log n)^{2/\kappa}\} $, and $k= O\{n^{\frac{\nu - 1}{\nu}} / (\log n)^2\}$, respectively.

Now we generalize the results in Theorems <ref> and <ref> to the model (<ref>). Besides <ref>–<ref>, we further make the following assumption on $\xb(\cdot)$ and the prior on $\betab$:

*  All $p$ components of $\xb(\cdot)$ are non-random functions in $\Scal$. The prior on $\betab$ is $N(\mub_{\betab}, \Sigmab_{\betab})$ and it is independent of the prior on $w(\cdot)$, which is $\text{GP}\{0,C_{\alphab}(\cdot,\cdot)\}$.

By the normality and joint independence in Assumption <ref>, it is straightforward to show that the mean function $\xb(\sbb)^T \betab + w(\sbb)$ has a GP prior $\text{GP}\{\xb(\cdot)^T \mub_{\betab}, \check C_{\alphab}(\cdot,\cdot)\}$, where the modified covariance kernel $\check C_{\alphab}$ is given by

\begin{align} \label{modified-kernel}
\check C_{\alphab}(\sbb_1, \sbb_2) &= \cov \left\{\xb(\sbb_1)^T \betab + w(\sbb_1), ~ \xb(\sbb_2)^T \betab + w(\sbb_2) \right\} \nonumber \\
&= \xb (\sbb_1)^T \Sigmab_{\betab}  \xb (\sbb_2) + C_{\alphab}(\sbb_1, \sbb_2),
\end{align}

for any $\sbb_1,\sbb_2\in \Scal$. With this modified covariance kernel, we have the following corollary:

If Assumption <ref> holds, and Assumptions <ref>–<ref> also hold with all $C_{\alphab}$ replaced by $\check C_{\alphab}$ defined in (<ref>), then the conclusions of Theorems <ref> and <ref> hold for the Bayes $L_2$-risk of the mean surface $\xb(\cdot)^T \betab + w(\cdot)$ in the model (<ref>).

§.§ Inference for spatial range autocorrelation

Besides the estimation of mean surface $w_0(\cdot)$, it is also important to have valid Bayesian inference for the spatial correlation. We demonstrate that the DISK estimate of the correlation structure is close to the true posterior correlation structure in the simulation example of Section <ref>; see Figure <ref>. In the following, we describe some heuristics; a detailed study is left for future research. To explain why our combination scheme preserves the spatial range autocorrelation, let us consider a special case of model (<ref>) where $\betab=0$, $\sigma^2$ and $\tau^2$ are known and fixed, and only the spatial range parameter $\phi$ is unknown and needs to be estimated. In this case, $\phi$ completely determines the correlation structure; therefore, for each data subset, the log likelihood function for the model in (<ref>) can be written as

\begin{align} \label{eq:prop1}
\log L_j(\phi) = k \log \ell_j(\phi) &= -\frac{km}{2} \log 2\pi - \frac{k}{2} \log \vert \Rb_j \vert - \frac{k}{2} \tr \left( \Rb_j^{-1} \yb_j \yb_j^T \right),
\end{align}

where $\ell_j(\phi)$ is the likelihood function without stochastic approximation, $\Rb_j$ is the shorthand for $\Rb_j(\phi)= \Cb_j(\phi) + \tau^2\Ib$ and $\Cb_j(\phi)$ is $\Cb_j(\alphab)$ with a fixed and known $\sigma^2$. If $\sigma^2$ is fixed at the true value, $\tau^2=0$, and $\phi$ is assigned a prior on the closed interval $[\phi_l, \phi_u]$ ($0<\phi_l<\phi_u<\infty$) that includes the true value $\phi_0$, then the maximum likelihood estimator (MLE) of $\phi$ over $\phi\in [\phi_l, \phi_u]$ is consistent and asymptotically normal for some commonly used covariance functions such as Matérn (33). Motivated by the existing frequentist theory on MLE, we conjecture that a Bernstein-von Mises (BvM) theorem under fixed $\sigma^2$ and $\tau^2$ possibly holds for $\phi$. Let $\Rb'_j(\phi) = \partial \Rb_j (\phi) / \partial \phi$ and $\Rb''_j(\phi) = \partial^2 \Rb_j(\phi) / \partial \phi^2$ be the matrices consisting of component-wise derivatives with respect to $\phi$. Let $\PP_{\phi_0}$ be the probability measure of the subset data $\yb_j$ for all $j=1,\ldots,k$ and also the full data $\yb$.

From (<ref>), we can derive the derivatives of log-likelihood function:

\begin{align}
\label{eq:prop2}
\frac{d \log L_j(\phi)}{d \phi} &= k\frac{d \log \ell_j(\phi)}{d \phi}  = - \frac{k}{2} \tr \left( \Rb_j^{-1} \Rb_j' \right) + \frac{k}{2} \tr(\Rb_j^{-1} \Rb_j' \Rb^{-1}_j \yb_j \yb_j^T ),
\nonumber \\
\frac{d^2 \log L_j(\phi)}{d\phi^2} &= k\frac{d^2 \log \ell_j(\phi)}{d \phi^2} =  - \frac{k}{2} \tr \left( \Rb_j^{-1} \Rb_j'' - \Rb_j^{-1} \Rb_j' \Rb_j^{-1} \Rb_j'   \right)  \nonumber \\
&\qquad \qquad \qquad ~~ - \frac{k}{2} \tr \{(2 \Rb^{-1}_j \Rb_j' \Rb_j^{-1} \Rb_j' \Rb^{-1}_j - \Rb_j^{-1} \Rb_j'' \Rb^{-1}_j ) \yb_j \yb_j^T\}.
\end{align}

At the true parameter $\phi=\phi_0$, the second derivative $d^2 \log \ell_j(\phi)/d \phi^2$ simplifies to $-\frac{1}{2}\tr \left\{ \Rb_j^{-1}(\phi_0) \Rb_j^{'}(\phi_0) \Rb_j^{-1}(\phi_0)\Rb_j^{'}(\phi_0)  \right\}$ because $\EE_{\phi_0}(\yb_j\yb_j^T)=\Rb_j$. In the special case where $\tau^2=0$ and $C_{\alphab}(\cdot,\cdot)$ is the Matérn covariance function with smoothness parameter $\nu$, one can derive from Theorem 2 of 33 that the asymptotic variance of the MLE of $\phi$ is $\phi_0^2/(2\nu^2)$. This implies that for the Matérn covariance function,

$$\lim_{m\to\infty} \frac{1}{2m}\tr \left\{ \Rb_j^{-1}(\phi_0) \Rb_j^{'}(\phi_0) \Rb_j^{-1}(\phi_0)\Rb_j^{'}(\phi_0)  \right\} = 2\nu^2/\phi_0^2.$$

Let $\Pi_{\phi,j}$ ($j=1,\ldots,k$) be the subset posteriors of $\phi$ after stochastic approximation. Let $\overline \Pi_{\phi,n}$ be the DISK posterior of $\phi$, which is the Wasserstein barycenter of $\Pi_{\phi,j}$'s and can be obtained by averaging the quantiles of $\Pi_{\phi,j}$'s over $j=1,\ldots,k$ (similar to $\overline \nu$ in (<ref>)). Let $\Pi_{\phi,n}$ be the full data posterior of $\phi$. Then, for the special case where the covariance function is Matérn and $\tau^2=0$, based on Theorem 2 of 33, we conjecture that under certain regularity conditions, one can follow the techniques used in 37 and show that both $\overline \Pi_{\phi,n}$ and $\Pi_{\phi,n}$ follow the BvM theorem and are asymptotically close to two normal distributions, each with the variance $2\nu^2/\phi_0^2$ but possibly different means. We conjecture that a similar result to 37 holds, in the sense that the means between the two normal limits are of the order $O_p(1/\sqrt{m})$, and hence $m^{1/2}\cdot W_2\left(\overline \Pi_{\phi,n}, \Pi_{\phi,n}\right)\to 0 $ as $m\to\infty$, in the $\PP_{\phi_0}$-probability if the data $\yb$ has the probability measure $\PP_{\phi_0}$. Furthermore, based on Theorem 3 of 33, we conjecture that with high $\PP_{\phi_0}$-probability, for any $\phi$ drawn from either the true posterior $\Pi_{\phi,n}$ or the DISK posterior $\overline \Pi_{\phi,n}$, the Gaussian process predictive variance at any testing location $\sbb^*\in \Scal$ is asymptotically close to the true Gaussian process predictive variance with the true range parameter $\phi_0$, in the sense that the ratio of these two predictive variances tends to 1 as $m\to\infty$. If $\sigma^2$ is unknown, it is well known that $\sigma^2$ and $\phi$ cannot be identified [70] but a microergodic parameter like $\sigma^2\phi^{2\nu}$ for the Matérn covariance functions can be identified. In this case, we conjecture that a similar BvM theorem can be shown for the microergodic parameter and similar conclusions about the equivalence of prediction variances can be obtained. We leave the thorough investigation on these issues for general covariance functions to future research.

*  The covariance function is continuous and at least three times differentiable in $\phi\in [\phi_l,\phi_u]$. The true parameter $\phi_0$ is an interior point of $(\phi_l,\phi_u)$.

*  $\underset{m \rightarrow \infty} {\lim} (2m)^{-1} \tr \left\{ \Rb_j^{-1}(\phi_0) \Rb_j^{'}(\phi_0) \Rb_j^{-1}(\phi_0)\Rb_j^{'}(\phi_0)  \right\} = \Ical_{\phi_0} $, where $0<\Ical_{\phi_0}< \infty$ is a constant that does not depend on $m$ and $k$.

*  There exists an envelope function $M(\yb_j)$, such that $m^{-1}|d^r \log \ell_j(\phi)/d \phi^r|\leq M(\yb_j)$ for $r=0,1,2,3$ and all $\phi\in [\phi_0-\delta_0,\phi_0+\delta_0]$ where $\delta_0>0$ is a small constant, $\EE_{\phi_0}[M(\yb_j)]<C_H$ for some constant $0<C_H<\infty$, for all $j=1,\ldots,k$.

*  For all $j=1,\ldots,k$, for any $\delta>0$, there exists an $\epsilon>0$ such that as $m\to\infty$, $\sup_{|\phi-\phi_0|\geq \delta} m^{-1} \left\{ \log \ell_j(\phi)-\log \ell_j(\phi_0)\right\} < -\epsilon$ holds with $\PP_{\phi_0}$-probability approaching 1.

*  The prior density function $\pi(\phi)$ is positive and continuous on $[\phi_l,\phi_u]$.

We have the following result about the Wasserstein-2 distance between the DISK posterior and the full data posterior.

Suppose that Model (<ref>) holds with $\betab=0$, $\sigma^2$ and $\tau^2$ being known and fixed. Let $\phi_0$ be the true value of the spatial range parameter. Suppose that Assumptions <ref>-<ref> hold. Then $m^{1/2}\cdot W_2\left(\overline \Pi_{\phi,n}, \Pi_{\phi,n}\right)\to 0 $ as $m\to\infty$ in $\PP_{\phi_0}$-probability.

Proposition <ref> implies that the DISK posterior of $\phi$ is close to the full data posterior of $\phi$ as $m\to\infty$, and the difference has the order $o_p(m^{-1/2})$. This establishes the validity of using the DISK posterior to make inference on the spatial range autocorrelation.

Suppose that Model (<ref>) holds with $\betab=0$, $\sigma^2$ and $\tau^2$ being known and fixed. Let $\phi_0$ be the true value of the spatial range parameter. Suppose that Assumptions <ref>-<ref> hold. Then for the subset data $\yb_j$, there exists a root $\hat\phi_j$ of the equation $d \log \ell_j(\hat\phi_j)/d \phi=0$ on $[\phi_l,\phi_u]$ that satisfies $\hat\phi_j \to \phi_0$ with $\PP_{\phi_0}$-probability approaching 1, as $m\to\infty$. Similarly, for the full data $\yb$, there exists a root $\hat\phi$ on $[\phi_l,\phi_u]$ of the equation $d \log \ell(\hat\phi)/d \phi=0$ that satisfies $\hat\phi \to \phi_0$ with $\PP_{\phi_0}$-probability approaching 1, as $m\to\infty$.

Furthermore, let $h = \sqrt{km}(\phi - \hat \phi_j)$. Then it follows that $\int_{-\infty}^{+\infty}(1+h^2)\Big|\pi_h\left(h\mid \yb_j \right) - \sqrt{\frac{\Ical_{\phi_0}}{2\pi}}e^{-\frac{\Ical_{\phi_0}}{2}h^2}\Big| dh \to 0$ as $m\to\infty$, for all $j=1,\ldots,k$, where $\pi_h\left(h\mid \yb_j \right) = (km)^{-1/2} \pi_{\phi}\left(\hat\phi_j + h/\sqrt{km} \mid \yb_j \right)$. Similarly, for the full data $\yb$, let $t = \sqrt{km}(\phi - \hat \phi)$. Then it follows that $\int_{-\infty}^{+\infty}(1+t^2)\Big|\pi_t\left(t\mid \yb \right) - \sqrt{\frac{\Ical_{\phi_0}}{2\pi}}e^{-\frac{\Ical_{\phi_0}}{2}t^2}\Big| dt \to 0$ as $m\to\infty$, where $\pi_t\left(t\mid \yb \right) = (km)^{-1/2} \pi_{\phi}\left(\hat\phi + t/\sqrt{km} \mid \yb \right)$.

The proof generally follows the proofs of Theorem 4.2 in <cit.> and Lemma 2 of 37. Even though the data are not i.i.d, our Assumptions C.3 and C.4 replace Assumptions A2 and A4 in <cit.>, and Assumptions 3 and 5 in 37. In the following, we only prove the results for each subset data $\yb_j$. The results for the full data $\yb$ can be shown similarly and are therefore omitted.

First we show that the sequence of $\hat\phi_j$ described in the lemma exists. We can derive the derivatives of log-likelihood function:

\begin{align}
\label{loglikderiv2}
\frac{d \log L_j(\phi)}{d \phi} &= k\frac{d \log \ell_j(\phi)}{d \phi}  = - \frac{k}{2} \tr \left( \Rb_j^{-1} \Rb_j' \right) + \frac{k}{2} \tr(\Rb_j^{-1} \Rb_j' \Rb^{-1}_j \yb_j \yb_j^T ),
\nonumber \\
\frac{d^2 \log L_j(\phi)}{d\phi^2} &= k\frac{d^2 \log \ell_j(\phi)}{d \phi^2} =  - \frac{k}{2} \tr \left( \Rb_j^{-1} \Rb_j'' - \Rb_j^{-1} \Rb_j' \Rb_j^{-1} \Rb_j'   \right)  \nonumber \\
&\qquad \qquad \qquad ~~ - \frac{k}{2} \tr \{(2 \Rb^{-1}_j \Rb_j' \Rb_j^{-1} \Rb_j' \Rb^{-1}_j - \Rb_j^{-1} \Rb_j'' \Rb^{-1}_j ) \yb_j \yb_j^T\}.
\end{align}

At $\phi=\phi_0$, we have

\begin{align*}
&\EE_{\phi_0}\left[\frac{d\log \ell_j(\phi)}{d\phi}\right] = \frac{1}{2}\tr\left\{\Rb_j^{-1}(\phi_0) \Rb_j'(\phi_0) \Rb^{-1}_j (\phi_0) \EE_{\phi_0}[\yb_j \yb_j^T] \right\} - \frac{1}{2} \tr \left\{ \Rb_j^{-1}(\phi_0) \Rb_j'(\phi_0) \right\} \\
&=  \frac{1}{2}\tr\left\{\Rb_j^{-1}(\phi_0) \Rb_j'(\phi_0) \Rb^{-1}_j(\phi_0) \Rb_j (\phi_0)\right\} - \frac{1}{2} \tr \left\{ \Rb_j^{-1}(\phi_0) \Rb_j'(\phi_0) \right\} =0.\\
&\EE_{\phi_0}\left[\frac{d^2\log \ell_j(\phi)}{d\phi^2}\right] = - \frac{1}{2} \tr \left\{\Rb_j^{-1} (\phi_0)\Rb_j'' (\phi_0)- \Rb_j^{-1} (\phi_0)\Rb_j'(\phi_0) \Rb_j^{-1}(\phi_0) \Rb_j' (\phi_0)  \right\} \nonumber \\
&- \frac{1}{2} \tr \left[\left\{2 \Rb^{-1}_j(\phi_0) \Rb_j'(\phi_0) \Rb_j^{-1}(\phi_0) \Rb_j'(\phi_0) \Rb^{-1}_j(\phi_0) - \Rb_j^{-1} (\phi_0)\Rb_j'' (\phi_0)\Rb^{-1}_j(\phi_0) \right\} \EE_{\phi_0}(\yb_j \yb_j^T)\right]\nonumber \\
&= -\frac{1}{2}\tr \left\{2 \Rb^{-1}_j(\phi_0) \Rb_j'(\phi_0) \Rb_j^{-1}(\phi_0) \Rb_j'(\phi_0) \Rb^{-1}_j(\phi_0)\right\}.
\end{align*}

By Assumption C.3, $d \log \ell_j(\phi)/d\phi$ and $d^2 \log \ell_j(\phi)/d\phi^2$ are both continuously differentiable in $[\phi_0-\delta_0,\phi_0+\delta_0]$. This implies that as $m\to\infty$, with large $\PP_{\phi_0}$-probability, there exists a root for the equation
$d \log \ell_j(\phi)/d\phi=0$ inside the neighborhood $[\phi_0-\delta,\phi_0+\delta]$ for any given $\delta>0$ that attains the maximum of $\log \ell(\theta)$. Denote the root by $\hat\phi_j$. By Assumption C.4, we have that for the given $\delta>0$, if $\phi$ is outside $[\phi_0-\delta,\phi_0+\delta]$, then there exists a constant $\epsilon>0$ such that for all sufficiently large $m$, with large $\PP_{\phi_0}$-probability, $\log \ell_j(\phi)< \log \ell_j(\phi_0)-m\epsilon \ll \log \ell_j(\phi_0) \leq \log \ell(\hat\phi_j) $, and hence the MLE on $[\phi_l,\phi_u]$ must be attained inside $[\phi_0-\delta,\phi_0+\delta]$ with large $\PP_{\phi_0}$-probability. This has proved that $\hat\phi_j\to \phi_0$ in $\PP_{\phi_0}$-probability.

Next, we prove the Bernstein-von Mises theorem for the subset posterior $\Pi_{\phi,j}$, where $j=1,\ldots,k$. Let

\begin{align*}
v(h) &= \log \ell_j\left\{\hat \phi_j +\frac{h}{(Km)^{1/2}}\right\} - \log \ell_j(\hat \phi_j), \\
C_v & = \int e^{k\cdot v(h)} \pi\left\{\hat \phi_j +\frac{h}{(k m)^{1/2}}\right\} dh.
\end{align*}

The $j$th subset likelihood function (after stochastic approximation) can be written as

$$\pi_{\phi} \left(\phi \mid\yb_j\right) = \frac{\exp\left\{k\log  \ell_j(\phi)\right\} \pi(\phi) d \phi} {\int \exp\left\{k\log \ell_j(\theta)\right\} \pi(\phi) d \phi}.$$

The induced $j$th subset posterior density of $h=(km)^{1/2}(\phi-\hat \phi_j)$ is

\begin{align*}
\pi_{h}(h\mid\yb_j) & =  \frac{1}{C_v}\exp\{kv(h)\}\pi\left\{\hat \phi_j + \frac{h}{(k m)^{1/2}}\right\}.
\end{align*}

Let $\mathcal{H} = \{h=(k m)^{1/2}(\phi-\hat\phi_j):\phi\in[\phi_l,\phi_u]\}$. Define

\begin{align*}
g(h) &= (1+h^2) \left[ e^{kv(h)} \pi\left\{\hat \phi_j +\frac{h}{(k m)^{1/2}}\right\} - \exp\left\{-\frac{1}{2}h^2\Ical_{\phi_0}\right\}\pi(\phi_0) \right].
\end{align*}

If we can show that $\int_{\mathcal{H}} |g(h)|d h \to 0$ as $m\to \infty$ in $\PP_{\phi_0}$-probability, then by Assumption C.5,

$$ C_v \to \int_{\mathcal{H} } \exp\left\{-\frac{1}{2}h^2\Ical_{\phi_0}\right\}\pi(\phi_0) d h  = (2\pi/\Ical_{\phi_0})^{1/2}\pi(\phi_0) >0,$$

as $m\to \infty$ in $\PP_{\phi_0}$-probability, and

\begin{align}\label{tv2exp}
& \int_{\mathcal{H} } (1+h^2)\left| \pi_h(h\mid \yb_j) - \frac{1}{(2\pi/\Ical_{\phi_0})^{1/2} }\exp\left\{-\frac{1}{2}h^2\Ical_{\phi_0}\right\}\right| dh \nonumber \\
={}& \int_{\mathcal{H} }(1+h^2) \left| \frac{ e^{kv(h)} \pi\left\{\hat \phi_j +\frac{h}{(Kk m)^{1/2}}\right\}}{C_v} - \frac{1}{(2\pi/\Ical_{\phi_0})^{1/2} }\exp\left\{-\frac{1}{2}h^2\Ical_{\phi_0}\right\}\right| dh \nonumber \\
\leq{}& \frac{1}{C_v} \int_{\mathcal{H}} |g(h)| dh +  \left|\frac{(2\pi/\Ical_{\phi_0})^{1/2} \pi(\phi_0)}{C_v}-1\right|\times \nonumber\\
& \int_{\mathcal{H} }  \frac{1+h^2}{(2\pi/\Ical_{\phi_0})^{1/2}}\exp\left\{-\frac{1}{2}h^2\Ical_{\phi_0}\right\}dh \to 0,
\end{align}

as $m\to\infty$ in $\PP_{\phi_0}$-probability, where the last step follows since the integral $\int_{-\infty}^{+\infty} \frac{1+h^2}{(2\pi/\Ical_{\phi_0})^{1/2}}\exp\left\{-\frac{1}{2}h^2\Ical_{\phi_0}\right\}dh$ is a finite number that only depends on $\Ical_{\phi_0}$. This is the exact conclusion of Lemma <ref> we need to show. Therefore, it suffices to prove $\int_{\mathcal{H}} |g(h)|d h \to 0$ as $m\to \infty$ in $\PP_{\phi_0}$-probability, under Assumptions C.1-C.5.

We break the integral of $g(h)$ into three parts: $A_1=\{h: |h|\geq \delta_1 (km)^{1/2}\}$, $A_2=\{h:\delta_2\leq |h|< \delta_1(k m)^{1/2}\}$, $A_3=\{h:|h|<\delta_2\}$, where the constants $\delta_1,\delta_2$ will be chosen later. Then

\begin{align}\label{gh}
\int_{\mathcal{H}} |g(h)| dh \leq \int_{A_1} |g(h)| dh + \int_{A_2} |g(h)| dh +  \int_{A_3} |g(h)| dh.
\end{align}

Additionally, assume that for every $j=1, \ldots, k$,

\begin{align}\label{eq:cond-schef}
\int  \frac{L_{j}(\hat \phi_{jm} + h / \sqrt{km})}{L_{j}(\hat \phi_{jm})} \pi(\hat \phi_{jm} + h / \sqrt{km}) d h \rightarrow
\int e^{-\frac{1}{2}h^2 i_{\phi_0}} \pi(\phi_0) dh
\end{align}

with probability 1.

We prove that under very general regularity conditions, $\overline \Pi_{\phi, n}$ and $\Pi_{\phi, n}$ are close to each other.

The model in (<ref>) implies that $\yb_j \sim N(\zero, \Rb_j)$, where $\Rb_j$ is a shorthand for $\Rb_j(\phi)$, and the stochastic-approximated log likelihood on subset $j$ is

\begin{align}
\label{eq:prop1}
\log L_j(\phi) = k \log \ell_j(\phi) &= -\frac{km}{2} \log 2\pi - \frac{k}{2} \log \vert \Rb_j \vert - \frac{k}{2} \tr \left\{ \Rb_j^{-1} \yb_j \yb_j^T \right\}.
\end{align}

Let $\Rb'_j$ and $\Rb''_j$ be the shorthand for $\Rb'_j(\phi)$ and $\Rb''_j(\phi)$. Differentiating the modified log likelihood in (<ref>) one and two times with respect to $\phi$ gives

\begin{align}
\label{eq:prop2}
\frac{\partial \log L_j(\phi)}{\partial \phi} &= \psi(\yb_j, \phi)   = - \frac{k}{2} \tr \left( \Rb_j^{-1} \Rb_j' \right) + \frac{k}{2} \tr(\Rb_j^{-1} \Rb_j' \Rb^{-1}_j \yb_j \yb_j^T ),
\nonumber \\
\frac{\partial^2 \log L_j(\phi)}{\partial \phi^2} &= \frac{\partial \psi(\yb_j, \phi)}{\partial \phi} =  - \frac{k}{2} \tr \left( \Rb_j^{-1} \Rb_j'' - \Rb_j^{-1} \Rb_j' \Rb_j^{-1} \Rb_j'   \right)  \nonumber \\
&\qquad \qquad \qquad ~~ - \frac{k}{2} \tr \{(2 \Rb^{-1}_j \Rb_j' \Rb_j^{-1} \Rb_j' \Rb^{-1}_j - \Rb_j^{-1} \Rb_j'' \Rb^{-1}_j ) \yb_j \yb_j^T\} .
\end{align}

If $\Ab_j = \Rb_j^{-1/2} \Rb_j' \Rb_j^{-1/2} $, where $\Rb_j^{1/2}$ is the symmetric square root of $\Rb_j$ and $\Rb_j^{-1/2}$ is the inverse of $\Rb_j^{1/2}$, and $\Bb_j = 2 \Ab_j^2 - \Rb_j^{-1/2} \Rb_j'' \Rb^{-1/2}_j$, then   (<ref>) reduces to

\begin{align}
\label{eq:prop3}
\frac{\partial \psi(\yb_j, \phi)}{\partial \phi} &=  \frac{k}{2} \tr \left\{  \Rb_j^{-1/2} \Ab_j^2 \Rb_j^{1/2}  -  \Rb_j^{-1} \Rb_j'' -
\Rb^{-1/2}_j  \Bb_j \Rb^{-1/2}_j \yb_j \yb_j^T \right\}  \nonumber\\
&=  \frac{k}{2} \tr ( \Ab_j^2)  - \frac{k}{2}  \tr (\Rb_j^{-1} \Rb_j'')
-  \frac{k}{2} (\Rb_j^{-1/2} \yb_j)^T \Bb_j (\Rb_j^{-1/2} \yb_j),\nonumber\\
- \frac{\partial \psi(\yb_j, \phi)}{\partial \phi} &=
- \frac{k}{2} \tr ( \Ab_j^2)  + \frac{k}{2}  \tr (\Rb_j^{-1} \Rb_j'')
+  \frac{k}{2} (\Rb_j^{-1/2} \yb_j)^T \Bb_j (\Rb_j^{-1/2} \yb_j), \nonumber\\
\left | \frac{\partial \psi(\yb_j, \phi)}{\partial \phi} \right |
&\leq \frac{k}{2}  \tr (\Rb_j^{-1} \Rb_j'')
+  \frac{k}{2} (\Rb_j^{-1/2} \yb_j)^T \Bb_j (\Rb_j^{-1/2} \yb_j) \equiv g(\yb_j).
\end{align}

The expectation of $g(\yb_j)$ in (<ref>) with respect to the true distribution of $\yb_j$  equals

\begin{align}
\label{eq:prop4}
\EE_{\phi_0} \left\{ g(\yb_j) \right\}
=  \frac{k}{2}  \tr (\Rb_j^{-1} \Rb_j'')
+  \frac{k}{2} \tr \left\{ \Rb_j^{-1/2} \Bb_j \Rb_j^{-1/2} \Rb_j(\phi_0)  \right\} < \infty,
\end{align}

which is finite based on our assumptions 1, 2, and 3 for any given $k,m$, and $n$. The bound on $\partial \psi(\yb_j, \phi) / \partial \phi$ in (<ref>) is also uniform for every value of $\phi \in [\phi_u, \phi_l]$; therefore, $-\EE_{\phi_0} \left\{  \partial \psi(\yb_j, \phi) / \partial \phi \right\}$ at $\phi=\phi_0$ exists using dominated convergence theorem and equals

\begin{align}
\label{eq:prop5}
- \EE_{\phi_0}  \frac{\partial \psi(\yb_j, \phi)}{\partial \phi} \bigg|_{\phi = \phi_0} = \frac{km}{2} \frac{1}{m} \tr \{\Rb^{-1}_j(\phi_0) \Rb_j'(\phi_0) \Rb_j^{-1}(\phi_0) \Rb_j'(\phi_0)\}.
\end{align}

Our assumption 4 implies that

\begin{align}
\label{eq:prop6}
\underset{m \rightarrow \infty} {\lim}  \frac{1}{km} \EE_{\phi_0}  \frac{\partial \psi(\yb_j, \phi)}{\partial \phi} \bigg|_{\phi = \phi_0} =  -i_{\phi_0}.
\end{align}

We now prove the density of the $j$th subset posterior $\phi$ is close (in total variation) to a Gaussian distribution centered at $\hat \phi_{jm}$ with covariance matrix $\{(km) i_{\phi_0}\}^{-1}$. Expand $\tilde \ell_j(\phi) = \log L_j(\phi)$ in the neighborhood $\hat \phi_{jm}$ to obtain that

\begin{align}
\label{eq:bv1}
\tilde \ell_{j}(\phi) &=
\tilde \ell_j(\hat \phi_{jm}) +
(\phi - \hat \phi_{jm}) \psi(\yb_j, \hat \phi_{jm}) + (\phi - \hat \phi_{jm})^2b_m(\phi),%  \nonumber \\
% &=
%   \tilde \ell_j(\hat \phi_m) + (\phi - \hat \phi_m)^2
%   \int_{0}^1 \int_{0}^1 v \partial \psi(\yb_j, \hat \phi_m + uv(\phi - \hat \phi_m)) / \partial \phi \nonumber \\
% &= \tilde \ell_j(\hat \phi_m) - (\phi - \hat \phi_m)^2
%   b_m, \quad b_m = - \int_{0}^1 \int_{0}^1 v \partial \psi(\yb_j, \hat \phi_m + uv(\phi - \hat \phi_m)) / \partial \phi
\end{align}

where $b_m(\phi) = \int_{0}^1 \int_{0}^1 v [\partial \psi\{\yb_j, \hat \phi_{jm} + uv(\phi - \hat \phi_{jm})\} / \partial \phi] du dv$.
Because $\hat \phi_{jm}$ is a strongly consistent sequence  of roots of the likelihood equation,  $\psi(\yb_j, \hat \phi_{jm}) = 0$ for a sufficiently large $m$ with probability 1. Define the local parameter $h$ as $h = \sqrt{km}(\phi - \hat \phi_{jm})$ (using $\sqrt{m}$ instead of $\sqrt{km}$ does not seem to work!) and substitute it in (<ref>) to obtain that

\begin{align}
\label{eq:bv2}
e^{\{\tilde \ell_{j}(\hat \phi_{jm} + h / \sqrt{km}) -
\tilde \ell_j(\hat \phi_{jm})\}} =
\frac{L_{j}(\hat \phi_{jm} + h / \sqrt{km})}{L_{j}(\hat \phi_{jm})} &= e^{  h^2 \frac{1}{km} b_m(\phi)} .
%                          (\phi - \hat \phi_{jm}) \psi(\yb_j, \hat \phi_{jm}) + (\phi - \hat \phi_{jm})^2b_m,%  \nonumber  &= \exp \left( - h^2 \frac{b_m}{m} \right)
\end{align}

If we show that ${b_m}/{(km)}$ converges to $-i_{\phi_0} / 2$ with probability 1, then

\begin{align}
\label{eq:bv3}
\frac{L_{j}(\hat \phi_{jm} + h / \sqrt{km})}{L_{j}(\hat \phi_{jm})} \pi(\hat \phi_{jm} + h / \sqrt{km})
\rightarrow e^{-\frac{1}{2} h^2 i(\phi_0)} \pi(\phi_0)
\end{align}

almost surely because $\hat \phi_{jm}$  is a strongly consistent sequence of roots of the likelihood equation, so $\hat \phi_{jm} \rightarrow \phi_0$ almost surely and $\pi(\phi)$ is continuous for every $\phi \in [\phi_l, \phi_u]$, including $\phi_0$. If we assume that (<ref>) holds, then an application of Scheffe's Theorem (<cit.>, Chapter 2, Theorem 2) implies that

\begin{align}
\label{eq:bv4}
\int | \pi_{km}(h \mid \yb_j) - N\{h \mid 0, i_{\phi_0}^{-1} \}| \, dh \rightarrow 0,
\end{align}

almost surely, where $\pi_{km}(h \mid \yb_j)$ is the density of $h$ obtained from the $j$th subset posterior density and $N\{h \mid 0, i_{\phi_0}^{-1} \}$ is the density of the Gaussian distribution with mean 0 and variance $i_{\phi_0}^{-1}$.

We finish the proof by showing that ${b_m}/{(km)}$ converges to $i_{\phi_0} / 2$ with probability 1. First, the inequality in (<ref>) implies that $(km)^{-1}\partial \psi(\yb_j, \phi) / \partial \phi$ converges uniformly to

\begin{align}
\label{eq:bv5}
\EE_{\phi_0} \frac{1}{km}\frac{\partial \psi(\yb_j, \phi)}{\partial \phi} =
\frac{1}{2m} \tr \{ \Ab_j^2 - \Rb_j^{-1} \Rb_j'' -  \Rb_j^{-1/2} \Bb_j \Rb_j^{-1/2} \Rb_j(\phi_0) \}
\end{align}

for every $\phi \in [\phi_l, \phi_u]$ as $m \rightarrow \infty$. We have assume that the covariance function used for defining $\Rb$ is continuous in $\phi$, so $\EE_{\phi_0} \{(km)^{-1}\partial \psi(\yb_j, \phi) / \partial \phi\}$ is continuous in $\phi$ because $\tr$ is a continuous function on positive semi-definite matrices. This implies that for every $\phi \in [\rho_l, \rho_u]$ and any $\epsilon > 0$, there exists a $\delta > 0$ such that for sufficiently large $m$

\begin{align}
\label{eq:bv6}
\vert \phi - \phi_0 \vert  < \delta \implies
\vert \EE_{\phi_0} \{(km)^{-1}\partial \psi(\yb_j, \phi) / \partial \phi\} + i_{\phi_0} \vert < \epsilon,
\end{align}

where we have used that $\EE_{\phi_0} \{(km)^{-1}\partial \psi(\yb_j, \phi) / \partial \phi\}$ equals $-i_{\phi_0}$ at $\phi = \phi_0$ for a sufficiently large $m$ from (<ref>).
Second, Theorem 16(a) in <cit.> can be extended to our case by choosing the parameter space as $[\rho_l, \rho_u]$, defining

\begin{align*}
\frac{1}{km} U(\yb_j, \phi) = \frac{1}{km} \left\{ \log \tilde \ell_j(\phi) - \log \tilde \ell_j(\phi_0) \right\},
\end{align*}

and using $g(\yb_j)$ in (<ref>) as the integrable function that upper bounds $U(\yb_j, \phi)$ in a neighborhood of $\phi_0$. This extension implies that there exists a $N_{\epsilon}$ such that

\begin{align}
\label{eq:bv7}
m > N_{\epsilon} \implies \underset{\phi \in [\phi_l, \phi_u]}{\sup} \big|
(km)^{-1}\partial \psi(\yb_j, \phi) / \partial \phi  - \EE_{\phi_0} \{(km)^{-1}\partial \psi(\yb_j, \phi) / \partial \phi\}
\big| < \epsilon
\end{align}

with probability 1. Assume that $N_{\epsilon}$ is large enough such that (<ref>) and (<ref>) hold simultaneously so that for every $m > N_{\epsilon}$ with probability 1

\begin{align}
\label{eq:bv8}
|(km)^{-1} b_m + i_{\phi_0}/2 | &=
\big| \int_{0}^1  \int_{0}^1 v \frac{\partial \psi(\yb_j, \hat \phi_m + uv \frac{h}{km})}{km \partial \phi}
du  dv + \int_{0}^1 v i_{\phi_0} dv  \big| \nonumber \\
\int_{0}^1  \int_{0}^1 v  \big| \frac{\partial \psi (\yb_j, \hat \phi_m + uv \frac{h}{km})}{km \partial \phi}
+ i_{\phi_0}  \big| du  dv \nonumber\\
\int_{0}^1  \int_{0}^1 v \underset{\phi \in [\phi_l, \phi_u]}{\sup} \bigg| \frac{\partial \psi(\yb_j, \phi)}{km \partial \phi} - \EE_{\phi_0}\frac{\partial \psi (\yb_j, \phi)}{km \partial \phi} \bigg| dudv + \nonumber\\
& ~~ \int_{0}^1  \int_{0}^1 v \bigg| \EE_{\phi_0}\frac{\partial \psi (\yb_j, \phi)}{km \partial \phi} + i_{\phi_0} \bigg| du dv
\leq \epsilon / 2 + \epsilon / 2 = \epsilon;
\end{align}

therefore, ${b_m}/{(km)}$ converges to $-i_{\phi_0} / 2$ with probability 1 and the proposition is proved.

We use this lemma to prove the following proposition which says that under general regularity assumptions, $\Pi_{\phi, n}$ and $\overline \Pi_{\phi, n}$ agree on inference on $\phi$ as $n \rightarrow \infty$.

$W_2^2(\Pi_{\phi, n}, \overline \Pi_{\phi, n}) = $

§ EXPERIMENTS

§.§ Simulation setup

We compare DISK with its competitors on synthetic data based on its performance in learning the process parameters, interpolating the unobserved  spatial surface, and predicting at new locations. This section presents three simulation studies. The first (Simulation 1) and second  (Simulation 2) simulations represent moderately large dataset with $12,025$ locations and the third (Simulation 3) simulation  analyzes a large dataset with $1,002,025$ locations. In first two simulations, we randomly select the data at $n=10^4$ locations for model fitting and $l=2025$ locations for predictions, while in Simulation 3 the training and test data are of size $n= 10^6$ and $l = 2025$.

For all three simulations, the response is simulated at $(n+l)$ locations as

\begin{align}
\label{eq:sim1}
y(\sbb_i) = \beta_0 + w_0(\sbb_i) + \epsilon_i, \quad \epsilon_i \sim N(0, \tau_0^2), \quad i = 1, \ldots, n+l.
\end{align}

Simulations 1 and 3 follow the data generation scheme described in [23]. Specifically, we set $\Dcal = [-2, 2] \times [-2, 2] \subset \RR^2$ and uniformly sample $(n+l)$ spatial locations $\sbb_i = (s_{i1}, s_{i2})$ in $\Dcal$ ($i = 1, \ldots, n+l$). For any $s \in [-2, 2]$,  define the function $f_0(s) = e^{-(s-1)^2} + e^{-0.8(s+1)^2} - 0.05 \sin \{8(s + 0.1)\}$ and set $w_0(\sbb_i)= - f_0(s_{i1})  f_0(s_{i2})$. Although the function $w_0(\cdot)$ simulated in this way is theoretically infinitely smooth,
the response surface simulated from (<ref>) exhibits complex local behavior, which is challenging to capture using spatial process-based models as we demonstrate later. Simulation 2 generates $\{w_0(\sbb_1),\ldots,w_0(\sbb_{n+l})\}$ from a GP(0,$\sigma_0^2\exp\{-\phi_0 \|\sbb_i-\sbb_j\|)\}$. This is a more familiar simulation scenario in the spatial context with the generated spatial surface being continuous everywhere but differentiable nowhere. As argued earlier, simulating data from an ordinary Gaussian process is expensive, hence we refrain from a massive simulation study under the scenario of Simulation 2. For  Simulations 1 and 3, the intercept $\beta_0$ and true error variance $\tau_0^2$ are set to $1$ and $0.01$, respectively. For Simulation 2, $\beta_0$, $\tau_0^2$, $\phi_0$ and $\sigma_0^2$ are taken to be $1$, $0.1$, $9$ and $1$ to keep long range spatial dependence and high spatial variance to nugget ratio. The inferential and predictive results for all simulations are based on $10$ replications. We compare DISK with a number of Bayesian and non-Bayesian spatial models in both simulations:

(i) LatticeKrig [43] using the  package in R [46]  with 3 resolutions [44];
(ii) nearest neighbor Gaussian process (NNGP) using the  package in R with the number of nearest neighbors (NN) as 5, 15, and 25 [14];
(iii) full-rank Gaussian process (GP) using the  package in R [17] with the full data;
(iv) modified predictive process (MPP) using the  package in R with the full data; and
(v) locally approximated Gaussian process (laGP) using the  package in R [23].

All five methods produce results in Simulation 1 and 2, but the first four methods fail due to numerical issues in Simulation 3. While laGP is not designed for full scale Bayesian inference, it is used as the benchmark for predictive point estimates and associated standard errors in (<ref>) due to its popularity in fitting computer models. For the spatial process-based Bayesian models, we employ (<ref>) with only an intercept $\beta$, putting a $N(0, 100)$ prior on $\beta$, a GP prior on $w(\cdot)$, and $\text{IG}(2, 0.1)$ prior on $\tau^2$, where IG$(a, b)$  is the Inverse-Gamma distribution with mean $b/(a-1)$ and variance $b / \{(a-1)^2 (a-2)\}$ for $a > 2$. In model fitting, we assume an exponential correlation in the random field given by $\cov\{w(\sbb),w(\sbb')\} = \sigma^2 e^{-\phi \| \sbb-\sbb' \|}, \; \sbb, \sbb' \in \Dcal$ and put IG(2, 2) prior on $\sigma^2$ and a uniform prior on $\phi$. The MPP prior on $w(\cdot)$ is fitted by setting the rank $r$ as 200 and 400, respectively, where the $r$ knots are selected randomly from the domain $\Dcal$. NNGP is chosen as a representative example of the current state-of-the-art Bayesian method for inference and predictions in massive spatial data.

The three-step DISK framework is applied using the full-rank GP and the low-rank MPP priors using the algorithm outlined in Section <ref>, yielding DISK (GP) and DISK (MPP) procedures, respectively. For all our simulations, DISK combines the subset marginal posteriors by averaging their quantiles, as described in Section <ref>, and we set $\xi = 10^{-4}$ in Equation (<ref>). We use consensus Monte Carlo (CMC; [51]), semiparametric density product (SDP; [42]) and meta kriging (MK; [25]) as representative  competitors for model-free subset posterior aggregation to highlight the advantages of DISK. Similar to DISK, these three approaches also operate in three steps. In steps 1 and 2, the MPP-based model in (<ref>) is fitted on every subset for CMC, SDP and MK using the  package. Unlike DISK, the methods of CMC, SDP and MK do not employ stochastic approximation in the sampling step. Third, we use  package with the default setting [41] for combining subset posterior MCMC samples in CMC and SDP, yielding CMC(MPP) and SDP (MPP) procedures respectively and  package for combining subset inferences in MK to yield MK (MPP). SDP (MPP) fails due to numerical issues when the posterior MCMC samples for predictions and the surface are combined. Further, Simulation 1 shows less than competitive performance for MK (MPP) in estimating the spatial surface and posterior distribution of the parameters, hence Simulations 2 and 3 do not use it as a competitor. Identical priors, covariance functions, ranks, and knots are used for the non-distributed process models and their distributed counterparts for a fair comparison.

All experiments are run on an Oracle Grid Engine cluster with 2.6GHz 16 core compute nodes. The non-distributed methods (LatticeKrig, GP, MPP, NNGP, and laGP) and the distributed methods (CMC, DISK, MK, and SDP) are allotted memory resources of 64GB and 16GB, respectively. Every MCMC sampling algorithm runs for 15,000 iterations, out of which the first 10,000 MCMC samples are discarded as burn-in MCMC samples and the rest of the chain is thinned by collecting every fifth MCMC sample. Convergence of the chains to their stationary distributions is confirmed using trace plots. All the interpolated spatial surfaces are obtained using the  package in .

We compare the quality of prediction and estimation of  spatial surface at predictive locations $\Scal^*=\{\sbb_1^*,\ldots,\sbb_l^*\}$. If $w_0(\sbb_{i'}^*)$ is the value of the spatial surface at $\sbb_{i'}^* \in \Scal^*$, the estimates of bias, variance, and Bayes $L_2$-risk in estimating $w_0(\cdot)$ are defined as

\begin{align}
\label{eq:sim2}
{\text{bias}}^2 = \frac{1}{l} \sum_{i'=1}^l \{\hat w(\sbb_{i'}^*) - w_0(\sbb_{i'}^*)\}^2, ~~ {\text{var}} = \frac{1}{l} \sum_{i'=1}^l \hat {\text{var}} \{w(\sbb_{i'}^*)\}, ~~ {L_2\text{-risk}} = {\text{bias}}^2 +  {\text{var}},
\end{align}

where $\hat w(\sbb_{i'}^*)$ and $\hat {\text{var}} \{w(\sbb_{i'}^*)\}$ denote the estimate of $w_0(\sbb_{i'}^*)$  obtained using any distributed or non-distributed methods and its variance, respectively. For sampling-based methods, we set $\hat w(\sbb_{i'}^*)$ and $\hat {\text{var}} \{w(\sbb_{i'}^*)\}$ to be the median and the variance of the posterior MCMC samples for $ w(\sbb_{i'}^*)$, respectively, for $i'=1, \ldots, l$. We also estimate the point-wise 95% credible or confidence intervals (CIs) of $w(\sbb^*_{i'})$  and predictive intervals (PIs) of $y(\sbb^*_{i'})$ for every $\sbb_{i'} \in \Scal^*$ and compare the CI and PI coverages and lengths for every method.  Also, the point predictive performance at the locations in  $\Scal^*$ are compared across competitors using the mean square prediction error (MSPE) defined as MSPE$= \sum_{i'=1}^l \{\hat y(\sbb_{i'}^*) - y(\sbb_{i'}^*)\}^2 / l$. Finally, we compare the performance of all the methods for parameter estimation using the posterior medians or point estimates and the 95% CIs for $\beta$, $\sigma^2$, $\tau^2$, and $\phi$.

§.§ Simulated data analysis

§.§.§ Simulation 1: moderately large spatial data.

We fit DISK (GP) for $k = 10, 20, 30, 40, 50$; CMC (MPP), SDP (MPP), MK (MPP) and DISK (MPP) for $k=10, 20$, along with other competitors. MPP using DISK, CMC, MK and SDP are fitted with the number of subsets below $20$ to ensure that the sample size in each subset $m$ is bigger than the number of knots $r$. Focusing on the estimation of $w_0(\sbb^*)$ for $\sbb^* \in \Scal^*$, CMC and DISK have smaller biases and larger variances than their non-distributed counterparts if $k \leq 30$ (Table <ref>). The DISK estimator's variance decreases and bias increases with increasing $k$, resulting in a decreasing $L_{2}-$risk in the estimation of $w_0$ initially and increasing after $k = 20$ (Figure <ref>), which empirically verifies Bayesian bias-variance trade-off revealed in our theory. The point-wise coverage of 95% CIs in DISK are similar to that of the non-distributed methods, except MPP and NNGP, for $k \leq 30$ and are above the nominal value for all $k$. On the other hand, coverage of CMC and MPP are below the nominal value for every $k$ and $r$ and NNGP fails to cover $w_0$ across all replications. The length of 95% CIs in DISK (GP) and DISK (MPP) are very close to that of their non-distributed version, whereas CMC's and NNGP's CIs greatly underestimate the posterior uncertainty. DISK (GP) and DISK (MPP) with $k=20$ are among the best performers, while MK (MPP) exhibits higher Bayes $L_2$-risk and wider 95% CIs' with different choices of $r$ and $k$. Since estimates of $w_0(\sbb^*)$ are not directly obtained for laGP or LatticeKrig from the  and  packages respectively, we subtract the estimated fixed effects from the predicted values to provide a rough estimate of $w_0(\sbb^*)$ for both of them. While such a strategy yields reasonable point estimates of $w_0(\sbb^*)$ from both of them, characterization of uncertainty of $w_0(\sbb^*)$ is perhaps less unreliable. Hence we refrain from interpreting coverage and length of laGP and LatticeKrig in Table <ref> any further, rather investigate the predictive performance of these methods vis-a-vis other competitors. All methods perform well in terms of MSPE and coverages and lengths of 95% PIs (Table <ref>). DISK (GP) and DISK (MPP) (with $k=10,20$) are more precise in the estimation of $\betab$ compared to their non-distributed versions. The most closely related competitor MK (MPP) shows competitive performance in terms of prediction, but much wider credible intervals for parameters, perhaps due to not accounting for stochastic approximation in subsets. Other competitors exhibit a similar degree of accuracy in estimating $\tau^2$, $\sigma^2$, and $\phi$ with the accuracy of parameter estimates and uncertainty quantification suffering beyond $k\geq 30$.

An interesting feature of our comparisons is the dramatic difference between the performance of DISK (MPP) and MPP with the same choices of knots (Figure <ref>). The performance of MPP in spatial surface and parameter estimation  using full data is sensitive to the choice of $r$, suffering greatly when $r=200$. Contrary to this, DISK (MPP) has substantially smaller $L_2$-risk with the same number of knots used in each subset and its  $w_0$ estimate is almost indistinguishable from the true spatial surface. The performance of MPP is unstable when $r/n$ is low due to poorly conditioned covariance matrix. While the poor performance of full data MPP is attributed to this fact, using $r$ knots in each subset of size $m$ results in relatively high $r/m$ ratio in each subset of DISK (MPP), which contributes to its strikingly superior performance.

While running NNGP with the  package, the inference from NNGP marginalized over $w(\sbb)$'s using the “response” option closely matches DISK in terms of inference on parameter estimates; however, $\wb^*$ cannot be estimated with “response” option, so we employ un-marginalized NNGP using the “sequential” option in the package. This results in severe auto-correlation among the latent variables and $\beta$, yielding NNGP's poor performance across all three choices of the number of nearest neighbors. In an ongoing work, this issue is addressed by employing conjugate gradient algorithms to estimate latent variables in NNGP.

The spatial surface $w_0$ at the locations in $\Scal^*$ for all the competing full Bayesian methods (except NNGP) in Simulation 1.
The 2.5%, 50%, and 97.5% quantile surfaces, respectively, represent pointwise quantiles of the posterior distribution for $w(\sbb^*)$ for every $\sbb^* \in \Scal^*$, where
the 50% quantile of $w(\sbb^*)$ is the estimate of $w_0(\sbb^*)$ and the 2.5% and 97.5% quantiles quantify uncertainty. The true spatial surface $w_0(\sbb^*), \sbb^* \in \Scal^*$ is in the first row and column. The remaining entries in the first column are the estimates of $w_0(\sbb^*), \sbb^* \in \Scal^*$ obtained using the full data posterior distributions with full-rank GP prior and MPP prior with $r=200, 400$. The remaining entries in the first row are the estimates of $w_0(\sbb^*), \sbb^* \in \Scal^*$ obtained using DISK with GP prior and $k=10, 30, 40$, respectively. All other entries provide point estimates and quantify uncertainty in inference on $w_0(\sbb^*), \sbb^* \in \Scal^*$ with $k=20$ and GP prior (second row), MPP prior with rank 200 (third row) and MPP prior with rank 400 (fourth row).

The empirical estimate of bias, variance, and Bayes $L_2$-risk in estimating the spatial surface $w_0$ at the locations in $\Scal^*$  in Simulation 1.
The GP is the theoretical benchmark and laGP is the state-of-the-art method for estimation. The solid lines represent averaged values and the dotted lines and arrows represent one standard deviation error for the full-GP over 10 simulation replications.

Inference on the values of spatial surface at the locations in $\Scal^*$ in Simulation 1. The numbers in parentheses are standard deviations over 10 simulation replications. The bias, variance, and Bayes $L_2$-risk in the estimation of $w_0$ are defined in (<ref>) and the coverage and credible intervals are calculated pointwise for the locations in $\Scal^*$

Bias$^2$     Variance    $L_2$-Risk     95% CI  Coverage     95% CI Length

laGP     0.0004 (0.0000)     0.0100 (0.0002)     0.0103 (0.0002)     1.0000 (0.0000)     0.3890 (0.0036)

LatticeKrig     0.0002 (0.0000)     0.0003 (0.0000)     0.0005 (0.0000)     0.9867 (0.0033)     0.0703 (0.0006)

GP     0.0077 (0.0049)     0.0103 (0.0002)     0.0180 (0.0049)     1.0000 (0.0002)     0.3943 (0.0036)

MPP ($r=200$)     0.3732 (0.3671)     0.0110 (0.0002)     0.3842 (0.3671)     0.0000 (0.0000)     0.4061 (0.0036)

MPP ($r=400$)     0.0623 (0.0369)     0.0105 (0.0002)     0.0727 (0.0370)     0.2946 (0.4662)     0.3976 (0.0037)

NNGP     5c|

NN$ = 5$     0.4213 (0.1373)     0.0021 (0.0002)     0.4233 (0.1373)     0.0000 (0.0000)     0.1778 (0.0079)

NN$ = 15$     0.4822 (0.0666)     0.0013 (0.0001)     0.4835 (0.0666)     0.0000 (0.0000)     0.1421 (0.0067)

NN$ = 25$     0.4887 (0.0668)     0.0013 (0.0001)     0.4900 (0.0668)     0.0000 (0.0000)     0.1398 (0.0032)

CMC (MPP)     5c|

$r=200$, $k=10$     0.0020 (0.0006)     0.0416 (0.0005)     0.0436 (0.0006)     0.8854 (0.0527)     0.1429 (0.0010)

$r=200$, $k=20$     0.0090 (0.0029)     0.0402 (0.0006)     0.0493 (0.0031)     0.1265 (0.1027)     0.1026 (0.0009)

CMC (MPP)     5c|

$r=400$, $k=10$     0.0031 (0.0013)     0.0424 (0.0006)     0.0455 (0.0014)     0.7710 (0.1315)     0.1398 (0.0013)

$r=400$, $k=20$     0.0013 (0.0005)     0.0409 (0.0006)     0.0422 (0.0009)     0.8255 (0.0987)     0.1005 (0.0009)

DISK (GP)     5c|

$k=10$     0.0012 (0.0007)     0.0160 (0.0007)     0.0332 (0.0018)     1.0000 (0.0000)     0.4971 (0.0104)

$k=20$     0.0008 (0.0005)     0.0106 (0.0004)     0.0221 (0.0005)     1.0000 (0.0000)     0.4041 (0.0070)

$k=30$     0.0055 (0.0016)     0.0089 (0.0001)     0.0232 (0.0015)     1.0000 (0.0000)     0.3694 (0.0026)

$k=40$     0.0136 (0.0019)     0.0085 (0.0001)     0.0306 (0.0018)     0.9946 (0.0048)     0.3612 (0.0026)

$k=50$     0.0255 (0.0021)     0.0086 (0.0001)     0.0427 (0.0020)     0.7949 (0.0572)     0.3626 (0.0022)

DISK (MPP)     5c|

$r=200$, $k=10$     0.0017 (0.0008)     0.0189 (0.0009)     0.0394 (0.0021)     1.0000 (0.0000)     0.5388 (0.0122)

$r=200$, $k=20$     0.0009 (0.0004)     0.0131 (0.0002)     0.0270 (0.0004)     1.0000 (0.0000)     0.4477 (0.0039)

DISK (MPP)     5c|

$r=400$, $k=10$     0.0015 (0.0008)     0.0177 (0.0007)     0.0369 (0.0017)     1.0000 (0.0000)     0.5211 (0.0099)

$r=400$, $k=20$     0.0007 (0.0004)     0.0118 (0.0002)     0.0243 (0.0003)     1.0000 (0.0000)     0.4253 (0.0031)

MK (MPP)     5c|

$r=200$, $k=10$     0.0410 (0.0224)     0.3777 (0.0293)     0.4196 (0.0373)     1.0000 (0.0000)     2.4180 (0.0897)

$r=200$, $k=20$     0.0234 (0.0081)     0.4064 (0.0190)     0.4298 (0.0209)     1.0000 (0.0000)     2.5139 (0.0628)

MK (MPP)     5c|

$r=400$, $k=10$     0.0247 (0.0113)     0.3871 (0.0200)     0.4118 (0.0225)     1.0000 (0.0000)     2.4576 (0.0641)

$r=400$, $k=20$     0.0143 (0.0003)     0.4254 (0.0148)     0.4398 (0.0167)     1.0000 (0.0000)     2.5794 (0.0473)

Parametric inference and prediction in Simulation 1. For parametric inference posterior medians are provided along with
the 95% credible intervals (CIs) in the parentheses, where available. Similarly mean squared prediction errors (MSPEs) along with length and coverage of 95% predictive intervals (PIs) are presented, where available. The upper and lower quantiles of 95% CIs and PIs are averaged over 10 simulation replications, with the numbers in parentheses for the last three columns denoting standard deviations across replications;  `-' indicates that the parameter estimate or prediction is not provided by the software or the competitor

$\beta$     $\sigma^2$     $\tau^2$     $\phi$     MSPE     Coverage     Length

Truth      1.00       -         0.01        -      -     -     -

laGP       -       -         -        -     0.010 (0.000)     0.94 (0.01)     0.39 (0.00)

LatticeKrig         -       -        -     -     0.010 (0.000)     0.95 (0.01)     0.39 (0.00)

GP     1.08 (0.50, 1.65)     0.12 (0.10, 0.14)     0.009 (0.009, 0.010)     0.115 (0.107, 0.135)     0.010 (0.000)     0.95 (0.01)     0.39 (0.00)

MPP ($r=200$)     1.56 (0.99, 2.15)     0.15 (0.13, 0.18)     0.008 (0.007, 0.008)     0.119 (0.110, 0.133)     0.010 (0.000)     0.95 (0.01)     0.41 (0.00)

MPP ($r=200$)     1.23 (0.61, 1.84)     0.16 (0.13, 0.19)     0.008 (0.008, 0.008)     0.120 (0.110, 0.148)     0.010 (0.000)     0.95 (0.01)     0.40 (0.00)

NNGP     7c|

NN$=5$     0.36 (0.36, 0.36)     0.29 (0.29, 0.29)     0.009 (0.009, 0.009)     0.123 (0.123, 0.123)     0.011 (0.000)     0.94 (0.01)     0.39 (0.02)

NN$=15$     0.31 (0.31, 0.31)     0.17 (0.17, 0.17)     0.009 (0.009, 0.009)     0.113 (0.113, 0.113)     0.010 (0.000)     0.95 (0.01)     0.40 (0.01)

NN$=25$     0.30 (0.30, 0.30)     0.16 (0.16, 0.16)     0.009 (0.009, 0.009)     0.112 (0.112, 0.112)     0.010 (0.000)     0.95 (0.01)     0.40 (0.01)

CMC (MPP)     7c|

$r=200$, $k=10$     0.99 (0.74, 1.23)     0.23 (0.22, 0.25)     0.006 (0.005, 0.006)     0.112 (0.107, 0.119)     0.011 (0.000)     0.50 (0.01)     0.14 (0.00)

$r=200$, $k=20$     1.09 (0.89, 1.28)     0.29 (0.27, 0.31)     0.007 (0.006, 0.007)     0.109 (0.106, 0.114)     0.011 (0.000)     0.38 (0.01)     0.10 (0.00)

CMC (MPP)     7c|

$r=400$, $k=10$     1.04 (0.77, 1.30)     0.26 (0.24, 0.28)     0.006 (0.006, 0.007)     0.108 (0.105, 0.114)     0.011 (0.000)     0.49 (0.01)     0.14 (0.00)

$r=400$, $k=20$     1.02 (0.82, 1.22)     0.29 (0.27, 0.32)     0.007 (0.007, 0.007)     0.112 (0.109, 0.119)     0.010 (0.000)     0.38 (0.01)     0.10 (0.00)

SDP (MPP)     7c|

$r=200$, $k=10$     0.98 (0.75, 1.23)     0.23 (0.22, 0.25)     0.006 (0.005, 0.006)     0.112 (0.106, 0.118)     -     -     -

$r=200$, $k=20$     1.08 (0.89, 1.27)     0.29 (0.27, 0.31)     0.007 (0.006, 0.007)     0.109 (0.106, 0.113)     -     -     -

SDP (MPP)     7c|

$r=400$, $k=10$      1.04 (0.79, 1.29)     0.26 (0.24, 0.28)     0.006 (0.006, 0.007)     0.109 (0.104, 0.113)     -     -     -

$r=400$, $k=20$      1.02 (0.83, 1.21)     0.30 (0.27, 0.32)     0.007 (0.007, 0.007)     0.113 (0.108, 0.118)     -     -     -

DISK (GP)     7c|

$k=10$     1.03 (0.80, 1.26)     0.21 (0.17, 0.24)     0.009 (0.008, 0.009)     0.124 (0.111, 0.147)     0.010 (0.000)     0.95 (0.01)     0.41 (0.00)

$k=20$     0.98 (0.82, 1.15)     0.22 (0.17, 0.26)     0.008 (0.008, 0.009)     0.142 (0.121, 0.179)     0.010 (0.000)     0.96 (0.01)     0.42 (0.00)

$k=30$     0.93 (0.80, 1.07)     0.20 (0.16, 0.24)     0.008 (0.008, 0.008)     0.171 (0.140, 0.219)     0.010 (0.000)     0.96 (0.00)     0.43 (0.00)

$k=40$     0.88 (0.78, 1.00)     0.18 (0.15, 0.23)     0.008 (0.007, 0.008)     0.201 (0.162, 0.252)     0.010 (0.000)     0.97 (0.01)     0.44 (0.00)

$k=50$     0.84 (0.75, 0.94)     0.17 (0.14, 0.21)     0.007 (0.007, 0.008)     0.231 (0.186, 0.285)     0.011 (0.000)     0.97 (0.00)     0.46 (0.00)

DISK (MPP)     7c|

$r=200$, $k=10$     1.03 (0.80, 1.27)     0.21 (0.18, 0.24)     0.009 (0.008, 0.009)     0.120 (0.109, 0.144)     0.010 (0.000)     0.97 (0.01)     0.44 (0.00)

$r=200$, $k=20$     0.98 (0.82, 1.16)     0.22 (0.17, 0.26)     0.008 (0.008, 0.009)     0.140 (0.119, 0.177)     0.010 (0.000)     0.97 (0.01)     0.46 (0.00)

DISK (MPP)     7c|

$r=400$, $k=10$     1.03 (0.80, 1.27)     0.21 (0.18, 0.24)     0.009 (0.008, 0.009)     0.119 (0.109, 0.143)     0.010 (0.000)     0.96 (0.01)     0.42 (0.00)

$r=400$, $k=20$     0.98 (0.82, 1.16)     0.22 (0.17, 0.26)     0.008 (0.008, 0.009)     0.140 (0.119, 0.181)     0.010 (0.000)     0.97 (0.01)     0.44 (0.00)

MK (MPP)     7c|

$r=200$, $k=10$     1.12 (0.03, 2.28)     0.41 (0.21, 0.58)     0.007 (0.004, 0.008)     0.058 (0.051, 0.127)     0.010 (0.000)     0.96 (0.01)     0.47 (0.00)

$r=200$, $k=20$     1.16 (0.04, 2.33)     0.46 (0.27, 0.66)     0.006 (0.004, 0.008)     0.060 (0.051, 0.112)     0.011 (0.000)     0.97 (0.01)     0.49 (0.00)

MK (MPP)     7c|

$r=400$, $k=10$     1.12 (0.07, 2.27)     0.41 (0.22, 0.60)     0.007 (0.005, 0.008)     0.068 (0.052, 0.129)     0.010 (0.000)     0.96 (0.00)     0.44 (0.00)

$r=400$, $k=20$     1.12 (-0.10, 2.38)     0.49 (0.26, 0.71)     0.007 (0.005, 0.009)     0.063 (0.050, 0.125)     0.010 (0.000)     0.97 (0.00)     0.48 (0.00)

§.§.§ Simulation 2: moderately large data with a rough spatial surface.

Covariance function of DISK-MPP as a function distance with 200 and 400 knots (r), respectively. The true covariance function is $\cov\{w_0(\sbb_i), w_0(\sbb_j)\} = \exp(-10 \|\sbb_i - \sbb_j \|_2)$.

Our second simulation example provides performance of DISK (MPP) with various choices of the number of knots when data are simulated from a Gaussian random field with nowhere differentiable  surface.  Based on the best results in Simulation 1, we use DISK (MPP) with $k=20$ and laGP as our only competitor. The bias, variance, and $L_2$-risk estimates of DISK (MPP) and laGP show a similar pattern as in Simulation 1 (Table <ref>). For both DISK (MPP) with $r=200$ and $r=400$ knots, the coverage of 95% CI turn out to be nominal. The posterior median of all parameters are close to the true value with 95% CIs covering the true value for $\beta,\sigma^2,\phi$, and the same is true for predictions and 95% predictive intervals; however, DISK (MPP) slightly overestimates $\tau^2$ (Table <ref>). This is expected given that the MPP prior applied to the full data tends to overestimate $\tau^2$.  The predictive inference of DISK shows desirable point prediction with precise characterization of uncertainty. We also emphasize that the full Bayesian inference from DISK (MPP) is computationally extremely efficient and takes less than 2 hours; see supplementary materials for detailed comparisons.

The DISK posterior of the covariance function from modified predictive process is plotted and contrasted with the true correlation function in Figure <ref>. While earlier articles [49] point out discrepancy between estimated correlation function of the modified predictive process and the true data generating correlation function, Figure <ref> shows that DISK posterior of the correlation function recovers the true correlation quite accurately. Moreover, the accuracy increases once we employ more knots in each subset. Table <ref> reveals that laGP offers better $L_2$-risk estimates with comparable quantification of uncertainties for  surface estimation; however, it needs to be emphasized that while running laGP using the available package, the $\beta$ parameter is fixed at its true value. In fact, as noted before, laGP does not offer posterior estimates of parameters which are readily available from DISK. Additionally, strong local variability in the generated nowhere differentiable true spatial surface naturally prefers a nearest neighbor approach over a low-rank approach. Since DISK (MPP) can be conceptualized as a computationally convenient alternative to MPP, a low-rank approach, with large number of knots, it is understandable that it may appear to be less effective in estimating the  surface than laGP in this case. Nevertheless, the model free nature of the DISK approach allows us to fit a nearest neighbor approach, including NNGP, on each subset to improve inference and expedite computations by multiple folds. We plan to carefully investigate this DISK feature elsewhere.

Inference on the values of spatial surface at the locations in $\Scal^*$ in Simulation 2.
The numbers in parentheses are standard deviations over 10 simulation replications. The bias, variance, and Bayes $L_2$-risk in the estimation of $w_0$ are defined in (<ref>) and the coverage and length of 95% credible intervals are calculated pointwise for the locations in $\Scal^*$

Bias$^2$     Variance    $L_2$-Risk     Coverage     Length

laGP                       0.4059 (0.0130)     0.4910 (0.0086)     0.8969 (0.0189)     0.9670 (0.0029)     2.7216 (0.0237)

DISK (MPP),$k=20,r=200$     0.8133 (0.0264)     0.8791 (0.0619)     1.6923 (0.0793)     0.9598 (0.0067)     3.6735 (0.1311)

DISK (MPP),$k=20,r=400$     0.7295 (0.0256)     0.8219 (0.0571)     1.5515 (0.0756)     0.9645 (0.0066)     3.5531 (0.1245)

Parameter estimates along with 95% credible intervals from DISK(MPP) with $r=200$ and $r=400$ in Simulation 2. For parametric inference posterior medians are provided along with the 95% credible intervals (CIs) in the parentheses , where available. Similarly mean squared prediction errors (MSPEs) along with length and coverage of 95% predictive intervals (PIs) are presented, where available. The upper and lower quantiles of 95% CIs and PIs are averaged over 10 simulation replications, with the numbers in parentheses for the last three columns denoting standard deviations across replications;  `-' indicates that the parameter estimate or prediction is not provided by the software or the competitor.

$\beta$     $\sigma^2$     $\tau^2$     $\phi$     MSPE     Coverage     Length     $\log_{10}$(Time)

Truth     1     1     0.1     9     -     -     -     -

laGP     -       -       -         -      0.50 (0.0120)     0.21 (0.0115)     0.38 (0.0003)     0.79 (0.0193)

DISK, $k=20,r=200$     1.00 (0.96, 1.03)     0.97 (0.92, 1.03)     0.14 (0.12, 0.17)     9.28 (8.94, 9.68)     0.90 (0.0210)     0.96 (0.0027)     3.97 (0.0773)     2.39 (0.0011)

DISK, $k=20,r=400$     1.00 (0.96, 1.03)     0.97 (0.92, 1.03)     0.14 (0.11, 0.17)     9.32 (8.96, 9.73)     0.82 (0.0194)     0.96 (0.0035)     3.86 (0.0717)     2.51 (0.0262)

§.§.§ Simulation 3: large spatial data.

Our ultimate goal is to apply DISK in massive data settings, so we evaluate its performance when $n=10^6$. As mentioned earlier, massive size of the data in Simulation 3 prohibits the fitting of models based on full-rank and low-rank GPs, including MPP, LatticeKrig, and NNGP due to numerical issues, leaving only laGP as a feasible competitor. Since Simulation 1 demonstrates similar performance of DISK (MPP) and DISK (GP) with DISK (MPP) having a smaller run time, we use only DISK (MPP) for comparisons with laGP in Simulation 3.
An identical three-step strategy for fitting DISK (MPP) is employed as in Simulation 1 but with $k=500$ and with $r=400$ and $r=600$. Notably, $r$ is increased from $200$ and $400$ in Simulation 1 to $400$ and $ 600$ to account for the larger subset size in this simulation, maintaining a high $r/m$ ratio.

The results for DISK's uncertainty quantification in parameter estimation and prediction agree with those observed in Simulation 1, but, unlike Simulation 1, DISK outperforms laGP in the estimation of $w_0(\sbb^*)$ for $\sbb^* \in \Scal^*$ for both choices of $r$ (Tables <ref> and <ref>). The point estimates of $\beta$ and $\tau^2$ are close to their true values with narrow 95% CIs. For $\tau^2$, the CI misses the truth, which is expected given that the full data GP in Simulation 1 underestimates $\tau^2$. The bias, variance, and Bayes $L_2$-risk of DISK (MPP) for both $r$s are lower than those of laGP. The coverages of 95% CIs for laGP and DISK (MPP) are the same but the lengths of 95% CIs in DISK (MPP) are smaller than those of laGP for both $r$s (see Table <ref>). We conclude that DISK (MPP) matches the performance of laGP in delivering predictive inference, while it outperforms laGP in terms of estimation of $w_0(\sbb^*)$ for $\sbb^* \in \Scal^*$.

Inference on the values of spatial surface at the locations in $\Scal^*$ in Simulation 3. The numbers in parentheses are standard deviations over 10 simulation replications. The bias, variance, and $L_2$-risk in the estimation of $w_0$ are defined in (<ref>) and the coverage and length of 95% credible intervals are calculated pointwise for the locations in $\Scal^*$.

Bias$^2$     Variance    $L_2$-Risk     Coverage     Length

laGP     0.0002 (0.0000)     0.0100 (0.0000)     0.0102 (0.0000)     1.0000 (0.0000)     0.3905 (0.0006)

DISK, $k=500$, $r=400$     0.0002 (0.0000)     0.0030 (0.0000)     0.0061 (0.0000)     1.0000 (0.0000)     0.2132 (0.0002)

DISK, $k=500$, $r=600$     0.0001 (0.0000)     0.0026 (0.0000)     0.0052 (0.0000)     1.0000 (0.0000)     0.1977 (0.0002)

Parametric inference and prediction in Simulation 3. For parametric inference posterior medians are provided along with
the 95% credible intervals (CIs) in the parentheses, where available. Similarly mean squared prediction errors (MSPEs) along with length and coverage of 95% predictive intervals (PIs) are presented, where available. The upper and lower quantiles of 95% CIs and PIs are averaged over 10 simulation replications, with the numbers in parentheses for the last three columns denoting standard deviations across replications;  `-' indicates that the parameter estimate or prediction is not provided by the software or the competitor.

$\beta$     $\sigma^2$     $\tau^2$     $\phi$

Truth      1.00       -         0.01        -

DISK (MPP) ($r=400$, $k=500$)      1.01 (0.98, 1.04)     0.16 (0.15, 0.17)     0.008 (0.008, 0.008)     0.13 (0.13, 0.13)

DISK (MPP) ($r=600$, $k=500$)      1.01 (0.98, 1.04)     0.16 (0.15, 0.17)     0.008 (0.008, 0.008)     0.13 (0.13, 0.13)

MSPE     Coverage     Length

laGP                              0.010 (0.0003)     0.94 (0.0040)     0.39 (0.0000)

DISK (MPP) ($r=400$, $k=500$)       0.009 (0.0003)     0.96 (0.0030)     0.41 (0.0000)

DISK (MPP) ($r=600$, $k=500$)       0.009 (0.0002)     0.95 (0.0040)     0.40 (0.0000)

§.§ Real data: Sea Surface Temperature data

A description of the evolution and dynamics of the SST is a key component of the study of the earth's climate. SST data  (in centigrade) from ocean samples have been collected by voluntary observing ships, buoys, and military and scientific cruises for decades. During the last 20 years or so, the SST database has been complemented by regular streams of remotely sensed observations from satellite orbiting the earth. A careful quantification of variability of SST data is important for climatological research, which includes determining the formation of sea breezes and sea fog and calibrating measurements from weather satellites [15]. A number of articles have appeared to address this issue in recent years; see [8, 36, 67].

We consider the problem of capturing the spatial trend and characterizing the uncertainties in the SST in the west coast of mainland U.S.A., Canada, and Alaska between $40^\circ$–$65^\circ$ north latitudes and $100^\circ$–$180^\circ$ west longitudes. The dataset is obtained from NODC World Ocean Database (<https://www.nodc.noaa.gov/OC5/WOD/pr_wod.html>). Due to our focus on spatial modeling, we ignore the temporal component. After screening the data for quality control, we choose a random subset of about $1,000,800$ spatial observations over the selected domain. From the selected observations, we randomly select $10^6$ observations as training data and the remaining observations are used to compare the performance of DISK and its competitors.  We replicate this setup ten times. The selected domain is large enough to allow considerable spatial variation in SST from north to south and provides an important first step in extending these models for analyzing global-scale SST database.

The SST data in the selected domain shows a clear decreasing trend in SST with increasing latitude (Figure <ref>). Based on this observation, we add latitude as a linear predictor in the univariate spatial regression model (<ref>) to explain the long-range directional variability in the SST.

The setup is identical to Simulation 2, except for the presence of the latitude predictor with the corresponding coefficient $\beta_1$.
We assign $N(0, 100)$ prior to $\beta_1$, and the remaining priors and DISK competitors are identical to those in Simulation 2. For application of CMC (MPP), DISK (MPP), and SDP (MPP), we follow the three-step strategy used in Simulation 1 but with $k=300$ and $r=400$. Due to the lack of ground truth for estimating $w_0(\sbb^*)$, we compare the four methods in terms of their inference on $\Omegab$ and prediction of $y(\sbb^*)$ for $\sbb^* \in \Scal^*$ in terms of MSPE and the length and coverage of 95% posterior PIs.

DISK (MPP) outperforms CMC (MPP) and SDP (MPP) in predictions while closely matching the results of laGP, the current state-of-the-art method for modeling massive spatial data.
The  50%, 2.5%, and 97.5% quantiles of the posterior distributions for $\Omegab$, $w(\sbb^*)$ and $y(\sbb^*)$ for every $\sbb^* \in \Scal^*$ are used for estimation and  uncertainty quantification. CMC (MPP), SDP (MPP) and DISK (MPP)  agree closely in their inference on $\Omegab$, but SDP (MPP) fails to provide any result for $w(\sbb^*)$ or $y(\sbb^*)$ due to the large size of $\Scal^*$ (Table <ref>). For every $\sbb^* \in \Scal^*$, CMC's and DISK's estimates of $w(\sbb^*)$ and $y(\sbb^*)$ agree closely, but CMC severely underestimates uncertainty in $w(\sbb^*)$ and $y(\sbb^*)$ (Figures <ref> and <ref> and Table <ref>). The pointwise predictive coverages of laGP and DISK match their nominal levels; however, the 95% posterior PIs of DISK are wider than those of laGP because DISK accounts for uncertainty due to the error term (Figure <ref> and Table <ref>).

As a whole, SST data analysis reinforces our findings on DISK as a computationally efficient, flexible, and fully Bayesian inferential tool.

Parametric inference and prediction in SST data. CMC, SDP, and DISK use  MPP-based  modeling with $r=400$ on $k=300$ subsets.
For parametric inference posterior medians are provided along with
The 95% credible intervals (CIs) in the parentheses, where available. Similarly mean squared prediction errors (MSPEs) along with length and coverage of 95% predictive intervals (PIs) are presented, where available. The upper and lower quantiles of 95% CIs and PIs are averaged over 10 simulation replications, with the numbers in parentheses for the last three columns denoting standard deviations across replications;  `-' indicates that the parameter estimate or prediction is not provided by the software or the competitor

$\beta_0$     $\beta_1$     $\sigma^2$     $\tau^2$     $\phi$

laGP     -     -     -     -     -

CMC     31.78 (31.19, 32.37)     -0.35 (-0.36, -0.34)     12.22 (11.78, 12.69)     0.110 (0.108, 0.112)     0.021 (0.020, 0.022)

SDP     31.67 (31.45, 31.82)     -0.35 (-0.36, -0.34)     14.42 (13.29, 14.80)     0.110 (0.108, 0.112)     0.021 (0.020, 0.022)

DISK     32.34 (31.74, 32.95)     -0.32 (-0.33, -0.31)     11.83 (11.23, 12.45)     0.184 (0.182, 0.185)     0.039 (0.037, 0.041)

MSPE     Coverage     Length

laGP     0.25 (0.00)     0.95 (0.00)     2.35 (0.00)

CMC     0.41 (0.00)     0.13 (0.00)     0.14 (0.00)

SDP     -     -     -

DISK     0.41 (0.00)     0.95 (0.00)     2.67 (0.00)

Predication of sea surface temperatures at the locations in $\Scal^*$. Negative longitude means degree west from Greenwich.
CMC and DISK use  MPP-based  modeling with $r=400$ on $k=300$ subsets and laGP uses the `nn' method. The 2.5%, 50%, and 97.5% quantile surfaces, respectively, represent pointwise quantiles of the posterior distribution for $y(\sbb^*)$ for every $\sbb^* \in \Scal^*$.

Interpolated spatial surface $w$ at the locations in $\Scal^*$. Negative longitude means degree west from Greenwich.
CMC and DISK use  MPP-based  modeling with $r=400$ on $k=300$ subsets. The 2.5%, 50%, and 97.5% quantile surfaces, respectively, represent pointwise quantiles of the posterior distribution for $w(\sbb^*)$ for every $\sbb^* \in \Scal^*$.

§ DISCUSSION

This article presents a novel distributed Bayesian approach for kriging with massive data using the divide-and-conquer technique.

We provide explicit upper bound on the number of subsets $k$ depending on the analytic properties of the spatial surface, so that the Bayes $L_2$-risk of the DISK posterior is nearly minimax optimal. We have confirmed this empirically via simulated and real data analyses, where DISK compares well with state-of-the-art methods. Additional theoretical results in the supplementary material shed light on the posterior convergence rate of the DISK posterior.

The simplicity and generality of the DISK framework enable scaling of any spatial model. For example, recent applications have confirmed that the NNGP prior requires modifications if scalability is desired for even a few millions of locations [18]. In future, we aim to scale ordinary NNGP and other multiscale approaches to tens of millions of locations with the DISK framework. Another important future work is to extend the DISK framework for scalable modeling of multiple correlated outcomes observed over massive number of locations.

This article focuses on developing the DISK framework for spatial modeling due to the motivating applications from massive geostatistical data. The DISK framework, however, is applicable to any mixed effects model where the random effects are assigned a GP prior, which includes Bayesian nonparametric regression using GP prior.

We plan to explore more general applications in the future with high dimensional covariates.

§ ACKNOWLEDGMENT

We thank Professor David B. Dunson of Duke University for inspiring many questions that we have addressed in this work and Professor Sudipto Banerjee of UCLA for helpful conversations. Cheng Li's research is supported by the National University of Singapore start-up grant (R155000172133). This material is based on research supported by the Office of Naval Research under Award Number ONR-BAA N00014-18-1-2741.

Supplementary Material for A Divide-and-Conquer Bayesian Approach to Large-Scale Kriging

§ PROOF OF THEOREMS IN SECTION 3.4

Recall that the spatial regression model with a GP prior considered in Section 2 is

\begin{align}
y(\sbb_i)= w(\sbb_i) + \epsilon(\sbb_i), \quad \epsilon(\sbb_i) \sim N(0, \tau^2), \quad i = 1, \ldots, n, \quad  w(\cdot) \sim \text{GP}\{0, C_{\alphab}(\cdot, \cdot) \}. \label{atheory-mdl}
\end{align}

Writing this model for the $n$ locations in $\Scal$ gives

\begin{align}
\label{eq:ap1}
\yb = \wb_0 + \epsilonb, \quad \epsilonb \mid \Scal \sim N(\zero, \tau^2 \Ib), \quad \yb \mid \Scal \sim N(\wb_0, \tau^2 \Ib),
\end{align}

where $\wb_0 = \{w_0(\sbb_1), \ldots, w_0(\sbb_n)\}$ and $\epsilonb = \{\epsilon(\sbb_1), \ldots, \epsilon(\sbb_n)\}$ are the true value of the residual spatial surface and white noise realized at the locations in $\Scal$. We can write the model in a similar format for each data subset. Let $\sbb \in \Dcal$ be a location, $w_0(\sbb)$ be the true value of the residual spatial surface, $\EE_{\sbb^*}$, $\EE_{0}$, $\EE_{\Scal}$, and $\EE_{0 \mid \Scal}$  respectively be the expectations with respect to the distributions of $\sbb^*$, $(\Scal, \yb)$, $\Scal$, and $\yb$ given $\Scal$.
If $\overline w(\sbb^*)$ is a random variable that follows the DISK posterior for estimating $w_0(\sbb^*)$, then $\overline w(\sbb^*)$ has the density $N(\overline m, \overline v)$, where

\begin{align}
\label{eq:ae3}
\overline m = \frac{1}{k} \sum_{j=1}^k \cb_{j,*}^T (\Cb_{j, j} + \tfrac{\tau^2}{k} \Ib)^{-1}\yb_j , \;
\overline v^{1/2} = \frac{1}{k}\sum_{j=1}^k v_j^{1/2}, \;
v_j = c_{*,*} - \cb_{j,*}^T (\Cb_{j,j} + \tfrac{\tau^2}{k} \Ib)^{-1} \cb_{j, *},
\end{align}

where $c_{*,*}=C_{\alphab}(\sbb^*,\sbb^*)$, and $\cb_{j,*}^T = \cb_j^T(\sbb^*) = [C_{\alphab}(\sbb_{j1}, \sbb^*), \ldots, C_{\alphab}(\sbb_{jm}, \sbb^*)]$. In the proofs below, without confusion, we use the notation $\cb_{j,*}$ and $\cb_j(\sbb^*)$ interchangeably.

The Bayes $L_2$-risk in estimating $w_0$ using the DISK posterior is defined as

\begin{align} \label{eq:ap2}
\EE_{0} \left[ \EE_{\sbb^*} \{\overline w(\sbb^*) - w_0(\sbb^*)\}^2 \right]  \overset{(i)}{=} \EE_{\Scal}  \int_{\Dcal} \EE_{0 \mid \Scal}\{\overline w(\sbb^*) - w_0(\sbb^*)\}^2 \PP_{\sbb}(d \sbb^*),
\end{align}

where $(i)$ follows from Fubini's theorem. Using bias-variance decomposition,

\begin{align*}
\EE_{0 \mid \Scal} \{ \overline w(\sbb^*) - w_0(\sbb^*)\}^2  &= \EE_{0 \mid \Scal} \left[\overline w(\sbb^*) - \EE_{0 \mid \Scal} \{\overline w(\sbb^*)\} +  \EE_{0 \mid \Scal} \{\overline w(\sbb^*)\} - w_0(\sbb^*) \right]^2\\
&= \left[ \EE_{0 \mid \Scal} \{\overline w(\sbb^*)\} - w_0(\sbb^*) \right]^2  + \EE_{0 \mid \Scal} \left[\overline w(\sbb^*) - \EE_{0 \mid \Scal} \{\overline w(\sbb^*)\}\right]^2   \\
&\equiv  \text{bias}_{0 \mid \Scal}^2\{\overline w(\sbb^*)\} + \var_{0 \mid \Scal}\{\overline w(\sbb^*)\} .
\end{align*}

If $\cb_j^T(\cdot) = [\cov\{w(\cdot), w(\sbb_{j1})\}, \ldots, \cov\{w(\cdot), w(\sbb_{jm})\}] = \{C_{\alphab}(\sbb_{j1}, \cdot), \ldots, C_{\alphab}(\sbb_{jm}, \cdot)\}$,  $\cb^T(\cdot) = \{\cb_1^T(\cdot), \ldots, \cb_k^T(\cdot)\}$, $\wb_{0j}^T = \{w_0(\sbb_{j1}), \ldots, w_0(\sbb_{jm})\}$, and $\wb_{0}^T = \{\wb^T_{01}, \ldots, \wb_{0k}^T\}$,
then the distribution of $\overline w(\sbb^*)$ in (<ref>) implies that

\begin{align*}
\EE_{0 \mid \Scal} \{\overline w(\sbb^*) \}  &= \frac{1}{k} \sum_{j=1}^k \cb_j^T (\sbb^*) \left(\Cb_{j, j} + \tfrac{\tau^2 }{k} I \right)^{-1} \wb_{0j} = \cb_j^T(\sbb^*) (k \Lb + \tau^2 \Ib)^{-1} \wb_0,\\
\var_{0 \mid \Scal}\{\overline w(\sbb^*)\} &= \var_{0 \mid \Scal} \left[ \EE \{\overline w(\sbb^*) \mid \yb\} \right] + \EE_{0 \mid \Scal} \left[ \var \{\overline w(\sbb^*) \mid \yb \} \right] \\
&= \tau^2 \cb^T(\sbb^*) (k \Lb + \tau^2 \Ib)^{-2} \cb(\sbb^*) + \overline v(\sbb^*),
\end{align*}

where $\Lb$ is a block-diagonal matrix with $\Cb_{1,1}, \ldots, \Cb_{k,k}$ along the diagonal;
therefore, the Bayes $L_2$-risk in (<ref>) can be decomposed into three parts:

\begin{align}
\label{eq:eap4}
\EE_{\sbb^*} \EE_{\Scal} \{  \cb_{*}^T(k \Lb + \tau^2 \Ib)^{-1}\wb_0 - w_0(\sbb^*)\}^2 + \tau^2 \EE_{\sbb^*} \EE_{\Scal} \left\{\cb^T_* (k \Lb + \tau^2 \Ib)^{-2} \cb_* \right\}  + \EE_{\sbb^*}\EE_{\Scal}  \overline v(\sbb^*),
\end{align}

which correspond to $\text{bias}^2$, ${\var}_{\text{mean}}$ and ${\var}_{\text{DISK}}$ in Theorem 1.

§.§ Proof of Theorem 1

The next three sections find upper bounds for each of the three terms in (<ref>). The conclusion of Theorem 1 follows directly by combining the three upper bounds.

§.§.§ An upper bound for the squared bias

Consider the squared-bias term in (<ref>). For ease of presentation, assume that $\{\sbb_1, \ldots, \sbb_n\}$ are relabeled to
$$\{\sbb_{11}, \ldots, \sbb_{1m}, \ldots, \sbb_{k1}, \ldots, \sbb_{km}\}$$
corresponding to the $k$ subsets. Define $\xi_{\sbb_{ji}}(\cdot) = C_{\alphab}(\sbb_{ji}, \cdot)$,

\begin{align}
\label{eq:14}