Highly Scalable Bayesian Geostatistical Modeling via Meshed Gaussian Processes on Partitioned Domains

Collecting large quantities of spatial and spatiotemporal data is now commonplace in many fields. In ecology and forestry, massive datasets are collected using satellite imaging and other remote sensing instruments such as LiDAR that periodically record high-resolution images. Unfortunately, clouds frequently obstruct the view resulting in large regions with missing information. Figure 1 shows this phenomenon in Normalized Difference Vegetation Index (NDVI) data from the Serengeti region. Filling such gaps in the data is an important goal as is quantifying uncertainty in predictions. This goal is achieved through stochastic modeling of the underlying phenomenon, which involves the specification of a spatial or spatiotemporal process characterizing dependence from a finite realization. Gaussian processes (GP) are a customary choice to characterize spatial dependence, but their implementation is notoriously burdened by their $O(n^{3})$ computational complexity. Consequently, intense research has been devoted in recent years to developing scalable models for large spatial datasets – see detailed reviews by Sun et al., (2011) and Banerjee, (2017).

Computational complexity can be reduced by considering low-rank models; among these, knot-based methods motivated by “kriging” ideas enjoy some optimality properties but oversmooth the estimates of spatial random effects unless the number of knots is large, and require corrections to avoid overestimation of the nugget (Banerjee et al.,, 2008; Cressie and Johannesson,, 2008; Banerjee et al.,, 2010; Guhaniyogi et al.,, 2011; Finley et al.,, 2012). Other methods reduce the computational burden by introducing sparsity in the covariance matrix; strategies include tapering (Furrer et al.,, 2006; Kaufman et al.,, 2008) or partitioning of the spatial domain into regions with a typical assumption of independence across regions (Sang and Huang,, 2012; Stein,, 2014). These can be improved by considering a recursive partitioning scheme, resulting in a multi-resolution approximation (MRA; Katzfuss, 2017). Other assumptions on conditional independence assumptions also have a good track record in terms of scalability to large spatial datasets: Gaussian random Markov random fields (GMRF; Rue and Held,, 2005), composite likelihood methods (Eidsvik et al.,, 2014), and neighbor-based likelihood approximations (Vecchia,, 1988) belong to this family.

The recent literature has witnessed substantial activity surrounding the so called Vecchia approximation (Vecchia,, 1988). This approximation can be regarded as a special case of the GMRF approximations with a simplified neighborhood structure motivated from a directed acyclic graphical (DAG) representation of a Gaussian process likelihood. Extensions leading to well-defined spatial processes to accommodate inference at arbitrary locations by extending the DAG representation to the entire domain include Nearest neighbor Gaussian processes (NNGP; Datta et al., 2016a ; Datta et al., 2016b ) and further generalizations by constructing DAGs over the augmented space of outcomes and spatial effects (Katzfuss and Guinness,, 2017). These approaches render computational scalability by introducing sparsity in the precision matrix. The DAG relies upon a specific topological ordering of the locations, which also determine the construction of neighborhood sets, and certain orderings tend to deliver improved performance of such models (Katzfuss and Guinness,, 2017; Guinness,, 2018).

When inference on the latent process is sought, Bayesian inference has the benefits of providing direct probability statements based upon the posterior distribution of the process. Inference based on asymptotic approximations are avoided, but there remain challenges in computing the posterior distribution given that inference is sought on a very high-dimensional parameter space (including the realizations of the latent process). One possibility, available for Gaussian first-stage likelihoods, is to work with a collapsed or marginalized likelihood by integrating out the spatial random effects. However, Gibbs samplers and other MCMC algorithms for the collapsed models can be inexorably slow and are impractical when data are in the millions. A sequential Gibbs sampler that updates the latent spatial effects (Datta et al., 2016a, ) is faster in updating the parameters but suffers from high autocorrelation and slow mixing. Another possibility emerges when interest lies in prediction or imputation of the outcome variable only and not the latent process. Here, a so called “response” model that models the outcome itself using an NNGP can be constructed. This model is much faster and enjoys superior convergence properties, but we lose inference on the latent process and its predictive performance tends to be inferior to the latent process model. Furthermore, these options are unavailable in non-Gaussian first-stage hierarchical models or when the focus is not uniquely on prediction. A detailed comparison of different approaches for computing Bayesian NNGP models is presented in Finley et al., (2019).

Our current contribution introduces a class of Meshed Gaussian Process (MGP) models for Bayesian hierarchical modeling of large spatial datasets. This class builds upon the aforementioned works that build upon Vecchia, (1988) and other DAG based models. The inferential focus remains within the context of massive spatial datasets over very large domains. We exploit the demonstrated benefits of the DAG based models, but we now adapt them to partitioned domains. We describe dependence across regions of a partitioned domain using a small, patterned DAG which we refer to as a mesh. Within each region, some locations are selected as reference and collectively mapped to a single node in the DAG. Relationships among nodes are governed by kriging ideas. In the resulting MGP, regions in the spatial domain depend on each other through the reference locations. Realizations at all other locations are assumed independent, conditional upon the reference locations. This construction leads to a valid standalone spatial process.

