# Grid-Parametrize-Split (GriPS) for Improved Scalable Inference in Spatial
Big Data Analysis

Michele Peruzzi111Department of Statistical Science, Duke University , Sudipto
Banerjee222Department of Biostatistics, UCLA Fielding School of Public Health,
David B. Dunson††footnotemark: , Andrew O. Finley333Department of Forestry,
Michigan State University

###### Abstract

Rapid advancements in spatial technologies including Geographic Information
Systems (GIS) and remote sensing have generated massive amounts of spatially
referenced data in a variety of scientific and data-driven industrial
applications. These advancements have led to a substantial, and still
expanding, literature on the modeling and analysis of spatially oriented big
data. In particular, Bayesian inferences for high-dimensional spatial
processes are being sought in a variety of remote-sensing applications
including, but not limited to, modeling next generation Light Detection and
Ranging (LiDAR) systems and other remotely sensed data. Massively scalable
spatial processes, in particular Gaussian processes (GPs), are being explored
extensively for the increasingly encountered big data settings. Recent
developments include GPs constructed from sparse Directed Acyclic Graphs
(DAGs) with a limited number of neighbors (parents) to characterize dependence
across the spatial domain. The DAG can be used to devise fast algorithms for
posterior sampling of the latent process, but these may exhibit pathological
behavior in estimating covariance parameters. While these issues are mitigated
by considering marginalized samplers that exploit the underlying sparse
precision matrix, these algorithms are slower, less flexible, and oblivious of
structure in the data. The current article introduces the Grid-Parametrize-
Split (GriPS) approach for conducting Bayesian inference in spatially oriented
big data settings by a combination of careful model construction and algorithm
design to effectuate substantial improvements in MCMC convergence. We
demonstrate the effectiveness of our proposed methods through simulation
experiments and subsequently undertake the modeling of LiDAR outcomes and
production of their predictive maps using G-LiHT and other remotely sensed
variables.

Keywords: Big data; Bayesian modeling; directed acyclic graphs; Gaussian
processes; sparse modeling.

##  1 Introduction

Geographic Information Systems (GIS) and related technologies such as remote
sensors, satellite imaging and portable devices that are capable of collecting
precise positioning information have spawned much interest in the production
of digitized geographic information, which, in turn, has led to the amassing
of massive amounts of spatial data. Rapidly increasing availability and
capabilities of remote sensing technologies and GIS have brought about
substantial developments in modeling and analyzing spatial data in diverse
scientific and data-driven industrial applications including, but not limited
to, natural and environmental sciences; economics; climate science; ecology;
forestry; and public health. With the abundance of spatial big data problems
in the sciences and engineering, GIS and spatial data science will likely
occupy a central place in the data revolution engulfing us.

With the scientific and technological communities embracing an era where open-
access data-rich environments provide extraordinary opportunities to
understand spatial and temporal complexity of processes, statisticians are
increasingly encountering inferential questions that require complex and
computationally challenging models. Analysts working with spatial technologies
customarily seek to (i) estimate the variations and associations among several
spatially referenced outputs; and (ii) predict outputs at new spatial
locations. A generic paradigm for spatial models posits a joint probability
law

$[\mbox{data}\,|\,\mbox{process}]\times[\mbox{process}\,|\,\mbox{parameters}]\times[\mbox{parameters}]\;,$
(1)

where each $[\cdot]$ represents a legitimate probability law. The above posits
an underlying “process” that generates the data. Modeling this process in a
rigorous and computationally efficient manner without sacrificing richness in
inference presents several challenges. How exactly the underlying process is
specified depends upon the data attributes relevant to inference. We focus on
point-referenced data, where spatial locations are indexed by coordinates on a
map. The “process” is modeled as a spatial random field over the domain of
interest and the observations are treated as a finite realization of this
random field. Endowing the finite-dimensional realizations of a random field
with a probability law that is multivariate normal leads to the pervasive
Gaussian process (GP).

The GP is, perhaps, the most conspicuous of process specifications and offers
flexibility and richness in modeling. However, GPs incur onerous computational
costs that make them impracticable for spatial big data analysis. The primary
burden stems from storage and computations involving the massive spatial
covariance matrices representing the dependencies inside the data. The
elements of these matrices are covariance functions, or kernels, depending on
parameters that describe the process. These covariance functions are nonlinear
and, depending upon the scientific hypotheses the model represents (e.g.,
nonstationary dependence), can be computationally burdensome. For irregularly
situated spatial locations, as is common in geostatistics, these covariance
matrices are typically dense and offer no exploitable structure to facilitate
computations. Even for a modestly large number of points by today’s norms
($\approx 50,000$+), the computational demands are prohibitive. Yet, spatial
analysts encounter data sets with locations numbering in the millions. The
size of such data sets easily exceed common CPU and memory configurations,
which means we need to rely upon statistical models to enable such analyses
within available computing resources.

There exists a very large literature on scalable GP approximations for
extending spatial and spatiotemporal inference to big data including, but not
limited to, low-rank methods (Quiñonero-Candela and Rasmussen, 2005; Snelson
and Ghahramani, 2007; Banerjee et al., 2008; Cressie and Johannesson, 2008),
covariance tapering (Furrer et al., 2006; Kaufman et al., 2008), domain
partitioning (Sang and Huang, 2012; Stein, 2014), divide-and-conquer methods
based on meta-posteriors (Guhaniyogi and Banerjee, 2018), recursive
partitioning and hierarchical matrices (Ambikasaran et al., 2016; Katzfuss,
2017; Huang and Sun, 2018; Geoga et al., 2020; Jurek and Katzfuss, 2020),
Gaussian random Markov random fields (GMRF; Rue and Held, 2005), composite
likelihood methods (Eidsvik et al., 2014), and local approximations (Gramacy
and Apley, 2015). We refer the reader to Sun et al. (2011), Banerjee (2017),
and Heaton et al. (2019) for in-depth reviews and comparisons of scalable
geostatistical methods.

One particular approach that has been garnering much attention in spatial and
spatiotemporal big data analysis stems from representing the joint probability
law of a finite realization of the spatial process using a directed acyclic
graph (DAG) constructed from the spatial locations as the nodes of the graph
with parents of each node as spatial neighbors. Vecchia (1988) used this idea
as a method to approximate the multivariate Gaussian likelihood function
corresponding to a realized GP over a finite but large number of locations.
Datta et al. (2016a) extends such finite-dimensional approximations to a well-
defined spatial process, the Nearest-Neighbor Gaussian Process (NNGP), over
continuous spatial domains by defining a set of reference locations over which
the DAG is constructed. The conditional independence of the outputs given
neighboring outputs that is imposed by the DAG yields a sparse precision
(inverse of covariance) matrix resulting in a sparse likelihood approximation
to the full GP likelihood (often called Vecchia’s approximation in the spatial
literature).

Alternative methods to build the DAG have been considered (Quiroz et al.,
2019; Katzfuss and Guinness, 2019). Recently, Peruzzi et al. (2020) introduced
a class of meshed Gaussian processes (MGPs) for massive spatial data sets.
MGPs are constructed by mapping the nodes of a pre-determined directed acyclic
graph (DAG) to groups of spatial locations. Scalability to big data follows
from special choices of DAGs; treed DAGs for multivariate outcomes can be
associated to recursive domain partitioning (Peruzzi and Dunson, 2020),
whereas DAGs with repeated patterns arise from domain tiling or axis-parallel
tessellations. In the latter case, Peruzzi et al. (2020) show that additional
computational speedups correspond to data on partly observed lattices when
sampling the latent process a posteriori.

However, inefficiencies associated with excessive data augmentation may arise
when the data are not regularly spaced. Bayesian inference about the latent
process and the parameters governing it can be carried out by collecting
dependent samples from their joint posterior distribution using Markov chain
Monte Carlo (MCMC) methods. Features of the posterior distribution for the
parameters of interest (e.g. their mean) are estimated using such dependent
samples. Efficiency of MCMC depends on (i) the computational complexity (in
terms of flops and time) of the algorithm used to generate a sample of size
$T$, and (ii) the effective sample size (ESS) of such a sample. The ESS (which
is typically $T^{*}<T$) is the size of an i.i.d. sample from the posterior
distribution with the same variance as the sample from the Markov chain. Let
$\\{\theta_{1},\dots,\theta_{T}\\}$ be the posterior samples of a parameter
$\theta$ obtained upon convergence of an MCMC algorithm. Then, ESS is found as
$T^{*}=T/\kappa(\theta)$, where $\kappa(\theta)$ is the autocorrelation time
associated to the sequence $\theta_{t}$, $t=1,\dots,T$. Therefore, ESS
measures the loss in efficiency due to the use of a Markov chain (Robert and
Casella, 2005). A slowly mixing Markov chain, which exhibits high
autocorrelations between successive samples, is linked to a lower ESS. Markov
chains that mix poorly require substantially longer time to produce Monte
Carlo samples that accurately characterize the posterior distribution.
Unfortunately, mixing time increases with the dimension of the state space
(see e.g. Dunson and Johndrow, 2020), which is particularly troubling in big
data settings.

Poor mixing of MCMC in spatial contexts negatively impacts the quality of
estimations and uncertainty quantification of the unknown covariance
parameters when sampling the high-dimensional latent process. Marginalized or
collapsed samplers exhibit considerable improvements in mixing because they
obviate sampling the high-dimensional spatial effects, which are subsequently
sampled _after_ the posterior samples of the collapsed model have been
obtained. Collapsed samplers are, however, effective only for Gaussian
likelihoods and tend to be much slower as they entail more expensive matrix
computations (Finley et al., 2019). Furthermore, the practical performance of
these algorithms in quantifying posterior uncertainty about the covariance
parameters may still be unsatisfactory. Improvements compared to (faster) non-
marginalized samplers may be small. Coupled with the observation that ESS is
frequently much higher for the latent process than it is for covariance
parameters when considering unmarginalized samplers, one may attribute poor
mixing to the absence of consistent estimators for the parameters in the
spatial covariance function. This has been studied for the Matérn covariance
kernel for in-fill or fixed-domain spatial asymptotics (Stein, 1990; Ying,
1991; Zhang, 2004), including in models with the nugget (Tang et al., 2020).
In the Bayesian context, inconsistent estimation implies that the effect of
the prior does not diminish even as we increase our sample size by populating
the fixed spatial domain with locations. It is also well known that model
parametrizations can impact performance of MCMC (Gelfand et al., 1995; Liu and
Wu, 1999; van Dyk and Meng, 2001; Yu and Meng, 2011; Kastner, 2014). The
connections between parametrizations and efficiency in spatial big data
contexts have not been analyzed, but we expect them to play an important role.
For example, a linear model of coregionalization (LMC; Matheron, 1982;
Wackernagel, 2003) for modeling multivariate spatial dependence in a Bayesian
regression model (see e.g. Zhang and Banerjee, 2020) is typically associated
with a “noncentered” parametrization (Papaspiliopoulos et al., 2007) that
inhibits posterior sampling efficiency in scenarios that are common in
geostatistical settings, e.g. when the measurement error variance is small
relative to the process variance, but is otherwise computationally convenient.

In this article, we combine three separate ideas to resolve the aforementioned
issues. First, we enforce a gridded (Gri) reference set on an MGP model to
ensure sizable scalability improvements for non-gridded multivariate outcomes
that may be misaligned. Second, we parametrize (P) the LMC via Matérn
covariances targeting explicit posterior sampling of the microergodic
parameter (Stein, 1999), which is well identified. Third, we overparametrize
the process covariance via a parameter splitting (S) strategy related to
half-t priors (Gelman, 2006). The latter two ideas are based on known results
but, to the best of our knowledge, their impact on multivariate spatial
regression models has not been considered. In the Matérn model, P and S imply
that the sill $\sigma^{2}$, range $1/\phi$, and smoothness $\nu$ all have
mixed centering. We then demonstrate that the combination of Grid,
Parametrize, and Split, leading to a GriPS regression, results in substantial
improvements in efficiency even when sampling very high-dimensional latent
effects.

Our proposed methodology and related software advancements are aimed for broad
application in global-scale public and private land surface monitoring efforts
that rely on next generation Light Detection and Ranging (LiDAR) systems and
other remotely sensed data. LiDAR provides high-dimensional characterization
of land surface structures at point-referenced locations, e.g., topography,
vegetation, ice/snow. Repeat measurement from these systems, and coupling with
other remotely sensed data, facilitate change detection and change
attribution, e.g., natural and human induced disturbances. Next generation
LiDAR systems include the National Science Foundations’s National Ecological
Observatory Network (NEON) Airborne Observation Platform (Kampe et al., 2010),
National Aeronautics and Space Administration (NASA) ICESAT-2 (ICESat-2,
2020), Global Ecosystem Dynamics Investigation (GEDI; GEDI 2020), and Goddard
LiDAR, Hyper-Spectral, and Thermal imager (G-LiHT, Cook et al. 2013). The GEDI
mission sensor, mounted on the International Space Station (ISS), began
collecting data in 2019 and, like ICESAT-2 and G-LiHT, is key to global-scale
mapping of forest change and quantifying forest loss contribution to
atmospheric CO2 concentrations. Unlike traditional LiDAR data collection
campaigns that exhaustively sample small regions, next generation LiDAR
systems such as ICESAT-2, G-LiHT, and GEDI cover large spatial extents by
sparsely sampling locations across the domain. While generating complete-
coverage high resolution maps of forest outcomes is still of interest, there
is also interest in creating predictive maps of the high-dimensional LiDAR
signals over non-sampled locations.

The remainder of the article proceeds as follows. Section 2 provides an
outline of a Bayesian multivariate spatial regression model and discusses the
major bottlenecks that detract from efficiency in MCMC samplers for current
state-of-the-art models in big data settings. We then introduce GriPS in
Section 3 to circumvent such bottlenecks. Section 4.1 outlines the posterior
sampling algorithm whose efficiency is demonstrated in Sections 4.2 and 5
through simulation experiments and subsequent analysis of LiDAR outcomes and
production of their predictive maps using G-LiHT and other remotely sensed
variables.

##  2 Posterior sampling (in)efficiency in spatial regressions

###  2.1 Bayesian multivariate spatial regression models

Let ${\cal D}\subset\Re^{d}$ denote our spatial domain and
$\boldsymbol{\ell}\in{\cal D}$ be a generic spatial location. A multivariate
spatial regression model describes the relationship between a vector of
spatial outputs and predictors while incorporating a multivariate spatial
process and white noise. Thus,

$\boldsymbol{y}(\boldsymbol{\ell})=\boldsymbol{X}(\boldsymbol{\ell})^{\top}\bbeta+\boldsymbol{w}(\boldsymbol{\ell})+\boldsymbol{\varepsilon}(\boldsymbol{\ell})\;,$
(2)

