Multi-Scale Process Modelling and Distributed Computation for Spatial Data

TL;DR

This paper introduces a scalable, multi-scale spatial modeling approach using Gaussian Markov random fields and parallel MCMC, enabling efficient analysis of large, complex, nonstationary spatial datasets like sea-surface temperature.

Contribution

It presents a novel multi-scale process model that combines nonstationarity and distributed computation for large spatial data, advancing beyond existing global models.

Findings

Successfully modeled sea-surface temperature with ~1 million observations.

Demonstrated improved scalability and flexibility over state-of-the-art methods.

Enabled distributed computing and data sharing for complex spatial models.

Abstract

Recent years have seen a huge development in spatial modelling and prediction methodology, driven by the increased availability of remote-sensing data and the reduced cost of distributed-processing technology. It is well known that modelling and prediction using infinite-dimensional process models is not possible with large data sets, and that both approximate models and, often, approximate-inference methods, are needed. The problem of fitting simple global spatial models to large data sets has been solved through the likes of multi-resolution approximations and nearest-neighbour techniques. Here we tackle the next challenge, that of fitting complex, nonstationary, multi-scale models to large data sets. We propose doing this through the use of superpositions of spatial processes with increasing spatial scale and increasing degrees of nonstationarity. Computation is facilitated through the use of Gaussian Markov random fields and parallel Markov chain Monte Carlo based on graph colouring. The resulting model allows for both distributed computing and distributed data. Importantly, it provides opportunities for genuine model and data scaleability and yet is still able to borrow strength across large spatial scales. We illustrate a two-scale version on a data set of sea-surface temperature containing on the order of one million observations, and compare our approach to state-of-the-art spatial modelling and prediction methods.

Tables

Table 1. Table 1 : Diagnostic results on the validation data 𝐙 v ( 1 ) superscript subscript 𝐙 𝑣 1 \mathbf{Z}_{v}^{(1)} (validation data in the vicinity of training data), 𝐙 v ( 2 ) superscript subscript 𝐙 𝑣 2 \mathbf{Z}_{v}^{(2)} (validation data not in the vicinity of training data), and 𝐙 v ( 3 ) superscript subscript 𝐙 𝑣 3 \mathbf{Z}_{v}^{(3)} (validation data in the 8 ∘ × 8 ∘ superscript 8 superscript 8 8^{\circ}\times 8^{\circ} box).

Data	Model	RMSPE	CRPS	IS90	Cov90
	FSA	0.353	0.185	1.659	0.915
	LTK	0.378	0.200	1.752	0.913
${𝐙}_v^{\left(1\right)}$	NNGP	0.345	0.180	1.632	0.913
	SPDE0	0.406	0.216	1.889	0.922
	SPDE1-indep	0.346	0.180	1.556	0.919
	SPDE2	0.345	0.176	1.479	0.907
	FSA	1.001	0.540	4.489	0.928
	LTK	0.763	0.430	4.065	0.959
${𝐙}_v^{\left(2\right)}$	NNGP	0.799	0.464	4.420	0.990
	SPDE0	0.699	0.377	3.170	0.928
	SPDE1-indep	0.728	0.382	3.278	0.935
	SPDE2	0.680	0.345	2.666	0.894
	FSA	0.460	0.325	3.570	0.993
	LTK	0.463	0.339	3.748	0.986
${𝐙}_v^{\left(3\right)}$	NNGP	0.515	0.422	5.000	0.996
	SPDE0	0.425	0.282	2.926	0.995
	SPDE1-indep	0.461	0.252	1.830	0.916
	SPDE2	0.367	0.210	1.798	0.956

Equations

$$({\kappa }_0^2-\Delta {)}^{{\alpha }_0/2}({\tau }_0{Y}_0(\text{s}))=W_0(\text{s}),\text{s}\in D,$$

$$({\kappa }_k^2(\text{s})-\Delta {)}^{{\alpha }_k/2}({\tau }_k(\text{s})Y_k(\text{s}))=W_k(\text{s}),\text{s}\in D,$$

$$\mathbf{Z}=\sum _{k=0}^K{\text{A}}_k{\mathbf{\eta }}_k+\mathbf{\epsilon },$$

$$\log {\tau }_k(\cdot )$$

$$\log {\kappa }_k(\cdot )$$

$$\log {\sigma }_{\epsilon }(\cdot )$$

$$\log {\sigma }_k(\cdot )$$

$$\log {\rho }_k(\cdot )$$

$$\frac{{\mathbf{\theta }}_X^{\{i\}}\mid {\omega }_X,{\lambda }_X\text{buildrel}\mathrm{IID}}{\sim \mathrm{Gau}({\omega }_X,{\lambda }_X),i=1,\dots ,r_X,}$$

$$\log {\sigma }_0\sim \mathrm{Gau}({\omega }_{{\sigma }_0},{\lambda }_{{\sigma }_0}),$$

$$\log {\rho }_0\sim \mathrm{Gau}({\omega }_{{\rho }_0},{\lambda }_{{\rho }_0}),$$

$${\mathbf{\theta }}_{\epsilon }\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}$$

$$\mathrm{pr}({\mathbf{\eta }}_0,{\mathbf{\theta }}_0\mid \mathrm{eelse})=$$

$$\times \mathrm{pr}({\mathbf{\theta }}_0\mid \mathrm{eelse}excl.{\mathbf{\eta }}_0),$$

$$\mathrm{pr}({\mathbf{\theta }}_0\mid \mathrm{eelse}excl.{\mathbf{\eta }}_0)$$

$$\propto \int \mathrm{pr}(\mathbf{Z}\mid {\mathbf{\eta }}_0,\mathrm{eelse})$$

$$\times \mathrm{pr}({\mathbf{\eta }}_0\mid {\mathbf{\theta }}_0)\text{d}{\mathbf{\eta }}_0\mathrm{pr}({\mathbf{\theta }}_0).$$