As a particular subclass of MGPs, we propose a novel partitioning and graph design based on domain tessellations. Unlike methods that build sparse DAGs by limiting dependence to $m$ nearest neighbors, our approach shapes the underlying DAG with a known, repeating pattern corresponding to the chosen tessellation geometry. The underlying sparse DAG enables scaling computations to large data settings and its known pattern guarantees the availability of block-parallel sampling schemes; furthermore, large computational savings can be achieved at no additional approximation cost if data are collected on patterned lattices. Finally, extensions to spatiotemporal and/or multivariate data are straightforward once a suitable covariance function has been defined. We use axis-parallel domain partitioning and the corresponding cubic DAG – resulting in cubic MGPs or Q-MGPs – to show substantial improvements in computational time and inferential performance relative to other models with data sizes ranging from the thousands to the several millions, for both spatial and spatiotemporal data and using multivariate spatial processes.

The present work may appear to share similarities with the block-NNGP model of Quiroz et al., (2019), who advocate building sparse DAGs on grouped locations based on their ordering and subsequent identification of $m$ “past” neighbors. Unlike block-NNGPs, our tessellated GPs consider the domain tessellation as generating the DAG; the number of parents of any node is thus fixed and depends on the geometry of the chosen tessellation rather than on a user-defined parameter. Inclusion of more distant locations in the parent set of any location will, therefore, not proceed by increasing the number of neighbors $m$, but rather by increasing the regions’ size and/or modifying their shape. Central to tessellated GPs is the idea of forcing a DAG with known coloring on the data, resulting in guaranteed efficiencies when recovering the latent spatial effects. This strategy is analogous in spirit to multi-resolution approximations (Katzfuss,, 2017; Gramacy and Lee,, 2008), which also force a DAG on the data, resulting in conditional independence patterns that are known in advance and that can be used to improve computations. However, while multi-resolution approximations are defined by branching graphs associated to recursive domain partitioning, tessellated GPs use a single domain partition, with each region connected in the DAG only to its immediate neighbors. Compared to treed graphs, tessellated GPs are associated to DAGs with fewer conditionally-independent groups and whose repeated patterns facilitate the identification of redundant matrix operations arising when one or more coordinate margins are gridded. We also note that while the idea of partitioning domains to create approximations is not new, construction of the DAG-based approximation over partitioned domains has received considerably less attention. Finally, our focus here is in developing tessellated GPs as a methodology that enables the efficient recovery of the latent spatial random effects and the Bayesian estimation of covariance parameters via MCMC; we are thus not focusing on alternative computational algorithms (see e.g. Finley et al., 2019), which have been developed for NNGPs but can nonetheless all be adapted to general MGP models.

The balance of this paper proceeds as follows. Section 2 introduces our general framework for hierarchical Bayesian modeling of spatial processes using networks of grouped spatial locations. The MGP is outlined in Section 3, where we provide a general, scalable computing algorithm in Section 3.1. Tessellation-based schemes and the specific case of Q-MGPs are outlined in Section 4, which highlights their properties and computational advantages. We illustrate the performance of our proposed approach in Section 5 using simulation experiments and an application on a massive dataset with millions of spatiotemporal locations. We conclude the paper with a discussion and pointers to further research. Supplementary material accompanying this manuscript as an Appendix is available online and contains further comparisons of Q-MGPs with several state-of-the-art methods for spatial data.

A $q\times 1$ spatial process assigns a probability law on $\{\boldsymbol{w}(\boldsymbol{\ell}):\boldsymbol{\ell}\in{\cal D}\}$, where $\boldsymbol{w}(\boldsymbol{\ell})$ is a $q\times 1$ random vector with elements $w_{i}(\boldsymbol{\ell})$ for $i=1,2,\ldots,q$. In the following general discussion we will not distinguish between spatial (${\cal D}\subset\Re^{d}$) and spatiotemporal domains (${\cal D}\subset\Re^{d+1}$), and denote spatial or spatiotemporal locations as $\boldsymbol{\ell},\boldsymbol{s}$, or $\boldsymbol{u}$.

For any finite set of spatial locations $\{\boldsymbol{\ell}_{1},\boldsymbol{\ell}_{2},\ldots,\boldsymbol{\ell}_{n_{{\cal L}}}\}={\cal L}\subset{\cal D}$ of size $n_{{\cal L}}$, let $P(\cdot)$ denote the probability law of the $n_{{\cal L}}q\times 1$ random vector $\boldsymbol{w}_{{\cal L}}=(\boldsymbol{w}(\boldsymbol{\ell}_{1})^{\top},\boldsymbol{w}(\boldsymbol{\ell}_{2})^{\top},\ldots,\boldsymbol{w}(\boldsymbol{\ell}_{n_{{\cal L}}})^{\top})^{\top}$ with probability density $p(\cdot)$. The joint density of $\boldsymbol{w}_{{\cal L}}$ can be expressed as a DAG (or a Bayesian network model) with respect to the ordered set of locations ${\cal L}$ as