where
$\boldsymbol{y}(\boldsymbol{\ell})=(y_{1}(\boldsymbol{\ell}),y_{2}(\boldsymbol{\ell}),\ldots,y_{q}(\boldsymbol{\ell}))^{\top}$
is the $q\times 1$ vector of outcomes at location $\boldsymbol{\ell}$,
$\boldsymbol{X}(\boldsymbol{\ell})^{\top}=\text{blockdiag}\\{\boldsymbol{x}_{i}(\boldsymbol{\ell})^{\top}\\}_{i=1}^{q}$
is a $q\times p$ matrix of spatially referenced predictors with each
$\boldsymbol{x}_{i}(\boldsymbol{\ell})$ being $p_{i}\times 1$ vector of
predictors corresponding to $y_{i}(\boldsymbol{\ell})$ at location
$\boldsymbol{\ell}$ and $p=\sum_{i=1}^{q}p_{i}$, $\bbeta$ is the $p\times 1$
vector of corresponding regression coefficients,
$\boldsymbol{w}(\boldsymbol{\ell})$ and
$\boldsymbol{\varepsilon}(\boldsymbol{\ell})$ are $q\times 1$ vectors of
spatial random effects and random noise with elements
$w_{i}(\boldsymbol{\ell})$ and $\epsilon_{i}(\boldsymbol{\ell})$,
respectively, for $i=1,2,\ldots,q$ such that
$\boldsymbol{\varepsilon}(\boldsymbol{\ell})\overset{iid}{\sim}N(0,\boldsymbol{D})$
with a diagonal $\boldsymbol{D}=\text{diag}(\tau^{2}_{1},\dots,\tau^{2}_{q})$.

We model the uncountable set
$\\{\boldsymbol{w}(\boldsymbol{\ell}):\boldsymbol{\ell}\in{\cal D}\\}$ as a
$q$-variate multivariate Gaussian process, denoted as
$\boldsymbol{w}(\boldsymbol{\ell})\sim
GP(\mathbf{0},\boldsymbol{C}_{\boldsymbol{\theta}}(\cdot,\cdot))$, where
$\boldsymbol{C}_{\boldsymbol{\theta}}(\cdot,\cdot)$ is a $q\times q$ matrix-
valued _cross-covariance_ function indexed by parameters $\boldsymbol{\theta}$
such that
$\boldsymbol{C}_{\boldsymbol{\theta}}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=[\text{cov}\\{w_{i}(\boldsymbol{\ell}),w_{j}(\boldsymbol{\ell}^{\prime})\\}]$
is the $q\times q$ matrix with $(i,j)$th element given by the covariance
between $w_{i}(\boldsymbol{\ell})$ and $w_{j}(\boldsymbol{\ell}^{\prime})$.
Henceforth we suppress $\boldsymbol{\theta}$ in the cross-covariance to
simplify notation. A valid cross-covariance satisfies: (i)
$\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=\boldsymbol{C}(\boldsymbol{\ell}^{\prime},\boldsymbol{\ell})^{\top}$;
and (ii)
$\sum_{i=1}^{n}\sum_{j=1}^{n}\boldsymbol{z}_{i}^{\top}\boldsymbol{C}(\boldsymbol{\ell}_{i},\boldsymbol{\ell}_{j})\boldsymbol{z}_{j}>0$
for any integer $n$ and any finite collection of points
$\\{\boldsymbol{\ell}_{1},\boldsymbol{\ell}_{2},\ldots,\boldsymbol{\ell}_{n}\\}$
and for all $\boldsymbol{z}_{i}\in\Re^{q}\setminus\\{\boldsymbol{0}\\}$ (see,
e.g., Genton and Kleiber, 2015). Over any finite set ${\cal L}\subset{\cal D}$
the process has a multivariate normal distribution $\boldsymbol{w}_{{\cal
L}}\sim N(0,\boldsymbol{C}_{{\cal L}})$, where $\boldsymbol{w}_{{\cal L}}$ is
the $qn_{{\cal L}}\times 1$ column vector and $\boldsymbol{C}_{{\cal L}}$ is
the $qn_{{\cal L}}\times qn_{{\cal L}}$ block matrix with the $q\times q$
matrix $\boldsymbol{C}(\boldsymbol{\ell}_{i},\boldsymbol{\ell}_{j})$ as its
$(i,j)$ block for $i,j=1,\dots,n_{{\cal L}}$.

Let ${\cal
T}=\\{\boldsymbol{\ell}_{1},\dots,\boldsymbol{\ell}_{n}\\}\subset{\cal D}$ be
the set of locations where at least one of our $q$ spatial outcomes has been
measured. We construct $\boldsymbol{y}_{{\cal
T}}=[\boldsymbol{y}(\boldsymbol{\ell}_{1})^{\top},\dots,\boldsymbol{y}(\boldsymbol{\ell}_{n})^{\top}]^{\top}$
as the $nq\times 1$ vector of $\boldsymbol{y}(\boldsymbol{\ell}_{i})$’s over
the $n$ locations, analogously define $\boldsymbol{w}_{{\cal T}}$ and
$\boldsymbol{\varepsilon}_{{\cal T}}$, and let $\boldsymbol{X}_{{\cal
T}}=[\boldsymbol{X}(\boldsymbol{\ell}_{1}):\ldots:\boldsymbol{X}(\boldsymbol{\ell}_{n})]^{\top}$.
A Bayesian hierarchical model is constructed from the above specifications as

$p(\bbeta,\boldsymbol{w}_{{\cal
T}},\boldsymbol{\tau},\boldsymbol{\theta}\,|\,\boldsymbol{y}_{{\cal
T}})\propto p(\boldsymbol{\tau},\boldsymbol{\theta})\times
N(\bbeta;\;\boldsymbol{m}_{\beta},\boldsymbol{V}_{\beta})\times
N(\boldsymbol{w}_{{\cal T}}\;;\mathbf{0},\boldsymbol{C}_{{\cal T}})\times
N(\boldsymbol{y}_{{\cal T}}\;;\boldsymbol{X}_{{\cal
T}}\bbeta+\boldsymbol{w},\boldsymbol{D}_{n})\;,$ (3)

where $\boldsymbol{\tau}=\\{\tau_{1}^{2},\dots,\tau_{q}^{2}\\}$ and
$\boldsymbol{D}_{n}=\text{blockdiag}(\\{\boldsymbol{D}\\}_{i=1}^{n})$. The
likelihood in (3) corresponds to
$\boldsymbol{y}=\boldsymbol{X}\bbeta+\boldsymbol{w}+\boldsymbol{\varepsilon}$,
where the specification over ${\cal T}$ is implicit.

Bayesian inference proceeds by sampling from the posterior distribution in (3)
using MCMC algorithms. Efficiency of MCMC algorithms can be measured by
effective sample size (ESS) per unit time. ESS is related to the
autocorrelation in the Markov chains and measures the size of an equivalent
i.i.d. sample, i.e., producing the same Monte Carlo error in ergodic averages,
to the dependent MCMC samples (Liu and Chen, 1995; Robert and Casella, 2005).
Scalable models effectively reduce the time for each MCMC iteration, but
issues adversely affecting efficient sampling persist. We consider, in
particular, the efficient inference for the cross-covariance parameters
$\boldsymbol{\theta}$ and the measurement error variances $\boldsymbol{\tau}$.

###  2.2 Parametric specifications for multivariate spatial processes

$\boldsymbol{y}$$\boldsymbol{w}$$\boldsymbol{\theta}$$\bbeta,\boldsymbol{\tau}$
(a) Centered parametrization

$\boldsymbol{y}$$\boldsymbol{\Lambda}$$\boldsymbol{\omega}$$\boldsymbol{\varphi}$$\bbeta,\boldsymbol{\tau}$
(b) Noncentered parametrization

Figure 1: Common parametrizations of Bayesian spatial regression models.

Modeling the latent $q$-variate process $\boldsymbol{w}(\cdot)$ in (2)
requires specifying a valid cross-covariance function. We consider a linear
model of coregionalization (LMC),

$\boldsymbol{w}(\boldsymbol{\ell})=\sum_{j=1}^{k}\boldsymbol{\lambda}_{j}\omega_{j}(\boldsymbol{\ell})=\boldsymbol{\Lambda}\boldsymbol{\omega}(\boldsymbol{\ell})$
(4)

where
$\boldsymbol{\Lambda}=[\boldsymbol{\lambda}_{1},\ldots,\boldsymbol{\lambda}_{k}]$
is a $q\times k$ full (column) rank matrix with $(i,j)$th entry
$\lambda_{ij}$, and each $\omega_{j}(\boldsymbol{\ell})$ is a univariate
spatial process with correlation function
$\rho_{j}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=\rho(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime};\phi_{j})$.
The $k\leq q$ components of $\boldsymbol{\omega}(\boldsymbol{\ell})$ are
independent of each other (hence
$\text{cov}\\{\omega_{j}(\boldsymbol{\ell}),\omega_{h}(\boldsymbol{\ell}^{\prime})\\}=0$
whenever $h\neq j$). This implies that
$\boldsymbol{\omega}(\boldsymbol{\ell})$ is a multivariate spatial process
with cross-correlation
$\boldsymbol{\rho}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime};\boldsymbol{\varphi})$
indexed by $\boldsymbol{\varphi}=\\{\phi_{1},\dots,\phi_{k}\\}$ and
$\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime};\boldsymbol{\theta})=\boldsymbol{\Lambda}\boldsymbol{\rho}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})\boldsymbol{\Lambda}^{\top}$.
Under stationarity, if $\|\boldsymbol{\ell}-\boldsymbol{\ell}^{\prime}\|=0$,
then
$\boldsymbol{C}(\mathbf{0};\boldsymbol{\theta})=\boldsymbol{\Lambda}\boldsymbol{\rho}(\mathbf{0};\boldsymbol{\varphi})\boldsymbol{\Lambda}^{\top}=\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{\top}$
since $\boldsymbol{\rho}(0;\boldsymbol{\varphi})=\boldsymbol{I}_{k}$.
Therefore, $\boldsymbol{\Lambda}$ identifies with the Cholesky factorization
of $\boldsymbol{C}(\mathbf{0})$ and is lower-triangular with positive diagonal
entries (see. e.g., Finley et al., 2008; Zhang and Banerjee, 2020, and
references therein for Bayesian LMC models). Model (2) becomes
$\boldsymbol{y}=\boldsymbol{X}\bbeta+(I_{n}\otimes\boldsymbol{\Lambda})\boldsymbol{\omega}+\boldsymbol{\varepsilon}$.
Integrating out
$\boldsymbol{w}=(\boldsymbol{I}_{n}\otimes\boldsymbol{\Lambda})\boldsymbol{\omega}$
from the Gaussian likelihood fetches the density
$p(\boldsymbol{y}\,|\,\boldsymbol{X},\bbeta,\boldsymbol{\tau},\boldsymbol{\theta})=N(\boldsymbol{y};\boldsymbol{X}\bbeta,\boldsymbol{C}_{\boldsymbol{\theta}}+\boldsymbol{D}_{n})$,
where
$\boldsymbol{C}_{\boldsymbol{\theta}}=(\boldsymbol{I}_{n}\otimes\boldsymbol{\Lambda})\boldsymbol{\rho}_{\boldsymbol{\varphi}}(\boldsymbol{I}_{n}\otimes\boldsymbol{\Lambda}^{\top})$
and $\boldsymbol{\theta}=\\{\boldsymbol{\Lambda},\boldsymbol{\varphi}\\}$,
irrespective of whether we take $\boldsymbol{w}$ or $\boldsymbol{\omega}$ as
the latent process. Otherwise, different hierarchical model parametrizations
can be considered (Papaspiliopoulos et al., 2007).

A centered parametrization (Figure 1(a)) implies that $\boldsymbol{\theta}$
and $\boldsymbol{y}$ are conditionally independent given $\boldsymbol{w}$. A
noncentered parametrization (Figure 1(b)) considers $\boldsymbol{\Lambda}$ and
$\boldsymbol{\omega}$ independent a priori. A fully-noncentered
parametrization can be built by writing
$\boldsymbol{w}=(I_{n}\otimes\boldsymbol{\Lambda})\boldsymbol{L}_{\boldsymbol{\varphi}}\boldsymbol{u}$
where
$\boldsymbol{L}_{\boldsymbol{\varphi}}\boldsymbol{L}_{\boldsymbol{\varphi}}^{\top}=\boldsymbol{\rho}_{\boldsymbol{\varphi}}$
and $\boldsymbol{u}$ is white noise—the graphical model reads
$\boldsymbol{u}\to\boldsymbol{y}$ and $\boldsymbol{\theta}\to\boldsymbol{y}$.
The difference between centered and noncentered parametrizations is best
appreciated via a simple example. If $q=1$, in the centered model we may let
$\boldsymbol{\Lambda}=\sigma>0$ and assign prior $\sigma^{2}\sim
InvG(\sigma^{2};a,b)$. Then, the full conditional posterior is
$p(\sigma^{2}\,|\,\boldsymbol{w},\boldsymbol{y})=p(\sigma^{2}\,|\,\boldsymbol{w})=InvG(\sigma^{2};a^{*},b^{*})$
where $a^{*}=a+n/2$ and
$b^{*}=b+\frac{1}{2}\boldsymbol{w}^{\top}\boldsymbol{C}^{-1}_{\boldsymbol{\theta}}\boldsymbol{w}$.
In the noncentered models, $\boldsymbol{\Lambda}=\lambda>0$ can be assigned
e.g. a truncated Gaussian prior $N(\lambda;m,v)\boldsymbol{1}(\lambda>0)$,
resulting in a truncated Gaussian full conditional posterior, with variance
$v^{*}=(v^{-1}+\boldsymbol{\omega}^{\top}\boldsymbol{\omega}/\tau^{2})^{-1}$,
mean
$m^{*}=v^{*}\boldsymbol{\omega}^{\top}(\boldsymbol{y}-\boldsymbol{X}\bbeta)/\tau^{2}$,
and truncation below zero. Updating $\boldsymbol{\varphi}$ in all models
typically proceeds via Metropolis steps, hence the fully noncentered model
provides no additional advantage and we will not consider it further.

Specifically, we expect the relative efficiency of the noncentered models to
be poor when the ratio between the noise variance and process variance is
small (Papaspiliopoulos et al., 2007). As this commonly occurs in geospatial
settings, the centered parametrization might lead to better MCMC performance.
However, noncentered models are advantageous as they involve $k$ a-priori
independent processes rather than $q$ dependent ones, leading to more
convenient updates for $\boldsymbol{\Lambda}$, especially when setting $k<q$
as in a spatial factor model (Ren and Banerjee, 2013; Taylor-Rodriguez et al.,
2019; Zhang and Banerjee, 2020).

###  2.3 Weak identifiability of Matérn covariance functions

The LMC of Section 2.2 requires the specification of the $k$ correlation
functions $\phi_{j}(\cdot,\cdot)$. The Matérn model is a parsimonious but
flexible choice in this context. For a univariate spatial process
$\omega(\cdot)\sim GP(0,C_{\boldsymbol{\theta}})$, one assumes that
$\text{cov}\\{\omega(\boldsymbol{\ell}),\omega(\boldsymbol{\ell}^{\prime})\\}=\sigma^{2}\rho(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime};\\{\phi,\nu\\})$,
where

$\displaystyle\rho(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime};\\{\phi,\nu\\})$
$\displaystyle=\frac{2^{1-\nu}}{\Gamma(\nu)}\phi^{\nu}\|\boldsymbol{\ell}-\boldsymbol{\ell}^{\prime}\|^{\nu}K_{\nu}\left(\phi\|\boldsymbol{\ell}-\boldsymbol{\ell}^{\prime}\|\right),$
(5)