$$\log \mathrm{pr}({\mathbf{\theta }}_0$$

$$-\frac{1}{2}\log |{\tilde{\text{Q}}}_0|-\frac{1}{2}{\tilde{\mathbf{Z}}}^{\top }{\text{Q}}_{\epsilon }\tilde{\mathbf{Z}}+\frac{1}{2}{\tilde{\mathbf{Z}}}^{\top }{\text{Q}}_{\epsilon }{\text{A}}_0{\tilde{\mathbf{\mu }}}_0$$

$$+\log \mathrm{pr}({\mathbf{\theta }}_0)+\text{const.},$$

$${\mathbf{\eta }}_0\mid \mathrm{eelse}\sim \mathrm{Gau}({\tilde{\mathbf{\mu }}}_0,{\tilde{\text{Q}}}_0^{-1}),$$

$${\mathbf{Z}}^{\mathcal{F}}$$

$$\mathrm{pr}({\mathbf{\eta }}_k^{\mathcal{T}},{\mathbf{\theta }}_k^{\mathcal{I}}\mid \mathrm{eelse})=$$

$$\times \mathrm{pr}({\mathbf{\theta }}_k^{\mathcal{I}}\mid \mathrm{eelse}\text{ excl. }{\mathbf{\eta }}_k^{\mathcal{T}}),$$

$$\mathrm{pr}({\mathbf{\theta }}_k^{\mathcal{I}}\mid$$

$$\propto \int \mathrm{pr}({\mathbf{Z}}^{\mathcal{F}}\mid {\mathbf{\eta }}_k^{\mathcal{T}},\mathrm{eelse})\mathrm{pr}({\mathbf{\eta }}_k^{\mathcal{T}}\mid {\mathbf{\eta }}_k^{\mathrm{rest}},{\mathbf{\theta }}_k^{\mathcal{I}},{\mathbf{\theta }}_k^{\mathrm{rest}})$$

$$\times \mathrm{pr}({\mathbf{\eta }}_k^{\mathrm{rest}}\mid {\mathbf{\theta }}_k^{\mathcal{I}},{\mathbf{\theta }}_k^{\mathrm{rest}})\text{d}{\mathbf{\eta }}_k^{\mathcal{T}}\mathrm{pr}({\mathbf{\theta }}_k^{\mathcal{I}}).$$

$$\begin{array}{r}\mathrm{pr}({\mathbf{\eta }}_k^{\mathrm{rest}}|\text{nolinebreak}{\mathbf{\theta }}_k^{\mathcal{I}},{\mathbf{\theta }}_k^{\mathrm{rest}})\approx \mathrm{pr}({\mathbf{\eta }}_k^{\mathrm{rest}}|\text{nolinebreak}{\tilde{\mathbf{\theta }}}_k^{\mathcal{I}},{\mathbf{\theta }}_k^{\mathrm{rest}}),\end{array}$$

$$\mathrm{pr}($$

$$\displaystyle\mathrel{\vbox{ \offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern 2.0pt\cr\sim\cr\kern-2.0pt\cr}}}\int\operatorname{\mathrm{pr}}(\mathbf{Z}^{{\mathcal{F}}}\mid{\bm{\eta}}_{k}^{\mathcal{T}},\mathrm{eelse})\operatorname{\mathrm{pr}}({\bm{\eta}}_{k}^{\mathcal{T}}\mid{\bm{\eta}}_{k}^{\mathrm{rest}},{\bm{\theta}}_{k}^{\mathcal{I}},{\bm{\theta}}_{k}^{\mathrm{rest}})\,\textrm{d}{\bm{\eta}}_{k}^{\mathcal{T}}$$

$$\times \mathrm{pr}({\mathbf{\theta }}_k^{\mathcal{I}}),$$

$${\mathbf{Z}}^{{\mathcal{F}}_1}$$

$$\mathrm{pr}($$

$$\displaystyle~{}\mathrel{\vbox{ \offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern 2.0pt\cr\sim\cr\kern-2.0pt\cr}}}~{}\operatorname{\mathrm{pr}}({\bm{\eta}}_{k}^{{\mathcal{T}}_{1}},{\bm{\theta}}_{k}^{{\mathcal{I}}_{1}}\mid\mathrm{eelse})\operatorname{\mathrm{pr}}({\bm{\eta}}_{k}^{{\mathcal{T}}_{2}},{\bm{\theta}}_{k}^{{\mathcal{I}}_{2}}\mid\mathrm{eelse})$$

$$=\mathrm{pr}({\mathbf{\eta }}_k^{{\mathcal{T}}_1}\mid \mathrm{eelse})\mathrm{pr}({\mathbf{\theta }}_k^{{\mathcal{I}}_1}\mid \mathrm{eelse}\text{ excl. }{\mathbf{\eta }}_k^{{\mathcal{T}}_1})$$

$$\times \mathrm{pr}({\mathbf{\eta }}_k^{{\mathcal{T}}_2}\mid \mathrm{eelse})\mathrm{pr}({\mathbf{\theta }}_k^{{\mathcal{I}}_2}\mid \mathrm{eelse}\text{ excl. }{\mathbf{\eta }}_k^{{\mathcal{T}}_2}).$$

$${\mathbf{\eta }}_k^{{\mathcal{T}}_1}\perp \perp \text{nolinebreak}{\mathbf{\eta }}_k^{{\mathcal{T}}_2}|\text{nolinebreak}{\mathbf{Z}}^{{\mathcal{F}}_1},{\mathbf{Z}}^{{\mathcal{F}}_2},{\mathbf{\eta }}_k^{\mathrm{rest}},\mathrm{eelse},$$

$$\mathbf{Z}|\text{nolinebreak}\mathbf{\eta },{\mathbf{\theta }}_{\epsilon }$$

$${\mathbf{Z}}^{{\mathcal{F}}_1}$$

Full text

**Keywords: ** Graph Colouring, Markov Chain Monte Carlo, Parallel Sampling, Spatial Statistics

1 Introduction

Large spatial/spatio-temporal data sets are now centre-stage in several of the environmental sciences such as meteorology and glaciology. Two popular tools available to the spatial statistician to deal with such data are the hierarchical model and the closely-related notion of conditional independence (Cressie and Wikle, 2011, Section 2.1.5). In a two layer, linear, Gaussian, data-process model, the widely adopted assumption that data are conditionally independent, given a low-dimensional underlying process, is sufficient for developing inferential algorithms that scale linearly with the dimension of the data (and quadratically with the dimension of the process). Several methods capitalise on this approach for the spatial or spatio-temporal analysis of big data; these include fixed rank kriging (Cressie and Johannesson, 2008), predictive processes (Banerjee et al., 2008), and a suite of approaches based on Gaussian Markov random field (GMRF) approximations to geostatistical models (e.g., Rue and Tjelmeland, 2002; Lindgren et al., 2011; Nychka et al., 2015). For spatio-temporal variants see, for example, Sahu and Mardia (2005); Dewar et al. (2009); Cressie and Wikle (2011) and Wikle et al. (2019).

In a two-layer data-process model, dimensionality reduction of the process needs to be done with care, as the space of the latent functions that can be reproduced is no larger than the span of the set of basis functions used for dimensionality reduction. If this space does not contain the true signal, then any residual observed signal variability will be, at best, attributed to other components in the hierarchical model, typically measurement error or fine-scale process variability. At worst, there is no such component and inferences are over-confident (e.g., Zammit-Mangion and Cressie, 2020). As such, low-dimensional models should only be used with big data if the analyst is confident that the true (hidden) process can be adequately reconstructed using the chosen basis-function representation.

As a result, model approximations have been developed that ensure that several basis functions can be used. Two popular ones include one based on multi-scale GMRFs (Nychka et al., 2015), and one based on a multi-resolution approximation to the process covariance function (Katzfuss, 2017). Unfortunately, models based on multi-scale GMRFs are currently hindered by the computational bottleneck of the required sparse Cholesky factorisations, while neither approach is well-suited to approximate highly nonstationary multi-scale processes. Rather, their motivation, and the associated inference methods that have been designed for them, are built on the premise that the underlying covariance function of the process is relatively simple; the approximations are made to be able to fit the model and predict with it when using large data sets. These methods work very well when this is indeed the case. The next challenge is therefore to use models and appropriate inference methods for when we have big data and for when the underlying process is multi-scale and highly nonstationary, as is the case in many applications.

The model we propose in Section 2 models the process $Y(\cdot)$ as a sum of several processes at various scales, where the degree of nonstationarity increases with the scale of the process. We approximate these processes using GMRFs so that the conditional dependence structure of the latent variables can be exploited for local processing and the use of a parallel Gibbs sampler, as discussed in Section 3. In Section 4 we give some practical guidelines on how to construct the model, and in Section 5 we demonstrate the two-scale variant on a data set containing on the order of one million observations using around one million basis functions, with the covariances constructed though several hundreds of parameters to allow for local nonstationarity. Section 6 concludes.

2 Multi-scale nonstationary spatial processes

In this section we detail the model through its hierarchical structure, adopting the terminology of Berliner (1996). The top layer in the hierarchy is the data model (Section 2.1), the middle layer is the process model (Section 2.2), and the bottom layer is the parameter model (Section 2.3). We often reference subsets of vectors and matrices through superscripts. Specifically, for a vector $\mathbf{X}$, $\mathbf{X}^{I_{1}}$ denotes the sub-vector of $\mathbf{X}$ constructed from the elements with indices in the set $I_{1}$. Sometimes we use subscripts to reference a specific sub-vector wherever using superscript notation would be tedious; for example $\mathbf{X}_{k}$ is the $k$th component, or $k$th sub-vector, of $\mathbf{X}$, where the component’s definition is context-specific. These ways of referencing can be used together; for example $\mathbf{X}_{k}^{\{i\}}$ refers to the $i$th element of the $k$th component of $\mathbf{X}$. When conditioning we sometimes use the symbol “$\mathrm{rest}$” to denote the remainder of a sub-vector, and its precise definition is to be taken from the context. Consider, for example, the $n$-dimensional vector $\mathbf{X}$. If we are considering the distribution of $\mathbf{X}^{{\mathcal{S}}}\mid\mathbf{X}^{\mathrm{rest}}$, where ${\mathcal{S}}\subset\{1,\dots,n\}$, then $\mathbf{X}^{\mathrm{rest}}\equiv\mathbf{X}^{\{1,\dots,n\}\backslash{\mathcal{S}}}$.

2.1 Data model

Assume one has measurements $\mathbf{Z}\equiv(Z_{1},\dots,Z_{m})^{\top}$ of some underlying zero-mean spatial process $Y(\cdot)$ on $D\subset\mathbb{R}^{2}$ or $D\subseteq\mathbb{S}^{2}$, where $\mathbb{S}^{2}$ is the surface of the sphere. Our data model is given by $Z_{l}=Y(\textbf{s}_{l})+\varepsilon_{l}$, for $l=1,\dots,m$, where $\varepsilon_{l}$ is Gaussian measurement-error that is independent of $Y(\cdot)$ and that has variance $\sigma^{2}_{\varepsilon}(\textbf{s}_{l})$, and $\{\textbf{s}_{1},\dots,\textbf{s}_{m}\}\subset D$ is a set of point-referenced measurement locations. In this article we assume that $\{\varepsilon_{l}\}$ are uncorrelated conditional on $\sigma^{2}_{\varepsilon}(\cdot)$, and that $m$ is large; in our case study in Section 5, $m$ is on the order of one million.

2.2 Process model

We define $Y(\cdot)$ to be the sum of $(K+1)$ independent Gaussian processes (GPs), $Y(\cdot)=Y_{0}(\cdot)+\sum_{k=1}^{K}Y_{k}(\cdot),$ where $Y_{0}(\cdot)$ is stationary with large spatial scale (i.e., a large range), and $Y_{k}(\cdot),k\geq 1,$ is nonstationary with spatial scale decreasing with increasing $k$. In our approach $K\geq 1$, and in the case study of Section 5, $K=1$.

Without loss of generality we assume that $Y_{k}(\cdot),k=0,\dots,K,$ have zero expectation. We model the GPs through the stochastic partial differential equations

[TABLE]

and, for $k=1,\dots,K$,

[TABLE]

where $W_{k}(\cdot),k=0,\dots,K,$ are Gaussian spatial white noise processes, $\kappa_{k}$ are spatial scale parameters, $\tau_{k}$ control the process variance, and $\operatorname{\mathrm{\Delta}}$ is the Laplacian. When $D\subset\mathbb{R}^{2}$, the solution to (1) is a stationary GP with Matérn covariance function with smoothness parameter $\nu_{0}=\alpha_{0}-1$, while that to (2) is a nonstationary process which, under some regularity assumptions on $\kappa_{k}(\cdot)$ and $\tau_{k}(\cdot)$, exhibits local Matérn behaviour with smoothness parameter $\nu_{k}=\alpha_{k}-1$ (Lindgren et al., 2011). Throughout, we assume that $\nu_{k},k=0,\dots,K,$ are known.

As in Lindgren et al. (2011) and Lindgren and Rue (2015), we project $Y_{k}(\cdot),k=0,1,\dots,K,$ onto a finite-dimensional basis through $Y_{k}(\cdot)=\textbf{a}_{k}(\cdot)^{\top}{\bm{\eta}}_{k}$, where $\textbf{a}_{k}(\cdot)$ is a vector of $r_{\eta_{k}}$ tent basis functions on a triangulation of $D$, $T_{\eta_{k}}$ say (Scherer, 2017, Section 12.6.2), and ${\bm{\eta}}_{k}\sim\mathrm{Gau}(\mathbf{0},\textbf{Q}_{k}^{-1})$, where $\textbf{Q}_{k}$ is a sparse precision matrix. (Throughout, we use $T_{\beta}$ to denote the triangulation associated with the random quantity $\bm{\beta}$). The triangulations $T_{\eta_{k}},k=0,\dots,K,$ become finer with increasing $k$ since $Y_{k}(\cdot)$ has spatial scale that decreases with $k$. Since each basis function corresponds to a vertex in the triangulation, the coefficients ${\bm{\eta}}_{k},k=0,\dots,K,$ which are independent across $k$, can be associated with the vertices of the basis functions they weight; this property will facilitate graph partitioning and colouring in Section 3.

After projecting onto our finite-dimensional subspace, we can write the measurements as

[TABLE]

where $\textbf{A}_{k}\equiv(\textbf{a}_{k}(\textbf{s}_{l}):l=1,\dots,m)^{\top}$, $\bm{\varepsilon}\sim\mathrm{Gau}(\mathbf{0},\textbf{Q}_{\varepsilon}^{-1}),$ and $\textbf{Q}_{\varepsilon}$ is diagonal with measurement-error precisions as non-zero elements. Collect all the unknown model parameters used to construct $\textbf{Q}_{\varepsilon}$ and $\textbf{Q}\equiv\textrm{bdiag}(\textbf{Q}_{0},\dots,\textbf{Q}_{K})$ into the vector ${\bm{\theta}}$. Let $\textbf{A}\equiv(\textbf{A}_{0},\dots,\textbf{A}_{K})$ and ${\bm{\eta}}\equiv({\bm{\eta}}_{0}^{\top},\dots,{\bm{\eta}}_{K}^{\top})^{\top}$. For this model, which is structurally identical to the LatticeKrig model (Nychka et al., 2015), ${\bm{\eta}}\mid\mathbf{Z},{\bm{\theta}}$ also has a sparse precision matrix since $\mathbf{Z}$ is point referenced. However, both $\mathbf{Z}$ and ${\bm{\eta}}$ can contain millions of elements so that storing, let alone computing with, A and this precision matrix, can become impractical. On the other hand GMRFs have the desirable property that they admit parallel-data and parallel-model conditional structures that can be used to facilitate computation. More details on how we exploit this property are given in Section 3.

2.3 Parameter model

As in Lindgren and Rue (2015) we do not consider the natural parameters $\tau_{0},\kappa_{0}$, and $\tau_{k}(\cdot),\kappa_{k}(\cdot),k=1,\dots,K,$ directly. Rather, we construct a parameter model for the process standard deviations $\sigma_{0}$ and $\{\sigma_{k}(\cdot)\}$, and the nominal ranges $\rho_{0}$ and $\{\rho_{k}(\cdot)\}$. In our setting, for $k=1,\dots,K$ (and similarly for $k=0$), these are related to the natural parameters through

[TABLE]

where $\operatorname{\mathrm{\Gamma}}(\cdot)$ is the gamma function.

We decompose the spatially-varying parameter $\log\sigma_{\varepsilon}(\cdot)$ as a weighted sum of $r_{\varepsilon}$ tent basis functions $\textbf{b}_{\varepsilon}(\cdot)$ constructed on a triangulation $T_{\varepsilon}$ on $D$. We also decompose $\log\sigma_{k}(\cdot)$ and $\log\rho_{k}(\cdot),k=1,\dots,K,$ as weighted sums of $r_{\theta_{k}}\equiv r_{\sigma_{k}}=r_{\rho_{k}}$ tent basis functions, $\textbf{b}_{k}(\cdot)$, constructed on a triangulation $T_{\theta_{k}}\equiv T_{\sigma_{k}}=T_{\rho_{k}}$ on $D$. Note that here, and below, we use the subscript $\theta_{k}$ to identify quantities such as basis functions and triangulations that are common to both $\sigma_{k}(\cdot)$ and $\rho_{k}(\cdot)$. Since the parameter fields at scale $k$ need to vary slowly relative to the process at that scale (Lindgren et al., 2011), the triangulation $T_{\theta_{k}}$ can be much coarser than $T_{\eta_{k}}$. These decompositions yield the following parameter models,

[TABLE]

for $k=1,\dots,K$, remembering that, for $k=0$, $\sigma_{0}$ and $\rho_{0}$ are spatially-invariant.

We assume that the elements of the parameter vectors ${\bm{\theta}}_{X},X\in\{\varepsilon,\sigma_{1},\rho_{1},\dots,\sigma_{K},\rho_{K}\}$, are independent and distributed according to Gaussian distributions. Specifically, we let

[TABLE]

where $\omega_{X}=\mathrm{E}({\bm{\theta}}_{X}^{\{i\}})$, $\lambda_{X}=\mathrm{var}({\bm{\theta}}_{X}^{\{i\}})$, and $r_{X}$ is the dimension of ${\bm{\theta}}_{X}$. Although we assume independence, a GMRF model for the parameters might also be used (see Monterrubio-Gómez et al., 2019), although this would introduce additional complexity to our updating algorithm of Section 3. The hyperparameters $\omega_{\rho_{k}}$ and $\lambda_{\rho_{k}}$ need to reflect a prior judgement that the process range decreases with increasing $k$; we provide some guidelines in Section 4 and an example of how this can be done in practice with a two-scale process in Section 5. Appropriate hyperparameters for constructing the prior distributions over ${\bm{\theta}}_{\varepsilon}$ and ${\bm{\theta}}_{\sigma_{k}}$ can usually be deduced from domain knowledge and exploratory data analysis (see Section 5.2).

We complete our parameter model by equipping the parameters appearing in the coarse component, $\rho_{0}$ and $\sigma_{0}$, with the prior lognormal distributions,

[TABLE]

where reasonable values of $\omega_{\sigma_{0}}$ and $\lambda_{\sigma_{0}}$ can generally be deduced from exploratory analysis of the data, while $\omega_{\rho_{0}}$ and $\lambda_{\rho_{0}}$ can generally be fixed using prior application knowledge (see Sections 4 and 5). As with $k>0$, we assume that $\sigma_{0}$ and $\rho_{0}$ are independent.

To summarise, the full set of parameters and the full conditional dependence structure of our model is given by the following directed acyclic graph (DAG):

[TABLE]

where the components of ${\bm{\theta}}_{\sigma_{k}}$ and $\bm{\theta}_{\rho_{k}}$ are mutually independent for $k=1,\dots,K$. The coefficients ${\bm{\eta}}_{0}$ construct a stationary process with spatially-invariant standard deviation and range, $\sigma_{0}$ and $\rho_{0}$, respectively. On the other hand, ${\bm{\eta}}_{1},\dots,{\bm{\eta}}_{K},$ construct nonstationary processes with spatially-varying standard deviations and ranges. For these nonstationary processes, the standard deviations and ranges are modelled through the coefficients ${\bm{\theta}}_{\sigma_{k}}$ and ${\bm{\theta}}_{\rho_{k}},~{}k=1,\dots,K$, where the individual elements are independent and their marginal distributions are pre-specified. The measurement errors $\bm{\varepsilon}$ are assumed to be spatially uncorrelated and are thus modelled using just the basis-function coefficients, ${\bm{\theta}}_{\varepsilon}$.

3 Inference

Our aim is to sample from the posterior distribution $\operatorname{\mathrm{pr}}({\bm{\eta}}_{0},\dots,{\bm{\eta}}_{K},{\bm{\theta}}_{\varepsilon},{\bm{\theta}}_{0},\dots,{\bm{\theta}}_{K}\mid\mathbf{Z})$, where ${\bm{\theta}}_{0}\equiv(\sigma_{0},\rho_{0})^{\top}$, and ${\bm{\theta}}_{k}\equiv({\bm{\theta}}_{\sigma_{k}}^{\top},{\bm{\theta}}_{\rho_{k}}^{\top})^{\top},k=1,\dots,K$. The value of the model is partly that it provides a flexible representation of spatial nonstationarity, and partly that its posterior distribution can be sampled from in a serial fashion using many Gibbs steps. Additionally, the Gibbs steps for ${\bm{\eta}}_{k}$ and ${\bm{\theta}}_{k},k=1,\dots,K,$ can be split and parallelised if they are too difficult to sample directly. We hence have a parallel MCMC scheme that can scale well with both data size and model complexity.

In Section 3.1 we discuss Gibbs sampler updates for the process and parameter coefficients in the first scale; in Section 3.2 the updates for the process and parameter coefficients in the other scales; in Section 3.3 a re-updating strategy for the process coefficients for $k>0$; in Section 3.4 the updates for the measurement-error variance parameters; and in Section 3.5 we summarise the sampler.

3.1 Updating the process and parameter coefficients for scale $k=0$

The process $Y_{0}(\cdot)$ captures large-scale effects, and hence its triangulation $T_{\eta_{0}}$ can be coarse relative to $T_{\eta_{k}}$ for $k=1,\dots,K$. Since ${\bm{\theta}}_{0}\equiv(\sigma_{0},\rho_{0})^{\top}$ and ${\bm{\eta}}_{0}$ can be expected to be highly correlated a posteriori (Knorr-Held and Rue, 2002) we sample them jointly in the Gibbs sampler, which substantially improves the mixing of the Markov chain. Specifically, we sample from the full conditional distribution

[TABLE]

by first sampling from $\operatorname{\mathrm{pr}}({\bm{\theta}}_{0}\mid\mathrm{eelse}~{}\mathrm{excl.}~{}{\bm{\eta}}_{0})$ and then from $\operatorname{\mathrm{pr}}({\bm{\eta}}_{0}\mid{\bm{\theta}}_{0},\mathrm{eelse})$, where we use ‘$\mathrm{eelse}$’ as shorthand for ‘everything else.’ Here, and in the rest of the article, we always sample using the most recently sampled values of the other variables.

The conditional distribution of ${\bm{\theta}}_{0}$ in (6) is a partially collapsed distribution, where ${\bm{\eta}}_{0}$ is integrated out:

[TABLE]

The integrand in (7) is Gaussian, and therefore the integral can be evaluated analytically to give

[TABLE]

where $\widetilde{\textbf{Q}}_{0}=\textbf{A}_{0}^{\top}\textbf{Q}_{\varepsilon}\textbf{A}_{0}+\textbf{Q}_{0},\widetilde{\bm{\mu}}_{0}=\widetilde{\textbf{Q}}_{0}^{-1}\textbf{A}_{0}^{\top}\textbf{Q}_{\varepsilon}\widetilde{\textbf{Z}},$ and $\widetilde{\mathbf{Z}}=\mathbf{Z}-\sum_{k>0}\textbf{A}_{k}{\bm{\eta}}_{k}$. The triangulation $T_{\eta_{0}}$ must be made sufficiently coarse that factorising $\textbf{Q}_{0}$, as well as algebraic operations of the form $\widetilde{\textbf{Q}}_{0}^{-1}\textbf{X}=\mathbf{Y}$, can be done on a single computing node. We provide more detailed guidelines on this in Section 4. All other operations can be computed in a distributed fashion from sums and products of smaller vectors and matrices corresponding to chunks of $\mathbf{Z},\{\textbf{A}_{k}\}$, and $\{{\bm{\eta}}_{k}\}$. There is therefore no theoretical limit on the data size and number of scales $K$ that will preclude sampling from this distribution in a reasonable time-frame given sufficient parallel computing resources. This is also true for all the sampling operations we outline below. The conditional distribution (8) does not have a recognisable distribution in ${\bm{\theta}}_{0}$, and so this update uses a Metropolis–Hastings step. The computational complexity of this step is dominated by the sparse Cholesky factorisation of $\widetilde{\textbf{Q}}_{0}$, which is approximately $O(r_{\eta_{0}}^{3/2})$ flops.

After updating ${\bm{\theta}}_{0}$ we update ${\bm{\eta}}_{0}$. The full conditional distribution of ${\bm{\eta}}_{0}$ is given by

[TABLE]

where $\widetilde{\bm{\mu}}_{0}$ and $\widetilde{\textbf{Q}}_{0}$ have already been computed for (8). A sample of (9) is therefore obtained ‘for free’ after the update of ${\bm{\theta}}_{0}$.

3.2 Updating the process and parameter coefficients for scales $k=1,\dots,K$

The update of ${\bm{\theta}}_{k}\equiv({\bm{\theta}}_{\sigma_{k}}^{\top},{\bm{\theta}}_{\rho_{k}}^{\top})^{\top}$ is more tricky than that of ${\bm{\theta}}_{0}$: The full conditional distribution of ${\bm{\theta}}_{k}$ has the same structure as that of ${\bm{\theta}}_{0}$ given in (7), but since for $k>0$ the triangulation $T_{\eta_{k}}$ is fine and $r_{\eta_{k}}$ is large, the elements of the integral cannot be evaluated in memory, and so a sequential strategy is necessary.

Let ${\mathcal{I}}$ index one element of ${\bm{\theta}}_{\sigma_{k}}$ (equivalently ${\bm{\theta}}_{\rho_{k}}$) and, with a slight abuse of notation, let ${\bm{\theta}}_{k}^{\mathcal{I}}$ denote the pair of elements $({\bm{\theta}}_{\sigma_{k}}^{\mathcal{I}},{\bm{\theta}}_{\rho_{k}}^{\mathcal{I}})^{\top}$. Let the set ${\mathcal{T}}$ be the effective process footprint of ${\bm{\theta}}_{k}^{\mathcal{I}}$, defined as the set of indices of ${\bm{\eta}}_{k}$ for which the (prior) conditional distribution of ${\bm{\eta}}_{k}^{\mathrm{rest}}$ excluding ${\bm{\eta}}_{k}^{\mathcal{T}}$, $\operatorname{\mathrm{pr}}({\bm{\eta}}_{k}^{\mathrm{rest}}\mathbin{|}\nolinebreak{\bm{\theta}}_{k}^{\mathcal{I}},{\bm{\theta}}_{k}^{\mathrm{rest}})$, is approximately independent of ${\bm{\theta}}_{k}^{\mathcal{I}}$ (see (11) below). It might be helpful to picture ${\bm{\eta}}_{k}^{\mathcal{T}}$ as a ‘halo’ around $\operatorname*{supp}(\textbf{b}_{k}^{{\mathcal{I}}})$, where $\operatorname*{supp}(\cdot)$ denotes a function’s support; see Fig. 1. Further, let ${\mathcal{F}}$ be the data footprint of ${\bm{\eta}}_{k}^{{\mathcal{T}}}$, that is, the set of indices of the data points that lie in the subset of the domain given by $\cup_{i\in{\mathcal{T}}}\operatorname*{supp}(\textbf{a}_{k}^{\{i\}})$. Under the partitioning of the coefficients $\{{\bm{\theta}}_{k}^{{\mathcal{I}}},{\bm{\theta}}_{k}^{\mathrm{rest}}\}$, $\{{\bm{\eta}}_{k}^{\mathcal{T}},{\bm{\eta}}_{k}^{\mathrm{rest}}\}$, and the data $\{\mathbf{Z}^{{\mathcal{F}}},\mathbf{Z}^{\mathrm{rest}}\}$, the relevant part of the full DAG is

[TABLE]

where we have used a dashed line to denote a conditional dependence relationship that is weak, and one that can potentially be ignored for practical purposes.

We jointly update $\{{\bm{\eta}}_{k}^{\mathcal{T}},{\bm{\theta}}_{k}^{\mathcal{I}}\}$ via the full conditional distribution

[TABLE]

by first sampling from $\operatorname{\mathrm{pr}}({\bm{\theta}}_{k}^{\mathcal{I}}\mid\mathrm{eelse}\textrm{~{}excl.~{}}{\bm{\eta}}_{k}^{\mathcal{T}})$ and subsequently from $\operatorname{\mathrm{pr}}({\bm{\eta}}_{k}^{\mathcal{T}}\mid{\bm{\theta}}_{k}^{\mathcal{I}},\mathrm{eelse})$. As with ${\bm{\theta}}_{0}$ we have that $\operatorname{\mathrm{pr}}({\bm{\theta}}_{k}^{\mathcal{I}}\mid\mathrm{eelse}\textrm{~{}excl.~{}}{\bm{\eta}}_{k}^{\mathcal{T}})$ is a partially collapsed distribution given by

[TABLE]

Since ${\bm{\eta}}_{k}$ is a Gaussian Markov random field it is computationally cheap to evaluate the conditional distribution $\operatorname{\mathrm{pr}}({\bm{\eta}}_{k}^{\mathcal{T}}\mid{\bm{\eta}}_{k}^{\mathrm{rest}},{\bm{\theta}}_{k}^{\mathcal{I}},{\bm{\theta}}_{k}^{\mathrm{rest}})$. However, it is computationally expensive to evaluate the ‘marginal’ distribution $\operatorname{\mathrm{pr}}({\bm{\eta}}_{k}^{\mathrm{rest}}\mid{\bm{\theta}}_{k}^{\mathcal{I}},{\bm{\theta}}_{k}^{\mathrm{rest}})$, as is required for this update, since it has a (large) precision matrix that is not sparse in general. We therefore resort to the approximation implied by the dashed line of the DAG in (10),

[TABLE]

where we have replaced ${\bm{\theta}}_{k}^{\mathcal{I}}$ by some value $\widetilde{\bm{\theta}}_{k}^{\mathcal{I}}$, which may depend on the value of ${\bm{\theta}}_{k}^{\mathrm{rest}}$. This approximation is motivated by the fact that, notionally, ${\bm{\theta}}_{k}^{\mathcal{I}}$ controls the variance parameters $\sigma_{k}(\cdot)$ and $\rho_{k}(\cdot)$ locally, and can therefore be expected to have diminishing effect on the probability of ${\bm{\eta}}_{k}^{\mathrm{rest}}$ by the marginalisation property of the Gaussian distribution. An attractive feature of this approximation is that it can be improved as much as needed, at the expense of increased computational cost, by increasing the size of the effective process footprint ${\mathcal{T}}$. The chosen effective process footprint should be one for which the approximation (11) is acceptable in practice. In our application, setting ${\mathcal{T}}$ to be the indices of the basis functions $\textbf{a}_{k}(\cdot)$ that have their and their neighbours’ maxima inside $\operatorname*{supp}(\textbf{b}_{k}^{{\mathcal{I}}})$ sufficed, but more conservative choices could be considered if needed; see Sect. 4.3.

Our approximation yields

[TABLE]

the log of which is identical in structure to (8) with an added term due to the presence of ${\bm{\eta}}_{k}^{\mathrm{rest}}$ in the conditional distribution of ${\bm{\eta}}_{k}^{\mathcal{T}}$ (Rue and Held, 2005, Sect. 2.2.5). This extra term does not alter the computational or memory complexity of the operations. In particular, the computational complexity is approximately $O(|{\mathcal{T}}^{3/2}|)$ flops, and thus sampling will be feasible as long as $|{\mathcal{T}}|$ has the same order of magnitude as $r_{\eta_{0}}$, which we have set to be sufficiently small to make it possible for these computations to be done on a single computing node.

As with ${\bm{\theta}}_{0}$, this distribution has no recognisable form and requires a Metropolis–Hastings step. Further, the computations for this conditional distribution can once again benefit from a distributed data framework. However, since there are $r_{\theta_{k}}$ of them (one for each parameter basis function), for the inference framework to be scaleable, the total computational time required by these updates must not depend on $r_{\theta_{k}}$. This is indeed the case under our approximation as we now show.

Take two sets of parameters ${\bm{\theta}}_{k}^{{\mathcal{I}}_{1}}$ and ${\bm{\theta}}_{k}^{{\mathcal{I}}_{2}}$ with associated effective process footprints ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$ for which ${\bm{\eta}}_{k}^{{\mathcal{T}}_{1}}\mathrel{\perp\!\!\perp}\nolinebreak{\bm{\eta}}_{k}^{{\mathcal{T}}_{2}}\mathbin{|}\nolinebreak{\bm{\eta}}_{k}^{\mathrm{rest}},{\bm{\theta}}_{k},$ and for which the associated data footprints ${\mathcal{F}}_{1}\cap{\mathcal{F}}_{2}=\emptyset$ (this condition is trivially satisfied in our context where $\mathbf{Z}$ are point referenced and the $\{\textbf{a}_{k}(\cdot)\}$ are tent basis functions). Partition ${\bm{\theta}}_{k}$ as $\{{\bm{\theta}}_{k}^{{\mathcal{I}}_{1}},{\bm{\theta}}_{k}^{{\mathcal{I}}_{2}},{\bm{\theta}}_{k}^{\mathrm{rest}}\}$, ${\bm{\eta}}_{k}$ as $\{{\bm{\eta}}_{k}^{{\mathcal{T}}_{1}},{\bm{\eta}}_{k}^{{\mathcal{T}}_{2}},{\bm{\eta}}_{k}^{\mathrm{rest}}\}$, and $\mathbf{Z}$ as $\{\mathbf{Z}^{{\mathcal{F}}_{1}},\mathbf{Z}^{{\mathcal{F}}_{2}},\mathbf{Z}^{\mathrm{rest}}\}$. Then, the relevant part of the DAG where we now omit the weak conditional dependencies for clarity, is

[TABLE]

From the moralised version (Lauritzen, 1996, Chapter 2) of (16), we immediately see that the coefficients $\{{\bm{\eta}}_{k}^{\mathrm{rest}},{\bm{\theta}}_{k}^{\mathrm{rest}}\}$ separate the two groups $\{{\bm{\eta}}_{k}^{{\mathcal{T}}_{1}},{\bm{\theta}}_{k}^{{\mathcal{I}}_{1}}\}$ and $\{{\bm{\eta}}_{k}^{{\mathcal{T}}_{2}},{\bm{\theta}}_{k}^{{\mathcal{I}}_{2}}\}$ so that

[TABLE]

Therefore, the update operations for ${\bm{\theta}}_{k}^{{\mathcal{I}}_{1}}$ and ${\bm{\theta}}_{k}^{{\mathcal{I}}_{2}}$ can be dispatched to different cores in a multicore computing architecture, and the updates can be done in parallel. We can identify the sets of parameters that can be updated in parallel by colouring the full conditional dependency graph of ${\bm{\theta}}_{k}$ under our chosen set of effective footprints, such that no conditionally dependent elements of ${\bm{\theta}}_{k}$ have the same colour. In a minimal setting, where we choose the effective process footprints as the subset of ${\bm{\eta}}_{k}$ for which the maxima of the associated basis functions and those of its Markov blanket lie in the support of $\textbf{b}_{k}^{\mathcal{I}}$, the conditional dependency graph is the edge-preserving bijection of the spatial graph $T_{\theta_{k}}$. With this choice of footprints we can then make use of the four-colour theorem (Gonthier, 2008) so that a sample of ${\bm{\theta}}_{k}$ can be done in exactly four steps, irrespective of the scale or data size. More conservative (i.e., larger) footprints will result in the need for more colours, and hence a lower degree of parallelisation, but still the number of colours will be independent of scale and data size, rendering this model and inferential technique scaleable. In our application we made use of the backtracking algorithm for colouring the graph (Bender and Wilf, 1985).

Consider now the joint update of the quantities $\{{\bm{\eta}}_{k}^{\mathcal{T}},{\bm{\theta}}_{k}^{\mathcal{I}}\}$. As in the joint update of $\{{\bm{\eta}}_{0},{\bm{\theta}}_{0}\}$, a sample from the full conditional distribution of ${\bm{\eta}}_{k}^{\mathcal{T}}$ is obtained ‘for free’ since it is Gaussian with a mean vector and a precision matrix that are also used when sampling ${\bm{\theta}}_{k}^{\mathcal{I}}$. One might be tempted to just skip this step, obtain a full sample of ${\bm{\theta}}_{k}$, and only then obtain a full sample of ${\bm{\eta}}_{k}$ using, for example, the re-updating strategy in Sect. 3.3. However, doing so will result in a collapsed Gibbs sampler that does not target the correct stationary distribution (Van Dyk and Park, 2008); see Appendix A for more details on why this is the case. It is therefore important that ${\bm{\eta}}_{k}^{\mathcal{T}}$ is sampled concurrently with ${\bm{\theta}}_{k}^{\mathcal{I}}$ although, as we discuss next, a re-updating of ${\bm{\eta}}_{k}$ is necessary in practice to improve convergence of the Markov chain.

3.3 Re-updating the process coefficients for scales $k=1,\dots,K$

When a block of ${\bm{\eta}}_{k}$ gets updated conditional on everything else it is pinned at its boundary, because ${\bm{\eta}}_{k}$ is locally smooth a priori and, in regions where data density is low, a posteriori as well. This results in slow mixing of the Markov chain, which was problematic in earlier versions of our implementation (see Appendix B for an illustration on a simple model). We now address this issue by re-sampling each ${\bm{\eta}}_{k}$ in a separate Gibbs step.

For each ${\bm{\eta}}_{k}$ we tile the domain $D$ into tiles so that each element of ${\bm{\eta}}_{k}$ is then associated with exactly one tile, specifically the tile in which the basis function it weights is a maximum. Each tile is no larger than what can be processed on a single node, and so will be associated with about $r_{\eta_{0}}$ basis functions, as per the update of ${\bm{\eta}}_{0}$. The tiles must also be large enough such that the elements of ${\bm{\eta}}_{k}$ associated with any two non-contiguous tiles (i.e., tiles that do not share a common boundary) are conditionally independent given the elements of ${\bm{\eta}}_{k}$ associated with the other tiles. These tiles and their neighbours can be used to establish a supergraph made up of blocks of ${\bm{\eta}}_{k}$, where each block (corresponding to one tile) is conditionally independent of the rest given its neighbours. This supergraph is then coloured such that no two neighbouring blocks have the same colour; for point-referenced data this will require at most four colours. The colouring also corresponds to one on the tiles where no two tiles which share a common border have the same colour.

As with the joint update of $\{{\bm{\eta}}_{k},{\bm{\theta}}_{k}\}$, this colouring allows us to develop a parallel sampler for ${\bm{\eta}}_{k}$. Specifically, consider two blocks of ${\bm{\eta}}_{k}$ that are associated with two tiles, with indices given by ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{2}$, respectively. When these two blocks are of the same colour, ${\bm{\eta}}_{k}^{{\mathcal{T}}_{1}}\mathrel{\perp\!\!\perp}\nolinebreak{\bm{\eta}}_{k}^{{\mathcal{T}}_{2}}\mathbin{|}\nolinebreak{\bm{\eta}}_{k}^{\mathrm{rest}}$ according to the structural zeros of $\textbf{Q}_{k}$. Let ${\mathcal{F}}_{1}$ be the data footprint of ${\bm{\eta}}_{k}^{{\mathcal{T}}_{1}}$. Similarly, let ${\mathcal{F}}_{2}$ be the data footprint of ${\bm{\eta}}_{k}^{{\mathcal{T}}_{2}}$. For point-referenced measurements it is straightforward to see that ${\mathcal{F}}_{1}\cap{\mathcal{F}}_{2}=\emptyset$ so that we can partition $\mathbf{Z}$ as $\{\mathbf{Z}^{{\mathcal{F}}_{1}},\mathbf{Z}^{{\mathcal{F}}_{2}},\mathbf{Z}^{\mathrm{rest}}\}$. The relevant part of the DAG for ${\bm{\eta}}_{k}$ and $\mathbf{Z}$ is therefore simply the first two lines of (16), from which we obtain

[TABLE]

and therefore ${\bm{\eta}}_{k}^{{\mathcal{T}}_{1}}$ and ${\bm{\eta}}_{k}^{{\mathcal{T}}_{2}}$ can be updated in parallel (Wilkinson, 2006). The log of the full conditional distributions are identical in form to (8) with the same extra term mentioned in Sect. 3.2, which is simple to compute. Therefore, each of these updates again requires only a small chunk of data, matrices, and subsets of samples of the other process coefficients, so that each of these updates can also be done via a distributed-data architecture. The computational complexity of the operations required to sample from ${\bm{\eta}}_{k}^{{\mathcal{T}}_{1}}$ and ${\bm{\eta}}_{k}^{{\mathcal{T}}_{2}}$ are approximately $O(|{\mathcal{T}}_{1}^{3/2}|)$ and $O(|{\mathcal{T}}_{2}^{3/2}|)$, respectively.

Even though this sampling strategy re-updates ${\bm{\eta}}_{k}$, this coloured-tile scheme has got its own boundaries close to which the process coefficients will experience poor mixing. We address this by using three tilings for each $k$, one after another in the Gibbs sampler. Each tiling is shifted by about a third of a tile, relative to the one before. This ensures that no element of ${\bm{\eta}}_{k}$ is always near a boundary when Gibbs sampling. We show three such tilings on $\mathbb{S}^{2}$ in Fig. 2. The resulting Gibbs sampler is an instance of a blocked sampler where the blocks are not disjoint; see, for example Jensen et al. (1995). Appendix B illustrates the benefit of an alternating tiling scheme on a simple model.

3.4 Updating the measurement-error variance parameters

Consider now the parameters ${\bm{\theta}}_{\varepsilon}$. The submodel for the Gibbs update of these parameters is made up of

[TABLE]

and the prior distribution (3) with $X=\varepsilon$, where $\textbf{Q}_{\varepsilon}^{-1}$ is diagonal with non-zero elements $[\exp(\textbf{b}_{\varepsilon}(\textbf{s}_{l})^{\top}{\bm{\theta}}_{\varepsilon})]^{2},l=1,\dots,m$. As with ${\bm{\eta}}_{k}$ and ${\bm{\theta}}_{k}$, elements of ${\bm{\theta}}_{\varepsilon}$ are spatially-indexed by the maxima of the basis functions they weight. Let ${\mathcal{I}}$ contain the index of one element of ${\bm{\theta}}_{\varepsilon}$. The data footprint ${\mathcal{F}}$ of ${\bm{\theta}}_{\varepsilon}^{\mathcal{I}}$ contains the indices of the data points that lie in $\operatorname*{supp}(\textbf{b}_{\varepsilon}^{{\mathcal{I}}})$. For two elements of ${\bm{\theta}}_{\varepsilon}$ with indices ${\mathcal{I}}_{1}$ and ${\mathcal{I}}_{2}$, respectively, and for which ${\mathcal{F}}_{1}\cap{\mathcal{F}}_{2}=\emptyset$, the relevant part of the conditional independence graph is

[TABLE]

because the components of $\mathbf{Z}$ are independent given ${\bm{\theta}}_{\varepsilon}$ and everything else, and the components of ${\bm{\theta}}_{\varepsilon}$ are IID in the prior.

As with ${\bm{\eta}}_{1}$, the condition ${\mathcal{F}}_{1}\cap{\mathcal{F}}_{2}=\emptyset$ is essential, but this is easily satisfied in our context since $\textbf{b}_{\varepsilon}(\cdot)$ are tent basis functions and the data are point-referenced. In particular, this condition is satisfied whenever ${\bm{\theta}}_{\varepsilon}^{{\mathcal{I}}_{1}}$ and ${\bm{\theta}}_{\varepsilon}^{{\mathcal{I}}_{2}}$ are not neighbours in the graph implied by $T_{\varepsilon}$. In this case, the full conditional dependence graph of ${\bm{\theta}}_{\varepsilon}$ is just an edge-preserving bijection of the spatial graph $T_{\varepsilon}$, with ${\bm{\theta}}_{\varepsilon}$ as the vertices. We therefore also colour it using four colours, and we update parameters of the same colour in parallel. The full conditional distribution,

[TABLE]

does not have a recognisable distribution, and so this update uses a Metropolis–Hastings step. Since the measurement errors are independent, computing the accept-reject ratio is computationally inexpensive.

3.5 Summary of Gibbs sampler

Algorithm 3.5 gives a summary of all the stages in one complete pass through the Gibbs sampler. Each sample is generated conditional on $\mathbf{Z}$ and the most recent samples of the other random variables, again denoted here as ‘$\mathrm{eelse}$’ for ‘everything else.’ When the conditional distribution to sample from is not available in closed form, the update is done via an accept-reject step; in our application in Sect. 5 we use adaptive random walk Metropolis proposals. In Algorithm 3.5 we use the superscript $*$ to denote intermediate quantities: these quantities are discarded in the final output of the Gibbs sampler.

Steps 1–3 in the Gibbs sampler use distributed data, while steps 2 and 3 also use parallel computation to obtain a complete sample in a fixed number of sweeps (usually 4). It is worth re-iterating that, despite the re-updating of ${\bm{\eta}}_{k}$ in step 2(b), updating ${\bm{\eta}}_{k}$ for each colour in step 2(a) is still required in order to ensure the correct stationary distribution is targeted when going through the colours in sequence (Van Dyk and Park, 2008); see Appendix A.

4 Implementation guidelines

The favourable computational properties of the updates in Sect. 3 can only be taken advantage of if sensible triangulations are used for the process and parameter decompositions. Further, the separation of $Y(\cdot)$ into separate scales introduces identifiability issues, akin to those seen in problems of spatial source separation (Nordhausen et al., 2015). Therefore, judicious construction of the triangulations and prior distributions on the nominal ranges (that help to separate the scales in spectral space) are required. Our informal guidelines in this section are based on point-referenced measurements that have relatively high signal-to-noise ratio; these guidelines may need to be slightly adjusted when this is not the case.

4.1 Triangulations

Process:

Each process $Y_{k}(\cdot)$ is associated with a triangulation $T_{\eta_{k}}$. Recall that, as $k$ increases, the process scale decreases, and therefore we seek triangulations that become finer with increasing $k$. Constructing appropriate triangulations on which to establish the multi-scale process is important. In Appendix C we empirically show the possible consequences of an inappropriate discretisation scheme for a two-scale process in one-dimensional space.

A key design criterion when constructing these triangulations is the size of largest matrix that can be efficiently factorised on a single computing node. Using hardware and linear algebra libraries current as of the year 2020, sparse precision matrices of dimension 50000 constructed from a second-order GMRF on a two-dimensional manifold (following appropriate fill-in reducing permutations) can be factorised extremely quickly (in about one second) on a single computing node. Since in our model ${\bm{\eta}}_{0}$ is updated as a whole (see (6)), a regular triangulation $T_{\eta_{0}}$ on $D$ should therefore be constructed such that $r_{\eta_{0}}$ is between $10^{4}$ and $10^{5}$.

At the other end of the spectrum, the finest process triangulation, $T_{\eta_{K}}$, should be fine enough that the size of the data footprint associated with each element of ${\bm{\eta}}_{K}$ is small, and preferably in the single digits. This guideline ensures that our multi-scale process model is able to resolve the highest frequency components that may be extracted from the data. If we do not do this, then our predictions at the finest scale are possibly over-smoothed (i.e., the model is underfitting).

The number of scales $K$ to use is a modelling choice which can be guided by spectral considerations. A good rule of thumb is to ensure that each scale is associated with signals with nominal ranges that approximately span an order of magnitude, and that the triangulations can adequately reproduce signals with these ranges. On a two-dimensional manifold, this rule of thumb implies that $|{\bm{\eta}}_{k}|\approx 100|{\bm{\eta}}_{k-1}|$, for $k=1,\dots,K$. In Sect. 5, we use two scales with $r_{\eta_{0}}$ on the order of $10^{4}$ and $r_{\eta_{1}}$ on the order of $10^{6}$.

Process parameters:

The update of ${\bm{\theta}}_{k}^{\mathcal{I}}$ in (15) requires factorising a square matrix with number of columns equal to the size of its effective process footprint, $|{\mathcal{T}}|$, which thus needs to be on the order of $r_{\eta_{0}}$. Therefore, since $T_{\eta_{k}}$ becomes finer with increasing $k$, this necessarily means that $T_{\theta_{k}}$ also needs to become finer with increasing $k$. Specifically, the number of basis functions in $\textbf{a}_{k}(\cdot)$ that have support intersecting $\textbf{b}_{k}^{\mathcal{I}}(\cdot)$, for each ${\mathcal{I}}$, needs to be on the order of $r_{\eta_{0}}$. This successive refining of $T_{\theta_{k}}$ with $k$ yields an attractive model where the degree of nonstationarity of the process increases with $k$.

Measurement-error standard deviation:

The coarseness of the triangulation $T_{\varepsilon}$ reflects a prior belief on how quickly (spatially) the measurement-error variance changes; this is application-specific.

4.2 Parameter prior distributions

Nominal ranges:

In Sect. 4.1 we associated $T_{\eta_{k}}$ with nominal ranges spanning approximately an order of magnitude, and this choice can be reflected in our prior distribution for the nominal ranges. We suggest that $\omega_{\rho_{0}}$ and $\lambda_{\rho_{0}}$ are set such that a low percentile of $\operatorname{\mathrm{pr}}(\rho_{0})$, the 2.5 percentile say, is no smaller than one to two times the average edge length in $T_{\eta_{0}}$. Similarly, we suggest that $\omega_{\rho_{k}}$ and $\lambda_{\rho_{k}}$ are set such that a low percentile of $\operatorname{\mathrm{pr}}(\exp({\bm{\theta}}_{\rho_{k}}^{\mathcal{I}}))$, for each ${\mathcal{I}}$, is no smaller than one to two times the average edge length in $T_{\eta_{k}}$. This choice is motivated by the Nyquist–Shannon criterion, which is also used for basis function placement by Zammit-Mangion et al. (2012). From a practical point of view, it encapsulates the fact that signal components with finer scales cannot be captured by this decomposition of $Y_{k}(\cdot)$ on $T_{\eta_{k}}$. The upper 2.5 percentile of $\operatorname{\mathrm{pr}}(\rho_{0})$, corresponding to the highest plausible nominal range of the multi-scale process, can usually be set based on application considerations. The upper 2.5 percentile of $\operatorname{\mathrm{pr}}(\exp({\bm{\theta}}_{\rho_{k}}^{\mathcal{I}})),k=1,2,\dots,K,$ for all ${\mathcal{I}}$, can be set to be close to, but higher than, the lower 2.5 percentie of $\operatorname{\mathrm{pr}}(\exp({\bm{\theta}}_{\rho_{k-1}}^{\mathcal{I}}))$. This ensures some spectral overlap across the scales.

Process standard deviations:

The prior distribution $\operatorname{\mathrm{pr}}(\sigma_{0})$ can be uninformative or based on results from an exploratory data analysis. Since in this model source separation is done via the prior distributions of $\rho_{k}$, the hyperparameters $\omega_{\sigma_{k}}$ and $\lambda_{\sigma_{k}}$ can be set such that $\operatorname{\mathrm{pr}}(\exp({\bm{\theta}}_{\sigma_{k}}^{{\mathcal{I}}}))$, for each ${\mathcal{I}}$, is relatively uninformative or identical to $\operatorname{\mathrm{pr}}(\sigma_{0})$.

Measurement-error standard deviation:

The hyperparameters $\omega_{\varepsilon}$ and $\lambda_{\varepsilon}$ encode the plausible range of values of measurement-error standard deviations, and their choice is therefore application-specific.

4.3 Effective process footprint

The effective process footprint of each ${\bm{\theta}}_{k}^{\mathcal{I}}$ is a design choice. In principle, the footprint needs to be one such that the approximation (11) is acceptable. For example, one might use a spatial buffer around $\operatorname*{supp}(\textbf{b}_{k}^{\mathcal{I}})$ of width equal to the upper 97.5 percentile of $\operatorname{\mathrm{pr}}(\exp({\bm{\theta}}_{\rho_{k}}^{\mathcal{I}}))$ to define the effective process footprint of ${\bm{\theta}}_{k}^{\mathcal{I}}$. If this footprint is so large that sampling of ${\bm{\theta}}_{k}^{\mathcal{I}}$ becomes prohibitive, then tighter prior distributions (and hence more scales in order to ensure spectral coverage everywhere) would need to be used if this approximation is of concern. Alternatively, one might simulate a number of Markov chains with different computationally inexpensive effective process footprints in parallel, and verify that the target distribution is indeed (approximately) independent of the chosen footprint size. As discussed in Sect. 3.2, in our application we have taken a minimal approach and used relatively small effective process footprints. This choice does not seem to have had a detrimental effect on the quality of the probabilistic predictions which, as we show next in Sect. 5, are considerably better than those obtained from competing models.

5 Case study: Spatial modelling and prediction of sea-surface temperature from VIIRS

A sophisticated distributed implementation on a multi-node cluster is required for the sampler to be used with data-set sizes on the order of tens or hundreds of millions and with models with many ($K>1$) scales. However, a two-scale model with approximately one million basis functions and one million data points can be implemented on a standard multi-core computing node using a straightforward implementation in R (R Core Team, 2019), whereby the “chunks” required for each parallel sample (implemented via the function mclapply) are loaded from disk when needed. This data size is still considered fairly large in spatial analyses, and we therefore restrict ourselves to this case in this article. An implementation that allows for orders of magnitude more data and basis functions is being developed and will be discussed elsewhere.

We compare inferences from the two-scale variant of the multi-scale spatial process model with those from several other spatial models and approximation methods, specifically a coarse-scale SPDE model; a suite of independent fine-scale SPDE models over spatially contiguous blocks of data; a single-scale model approximated using a full-scale approximation (FSA); a stationary multi-scale process model constructed using a relatively small number of basis functions, and a single-scale process model using the approximations via a nearest-neighbour Gaussian process approximation. These models and approximations were chosen to show that several aspects of the proposed multi-scale model, namely the use of multiple scales, the use of a large number of basis functions, and the ability to model nonstationarity, are crucial to get good predictions in a typical big-data environmental spatial data analysis.

Reproducible code and additional material, such as MCMC trace plots, are available from https://github.com/andrewzm/SpatialMultiScale.

5.1 Data

As data we used global sea-surface temperature (SST) obtained from the Visible Infrared Imaging Radiometer Suite (VIIRS) on board the Suomi National Polar-orbiting Partnership (Suomi NPP) weather satellite on October 14 2014 (Cao et al., 2013). We sampled one million data points at random from the complete data set of 238 million data points, intentionally leaving data out from an $8^{\circ}\times 8^{\circ}$ box centred at ($155^{\circ}$W, $0^{\circ}$N) in the Pacific Ocean, and used these as observation/training data. We then sampled another one million data points at random from the remainder and used these for assessing the quality of our predictions. After discarding data with invalid coordinate entries, and data points that fell outside our constructed meshes (see Sect. 5.2) we were left with about 900000 data points in each of the training and validation data sets.

We then detrended the training and validation data sets using a linear model with an intercept, the latitude coordinate, and the square of the latitude coordinate, as covariates. The residuals from the training data are what we call $\mathbf{Z}$ in Sect. 2, while we denote the residuals from the validation data set as $\mathbf{Z}_{v}$. In Fig. 3 we show the residuals $\mathbf{Z}$ on the globe together with the $8^{\circ}\times 8^{\circ}$ box (left panel) and a zoomed-in view of these residuals around Papua New Guinea (right panel). We can see from the figure that the spatial data are very irregularly-sampled in space, and that the distance from anywhere in the ocean to the nearest data point can range from a few kms to thousands of kms. The data also reveal both large-scale and small-scale features that need to be modelled.

5.2 Model Setup

We constructed the Delaunay triangulations on $\mathbb{S}^{2}$ with a coarse boundary around land masses (used in a study by Simpson et al., 2016) using the inla.mesh.create function in INLA, and fixed $\nu_{0}=\nu_{1}=1$. We determined the approximate number of vertices in the triangulations as follows, using the guidelines in Sect. 4. First, since it is required to factorise $\textbf{Q}_{0}$ we made sure that $T_{\eta_{0}}$ contains only a few tens of thousands of vertices, in our case it contains $r_{\eta_{0}}=38274$ vertices. Second, since $Y_{1}(\cdot)$ models the small-scale variability and we are only considering two scales, we require $T_{\eta_{1}}$ to be spatially dense with respect to the observed data. In our case we constructed a triangulation $T_{\eta_{1}}$ with 942349 vertices. Using this triangulation, the proportion of elements in ${\bm{\eta}}_{1}$ which have a data footprint size in the single digits is 92%. Finally, we want the nonstationarity in the fine-scale process to vary smoothly in space, and we therefore let $T_{\varepsilon}=T_{\theta_{1}}$ have only a few hundred vertices, in our case 205 vertices. Using this triangulation the largest effective process footprint is of size 27198, which is on the same order as $r_{\eta_{0}}$ as desired. The triangulations on a small region of the Pacific Ocean are shown in Fig. 4.

In order to re-sample ${\bm{\eta}}_{1}$ as discussed in Sect. 3.3, we constructed three partitionings of ${\bm{\eta}}_{1}$ from three spatial tilings of $D$. The first tiling was done using the third resolution of an Icosahedral Snyder Equal Area Aperture 3 Hexagonal (ISEA3H) grid on the surface of the sphere (Sahr, 2008). The second and third tilings were then done by shifting the ISEA3H grid north and east, respectively; see Fig. 2. To improve mixing of the Markov chain, spatial tiles that contained less than 100 data points or less than 200 basis functions at scale $k=1$, were merged with neighbouring tiles. The three tilings and colourings are shown in Fig. 5.

We constructed the prior distributions following the guidelines of Sect. 4. As prior distribution for $\rho_{0}$ we used a lognormal distribution such that $(\rho_{0,0.025},\rho_{0,0.975})=(300,10000)$ km, where $\rho_{0,q}$ denotes the $q$th quantile of the distribution. Our choice of $\rho_{0,0.025}$ stems from us having a maximum edge length of 150 km when constructing $T_{\eta_{0}}$, while our choice of $\rho_{0,0.975}$ presents a soft maximum on what we believe the spatial correlations in SSTs are (10000 km equates to approximately half the greatest east-west span of the Pacific Ocean). On the other hand, as prior distribution for ${\bm{\theta}}_{\rho_{1}}^{\{i\}},i=1,\dots,r_{\rho_{1}},$ we used a Gaussian distribution such that $(\exp({\bm{\theta}}_{\rho_{1},0.025}^{\{i\}}),\exp({\bm{\theta}}_{\rho_{1},0.975}^{\{i\}}))=(30,900)$ km. The lower quantile was selected on the basis that we used 30 km as our maximum edge size when constructing $T_{\eta_{1}}$, while the upper quantile was selected to ensure that there is sufficient spectral overlap between $Y_{0}(\cdot)$ and $Y_{1}(\cdot)$. As for the marginal process standard deviation, we used prior distributions such that $(\sigma_{0,0.025},\sigma_{0,0.975})=(\exp({\bm{\theta}}_{\sigma_{k},0.025}^{\{i\}}),\exp({\bm{\theta}}_{\sigma_{k},0.975}^{\{i\}}))=(0.1,6)^{\circ}$C, where our upper quantile was selected from the fact that the empirical standard deviation of the SST residuals is approximately equal to 2∘C. Finally, as prior for ${\bm{\theta}}_{\varepsilon}$ we used a Gaussian distribution such that $(\exp({\bm{\theta}}_{\varepsilon,0.025}^{\{i\}}),\exp({\bm{\theta}}_{\varepsilon,0.975}^{\{i\}}))=(0.5,5)^{\circ}$C, which reflects our prior belief that it is unlikely that the measurement error standard deviations are less than 0.5∘C or larger than $5^{\circ}$C.

As initial values for our parameters in the MCMC algorithm we used samples of ${\bm{\theta}}_{0},{\bm{\theta}}_{1},$ and ${\bm{\theta}}_{\varepsilon}$ drawn from their respective prior distributions. As initial values for ${\bm{\eta}}_{0}$ we used least-squares estimates from the model $\mathbf{Z}=\textbf{A}_{0}{\bm{\eta}}_{0}+\bm{\varepsilon}$, perturbed using independent samples from a Gaussian distribution with mean zero and variance 0.04. As initial values for ${\bm{\eta}}_{1}$ we used independent samples from a Gaussian distribution with mean zero and variance 0.01. Our adaptive random-walk Metropolis proposals were initialised to have variance 0.001; adaptation was done every 30 iterations for the first 2000 iterations.

The Gibbs sampler algorithm in Sect. 3.5 was run for 10000 iterations. The first 5000 samples were discarded as burn-in and the remaining 5000 were thinned by a factor of 50, to yield 100 samples from the posterior distribution over all the states and parameters. It took approximately 5 days of computing time to obtain 10000 samples using our implementation.

The long run-time of our algorithm made a detailed test of MCMC convergence impractical. Instead, we made an informal test of convergence by analysing samples of the process $Y(\cdot)$ evaluated at approximately 350000 points regularly interspersed throughout the ocean. Specifically, we computed the effective sample sizes (see Eq. 11.8 of Gelman et al., 2013) of the thinned chain at each of these locations, and compared the empirical density of these sizes to that constructed from independent samples from a Gaussian distribution, and to that constructed from traces that are first-order autoregressive processes with auto-correlation coefficient 0.1. In Fig. 6 we see that the empirical density is very sensitive to auto-correlation, and that the effective sample sizes from our sampler are virtually indistinguishable from those computed from independent samples. This indicates that most likely we do indeed have independent samples from the posterior distribution. Visual inspection of a random selection of trace plots also suggested the absence of auto-correlation.

Traces of other model components were not all as uncorrelated: While, for example, traces for ${\bm{\eta}}_{0}$, ${\bm{\eta}}_{1}$, and ${\bm{\theta}}_{\varepsilon}$, had median effective sample sizes of 87, 90, and 90, respectively, those for ${\bm{\theta}}_{\sigma_{1}}$ and ${\bm{\theta}}_{\rho_{1}}$ had median effective sample sizes of 52 and 46, respectively. While the latter low effective sample sizes imply that we should be cautious when making inferential statements on the nominal ranges and variances, they do not necessarily suggest that our process predictions could benefit from a longer run. This is a consequence of weak identifiability of the process model we are using; see Besag et al. (1995); Eberly and Carlin (2000) and Gelfand et al. (2001) for more discussion. Moreover, the purpose of our approach is to predict, and for prediction the ultimate goal is to maximise out-of-sample performance. In this respect the key finding is that our approach is well-calibrated out-of-sample, and that it considerably outperforms other approaches, as we describe in Sect. 5.4.

5.3 Other models

We compared the predictions of the proposed model (SPDE2) to several other models and approximations commonly used in applications containing large spatial data sets. We describe these in more detail below.

Global coarse-scale model (SPDE0):

We used INLA to fit the SPDE (Lindgren et al., 2011) in (1) to the entire data set via a GMRF on the coefficients of the tent basis functions on $T_{\eta_{0}}$. We let $\nu=1$, and used the same prior distributions for $\rho_{0}$ and $\sigma_{0}$ as we did for our two-scale spatial model in Sect. 5.2. We compare to this model to demonstrate the benefit of using a second, small-scale process when modelling such large spatial data sets. It took approximately 30 minutes of computing time to fit the model and obtain predictions with SPDE0. The nominal range for this model was estimated to be approximately 800 km.

Multiple spatially-independent small-scale processes (SPDE1-indep):

We used INLA to independently fit SPDEs with spatially-invariant parameters to data that lie in each of the tiles shown in Fig. 5, top panel. For each set of tile indices ${\mathcal{T}}$ in our model SPDE2, we generated a fine-resolution mesh that has a similar number of basis functions as $|{\mathcal{T}}|$, and a few more around the boundary of the partition to reduce boundary effects. We let $\nu=1$, and used the same prior distributions for (the now spatially-invariant) $\sigma$ and $\rho$ in each tile as we did for our two-scale process model in Sect. 5.2. These independent models were used to predict the noisy process at validation data locations that lie within their associated tile. Results from this model are used to demonstrate the benefit of having a large-scale process that can borrow strength across large spatial scales when modelling large spatial data. It took approximately 6 hours to sequentially fit the models and obtain predictions with SPDE1.

Multi-scale stationary process with a relatively small number of basis functions via LatticeKrig (LTK):

We used LatticeKrig (Nychka et al., 2015) to fit a 3-resolution LatticeKrig model on the cylinder comprising a total of 137577 basis functions. Note that many of these basis functions (approximately one third) are over land. We fitted three LTK models with the parameter a.wght, which dictates the spatial range, fixed to 4.01, 5.01, and 6.01, respectively. In Sect. 5.4 we only show results for the case a.wght = 6.01, which provided slightly better predictions than the other two cases (all cases gave very similar predictions). LTK models are generally limited to $r=100000$ to $200000$ basis functions since they require factorisation of precision matrices of dimension $r\times r$ when being fit. Results from this model are used to demonstrate that hundreds of thousands of basis functions are still likely to be insufficient when modelling large spatial data, despite the use of multiple scales. Fitting and predicting (via 30 conditional simulations) with LTK required approximately one day of computing time.

Single-scale process approximated using a full-scale approximation (FSA):

The full-scale approximation (Sang and Huang, 2012) to a Matérn covariance function with smoothness $\nu=1$ was implemented on $\mathbb{R}^{2}$ using 100 knots randomly placed on the surface of the sphere before projection onto the plane. The number of blocks on which to approximate the residual field was set to 9000 and data were attributed to each of these blocks using a k-means clustering algorithm on the lon–lat coordinates of the data. This choice ensured that the number of data points in each block was computationally feasible at about 100 data points per block. It took approximately two days to fit and predict using the FSA. The nominal range for this model was estimated to be approximately 850 km.

Single-scale process approximated using nearest neighbours (NNGP):

We used the R package spNNGP to implement a conjugate version of the nearest-neighbour Gaussian process (NNGP, Finley et al., 2019), where NNGPs are fitted for several plausible values of the latent-process range and variance, and cross-validation is used to select the best of these models for prediction. In each of these models an inverse-gamma prior distribution is used for the observation measurement-error variance, which yields closed-form predictive distributions that are quick to evaluate. In our implementation we set the number of neighbours to 15, the covariance function to a Matérn covariance function with $\nu=1$, and used a fine grid for the latent-process range and variance, ensuring that an interior point of this grid was selected as optimal. It took approximately three hours to find the optimal NNGP model through cross-validation and to predict using the optimal model. The nominal range for this model was estimated to be approximately 700 km.

Results from the FSA and NNGP are used to show that single-scale process models, which use local methods when predicting, are not competitive against models like SPDE2 which are able to borrow strength across large spatial scales when predicting over large gaps.

5.4 Results

In Fig. 7 we show predictions (posterior means) and prediction standard errors (posterior standard errors) for $Y(\cdot)$ using the two-scale spatial-process model SPDE2. The detail in these maps is very high – the one million basis functions effectively model the SST at a 30 $\times$ 30 km resolution globally. One should be cautious when interpreting the spatially-varying parameters (see Sect. 5.2), however it is evident that the posterior expectations of ${\bm{\theta}}_{\varepsilon},{\bm{\theta}}_{\sigma_{1}}$, and ${\bm{\theta}}_{\rho_{1}}$ vary considerably. In particular, the elements in $\mathrm{E}(\exp({\bm{\theta}}_{\sigma_{1}})\mid\mathbf{Z})$ had an interquartile range of $(0.22,0.85)^{\circ}$C while those in $\mathrm{E}(\exp({\bm{\theta}}_{\rho_{1}})\mid\mathbf{Z})$ had an interquartile range of $(40,220)$ km. This suggests that there is ample nonstationarity at small scales that is being captured by SPDE2. These posterior expectations can also be used to interpret the behaviour of the process at each scale. For example, the posterior expectations of the nominal ranges at $k=1$ (ranging from a few dozen to a few hundreds of kms) are smaller than $\mathrm{E}(\rho_{0}\mid\mathbf{Z})\approx 1700$ km, as expected due to the imposed prior distributions and triangulations. Also, $\mathrm{E}(\sigma_{0}\mid\mathbf{Z})=1.84^{\circ}$C, and therefore there is more power being allocated to the signal at the coarse scale than at the fine scale. Again, this was expected since power spectral densities tend toward zero at large frequencies. More intuition into the behaviour of the process at the separate scales can be obtained by visualising predictions of $Y_{0}(\cdot)$ and $Y_{1}(\cdot)$ separately; we show one such plot centred on the Brazil–Malvinas confluence in our online code repository.

We compare the predictions from all the models in terms of the root-mean-squared prediction error (RMSPE), continuous-ranked probability score (CRPS), 90% interval score (IS90), and the 90% coverage (Cov90) at the validation-data locations (Gneiting and Raftery, 2007). To assess the models’ ability to predict and quantify uncertainty correctly at validation-data locations that are in regions of both dense and sparse training data, we split the validation data set into three: $\mathbf{Z}_{v}^{(1)}$ contains all the validation data outside the $8^{\circ}\times 8^{\circ}$ box in the Pacific Ocean (see Sect. 5.1) that are in the vicinity of training data (specifically, a $1.5^{\circ}\times 1.5^{\circ}$ box centred on each datum in $\mathbf{Z}_{v}^{(1)}$ contains at least one training datum within it); $\mathbf{Z}_{v}^{(2)}$, which contains the validation data not in the vicinity of training data and outside the $8^{\circ}\times 8^{\circ}$ box; and $\mathbf{Z}_{v}^{(3)}$ that contains the validation data within the $8^{\circ}\times 8^{\circ}$ box. The dimensions of $\mathbf{Z}_{v}^{(1)},\mathbf{Z}_{v}^{(2)}$, and $\mathbf{Z}_{v}^{(3)}$ were $897282,293$, and $6964$, respectively.

Diagnostics of the marginal prediction distributions from each model are provided in Table 1. In the first six rows of Table 1 we show diagnostics for when prediction locations are in the vicinity of data points. As expected, the models that rely on a relatively small number of basis functions (SPDE0 and LTK) give relatively poor predictions. On the other hand, the use of a high number of basis functions (SPDE2 and SPDE1-indep) and the single-scale models that use local methods when predicting (FSA and NNGP) give better predictions. All methods give the correct coverage. What we deduce from these results is well-known: When doing spatial prediction on large data sets there is little benefit in using a global model when predicting at locations close to observed data (unless the signal-to-noise-ratio is very low; see, e.g., Zammit-Mangion et al., 2018) if the target of the spatial analysis is pointwise process prediction. There is also likely to be little benefit to be gained from using nonstationary process models since the inferences are predominantly data-driven and not model-driven. Rather, it is crucial to use a latent process model that is able to reproduce the fine-scale variation in the underlying process.

In the second six rows of Table 1 we show diagnostics for when prediction locations are not in the vicinity (through the $1.5^{\circ}\times 1.5^{\circ}$ boxes described earlier) of observations. The models LTK and SPDE0 are now not compromised as much for their use of relatively few basis functions. On the other hand, since FSA and NNGP utilise only a few data points when predicting, they begin to give poorer predictions. SPDE1-indep performs surprisingly well in this scenario, but slightly worse than SPDE0. All methods performed worse than SPDE2. Finally, in the last six rows of Table 1 we show diagnostics for validation data inside the $8^{\circ}\times 8^{\circ}$ lon-lat box, while in Fig. 8 we show the predictions at the data locations. We see that predicting using ‘nearest data’ generates unwanted artefacts in the case of the FSA and the NNGP. We expected LTK to perform well in this case; however, it appears to have over-fit at the region’s boundary, resulting in a large negative prediction in the region’s interior. The two-scale model outperforms the others by quite some margin in terms of RMSPE and CRPS since it is able to correctly predict the region of low SST anomaly traversing this box from west to east, while resolving fine scales on the region’s boundary. Interestingly all models except SPDE1-indep, and to a lesser extent SPDE2, are largely underconfident in this unobserved area.

6 Conclusion

The majority of models and inference techniques developed for analysing massive spatial data sets are not designed to exploit the multi-scale, nonstationary nature of the underlying process. The result is the widespread use of techniques that are able to provide good predictions at the fine scale (e.g., NNGP) or at the coarse scale when the number of basis functions for decomposition is capped for computational reasons, but not both. In this article we have proposed a multi-scale model where the degree of nonstationarity increases with the scale of the component processes, together with an approximate inference algorithm that is scaleable with both data and model size. We have shown that it outperforms other state-of-the-art approaches that are amenable to big data scenarios, but that are unable to capture the complexities of the underlying process due to either the use of a few basis functions, or the use of local methods that only take into account a small subset of the data when predicting. We conclude that using information at multiple scales is important for accurate prediction in environmental applications such as the analysis of SST, particularly in regions where remote sensing data is sparse. We used SPDEs to model the processes at each scale, but other spatial GMRFs (such as conditional autoregressive models) can be used instead. However, we anticipate that constructing appropriate grids, modelling the nonstationary structure, and eliciting interpretable prior distributions for the parameters would become more challenging.

There are several avenues for future work. First, while the approximation (11) is a sensible one under the guidelines of Sect. 4, we do not provide any theoretical bound on the effect of the approximation on the target distribution, which is not straightforward to derive. Second, while our implementation in Sect. 5 used chunked-up data and models when doing inference, these were saved and loaded to disk when sampling via graph colouring. Sophisticated implementations include one where static objects (e.g., data and basis-function matrices) reside permanently in memory on dedicated nodes (Katzfuss and Hammerling, 2017), while communication of selected samples between the nodes is done efficiently via message passing (e.g., using MPI). Such an implementation would speed up the Gibbs sampler, and is necessary for the consideration of higher scales and data sets that are one to two orders of magnitude larger than that considered in Sect. 5. Third, there is considerable interest in making inference with large non-Gaussian spatial data. Parallelisation is still possible for this case, but our computational framework would need to be modified slightly since the process coefficients cannot be sampled directly when the data is non-Gaussian. Specifically, when carrying out an update corresponding to a colour, the process coefficients would have to be proposed and then accepted or rejected jointly with the parameter coefficients; see, for example, Knorr-Held and Rue (2002).

Acknowledgements

We thank Yuliya Marchetti for providing the sea-surface temperature data set, Bohai Zhang for providing the MATLAB code for the implementation of the FSA, Matt Moores for general discussions on improving MCMC mixing, and Quan Vu for providing comments on an early version of this manuscript. AZ–M was supported by the Australian Research Council (ARC) Discovery Early Career Research Award, DE180100203.

Appendix A Targeted distribution of the Markov chain

In Sect. 3.2 we asserted that it is important to update ${\bm{\eta}}_{k}$ concurrently with ${\bm{\theta}}_{k}$, even though ${\bm{\eta}}_{k}$ is later re-updated (Sect. 3.3). If one does not do this, ${\bm{\eta}}_{k}^{\mathrm{rest}}$ in (17) would be ‘out of sync’ with the updated parameters ${\bm{\theta}}_{k}$; as a consequence, when ${\bm{\eta}}_{k}^{{\mathcal{T}}_{1}}$ or ${\bm{\eta}}_{k}^{{\mathcal{T}}_{2}}$ is updated in (17), an incorrect distribution would be targeted. This phenomenon occurs when marginalising (termed ‘trimming’ by Van Dyk and Park, 2008) in Gibbs samplers.

We show the importance of resampling on a very simple spatial model, where we have two sets of parameters, ${\bm{\theta}}_{1}$ and ${\bm{\theta}}_{2}$, and two sets of basis-function coefficients, ${\bm{\eta}}_{1}$ and ${\bm{\eta}}_{2}$. In what follows we omit conditioning on the data $\mathbf{Z}$, since this is implicit in all the distributions. We denote the target (posterior) distribution as $\operatorname{\mathrm{pr}}_{0}({\bm{\Theta}})$ where ${\bm{\Theta}}\equiv\{{\bm{\theta}}_{1},{\bm{\theta}}_{2},{\bm{\eta}}_{1},{\bm{\eta}}_{2}\}$. In MCMC we seek a transition kernel $q({\bm{\Theta}}^{\prime}\mid{\bm{\Theta}})$ such that

[TABLE]

If (A.20) holds, then we say that the Markov chain preserves the target distribution, $\operatorname{\mathrm{pr}}_{0}$.

In a vanilla Gibbs sampler, one constructs the transition kernel from full conditional distributions of the target distribution:

[TABLE]

Successive updating of the parameters in this fashion preserves the target distribution. To see this, substitute (A.21) in (A.20) to obtain

[TABLE]

as required. Now, this vanilla sampler does not mix well due to the correlation a posteriori between ${\bm{\eta}}_{i}$ and ${\bm{\theta}}_{i}$, $i=1,2$ (Knorr-Held and Rue, 2002). Since our model is in a linear, Gaussian, setting, one might be tempted to instead use the following transition kernel:

[TABLE]

A similar treatment to the vanilla Gibbs case reveals that this only targets the correct distribution if $\operatorname{\mathrm{pr}}_{0}({\bm{\theta}}_{1},{\bm{\theta}}_{2},{\bm{\eta}}_{1},{\bm{\eta}}_{2})=\operatorname{\mathrm{pr}}_{0}({\bm{\theta}}_{1},{\bm{\theta}}_{2})\operatorname{\mathrm{pr}}_{0}({\bm{\eta}}_{1},{\bm{\eta}}_{2})$, which is almost certainly not the case in our spatial models. Therefore, (A.22) is not a target-preserving kernel. In our MCMC scheme this is important: updating ${\bm{\theta}}_{k}$ in blocks and subsequently updating ${\bm{\eta}}_{k}$ in blocks will not yield samples from the posterior distribution.

A kernel which preserves the target can be constructed by updating ${\bm{\eta}}_{1}$ and ${\bm{\eta}}_{2}$ twice, with the intermediate quantities then discarded. As in Algorithm 3.5, denote these intermediate quantities as ${\bm{\eta}}_{1}^{*}$ and ${\bm{\eta}}_{2}^{*}$, respectively, and consider the transition kernel

[TABLE]

This kernel preserves the target (posterior) distribution since

[TABLE]

as required.

Appendix B Simulation experiment illustrating the benefit of alternating tilings

Consider a GMRF ${\bm{\eta}}\sim\mathrm{Gau}(\mathbf{0},\textbf{Q}^{-1})$, where Q is the sparse tridiagonal matrix

[TABLE]

where omitted entries are zero, $\phi$ is a length-scale parameter, and $\sigma^{2}_{v}$ a variance parameter. We consider the case where ${\bm{\eta}}\in\mathbb{R}^{n}$ with $n=99$, and consider the following tilings,

[TABLE]

Here, we compare Markov chain behaviour when sampling $N$ times from this distribution using the two samplers, Sampler 1 and Sampler 2, detailed in Algorithms B and B, respectively. Sampler 1 is a blocked Gibbs sampler which samples ${\bm{\eta}}$ using only ${\mathcal{T}}_{11}$ and ${\mathcal{T}}_{12}$, while Sampler 2 changes the tiling used, $\{{\mathcal{T}}_{11},{\mathcal{T}}_{12}\}$ or $\{{\mathcal{T}}_{21},{\mathcal{T}}_{22},{\mathcal{T}}_{23}\}$, at each iteration.

In our simulation experiment we simulated ${\bm{\eta}}$ using $\phi=0.9$ and $\sigma^{2}_{v}=0.2$, and let $N=10000$. We then generated two Markov chains, one using Sampler 1 and another using Sampler 2, and for each chain took the last 5000 samples and applied a thinning factor of 2. In Fig. 9 we show the empirical auto-correlation functions computed from the trace plots of ${\bm{\eta}}^{\{49\}}$ from both samplers. We see that samples of ${\bm{\eta}}^{\{49\}}$ from Sampler 1 are highly correlated due to the proximity of this variable to the tiling boundary. Samples of ${\bm{\eta}}^{\{49\}}$ from Sampler 2 are virtually uncorrelated. Hence, a system of shifting tiles in a Gibbs sampler for spatial GMRFs (as done in Algorithm 3.5) can virtually eliminate any auto-correlation that may appear due to tile boundaries. Note that a thinning factor greater than the number of tilings needs to be used to effectively remove any auto-correlation.

Appendix C Simulation experiment illustrating the sensitivity of the predictions to the chosen basis functions

In this section we conduct a simulation experiment that demonstrates the effect of a coarse discretisation, and the corresponding basis-function representation, on the prediction performance of the multi-scale model. The experiment is done in a one-dimensional, two-scale-process setting.

Consider a 1D Gaussian process on $D=[0,1]$, which has as covariance function a sum of two exponential covariance functions, $C_{0}(\cdot)$ and $C_{1}(\cdot)$. The exponential covariance functions have as parameters the variances $\sigma^{2}_{k}$ and ranges $\tau_{k}$, $k=0,1$, and are given by

[TABLE]

We model the process of interest $Y(\cdot)$ as a sum of the two processes $Y_{k}(\cdot)=\textbf{a}_{k}(\cdot)^{\top}{\bm{\eta}}_{k},k=0,1$, where $\textbf{a}_{k}(\cdot)$ are basis functions and ${\bm{\eta}}_{k}$ are basis-function coefficients. Now, let the basis functions $\textbf{a}_{k}(\cdot),k=0,1,$ be piecewise constants on a regular partitioning of $D$. The basis functions have width $\delta_{0}$ and $\delta_{1}$, respectively, where $\delta_{1}<\delta_{0}$. We are interested in what the effect of a poor choice for $\delta_{0}$ and $\delta_{1}$ is when predicting $Y(\cdot)$ from noisy observations $\mathbf{Z}$.

If $\delta_{k}\rightarrow 0$, then the $k$th scale of the original process is reconstructed exactly if we let ${\bm{\eta}}_{k}\sim\mathrm{Gau}(\mathbf{0},\textbf{Q}_{k}^{-1}),k=0,1,$ where $\textbf{Q}_{k}$ is the (sparse) tridiagonal matrix given by (B.23) with $\sigma^{2}_{v}$ replaced with $\sigma^{2}_{v,k}$ and $\phi$ replaced with $\phi_{k}$. Here, $\phi_{k}=\exp(-\delta_{k}/\tau_{k})$ and $\sigma^{2}_{v,k}=\sigma^{2}_{k}(1-\phi_{k}^{2})$. In a simulation environment where we have access to $\sigma^{2}_{k},\tau_{k},k=0,1$, we can therefore obtain accurate GMRF representations for our processes at the individual scales for when $\delta_{k}$ is small. We can then also see what happens as $\delta_{k}$ grows (this corresponds to coarsening a triangulation in 2D).

In our simulation environment we fixed $\tau_{0}=0.4,\tau_{1}=0.04,\sigma^{2}_{0}=1$ and $\sigma^{2}_{1}=0.05$ and conducted 100 Monte Carlo simulations, where in each simulation we did the following:

Randomly established an ‘unobserved’ region in $D$, $D_{\operatorname*{gap}}$ say, where $|D_{\operatorname*{gap}}|=0.2$. 2. 2.

Generated 1100 observations on $D$, with 1000 in $D\backslash D_{\operatorname*{gap}}$ and 100 in $D_{\operatorname*{gap}}$ with measurement-error variance $\sigma^{2}_{\varepsilon}=0.0002$. Five hundred of those in $D\backslash D_{\operatorname*{gap}}$ were used as training data $\mathbf{Z}$, the remaining 500 in $D\backslash D_{\operatorname*{gap}}$ as validation data $\mathbf{Z}_{v}^{(1)}$, and those 100 in $D_{\operatorname*{gap}}$ as validation data $\mathbf{Z}_{v}^{(2)}$. 3. 3.

For various values of $\delta_{0}$ and $\delta_{1}$, we constructed $\textbf{Q}_{0}$ and $\textbf{Q}_{1}$ according to (B.23) and used these, the true measurement-error variance, and $\mathbf{Z}$, to predict $Y(\cdot)$ at all validation data locations. 4. 4.

Computed the RMSPE and CRPS at the validation locations.

Each Monte Carlo simulation provided us with an RMSPE and a CRPS corresponding to a combination of $\{\delta_{0},\delta_{1}\}$. We then averaged over the 100 simulations to provide averaged RMSPEs and CRPSs corresponding to each combination of $\{\delta_{0},\delta_{1}\}$. This experiment allows us to analyse the detrimental effect of a large $\delta_{0}$ or $\delta_{1}$ on our predictions.

The results from this experiment are summarised in Fig. 10. The figure clearly shows that the RMSPE considerably increases in regions where data is dense ($\mathbf{Z}_{v}^{(1)}$, left panels) and $\delta_{1}$ is large; on the other hand $\delta_{0}$ does not have much of an effect in these regions. The situation is reversed in regions where data is missing in large contiguous blocks ($\mathbf{Z}_{v}^{(2)}$, right panels). Here $\delta_{1}$ does not play a big role while $\delta_{0}$ does. When doing simple kriging (with the exact model) the mean RMSPEs were 0.07 and 0.41, respectively, while the mean CRPSs were 0.033 and 0.24, respectively, which are relatively close to what was obtained with the smallest values we chose for $\delta_{0}$ and $\delta_{1}$. Therefore, the way in which we discretise both scales are important, and both $\delta_{0}$ and $\delta_{1}$ should be made as small as needed for their respective scale; in this experiment $\delta_{1}=0.001$ and $\delta_{0}=0.01$ are suitable choices. In practice, the coarseness of the grid (or triangulation in 2D) will be determined through computational considerations. Fortunately, we see that predictive performance deteriorates at an reasonably slow rate as the discretisations get coarser and coarser. A detailed analysis taking into account convergence rates of finite-element approximations might be needed for an in-depth analysis.