$\displaystyle p(\boldsymbol{w}_{{\cal L}})$

$\displaystyle=\prod_{i=1}^{n_{{\cal L}}}p(\boldsymbol{w}(\boldsymbol{\ell}_{i})\,|\,\boldsymbol{w}(\boldsymbol{\ell}_{1}),\ldots,\boldsymbol{w}(\boldsymbol{\ell}_{i-1})),$

(1)

where the conditional set for each $\boldsymbol{w}(\boldsymbol{\ell}_{i})$ can be interpreted as the set of its parents in a large, dense Bayesian network. Defining a simplified valid joint density on ${\cal L}$ by reducing the size of the conditioning sets is a popular strategy for fast likelihood approximations in the context of large spatial datasets. One typically limits dependence to “past” neighboring locations with respect to the ordering in (1) (Vecchia,, 1988; Gramacy and Apley,, 2015; Stein et al.,, 2004; Datta et al., 2016a, ; Katzfuss and Guinness,, 2017). The neighbors are defined and fixed and model performance may benefit from the addition of some distant locations (Stein et al.,, 2004). The ordering in ${\cal L}$ is also fixed and inferential performance may benefit from the use of some fixed permutations (Guinness,, 2018). The result of shrinking the conditional sets to a smaller set of neighbors from the past yields a sparse DAG or Bayesian network, which yields potentially massive computational gains.

We proceed in a similar manner, but instead of defining a sparse DAG at the level of each individual location, we map entire groups of locations to nodes in a much smaller graph; the same graph will be used to model the dependence between any location in the spatial domain and, therefore, to define a spatial process. Let ${\cal P}=\{{\cal D}_{1},\dots,{\cal D}_{M}\}$ be a partition of ${\cal D}$ into $M$ mutually exclusive subsets so that ${\cal D}=\cup_{i=1}^{M}{\cal D}_{i}$ and ${\cal D}_{i}\cap{\cal D}_{j}=\emptyset$ whenever $i\neq j$. Similar to the nomenclature in the NNGP, we fix a reference set ${\cal S}=\{\boldsymbol{s}_{1},\dots,\boldsymbol{s}_{n_{{\cal S}}}\}\subset{\cal D}$, which itself is partitioned using ${\cal P}$ by letting ${\cal S}_{j}={\cal D}_{j}\cap{\cal S}$. The set of non-reference locations is similarly partitioned with ${\cal U}_{j}={\cal D}_{j}\setminus{\cal S}_{j}$ so that ${\cal D}_{j}={\cal S}_{j}\cup{\cal U}_{j}$ for each $j=1,2,\ldots,M$. We now construct a DAG to model dependence within and between ${\cal S}$ and ${\cal U}$. Let ${\cal G}=\{\boldsymbol{V},\boldsymbol{E}\}$ be a graph with nodes $\boldsymbol{V}=\boldsymbol{A}\cup\boldsymbol{B}$, where we refer to $\boldsymbol{A}=\{\boldsymbol{a}_{1},\dots,\boldsymbol{a}_{M}\}$ as the reference nodes and to $\boldsymbol{B}=\{\boldsymbol{b}_{1},\dots,\boldsymbol{b}_{M}\}$ as the non-reference, or simply “other”, nodes. Let $\boldsymbol{A}\cap\boldsymbol{B}=\emptyset$. We introduce a map $\eta:{\cal D}\rightarrow\boldsymbol{V}$ such that

$\displaystyle\eta(\boldsymbol{\ell})$