and $K_{\nu}$ is the modified Bessel function of the second kind of order
$\nu$. We know that $\boldsymbol{\theta}=\\{\sigma^{2},\phi,\nu\\}$ are weakly
identified (Banerjee et al., 2014). In fixed-domain asymptotics, only the
microergodic parameter $\\{\sigma^{2}\phi^{2\nu}\\}$ can be estimated
consistently (Stein, 1990; Ying, 1991; Zhang, 2004; Tang et al., 2020). In
practice, predictions at unknown locations are less sensitive to the
individual parameters’ value, since $\sigma^{2}\phi^{2\nu}$ matters more for
spatial interpolation. Nevertheless, weak identifiability is expected to
adversely impact MCMC mixing and convergence for $\\{\sigma^{2},\phi,\nu\\}$.
If one seeks to estimate the individual parameters and quantify uncertainty,
informative priors become necessary. It is not uncommon in these contexts to
fix some parameters to reasonable values, constrain the prior on a grid, or
perform greedy searches for optimal values through cross-validation (Shirota
et al., 2019; Banerjee, 2020).

###  2.4 Pitfalls of marginalization

Spatial dependence of the outcomes in (2) is modeled using the latent process.
Integrating out the latent effects from (3) yields the marginal or collapsed
Bayesian model

$p(\bbeta,\boldsymbol{\tau},\boldsymbol{\theta}\,|\,\boldsymbol{y})\propto
p(\boldsymbol{\tau},\boldsymbol{\theta})\times
N(\bbeta;\;\boldsymbol{\mu}_{\beta},\boldsymbol{V}_{\beta})\times
N(\boldsymbol{y};\;\boldsymbol{X}\bbeta,\boldsymbol{C}_{\boldsymbol{\theta}}+\boldsymbol{D}_{n})\;.$
(6)

Now we can sample the unknown parameters
$\\{\bbeta,\boldsymbol{\tau},\boldsymbol{\theta}\\}$ from (6). This is
appealing because we avoid sampling the high dimensional $\boldsymbol{w}$, a
step that is known to significantly encumber MCMC mixing and convergence for
parameters that are highly correlated to $\boldsymbol{w}$ a posteriori. Even
without $\boldsymbol{w}$ in (6), we can sample from
$p(\boldsymbol{w}\,|\,\boldsymbol{y})$ by drawing one $\boldsymbol{w}$ from
$p(\boldsymbol{w}\,|\,\bbeta,\boldsymbol{\tau},\boldsymbol{\theta},\boldsymbol{y})$,
which is Gaussian, for each sampled value of
$\\{\bbeta,\boldsymbol{\tau},\boldsymbol{\theta}\\}$ from (6).

Sampling from (6), however, entails
$(\boldsymbol{C}_{\boldsymbol{\theta}}+\boldsymbol{D}_{n})^{-1}$ and its
determinant, which becomes prohibitive in big data settings. To further
elucidate, the chief computational burden arises from the Cholesky
decomposition of
$\boldsymbol{K}_{\boldsymbol{\theta}}=\boldsymbol{C}_{\boldsymbol{\theta}}+\boldsymbol{D}_{n}$
or its inverse. Once the cholesky decomposition has been computed, say
$\texttt{chol}(\boldsymbol{K}_{\boldsymbol{\theta}})=\boldsymbol{L}\boldsymbol{F}\boldsymbol{L}^{\top}$,
where $\boldsymbol{L}=[l_{ii}]$ is unit lower triangular (with $l_{ii}=1$) and
$\boldsymbol{F}=\mbox{diag}[f_{ii}]$ is diagonal with $f_{ii}>0$, then

$\log
N(\boldsymbol{y};\;\boldsymbol{X}\bbeta,\boldsymbol{C}_{\boldsymbol{\theta}}+\boldsymbol{D}_{n})=\mbox{constant}-\frac{1}{2}\sum_{i=1}^{n}\log\left(f_{ii}\right)-\frac{1}{2}\sum_{i=1}^{n}\frac{u_{i}^{2}}{f_{ii}}\;,$
(7)

where $u_{i}$ is the $i$-th element of the vector $\boldsymbol{u}$ that solves
the lower-triangular system
$\boldsymbol{L}\boldsymbol{u}=\boldsymbol{y}-\boldsymbol{X}\bbeta$. Given the
Cholesky decomposition, computing (7) is cheap (of $O(nq)$). However, the
Cholesky decomposition itself involves $O(n^{3}q^{3})$ floating point
operations (flops) and $O(n^{2}q^{2})$ in storage complexity and, hence, is
impracticable for very large values of $nq$.

One option to improve computational speed is to replace
$\boldsymbol{w}(\cdot)$ with a low-rank process
$\widetilde{\boldsymbol{w}}(\cdot)$ with covariance function
$\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=\boldsymbol{C}_{\boldsymbol{\theta}}(\boldsymbol{\ell},{\cal
S})\boldsymbol{C}_{\boldsymbol{\theta}}^{-1}({\cal
S})\boldsymbol{C}_{\boldsymbol{\theta}}({\cal S},\boldsymbol{\ell}^{\prime})$,
where ${\cal S}$ is a set of $n_{k}$ knots (Banerjee et al., 2008). In this
case, the Sherman-Woodbury-Morrison (SWM) identity and the matrix determinant
lemma on the integrated covariance
$\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}+\boldsymbol{D}_{n}=\boldsymbol{C}_{\boldsymbol{\theta}}(\boldsymbol{\ell},{\cal
S})\boldsymbol{C}_{\boldsymbol{\theta}}^{-1}({\cal
S})\boldsymbol{C}_{\boldsymbol{\theta}}({\cal
S},\boldsymbol{\ell}^{\prime})+\boldsymbol{D}_{n}$ imply that one only needs
to compute $\boldsymbol{C}_{\boldsymbol{\theta}}^{-1}({\cal S})$ which is of
size $n_{k}\times n_{k}$ rather than $n\times n$. Directly computing the SWM
identity is not optimal and can be avoided using
$\texttt{chol}(\boldsymbol{C}_{\boldsymbol{\theta}}({\cal S}))$ (see Banerjee,
2017, for details).

Greater scalability is achieved when $\boldsymbol{w}(\cdot)$ is replaced by
any DAG-based process such as an NNGP or an MGP. Here,
$\boldsymbol{C}_{\boldsymbol{\theta}}$ is replaced by
$\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}$ such that
$\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}^{-1}$ is sparse.
Computations will ultimately depend on
$(\boldsymbol{C}_{\boldsymbol{\theta}}^{-1}+\boldsymbol{D}_{n}^{-1})^{-1}$ and
its determinant. This relies on sparse libraries (e.g. sparse Cholesky) which,
unfortunately, can impede iterative computations in big data analysis and may
be circumvented in some special circumstances (Peruzzi and Dunson, 2020).
Alternatively, one can directly model the responses as a sparse GP, in which
case the original cross-covariance function
$\boldsymbol{C}_{\boldsymbol{\theta}}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})$
is first replaced by
$\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})_{\boldsymbol{\theta}}+\delta_{\boldsymbol{\ell}=\boldsymbol{\ell}^{\prime}}\\{\boldsymbol{D}\\}$,
then the resulting full model’s precision
$\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}^{-1}$ is directly replaced
by a sparse counterpart
$\widetilde{\widetilde{\boldsymbol{C}}}_{\boldsymbol{\theta}}^{-1}$. Now one
does not need sparse Cholesky libraries as
$\widetilde{\widetilde{\boldsymbol{C}}}_{\boldsymbol{\theta}}^{-1}$ and its
determinant are computed directly from the DAG Finley et al. (2019). This
model does not include latent effects, but Vecchia approximations on the
latent effects tend to produce inference closer to the full GP model than the
response model (Katzfuss and Guinness, 2019). While marginalization avoids
some issues related to model parametrization, it does not address the weak
identifiability of the Matérn parameters in (5). Finally, the marginal and
response models are only feasible when all outcomes are Gaussian. Sampling
$\boldsymbol{w}$ a posteriori is, therefore, more generally applicable.

###  2.5 QMGPs at irregularly-spaced locations

Let us return to the Bayesian hierarchical model in (3) that explicitly
includes $\boldsymbol{w}$. While $\bbeta$ and $\boldsymbol{w}$ can be directly
updated using Gibbs steps from their (Gaussian) full conditional distributions
and $\boldsymbol{\tau}$ can also be updated directly if we assign independent
conjugate distributions (Inverse-Gamma) for each $\tau_{i}^{2}$, the cross-
covariance parameters $\boldsymbol{\theta}$ do not render closed-form full
conditional distributions. Therefore, we update $\boldsymbol{\theta}$ using
some stationarity-preserving MCMC method such as Metropolis random walk of
Hamiltonian Monte Carlo. These entail computing, up to a proportionality
constant, a density proportional to
$p(\boldsymbol{w}\,|\,\boldsymbol{\theta})=N(\boldsymbol{w};\;\mathbf{0},\boldsymbol{C}_{\boldsymbol{\theta}})\times
p(\boldsymbol{\theta})$, where $p(\boldsymbol{\theta})$ is any proper prior
distribution for $\boldsymbol{\theta}$. Unfortunately, once again, we are
faced with computing the inverse and determinant of a high-dimensional
covariance matrix.

Scalable methods accrue significant benefits in computing
$p(\boldsymbol{w}\,|\,\boldsymbol{\theta})$. In particular, a “cubic” meshed
GP (QMGP; Peruzzi et al., 2020) produces a reversal in the balance of
computational costs: evaluations of
$p(\boldsymbol{w}\,|\,\boldsymbol{\theta})$ become faster than sampling
$\boldsymbol{w}$—usually they are about 10 times slower. In order to achieve
such speedups, one requires that the so-called reference set ${\cal S}$ be
defined on a regular grid. ${\cal S}$ is akin to the set of “knots” or
“sensors” in low rank methods, whose choice is not a trivial matter, see e.g.
Krause et al. (2008) or Guhaniyogi et al. (2011). Standard QMGPs are
associated with MCMC efficiency losses when a gridded reference set does not
overlap with the measurement locations. In this case the size of the latent
effects might greatly exceed the sample size, with negative impacts on MCMC
mixing and convergence of the centered parameters. We clarify this point after
building the QMGP that we subsequently use in the remainder of the article.

For any given finite set ${\cal L}$ of spatial locations and a fixed reference
set ${\cal S}$, let ${\cal U}={\cal L}\setminus{\cal S}$ denote the non-
reference locations. We partition the domain using an axis-parallel
tessellation $\boldsymbol{T}=\\{D_{j}\\}_{j=1}^{M}$ such that ${\cal
D}=\cup_{j=1}^{M}D_{j}$, where $M=\prod_{i=1}^{d}M_{i}$. We then partition the
reference and non-reference sets into ${\cal S}=\cup_{i=1}^{M}S_{i}$ and
${\cal U}=\cup_{i=1}^{M}U_{i}$, where $S_{i}={\cal S}\cap D_{i}$ and
$U_{i}={\cal U}\cap D_{i}$. Let the DAG ${\cal
G}=\\{\boldsymbol{V},\boldsymbol{E}\\}$ have nodes
$\boldsymbol{V}=\boldsymbol{A}\cup\boldsymbol{B}$ where $\boldsymbol{A}$ and
$\boldsymbol{B}$ are the $2M$ reference and non-reference nodes, respectively.
$\boldsymbol{T}$ and ${\cal G}$ are linked via $\eta:{\cal
D}\to\boldsymbol{V}$ which maps any reference location to a node in
$\boldsymbol{A}$ and any non-reference location to a node in $\boldsymbol{B}$.
The parent set of nodes for $\boldsymbol{v}\in\boldsymbol{V}$ is denoted as
$\text{Pa}[{\boldsymbol{v}}]$, and
$\boldsymbol{E}=\\{\text{Pa}[{\boldsymbol{v}}]\\}_{\boldsymbol{v}\in\boldsymbol{V}}$.
By construction, we let $\text{Pa}[{\boldsymbol{v}}]\subset\boldsymbol{A}$ if
$\boldsymbol{v}\in\boldsymbol{B}$, i.e. non-reference nodes have no children.
Furthermore, in the cubic MGP’s specification of ${\cal G}$ we have
$|\text{Pa}[{\boldsymbol{v}}]|\leq d$ and
$|\text{Ch}[{\boldsymbol{v}}]\cap\boldsymbol{A}|\leq d$ whenever
$\boldsymbol{v}\in\boldsymbol{A}$, where $d$ is the dimension of ${\cal D}$.
For non-reference nodes, let $\eta^{-1}(U_{i})=\boldsymbol{b}_{i}$ and
$\eta^{-1}(S_{i})=\boldsymbol{a}_{i}$. Then, we specify
$\text{Pa}[{\boldsymbol{b}_{i}}]=\\{\boldsymbol{a}_{i}\\}$—non-reference
locations’ parent set comprises reference locations in the same domain
partition. These specifications for ${\cal G}$ imposes a coloring on ${\cal
G}$ enabling parallel sampling of the latent random effects (see, e.g.,
Gonzalez et al., 2011), independently of observed locations ${\cal T}$ and
$d$. At ${\cal L}$, the new process has density $p^{*}$,

$\displaystyle p^{*}(\boldsymbol{w}_{{\cal L}})$
$\displaystyle=\prod_{i=1}^{M}p(\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[{i}]})\prod_{\boldsymbol{\ell}\in
U_{i}}p(\boldsymbol{w}(\boldsymbol{\ell})\,|\,\boldsymbol{w}_{i})\;,$ (8)

where
$\boldsymbol{w}_{i}=\boldsymbol{w}(\eta^{-1}(\boldsymbol{v}_{i}))=\boldsymbol{w}(S_{i})=(\boldsymbol{w}(\boldsymbol{\ell}_{1})^{\top},\dots,\boldsymbol{w}(\boldsymbol{\ell}_{n_{S_{i}}})^{\top})^{\top}$,
$\boldsymbol{w}_{[{i}]}=\boldsymbol{w}(\eta^{-1}(\text{Pa}[{\boldsymbol{v}_{i}}]))$,
using $\eta^{-1}$ to find all reference locations mapped to the same node. We
assume conditional independence at non-reference locations in this article. A
Gaussian base process leads to (8) as a product of multivariate Gaussian
densities; $p^{*}(\boldsymbol{w}_{{\cal L}})$ is then Gaussian with a sparse
precision matrix – refer to Peruzzi et al. (2020) for more details.

For an observed location $\boldsymbol{\ell}\in{\cal T}\cap D_{i}$, the
standard MGP regression model is

$\begin{split}\boldsymbol{y}(\boldsymbol{\ell})&=\boldsymbol{X}(\boldsymbol{\ell})^{\top}\bbeta+\boldsymbol{w}(\boldsymbol{\ell})+\boldsymbol{\varepsilon}(\boldsymbol{\ell})\;;\quad\boldsymbol{\varepsilon}(\boldsymbol{\ell})\sim
N(\mathbf{0},\boldsymbol{D})\;;\\\
\boldsymbol{w}(\boldsymbol{\ell})&=\boldsymbol{H}_{\boldsymbol{\ell}}\boldsymbol{w}_{[\boldsymbol{\ell}]}+\boldsymbol{\xi}(\boldsymbol{\ell})\;;\quad\boldsymbol{\xi}(\boldsymbol{\ell})\sim
N(\mathbf{0},\boldsymbol{R}_{\boldsymbol{\ell}})\;,\end{split}$ (9)