$\displaystyle=\left\{\begin{array}[]{l}\boldsymbol{a}_{j}\in\boldsymbol{A}\;\mbox{ if }\;\boldsymbol{\ell}\in{\cal S}_{j}\\ \boldsymbol{b}_{j}\in\boldsymbol{B}\;\mbox{ if }\;\boldsymbol{\ell}\in{\cal U}_{j}\end{array}\right.\;.$

(4)

This surjective many-to-one map links each location in ${\cal S}_{j}$ and ${\cal U}_{j}$ to a node in ${\cal G}$. The edges connecting nodes in ${\cal G}$ are $\boldsymbol{E}=\{\text{Pa}[{\boldsymbol{v}_{1}}],\dots,\text{Pa}[{\boldsymbol{v}_{2M}}]\}$ where $\text{Pa}[{\boldsymbol{v}}]\subset\boldsymbol{V}$ denotes the set of parents of any $\boldsymbol{v}\in\boldsymbol{V}$ and, hence, identifies the directed edges pointing to $\boldsymbol{v}$. We let $\mathcal{G}$ be acyclic, i.e., there is no chain $\{\boldsymbol{v}_{i_{1}}\to\boldsymbol{v}_{i_{2}}\to\cdots\to\boldsymbol{v}_{i_{t}}\}$ of elements of $\boldsymbol{V}$ such that $\boldsymbol{v}_{i_{j}}\in\text{Pa}[{\boldsymbol{v}_{i_{j+1}}}]$ and $\boldsymbol{v}_{i_{j+1}}\in\text{Pa}[{\boldsymbol{v}_{i_{1}}}]$. Crucially, we assume that $\text{Pa}[{\boldsymbol{v}}]\subset\boldsymbol{A}$ for all $\boldsymbol{v}\in\boldsymbol{V}$, i.e., that only reference nodes have children, to distinguish the reference nodes $\boldsymbol{A}$ from the other nodes $\boldsymbol{B}$. Apart from the assumption that $\boldsymbol{a}_{j}\in\text{Pa}[{\boldsymbol{b}_{j}}]$, we refrain from defining the parents of a node, thereby retaining flexibility. In general, however, all locations in ${\cal U}_{j}$ will share the same parent set. In Section 4 we will consider meshes associated to domain tessellations.

Consider the enumeration ${\cal S}_{i}=\{\boldsymbol{s}_{i_{1}},\dots,\boldsymbol{s}_{i_{n_{i}}}\}$, where $\{i_{1},i_{2},\ldots,i_{n_{i}}\}\subset\{1,2,\ldots,n_{{\cal S}}\}$, and let $\boldsymbol{w}_{i}=(\boldsymbol{w}(\boldsymbol{s}_{i_{1}})^{\top},\boldsymbol{w}(\boldsymbol{s}_{i_{2}})^{\top},\ldots,\boldsymbol{w}(\boldsymbol{s}_{i_{n_{i}}})^{\top})^{\top}$ be the $n_{i}q\times 1$ random vector listing elements of $\boldsymbol{w}(\boldsymbol{s})$ for each $\boldsymbol{s}\in{\cal S}_{i}$. We now rewrite (1) as a product of $M$ conditional densities

$\displaystyle p(\boldsymbol{w}_{{\cal S}})$

$\displaystyle=p(\boldsymbol{w}_{1},\boldsymbol{w}_{2},\ldots,\boldsymbol{w}_{M})=\prod_{i=1}^{M}p(\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{1},\ldots,\boldsymbol{w}_{i-1}).$

(5)

The conditioning sets are then reduced based on the graph ${\cal G}$:

$\displaystyle\widetilde{p}(\boldsymbol{w}_{{\cal S}})$

$\displaystyle=\prod_{i=1}^{M}p(\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[{i}]})\;,$

(6)

where we denote $\boldsymbol{w}_{[{i}]}=\{\boldsymbol{w}_{j}:\boldsymbol{a}_{j}\in\text{Pa}[{\boldsymbol{a}_{i}}]\}$, and $\text{Pa}[{\boldsymbol{a}_{i}}]\subset\{\boldsymbol{a}_{1},\dots,\boldsymbol{a}_{i-1}\}\subset\boldsymbol{A}$. This is a proper multivariate joint density since the graph is acyclic (Lauritzen,, 1996). It is instructive to note how the above approximation behaves when the size of the parent set shrinks, for a given domain partitioning scheme. To this end, we adapt a result in Banerjee, (2020) and show that sparser DAGs correspond to a larger Kullback-Leibler (KL) divergence from the base density $p$. This result has been proved earlier for Gaussian likelihoods by Guinness, (2018), but the argument given below is free of distributional assumptions and is linked to the submodularity of entropy and the “information never hurts” principle (see e.g. Cover and Thomas,, 1991).

Consider random vector $\boldsymbol{w}$ and some partition of the domain ${\cal P}$ corresponding to nodes $\boldsymbol{V}=\{\boldsymbol{v}_{1},\dots,\boldsymbol{v}_{M}\}$ via map $\eta$. Let the base process correspond to graph ${\cal G}_{0}=\{\boldsymbol{V},\boldsymbol{E}_{0}\}$ where $\boldsymbol{E}_{0}=\{\text{Pa}_{0}[{\boldsymbol{v}_{1}}],\dots,\text{Pa}_{0}[{\boldsymbol{v}_{M}}]\}$. Then, let ${\cal G}_{1}=\{\boldsymbol{V},\boldsymbol{E}_{1}\}$ where $\boldsymbol{E}_{1}=\{\text{Pa}_{1}[{\boldsymbol{v}_{1}}],\dots,\text{Pa}_{1}[{\boldsymbol{v}_{M}}]\}$ and $\text{Pa}_{1}[{\boldsymbol{v}_{i}}]\subseteq\text{Pa}_{0}[{\boldsymbol{v}_{i}}]$ for all $i\in\{1,\dots,M\}$. Finally construct ${\cal G}_{2}=\{\boldsymbol{V},\boldsymbol{E}_{2}\}$ by letting $\text{Pa}_{2}[{\boldsymbol{v}_{i^{*}}}]=\text{Pa}_{1}[{\boldsymbol{v}_{i^{*}}}]\setminus\{\boldsymbol{v}^{*}\}$ for some $\boldsymbol{v}^{*}\in\text{Pa}_{1}[{\boldsymbol{v}_{i^{*}}}]$. In other words, graph ${\cal G}_{2}$ is obtained by removing the directed edge $\boldsymbol{v}^{*}\to\boldsymbol{v}_{i^{*}}$ from ${\cal G}_{1}$. We approximate $p$ using densities $p_{1}$ and $p_{2}$ based on ${\cal G}_{1}$ and ${\cal G}_{2}$, respectively, obtaining