where, denoting
$\boldsymbol{C}_{\boldsymbol{\ell}}=\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell})$,
the set of parent locations as
$B=\eta^{-1}(\text{Pa}[{\eta(\boldsymbol{\ell})}])$ and
$\boldsymbol{C}_{[\boldsymbol{\ell}]}=\boldsymbol{C}(B,B)$, we have
$\boldsymbol{H}_{\boldsymbol{\ell}}=\boldsymbol{C}_{\boldsymbol{\ell},[\boldsymbol{\ell}]}\boldsymbol{C}_{[\boldsymbol{\ell}]}^{-1}$
and
$\boldsymbol{R}_{\boldsymbol{\ell}}=\boldsymbol{C}_{\boldsymbol{\ell}}-\boldsymbol{C}_{\boldsymbol{\ell},[\boldsymbol{\ell}]}\boldsymbol{C}_{[\boldsymbol{\ell}]}^{-1}\boldsymbol{C}_{[\boldsymbol{\ell}],\boldsymbol{\ell}}$.
Since evaluating $p(\boldsymbol{w}\,|\,\boldsymbol{\theta})$ is required in
posterior sampling of $\boldsymbol{\theta}$, the primary computing cost
depends on $\boldsymbol{C}_{[\boldsymbol{\ell}]}^{-1}$. The number of such
inverses and their sizes depend on ${\cal T}$, ${\cal S}$, and ${\cal U}$, for
a given ${\cal G}$ and $\boldsymbol{T}$. If ${\cal T}$ are gridded, we let
${\cal S}\supset{\cal T}$ be the observed grid plus the missing grid
locations, and ${\cal U}=\emptyset$. Only $O(d)$ larger matrix inverses are
needed, rather than $M$, and $\boldsymbol{w}$ is sampled at $|{\cal S}|$
locations.

In fact, with stationary cross-covariances, for a given node
$\boldsymbol{v}_{i}$ one can find “prototypes” $\boldsymbol{v}_{i^{\\#}}$ and
$\boldsymbol{v}_{i^{*}}$ such that
$\boldsymbol{C}_{i}=\boldsymbol{C}_{i^{\\#}}$ and
$\boldsymbol{C}_{[i]}=\boldsymbol{C}_{[i^{*}]}$, and the number of prototype
nodes is $O(d)$. One computes the inverses for the prototypes and copies the
result across the whole domain. If ${\cal T}$ are not gridded, then we still
let ${\cal S}$ be a grid, again leading to $O(d)$ large matrix solvers.
However, now ${\cal T}={\cal U}$ and we sample $\boldsymbol{w}$ at $|{\cal
S}|+|{\cal U}|$ locations. Since $n_{{\cal S}}=|{\cal S}|\approx n$ to avoid
over-smoothing the spatial surface, we would be using about twice as many
latent variables as we have observations. Updates of the centered parameters
proceed by evaluating $p(\boldsymbol{w}_{{\cal
S}}\,|\,\boldsymbol{\theta})\times p(\boldsymbol{w}_{{\cal
U}}\,|\,\boldsymbol{w}_{{\cal S}},\boldsymbol{\theta})$, which implies that
the move from $\boldsymbol{\theta}$ to $\boldsymbol{\theta}^{\prime}$ must
agree with the current high-dimensional sample of the latent process at both
reference and non-reference locations. Consequently, MCMC convergence will be
slow. Clearly, the computational speedup of gridded reference sets may not
compensate for such inefficiencies.

##  3 Grid, Parametrize, Split

We now generalize the computational advantages of QMGPs to non-gridded data,
and rebuild the linear model of coregionalization to improve the overall
sampling efficiency. We assume, without loss of generality, that the observed
locations ${\cal T}$ are irregularly spaced. We target raw speed in the first
subsection, efficiency in the next.

###  3.1 Reference set on a Grid

We build the reference set ${\cal S}$ as a grid, regardless of ${\cal T}$. If
${\cal D}=[0,1]^{d}$, we may take ${\cal
S}=\left\\{\left(\frac{j_{1}}{N_{1}},\dots,\frac{j_{d}}{N_{d}}\right):(j_{1},\dots,j_{d})\in\\{1,\dots,N_{d}\\}^{d}\right\\}$,
so the number of reference locations is $n_{S}=\prod_{j=1}^{d}N_{j}$. From (9)
we can partly marginalize the latent process at the observed locations. For a
location $\boldsymbol{\ell}\in D_{i}$ and conditioning on realizations of the
latent process at ${\cal S}$, we have

$\boldsymbol{y}(\boldsymbol{\ell})=\boldsymbol{X}(\boldsymbol{\ell})^{\top}\bbeta+\boldsymbol{H}_{\boldsymbol{\ell}}\boldsymbol{w}_{[\boldsymbol{\ell}]}+\boldsymbol{\varepsilon}^{\prime}(\boldsymbol{\ell})\;;\quad\boldsymbol{\varepsilon}^{\prime}(\boldsymbol{\ell})\sim
N(\mathbf{0},\boldsymbol{D}+\boldsymbol{R}_{\boldsymbol{\ell}})\;,$ (10)

where $\boldsymbol{w}_{[\boldsymbol{\ell}]}$ are process realizations at
locations $\eta^{-1}(\text{Pa}[{\boldsymbol{b}_{i}}])$ since
$\boldsymbol{\ell}\in D_{i}$ and, hence,
$\eta(\boldsymbol{\ell})=\boldsymbol{b}_{i}$,
$\boldsymbol{H}_{\boldsymbol{\ell}}$ and $\boldsymbol{R}_{\boldsymbol{\ell}}$
are $q\times n_{S_{i}}$ and $q\times q$ matrices, respectively. For a given
$\boldsymbol{w}_{{\cal S}}$, if $\boldsymbol{\ell}\in{\cal S}$ then
$\boldsymbol{R}_{\boldsymbol{\ell}}=\mathbf{0}$ and
$\boldsymbol{H}_{\boldsymbol{\ell}}\boldsymbol{w}_{[\boldsymbol{\ell}]}=\boldsymbol{w}_{\boldsymbol{\ell}}$.
In other words, (10) generalizes QMGPs to nongridded data without sampling at
too many locations. Furthermore, while no issue arise if ${\cal T}$ and ${\cal
S}$ partly overlap, in general they are exclusive.

Partial marginalization of the latent process removes the need to sample at
non-reference locations. However, unlike full marginalization, the advantages
induced by having ${\cal S}$ on a grid are maintained. Moreover, choosing a
grid ${\cal S}$ such that $n_{{\cal S}}\ll n$ leads to interpreting (10) as a
“meshed” extension of the bias-adjusted predictive process (MPP; Finley et
al., 2009). While computational considerations make it advantageous to choose
${\cal S}$ on a grid, (10) does not preclude alternatives. For example,
typical ${\cal T}={\cal S}$ implies prediction locations do are not in the
reference set, but here—assuming $\boldsymbol{X}(\boldsymbol{\ell})$ are
available—we can choose ${\cal S}$ as the prediction set. Then, sampling
$\boldsymbol{w}_{{\cal S}}$ yields the desired predictions. This feature can
be used to produce model-based heat maps (i.e. raster images) of the latent
process at any resolution. Finally, while partial marginalization is available
for Gaussian first stage models, we expect analogous strategies to be
similarly advantageous otherwise, taking advantage of conditional independence
of all outcomes on the latent process.

###  3.2 Parametrize a Matérn LMC, Split its variance

This subsection does not assume $\boldsymbol{w}$ is a QMGP. After letting
$\widetilde{\boldsymbol{\omega}}(\boldsymbol{\ell})=(\widetilde{\omega}_{1}(\boldsymbol{\ell}),\dots,\widetilde{\omega}_{k}(\boldsymbol{\ell}))^{\top}$
and $\widetilde{\omega}_{j}(\boldsymbol{\ell})\sim
GP(0,\widetilde{\boldsymbol{C}}_{j})$, where

$\displaystyle\widetilde{\boldsymbol{C}}_{j}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})$
$\displaystyle=\frac{\rho(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime};\\{\phi_{j},\nu_{j}\\})}{\phi_{j}^{2\nu_{j}}},$
(11)

we write the LMC as
$\boldsymbol{w}(\boldsymbol{\ell})=\boldsymbol{\Lambda}\boldsymbol{\omega}(\boldsymbol{\ell})=\boldsymbol{\Lambda}\boldsymbol{Q}\boldsymbol{Q}^{-1}\boldsymbol{\omega}(\boldsymbol{\ell})=\widetilde{\boldsymbol{\Lambda}}\widetilde{\boldsymbol{\omega}}(\boldsymbol{\ell})$,
where
$\boldsymbol{Q}=\text{diag}\\{\phi_{1}^{\nu_{1}},\dots,\phi_{k}^{\nu_{k}}\\}$.
$\widetilde{\boldsymbol{\Lambda}}=\boldsymbol{\Lambda}\boldsymbol{Q}$ is thus
a lower-triangular rescaling matrix whose $j$th diagonal element is
$\lambda_{jj}\phi_{j}^{\nu_{j}}$, which is the square root of the microergodic
parameter (Stein, 1999)—its better identifiability impacts MCMC efficiency
positively. The immediate consequence of this change is that posterior samples
of $\boldsymbol{\Lambda}$ will only be obtained via post-processing of MCMC
output—we will instead directly sample $\widetilde{\boldsymbol{\Lambda}}$ and
$\widetilde{\boldsymbol{\omega}}$.

The LMC is typically estimated as a noncentered parametrization, which as we
mentioned is associated with poor mixing when the nuggets are small. We
address this issue by parameter expansion (or splitting) of
$\lambda_{jj}^{2}$: denoting
$\boldsymbol{\Psi}=\\{\lambda_{11},\dots,\lambda_{kk}\\}$, we find
$\boldsymbol{L}=\boldsymbol{\Lambda}\boldsymbol{\Psi}^{-1}$ as a lower
triangular matrix with units along the diagonal. Letting
$\boldsymbol{J}=\text{diag}\\{\sigma_{1}^{-1},\dots,\sigma_{k}^{-1}\\}$, we
write the LMC as
$\boldsymbol{w}(\boldsymbol{\ell})=\boldsymbol{\Lambda}\boldsymbol{\omega}(\boldsymbol{\ell})=\boldsymbol{L}\boldsymbol{\Psi}\boldsymbol{\omega}(\boldsymbol{\ell})$.
The split follows by writing
$\boldsymbol{w}(\boldsymbol{\ell})=\boldsymbol{\Lambda}\boldsymbol{J}\boldsymbol{J}^{-1}\boldsymbol{\omega}(\boldsymbol{\ell})=\overline{\boldsymbol{\Lambda}}\overline{\boldsymbol{\omega}}(\boldsymbol{\ell})$,
where
$\overline{\boldsymbol{\Lambda}}=\boldsymbol{L}\boldsymbol{\Psi}\boldsymbol{J}$,
$\overline{\omega}_{j}(\boldsymbol{\ell})\sim GP(0,\sigma^{2}_{j}\rho_{j})$
and $\rho_{j}$ is as above. In practice, we sample
$\overline{\boldsymbol{\omega}}$ alongside $\overline{\boldsymbol{\Lambda}}$,
whose $j$th diagonal element is
$\overline{\lambda}_{jj}=\lambda_{jj}/\sigma_{j}$. The split entails an
overparametrization of the variances as we assign a prior on
$\overline{\boldsymbol{\Lambda}}$ which is independent of the prior on
$\boldsymbol{J}$. If $q=1$, we may assign a truncated Gaussian prior on
$\overline{\lambda}$ centered at zero with truncation below zero, and an
inverse Gamma on $\sigma^{2}$, to get a half-t for
$\lambda=\overline{\lambda}\sigma$ (Gelman, 2006). In the context of Bayesian
factor models, a similar parameter expansion was introduced in Ghosh and
Dunson (2009).

##  4 Bayesian GriPS regression

We build the model for $\boldsymbol{w}$ by combining parametrization and
splitting strategies:

$\begin{split}\boldsymbol{w}(\boldsymbol{\ell})&=\boldsymbol{\Lambda}\boldsymbol{\omega}(\boldsymbol{\ell})=\boldsymbol{\Lambda}\boldsymbol{Q}\boldsymbol{J}\boldsymbol{J}^{-1}\boldsymbol{Q}^{-1}\boldsymbol{\omega}(\boldsymbol{\ell})=\boldsymbol{{\cal
A}}\boldsymbol{r}(\boldsymbol{\ell}),\end{split}$ (12)

where $\boldsymbol{{\cal A}}=\boldsymbol{\Lambda}\boldsymbol{Q}\boldsymbol{J}$
is a $q\times k$ matrix with zeros in the upper triangular part, the
$k$-variate process is $\boldsymbol{r}(\cdot)\sim
QMGP(\mathbf{0},\boldsymbol{K}(\cdot))$ such that
$\boldsymbol{r}(\boldsymbol{\ell})=(r_{1}(\boldsymbol{\ell}),\dots,r_{k}(\boldsymbol{\ell}))^{\top}$
and, for each $j=1,\dots,k$, $r_{j}(\boldsymbol{\ell})\sim
QMGP\left(0,K_{j}(\cdot)\right)$ where
$K_{j}(\cdot)=\sigma^{2}_{j}\frac{\rho_{j}(\cdot)}{\phi_{j}^{2\nu}}$ with
$\rho_{j}(\cdot)$ as the Matérn correlation in (5). All independent QMGPs are
built on the same domain partition $\boldsymbol{T}$ and the same graph ${\cal
G}$.

$\boldsymbol{y}$$\boldsymbol{{\cal
A}}$$\boldsymbol{r}$$\boldsymbol{\varphi}$$\bbeta,\boldsymbol{\tau}$$\boldsymbol{{\cal
A}}=\boldsymbol{\Lambda}\boldsymbol{Q}\boldsymbol{J}=\boldsymbol{\Lambda}\left[\begin{smallmatrix}\phi_{1}^{\nu_{1}}&&\\\
&\ddots&\\\
&&\phi_{k}^{\nu_{k}}\end{smallmatrix}\right]\left[\begin{smallmatrix}\sigma_{1}^{-1}&&\\\
&\ddots&\\\
&&\sigma_{k}^{-1}\end{smallmatrix}\right]$$\boldsymbol{\varphi}=\\{\sigma^{2}_{j},\phi_{j},\nu_{j}\\}_{j=1}^{k}$
Figure 2: DAG of the Bayesian spatial regression model using a PS
parametrization. The relationship with a standard parametrization is
summarised on the right.

The GriPS spatial regression model then becomes

$\begin{split}\boldsymbol{y}(\boldsymbol{\ell})&=\boldsymbol{X}(\boldsymbol{\ell})^{\top}\bbeta+\boldsymbol{{\cal
Z}}_{\boldsymbol{\ell}}\boldsymbol{r}_{[\boldsymbol{\ell}]}+\boldsymbol{\varepsilon}^{\prime}(\boldsymbol{\ell})\;;\quad\boldsymbol{\varepsilon}^{\prime}(\boldsymbol{\ell})\sim
N(\mathbf{0},\boldsymbol{D}+\boldsymbol{\Sigma}_{\boldsymbol{\ell}}),\end{split}$
(13)

where $\boldsymbol{{\cal Z}}_{\boldsymbol{\ell}}=\boldsymbol{{\cal
A}}\boldsymbol{H}_{\boldsymbol{\ell}}$,
$\boldsymbol{\Sigma}_{\boldsymbol{\ell}}=\boldsymbol{{\cal
A}}\boldsymbol{R}_{\boldsymbol{\ell}}\boldsymbol{{\cal A}}^{\top}$. Writing
(13) for locations ${\cal T}$, we get

$\begin{split}\boldsymbol{y}&=\boldsymbol{X}\bbeta+\boldsymbol{{\cal
Z}}\boldsymbol{r}+\boldsymbol{\varepsilon}^{\prime}\;;\quad\boldsymbol{\varepsilon}^{\prime}\sim
N(\mathbf{0},\boldsymbol{D}_{n}+\boldsymbol{\Sigma})\;,\end{split}$ (14)

where $\otimes$ denotes the Kronecker product, $\boldsymbol{{\cal
Z}}=(\boldsymbol{I}_{n}\otimes\boldsymbol{{\cal A}})\boldsymbol{{\cal H}}$,
and $\boldsymbol{\Sigma}=(\boldsymbol{I}_{n}\otimes\boldsymbol{{\cal
A}})\boldsymbol{{\cal R}}(I_{n}\otimes\boldsymbol{{\cal A}}^{\top})$. In
addition to the terms already appearing in (2), we now explicitly have
$\boldsymbol{r}$ representing the $nk\times 1$ vector of latent effects at the
gridded reference set ${\cal S}$. Furthermore, $\boldsymbol{{\cal H}}$ is a
sparse $nk\times n_{{\cal S}}k$ matrix whose $i$th row block is filled with
$\boldsymbol{H}_{\boldsymbol{\ell}_{i}}$, with zeroes at columns corresponding
to locations in ${\cal S}$ not in the parent set of $\boldsymbol{\ell}_{i}$.
Since the $k$ margins of the $\boldsymbol{r}$ process are independent, the
permutation of the columns of $\boldsymbol{H}_{\boldsymbol{\ell}_{i}}$ that
arranges it according to the process margins is a block-diagonal matrix whose
$j$th (of the total $k$) row block is
$K_{j}(\boldsymbol{\ell},[\boldsymbol{\ell}])K_{j}^{-1}([\boldsymbol{\ell}])$.
Similarly, $\boldsymbol{{\cal R}}$ is a diagonal $nk\times nk$ matrix whose
$i$th block is $\boldsymbol{R}_{\boldsymbol{\ell}_{i}}$, itself diagonal with
$j$th diagonal element
$K_{j}(\boldsymbol{\ell},\boldsymbol{\ell})-K_{j}(\boldsymbol{\ell},[\boldsymbol{\ell}])K_{j}^{-1}([\boldsymbol{\ell}])K_{j}([\boldsymbol{\ell}],\boldsymbol{\ell})$.
We find $\boldsymbol{\Sigma}$ as a block-diagonal matrix of dimension
$nq\times nq$ with $q\times q$ blocks—the $i$th block is $\boldsymbol{{\cal
A}}\boldsymbol{R}_{\boldsymbol{\ell}_{i}}\boldsymbol{{\cal A}}^{\top}$.

We complete the Bayesian model by assigning multivariate Gaussian prior
distributions to $\bbeta$ and the lower triangular elements of
$\boldsymbol{{\cal A}}$, and independent zero-mean Gaussian distributions
truncated below zero for the diagonal elements of $\boldsymbol{{\cal A}}$. We
also assign independent Inverse Gamma priors for each $\sigma^{2}_{j}$,
$j=1,\dots,k$ and each $\tau^{2}_{i}$, $i=1,\dots,q$. Figure 2 summarises the
PS parametrization—GriPS regressions are noncentered on $\boldsymbol{{\cal
A}}$ and centered on
$\boldsymbol{\varphi}=\\{\sigma^{2}_{j},\phi_{j},\nu_{j}\\}_{j=1}^{k}$.
Ultimately, we are interested in $\boldsymbol{\Lambda}=\boldsymbol{{\cal
A}}\boldsymbol{J}^{-1}\boldsymbol{Q}^{-1}$, meaning that a PS parametrization
is simultaneously centered and noncentered on $\boldsymbol{\Lambda}$ since
$\boldsymbol{J}$ and $\boldsymbol{Q}$ include the centered parameters. Next,
we outline the sampling algorithms for estimation and prediction.

###  4.1 Estimation and prediction

Updating $\boldsymbol{r}$ proceeds following ${\cal G}$. Taking reference node
$\boldsymbol{v}_{i}\in\boldsymbol{A}$, the set $T_{i}$ are locations having
$S_{i}$ as parent set, hence we can write
$\boldsymbol{y}_{T_{i}}=\boldsymbol{X}_{T_{i}}\bbeta+\boldsymbol{{\cal
Z}}_{T_{i}}\boldsymbol{r}_{i}+\boldsymbol{\varepsilon}_{T_{i}}$ and
$\boldsymbol{\varepsilon}_{T_{i}}\sim
N(\mathbf{0},\boldsymbol{D}_{T_{i}}+\boldsymbol{\Sigma}_{T_{i}})$. Denote
$\widetilde{\boldsymbol{y}}_{T_{i}}=\boldsymbol{I}_{T_{i}}(\boldsymbol{y}_{T_{i}}-\boldsymbol{X}_{T_{i}}\bbeta)$
and
$\widetilde{\boldsymbol{D}}_{T_{i}}=\boldsymbol{I}_{T_{i}}(\boldsymbol{D}_{T_{i}}+\boldsymbol{\Sigma}_{T_{i}})\boldsymbol{I}_{T_{i}}$,
where $\boldsymbol{I}_{T_{i}}$ is a diagonal matrix with zeros at elements
corresponding to missing outcomes in $T_{i}$, ones elsewhere. Calling
$\text{Ch}[{\boldsymbol{v}_{i}}]=\\{\boldsymbol{v}_{j}:\boldsymbol{v}_{i}\in\text{Pa}[{\boldsymbol{v}_{j}}]\\}$,
we partition
$\boldsymbol{H}_{j}=[\boldsymbol{H}_{i}:\boldsymbol{H}_{[j]\setminus i}]$ for
each $\boldsymbol{v}_{j}\in\text{Ch}[{\boldsymbol{v}_{i}}]$, and similarly
$\boldsymbol{r}_{[j]}=(\boldsymbol{r}_{i}^{\top},\boldsymbol{r}_{[j]\setminus
i}^{\top})^{\top}$, where the $[j]\setminus i$ notation refers the set
$\text{Pa}[{\boldsymbol{v}_{j}}]\setminus\\{\boldsymbol{v}_{i}\\}$. Then, the
full conditional distribution of $\boldsymbol{r}_{i}$ is Gaussian, with mean
$\boldsymbol{g}_{i}$ and covariance $\boldsymbol{G}_{i}$ such that:

$\displaystyle\boldsymbol{G}_{i}^{-1}$ $\displaystyle=\boldsymbol{{\cal
Z}}_{T_{i}}^{\top}\widetilde{\boldsymbol{D}}_{T_{i}}^{-1}\boldsymbol{{\cal
Z}}_{T_{i}}+\boldsymbol{R}_{i}^{-1}+\sum\limits_{j\in\text{Ch}[{\boldsymbol{v}_{i}}]}\boldsymbol{H}_{[j]\setminus
i}^{\top}\boldsymbol{R}_{j}^{-1}\boldsymbol{H}_{[j]\setminus i}$
$\displaystyle\boldsymbol{G}_{i}\boldsymbol{g}_{i}$
$\displaystyle=\boldsymbol{R}_{i}^{-1}\boldsymbol{H}_{i}\boldsymbol{r}_{[i]}+\boldsymbol{{\cal
Z}}_{T_{i}}^{\top}\widetilde{\boldsymbol{D}}_{T_{i}}^{-1}\widetilde{\boldsymbol{y}}_{T_{i}}+\sum\limits_{j\in\text{Ch}[{\boldsymbol{v}_{i}}]}\boldsymbol{H}_{[j]\setminus
i}^{\top}\boldsymbol{R}_{j}^{-1}\boldsymbol{r}_{[j]\setminus i}.$

Updating $\bbeta$ proceeds via its full conditional distribution, which is
$N(\bbeta;\boldsymbol{m}_{\bbeta}^{(n)},\boldsymbol{V}_{\bbeta}^{(n)})$, where
$\boldsymbol{V}_{\bbeta}^{{}^{(n)}-1}=\boldsymbol{V}_{\bbeta}^{-1}+\boldsymbol{X}^{\top}(\boldsymbol{D}_{n}+\boldsymbol{\Sigma})^{-1}\boldsymbol{X}$
and
$\boldsymbol{V}_{\bbeta}^{(n)}\boldsymbol{m}_{\bbeta}^{(n)}=\boldsymbol{V}_{\bbeta}^{-1}\boldsymbol{m}_{\bbeta}+\boldsymbol{X}^{\top}(\boldsymbol{D}_{n}+\boldsymbol{\Sigma})^{-1}(\boldsymbol{y}-\boldsymbol{{\cal
Z}}\boldsymbol{r})$. Updates for other unknowns may proceed via Metropolis
steps: for $\boldsymbol{{\cal A}}$ and $\boldsymbol{\tau}$, the centered
parametrization results in being target density being the prior
$\pi(\boldsymbol{{\cal A}})\pi(\boldsymbol{\tau})$, multiplied by the
following product of multivariate Gaussians

$\prod\limits_{i\in{\cal
T}}\left.N(\boldsymbol{y}_{\boldsymbol{\ell}_{i}};\boldsymbol{X}_{\boldsymbol{\ell}_{i}}\bbeta+\boldsymbol{{\cal
Z}}_{\boldsymbol{\ell}_{i}}\boldsymbol{r}_{\boldsymbol{\ell}_{i}},\boldsymbol{D}+\boldsymbol{\Sigma}_{\boldsymbol{\ell}_{i}})\right|_{\text{obs}},$
(15)

where we restrict each term to the observed outcomes, i.e. each term is an
$l$-dimensional Gaussian corresponding to the $l\leq q$ observed outcomes at
spatial location $\boldsymbol{\ell}_{i}$.

For $\boldsymbol{\varphi}=\\{\sigma^{2}_{j},\phi_{j},\nu_{j}\\}_{j=1}^{k}$,
the Metropolis update following the centered parametrization targets the
density obtained by multiplying the prior $\pi(\boldsymbol{\varphi})$ to
$\prod_{i}^{M}N(\boldsymbol{r}_{i};\boldsymbol{H}_{i}\boldsymbol{r}_{[i]},\boldsymbol{R}_{i})$,
which translates (8) for the GriPS coregionalization model. In our
implementation, we use the robust adaptive Metropolis algorithm of Vihola
(2012) to target an acceptance rate of about 0.23. We note that it is
computationally easier to update $\boldsymbol{{\cal A}}$, $\boldsymbol{\tau}$
and $\boldsymbol{\sigma}^{2}$ than $\boldsymbol{\phi}$ and $\boldsymbol{\nu}$
as the latter two require recomputing the correlations. Finally, if ${\cal T}$
is an incomplete grid and we choose ${\cal S}\supset{\cal T}$ as the grid
which “fills the gaps” in ${\cal T}$, then the regression error is
$\boldsymbol{\varepsilon}\sim N(\mathbf{0},\boldsymbol{D}_{n})$. In that case,
with the prior distributions we have mentioned earlier, $\boldsymbol{{\cal
A}}$ and $\boldsymbol{\tau}$ can be updated via their full conditionals—refer
to the supplementary material for details. If ${\cal S}\neq{\cal T}$, then
those full conditionals are the correct updates for posterior sampling of
model $\boldsymbol{y}=\boldsymbol{X}\bbeta+\boldsymbol{{\cal
Z}}\boldsymbol{r}+\widetilde{\boldsymbol{\varepsilon}}$ where
$\widetilde{\boldsymbol{\varepsilon}}\sim N(0,\boldsymbol{D}_{n})$, which is
not (14) as $\boldsymbol{\Sigma}$ is omitted and is instead akin to a meshed
version of the “unmodified” predictive process of Banerjee et al. (2008); we
expect the performance of the two models to be very similar when ${\cal S}$ is
large and ${\cal S}\supset\widetilde{{\cal S}}\approx{\cal T}$. All
applications in this article use Metropolis updates for (14) on nongridded
data.

###  4.2 Applications on synthetic data

We illustrate the performance of GriPS in estimation and prediction in
univariate spatial regressions using simulated data. Each generated data set
includes a realization of a Gaussian process with Matérn covariance (5),
fixing the smoothness $\nu$ to its known value (either $\nu=0.5$ or
$\nu=1.5$). The outcomes are observed at irregularly-spaced locations in
$[0,1]^{2}$. The total number of observed locations is $n=5000$. Our main goal
is the efficient estimation of the posterior distribution of
$\sigma^{2},\phi,\tau^{2}$. Predictive performance is used to validate our
methodology. We use a grid of $n_{\text{out}}=10000$ spatial locations as our
test set. We set ${\cal S}$ in GriPS models to overlap with the test set—we
are thus sampling the spatial process at twice as many locations as there are
observations. For $\nu=0.5$, we generate 50 datasets for each combination of
$\sigma^{2}\in\\{1,5\\}$ and $\phi\in\\{0.5,5\\}$, and similarly for $\nu=1.5$
but with $\phi\in\\{5,25\\}$, for a total of 400 datasets each of size
$n_{\text{all}}=15000$.

|  | $\nu=0.5$ | $\nu=1.5$
---|---|---|---
Efficiency | ESS | ESS/s | Relative | ESS | ESS/s | Relative
QMGP GriPS | $\phi$ | 245.46 | 12.78 | 52.11 | 125.40 | 6.49 | 50.59
$\sigma^{2}$ | 248.68 | 12.95 | 49.04 | 224.83 | 11.68 | 97.44
$\tau^{2}$ | 116.04 | 6.04 | 1.38 | 231.46 | 12.03 | 6.68
QMGP latent | $\phi$ | 4.46 | 0.19 | 0.78 | 4.84 | 0.28 | 2.16
$\sigma^{2}$ | 5.88 | 0.25 | 0.95 | 7.34 | 0.42 | 3.51
$\tau^{2}$ | 206.54 | 8.78 | 2.00 | 356.89 | 20.50 | 11.39
NNGP response | $\phi$ | 11.30 | 0.25 | 1.00 | 23.00 | 0.13 | 1.00
$\sigma^{2}$ | 12.17 | 0.26 | 1.00 | 21.45 | 0.12 | 1.00
$\tau^{2}$ | 201.44 | 4.38 | 1.00 | 310.28 | 1.80 | 1.00
NNGP latent | $\phi$ | 5.81 | 0.09 | 0.38 | 3.88 | 0.02 | 0.16
$\sigma^{2}$ | 8.39 | 0.13 | 0.51 | 7.56 | 0.04 | 0.33
$\tau^{2}$ | 207.61 | 3.33 | 0.76 | 71.34 | 0.37 | 0.21

Table 1: MCMC efficiency of GriPS in the estimation of Matérn parameters in
univariate settings. For posterior samples of size 2500, we report effective
sample size—absolute, per unit time, and relative to the NNGP response model.
Figure 3: Effective sample size (log-scaled) in estimating $\sigma^{2}$,
$\phi$, $\tau^{2}$ at different Matérn smoothness ($\nu=0.5$ or $\nu=1.5$),
from an MCMC sample of size 2500.

|  | $\nu=0.5$ | $\nu=1.5$
---|---|---|---
Estimation | RMSE | M%E | UQ: Covg | RMSE | M%E | UQ: Covg
QMGP GriPS | $\phi$ | 1.1705 | 0.7463 | 0.8100 | 1.9102 | 0.1176 | 0.6850
$\sigma^{2}$ | 1.8981 | 0.3764 | 0.8350 | 1.8311 | 0.3332 | 0.8600
$\tau^{2}$ | 0.0048 | 0.0628 | 0.9450 | 0.0015 | 0.0239 | 0.9600
QMGP latent | $\phi$ | 1.7042 | 0.8251 | 0.6450 | 1.9753 | 0.1581 | 0.8600
$\sigma^{2}$ | 4.3575 | 0.9618 | 0.6400 | 4.9800 | 0.6827 | 0.8500
$\tau^{2}$ | 0.0069 | 0.0955 | 0.7150 | 0.0018 | 0.0282 | 0.9200
NNGP response | $\phi$ | 1.4851 | 1.0549 | 0.5500 | 1.5224 | 0.1494 | 0.7450
$\sigma^{2}$ | 1.9187 | 0.4221 | 0.5450 | 1.4961 | 0.2743 | 0.7150
$\tau^{2}$ | 0.0050 | 0.0649 | 0.8650 | 0.0015 | 0.0227 | 0.9500
NNGP latent | $\phi$ | 1.2829 | 0.9697 | 0.5150 | 7.6886 | 1.0344 | 0.3700
$\sigma^{2}$ | 2.0355 | 0.4059 | 0.5100 | 2.6931 | 0.5092 | 0.4400
$\tau^{2}$ | 0.0059 | 0.0754 | 0.7950 | 2.9921 | 27.3383 | 0.3350
SPDE INLA | $\phi$ | 6.1887 | 6.0161 | 0.0000 |  |  |
$\sigma^{2}$ | 2.5137 | 0.6386 | 0.0050 | unsupported in R-INLA
$\tau^{2}$ | 0.0520 | 0.7080 | 0.0050 |  |  |

Table 2: Accuracy in estimating the Matérn parameters
$\phi,\sigma^{2},\tau^{2}$ in univariate settings: root mean square error
(RMSE), mean absolute percentage error (M%E), and uncertainty quantification
in terms of coverage (UQ: Covg, with target 95%), averaged across 400
datasets.

| $\nu=0.5$ | $\nu=1.5$
---|---|---
Prediction | RMSPE | CRPS | UQ: Covg | RMSPE | CRPS | UQ: Covg
QMGP-GriPS | 0.3489 | 0.1960 | 0.9445 | 0.2651 | 0.1490 | 0.9478
QMGP latent | 0.3512 | 0.1967 | 0.9509 | 0.2677 | 0.1503 | 0.9438
NNGP response | 0.3474 | 0.1955 | 0.9464 | 0.2709 | 0.1523 | 0.9350
NNGP latent | 0.3462 | 0.1945 | 0.9476 | 0.8003 | 0.5244 | 0.9723
SPDE-INLA | 0.3462 | 0.1948 | 0.9602 | unsupported in R-INLA

Table 3: Accuracy in predictions at 10000 unobserved locations, averaged
across 400 datasets, measured as root mean square prediction error (RMSPE),
continuous ranked probability score (CRPS), and uncertainty quantification
measured as coverage of 95% prediction intervals.

QMGPs based on GriPS are compared with unmodified “latent” QMGPs and NNGPs,
both of which also sample the latent process, but without reparametrization or
splitting. We also consider a “response” NNGP model which avoids sampling the
latent process, and a stochastic partial differential equations approach
(Lindgren et al., 2011) via integrated nested Laplace approximations (INLA;
Rue et al., 2009) that avoids MCMC altogether. The SPDE-INLA approach scales
to big data via a sparse GMRF approximation of the Gaussian field which
requires the user to define a mesh of points. We fix the mesh at the test
set—also the reference set in GriPS. Carefully assessing MCMC convergence for
each of the 400 datasets is impractical. Hence, we fix the number of
iterations to 5000 for all MCMC-based methods. In these settings, GriPS is
almost always faster in computing 5000 MCMC iterations than the non sampling-
based INLA SPDE method on the same task. We provide more details on the
implemented methods in the supplementary material, along with additional
figures and summaries of computing time for each approach.

Sampling efficiency in terms of effective sample size per unit time is shown
in Table 1. As expected, the “latent” NNGP method is less efficient than its
“response” counterpart. The latter displays performance similar to the
unmodified QMGP model. On the other hand, the GriPS method is massively more
efficient than all other methods when estimating $\sigma^{2}$ and $\phi$—the
relative efficiency improvements are between one and two orders of magnitude.

These efficiency improvements are reflected in the accuracy of estimations, as
shown in Table 2. In particular, GriPS provide more reliable uncertainty
quantification as evidenced by the better coverage of the posterior 95%
credible intervals (uncertainty intervals for all 400 data sets are reported
in the supplementary material). Furthermore, GriPS-based methods most
accurately estimate the exponential covariance parameters, with the smoother
$\nu=1.5$ slightly favoring an NNGP response model. However, in the latter
scenario GriPS still achieves the best predictive performance, as shown in
Table 3.

We have, thus far, considered unmodified QMGPs in our comparison, which use
${\cal S}={\cal T}$ with irregularly spaced locations. As a consequence, the
net effect of the PS parametrization alone is unclear. We address this matter
in the supplementary material by considering gridded data, in which case all
QMGP models use the same reference set and thus all performance differences
will be attributable to the PS parametrization. The results of our analysis
show that the efficiency improvements due to PS are even larger than shown in
this section, and QMGP-GriPS is shown to be more than two orders of magnitude
more efficient than an NNGP model of the response in estimating $\sigma^{2}$
and $\phi$.

##  5 Alaskan LiDAR data

Figure 4: At the top: tree cover in 2000 (left) and forest loss year (right)
over the entire study area. At the bottom: the two outcomes (p.90 and f.cvr)
shown at a smaller subregion. The strips along which the outcomes are recorded
appear as thin lines in the top maps.

We consider a large set of remotely sensed light detection and ranging (LiDAR)
data, collected using the National Aeronautics and Space Administration (NASA)
Goddard’s LiDAR, Hyperspectral, and Thermal (G-LiHT) Airborne Imager (Cook et
al., 2013) sensor in Summer 2014 over a large region in Alaska. The data set
includes observations of forest canopy height (labeled p.90 and measured in
meters) and fractional forest cover (f.cvr) at 12,270,334 spatial locations,
which we seek to jointly characterize in a multivariate regression model.
Owing to the low flying altitude (355m) and narrow field of view (30°) of the
sensor, the observed data locations are not uniformly spread out in the study
area. Instead, they are highly concentrated along swaths less than 200m wide
and separated by 9km of unobserved areas. This data set has been considered in
the context of univariate spatial regressions in Finley et al. (2019), to
which we refer for additional details.

We select a sub-region of the data of size $150\text{km}\times 250\text{km}$
covering about 31.2% of the observed data, and sample 2.5 million spatial
locations uniformly at random in this area to build the data set we consider.
The resulting study area includes 17 strips of observed data. The effective
dimension of our data is thus $nq=5\times 10^{6}$ accounting for $q=2$. We
seek to make predictions of both outcomes on the full sub-region (including
the wide areas separating the observed strips) based on the spatial regression
model (14) and using two covariates that are available across the entire study
area as shown in Figure 4. The two covariates are: (i) percent forest cover
measured in the year 2000 (labeled tc); and (ii) an indicator variable of
forest loss (loss) that gives the number of years since 2000 that a location
had a forest-destroying disturbance such as a fire. Outcomes and predictors
are shown in Figure 4.

We build the reference set ${\cal S}$ for GriPS by taking a grid (colored in
orange in Figure 5) of regularly-spaced locations. Its size is
$n_{\text{orange}}=1.52\times 10^{6}$ which we use to cover the selected study
area. However, this grid is not dense enough to accurately characterize
spatial dependence since only a small portion of the reference locations are
near the observed outcomes. For this reason, we augment ${\cal S}$ with an
additional patterned grid (colored in blue in Figure 5)—of size
$n_{\text{blue}}=5.44\times 10^{5}$—that is built to cover the 17 strips, but
not the unobserved areas. Rather than aiming for perfectly overlapping this
grid with the data, we create a regular pattern of blue strips by equally
spacing their centers. The slight off-centering of some data swaths along this
pattern motivates us to expand each grid strip symmetrically in both
directions.

Figure 5: Construction of ${\cal S}$ by overlaying a grid of regularly-spaced
locations (orange) to a patterned grid (blue); alignment with the data swaths
is achieved by rotating ${\cal S}$.

We implement the QMGP model on ${\cal S}$ by partitioning the coordinates into
$33$ and $500$ intervals, respectively. This patterned setup ensures that only
a few matrix inverses have to be computed at each iteration. Given that
$\text{p.90}>0$ and $0<\text{f.cvr}<1$, we fit (14) on
$\widetilde{\text{p.90}}=\log\\{\text{p.90}\\}$ and analogously take
$\widetilde{\text{f.cvr}}=\log\\{\text{f.cvr}/(1-\text{f.cvr})\\})$. We ran
MCMC for a total of 20,000 iterations, of which we discarded the first 5,000
as burn-in. Predictions on a test set made of 10,000 left-out observations are
based on a thinned chain consisting of a total of 500 samples—memory
constraints inhibit storage of a significantly larger number of samples—and
shown in Table 5. Full prediction maps over the whole study areas along with
uncertainty bands are shown in Figure 6. All predictions and their summaries
have been performed after inverse-transforming all posterior samples from
MCMC.