$$\frac{p_{1}(\boldsymbol{w})}{p_{2}(\boldsymbol{w})}=\prod_{i=1}^{M}\frac{p(\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[{i}]_{1}})}{p(\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[{i}]_{2}})}=\frac{p(\boldsymbol{w}_{i^{*}}\,|\,\boldsymbol{w}_{[{i^{*}}]_{1}})}{p(\boldsymbol{w}_{i^{*}}\,|\,\boldsymbol{w}_{[{i^{*}}]_{2}})}.$$

(7)

Considering the Kullback-Leibler divergence of each density from $p$, and denoting $\boldsymbol{V}^{*}=\boldsymbol{V}\setminus\{\{i^{*}\}\cup\text{Pa}_{1}[{i^{*}}]\}$, we find

	$\displaystyle KL(p_{2}\|p)-KL(p_{1}$	$\displaystyle\|p)=\int\left\{\log\left(\frac{p(\boldsymbol{w})}{p_{2}(\boldsymbol{w})}\right)-\log\left(\frac{p(\boldsymbol{w}}{p_{1}(\boldsymbol{w}}\right)\right\}p(\boldsymbol{w})d\boldsymbol{w}$		(8)
		$\displaystyle=\int\log\left(\frac{p_{1}(\boldsymbol{w})}{p_{2}(\boldsymbol{w})}\right)p(\boldsymbol{w})d\boldsymbol{w}=\int\log\left(\frac{p(\boldsymbol{w}_{i^{*}}\,|\,\boldsymbol{w}_{[{i^{*}}]_{1}})}{p(\boldsymbol{w}_{i^{*}}\,|\,\boldsymbol{w}_{[{i^{*}}]_{2}})}\right)p(\boldsymbol{w})d\boldsymbol{w}$
		$\displaystyle=\int\log\left(\frac{p(\boldsymbol{w}_{i^{*}}\,|\,\boldsymbol{w}_{[{i^{*}}]_{1}})}{p(\boldsymbol{w}_{i^{*}}\,|\,\boldsymbol{w}_{[{i^{*}}]_{2}})}\right)p(\boldsymbol{w}_{i^{*}},\boldsymbol{w}_{[{i^{*}}]_{1}})d\boldsymbol{w}_{i^{*}}d\boldsymbol{w}_{[{i^{*}}]_{1}}$
		$\displaystyle=\int\left\{\int\log\left(\frac{p(\boldsymbol{w}_{i^{*}}\,|\,\boldsymbol{w}_{[{i^{*}}]_{1}})}{p(\boldsymbol{w}_{i^{*}}\,|\,\boldsymbol{w}_{[{i^{*}}]_{2}})}\right)p(\boldsymbol{w}_{i^{*}}\,|\,\boldsymbol{w}_{[{i^{*}}]_{1}})d\boldsymbol{w}_{i^{*}}\right\}p(\boldsymbol{w}_{[{i^{*}}]_{1}})d\boldsymbol{w}_{[{i^{*}}]_{1}}\geq 0,$

where we use (7), the fact that $\boldsymbol{V}^{*}$ and $\{i^{*}\}\cup\text{Pa}_{1}[{i^{*}}]$ are disjoint, and Jensen’s inequality. This result implies that larger parent sets are preferrable as they correspond to better approximations to the full model; the choice of sparser graphs will be driven by computational considerations – see Section 3.2.

We construct the spatial process over arbitrary locations by enumerating other locations as ${\cal U}=\{\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{n_{{\cal U}}}\}\subset{\cal D}\setminus{\cal S}$ and extending (6) to the non-reference locations. Given the partition of ${\cal U}$ defined earlier with components ${\cal U}_{j}$ for $j=1,2,\ldots,M$, for each $\boldsymbol{u}\in{\cal U}_{j}$ we set $\eta(\boldsymbol{u})=\boldsymbol{b}_{j}$ and recall that $\text{Pa}[{\boldsymbol{b}_{i}}]\subset\boldsymbol{A}$ by construction. For each $i=1,\dots,n_{{\cal U}}$, we denote $\boldsymbol{w}_{[{\boldsymbol{u}_{i}}]}=\{\boldsymbol{w}_{j}:\boldsymbol{a}_{j}\in\text{Pa}[{\eta(\boldsymbol{u}_{i})}]\}\subset\boldsymbol{w}_{{\cal S}}$ and define the conditional density of $\boldsymbol{w}_{{\cal U}}$ given $\boldsymbol{w}_{{\cal S}}$ as

$\displaystyle\widetilde{p}(\boldsymbol{w}_{{\cal U}}\mid\boldsymbol{w}_{{\cal S}})$

$\displaystyle=\prod_{\boldsymbol{u}_{i}\in{\cal U}}p(\boldsymbol{w}({\boldsymbol{u}_{i}})\mid\boldsymbol{w}_{[{\boldsymbol{u}_{i}}]})=\prod_{j=1}^{M}p(\boldsymbol{w}_{{\cal U}_{j}}\mid\boldsymbol{w}_{[{\boldsymbol{b}_{j}}]}).$

(9)

Therefore, for any finite subset of spatial locations $\mathcal{L}\subset\mathcal{D}$ we can let ${\cal U}=\mathcal{L}\setminus{\cal S}$ and obtain

$\displaystyle\widetilde{p}(\boldsymbol{w}_{\mathcal{L}})$

$\displaystyle=\int\widetilde{p}(\boldsymbol{w}_{{\cal U}}\mid\boldsymbol{w}_{{\cal S}})\widetilde{p}(\boldsymbol{w}_{{\cal S}})\prod_{\boldsymbol{s}_{i}\in{\cal S}\setminus\mathcal{L}}d(\boldsymbol{w}(\boldsymbol{s}_{i}))\;.$

We show (see Appendix A, available online) that this is a well-defined process by verifying the Kolmogorov consistency conditions. This new process can be built starting from a base process, a fixed reference set, domain partition ${\cal P}$ and a graph ${\cal G}$. Next, we elucidate with Gaussian processes.

Let $\{\boldsymbol{w}(\boldsymbol{\ell}):\boldsymbol{\ell}\in{\cal D}\}$ be a $q$-variate multivariate Gaussian process, denoted as $\boldsymbol{w}(\boldsymbol{\ell})\sim GP(\mathbf{0},\boldsymbol{C}(\cdot,\cdot\mid\boldsymbol{\theta}))$. The cross-covariance $\boldsymbol{C}(\cdot,\cdot\mid\boldsymbol{\theta})$ indexed by parameters $\boldsymbol{\theta}$ is a function $\boldsymbol{C}:{\cal D}\times{\cal D}\rightarrow{\cal M}_{q\times q}$, where ${\cal M}_{q\times q}$ is a subset of $\Re^{q\times q}$ (the space of all $q\times q$ real matrices) such that the $(i,j)$-th entry of $\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}\mid\boldsymbol{\theta})$ evaluates the covariance between the $i$-th and $j$-th elements of $\boldsymbol{w}(\boldsymbol{\ell})$ at $\boldsymbol{\ell}$ and $\boldsymbol{\ell}^{\prime}$, respectively, i.e., $\mbox{cov}(w_{i}(\boldsymbol{\ell}),w_{j}(\boldsymbol{\ell}^{\prime}))$. We omit dependence on $\boldsymbol{\theta}$ to simplify notation. The cross-covariance function itself needs to be neither symmetric nor positive-definite, but must satisfy the following two properties: (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 Genton and Kleiber, (2015) for a review of cross-covariance functions for multivariate processes. The (partial) realization of the multivariate process over any finite set ${\cal L}$ 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}}$.

We construct the MGP from a base, or parent, multivariate GP for $\boldsymbol{w}(\boldsymbol{\ell})$ and then, using the graph ${\cal G}$ defined in Section 2, represent the joint density at the reference set ${\cal S}$ as

$\displaystyle\widetilde{p}(\boldsymbol{w}_{{\cal S}})$

$\displaystyle=\prod_{j=1}^{M}N(\boldsymbol{w}_{j}\mid\boldsymbol{H}_{j}\boldsymbol{w}_{[{j}]},\boldsymbol{R}_{j}),$

(10)

where $\boldsymbol{H}_{1}=\boldsymbol{O}_{n_{1}\times 1}$, $\boldsymbol{R}_{1}=\boldsymbol{C}_{{\cal S}_{j}}$ and for $j>1$, $\boldsymbol{H}_{j}=\boldsymbol{C}_{{\cal S}_{j},{\cal S}_{[{j}]}}\boldsymbol{C}^{-1}_{{\cal S}_{[{j}]}}$ and $\boldsymbol{R}_{j}=\boldsymbol{C}_{{\cal S}_{j}}-\boldsymbol{C}_{{\cal S}_{j},{\cal S}_{[{j}]}}\boldsymbol{C}^{-1}_{{\cal S}_{[{j}]}}\boldsymbol{C}_{{\cal S}_{[{j}]},{\cal S}_{j}}$. The resulting joint density $\widetilde{p}(\boldsymbol{w}_{{\cal S}})$ is multivariate normal with covariance $\widetilde{\boldsymbol{C}}_{{\cal S}}$ and a precision matrix $\widetilde{\boldsymbol{C}}^{-1}_{{\cal S}}$. The precision matrix for Gaussian graphical models is easily derived using customary linear model representations for each conditional regression. Consider the DAG in (6). Each $\boldsymbol{w}_{i}$ is $n_{i}q\times 1$ and let $J_{i}=|\text{Pa}[{\boldsymbol{a}_{i}}]|$ be the number of parents for $\boldsymbol{a}_{i}$ in the graph ${\cal G}$. Furthermore, let $\boldsymbol{C}_{i,j}$ be the $n_{i}q\times n_{j}q$ covariance matrix between $\boldsymbol{w}_{i}$ and $\boldsymbol{w}_{j}$, $\boldsymbol{C}_{i,[i]}$ be the $n_{i}q\times J_{i}q$ covariance matrix between $\boldsymbol{w}_{i}$ and $\boldsymbol{w}_{[i]}$, and $\boldsymbol{C}_{[i],[i]}$ be the $J_{i}q\times J_{i}q$ covariance matrix between $\boldsymbol{w}_{[{i}]}$ and itself. Representing each conditional density in (6) as a linear regression on $\boldsymbol{w}_{i}$, we get