We implement GriPS on two alternative partitions of the spatial domain. The
first partitioning scheme tessellates the domain into 33 intervals along the
horizontal and 500 along the vertical axis, resulting in a coarse partitioning
of ${\cal S}$. The second uses the same partitioning for the horizontal axis,
but $1000$ intervals for the vertical axis. Prediction results from the two
implementations are reported in Table 5. All other figures and tables refer to
the coarse partition ($33\times 500$) model. We observe that the coarser
partitioning strategy is conservative as it results in relatively large
regions including up to 140 spatial locations. Finer partitioning leads to
faster computations at a minor cost in terms of predictive accuracy. Finally,
Table 4 provides marginal summaries of the posterior distribution of the
unknown covariance and regression parameters. Marginal posterior densities and
the full Markov chains are available in the Appendix. We find that loss has a
negative impact on both outcomes even after accounting for spatial variability
via the bivariate latent process. The two margins are positively associated as
evidenced by the positive value of the off-diagonal term of
$\boldsymbol{\Lambda}$ (their correlation at a distance $\boldsymbol{h}=0$ is
estimated to be about $\sim 0.686$).

Figure 6: Predictions of p.90 and f.cvr and uncertainty quantification in
terms of the width of 95% credible intervals.

| $\boldsymbol{\theta}$ | $\lambda_{\text{p.90}}$ | $\phi_{\text{p.90}}$ | $\lambda_{\text{f.cvr}}$ | $\phi_{\text{f.cvr}}$ | $\lambda_{\text{p.90},\text{f.cvr}}$
---|---|---|---|---|---
2.5% quantile | 0.5085 | 232.5055 | 1.1588 | 357.3680 | 1.0777
mean | 0.5213 | 246.4702 | 1.1798 | 376.1024 | 1.1118
97.5% quantile | 0.5354 | 259.7343 | 1.2012 | 395.0176 | 1.1488
| $\bbeta$ and $\boldsymbol{\tau}^{2}$ | $\beta_{\text{tc}}$ | $\beta_{\text{tc}}$ | $\beta_{\text{loss}}$ | $\beta_{\text{loss}}$ | $\tau^{2}_{\text{p.90}}$ | $\tau^{2}_{\text{f.cvr}}$
---|---|---|---|---|---|---
(on p.90) | (on f.cvr) | (on p.90) | (on f.cvr)
2.5% quantile | 0.0032 | 0.0162 | -0.0018 | -0.0192 | 0.0157 | 0.1325
mean | 0.0032 | 0.0163 | -0.0017 | -0.0187 | 0.0163 | 0.1430
97.5% quantile | 0.0033 | 0.0164 | -0.0015 | -0.0182 | 0.0171 | 0.1533

Table 4: Estimated posterior summaries for unknown covariance parameters $\boldsymbol{\theta}$, regression coefficients $\bbeta$, and measurement errors $\boldsymbol{\tau}^{2}$ for the GriPS regression on the Alaskan LiDAR data set. |  $300\times 500$ partition |  $33\times 1000$ partition
---|---|---
| p.90 | f.cvr | p.90 | f.cvr
RMSPE | 1.7045 | 0.1181 | 1.7320 | 0.1197
Covg (95%) | 0.9396 | 0.9384 | 0.9400 | 0.9401
Time/iteration | 8.18s | 3.98s
Num. iterations | 20,000 | 30,000
Total time | 45.46h | 33.17h
Table 5: Prediction results over the test set, and timing, of the two
implemented models.

##  6 Discussion

We have introduced a Grid-Parametrize-Split (GriPS) strategy for resolving
inefficiencies in MCMC sampling of spatial regression models and demonstrated
its validity as an alternative to methods not based on sampling the latent
process a posteriori. Our development of the methodology and the subsequent
numerical experiments cogently illustrate the effectiveness of GriPS in
conducting spatial big data analysis. Our focus on MGPs, Matérn covariances
and LMCs, and Gaussian outcomes suggests several possible directions for
future research. On one hand, we expect our insight to be applicable more
generally to Bayesian estimation of spatial regressions with latent processes
(not necessarily MGPs). Moreover, research may focus on efficient Bayesian
parametrizations of covariance models alternative to the Matérn and/or the LMC
(see, e.g., Gneiting et al., 2010; Genton and Kleiber, 2015; Zhang and
Banerjee, 2020), and on applying similar ideas to non Gaussian outcomes.

Other natural extensions include adapting DAG-based processes to spatial-
temporal data analysis. Challenges specific to spatial-temporal domains
include the identifying “neighbors” over space and time as the scale of
associations along space significantly differ from that along time. Dynamic
NNGPs (Datta et al., 2016b) resolve this problem by estimating neighbors
instead of fixing them. QMGPs have been implemented on spatiotemporal domains
(Peruzzi et al., 2020), but efficient estimation of the covariance parameters
via MCMC remains challenging. Therefore, extensions of GriPS to spatial-
temporal domains could circumvent some of these problems and can, possibly,
provide easier and more effective analysis of spatial-temporal big data.

## Funding

M.P. and D.B.D. have received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 856506), and grant R01ES028804 of the United States
National Institutes of Health (NIH). A.O.F. received funding from the U.S.
National Science Foundation (NSF) grants DMS-1916395 and EF-1253225, and NASA
Carbon Monitoring System program. S.B. received funding from the National
Science Foundation (NSF) through grants NSF/DMS 1916349 and NSF/IIS 1562303,
and from the National Institute of Environmental Health Sciences (NIEHS)
through grants R01ES030210 and R01ES027027.