$$\boldsymbol{w}_{1}=\boldsymbol{\omega}_{1}\sim N(\boldsymbol{0},\boldsymbol{R}_{1})\;;\quad\boldsymbol{w}_{i}=\sum_{\{j:\boldsymbol{a}_{j}\in\text{Pa}[{\boldsymbol{a}_{i}}]\}}\boldsymbol{H}_{ij}\boldsymbol{w}_{j}+\boldsymbol{\omega}_{i}\;,\;\;i=2,3,\ldots,M\;,$$

(11)

where each $\boldsymbol{H}_{ij}$ is an $n_{i}q\times n_{j}q$ is a coefficient matrix representing the multivariate regression of $\boldsymbol{w}_{j}$ given $\boldsymbol{w}_{[i]}$, $\boldsymbol{\omega}_{i}\stackrel{{\scriptstyle ind}}{{\sim}}N(\boldsymbol{0},\boldsymbol{R}_{i})$ for $i=1,2,\ldots,M$, and each $\boldsymbol{R}_{i}$ is an $n_{i}q\times n_{i}q$ residual covariance matrix. We set $\boldsymbol{H}_{ii}=\boldsymbol{O}$ and $H_{ij}=\boldsymbol{O}$, where $\boldsymbol{O}$ is the matrix of zeros, whenever $j\in\{j:\boldsymbol{a}_{j}\notin\text{Pa}[{\boldsymbol{a}_{i}}]\}$. For $j\in\{j:\boldsymbol{a}_{j}\in\text{Pa}[{\boldsymbol{a}_{i}}]\}$, let $\{j_{1},j_{2},\ldots,j_{J_{i}}\}$ be the indices in $\text{Pa}[{\boldsymbol{a}_{i}}]$ and let $\boldsymbol{H}_{i,[i]}=\left[\boldsymbol{H}_{i,j_{1}},\boldsymbol{H}_{i,j_{2}},\ldots,\boldsymbol{H}_{i,j_{J_{i}}}\right]$ be the $n_{i}q\times(\sum_{k=1}^{J_{i}}n_{j_{k}})q$ block matrix formed by stacking $\boldsymbol{H}_{i,j_{k}}$ side by side for each $\boldsymbol{a}_{j_{k}}\in\text{Pa}[{\boldsymbol{a}_{i}}]$. Since $\mbox{E}[\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[i]}]=\boldsymbol{H}_{i,[i]}\boldsymbol{w}_{[i]}=\boldsymbol{C}_{i,[i]}\boldsymbol{C}_{[i][i]}^{-1}\boldsymbol{w}_{[i]}$, we obtain $\boldsymbol{H}_{i,[i]}=\boldsymbol{C}_{i,[i]}\boldsymbol{C}_{[i][i]}^{-1}$ and each $\boldsymbol{H}_{ij_{K}}$ can be obtained from the respective submatrix of $\boldsymbol{H}_{i[i]}$. We also obtain $\boldsymbol{R}_{i}=\mbox{var}\{\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[i]}\}=\boldsymbol{C}_{i,i}-\boldsymbol{C}_{i,[i]}\boldsymbol{C}_{[i][i]}^{-1}\boldsymbol{C}_{[i],i}$. Therefore, all the $\boldsymbol{H}_{ij}$’s and $\boldsymbol{R}_{i}$’s can be computed from the base cross-covariance function.

The distribution of $\boldsymbol{w}=[\boldsymbol{w}_{1}^{\top},\boldsymbol{w}_{2}^{\top},\ldots,\boldsymbol{w}_{M}^{\top}]^{\top}$ can be obtained by noting that $\boldsymbol{w}=\boldsymbol{H}\boldsymbol{w}+\boldsymbol{\omega}$, where $\boldsymbol{H}=\{\boldsymbol{H}_{ij}\}$ is the $(\sum_{i=1}^{M}n_{i}q)\times(\sum_{i=1}^{M}n_{i}q)$ block matrix with $\{\boldsymbol{H}_{ij}\}$ as $(i,j)$-th block. Therefore, $\widetilde{\boldsymbol{C}}_{{\cal S}}=\mbox{var}(\boldsymbol{w})=(\boldsymbol{I}-\boldsymbol{H})^{-1}\boldsymbol{R}(\boldsymbol{I}-\boldsymbol{H})^{-\top}$, where $\boldsymbol{R}$ is block-diagonal with $\boldsymbol{R}_{i}$ as the $(i,i)$-th block. Note that $\boldsymbol{I}-\boldsymbol{H}$ is block lower-triangular with $1$’s on the diagonal, hence non-singular. Also, the precision matrix $\widetilde{\boldsymbol{C}}^{-1}_{{\cal S}}=(\boldsymbol{I}-\boldsymbol{H})^{\top}\boldsymbol{R}^{-1}(\boldsymbol{I}-\boldsymbol{H})$ is sparse because of $\boldsymbol{H}_{ij}=\boldsymbol{O}$ whenever $\boldsymbol{a}_{j}\notin\text{Pa}[{\boldsymbol{a}_{i}}]$. Block-sparsity of $\widetilde{\boldsymbol{C}}^{-1}_{{\cal S}}$ can be induced by building $\mathcal{G}$ with few, carefully placed directed edges among nodes in $\boldsymbol{A}$; Appendix B, available online, contains a more in-depth treatment. We extend (10) to the collection of non-reference locations ${\cal U}\subset{\cal D}\setminus{\cal S}$:

$\displaystyle\widetilde{p}(\boldsymbol{w}_{{\cal U}}\mid\boldsymbol{w}_{{\cal S}})$

$\displaystyle=\prod_{j=1}^{M}N(\boldsymbol{w}_{{\cal U}_{j}}\mid\boldsymbol{H}_{{\cal U}_{j}}\boldsymbol{w}_{[{\boldsymbol{b}_{j}}]},\boldsymbol{R}_{{\cal U}_{j}})=N(\boldsymbol{w}_{{\cal U}}\mid\boldsymbol{H}_{{\cal U}}\boldsymbol{w}_{{\cal S}},\boldsymbol{R}_{{\cal U}}),$

(12)

where $\boldsymbol{H}_{{\cal U}_{j}}=\boldsymbol{C}_{{\cal U}_{j},{\cal S}_{[{\boldsymbol{b}_{j}}]}}\boldsymbol{C}^{-1}_{{\cal S}_{[{\boldsymbol{b}_{j}}]}}$ and $\boldsymbol{R}_{{\cal U}_{j}}=\boldsymbol{C}_{{\cal U}_{j}}-\boldsymbol{C}_{{\cal U}_{j},{\cal S}_{[{\boldsymbol{b}_{j}}]}}\boldsymbol{C}^{-1}_{{\cal S}_{[{\boldsymbol{b}_{j}}]}}\boldsymbol{C}_{{\cal S}_{[{\boldsymbol{b}_{j}}]},{\cal U}_{j}}$, analogously to (10), while $\boldsymbol{H}_{{\cal U}}$ and $\boldsymbol{R}_{{\cal U}}$ are analogous to $\boldsymbol{H}_{{\cal S}}$ and $\boldsymbol{R}_{{\cal S}}$. Clearly, given that all the $\widetilde{p}$ densities are Gaussian, all finite dimensional distributions will also be Gaussian. We have constructed a Gaussian process with the following cross-covariance function for any two locations $\boldsymbol{\ell}_{1},\boldsymbol{\ell}_{2}\in{\cal D}$

$\displaystyle Cov_{\widetilde{p}}(\boldsymbol{w}(\boldsymbol{\ell}_{1}),\boldsymbol{w}(\boldsymbol{\ell}_{2}))$

$\displaystyle=\begin{cases}\widetilde{\boldsymbol{C}}_{\boldsymbol{s}_{i},\boldsymbol{s}_{j}}&\text{if\ }\boldsymbol{\ell}_{1}=\boldsymbol{s}_{i},\boldsymbol{\ell}_{2}=\boldsymbol{s}_{j}\text{ and }\boldsymbol{s}_{i},\boldsymbol{s}_{j}\in{\cal S}\\ \boldsymbol{H}_{\boldsymbol{\ell}_{1}}\widetilde{\boldsymbol{C}}_{{\cal S}_{[{\boldsymbol{\ell}_{1}}]},\boldsymbol{s}_{j}}&\text{if\ }\boldsymbol{\ell}_{1}\in{\cal D}\setminus{\cal S},\boldsymbol{\ell}_{2}=\boldsymbol{s}_{j}\text{ and }\boldsymbol{s}_{j}\in{\cal S}\\ \delta_{(\boldsymbol{\ell}_{1}=\boldsymbol{\ell}_{2})}\boldsymbol{R}_{\boldsymbol{\ell}_{1}}+\boldsymbol{H}_{\boldsymbol{\ell}_{1}}\widetilde{\boldsymbol{C}}_{{\cal S}_{[{\boldsymbol{\ell}_{1}}]},{\cal S}_{[{\boldsymbol{\ell}_{2}}]}}\boldsymbol{H}_{\boldsymbol{\ell}_{2}}^{\top}&\text{otherwise.}\end{cases}$

For a given base Gaussian covariance function $\boldsymbol{C}$, domain partitioning ${\cal P}$, mesh ${\cal G}$, and reference set ${\cal S}$, we denote the corresponding meshed Gaussian process as $\text{MGP}({\cal G},{\cal P},{\cal S},\boldsymbol{C})$.