##  Appendix A Full conditional updates of $\boldsymbol{{\cal A}}$ and
$\boldsymbol{\tau}^{2}$

If the data are gridded or we sample $\boldsymbol{r}$ at non-reference
locations, full conditionals for $\boldsymbol{{\cal A}}$ and
$\boldsymbol{\tau}$ are available. Let $\boldsymbol{Y}$ be the $n\times q$
matrix whose $j$th column is
$\boldsymbol{Y}_{j}=\boldsymbol{y}_{j}-\boldsymbol{X}_{j}\bbeta_{j}$, i.e. the
$j$th outcome at locations ${\cal T}$, minus the static linear predictor for
the same outcome. Similarly denote as $\boldsymbol{F}$ the $n\times k$ matrix
whose $j$th column is the $j$th margin of the $\boldsymbol{r}$ process. We can
then write $\boldsymbol{Y}=\boldsymbol{F}\boldsymbol{{\cal
A}}^{\top}+\boldsymbol{E}$ where $\boldsymbol{E}$ is a $n\times q$ matrix such
that $\text{vec}\\{\boldsymbol{E}\\}\sim
N(\mathbf{0},I_{n}\otimes\boldsymbol{D})$, with $\text{vec}\\{\cdot\\}$ the
vectorization operator. For outcome $j$, let $\boldsymbol{F}_{(j)}$ subset
$\boldsymbol{F}$ to columns corresponding to the nonzero elements in the $j$th
row of $\boldsymbol{{\cal A}}$, denoted as $\boldsymbol{{\cal A}}_{j\setminus
0}$. Then, with
$\boldsymbol{D}=\text{diag}\\{\tau_{1}^{2},\dots,\tau_{q}^{2}\\}$ the full
conditional distribution of $\tau_{j}^{2}$ is $InvGamma(a_{j}^{*},b^{*}_{j})$
whose parameters are $a_{j}^{*}=a+n_{j}/2$ and
$b^{*}_{j}=b+\frac{1}{2}\hat{\boldsymbol{E}}_{j}^{\top}\hat{\boldsymbol{E}}_{j}$,
where
$\hat{\boldsymbol{E}}_{j}=\boldsymbol{Y}_{j,\text{obs}}-\boldsymbol{F}_{(j),\text{obs}}\boldsymbol{{\cal
A}}_{j\setminus 0}^{\top}$ and $\boldsymbol{Y}_{j,\text{obs}}$ is the
$n_{j}\times 1$ vector collecting the observed elements of
$\boldsymbol{Y}_{j}$. Finally, $\boldsymbol{{\cal A}}_{j\setminus 0}$ has
Gaussian full conditional distribution with precision
$\boldsymbol{V}_{\boldsymbol{{\cal
A}}_{j}}^{{}^{(n)}-1}=\boldsymbol{V}_{\boldsymbol{{\cal
A}}_{j}}^{-1}+\boldsymbol{F}_{(j),\text{obs}}^{\top}(I_{n_{j}}\otimes\boldsymbol{D}^{-1})\boldsymbol{F}_{(j),\text{obs}}$
and mean $\boldsymbol{m}_{\boldsymbol{{\cal
A}}_{j}}^{(n)}=\boldsymbol{V}_{\boldsymbol{{\cal
A}}_{j}}^{{}^{(n)}}\boldsymbol{F}_{(j),\text{obs}}^{\top}(I_{n_{j}}\otimes\boldsymbol{D}^{-1})\boldsymbol{Y}_{j,\text{obs}}$.

##  Appendix B Comparisons on synthetic gridded data

We compare QMGP-GriPS to unmodified QMGPs and NNGP models of the response on
gridded data. Since the data grid can be directly used as reference set, we
can write
$\boldsymbol{y}(\boldsymbol{\ell})=\boldsymbol{X}^{\top}\bbeta+\boldsymbol{w}(\boldsymbol{\ell})+\boldsymbol{\varepsilon}(\boldsymbol{\ell})$
where $\boldsymbol{\varepsilon}(\boldsymbol{\ell})\sim N(0,\boldsymbol{D})$
and we will therefore not need to modify QMGPs. We thus relabel GriPS as QMGP-
PS in this section.

Figure 7: Effective sample size in the gridded synthetic data application.

While the overall simulation design is the same as in the main article, here
we sample the univariate Gaussian process with Matérn covariance at $n=15625$
spatial locations on a $125\times 125$ grid using ${\cal D}=[0,1]^{2}$ as our
spatial domain. Then, all process realizations inside the $[0.45,0.55]^{2}$
region and at about $10\%$ of the remaining locations (chosen uniformly at
random) are removed from the data and considered as test set for predictions.
All other design choices remain the same, implying that once again our
comparisons are based on 400 datasets.

We report MCMC efficiency summaries in Table 6 and estimation and prediction
results in Table 7. Compared to the synthetic data analysis in the main
article and relative to unmodified QMGPs, there is only a slight improvement
in QMGP-PS relative to QMGP-GriPS. This is a positive result for GriPS since
it implies the loss in performance when choosing ${\cal S}\neq{\cal T}$ is
only minor, and attributes most of the improvements to the PS parametrization.
Furthermore, the simplifications arising in QMGPs from gridded data result in
much faster computations. Relative to NNGP models of the response, we show
here an additional 2-fold (or larger) improvement in MCMC efficiency, compared
to the nongridded data setting. Estimation and prediction accuracy shown in
Table 7 largely overlap with the nongridded data setting.

|  | $\nu=0.5$ | $\nu=1.5$
---|---|---|---
Efficiency | ESS | ESS/s | Relative | ESS | ESS/s | Relative
QMGP -PS | $\phi$ | 270.13 | 24.89 | 148.61 | 142.59 | 13.15 | 222.56
$\sigma^{2}$ | 242.91 | 22.39 | 117.32 | 219.82 | 20.27 | 341.84
$\tau^{2}$ | 101.94 | 9.36 | 3.04 | 618.25 | 56.96 | 53.44
QMGP latent | $\phi$ | 3.96 | 0.28 | 1.69 | 3.87 | 0.28 | 4.70
$\sigma^{2}$ | 4.97 | 0.36 | 1.87 | 5.58 | 0.40 | 6.75
$\tau^{2}$ | 149.04 | 10.67 | 3.46 | 215.96 | 15.48 | 14.52
NNGP response | $\phi$ | 7.84 | 0.17 | 1.00 | 10.45 | 0.06 | 1.00
$\sigma^{2}$ | 8.93 | 0.19 | 1.00 | 10.48 | 0.06 | 1.00
$\tau^{2}$ | 144.05 | 3.08 | 1.00 | 187.57 | 1.07 | 1.00

Table 6: MCMC efficiency in estimating the Matérn parameters in univariate
gridded data settings. For posterior samples of size 2500, we report effective
sample size—absolute, per unit time, and relative to the NNGP response model.

|  |  | $\nu=0.5$ | $\nu=1.5$
---|---|---|---
Estimation | RMSE | M%E | UQ: Covg | RMSE | M%E | UQ: Covg
QMGP -PS | $\phi$ | 1.0655 | 0.7269 | 0.8650 | 1.4607 | 0.0925 | 0.7200
$\sigma^{2}$ | 1.7736 | 0.3866 | 0.8600 | 1.9217 | 0.3385 | 0.8100
$\tau^{2}$ | 0.0047 | 0.0584 | 0.7600 | 0.0017 | 0.017 | 0.8550
QMGP latent | $\phi$ | 1.6210 | 1.0024 | 0.5750 | 4.0534 | 0.3468 | 0.4700
$\sigma^{2}$ | 2.9009 | 0.6010 | 0.5700 | 2.8427 | 0.7441 | 0.4550
$\tau^{2}$ | 0.0056 | 0.0696 | 0.6850 | 1.3032 | 8.0071 | 0.7400
NNGP response | $\phi$ | 1.3121 | 1.0352 | 0.5550 | 1.7365 | 0.1765 | 0.7300
$\sigma^{2}$ | 1.8878 | 0.4124 | 0.5350 | 1.5838 | 0.3094 | 0.7000
$\tau^{2}$ | 0.0034 | 0.0424 | 0.8500 | 0.0009 | 0.0137 | 0.8800
|  | $\nu=0.5$ | $\nu=1.5$
---|---|---
Prediction | RMSPE | CRPS | UQ: Covg | RMSPE | CRPS | UQ: Covg
QMGP-PS | 0.3373 | 0.1886 | 0.9502 | 0.2670 | 0.1461 | 0.9490
QMGP latent | 0.3374 | 0.1887 | 0.9506 | 0.4816 | 0.2607 | 0.9554
NNGP response | 0.3432 | 0.1914 | 0.9778 | 0.2789 | 0.1512 | 0.9782

Table 7: Top: accuracy in estimating the Matérn parameters $\phi,\sigma^{2}$
plus $\tau^{2}$, in univariate settings—root mean square error (RMSE), mean
absolute percentage error (M%E), and uncertainty quantification in terms of
coverage (UQ: Covg, with target 95%), averaged across 400 datasets. Bottom:
accuracy in predictions—root mean square prediction error (RMSPE), continuous
ranked probability score (CRPS), and 95% coverage.

##  Appendix C Details on implemented models

In the synthetic data applications, the prior distributions on $\tau^{2}$ was
chosen as an inverse Gamma with parameters $2.01$ and $1$ for QMGP-GriPS
latent QMGPs, and both NNGP models. $\pi(\sigma^{2})=\text{Inv.Gamma}(2.01,1)$
for latent QMGPs and NNGPs, whereas in GriPS this parameter is expanded into
$\sigma^{2}$ with prior $\pi(\sigma^{2})=\text{Inv.Gamma}(\epsilon,\epsilon)$
with $\epsilon$ small (we set it to $10^{-3}$), and
$\lambda\sim\pi(\lambda)=N_{\lambda>0}(0,1)$. The SPDE-INLA approach was
implemented for $\nu=0.5$ following Krainski et al. (2019), Chapter 2, setting
$\alpha=\nu+1=1.5$. Setting $\nu=1.5$ was unsupported by inla.spde2.pcmatern
in R-INLA. The priors on $\sigma^{2},\phi$ were such that
$Pr(\sigma>2.5)=0.01$ and $Pr(\phi<0.1)=0.01$. The mesh for INLA was fixed at
a $10000$-sized grid which overlapped with the spatial locations of the test
set. In other words we used for INLA the same mesh as the reference set in the
QMGP-GriPS method. Smaller meshes are associated with faster computations in
INLAs but negatively impact predictions and estimation accuracy (especially
about the nugget $\tau^{2}$). All models were run in a Ubuntu 20.04
workstation based on a 12-core, 24-thread AMD Ryzen 9 3900 CPU and 128GB of
memory; R version 4.0.3 was linked to Intel MKL 2019.5-075 which includes
efficient LAPACK (Anderson et al., 1999) and BLAS (Blackford et al., 2002)
libraries for linear algebra. GriPS-based QMGPs were run on 12 threads using
OpenMP (Dagum and Menon, 1998), all other models on 22 threads. Figure 8
summarises the compute times for each model, whereas Figure 9 reports all 95%
posterior intervals for the estimation of the unknown covariance parameters.
Effective sample size reported in the tables throughout the article was
calculated using the effectiveSize function of R package coda (Plummer et al.,
2006).

Figure 8: Compute time for each model. The top figure refers to the analysis
and models implemented in the main article, the bottom to those in this
supplement.

Figure 9: Coverage of posterior 95% intervals in all 400 data sets in the
synthetic data application appearing in the main article, and for all models.
Orange lines represent the true values; blue lines are 95% (marginal)
intervals that include the true value, grey lines are intervals that do not
include the true value. Intervals are calculated using MCMC samples of size
2500. SPDE-INLA intervals are extracted from the approximated marginal
posterior distributions.

##  Appendix D Additional details on the LiDAR data application

Figures 10 and 11 add to the posterior summaries reported in the main article.

Figure 10: Markov chains for the unknown parameters in the regression model
for the LiDAR data. Figure 11: Estimated marginal posterior densities of the
covariance parameters
$\boldsymbol{\theta}=\\{\boldsymbol{\Lambda},\boldsymbol{\phi}\\}$ and
regression parameters $\bbeta$ and $\boldsymbol{\tau}^{2}$.

## References

  * Ambikasaran et al. (2016) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., and O’Neil, M. (2016), “Fast Direct Methods for Gaussian Processes,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252–265, doi:10.1109/TPAMI.2015.2448083.
  * Anderson et al. (1999) Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., and Sorensen, D. (1999), LAPACK Users’ Guide, Philadelphia, PA: Society for Industrial and Applied Mathematics, 3rd ed.
  * Banerjee (2017) Banerjee, S. (2017), “High-Dimensional Bayesian Geostatistics,” Bayesian Analysis, 12, 583–614, doi:10.1214/17-BA1056R.
  * Banerjee (2020) — (2020), “Modeling Massive Spatial Datasets Using a Conjugate Bayesian Linear Modeling Framework,” Spatial Statistics, in press, doi:10.1016/j.spasta.2020.100417.
  * Banerjee et al. (2014) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2014), Hierarchical Modeling and Analysis for Spatial Data, 2nd edition, Chapman and Hall/CRC, doi:10.1201/b17115.
  * Banerjee et al. (2008) Banerjee, S., Gelfand, A. E., Finley, A. O., and Sang, H. (2008), “Gaussian predictive process models for large spatial data sets,” Journal of the Royal Statistical Society, Series B, 70, 825–848, doi:10.1111/j.1467-9868.2008.00663.x.
  * Blackford et al. (2002) Blackford, L. S., Petitet, A., Pozo, R., Remington, K., Whaley, R. C., Demmel, J., Dongarra, J., Duff, I., Hammarling, S., Henry, G., et al. (2002), “An updated set of basic linear algebra subprograms (BLAS),” ACM Transactions on Mathematical Software, 28, 135–151.
  * Cook et al. (2013) Cook, B. D., Corp, L. A., Nelson, R. F., Middleton, E. M., Morton, D. C., McCorkel, J. T., Masek, J. G., Ranson, K. J., Ly, V., and Montesano, P. M. (2013), “NASA Goddard’s LiDAR, Hyperspectral and Thermal (G-LiHT) Airborne Imager,” Remote Sensing, 5, 4045–4066, doi:10.3390/rs5084045.
  * Cressie and Johannesson (2008) Cressie, N. and Johannesson, G. (2008), “Fixed Rank Kriging for Very Large Spatial Data Sets,” Journal of the Royal Statistical Society, Series B, 70, 209–226, doi:10.1111/j.1467-9868.2007.00633.x.
  * Dagum and Menon (1998) Dagum, L. and Menon, R. (1998), “OpenMP: an industry standard API for shared-memory programming,” Computational Science & Engineering, IEEE, 5, 46–55.
  * Datta et al. (2016a) Datta, A., Banerjee, S., Finley, A. O., and Gelfand, A. E. (2016a), “Hierarchical Nearest-Neighbor Gaussian Process Models for Large Geostatistical Datasets,” Journal of the American Statistical Association, 111, 800–812, doi:10.1080/01621459.2015.1044091.
  * Datta et al. (2016b) Datta, A., Banerjee, S., Finley, A. O., Hamm, N. A. S., and Schaap, M. (2016b), “Nonseparable dynamic nearest neighbor Gaussian process models for large spatio-temporal data with an application to particulate matter analysis,” The Annals of Applied Statistics, 10, 1286–1316, doi:10.1214/16-AOAS931.
  * Dunson and Johndrow (2020) Dunson, D. and Johndrow, J. E. (2020), “The Hastings algorithm at fifty,” Biometrika, 107, 1–23, doi:10.1093/biomet/asz066.
  * Eidsvik et al. (2014) Eidsvik, J., Shaby, B. A., Reich, B. J., Wheeler, M., and Niemi, J. (2014), “Estimation and Prediction in Spatial Models With Block Composite Likelihoods,” Journal of Computational and Graphical Statistics, 23, 295–315, doi:10.1080/10618600.2012.760460.
  * Finley et al. (2008) Finley, A. O., Banerjee, S., Ek, A. R., and McRoberts, R. E. (2008), “Bayesian multivariate process modeling for prediction of forest attributes,” Journal of Agricultural, Biological, and Environmental Statistics, 13, 60, doi:10.1198/108571108X273160.
  * Finley et al. (2019) Finley, A. O., Datta, A., Cook, B. D., Morton, D. C., Andersen, H. E., and Banerjee, S. (2019), “Efficient Algorithms for Bayesian Nearest Neighbor Gaussian Processes,” Journal of Computational and Graphical Statistics, 28, 401–414, doi:10.1080/10618600.2018.1537924.
  * Finley et al. (2009) Finley, A. O., Sang, H., Banerjee, S., and Gelfand, A. E. (2009), “Improving the performance of predictive process modeling for large datasets,” Computational Statistics and Data Analysis, 53, 2873–2884, doi:10.1016/j.csda.2008.09.008.
  * Furrer et al. (2006) Furrer, R., Genton, M. G., and Nychka, D. (2006), “Covariance Tapering for Interpolation of Large Spatial Datasets,” Journal of Computational and Graphical Statistics, 15, 502–523, doi:10.1198/106186006X132178.
  * GEDI (2020) GEDI (2020), “Global Ecosystem Dynamics Investigation LiDAR,” https://gedi.umd.edu.
  * Gelfand et al. (1995) Gelfand, A. E., Sahu, S. K., and Carlin, B. P. (1995), “Efficient Parametrisations for Normal Linear Mixed Models,” Biometrika, 82, 479–488, doi:10.2307/2337527.
  * Gelman (2006) Gelman, A. (2006), “Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper),” Bayesian Analysis, 1, 513–533, doi:10.1214/06-BA117A.
  * Genton and Kleiber (2015) Genton, M. G. and Kleiber, W. (2015), “Cross-Covariance Functions for Multivariate Geostatistics,” Statistical Science, 30, 147–163, doi:10.1214/14-STS487.
  * Geoga et al. (2020) Geoga, C. J., Anitescu, M., and Stein, M. L. (2020), “Scalable Gaussian Process Computations Using Hierarchical Matrices,” Journal of Computational and Graphical Statistics, 29, 227–237, doi:10.1080/10618600.2019.1652616.
  * Ghosh and Dunson (2009) Ghosh, J. and Dunson, D. B. (2009), “Default Prior Distributions and Efficient Posterior Computation in Bayesian Factor Analysis,” Journal of Computational and Graphical Statistics, 18, 306–320, doi:10.1198/jcgs.2009.07145.
  * Gneiting et al. (2010) Gneiting, T., Kleiber, W., and Schlather, M. (2010), “Matérn Cross-Covariance Functions for Multivariate Random Fields,” Journal of the American Statistical Association, 105, 1167–1177, doi:10.1198/jasa.2010.tm09420.
  * Gonzalez et al. (2011) Gonzalez, J., Low, Y., Gretton, A., and Guestrin, C. (2011), “Parallel Gibbs Sampling: From Colored Fields to Thin Junction Trees,” in Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, eds. Gordon, G., Dunson, D., and Dudík, M., Fort Lauderdale, FL, USA: PMLR, vol. 15 of Proceedings of Machine Learning Research, pp. 324–332.
  * Gramacy and Apley (2015) Gramacy, R. B. and Apley, D. W. (2015), “Local Gaussian Process Approximation for Large Computer Experiments,” Journal of Computational and Graphical Statistics, 24, 561–578, doi:10.1080/10618600.2014.914442.
  * Guhaniyogi and Banerjee (2018) Guhaniyogi, R. and Banerjee, S. (2018), “Meta-Kriging: Scalable Bayesian Modeling and Inference for Massive Spatial Datasets,” Technometrics, 60, 430–444, doi:10.1080/00401706.2018.1437474.
  * Guhaniyogi et al. (2011) Guhaniyogi, R., Finley, A. O., Banerjee, S., and Gelfand, A. E. (2011), “Adaptive Gaussian predictive process models for large spatial datasets,” Environmetrics, 22, 997–1007, doi:10.1002/env.1131.
  * Heaton et al. (2019) Heaton, M. J., Datta, A., Finley, A. O., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, F., Gramacy, R. B., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D. W., Sun, F., and Zammit-Mangion, A. (2019), “A Case Study Competition Among Methods for Analyzing Large Spatial Data,” Journal of Agricultural, Biological and Environmental Statistics, 24, 398–425, doi:10.1007/s13253-018-00348-w.
  * Huang and Sun (2018) Huang, H. and Sun, Y. (2018), “Hierarchical Low Rank Approximation of Likelihoods for Large Spatial Datasets,” Journal of Computational and Graphical Statistics, 27, 110–118, doi:10.1080/10618600.2017.1356324.
  * ICESat-2 (2020) ICESat-2 (2020), “Ice, Cloud, and Land Elevation Satellite-2,” https://icesat-2.gsfc.nasa.gov/.
  * Jurek and Katzfuss (2020) Jurek, M. and Katzfuss, M. (2020), “Hierarchical sparse Cholesky decomposition with applications to high-dimensional spatio-temporal filtering,” arXiv:2006.16901.
  * Kampe et al. (2010) Kampe, T. U., Johnson, B. R., Kuester, M. A., and Keller, M. (2010), “NEON: the first continental-scale ecological observatory with airborne remote sensing of vegetation canopy biochemistry and structure,” Journal of Applied Remote Sensing, 4, 1–24, doi:10.1117/1.3361375.
  * Kastner (2014) Kastner, Gregor Frühwirth-Schnatter, S. (2014), “Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility models,” Computational Statistics and Data Analysis, 76, 408–423, doi:10.1016/j.csda.2013.01.002.
  * Katzfuss (2017) Katzfuss, M. (2017), “A Multi-Resolution Approximation for Massive Spatial Datasets,” Journal of the American Statistical Association, 112, 201–214, doi:10.1080/01621459.2015.1123632.
  * Katzfuss and Guinness (2019) Katzfuss, M. and Guinness, J. (2019), “A general framework for Vecchia approximations of Gaussian processes,” arXiv:1708.06302.
  * Kaufman et al. (2008) Kaufman, C. G., Schervish, M. J., and Nychka, D. W. (2008), “Covariance Tapering for Likelihood-Based Estimation in Large Spatial Data Sets,” Journal of the American Statistical Association, 103, 1545–1555, doi:10.1198/016214508000000959.
  * Krainski et al. (2019) Krainski, E. T., Gómez-Rubio, V., Bakka, H., Lenzi, A., Castro-Camilo, D., Simpson, D., Lindgren, F., and Rue, H. (2019), Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA, CRC Press/Taylor and Francis Group.
  * Krause et al. (2008) Krause, A., Singh, A., and Guestrin, C. (2008), “Near-Optimal Sensor Placements in Gaussian Processes: Theory, Efficient Algorithms and Empirical Studies,” Journal of Machine Learning Research, 8, 235–284, http://www.jmlr.org/papers/v9/krause08a.html.
  * Lindgren et al. (2011) Lindgren, F., Rue, H., and Lindström, J. (2011), “An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach,” Journal of the Royal Statistical Society: Series B, 73, 423–498, doi:10.1111/j.1467-9868.2011.00777.x.
  * Liu and Chen (1995) Liu, J. S. and Chen, R. (1995), “Blind Deconvolution via Sequential Imputations,” Journal of the American Statistical Association, 90, 567–576, doi:10.1080/01621459.1995.10476549.
  * Liu and Wu (1999) Liu, J. S. and Wu, Y. N. (1999), “Parameter Expansion for Data Augmentation,” Journal of the American Statistical Association, 94, 1264–1274, doi:10.1080/01621459.1999.10473879.
  * Matheron (1982) Matheron, G. (1982), “Pour une analyse krigeante des données régionalisées,” Technical report N.732, Centre de Géostatistique.
  * Papaspiliopoulos et al. (2007) Papaspiliopoulos, O., O., R. G., and Sköld, M. (2007), “A General Framework for the Parametrization of Hierarchical Models,” Statistical Science, 22, 59–73, doi:10.1214/088342307000000014.
  * Peruzzi et al. (2020) Peruzzi, M., Banerjee, S., and Finley, A. O. (2020), “Highly Scalable Bayesian Geostatistical Modeling via Meshed Gaussian Processes on Partitioned Domains,” Journal of the American Statistical Association, in press. doi:10.1080/01621459.2020.1833889.
  * Peruzzi and Dunson (2020) Peruzzi, M. and Dunson, D. B. (2020), “Spatial multivariate trees for big data Bayesian regression,” arXiv:2012.00943.
  * Plummer et al. (2006) Plummer, M., Best, N., Cowles, K., and Vines, K. (2006), “CODA: Convergence Diagnosis and Output Analysis for MCMC,” R News, 6, 7–11.
  * Quiroz et al. (2019) Quiroz, Z. C., Prates, M. O., and Dey, D. K. (2019), “Block Nearest Neighboor Gaussian processes for large datasets,” arXiv:1604.08403.
  * Quiñonero-Candela and Rasmussen (2005) Quiñonero-Candela, J. and Rasmussen, C. E. (2005), “A Unifying View of Sparse Approximate Gaussian process Regression,” Journal of Machine Learning Research, 6, 1939–1959, https://www.jmlr.org/papers/volume6/quinonero-candela05a/quinonero-candela05a.pdf.
  * Ren and Banerjee (2013) Ren, Q. and Banerjee, S. (2013), “Hierarchical Factor Models for Large Spatially Misaligned Data: A Low-Rank Predictive Process Approach,” Biometrics, 69, 19–30, doi:10.1111/j.1541-0420.2012.01832.x.
  * Robert and Casella (2005) Robert, C. P. and Casella, G. (2005), Monte Carlo Statistical Methods (Springer Texts in Statistics), Berlin, Heidelberg: Springer-Verlag.
  * Rue and Held (2005) Rue, H. and Held, L. (2005), Gaussian Markov Random Fields: Theory and Applications, Chapman & Hall/CRC, doi:10.1007/978-3-642-20192-9.
  * Rue et al. (2009) Rue, H., Martino, S., and Chopin, N. (2009), “Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations,” Journal of the Royal Statistical Society: Series B, 71, 319–392, doi:10.1111/j.1467-9868.2008.00700.x.
  * Sang and Huang (2012) Sang, H. and Huang, J. Z. (2012), “A full scale approximation of covariance functions for large spatial data sets,” Journal of the Royal Statistical Society, Series B, 74, 111–132, doi:10.1111/j.1467-9868.2011.01007.x.
  * Shirota et al. (2019) Shirota, S., Finley, A. O., Cook, B. D., and Banerjee, S. (2019), “Conjugate Nearest Neighbor Gaussian Process Models for Efficient Statistical Interpolation of Large Spatial Data,” arXiv:1907.10109.
  * Snelson and Ghahramani (2007) Snelson, E. and Ghahramani, Z. (2007), “Local and global sparse Gaussian process approximations,” in Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics, vol. 2 of Proceedings of Machine Learning Research, pp. 524–531, http://proceedings.mlr.press/v2/snelson07a.html.
  * Stein (1990) Stein, M. (1990), “Uniform Asymptotic Optimality of Linear Predictions of a Random Field Using an Incorrect Second-Order Structure,” The Annals of Statistics, 18, 850–872, doi:10.1214/aos/1176347629.
  * Stein (1999) Stein, M. L. (1999), Interpolation of spatial data: Some theory for Kriging, Springer Series in Statistics, New York: Springer-Verlag, doi:10.1007/978-1-4612-1494-6.
  * Stein (2014) — (2014), “Limitations on low rank approximations for covariance matrices of spatial data,” Spatial Statistics, 8, 1–19, doi:doi:10.1016/j.spasta.2013.06.003.
  * Sun et al. (2011) Sun, Y., Li, B., and Genton, M. (2011), “Geostatistics for large datasets,” in Advances and Challenges in Space-time Modelling of Natural Events, eds. Montero, J., Porcu, E., and Schlather, M., Berlin Heidelberg: Springer-Verlag, pp. 55–77, doi:10.1007/978-3-642-17086-7.
  * Tang et al. (2020) Tang, W., Zhang, L., and Banerjee, S. (2020), “On identifiability and consistency of the nugget in Gaussian spatial process models,” arXiv:1908.05726.
  * Taylor-Rodriguez et al. (2019) Taylor-Rodriguez, D., Finley, A. O., Datta, A., Babcock, C., Andersen, H. E., Cook, B. D., Morton, D. C., and Banerjee, S. (2019), “Spatial factor models for high-dimensional and large spatial data: An application in forest variable mapping,” Statistica Sinica, 29, 1155–1180, doi:10.5705/ss.202018.0005.
  * van Dyk and Meng (2001) van Dyk, D. A. and Meng, X.-L. (2001), “The Art of Data Augmentation,” Journal of Computational and Graphical Statistics, 10, 1–50, doi:10.1198/10618600152418584.
  * Vecchia (1988) Vecchia, A. V. (1988), “Estimation and Model Identification for Continuous Spatial Processes,” Journal of the Royal Statistical Society, Series B, 50, 297–312, doi:10.1111/j.2517-6161.1988.tb01729.x.
  * Vihola (2012) Vihola, M. (2012), “Robust adaptive Metropolis algorithm with coerced acceptance rate,” Statistics and Computing, 22, 997–1008, doi:10.1007/s11222-011-9269-5.
  * Wackernagel (2003) Wackernagel, H. (2003), Multivariate Geostatistics: An Introduction with Applications, Springer, Berlin, doi:10.1007/978-3-662-05294-5.
  * Ying (1991) Ying, Z. (1991), “Asymptotic properties of a maximum likelihood estimator with data from a Gaussian process,” Journal of Multivariate Analysis, 36, 280–296, doi:10.1016/0047-259X(91)90062-7.
  * Yu and Meng (2011) Yu, Y. and Meng, X.-L. (2011), “To Center or Not to Center: That Is Not the Question—An Ancillarity–Sufficiency Interweaving Strategy (ASIS) for Boosting MCMC Efficiency,” Journal of Computational and Graphical Statistics, 20, 531–570, doi:10.1198/jcgs.2011.203main.
  * Zhang (2004) Zhang, H. (2004), “Inconsistent Estimation and Asymptotically Equal Interpolations in Model-Based Geostatistics,” Journal of the American Statistical Association, 99, 250–261, doi:10.1198/016214504000000241.
  * Zhang and Banerjee (2020) Zhang, L. and Banerjee, S. (2020), “Spatial Factor Modeling: A Bayesian Matrix-Normal Approach for Misaligned Data,” arXiv:2006.00595.