# Spatial Multivariate Trees
for Big Data Bayesian Regression

Michele Peruzzi111Department of Statistical Science, Duke University David B.
Dunson††footnotemark:

###### Abstract

High resolution geospatial data are challenging because standard
geostatistical models based on Gaussian processes are known to not scale to
large data sizes. While progress has been made towards methods that can be
computed more efficiently, considerably less attention has been devoted to big
data methods that allow the description of complex relationships between
several outcomes recorded at high resolutions by different sensors. Our
Bayesian multivariate regression models based on spatial multivariate trees
(SpamTrees) achieve scalability via conditional independence assumptions on
latent random effects following a treed directed acyclic graph. Information-
theoretic arguments and considerations on computational efficiency guide the
construction of the tree and the related efficient sampling algorithms in
imbalanced multivariate settings. In addition to simulated data examples, we
illustrate SpamTrees using a large climate data set which combines satellite
data with land-based station data. Source code is available at
github.com/mkln/spamtree.

Keywords: Directed acyclic graph, Gaussian process, Geostatistics,
Multivariate regression, Markov chain Monte Carlo, Multiscale/multiresolution.

##  1 Introduction

It is increasingly common in the natural and social sciences to amass large
quantities of geo-referenced data. Researchers seek to use these data to
understand phenomena and make predictions via interpretable models that
quantify uncertainty taking into account the spatial and temporal dimensions.
Gaussian processes (GP) are flexible tools that can be used to characterize
spatial and temporal variability and quantify uncertainty, and considerable
attention has been devoted to developing GP-based methods that overcome their
notoriously poor scalability to large data. The literature on scaling GPs to
big data is now extensive. We mention low-rank methods (Quiñonero-Candela and
Rasmussen, , 2005; Snelson and Ghahramani, , 2007; Banerjee et al., , 2008;
Cressie and Johannesson, , 2008); their extensions (Low et al., , 2015;
Ambikasaran et al.,, 2016; Huang and Sun, , 2018; Geoga et al., , 2020);
methods that exploit special structure or simplify the representation of
multidimensional inputs—for instance, a Toeplitz structure of the covariance
matrix scales GPs to big time series data, and tensor products of scalable
univariate kernels can be used for multidimensional inputs (Gilboa et al., ,
2015; Moran and Wheeler, , 2020; Loper et al., , 2020). These methods may be
unavailable or perform poorly in geostatistical settings, which focus on
small-dimensional inputs, i.e. the spatial coordinates plus time. In these
scenarios, low-rank methods oversmooth the spatial surface (Banerjee et al., ,
2010), Toeplitz-like structures are typically absent, and so-called separable
covariance functions obtained via tensor products poorly characterize spatial
and temporal dependence. To overcome these hurdles, one can use covariance
tapering and domain partitioning (Furrer et al., , 2006; Kaufman et al., ,
2008; Sang and Huang, , 2012; Stein, , 2014; Katzfuss, , 2017) or composite
likelihood methods and sparse precison matrix approximations (Vecchia, , 1988;
Rue and Held, , 2005; Eidsvik et al., , 2014); refer to Sun et al.,  (2011),
Banerjee,  (2017), Heaton et al.,  (2019) for reviews of scalable
geostatistical methods.

Figure 1: Observed data of Section 4.2. Missing outcomes are in magenta. GHCN
data are much more sparsely observed compared to satellite imaging from MODIS.

Additional difficulties arise in multivariate (or multi-output) regression
settings. Multivariate geostatistical data are commonly misaligned, i.e.
observed at non-overlapping spatial locations (Gelfand et al., , 2010). Figure
1 shows several variables measured at non-overlapping locations, with one
measurement grid considerably sparser than the others. This issue can be
solved by modeling cross-dependence of the outputs via latent spatial random
effects thought of as a realization of an underlying GP and embedded in a
larger hierarchical model.

Unfortunately, GP approximations that do not correspond to a valid stochastic
process may inaccurately characterize uncertainty, as the models used for
estimation and interpolation may not coincide. Rather than seeking
approximations to the full GP, one can develop valid standalone spatial
processes by introducing conditional independence across spatial locations as
prescribed by a sparse directed acyclic graph (DAG). These models are
advantageous because they lead to scalability by construction; in other words,
posterior computing algorithms for these methods can be interpreted as
approximate algorithms for the full GP, but also as exact algorithms for the
standalone process.

This family of method includes nearest-neighbor Gaussian processes, which
limit dependence to a small number of neighboring locations (NNGP; Datta et
al., 2016a ; Datta et al., 2016b ), and block-NNGPs (Quiroz et al., , 2019).
There is a close relation between DAG structure and computational performance
of NNGPs: some orderings may be associated to improved approximations
(Guinness, , 2018), and graph coloring algorithms (Molloy and Reed, , 2002;
Lewis, , 2016) can be used for parallel Gibbs sampling. Inferring ordering or
coloring can be problematic when data are in the millions, but these issues
can be circumvented by forcing DAGs with known properties onto the data; in
meshed GPs (MGPs; Peruzzi et al., , 2020), patterned DAGs associated to domain
tiling are associated to more efficient sampling of the latent effects.
Alternative so-called multiscale or multiresolution methods correspond to DAGs
with hierarchical node structures (trees), which are typically coupled with
recursive domain partitioning; in this case, too, efficiencies follow from the
properties of the chosen DAG. There is a rich literature on Gaussian processes
and recursive partitioning, see e.g Ferreira and Lee,  (2007); Gramacy and
Lee,  (2008); Fox and Dunson,  (2012); in geospatial contexts, in addition to
the GMRF-based method of Nychka et al.,  (2015), multi-resolution
approximations (MRA; Katzfuss, , 2017) replace an orthogonal basis
decomposition with approximations based on tapering or domain partitioning and
also have a DAG interpretation (Katzfuss and Guinness, , 2019).

Considerably less attention has been devoted to process-based methods that
ensure scalability in multivariate contexts, with the goal of modeling the
spatial and/or temporal variability of several variables jointly via flexible
cross-covariance functions (Genton and Kleiber, , 2015). When scalability of
GP methods is achieved via reductions in the conditioning sets, including more
distant locations is thought to aid in the estimation of unknown covariance
parameters (Stein et al., , 2004). However, the size of such sets may need to
be reduced excessively when outcomes are not of very small dimension. One
could restrict spatial coverage of the conditioning sets, but this works best
when data are not misaligned, in which case all conditioning sets will include
outcomes from all margins; this cannot be achieved for misaligned data,
leading to pathological behavior. Alternatively, one can model the
multivariate outcomes themselves as a DAG; however this may only work on a
case-by-case basis. Similarly, recursive domain partitioning strategies work
best for data that are measured uniformly in space as this guarantees
similarly sized conditioning sets; on the contrary, recursive partitioning
struggles in predicting the outcomes at large unobserved areas as they tend to
be associated to the small conditioning sets making up the coarser scales or
resolutions.

In this article, we solve these issues by introducing a Bayesian regression
model that encodes spatial dependence as a latent spatial multivariate tree
(SpamTree); conditional independence relations at the reference locations are
governed by the branches in a treed DAG, whereas a map is used to assign all
non-reference locations to leaf nodes of the same DAG. This assignment map
controls the nature and the size of the conditioning sets at all locations;
when severe restrictions on the reference set of locations become necessary
due to data size, this map is used to improve estimation and predictions and
overcome common issues in standard nearest-neighbor and recursive partition
methods while maintaining the desirable recursive properties of treed DAGs.
Unlike methods based on defining conditioning sets based solely on spatial
proximity, SpamTrees scale to large data sets without excessive reduction of
the conditioning sets. Furthermore, SpamTrees are less restrictive than
methods based on recursive partitioning and can be built to guarantee
similarly-sized conditioning sets at all locations.

The present work adds to the growing literature on spatial processes defined
on DAGs by developing a method that targets efficient computations of Bayesian
multivariate spatial regression models. SpamTrees share similarities with MRAs
(Katzfuss, , 2017); however, while MRAs are defined as a basis function
expansion, they can be represented by a treed graph of a SpamTree with full
“depth” as defined later (the DAG on the right of Figure 2), in univariate
settings, and “response” models. All these restrictions are relaxed in this
article. In considering spatial proximity to add “leaves” to our treed graph,
our methodology also borrows from nearest-neighbor methods (Datta et al.,
2016a, ). However, while we use spatial neighbors to populate the conditioning
sets for non-reference locations, the same cannot be said about reference
locations for which the treed graph is used instead. Our construction of the
SpamTree process also borrows from MGPs on tessellated domains (Peruzzi et
al., , 2020); however, the treed DAG we consider here induces markedly
different properties on the resulting spatial process owing to its recursive
nature. Finally, a contribution of this article is in developing self-
contained sampling algorithms which, based on the graphical model
representation of the model, will not require any external libraries.

The article builds SpamTrees as a standalone process based on a DAG
representation in Section 2. A Gaussian base process is considered in Section
3 and the resulting properties outlined, along with sampling algorithms.
Simulated data and real-world applications are in Section 4; we conclude with
a discussion in Section 5. The Appendix provides more in-depth treatment of
several topics and additional algorithms.

##  2 Spatial Multivariate Trees

Figure 2: Three SpamTrees on $M=4$ levels with depths $\delta=1$ (left),
$\delta=3$ (center), and $\delta=4$ (right). Nodes are represented by circles,
with branches colored in brown and leaves in green.

Consider a spatial or spatiotemporal domain ${\cal D}$. With the temporal
dimension, we have ${\cal D}\subset\Re^{d}\times[0,\infty)$, otherwise ${\cal
D}\subset\Re^{d}$. A $q$-variate spatial process is defined as an uncountable
set of random variables
$\\{\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$, paired with a
probability law $P$ defining the joint distribution of any finite sample from
that set. Let
$\\{\boldsymbol{\ell}_{1},\boldsymbol{\ell}_{2},\ldots,\boldsymbol{\ell}_{n_{{\cal
L}}}\\}={\cal L}\subset{\cal D}$ be of size $n_{{\cal L}}$. 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}$ has joint density $p(\boldsymbol{w}_{{\cal L}})$. After
choosing an arbitrary order of the locations, $p(\boldsymbol{w}_{{\cal
L}})=\prod_{i=1}^{n_{{\cal
L}}}p(\boldsymbol{w}(\boldsymbol{\ell}_{i})\,|\,\boldsymbol{w}(\boldsymbol{\ell}_{1}),\ldots,\boldsymbol{w}(\boldsymbol{\ell}_{i-1})),$
where the conditioning set for each $\boldsymbol{w}(\boldsymbol{\ell}_{i})$
can be interpreted as the set of nodes that have a directed edge towards
$\boldsymbol{w}(\boldsymbol{\ell}_{i})$ in a DAG. Some scalable spatial
processes result from reductions in size of the conditioning sets, following
one of several proposed strategies (Vecchia, , 1988; Stein et al., , 2004;
Gramacy and Apley, , 2015; Datta et al., 2016a, ; Katzfuss and Guinness, ,
2019; Peruzzi et al., , 2020). Accordingly,

$\displaystyle p(\boldsymbol{w}_{{\cal L}})$
$\displaystyle=\prod_{i=1}^{n_{{\cal
L}}}p(\boldsymbol{w}(\boldsymbol{\ell}_{i})\,|\,\boldsymbol{w}(\text{Pa}[{\boldsymbol{\ell}_{i}}])),$
(1)

where $\text{Pa}[{\boldsymbol{\ell}_{i}}]$ is the set of spatial locations
that correspond to directed edges pointing to $\boldsymbol{\ell}_{i}$ in the
DAG. If $\text{Pa}[{\boldsymbol{\ell}_{i}}]$ is of size $J$ or less for all
$i=1,\dots,n_{{\cal L}}$, then
$\boldsymbol{w}(\text{Pa}[{\boldsymbol{\ell}_{i}}]))$ is of size $Jq$. Methods
that rely on reducing the size of parent sets are thus negatively impacted by
the dimension $q$ of the multivariate outcome; if $q$ is not very small,
reducing the number of parent locations $J$ may be insufficient for scalable
computations. As an example, an NNGP model has
$\text{Pa}[{\boldsymbol{\ell}_{i}}]=N(\boldsymbol{\ell}_{i})$, where
$N(\cdot)$ maps a location in the spatial domain to its neighbor set. It is
customary in practice to consider $Jq=m\leq 20$ for accurate and scalable
estimation and predictions in univariate settings, but this may be restrictive
in some multivariate settings as one must reduce $J$ to maintain similar
computing times, possibly harming estimation and prediction accuracy.

We represent the $i$th component of the $q\times 1$ vector
$\boldsymbol{w}(\boldsymbol{\ell})$ as $w(\boldsymbol{\ell},\xi_{i})$, where
$\xi_{i}=(\xi_{i1},\dots,\xi_{ik})^{\top}\in\Xi$ for some $k$ and $\Xi$ serves
as the $k$-dimensional latent spatial domain of variables. The $q$-variate
process $\boldsymbol{w}(\boldsymbol{\ell})$ is thus recast as
$\\{w(\boldsymbol{\ell},\xi):(\boldsymbol{\ell},\xi)\in{\cal D}\times\Xi\\}$,
with $\xi$ representing the latent location in the domain of variables. We can
then write (1) as

$\displaystyle p(\boldsymbol{w}_{{\cal L}^{*}})$
$\displaystyle=\prod_{i=1}^{n_{{\cal
L}^{*}}}p(w(\boldsymbol{\ell}^{*}_{i})\,|\,w(\text{Pa}[{\boldsymbol{\ell}^{*}_{i}}])),$
(2)

where ${\cal L}^{*}=\\{\boldsymbol{\ell}^{*}_{i}\\}_{i=1}^{n_{{\cal L}^{*}}}$,
$\boldsymbol{\ell}^{*}_{i}\in{\cal D}\times\Xi={\cal D}^{*}$, and $w(\cdot)$
is a univariate process on the expanded domain ${\cal D}^{*}$. This
representation is useful as it provides a clearer accounting of the assumed
conditional independence structure of the process in a multivariate context.

###  2.1 Constructing spatial multivariate DAGs

SpamTrees are defined by (i) a treed DAG ${\cal G}$ with branches and leaves
on $M$ levels and with depth $\delta\leq M$; (ii) a reference set of locations
${\cal S}$; (iii) a cherry picking map. The graph is ${\cal
G}=\\{\boldsymbol{V},\boldsymbol{E}\\}$ where the nodes are
$\boldsymbol{V}=\\{\boldsymbol{v}_{1},\dots,\boldsymbol{v}_{m_{V}}\\}=\boldsymbol{A}\cup\boldsymbol{B}$,
$\boldsymbol{A}\cap\boldsymbol{B}=\emptyset$. The reference or branch nodes
are
$\boldsymbol{A}=\\{\boldsymbol{a}_{1},\dots,\boldsymbol{a}_{m_{A}}\\}=\boldsymbol{A}_{0}\cup\boldsymbol{A}_{1}\cup\dots\cup\boldsymbol{A}_{M-1}$,
where
$\boldsymbol{A}_{i}=\\{\boldsymbol{a}_{i,1},\dots,\boldsymbol{a}_{i,m_{i}}\\}$
for all $i=0,\dots,M-1$ and with
$\boldsymbol{A}_{i}\cap\boldsymbol{A}_{j}=\emptyset$ if $i\neq j$. The non-
reference or leaf nodes are
$\boldsymbol{B}=\\{\boldsymbol{b}_{1},\dots,\boldsymbol{b}_{m_{B}}\\}$,
$\boldsymbol{A}\cap\boldsymbol{B}=\emptyset$. We also denote
$\boldsymbol{V}_{r}=\boldsymbol{A}_{r}$ for $r=0,\dots,M-1$ and
$\boldsymbol{V}_{M}=\boldsymbol{B}$. The edges are
$\boldsymbol{E}=\\{\text{Pa}[{\boldsymbol{v}}]\subset\boldsymbol{V}:\boldsymbol{v}\in\boldsymbol{V}\\}$
and similarly
$\text{Ch}[{\boldsymbol{v}}]=\\{\boldsymbol{v}^{\prime}\in\boldsymbol{V}:\boldsymbol{v}\in\text{Pa}[{\boldsymbol{v}^{\prime}}]\\}$.
The reference set ${\cal S}$ is partitioned in $M$ levels starting from zero,
and each level is itself partitioned into reference subsets: ${\cal
S}=\cup_{r=0}^{M-1}{\cal S}_{r}=\cup_{r=0}^{M-1}\cup_{i=1}^{m_{i}}S_{ri}$,
where $S_{ri}\cap S_{r^{\prime}i^{\prime}}=\emptyset$ if $r\neq r^{\prime}$ or
$i\neq i^{\prime}$ and its complement set of non-reference or other locations
${\cal U}={\cal D}^{*}\setminus{\cal S}$. The cherry picking map is
$\eta:{\cal D}^{*}\rightarrow\boldsymbol{V}$ and assigns a node (and therefore
all the edges directed to it in ${\cal G}$) to any location in the domain,
following a user-specified criterion.

####  2.1.1 Branches and leaves

For a given $M$ and a depth $\delta\leq M$, we impose a treed structure on
${\cal G}$ by assuming that if $\boldsymbol{v}\in\boldsymbol{A}_{i}$ and
$i>M-\delta=M_{\delta}$ then there exists a sequence of nodes
$\\{\boldsymbol{v}_{r_{M_{\delta}}},\dots,\boldsymbol{v}_{r_{i-1}}\\}$ such
that $\boldsymbol{v}_{r_{j}}\in\boldsymbol{A}_{j}$ for
$j=M_{\delta},\dots,i-1$ and
$\text{Pa}[{\boldsymbol{v}}]=\\{\boldsymbol{v}_{r_{M_{\delta}}},\boldsymbol{v}_{r_{1}},\dots,\boldsymbol{v}_{r_{j-1}}\\}$.
If $i\leq M-\delta=M_{\delta}$ then
$\text{Pa}[{\boldsymbol{v}}]=\\{\boldsymbol{v}_{i-1}\\}$ with
$\boldsymbol{v}_{i-1}\in\boldsymbol{A}_{i-1}$. $\boldsymbol{A}_{0}$ is the
tree root and is such that $\text{Pa}[{\boldsymbol{v}_{0}}]=\emptyset$ for all
$\boldsymbol{v}_{0}\in\boldsymbol{A}_{0}$. The depth $\delta$ determines the
number of levels of ${\cal G}$ (from the top) across which the parent sets are
nested. Choosing $\delta=1$ implies that all nodes have a single parent;
choosing $\delta=M$ implies fully nested parent sets (i.e. if
$\boldsymbol{v}_{i}\in\text{Pa}[{\boldsymbol{v}_{j}}]$ then
$\text{Pa}[{\boldsymbol{v}_{i}}]\subset\text{Pa}[{\boldsymbol{v}_{j}}]$ for
all $\boldsymbol{v}_{i},\boldsymbol{v}_{j}\in\boldsymbol{V}$). The $m_{i}$
elements of $\boldsymbol{A}_{i}$ are the branches at level $i$ of ${\cal G}$
and they have $i-M_{\delta}$ parents if the current level $i$ is above the
depth level $M_{\delta}$ and 1 parent otherwise. We refer to terminal branches
as nodes $\boldsymbol{v}\in\boldsymbol{A}$ such that
$\text{Ch}[{\boldsymbol{v}}]\subset\boldsymbol{B}$. For all choices of
$\delta$,
$\boldsymbol{v}\in\boldsymbol{A}_{i},\boldsymbol{v}^{\prime}\in\boldsymbol{A}_{j}$
and $\boldsymbol{v}\in\text{Pa}[{\boldsymbol{v}^{\prime}}]$ implies $i<j$;
this guarantees acyclicity.

As for the leaves, for all $\boldsymbol{v}\in\boldsymbol{B}$ we assume
$\text{Pa}[{\boldsymbol{v}}]=\\{\boldsymbol{v}_{r_{M_{\delta}}},\dots,\boldsymbol{v}_{r_{k}}\\}$
for some integer sequence $\\{r_{M_{\delta}},\dots,r_{k}\\}$ and
$\boldsymbol{v}_{r_{i}}\in\boldsymbol{A}_{i}$ with $i\geq M_{\delta}$. We
allow the existence of multiple leaves with the same parent set, i.e. there
can be $k$ and $\boldsymbol{b}_{i_{1}},\dots,\boldsymbol{b}_{i_{k}}$ such that
for all $i_{2},\dots,i_{k}$,
$\text{Pa}[{\boldsymbol{b}_{i_{h}}}]=\text{Pa}[{\boldsymbol{b}_{i_{1}}}]$.
Acyclicity of ${\cal G}$ is maintained as leaves are assumed to have no
children. Figure 2 represents the graph associated to SpamTrees with different
depths.

####  2.1.2 Cherry picking via $\eta(\cdot)$

The link between ${\cal G}$, ${\cal S}$ and ${\cal U}$ is established via the
map $\eta:{\cal D}^{*}\rightarrow\boldsymbol{V}$ which associates a node in
${\cal G}$ to any location $\boldsymbol{\ell}^{*}$ in the expanded domain
${\cal D}^{*}$:

$\displaystyle\eta(\boldsymbol{\ell}^{*})$
$\displaystyle=\left\\{\begin{array}[]{l}\eta_{A}(\boldsymbol{\ell}^{*})=\boldsymbol{a}_{ri}\in\boldsymbol{A}_{r}\;\mbox{
if }\;\boldsymbol{\ell}^{*}\in S_{ri},\\\
\eta_{B}(\boldsymbol{\ell}^{*})=\boldsymbol{b}\in\boldsymbol{B}\;\mbox{ if
}\;\boldsymbol{\ell}^{*}\in{\cal U}.\end{array}\right.$ (5)

This is a many-to-one map; note however that all locations in $S_{ij}$ are
mapped to $\boldsymbol{a}_{ij}$: by calling
$\eta(X)=\\{\eta(\boldsymbol{\ell}^{*}):\boldsymbol{\ell}^{*}\in X\\}$ then
for any $i=0,\dots,M-1$ and any $j=1,\dots,m_{i}$ we have
$\eta(S_{ij})=\eta_{A}(S_{ij})=\boldsymbol{a}_{ij}$. SpamTrees introduce
flexibility by cherry picking the leaves, i.e. using $\eta_{B}:{\cal
U}\rightarrow\boldsymbol{B}$, the restriction of $\eta$ to ${\cal U}$. Since
each leaf node $\boldsymbol{b}_{j}$ determines a unique path in ${\cal G}$
ending in $\boldsymbol{b}_{j}$, we use $\eta_{B}$ to assign a convenient
parent set to $w(\boldsymbol{u})$, $\boldsymbol{u}\in{\cal U}$, following some
criterion.

For example, suppose that $\boldsymbol{u}=(\boldsymbol{\ell},\xi_{s})$ meaning
that $w(\boldsymbol{u})=w(\boldsymbol{\ell},\xi_{s})$ is the realization of
the $s$-th variable at the spatial location $\boldsymbol{\ell}$, and we wish
to ensure that $\text{Pa}[{w(\boldsymbol{u})}]$ includes realizations of the
same variable. Denote
$\boldsymbol{T}=\\{\boldsymbol{v}\in\boldsymbol{A}:\text{Ch}[{\boldsymbol{v}}]\subset\boldsymbol{B}\\}$
as the set of terminal branches of ${\cal G}$. Then we find
$(\boldsymbol{\ell},\xi_{s})_{\text{opt}}=\arg\min_{(\boldsymbol{\ell}^{\prime},\xi^{\prime}=\xi_{s})\in\eta_{A}^{-1}(\boldsymbol{T})}d(\boldsymbol{\ell}^{\prime},\boldsymbol{\ell})$
where $d(\cdot,\cdot)$ is the Euclidean distance. Since
$(\boldsymbol{\ell},\xi_{s})_{\text{opt}}\in S_{ij}$ for some $i,j$ we have
$\eta_{A}((\boldsymbol{\ell},\xi_{s})_{\text{opt}})=\boldsymbol{a}_{ij}$. We
then set $\eta_{B}(\boldsymbol{u})=\boldsymbol{b}_{k}$ where
$\text{Pa}[{\boldsymbol{b}_{k}}]=\\{\boldsymbol{a}_{ij}\\}$. In a sense
$\boldsymbol{a}_{ij}$ is the terminal node nearest to $\boldsymbol{u}$; having
defined $\eta_{B}$ in such a way forces the parent set of any location to
include at least one realization of the process from the same variable. There
is no penalty in using ${\cal D}^{*}={\cal D}\times\Xi$ as we can write
$p(\boldsymbol{w}(\boldsymbol{u})\,|\,\text{Pa}[{\boldsymbol{w}(\boldsymbol{u})}])=p(\boldsymbol{w}((\boldsymbol{\ell},\xi_{1}),\dots,(\boldsymbol{\ell},\xi_{q}))\,|\,\text{Pa}[{\boldsymbol{w}(\boldsymbol{u})}])=\prod_{s=1}^{q}p(w(\boldsymbol{\ell},\xi_{s})\,|\,w(\boldsymbol{\ell},\xi_{1}),\dots,w(\boldsymbol{\ell},\xi_{s-1}),\text{Pa}[{\boldsymbol{w}(\boldsymbol{\ell})}]))$,
which also implies that the size of the parent set may depend on the variable
index. Assumptions of conditional independence across variables can be encoded
similarly. Also note that any specific choice of $\eta_{B}$ induces a
partition on ${\cal U}$; let $U_{j}=\\{\boldsymbol{u}\in{\cal
U}:\eta_{B}(\boldsymbol{u})=\boldsymbol{b}_{j}\\}$, then clearly ${\cal
U}=\cup_{j=1}^{m_{U}}U_{j}$ with $U_{i}\cap U_{j}=\emptyset$ if $i\neq j$.
This partition does not necessarily correspond to the partitioning scheme used
on ${\cal S}$. $\eta_{B}$ may by designed to ignore part of the tree and
result in $m_{U}<m_{B}$. However, we can just drop the unused leaves from
${\cal G}$ and set $\text{Ch}[{\boldsymbol{a}}]=\emptyset$ for terminal nodes
whose leaf is inactive, resulting in $m_{U}=m_{B}$. We will thus henceforth
assume that $m_{U}=m_{B}$ without loss of generality.

###  2.2 SpamTrees as a standalone spatial process

We define a valid joint density for any finite set of locations in ${\cal
D}^{*}$ satisfying the Kolmogorov consistency conditions in order to define a
valid process. We approach this problem analogously to Datta et al., 2016a and
Peruzzi et al.,  (2020). Enumerate each of the $m_{S}$ reference subsets as
$S_{i}=\\{\boldsymbol{s}_{i_{1}},\dots,\boldsymbol{s}_{i_{n_{i}}}\\}$ where
$\\{i_{1},\ldots,i_{n_{i}}\\}\subset\\{1,\ldots,n_{{\cal S}}\\}$, and each of
the $m_{U}$ non-reference subsets as
$U_{i}=\\{\boldsymbol{u}_{i_{1}},\dots,\boldsymbol{u}_{i_{n_{i}}}\\}$ where
$\\{i_{1},\ldots,i_{n_{i}}\\}\subset\\{1,\ldots,n_{{\cal U}}\\}$. Then
introduce ${\cal V}=\\{V_{1},\dots,V_{m_{V}}\\}$ where $m_{V}=m_{S}+m_{U}$ and
$V_{i}=S_{i}$ for $i=1,\dots,m_{S}$, $V_{m_{S}+i}=U_{i}$ for
$i=1,\dots,m_{U}$. Then take
$\boldsymbol{w}_{i}=(w(\boldsymbol{\ell}_{i_{1}}),\ldots,w(\boldsymbol{\ell}_{i_{n_{i}}}))^{\top}$
as the $n_{i}\times 1$ random vector with elements of $w(\boldsymbol{\ell})$
for each $\boldsymbol{\ell}\in V_{i}$. Denote
$\boldsymbol{w}_{[{i}]}=\boldsymbol{w}(\eta^{-1}(\text{Pa}[{\boldsymbol{v}_{i}}]))$.
Then

$\displaystyle\widetilde{p}(\boldsymbol{w}_{{\cal
S}})=\widetilde{p}(\boldsymbol{w}_{1},\dots,\boldsymbol{w}_{m_{S}})=\prod_{r=0}^{M-1}\prod_{i:\\{\boldsymbol{v}_{i}\in\boldsymbol{A}_{r}\\}}p(\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[{i}]})$
$\displaystyle\quad\widetilde{p}(\boldsymbol{w}_{{\cal
U}}\mid\boldsymbol{w}_{{\cal
S}})=\prod_{i:\\{\boldsymbol{v}_{i}\in\boldsymbol{B}\\}}p(\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[{i}]})$
(6) $\displaystyle\widetilde{p}(\boldsymbol{w}_{{\cal
S}})\widetilde{p}(\boldsymbol{w}_{{\cal U}}\mid\boldsymbol{w}_{{\cal
S}})=\prod_{r=0}^{M-1}\prod_{i:\\{\boldsymbol{v}_{i}\in\boldsymbol{A}_{r}\\}}p(\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[{i}]})$
$\displaystyle\prod_{i:\\{\boldsymbol{v}_{i}\in\boldsymbol{B}\\}}p(\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[{i}]})$

which is a proper multivariate joint density since ${\cal G}$ is acyclic
(Lauritzen, , 1996). All locations inside $U_{j}$ always share the same parent
set, but a parent set is not necessarily unique to a single $U_{j}$. This
includes as a special case a scenario in which one can assume

$\displaystyle\widetilde{p}(\boldsymbol{w}_{{\cal U}}\mid\boldsymbol{w}_{{\cal
S}})=\prod_{j=1}^{m_{U}}\prod_{i=1}^{|U_{j}|}p(w(\boldsymbol{u}_{i})\,|\,\boldsymbol{w}(\eta^{-1}(\text{Pa}[{\boldsymbol{b}_{j}}])));$
(7)

in this case each location corresponds to a leaf node. To conclude the
construction, for any finite subset of spatial locations
$\mathcal{L}\subset\mathcal{D}$ we can let ${\cal U}=\mathcal{L}\setminus{\cal
S}$ and obtain

$\widetilde{p}(\boldsymbol{w}_{\mathcal{L}})=\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})),$

leading to a well-defined process satisfying the Kolmogorov conditions (see
Appendix A).

####  2.2.1 Positioning of spatial locations in conditioning sets

In spatial models based on sparse DAGs, larger conditioning sets yield
processes that are closer to the base process $p$ in terms of Kullback-Leibler
divergence (Banerjee, , 2020; Peruzzi et al., , 2020), denoted as
$KL(\cdot\|p)$. The same results cannot be applied directly to SpamTrees given
the treed structure of the DAG. For a given ${\cal S}$, we consider the
distinct but related issues of placing individual locations into reference
subsets (1) at different levels of the treed hierarchy; (2) within the same
level of the hierarchy.

######  Proposition 2.1.

Suppose ${\cal S}={\cal S}_{0}\cup{\cal S}_{1}$ where $S_{0}\cap
S_{1}=\emptyset$ and ${\cal S}_{1}=S_{11}\cup S_{12}$, $S_{11}\cap
S_{12}=\emptyset$. Take $\boldsymbol{s}^{*}\notin{\cal S}$. Consider the graph
${\cal
G}=\\{\boldsymbol{V}=\\{\boldsymbol{v}_{0},\boldsymbol{v}_{1},\boldsymbol{v}_{2}\\},\boldsymbol{E}=\\{\boldsymbol{v}_{0}\to\boldsymbol{v}_{1},\boldsymbol{v}_{0}\to\boldsymbol{v}_{2}\\}\\}$;
denote as $p_{0}$ the density of a SpamTree using $\eta({\cal
S}_{0}\cup\\{\boldsymbol{s}^{*}\\})=\boldsymbol{v}_{0}$,
$\eta(S_{11})=\boldsymbol{v}_{1}$ and $\eta(S_{12})=\boldsymbol{v}_{2}$,
whereas let $p_{1}$ be the density of a SpamTree with $\eta({\cal
S}_{0})=\boldsymbol{v}_{0}$,
$\eta(S_{11}\cup\\{\boldsymbol{s}^{*}\\})=\boldsymbol{v}_{1}$ and
$\eta(S_{12})=\boldsymbol{v}_{2}$. Then $KL(p_{1}\|p)-KL(p_{0}\|p)>0.$

The proof proceeds by an “information never hurts” argument (Cover and Thomas,
, 1991). Denote ${\cal S}^{*}={\cal S}\cup\\{\boldsymbol{s}^{*}\\}$,
$\boldsymbol{w}^{*}=\boldsymbol{w}_{{\cal S}^{*}}$,
$w^{*}=w(\boldsymbol{s}^{*})$ and
$\boldsymbol{w}_{j}^{*}=(\boldsymbol{w}_{j}^{\top},w^{*})^{\top}$. Then

$\displaystyle p_{0}(\boldsymbol{w}^{*})$
$\displaystyle=p(\boldsymbol{w}_{0}^{*})p(\boldsymbol{w}_{1}\,|\,\boldsymbol{w}_{0}^{*})p(\boldsymbol{w}_{2}\,|\,\boldsymbol{w}_{0}^{*})=p(\boldsymbol{w}_{0})p(w^{*}\,|\,\boldsymbol{w}_{0})p(\boldsymbol{w}_{1}\,|\,\boldsymbol{w}_{0},w^{*})p(\boldsymbol{w}_{2}\,|\,\boldsymbol{w}_{0}^{*})$
$\displaystyle p_{1}(\boldsymbol{w}^{*})$
$\displaystyle=p(\boldsymbol{w}_{0})p(\boldsymbol{w}_{1}^{*}\,|\,\boldsymbol{w}_{0})p(\boldsymbol{w}_{2}\,|\,\boldsymbol{w}_{0})=p(\boldsymbol{w}_{0})p(w^{*}\,|\,\boldsymbol{w}_{0})p(\boldsymbol{w}_{1}\,|\,\boldsymbol{w}_{0},w^{*})p(\boldsymbol{w}_{2}\,|\,\boldsymbol{w}_{0}),$

therefore
$p_{0}(\boldsymbol{w}^{*})/p_{1}(\boldsymbol{w}^{*})=p(\boldsymbol{w}_{2}\,|\,\boldsymbol{w}_{0}^{*})/p(\boldsymbol{w}_{2}\,|\,\boldsymbol{w}_{0})$;
then by Jensen’s inequality

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

Intuitively, this shows that there is a penalty associated to positioning
reference locations at higher levels of the treed hierarchy. Increasing the
size of the reference set at the root augments the conditioning sets at all
its children; since this is not true when the increase is at a branch level,
the density $p_{0}$ is closer to $p$ than $p_{1}$. In other words there is a
cost of branching in ${\cal G}$ which must be justified by arguments related
to computational efficiency. The above proposition also suggests populating
near-root branches with locations of sparsely-observed outcomes. Not doing so
in highly imbalanced settings may result in possibly too restrictive spatial
conditional independence assumptions.

######  Proposition 2.2.

Consider the same setup as Proposition 2.1 and let $p_{2}$ be the density of a
SpamTree such that
$\eta(S_{12}\cup\\{\boldsymbol{s}^{*}\\})=\boldsymbol{v}_{2}$. Let $H_{p}$ be
the conditional entropy of base process $p$. Then
$H_{p}(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{2})<H_{p}(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{1})$
implies $KL(p_{2}\|p)<KL(p_{1}\|p)$.

The density of the new model is

$p_{2}(\boldsymbol{w}^{*})=p(\boldsymbol{w}_{0})p(\boldsymbol{w}_{1}\,|\,\boldsymbol{w}_{0})p(\boldsymbol{w}_{2}^{*}\,|\,\boldsymbol{w}_{0})=p(\boldsymbol{w}_{0})p(\boldsymbol{w}_{1}\,|\,\boldsymbol{w}_{0})p(\boldsymbol{w}_{2}\,|\,\boldsymbol{w}_{0})p(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{2}).$

Then, noting that
$p(\boldsymbol{w}^{*}_{1}\,|\,\boldsymbol{w}_{0})=p(\boldsymbol{w}_{1}\,|\,\boldsymbol{w}_{0})p(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{1})$,
we get
$\frac{p_{1}(\boldsymbol{w}^{*})}{p_{2}(\boldsymbol{w}^{*})}=\frac{p(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{1})}{p(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{2})}$
and

$\displaystyle KL(p_{2}\|p)-KL(p_{1}$ $\displaystyle\|p)=\int\log
p(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{1})p(\boldsymbol{w}^{*})d\boldsymbol{w}^{*}-\int\log
p(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{2})p(\boldsymbol{w}^{*})d\boldsymbol{w}^{*}$
$\displaystyle=H_{p}(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{2})-H_{p}(w^{*}\,|\,\boldsymbol{w}_{0},\boldsymbol{w}_{1}).$

This result suggests placing a new reference location $\boldsymbol{s}^{*}$ in
the reference subset least uncertain about the realization of the process at
$\boldsymbol{s}^{*}$. We interpret this as justifying recursive domain
partitioning on ${\cal S}$ in spatial contexts in which local spatial clusters
of locations are likely less uncertain about process realization in the same
spatial region. In the remainder of this article, we will consider a given
reference set ${\cal S}$ which typically will be based on a subset of observed
locations; the combinatorial problem of selecting an optimal ${\cal S}$ (in
some sense) is beyond the scope of this article. If ${\cal S}$ is not
partitioned, it can be considered as a set of knots or “sensors” and one can
refer to a large literature on experimental design and optimal sensor
placement (see e.g. Krause et al., , 2008, and references therein). It might
be possible to extend previous work on adaptive knot placement (Guhaniyogi et
al., , 2011), but this will come at a steep cost in terms of computational
performance.

##  3 Bayesian spatial regressions using SpamTrees

Suppose we observe an $l$-variate outcome at spatial locations
$\boldsymbol{\ell}\in{\cal D}\subset\Re^{d}$ which we wish to model using a
spatially-varying regression model:

$\displaystyle y_{j}(\boldsymbol{\ell})$
$\displaystyle=\boldsymbol{x}_{j}(\boldsymbol{\ell})^{\top}\bbeta_{j}+\sum_{k}z_{jk}(\boldsymbol{\ell})w(\boldsymbol{\ell},\boldsymbol{\xi}_{k})+\varepsilon_{j}(\boldsymbol{\ell}),\quad
j=1,\dots,l,$ (9)

where $y_{j}(\boldsymbol{\ell})$ is the $j$-th point-referenced outcome at
$\boldsymbol{\ell}$, $\boldsymbol{x}_{j}(\boldsymbol{\ell})$ is a $p_{j}\times
1$ vector of spatially referenced predictors linked to constant coefficients
$\bbeta_{j}$,
$\varepsilon_{j}(\boldsymbol{\ell})\overset{iid}{\sim}N(0,\tau^{2}_{j})$ is
the measurement error for outcome $j$, and $z_{jk}(\boldsymbol{\ell})$ is the
$k$-th (of $q$) covariates for the $j$-th outcome modeled with spatially-
varying coefficient $w(\boldsymbol{\ell},\boldsymbol{\xi}_{k})$,
$\boldsymbol{\ell}\in{\cal D}$, $\boldsymbol{\xi}_{k}\in\Xi$. This coefficient
$w(\boldsymbol{\ell},\boldsymbol{\xi}_{k})$ corresponds to the $k$-th margin
of a $q$-variate Gaussian process
$\\{\boldsymbol{w}(\boldsymbol{\ell}):\boldsymbol{\ell}\in{\cal D}\\}$ denoted
as $\boldsymbol{w}(\boldsymbol{\ell})\sim
GP(\mathbf{0},\boldsymbol{C}_{\boldsymbol{\theta}}(\cdot,\cdot))$ with cross-
covariance $\boldsymbol{C}_{\boldsymbol{\theta}}$ indexed by unknown
parameters $\boldsymbol{\theta}$ which we omit in notation for simplicity. A
valid cross-covariance function is defined as
$\boldsymbol{C}_{\boldsymbol{\theta}}:{\cal D}\times{\cal D}\rightarrow{\cal
M}_{q\times q}$, where ${\cal M}_{q\times q}$ is a subset of the space of all
$q\times q$ real matrices $\Re^{q\times q}$. It must satisfy
$\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=\boldsymbol{C}(\boldsymbol{\ell}^{\prime},\boldsymbol{\ell})^{\top}$
for any two locations $\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}\in{\cal
D}$, and
$\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 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}\\}$.

We replace the full GP with a Gaussian SpamTree for scalable computation
considering the $q$-variate multivariate Gaussian process
$\boldsymbol{w}(\cdot)$ as the base process. Since the $(i,j)$-th entry of
$\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})$ is
$\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})_{i,j}=\mbox{Cov}(w_{i}(\boldsymbol{\ell}),w_{j}(\boldsymbol{\ell}^{\prime}))$,
i.e. the covariance between the $i$-th and $j$-th elements of
$\boldsymbol{w}(\boldsymbol{\ell})$ at $\boldsymbol{\ell}$ and
$\boldsymbol{\ell}^{\prime}$, we can obtain a covariance function on the
augmented domain $\boldsymbol{C}^{*}:{\cal D}^{*}\times{\cal
D}^{*}\rightarrow\Re$ as
$\boldsymbol{C}^{*}((\boldsymbol{\ell},\boldsymbol{\xi}),(\boldsymbol{\ell}^{\prime},\boldsymbol{\xi}^{\prime}))=\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})_{i,i^{\prime}}$
where $\boldsymbol{\xi}$ and $\boldsymbol{\xi}^{\prime}$ are the locations in
$\Xi$ of variables $i$ and $j$, respectively. Apanasovich and Genton,  (2010)
use a similar representation to build valid cross-covariances based on
existing univariate covariance functions; their approach amounts to
considering $\boldsymbol{\xi}$ or
$\|\boldsymbol{\xi}-\boldsymbol{\xi}^{\prime}\|$ as a parameter to be
estimated. Our approach can be based on any valid cross-covariance as we may
just set $\Xi=\\{1,\dots,q\\}$. Refer to e.g. Genton and Kleiber,  (2015) for
an extensive review of cross-covariance functions for multivariate processes.
Moving forward, we will not distinguish between $\boldsymbol{C}^{*}$ and
$\boldsymbol{C}$. The linear multivariate spatially-varying regression model
(9) allows the $l$ outcomes to be observed at different locations; we later
consider the case $l=q$ and $\boldsymbol{Z}(\boldsymbol{\ell})=I_{q}$
resulting in a multivariate space-varying intercept model.

###  3.1 Gaussian SpamTrees

Enumerate the set of nodes as
$\boldsymbol{V}=\\{\boldsymbol{v}_{1},\dots,\boldsymbol{v}_{m_{V}}\\}$,
$m_{V}=m_{S}+m_{U}$ and denote
$\boldsymbol{w}_{i}=w(\eta^{-1}(\boldsymbol{v}_{i}))$, $\boldsymbol{C}_{ij}$
as the $n_{i}\times n_{j}$ covariance matrix between $\boldsymbol{w}_{i}$ and
$\boldsymbol{w}_{j}$, $\boldsymbol{C}_{i,[i]}$ the $n_{i}\times J_{i}$
covariance matrix between $\boldsymbol{w}_{i}$ and $\boldsymbol{w}_{[{i}]}$,
$\boldsymbol{C}_{i}$ the $n_{i}\times n_{i}$ covariance matrix between
$\boldsymbol{w}_{i}$ and itself, and $\boldsymbol{C}_{[i]}$ the $J_{i}\times
J_{i}$ covariance matrix between $\boldsymbol{w}_{[{i}]}$ and itself. A base
Gaussian process induces $\widetilde{p}(\boldsymbol{w}_{{\cal
S}})=\prod_{j:\\{\boldsymbol{v}_{j}\in\boldsymbol{A}\\}}N(\boldsymbol{w}_{j}\mid\boldsymbol{H}_{j}\boldsymbol{w}_{[{j}]},\boldsymbol{R}_{j})$,
where

$\boldsymbol{H}_{j}=\boldsymbol{C}_{j,[j]}\boldsymbol{C}^{-1}_{[j]}\quad\text{and}\quad\boldsymbol{R}_{j}=\boldsymbol{C}_{j}-\boldsymbol{C}_{j,[j]}\boldsymbol{C}^{-1}_{[j]}\boldsymbol{C}_{[j],j},$
(10)

implying that the joint density $\widetilde{p}(\boldsymbol{w}_{{\cal S}})$ is
multivariate normal with covariance $\widetilde{\boldsymbol{C}}_{{\cal S}}$
and precision matrix $\widetilde{\boldsymbol{C}}^{-1}_{{\cal S}}$. At ${\cal
U}$ we have $\widetilde{p}(\boldsymbol{w}_{{\cal U}}\mid\boldsymbol{w}_{{\cal
S}})=\prod_{j:\\{\boldsymbol{v}_{j}\in\boldsymbol{B}\\}}N(\boldsymbol{w}_{j}\mid\boldsymbol{H}_{j}\boldsymbol{w}_{[{j}]},\boldsymbol{R}_{j})$,
where $\boldsymbol{H}_{j}$ and $\boldsymbol{R}_{j}$ are as in (10). All
quantities can be computed using the base cross-covariance function. Given
that the $\widetilde{p}$ densities are Gaussian, so will be the finite
dimensional distributions.

The treed graph ${\cal G}$ leads to properties which we analyze in more detail
in Appendix B and summarize here. For two nodes
$\boldsymbol{v}_{i},\boldsymbol{v}_{j}\in\boldsymbol{V}$ denote the common
descendants as
$\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=(\\{\boldsymbol{v}_{i}\\}\cup\text{Ch}[{\boldsymbol{v}_{i}}])\cap(\\{\boldsymbol{v}_{j}\\}\cup\text{Ch}[{\boldsymbol{v}_{j}}])$.
If $\boldsymbol{v}_{i}\in\text{Pa}[{\boldsymbol{v}_{j}}]$ denote
$\boldsymbol{H}_{i\to j}$ and $\boldsymbol{H}_{\setminus i\to j}$ as the
matrix obtained by subsetting $\boldsymbol{H}_{j}$ to columns corresponding to
$\boldsymbol{v}_{i}$, or to
$\text{Pa}[{\boldsymbol{v}_{j}}]\setminus\\{\boldsymbol{v}_{i}\\}$,
respectively. Similarly define $\boldsymbol{w}_{[i\to j]}=\boldsymbol{w}_{i}$
and $\boldsymbol{w}_{[\setminus i\to j]}$. As a special case, if the tree
depth is $\delta=1$ and
$\\{\boldsymbol{v}_{j}\\}=\text{Pa}[{\boldsymbol{v}_{i}}]$ then
$\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=\\{\boldsymbol{v}_{i}\\}$,
$\boldsymbol{H}_{i\to j}=\boldsymbol{H}_{j}$, and $\boldsymbol{w}_{[i\to
j]}=\boldsymbol{w}_{[{j}]}$. Define $\boldsymbol{\mathcal{H}}$ as the matrix
whose $(i,j)$ block is
$\boldsymbol{\mathcal{H}}_{ij}=\boldsymbol{O}_{n_{i}\times n_{j}}$ if
$\boldsymbol{v}_{j}\notin\text{Pa}[{\boldsymbol{v}_{i}}]$, and otherwise
$\boldsymbol{\mathcal{H}}_{ij}=\boldsymbol{H}_{j\to i}$.

####  3.1.1 Precision matrix

The $(i,j)$ block of the precision matrix at both reference and non-reference
locations $\widetilde{\boldsymbol{C}}^{-1}$ is denoted by
$\widetilde{\boldsymbol{C}}^{-1}(i,j)$, with $i,j=1,\dots,m_{V}$ corresponding
to nodes $\boldsymbol{v}_{i},\boldsymbol{v}_{j}\in\boldsymbol{V}$ for some
$i,j$; it is nonzero if
$\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=\emptyset$, otherwise:

$\displaystyle\widetilde{\boldsymbol{C}}^{-1}(i,j)$
$\displaystyle=\sum\limits_{\boldsymbol{v}_{k}\in\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})}(\boldsymbol{I}_{ki}-\boldsymbol{H}_{i\to
k})^{\top}\boldsymbol{R}_{k}^{-1}(\boldsymbol{I}_{kj}-\boldsymbol{H}_{j\to
k})$ (11)
$\displaystyle=\sum\limits_{\boldsymbol{v}_{k}\in\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})}(\boldsymbol{I}_{ki}-\boldsymbol{\mathcal{H}}_{ki})^{\top}\boldsymbol{R}_{k}^{-1}(\boldsymbol{I}_{kj}-\boldsymbol{\mathcal{H}}_{kj}),$

where $\boldsymbol{I}_{ij}$ is the $(i,j)$ block of an identity matrix with
$n_{{\cal S}}+n_{{\cal U}}$ rows and is nonzero if and only if $i=j$. We thus
obtain that the number of nonzero elements of
$\widetilde{\boldsymbol{C}}^{-1}$ is

$\displaystyle\text{nnz}(\widetilde{\boldsymbol{C}}^{-1})$
$\displaystyle=\sum_{i=1}^{m_{V}}\left(2n_{i}J_{i}+n_{i}^{2}\boldsymbol{1}\\{\boldsymbol{v}_{i}\in\boldsymbol{V}\\}\right),$
(12)

where $n_{i}=|\eta^{-1}(\boldsymbol{v}_{i})|$,
$J_{i}=|\eta^{-1}(\text{Pa}[{\boldsymbol{v}_{i}}])|$, and by symmetry
$(\widetilde{\boldsymbol{C}}^{-1}(i,j))^{\top}=\widetilde{\boldsymbol{C}}^{-1}(j,i)$.

If $\delta>1$, the size of $\boldsymbol{C}_{[i]}$ is larger for nodes
$\boldsymbol{v}_{i}$ at levels of the treed hierarchy farther from
$\boldsymbol{A}_{M_{\delta}}$. However suppose
$\boldsymbol{v}_{i},\boldsymbol{v}_{j}$ are such that
$\text{Pa}[{\boldsymbol{v}_{j}}]=\\{\boldsymbol{v}_{i}\\}\cup\text{Pa}[{\boldsymbol{v}_{i}}]$.
Then computing $\boldsymbol{C}_{[j]}^{-1}$ proceeds more cheaply by
recursively applying the following:

$\displaystyle\boldsymbol{C}_{[j]}^{-1}=\begin{bmatrix}\boldsymbol{C}_{[i]}^{-1}+\boldsymbol{H}_{i}^{\top}\boldsymbol{R}_{i}^{-1}\boldsymbol{H}_{i}&-\boldsymbol{H}_{i}^{\top}\boldsymbol{R}_{i}^{-1}\\\
-\boldsymbol{R}_{i}^{-1}\boldsymbol{H}_{i}&\boldsymbol{R}_{i}^{-1}\end{bmatrix}.$
(13)

####  3.1.2 Induced covariance

Define a path from $\boldsymbol{v}_{k}$ to $\boldsymbol{v}_{j}$ as ${\cal
P}_{k\to j}=\\{\boldsymbol{v}_{i_{1}},\dots,\boldsymbol{v}_{i_{r}}\\}$ where
$\boldsymbol{v}_{i_{1}}=\boldsymbol{v}_{k}$,
$\boldsymbol{v}_{i_{r}}=\boldsymbol{v}_{j}$, and
$\boldsymbol{v}_{i_{h}}\in\text{Pa}[{\boldsymbol{v}_{i_{h+1}}}]$. The longest
path $\widetilde{{\cal P}}_{k\to j}$ is such that if
$\boldsymbol{v}_{k}\in\boldsymbol{A}_{r_{k}}$ and
$\boldsymbol{v}_{j}\in\boldsymbol{A}_{r_{j}}$ then $|\widetilde{{\cal
P}}_{k\to j}|=r_{j}-r_{k}+1$. The shortest path $\bar{{\cal P}}_{k\to j}$ is
the path from $\boldsymbol{v}_{k}$ to $\boldsymbol{v}_{j}$ with minimum number
of steps. We denote the longest path from the root to $\boldsymbol{v}_{j}$ as
$\widetilde{{\cal P}}_{0\to j}$; this corresponds to the full set of ancestors
of $\boldsymbol{v}_{j}$, and
$\text{Pa}[{\boldsymbol{v}_{j}}]\subset\widetilde{{\cal P}}_{0\to j}$. For two
nodes $\boldsymbol{v}_{i}$ and $\boldsymbol{v}_{j}$ we have
$(\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}])\subset(\widetilde{{\cal
P}}_{0\to i}\cap\widetilde{{\cal P}}_{0\to j})$. We define the concestor
between $\boldsymbol{v}_{i}$ and $\boldsymbol{v}_{j}$ as
$\text{con}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=\arg\max_{\boldsymbol{v}_{k}\in\boldsymbol{V}}\\{k:{\cal
P}_{k\to i}\cap{\cal P}_{k\to j}\neq\emptyset\\}$ i.e. the last common
ancestor of the two nodes.

Take the path $\widetilde{{\cal P}}_{M_{\delta}\to j}$ in ${\cal G}$ from a
node at $\boldsymbol{A}_{M_{\delta}}$ leading to $\boldsymbol{v}_{j}$. After
defining the cross-covariance function
$\boldsymbol{K}_{i}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=\boldsymbol{C}_{\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}}-\boldsymbol{C}_{\boldsymbol{\ell},[i]}\boldsymbol{C}^{-1}_{[i]}\boldsymbol{C}_{[i],\boldsymbol{\ell}^{\prime}}$
and denoting
$\boldsymbol{K}_{i}(\boldsymbol{\ell},s)=\boldsymbol{K}_{i}(\boldsymbol{\ell},\eta^{-1}(\boldsymbol{v}_{s}))$
we can write

$\displaystyle\boldsymbol{w}_{j}=\sum\limits_{s=i_{M_{\delta}}}^{i_{r-1}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{e}_{j},$
(14)

where for $s>i_{M_{\delta}}$ the $\boldsymbol{e}_{s}$ are independent zero-
mean GPs with covariance $K_{s}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})$
and we set
$K_{i_{M_{\delta}}}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})$
and $\boldsymbol{e}_{i_{M_{\delta}}}=\boldsymbol{w}_{i_{M_{\delta}}}\sim
N(0,\boldsymbol{C}_{i_{M_{\delta}}})$. Take two locations
$\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}$ such that
$\boldsymbol{v}_{i}=\eta(\boldsymbol{\ell}),\boldsymbol{v}_{j}=\eta(\boldsymbol{\ell}^{\prime})$
and let
$\boldsymbol{v}_{z}=\text{con}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})$; if
$\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]\neq\emptyset$
then the above leads to

$\displaystyle
Cov_{\widetilde{p}}(\boldsymbol{w}(\boldsymbol{\ell}),\boldsymbol{w}(\boldsymbol{\ell}^{\prime}))$
$\displaystyle=\sum\limits_{s\in\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]}\boldsymbol{K}_{s}(\boldsymbol{\ell},s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{K}_{s}(s,\boldsymbol{\ell}^{\prime})+\boldsymbol{1}\\{\boldsymbol{\ell}=\boldsymbol{\ell}^{\prime}\\}\boldsymbol{K}_{j}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}),$
(15)

where
$\boldsymbol{K}_{z}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=\boldsymbol{C}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})$.
If
$\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]=\emptyset$
take the shortest paths $\bar{{\cal P}}_{z\to i}=\\{i_{1},\dots,i_{r_{i}}\\}$
and $\bar{{\cal P}}_{z\to j}=\\{j_{1},\dots,j_{r_{j}}\\}$; setting
$\boldsymbol{F}_{i_{h}}=\boldsymbol{C}_{i_{h},i_{h-1}}\boldsymbol{C}_{i_{h-1}}^{-1}$
we get

$\displaystyle
Cov_{\widetilde{p}}(\boldsymbol{w}(\boldsymbol{\ell}),\boldsymbol{w}(\boldsymbol{\ell}^{\prime}))$
$\displaystyle=\boldsymbol{F}_{i_{r_{i}}}\cdots\boldsymbol{F}_{i_{1}}\boldsymbol{C}_{z}\boldsymbol{F}_{j_{1}}^{\top}\cdots\boldsymbol{F}_{j_{r_{j}}}^{\top}.$
(16)

In particular if $\delta=M$ then
$\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]\neq\emptyset$
for all $i,j$ and only (15) is used, whereas if $\delta=1$ then the only
scenario in which (15) holds is
$\\{\boldsymbol{v}_{z}\\}=\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]$
in which case the two are equivalent. In univariate settings, the special case
in which $\delta=M$, and hence $M_{\delta}=0$, leads to an interpretation of
(14) as a basis function decomposition; considering all leaf paths ${\cal
P}_{j}$ for $\boldsymbol{v}_{j}\in\boldsymbol{B}$, this leads to an MRA
(Katzfuss, , 2017; Katzfuss and Gong, , 2019). On the other hand, keeping
other parameters constant, $\delta<M$ and in particular $\delta=1$ may be
associated to savings in computing cost, leading to a trade-off between graph
complexity and size of reference subsets; see Appendix B.5.

####  3.1.3 Block-sparse Cholesky decompositions

In recent work Jurek and Katzfuss,  (2020) consider sparse Cholesky
decompositions of covariance and precision matrices for treed graphs
corresponding to the case $\delta=M$ above in the context of space-time
filtering; their methods involve sparse Cholesky routines on reverse orderings
of $\widetilde{\boldsymbol{C}}^{-1}$ at the level of individual locations. In
doing so, the relationship between Cholesky decompositions and ${\cal G}$,
$\widetilde{\boldsymbol{C}}^{-1}$ and the block structure in ${\cal S}$
remains somewhat hidden, and sparse Cholesky libraries are typically
associated to bottlenecks in MCMC algorithms. However we note that a
consequence of (11) is that it leads to a direct algorithm, for any $\delta$,
for the block-decomposition of any symmetric positive-definite matrix
$\boldsymbol{\Lambda}$ conforming to ${\cal G}$, i.e. with the same block-
sparse structure as $\widetilde{\boldsymbol{C}}^{-1}$. This allows us to write
$\boldsymbol{\Lambda}=(\boldsymbol{I}-\boldsymbol{L})^{\top}\boldsymbol{D}(\boldsymbol{I}-\boldsymbol{L})$
where $\boldsymbol{I}$ is the identity matrix, $\boldsymbol{L}$ is block lower
triangular with the same block-sparsity pattern as $\boldsymbol{\mathcal{H}}$
above, and $\boldsymbol{D}$ is block diagonal symmetric positive-definite. In
Appendix B.2.3 we outline Algorithm 4 which (i) makes direct use of the
structure of ${\cal G}$, (ii) computes the decomposition at blocks of
reference and non-reference locations, and (iii) requires no external sparse
matrix library, in particular no sparse Cholesky solvers. Along with Algorithm
5 for the block-computation of $(\boldsymbol{I}-\boldsymbol{L})^{-1}$, it can
be used to compute
$\boldsymbol{\Lambda}^{-1}=(\widetilde{\boldsymbol{C}}^{-1}+\boldsymbol{\Sigma})^{-1}$
where $\boldsymbol{\Sigma}$ is a block-diagonal matrix; it is thus useful in
computing the Gaussian integrated likelihood.

###  3.2 Estimation and prediction

We introduce notation to aid in obtaining the full conditional distributions.
Write (9) as

$\displaystyle\boldsymbol{y}(\boldsymbol{\ell})$
$\displaystyle=\boldsymbol{X}(\boldsymbol{\ell})\bbeta+\boldsymbol{Z}(\boldsymbol{\ell})\boldsymbol{w}(\boldsymbol{\ell})+\boldsymbol{\varepsilon}(\boldsymbol{\ell}),$
(17)

where
$\boldsymbol{y}(\boldsymbol{\ell})=(\\{y_{j}(\boldsymbol{\ell})\\}_{j=1}^{l})^{\top}$,
$\boldsymbol{\varepsilon}(\boldsymbol{\ell})=(\\{\varepsilon_{j}(\boldsymbol{\ell})\\}_{j=1}^{l})^{\top}\sim
N(\mathbf{0},\boldsymbol{D}_{\tau})$,
$\boldsymbol{D}_{\tau}=\text{diag}(\tau^{2}_{1},\dots,\tau^{2}_{l})$,
$\boldsymbol{X}(\boldsymbol{\ell})=\text{blockdiag}\\{\boldsymbol{x}_{j}(\boldsymbol{\ell})^{\top},j=1,\dots,l\\}$,
$\bbeta=(\bbeta_{p_{1}}^{\top},\dots,\bbeta_{p_{j}}^{\top})^{\top}$. The
$l\times q$ matrix
$\boldsymbol{Z}(\boldsymbol{\ell})=(\boldsymbol{z}_{j}(\boldsymbol{\ell})^{\top},j=1,\dots,l)$
with
$\boldsymbol{z}_{j}(\boldsymbol{\ell})^{\top}=(z_{jk}(\boldsymbol{\ell}),k=1,\dots,q)$
acts a design matrix for spatial location $\boldsymbol{\ell}$. Collecting all
locations along the $j$-th margin, we build ${\cal
T}_{j}=\\{\boldsymbol{\ell}_{1}^{(j)},\dots,\boldsymbol{\ell}_{N_{j}}^{(j)}\\}$
and ${\cal T}=\cup_{j}{\cal T}_{j}$. We then call
$\boldsymbol{y}^{(j)}=(y_{j}(\boldsymbol{\ell}_{1}^{(j)}),\dots,y_{j}(\boldsymbol{\ell}_{N_{j}}^{(j)}))^{\top}$
and $\boldsymbol{\varepsilon}^{(j)}$ similarly,
$\boldsymbol{X}^{(j)}=(\boldsymbol{x}_{j}(\boldsymbol{\ell}_{1}^{(j)}),\dots,\boldsymbol{x}_{j}(\boldsymbol{\ell}_{N_{j}}^{(j)}))^{\top}$,
$\boldsymbol{w}^{(j)}=(\boldsymbol{w}(\boldsymbol{\ell}_{1}^{(j)},\boldsymbol{\xi})^{\top},\dots,\boldsymbol{w}(\boldsymbol{\ell}_{N_{j}}^{(j)},\boldsymbol{\xi})^{\top})^{\top}$
and
$\boldsymbol{Z}^{(j)}=\text{blockdiag}\\{\boldsymbol{z}_{j}(\boldsymbol{\ell}_{s}^{(j)})^{\top}\\}_{s=1}^{N_{j}}$.
The full observed data are $\boldsymbol{y},\boldsymbol{X},\boldsymbol{Z}$.
Denoting the number of observations as $n=\sum_{j=1}^{l}N_{j}$,
$\boldsymbol{Z}$ is thus a $n\times qn$ block-diagonal matrix, and similarly
$\boldsymbol{w}$ is a $qn\times 1$ vector. We introduce the diagonal matrix
$\boldsymbol{D}_{n}$ such that
$\text{diag}(\boldsymbol{D}_{n})=(\tau^{2}_{1}\boldsymbol{1}_{N_{1}}^{\top},\dots,\tau^{2}_{l}\boldsymbol{1}_{N_{l}}^{\top})^{\top}$.

By construction we may have $\eta(S_{i})=\boldsymbol{v}_{i}$ and
$\eta(S_{j})=\boldsymbol{v}_{j}$ such that
$(\boldsymbol{\ell},\boldsymbol{\xi})\in S_{i}$ and
$(\boldsymbol{\ell}^{\prime},\boldsymbol{\xi}^{\prime})\in S_{j}$ where
$\boldsymbol{\ell}^{\prime}=\boldsymbol{\ell}$,
$\boldsymbol{\xi}\neq\boldsymbol{\xi}^{\prime}$ and similarly for non-
reference subsets. Suppose ${\cal A}\subset{\cal D}\times\Xi$ is a generic
reference or non-reference subset. We denote $\bar{{\cal A}}\subset{\cal
D}\times\Xi$ as the set of all combinations of spatial locations of ${\cal A}$
and variables i.e. $\bar{{\cal A}}={\cal A}\big{|}_{{\cal D}}\times{\cal
A}\big{|}_{\Xi}$ where ${\cal A}\big{|}_{{\cal D}}\subset{\cal D}$ is the set
of unique spatial locations in ${\cal A}$ and ${\cal A}\big{|}_{\Xi}$ are the
unique latent variable coordinates. By subtraction we find ${\cal
A}_{-}=\bar{{\cal A}}\setminus{\cal A}$ as the set of locations whose spatial
location is in ${\cal A}$ but whose variable is not. Let
$\boldsymbol{y}(\bar{{\cal A}})=\boldsymbol{y}({\cal
A})=(\\{\boldsymbol{y}(\boldsymbol{\ell}),\boldsymbol{\ell}\in{\cal
A}\big{|}_{{\cal D}}\\})^{\top}$, $\boldsymbol{X}(\bar{{\cal
A}})=\boldsymbol{X}({\cal
A})=\text{blockdiag}\\{\boldsymbol{X}(\boldsymbol{\ell})^{\top},\boldsymbol{\ell}\in{\cal
A}\big{|}_{{\cal D}}\\}$; values corresponding to unobserved locations will be
dealt with by defining $\widetilde{\boldsymbol{D}}_{n}({\cal A})$ as the
diagonal matrix obtained from $\boldsymbol{D}_{n}$ by replacing unobserved
outcomes with zeros. Denote $\boldsymbol{Z}(\bar{{\cal
A}})=\text{blockdiag}\\{\boldsymbol{Z}(\boldsymbol{\ell}),\boldsymbol{\ell}\in{\cal
A}\big{|}_{{\cal D}}\\}$ and $\boldsymbol{w}(\bar{{\cal A}})$ similarly. If
${\cal A}$ includes $L$ unique spatial locations then
$\boldsymbol{y}(\bar{{\cal A}})$ is a $L\ l\times 1$ vector and
$\boldsymbol{X}({\cal A})$ is a $L\ l\times pl$ matrix. In particular,
$\boldsymbol{Z}(\bar{{\cal A}})$ is a $L\ l\times Lql$ matrix; the subset of
its columns with locations in ${\cal A}$ is denoted as $\boldsymbol{Z}({\cal
A})$ whereas at other locations we get $\boldsymbol{Z}({\cal A}_{-})$. We can
then separate the contribution of $\boldsymbol{w}({\cal A})$ to
$\boldsymbol{y}({{\cal A}})$ from the contribution of $\boldsymbol{w}({\cal
A}_{-})$ by writing $\boldsymbol{y}({{\cal A}})=\boldsymbol{X}({{\cal
A}})\bbeta+\boldsymbol{Z}({\cal A}_{-})\boldsymbol{w}({\cal
A}_{-})+\boldsymbol{Z}({\cal A})\boldsymbol{w}({\cal
A})+\boldsymbol{\varepsilon}({{\cal A}})$, using which we let
$\widetilde{\boldsymbol{y}}({{\cal A}})=\boldsymbol{y}({{\cal
A}})-\boldsymbol{X}({{\cal A}})\bbeta-\boldsymbol{Z}({\cal
A}_{-})\boldsymbol{w}({\cal A}_{-})$.

With customary prior distributions $\bbeta\sim
N(\mathbf{0},\boldsymbol{V}_{\beta})$ and $\tau^{2}_{j}\sim
Inv.Gamma(a_{\tau},b_{\tau})$ along with a Gaussian SpamTree prior on
$\boldsymbol{w}$, we obtain the posterior distribution as

$\displaystyle
p(\boldsymbol{w},\bbeta,\\{\tau_{j}^{2}\\}_{j=1}^{l},\boldsymbol{\theta}\,|\,\boldsymbol{y})$
$\displaystyle\propto
p(\boldsymbol{y}\,|\,\boldsymbol{w},\bbeta,\\{\tau_{j}^{2}\\}_{j=1}^{l})p(\boldsymbol{w}\,|\,\boldsymbol{\theta})p(\boldsymbol{\theta})p(\bbeta)\prod_{j=1}^{l}p(\tau^{2}_{j}).$
(18)

We compute the full conditional distributions of unknowns in the model, save
for $\boldsymbol{\theta}$; iterating sampling from each of these distributions
corresponds to a Gibbs sampler which ultimately leads to samples from the
posterior distribution above.

####  3.2.1 Full conditional distributions

The full conditional distribution for $\bbeta$ is Gaussian with covariance
$\boldsymbol{\Sigma}^{*}_{\bbeta}=(\boldsymbol{V}_{\bbeta}^{-1}+\boldsymbol{X}^{\top}\boldsymbol{D}_{n}^{-1}\boldsymbol{X})^{-1}$
and mean
$\mu^{*}_{\bbeta}=\boldsymbol{\Sigma}_{\bbeta}\boldsymbol{X}^{\top}\boldsymbol{D}_{n}^{-1}(\boldsymbol{y}-\boldsymbol{Z}\boldsymbol{w})$.
For $j=1,\dots,l$,
$p(\tau^{2}_{j}\,|\,\bbeta,\boldsymbol{y},\boldsymbol{w})=Inv.Gamma(a_{\tau,j}^{*},b^{*}_{\tau,j})$
where $a_{\tau,j}^{*}=a_{\tau}+N_{j}/2$ and
$b^{*}_{\tau,j}=b_{\tau}+\frac{1}{2}\boldsymbol{E}^{(j)\top}\boldsymbol{E}^{(j)}$
with
$\boldsymbol{E}^{(j)}=\boldsymbol{y}^{(j)}-\boldsymbol{X}^{(j)}\bbeta_{j}-\boldsymbol{Z}^{(j)}\boldsymbol{w}^{(j)}$.

Take a node $\boldsymbol{v}_{i}\in\boldsymbol{V}$. If
$\boldsymbol{v}_{i}\in\boldsymbol{A}$ then
$\eta^{-1}(\boldsymbol{v}_{i})=S_{i}$ and for
$\boldsymbol{v}_{j}\in\text{Ch}[{\boldsymbol{v}_{i}}]$ denote
$\widetilde{\boldsymbol{w}}_{j}=\boldsymbol{w}_{j}-\boldsymbol{H}_{\setminus
i\to j}\boldsymbol{w}_{[\setminus i\to j]}$. The full conditional distribution
of $\boldsymbol{w}_{i}$ is $N(\boldsymbol{\mu}_{i},\boldsymbol{\Sigma}_{i})$,
where

$\displaystyle\boldsymbol{\Sigma}_{i}^{-1}=\boldsymbol{Z}(S_{i})^{\top}\boldsymbol{D}_{n}(S_{i})^{-1}$
$\displaystyle\boldsymbol{Z}(S_{i})+\boldsymbol{R}_{i}^{-1}+\boldsymbol{F}^{(c)}_{i}$
(19)
$\displaystyle\boldsymbol{\Sigma}_{i}^{-1}\boldsymbol{\mu}_{i}=\boldsymbol{Z}(S_{i})^{\top}\boldsymbol{D}_{n}(S_{i})^{-1}$
$\displaystyle\widetilde{\boldsymbol{y}}(S_{i})+\boldsymbol{R}_{i}^{-1}\boldsymbol{H}_{i}\boldsymbol{w}_{[i]}+\boldsymbol{m}^{(c)}_{i}$
$\displaystyle\boldsymbol{F}^{(c)}_{i}=\sum_{j:\\{\boldsymbol{v}_{j}\in\text{Ch}[{\boldsymbol{v}_{i}}]\\}}\boldsymbol{H}_{i\to
j}^{\top}\boldsymbol{R}_{j}^{-1}\boldsymbol{H}_{i\to j}$
$\displaystyle\qquad\boldsymbol{m}^{(c)}_{i}=\sum_{j:\\{\boldsymbol{v}_{j}\in\text{Ch}[{\boldsymbol{v}_{i}}]\\}}\boldsymbol{H}_{i\to
j}^{\top}\boldsymbol{R}_{j}^{-1}\widetilde{\boldsymbol{w}}_{j}$

If $\boldsymbol{v}_{i}\in\boldsymbol{B}$ instead
$\boldsymbol{\Sigma}_{i}=(\boldsymbol{Z}(U_{i})^{\top}\boldsymbol{D}_{n}(U_{i})^{-1}\boldsymbol{Z}(U_{i})+\boldsymbol{R}_{i})^{-1}$
and
$\boldsymbol{\mu}_{i}=\boldsymbol{\Sigma}_{i}(\boldsymbol{Z}(U_{i})^{\top}\boldsymbol{D}_{n}(U_{i})^{-1}\widetilde{\boldsymbol{y}}(U_{i})+\boldsymbol{R}_{i}^{-1}\boldsymbol{H}_{i}\boldsymbol{w}_{[i]})$.
Sampling of $\boldsymbol{w}$ at nodes at the same level $r$ proceeds in
parallel given the assumed conditional independence structure in ${\cal G}$.
It is thus essential to minimize the computational burden at levels with a
small number of nodes to avoid bottlenecks. In particular computing
$\boldsymbol{F}^{(c)}_{i}$ and $\boldsymbol{m}^{(c)}_{i}$ can become expensive
at the root when the number of children is very large. In Algorithm 3 we show
that one can efficiently sample at a near-root node $\boldsymbol{v}_{i}$ by
updating $\boldsymbol{F}^{(c)}_{i}$ and $\boldsymbol{m}^{(c)}_{i}$ via
message-passing from the children of $\boldsymbol{v}_{i}$.

Initialize: $\text{\Large$\ell$}=0$;

for  _$r\in\\{0,\dots,M\\}$_ do

for  _$j:\\{\boldsymbol{v}_{j}\in\boldsymbol{V}_{r}\\}$_ do // [parallel for]

Compute
$\boldsymbol{R}_{j}^{-1}=(\boldsymbol{C}_{j}-\boldsymbol{C}_{j,[j]}\boldsymbol{C}_{[j]}^{-1}\boldsymbol{C}_{[j],j})^{-1}$
and $|\boldsymbol{R}_{j}^{-1}|$;

$\text{\Large$\ell$}=\text{\Large$\ell$}+\frac{1}{2}\log|\boldsymbol{R}_{j}^{-1}|-\frac{1}{2}(\boldsymbol{w}_{j}-\boldsymbol{H}_{j}\boldsymbol{w}_{[j]})^{\top}\boldsymbol{R}_{j}^{-1}(\boldsymbol{w}_{j}-\boldsymbol{H}_{j}\boldsymbol{w}_{[j]})$;

if  _$\text{Ch}[{\boldsymbol{v}_{j}}]\neq\emptyset$_ then

Identify $\boldsymbol{v}_{i}\in\text{Ch}[{\boldsymbol{v}_{j}}]$ such that
$\boldsymbol{v}_{i}\in\boldsymbol{V}_{r+1}$;

Compute and store $\boldsymbol{C}_{[i]}^{-1}$ (possibly via (13));

Result: $\exp(\text{\Large$\ell$})\propto
p(\boldsymbol{w}\,|\,\boldsymbol{\theta})=\prod_{i}N(\boldsymbol{w}_{i}\,|\,\boldsymbol{H}_{i}\boldsymbol{w}_{[i]},\boldsymbol{R}_{i})$.

Algorithm 1 Computing $p(\boldsymbol{w}\,|\,\boldsymbol{\theta})$.

Input:$\boldsymbol{C}_{[j]}$ for all $j$ from Algorithm 1;

$\boldsymbol{W}_{e}=\bigcup\limits_{r\text{ is even}}\boldsymbol{V}_{r}$;
$\boldsymbol{W}_{o}=\bigcup\limits_{r\text{ is odd}}\boldsymbol{V}_{r}$;

for  _$i\in\\{e,o\\}$_ do

for  _$j:\\{\boldsymbol{v}_{j}\in\boldsymbol{W}_{i}\\}$_ do // [parallel for]

Sample $\boldsymbol{w}_{j}\sim
N(\boldsymbol{\mu}_{j},\boldsymbol{\Sigma}_{j})$ using (19);

Let $\text{Pa}[{\boldsymbol{v}_{j}}]=\\{\boldsymbol{v}_{p}\\}$, then
$\boldsymbol{m}^{(c)}_{p}=\boldsymbol{H}^{\top}_{j}\boldsymbol{R}^{-1}_{j}\boldsymbol{w}_{j}$
and
$\boldsymbol{F}^{(c)}_{p}=\boldsymbol{H}^{\top}_{j}\boldsymbol{R}^{-1}_{j}\boldsymbol{H}_{j}$;

Result: sample from
$p(\boldsymbol{w}_{j}\,|\,\boldsymbol{w}_{-j},\boldsymbol{y},\bbeta,\boldsymbol{\theta},\boldsymbol{\tau})$
for all $\boldsymbol{v}_{j}\in\boldsymbol{V}$.

Algorithm 2 Sampling from the full conditional distribution of
$\boldsymbol{w}_{i}$ when $\delta=1$.

Input : $\boldsymbol{C}_{[j]}$ for all $j$ from Algorithm 1

Initialize: for all $i$, $\boldsymbol{m}^{(c)}_{i}=\mathbf{0}_{n_{i}\times 1}$
and $\boldsymbol{F}^{(c)}_{i}=\boldsymbol{O}_{n_{i}\times n_{i}}$;

for  _$r\in\\{M,\dots,0\\}$_ do

for  _$j:\\{\boldsymbol{v}_{j}\in\boldsymbol{V}_{r}\\}$_ do // [parallel for]

Sample $\boldsymbol{w}_{j}\sim
N(\boldsymbol{\mu}_{j},\boldsymbol{\Sigma}_{j})$ using (19);

for  _$p:\\{\boldsymbol{v}_{p}\in\text{Pa}[{\boldsymbol{v}_{j}}]\\}$_ do

$\boldsymbol{m}^{(c)}_{p}=\boldsymbol{m}^{(c)}_{p}+\boldsymbol{H}^{\top}_{p\to
j}\boldsymbol{R}^{-1}_{j}\boldsymbol{w}_{j}$;

$\boldsymbol{F}^{(c)}_{p}=\boldsymbol{F}^{(c)}_{p}+\boldsymbol{H}^{\top}_{p\to
j}\boldsymbol{R}^{-1}_{j}\boldsymbol{H}_{p\to j}$;

Result: sample from
$p(\boldsymbol{w}_{j}\,|\,\boldsymbol{w}_{-j},\boldsymbol{y},\bbeta,\boldsymbol{\theta},\boldsymbol{\tau})$
for all $\boldsymbol{v}_{j}\in\boldsymbol{V}$.

Algorithm 3 Sampling from the full conditional distribution of
$\boldsymbol{w}_{j}$ when $\delta=M$.

####  3.2.2 Update of $\boldsymbol{\theta}$

The full conditional distribution of $\boldsymbol{\theta}$—which may include
$\boldsymbol{\xi}_{j}$ for $j=1,\dots,q$ or equivalently
$\delta_{ij}=\|\boldsymbol{\xi}_{i}-\boldsymbol{\xi}_{j}\|$ if the chosen
cross-covariance function is defined on a latent domain of variables—is not
available in closed form and sampling a posteriori can proceed via Metropolis-
Hastings steps which involve accept/reject steps with acceptance probability
$\alpha=\min\\{1,\frac{p(\boldsymbol{w}\,|\,\boldsymbol{\theta}^{\prime})p(\boldsymbol{\theta}^{\prime})q(\boldsymbol{\theta}\,|\,\boldsymbol{\theta}^{\prime})}{p(\boldsymbol{w}\,|\,\boldsymbol{\theta})p(\boldsymbol{\theta})q(\boldsymbol{\theta}^{\prime}\,|\,\boldsymbol{\theta})}\\}$.
In our implementation, we adaptively tune the standard deviation of the
proposal distribution via the robust adaptive Metropolis algorithm (RAM;
Vihola, , 2012). In these settings, unlike similar models based on DAG
representations such as NNGPs and MGPs, direct computation via
$p(\boldsymbol{w}\,|\,\boldsymbol{\theta})=\prod_{i}N(\boldsymbol{w}_{i}\,|\,\boldsymbol{H}_{i}\boldsymbol{w}_{[i]},\boldsymbol{R}_{i})$
is inefficient as it requires computing $\boldsymbol{C}^{-1}_{[i]}$ whose size
grows along the hierarchy in ${\cal G}$. We thus outline Algorithm 1 for
computing $p(\boldsymbol{w}\,|\,\boldsymbol{\theta})$ via (13). As an
alternative we can perform the update using ratios of
$p(\boldsymbol{y}\,|\,\bbeta,\boldsymbol{\theta},\boldsymbol{\tau})=\int
p(\boldsymbol{y}\,|\,\boldsymbol{w},\bbeta,\boldsymbol{\tau})p(\boldsymbol{w}\,|\,\boldsymbol{\theta})d\boldsymbol{w}=N(\boldsymbol{y}\,|\,\boldsymbol{X}\bbeta,\boldsymbol{Z}\widetilde{\boldsymbol{C}}\boldsymbol{Z}^{\top}+\boldsymbol{D}_{n})$
using Algorithms 4 and 5 outlined in Appendix B.2.3 which require no sparse
matrix library.

####  3.2.3 Graph coloring for parallel sampling

An advantage of the treed structure of ${\cal G}$ is that it leads to fixed
graph coloring associated to parallel Gibbs sampling; no graph coloring
algorithms are necessary (see e.g. Molloy and Reed, , 2002; Lewis, , 2016).
Specifically, if $\delta=M$ (full depth) then there is a one to one
correspondence between the $M+1$ levels of ${\cal G}$ and graph colors, as
evidenced by the parallel blocks in Algorithms 1 and 3. In the case
$\delta=1$, ${\cal G}$ is associated to only two colors alternating the odd
levels with the even ones. This is possible because the Markov blanket of each
node at level $r$, with $r$ even, only includes nodes at odd levels, and vice-
versa.

####  3.2.4 Prediction of the outcome at new locations

The Gibbs sampling algorithm will iterate across the above steps and, upon
convergence, will produce samples from
$p(\bbeta,\\{\tau^{2}_{j}\\}_{j=1}^{q},\boldsymbol{w}\mid\boldsymbol{y})$. We
obtain posterior predictive inference at arbitrary $\boldsymbol{\ell}\in{\cal
D}$ by evaluating $p(\boldsymbol{y}(\boldsymbol{\ell})\,|\,\boldsymbol{y})$.
If $\boldsymbol{\ell}\in{\cal S}\cup{\cal U}$, then we draw one sample of
$\boldsymbol{y}(\boldsymbol{\ell})\sim
N(\boldsymbol{X}(\boldsymbol{\ell})^{\top}\bbeta+\boldsymbol{Z}(\boldsymbol{\ell})^{\top}\boldsymbol{w}(\boldsymbol{\ell}),\boldsymbol{D}_{n}(\boldsymbol{\ell}))$
for each draw of the parameters from
$p(\bbeta,\\{\tau^{2}_{j}\\}_{j=1}^{l},\boldsymbol{w}\mid\boldsymbol{y})$.
Otherwise, considering that
$\eta(\boldsymbol{\ell})=\boldsymbol{v}_{j}\in\boldsymbol{B}$ for some $j$,
with parent nodes $\text{Pa}[{\boldsymbol{v}_{j}}]$, we sample
$\boldsymbol{w}(\boldsymbol{\ell})$ from the full conditional
$N(\boldsymbol{\mu}_{\boldsymbol{\ell}}^{*},\boldsymbol{\Sigma}_{\boldsymbol{\ell}}^{*})$,
where
$\boldsymbol{\Sigma}_{\boldsymbol{\ell}}^{*}=(\boldsymbol{Z}(\boldsymbol{\ell})\boldsymbol{D}_{n}(\boldsymbol{\ell})^{-1}\boldsymbol{Z}(\boldsymbol{\ell})^{\top}+\boldsymbol{R}_{\boldsymbol{\ell}}^{-1})^{-1}$
and
$\boldsymbol{\mu}_{\boldsymbol{\ell}}^{*}=\boldsymbol{\Sigma}_{\boldsymbol{\ell}}^{*}(\boldsymbol{Z}(\boldsymbol{\ell})\boldsymbol{D}^{-1}(\boldsymbol{y}(\boldsymbol{\ell})-\boldsymbol{X}(\boldsymbol{\ell})^{\top}\bbeta)+\boldsymbol{R}_{\boldsymbol{\ell}}^{-1}\boldsymbol{H}_{\boldsymbol{\ell}}\boldsymbol{w}_{[{j}]})$,
then draw $\boldsymbol{y}(\boldsymbol{\ell})\sim
N(\boldsymbol{X}(\boldsymbol{\ell})^{\top}\bbeta+\boldsymbol{Z}(\boldsymbol{\ell})^{\top}\boldsymbol{w}(\boldsymbol{\ell}),\boldsymbol{D}_{n})$.

####  3.2.5 Computing and storage cost

The update of $\tau_{j}^{2}$ and $\bbeta$ can be performed at a minimal cost
as typically $p=\sum_{j=1}^{l}p_{j}$ is small; almost all the computation
budget must be dedicated to computing
$p(\boldsymbol{w}\,|\,\boldsymbol{\theta})$ and sampling
$p(\boldsymbol{w}\,|\,\boldsymbol{y},\bbeta,\boldsymbol{\tau}^{2})$. Assume
that reference locations are all observed ${\cal S}\subset{\cal T}$ and that
all reference subsets have the same size i.e. $|S_{i}|=N_{s}$ for all $i$. We
show in Appendix B.5 that the cost of computing SpamTrees is $O(nN_{s}^{2})$.
As a result, SpamTrees compare favorably to other models specifically in not
scaling with the cube of the number of samples. $\delta$ does not impact the
computational order, however, compared to $\delta=M$, choosing $\delta=1$
lowers the cost by a factor of $M$ or more. For a fixed reference set
partition and corresponding nodes, choosing larger $\delta$ will result in
stronger dependence between leaf nodes and nodes closer to the root, and this
typically corresponds to leaf nodes being assigned conditioning sets that span
larger distances in space. The computational speedup corresponding to choosing
$\delta=1$ can effectively be traded for a coarser partitioning of ${\cal S}$,
resulting in large conditioning sets that are more local to the leaves.

##  4 Applications

We consider Gaussian SpamTrees for the multivariate regression model (17).
Consider the spatial locations
$\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}\in{\cal D}$ and the locations of
variables $i$ and $j$ in the latent domain of variables
$\boldsymbol{\xi}_{i},\boldsymbol{\xi}_{j}\in\Xi$, then denote
$\boldsymbol{h}=\|\boldsymbol{\ell}-\boldsymbol{\ell}^{\prime}\|$,
$\Delta=\delta_{ij}=\|\boldsymbol{\xi}_{i}-\boldsymbol{\xi}_{j}\|$, and

$\displaystyle
C(\boldsymbol{h},\Delta)=\frac{\exp\left\\{-\phi\|\boldsymbol{h}\|/\exp\left\\{\frac{1}{2}\beta\log(1+\alpha\Delta)\right\\}\right\\}}{\exp\left\\{\beta\log(1+\alpha\Delta)\right\\}}.$

For $j=1,\dots,q$ we also introduce
$C_{j}(\boldsymbol{h})=\exp\left\\{-\phi_{j}\|\boldsymbol{h}\|\right\\}$. A
non-separable cross-covariance function for a multivariate process can be
defined as

$\displaystyle\text{Cov}(w(\boldsymbol{\ell},\boldsymbol{\xi}_{i}),w(\boldsymbol{\ell}^{\prime},\boldsymbol{\xi}_{j}))=\boldsymbol{C}_{ij}(\boldsymbol{h})$
$\displaystyle=\begin{cases}\sigma^{2}_{i1}C(\boldsymbol{h},\delta_{ij})+\sigma^{2}_{i2}C_{i}(\boldsymbol{h})&\text{if
}i=j\\\ \sigma_{i1}\sigma_{j1}C(\boldsymbol{h},\delta_{ij})&\text{if }i\neq
j,\end{cases}$ (20)

which is derived from eq. (7) of Apanasovich and Genton,  (2010); locations of
variables in the latent domain are unknown, therefore
$\boldsymbol{\theta}=\\{\sigma_{i1},\sigma_{i2},\phi_{i}\\}_{i=1,\dots,q}\cup\\{\delta_{ij}\\}_{i=1,\dots,q}^{j<i}\cup\\{\alpha,\beta,\phi\\}$
for a total of $3q+q(q-1)/2+3$ unknown parameters.

###  4.1 Synthetic data

In this section we focus on bivariate outcomes ($q=2$) using the cross-
covariance (20). For simplicity we set $\delta_{21}=1,\alpha=1,\beta=1$. For
each combination of $\sigma_{ij}\in\\{1,2\\}$ for $i,j\in\\{1,2\\}$,
$\phi_{1}=\phi_{2}\in\\{0.1,1,10\\}$, and $\phi\in\\{0.1,1,10\\}$ we generate
25 data sets. Considering model (17), we set $\bbeta=\mathbf{0}$,
$\boldsymbol{Z}=I_{q}$ and take the sampling locations as a regular grid of
size $70\times 70$ for a total of 4,900 spatial locations.

Figure 3: Left half: Full data set – a bivariate outcome is generated on 4,900
spatial locations. Right half: Observed data set – the training sample is
built via subsampling each outcome at a smaller set of locations.

We simulate the spatial effects by sampling the full GP; the nuggets for the
two outcomes are set to $\tau^{2}_{1}=0.01$ and $\tau^{2}_{2}=0.1$. As a
result, we obtain 900 data sets, each simulating a realization of a bivariate
outcome. We mimick real-world data measured irregularly in space by replacing
the outcomes with missing values at $\approx 80\%$ of the spatial locations
chosen uniformly at random and independently across the two margins of the
outcome. In order to replicate the occurrence of regions with more sparsely
observed outcomes, we replace outcomes with missing values at $\approx 99\%$
of spatial locations inside small circular areas whose center is chosen
uniformly at random in $[0,1]^{2}$. Figure 3 shows one of the resulting 900
data sets which we use to evaluate the performance of SpamTrees. We consider
multivariate SpamTrees with $\delta=M$ and $\delta=1$ and compare them with
multivariate cubic meshed GPs (Q-MGPs; Peruzzi et al., , 2020), integrated
nested Laplace approximations (INLA; Rue et al., , 2009) implemented via
R-INLA using a $15\times 15$ grid, a low-rank multivariate GP method (labeled
lowrank) on 25 knots obtained via SpamTrees by setting $M=1$ with no domain
partitioning, and an independent partitioning GP method (labeled ind-part)
implemented by setting $M=1$ and partitioning the domain into 25 regions.
Refer e.g. to Heaton et al.,  (2019) for an overview of low-rank and
independent partitioning methods. We also include results from a non-spatial
regression using Bayesian additive regression trees (BART; Chipman et al., ,
2010). Each method was setup to target a compute time of 15 seconds for each
data set. The total runtime for the 900 data sets thus amounted to about 26
hours.

Tables 1 and 2 summarise the results across all 900 data sets. All Bayesian
methods based on latent GPs exhibit very good coverage; in these simulated
scenarios, SpamTrees exhibit comparatively lower out-of-sample prediction
errors. We highlight that the construction of DAG-based Bayesian methods for
spatial regression depends on the underlying covariance function or kernel;
comparisons of different covariance specifications across different methods is
beyond the scope of this article. Additional implementation details and
figures can be found in Appendix C.1.1.

Model | Cov. (95%) | RMSE$(\boldsymbol{y})$ | MAE$(\boldsymbol{y})$
---|---|---|---
SpamTrees $\delta=M$ | 96.14 | 0.7171 | 0.5144
inla | 92.83 | 0.7521 | 0.5742
q-mgp | 95.87 | 0.7529 | 0.5410
SpamTrees $\delta=1$ | 96.31 | 0.8027 | 0.5737
bart | 92.61 | 0.8820 | 0.6895
lowrank | 96.31 | 1.0230 | 0.7676
ind-part | 95.97 | 1.0796 | 0.8198

Table 1: Prediction performance on multivariate synthetic data: average
coverage of 95% prediction intervals, root mean square error (RMSE), and mean
absolute error in prediction (MAE), over 900 data sets, sorted by lowest RMSE

Model | RMSE$(\boldsymbol{y})$ | MAE$(\boldsymbol{y})$
---|---|---
SpamTrees $\delta=M$ | 2.2882 | 1.6873
lowrank | 2.6953 | 1.4895
q-mgp | 2.8944 | 1.9630
ind-part | 2.9604 | 1.6079
SpamTrees $\delta=1$ | 3.5883 | 2.3322

Table 2: RMSE and MAE in the estimation of $\boldsymbol{\theta}$ from
covariance function (20), averaged over 900 data sets.

###  4.2 Climate data: MODIS-TERRA and GHCN

Climate data are collected from multiple sources in large quantities; when
originating from satellites and remote sensing, they are typically collected
at high spatial and relatively low temporal resolution. Atmospheric and land-
surface products are obtained via post-processing of satellite imaging, and
their quality is negatively impacted by cloud cover and other atmospheric
disturbances. On the other hand, data from a relatively small number of land-
based stations is of low spatial but high temporal resolution. An advantage of
land-based stations is that they measure phenomena related to atmospheric
conditions which cannot be easily measured from satellites (e.g. precipitation
data, depth of snow cover).

We consider the joint analysis of five spatial outcomes collected from two
sources. First, we consider Moderate Resolution Imaging Spectroradiometer
(MODIS) data from the Terra satellite which is part of the NASA’s Earth
Observing System. Specifically, data product MOD11C3 v. 6 provides monthly
Land Surface Temperature (LST) values in a 0.05 degree latitude/longitude grid
(the Climate Modeling Grid or CMG). The monthly data sets cover the whole
globe from 2000-02-01 and consist of daytime and nighttime LSTs, quality
control assessments, in addition to emissivities and clear-sky observations.
The second source of data is the Global Historical Climatology Network (GHCN)
database which includes climate summaries from land surface stations across
the globe subjected to common quality assurance reviews. Data are published by
the National Centers of Environmental Information (NCEI) of the National
Oceanic and Atmospheric Administration (NOAA) at several different temporal
resolutions; daily products report five core elements (precipitation,
snowfall, snow depth, maximum and minimum temperature) in addition to several
other measurements.

We build our data set for analysis by focusing on the continental United
States in October, 2018. The MODIS data correspond to 359,822 spatial
locations. Of these, 250,874 are collected at the maximum reported quality; we
consider all remaining 108,948 spatial locations as missing, and extract (1)
daytime LST (LST_Day_CMG), (2) nighttime LST (LST_Night_CMG), (3) number of
days with clear skies (Clear_sky_days), (4) number of nights with clear skies
(Clear_sky_nights). From the GHCN database we use daily data to obtain monthly
averages for precipitation (PRCP), which is available at 24,066 spatial
locations corresponding to U.S. weather stations; we log-transform PRCP. The
two data sources do not share measurement locations as there is no overlap
between measurement locations in MODIS and GHCN, with the latter data being
collected more sparsely—this is a scenario of complete spatial misalignment.
For this reason we build SpamTrees favoring placement of GHCN locations at
root nodes following Proposition 2.1. Further implementation details are
outlined at Appendix C.

From the resulting data set of size $n=$1,027,562 we remove all observations
in a large $3\times 3$ degree area in the central U.S. (from -100W to -97W and
from 35N to 38N, i.e. the red area of Figure 4) to build a test set on which
we calculate coverage, MAE and RMSE of the predictions. We implement SpamTrees
on the covariance function (20).

Figure 5 maps the predictions at all locations and the corresponding posterior
uncertainties. Comparisons with other methods are difficult due to data size
and complete misalignment. We implemented a tessellated MGP with the same
covariance function and targeting similar computing time (predictive
performance is reported in Appendix C.1.2); SpamTrees displayed lower
prediction errors and better coverage as seen in Table 3. We report posterior
summaries of $\boldsymbol{\theta}$ in Appendix C.1.2. Opposite signs of
$\sigma_{i1}$ and $\sigma_{j1}$ for pairs of variables $i,j\in\\{1,\dots,q\\}$
imply a negative relationship; however, the degree of spatial decay of these
correlations is different for each pair as prescribed by the latent distances
in the domain of variables $\delta_{ij}$. Figure 6 depicts the resulting
cross-covariance function for three pairs of variables.

Figure 4: Prediction area | Measure | Clear_sky_days | Clear_sky_nights | LST_Day_CMG | LST_Night_CMG | PRCP
---|---|---|---|---|---
95% Coverage | 0.9798 | 0.9894 | 1.0000 | 0.9993 | 0.9717
MAE | 1.2824 | 1.3029 | 0.9686 | 0.8440 | 0.3517
RMSE | 1.6114 | 1.6214 | 1.2547 | 1.0764 | 0.5168

|  $n=$1,014,017 | Total iterations: 30,000 | Total time: 16.14h | Average time/iteration: 1.9s
---|---|---|---

Table 3: Prediction results of SpamTrees over the $3\times 3$ degree area
shown in Figure 4 Figure 5: Predicted values of the outcomes at all locations
(top row) and associated 95% uncertainty (bottom row), with darker spots
corresponding to wider credible intervals. Figure 6: Given the latent
dimensions $\delta_{ij}$, the color-coded lines represent
$C(\boldsymbol{h},\delta_{ij})$ whereas
$\boldsymbol{C}_{ij}(\boldsymbol{h})=\sigma_{1i}\sigma_{1j}C(\boldsymbol{h},\delta_{ij})$
is shown as a dashed grey line.

##  5 Discussion

In this article, we introduced SpamTrees for Bayesian spatial multivariate
regression modeling and provided algorithms for scalable estimation and
prediction. SpamTrees add significantly to the class of methods for regression
in spatially-dependent data settings. We have demonstrated that SpamTrees
maintain accurate characterization of spatial dependence and scalability even
in challenging settings involving multivariate data that are spatially
misaligned. Such complexities create problems for competing approaches,
including recent DAG-based approaches ranging from NNGPs to MGPs.

One potential concern is the need for users to choose a tree, and in
particular specify the number of locations associated to each node and the
multivariate composition of locations in each node. Although one can
potentially estimate the tree structure based on the data, this would
eliminate much of the computational speedup. We have provided theoretical
guidance based on KL divergence from the full GP and computational cost
associated to different tree structures. This and our computational
experiments lead to practical guidelines that can be used routinely in tree
building. Choosing a tree is simpler than the common task of choosing a neural
network architecture in deep learning, and provides a useful degree of user-
input to refine and improve upon an approach.

We have focused on sampling algorithms for the latent effects because they
provide a general blueprint which may be used for posterior computations in
non-Gaussian outcome models; efficient algorithms for non-Gaussian big
geostatistical data sets are currently lacking and are the focus of ongoing
research. Including time as a dimension is challenging and care must be taken
when building a sparse DAG to avoid unreasonable assumptions on temporal
dependence. For these reasons, future research may be devoted to building
sparse DAG methods combining the advantages of treed structures with e.g.
Markov-type assumptions of conditional independence.

### Acknowledgements

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 856506). This project was also partially funded by grant
R01ES028804 of the United States National Institutes of Health (NIH).

Appendix

##  Appendix A Kolmogorov consistency conditions for SpamTrees

We adapt results from Datta et al., 2016a and Peruzzi et al.,  (2020). Let
$w(\boldsymbol{s}),\boldsymbol{s}\in\mathcal{D}^{*}$ be the univariate
representation in the augmented domain of the multivariate base process
$\\{\boldsymbol{w}(\boldsymbol{\ell}),\boldsymbol{\ell}\in{\cal
D}\subset\Re^{d}\\}$. Fix the reference set ${\cal S}\subset{\cal D}^{*}$ and
let ${\cal
L}=\\{\boldsymbol{\ell}_{1},\dots,\boldsymbol{\ell}_{n}\\}\subset{\cal D}^{*}$
and ${\cal U}={\cal L}\setminus{\cal S}$. Then

$\displaystyle\int\widetilde{p}(\boldsymbol{w}_{{\cal
L}})\prod_{\boldsymbol{\ell}_{i}\in{\cal L}}dw(\boldsymbol{\ell}_{i})$
$\displaystyle=\int\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}}dw(\boldsymbol{s}_{i})\prod_{\boldsymbol{\ell}_{i}\in{\cal
L}}dw(\boldsymbol{\ell}_{i})$
$\displaystyle=\int\widetilde{p}(\boldsymbol{w}_{{\cal
S}})\left(\int\widetilde{p}(\boldsymbol{w}_{{\cal U}}\mid\boldsymbol{w}_{{\cal
S}})\prod_{\boldsymbol{\ell}_{i}\in{\cal
U}}dw(\boldsymbol{\ell}_{i})\right)\prod_{\boldsymbol{\ell}_{i}\in{\cal
S}}dw(\boldsymbol{\ell}_{i})=1,$

hence $\widetilde{p}(\boldsymbol{w}_{{\cal L}})$ is a proper joint density. To
verify the Kolmogorov consistency conditions, take the permutation ${\cal
L}_{\pi}=\\{\boldsymbol{\ell}_{\pi(1)},\dots,\boldsymbol{\ell}_{\pi(n)}\\}$
and call ${\cal U}_{\pi}={\cal L}_{\pi}\setminus{\cal S}$. Clearly ${\cal
U}_{\pi}={\cal L}_{\pi}\setminus{\cal S}={\cal L}\setminus{\cal S}={\cal U}$
and similarly ${\cal S}\setminus{\cal L}_{\pi}={\cal S}\setminus{\cal L}$ so
that

$\displaystyle\widetilde{p}(\boldsymbol{w}_{{\cal L}_{\pi}})$
$\displaystyle=\int\widetilde{p}(\boldsymbol{w}_{{\cal
U}_{\pi}}\mid\boldsymbol{w}_{{\cal S}})\widetilde{p}(\boldsymbol{w}_{{\cal
S}})\prod_{\boldsymbol{s}_{i}\in{\cal S}\setminus{\cal
L}_{\pi}}dw(\boldsymbol{s}_{i})$
$\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{\cal
L}}dw(\boldsymbol{s}_{i})=\widetilde{p}(\boldsymbol{w}_{\cal L})$

implying

$\displaystyle\widetilde{p}(\boldsymbol{w}(\boldsymbol{\ell}_{1}),\dots,\boldsymbol{w}(\boldsymbol{\ell}_{n}))$
$\displaystyle=\widetilde{p}(\boldsymbol{w}(\boldsymbol{\ell}_{\pi(1)}),\dots,\boldsymbol{w}(\boldsymbol{\ell}_{\pi(n)})).$

Next, take a new location location $\boldsymbol{\ell}_{0}\in{\cal D}^{*}$.
Call ${\cal L}_{1}={\cal L}\cup\\{\boldsymbol{\ell}_{0}\\}$. We want to show
that $\int\widetilde{p}(\boldsymbol{w}_{{\cal
L}_{1}})dw(\boldsymbol{\ell}_{0})=\widetilde{p}(\boldsymbol{w}_{{\cal L}})$.
If $\boldsymbol{\ell}_{0}\in{\cal S}$ then ${\cal L}_{1}\setminus{\cal
S}={\cal L}\setminus{\cal S}={\cal U}$ and hence

$\displaystyle\int\widetilde{p}(\boldsymbol{w}_{{\cal
L}_{1}})dw(\boldsymbol{\ell}_{0})$
$\displaystyle=\int\left(\widetilde{p}(\boldsymbol{w}_{{\cal
L}_{1}\setminus{\cal S}}\mid\boldsymbol{w}_{{\cal
S}})\widetilde{p}(\boldsymbol{w}_{{\cal S}})\prod_{\boldsymbol{s}_{i}\in{\cal
S}\setminus{\cal
L}_{1}}dw(\boldsymbol{s}_{i})\right)dw(\boldsymbol{\ell}_{0})$
$\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{\cal
L}}dw(\boldsymbol{s}_{i})=\widetilde{p}(\boldsymbol{w}_{{\cal L}}).$

If $\boldsymbol{\ell}_{0}\notin{\cal S}$ we have

$\displaystyle\int\widetilde{p}(\boldsymbol{w}_{{\cal
L}_{1}})dw(\boldsymbol{\ell}_{0})$
$\displaystyle=\int\left(\int\widetilde{p}(\boldsymbol{w}_{{\cal
L}_{1}\setminus{\cal S}}\mid\boldsymbol{w}_{{\cal
S}})\widetilde{p}(\boldsymbol{w}_{{\cal S}})\prod_{\boldsymbol{s}_{i}\in{\cal
S}\setminus{\cal
L}_{1}}dw(\boldsymbol{s}_{i})\right)dw(\boldsymbol{\ell}_{0})$
$\displaystyle=\int\left(\int\widetilde{p}(\boldsymbol{w}_{{\cal
L}\setminus{\cal S}\cup\\{\boldsymbol{\ell}_{0}\\}}\mid\boldsymbol{w}_{{\cal
S}})\widetilde{p}(\boldsymbol{w}_{{\cal S}})\prod_{\boldsymbol{s}_{i}\in{\cal
S}\setminus{\cal L}}dw(\boldsymbol{s}_{i})\right)dw(\boldsymbol{\ell}_{0})$
$\displaystyle=\int\left(\int\widetilde{p}(\boldsymbol{w}_{\\{\boldsymbol{\ell}_{0}\\}}\mid\boldsymbol{w}_{{\cal
L}\setminus{\cal S}},\boldsymbol{w}_{{\cal
S}})\widetilde{p}(\boldsymbol{w}_{{\cal L}\setminus{\cal
S}}\mid\boldsymbol{w}_{{\cal S}})\widetilde{p}(\boldsymbol{w}_{{\cal
S}})\prod_{\boldsymbol{s}_{i}\in{\cal S}\setminus{\cal
L}}dw(\boldsymbol{s}_{i})\right)dw(\boldsymbol{\ell}_{0})$
$\displaystyle=\int\widetilde{p}(\boldsymbol{w}_{{\cal L}\setminus{\cal
S}}\mid\boldsymbol{w}_{{\cal S}})\widetilde{p}(\boldsymbol{w}_{{\cal
S}})\prod_{\boldsymbol{s}_{i}\in{\cal S}\setminus{\cal
L}}dw(\boldsymbol{s}_{i})\int\widetilde{p}(\boldsymbol{w}_{\\{\boldsymbol{\ell}_{0}\\}}\mid\boldsymbol{w}_{{\cal
S}})dw(\boldsymbol{\ell}_{0})$
$\displaystyle=\int\widetilde{p}(\boldsymbol{w}_{{\cal L}\setminus{\cal
S}}\mid\boldsymbol{w}_{{\cal S}})\widetilde{p}(\boldsymbol{w}_{{\cal
S}})\prod_{\boldsymbol{s}_{i}\in{\cal S}\setminus{\cal
L}}dw(\boldsymbol{s}_{i})$ $\displaystyle=\widetilde{p}(\boldsymbol{w}_{{\cal
L}}).$

##  Appendix B Properties of Gaussian SpamTrees

Consider the treed graph ${\cal G}$ of a SpamTree. In this section, we make no
distinction between reference and non-reference nodes, and instead label
$\boldsymbol{V}_{i}=\boldsymbol{A}_{i}$ for $i=0,\dots,M-1$ and
$\boldsymbol{V}_{M}=\boldsymbol{B}$ so that
$\boldsymbol{V}=\\{\boldsymbol{A},\boldsymbol{B}\\}=\\{\boldsymbol{V}_{0},\dots,\boldsymbol{V}_{M-1},\boldsymbol{V}_{M}\\}$
and the $\boldsymbol{V}_{M}$ are the leaf nodes. Each $\boldsymbol{w}_{i}$ is
$n_{i}\times 1$ and corresponds to $\boldsymbol{v}_{i}\in\boldsymbol{V}_{r}$
for some $r=0,\dots,M$ so that
$\text{Pa}[{\boldsymbol{v}_{i}}]=\\{\boldsymbol{v}_{j_{1}},\dots,\boldsymbol{v}_{j_{r}}\\}$
for some sequence $\\{j_{1},\dots,j_{r}\\}$, and
$\eta^{-1}(\text{Pa}[{\boldsymbol{v}_{i}}])=\\{S_{j_{1}},\dots,S_{j_{r}}\\}$.
Denote the $h$-th parent of $\boldsymbol{v}_{i}$ as
$\text{Pa}[{\boldsymbol{v}_{i}}](h)$.

###  B.1 Building the precision matrix

We can represent each conditional density
$N(\boldsymbol{w}_{i}\mid\boldsymbol{H}_{i}\boldsymbol{w}_{[{i}]},\boldsymbol{R}_{i})$
as a linear regression on $\boldsymbol{w}_{i}$:

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

where each $\boldsymbol{h}_{ij}$ is an $n_{i}\times n_{j}$ coefficient matrix
representing the regression of $\boldsymbol{w}_{i}$ given
$\boldsymbol{w}_{[i]}$, $\boldsymbol{\omega}_{i}\stackrel{{\scriptstyle
ind}}{{\sim}}N(\boldsymbol{0},\boldsymbol{R}_{i})$ for $i=0,1,\ldots,M$, and
each $\boldsymbol{R}_{i}$ is an $n_{i}\times n_{i}$ residual covariance
matrix. We set $\boldsymbol{h}_{ii}=\boldsymbol{O}$ and
$\boldsymbol{h}_{ij}=\boldsymbol{O}$, where $\boldsymbol{O}$ is the matrix of
zeros, whenever $j\notin\\{j_{1},\dots,j_{r}\\}$. Using this representation,
we have
$\boldsymbol{H}_{i}=\left[\boldsymbol{h}_{i,j_{1}},\boldsymbol{h}_{i,j_{2}},\ldots,\boldsymbol{h}_{i,j_{r}}\right]$,
which is an $n_{i}\times J_{i}$ block matrix formed by stacking
$\boldsymbol{h}_{i,j_{k}}$ side by side for $k=1,\dots,r$. Since
$\mbox{E}[\boldsymbol{w}_{i}\,|\,\boldsymbol{w}_{[i]}]=\boldsymbol{H}_{i}\boldsymbol{w}_{[i]}=\boldsymbol{C}_{i,[i]}\boldsymbol{C}_{[i]}^{-1}\boldsymbol{w}_{[i]}$,
we obtain
$\boldsymbol{H}_{i}=\boldsymbol{C}_{i,[i]}\boldsymbol{C}_{[i]}^{-1}$. 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]}^{-1}\boldsymbol{C}_{[i],i}$,
hence all $\boldsymbol{H}_{i}$’s, $\boldsymbol{h}_{ij}$’s, and
$\boldsymbol{R}_{i}$’s can be computed from the base covariance function.

In order to continue building the precision matrix, define the block matrix
$\boldsymbol{\mathcal{H}}=\\{\boldsymbol{h}_{ij}\\}$. We can write

$\boldsymbol{h}_{ij}=\begin{cases}\boldsymbol{O}&\text{ if
}\boldsymbol{v}_{j}\notin\text{Pa}[{\boldsymbol{v}_{i}}]\\\
(\boldsymbol{C}_{i,[i]}\boldsymbol{C}^{-1}_{[i]})(\cdot,h)=\boldsymbol{H}_{i}(\cdot,h)&\text{
if
}\boldsymbol{v}_{j}=\boldsymbol{v}_{j_{h}}\in\text{Pa}[{\boldsymbol{v}_{i}}],\end{cases}$
(22)

where $(\cdot,h)$ refers to the $h$-th block column. More compactly using the
indicator function $\mathbf{1}\\{\cdot\\}$ we have
$\boldsymbol{h}_{ij}=\mathbf{1}\\{\exists
h:\boldsymbol{v}_{j}=\text{Pa}[{\boldsymbol{v}_{i}}](h)\\}(\boldsymbol{C}_{i,[i]}\boldsymbol{C}^{-1}_{[i]})(\cdot,h)$.
If we stack all the $\boldsymbol{h}_{ik}$ horizontally for
$k=0,\dots,M_{S}-1$, we obtain the $n_{i}\times n$ matrix
$\boldsymbol{\mathcal{H}}(i,\cdot)$, which is $i$-th block row of
$\boldsymbol{\mathcal{H}}$. Intuitively, $\boldsymbol{\mathcal{H}}(i,\cdot)$
is a sparse matrix with the coefficients linking the full $\boldsymbol{w}$ to
$\boldsymbol{w}_{i}$, with zero blocks at locations whose corresponding node
is $\boldsymbol{v}_{j}\notin\text{Pa}[{\boldsymbol{v}_{i}}]$. The $i$th block-
row of $\boldsymbol{\mathcal{H}}$ is of size $n_{i}\times n$ but only has $r$
non-empty sub-blocks, with sizes $n_{i}\times n_{j}$ for
$j\in\\{j_{1},\dots,j_{r}\\}$, respectively. Instead, $\boldsymbol{H}_{i}$ is
a dense matrix obtained by dropping all the zero-blocks from
$\boldsymbol{\mathcal{H}}(i,\cdot)$, and stores the coefficients linking
$\boldsymbol{w}_{[{i}]}$ to $\boldsymbol{w}_{i}$. The two are linked as
$\boldsymbol{H}_{i}\boldsymbol{w}_{[{i}]}=\boldsymbol{\mathcal{H}}(i,\cdot)\boldsymbol{w}$.

Since
$\boldsymbol{w}=\boldsymbol{\mathcal{H}}\boldsymbol{w}+\boldsymbol{\omega}$,
$\widetilde{\boldsymbol{C}}=\mbox{var}(\boldsymbol{w})=(\boldsymbol{I}-\boldsymbol{\mathcal{H}})^{-1}\boldsymbol{R}(\boldsymbol{I}-\boldsymbol{\mathcal{H}})^{-\top}$,
where $\boldsymbol{R}=\text{blockdiag}\\{\boldsymbol{R}_{i}\\}$ and
$\boldsymbol{I}-\boldsymbol{\mathcal{H}}$ is block lower-triangular with unit
diagonal, hence non-singular. We find the precision matrix as
$\widetilde{\boldsymbol{C}}^{-1}=(\boldsymbol{I}-\boldsymbol{\mathcal{H}})^{\top}\boldsymbol{R}^{-1}(\boldsymbol{I}-\boldsymbol{\mathcal{H}})$.

###  B.2 Properties of $\widetilde{\boldsymbol{C}}^{-1}$

When not ambiguous, we use the notation $\boldsymbol{X}_{ij}$ to denote the
$(i,j)$ block of $\boldsymbol{X}$. An exception to this is the $(i,j)$ block
of $\widetilde{\boldsymbol{C}}^{-1}$ which we denote as
$\widetilde{\boldsymbol{C}}^{-1}(i,j)$. In SpamTrees,
$\widetilde{\boldsymbol{C}}^{-1}(i,j)$ is nonzero if $i=j$ or if the
corresponding nodes $\boldsymbol{v}_{i}$ and $\boldsymbol{v}_{j}$ are
connected in the moral graph $\mathcal{G}^{m}$, which is an undirected graph
based on ${\cal G}$ in which an edge connects all nodes that share a child.
This means that either (1)
$\boldsymbol{v}_{i}\in\text{Pa}[{\boldsymbol{v}_{j}}]$ or vice-versa, or (2)
there exists $\boldsymbol{a}^{*}$ such that
$\\{\boldsymbol{v}_{i},\boldsymbol{v}_{j}\\}\subset\text{Pa}[{\boldsymbol{a}^{*}}]$.
In SpamTrees, ${\cal G}^{m}={\cal G}$. In fact, suppose there is a node
$\boldsymbol{a}^{*}\in\boldsymbol{v}_{r^{*}}$ such that
$\boldsymbol{a}^{*}\in\text{Ch}[{\boldsymbol{v}_{j}}]\cap\text{Ch}[{\boldsymbol{v}_{k}}]$,
where $\boldsymbol{v}_{j}\in\boldsymbol{v}_{r_{j}}$ and
$\boldsymbol{v}_{k}\in\boldsymbol{v}_{r_{k}}$. By definition of ${\cal G}$
there exists a sequence $\\{i_{1},\dots,i_{r^{*}}\\}$ such that
$\text{Pa}[{\boldsymbol{a}^{*}}]=\\{\boldsymbol{v}_{i_{1}},\dots,\boldsymbol{v}_{i_{r^{*}}}\\}\supset\\{\boldsymbol{v}_{j},\boldsymbol{v}_{k}\\}$,
and furthermore
$\text{Pa}[{\boldsymbol{v}_{i_{h}}}]=\\{\boldsymbol{v}_{i_{1}},\dots,\boldsymbol{v}_{i_{h-1}}\\}$
for $h\leq r^{*}$. This implies that if $j=k$ then
$\boldsymbol{v}_{j}=\boldsymbol{v}_{k}$, whereas if $j>k$ then
$\boldsymbol{v}_{k}\in\text{Pa}[{\boldsymbol{v}_{j}}]$, meaning that no
additional edge is necessary to build ${\cal G}^{m}$.

####  B.2.1 Explicit derivation of $\widetilde{\boldsymbol{C}}^{-1}(i,j)$

Denote
$\boldsymbol{R}^{-1}=\boldsymbol{R}^{-\frac{1}{2}}\boldsymbol{R}^{-\frac{\top}{2}}$,
$\boldsymbol{U}=(\boldsymbol{I}-\boldsymbol{\mathcal{H}})^{\top}\boldsymbol{R}^{-\frac{1}{2}}$,
and define the “common descendants” as
$\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=(\\{\boldsymbol{v}_{i}\\}\cup\text{Ch}[{\boldsymbol{v}_{i}}])\cap(\\{\boldsymbol{v}_{j}\\}\cup\text{Ch}[{\boldsymbol{v}_{j}}])$.
Then consider
$\boldsymbol{a}_{i}\in\boldsymbol{A},\boldsymbol{v}_{j}\in\boldsymbol{V}$ such
that $\boldsymbol{a}_{i}\in\text{Pa}[{\boldsymbol{v}_{j}}]$ and denote as
$\boldsymbol{H}_{i\to j}$ the matrix obtained by subsetting
$\boldsymbol{H}_{j}$ to columns corresponding to $\boldsymbol{a}_{i}$ and note
that $\boldsymbol{H}_{i\to j}=\boldsymbol{\mathcal{H}}_{ji}$. The $(i,j)$
block of $\boldsymbol{U}$ is then

$\boldsymbol{U}_{ij}=\begin{cases}\boldsymbol{O}_{n_{i}\times n_{j}}&\text{ if
}\boldsymbol{v}_{j}\notin\text{Ch}[{\boldsymbol{v}_{i}}]\\\
\boldsymbol{I}_{n_{i}\times n_{i}}&\text{ if }i=j\\\
(\boldsymbol{I}_{ji}-\boldsymbol{H}_{i\to
j})^{\top}\boldsymbol{R}_{j}^{-\frac{1}{2}}&\text{ if
}\boldsymbol{v}_{j}\in\text{Ch}[{\boldsymbol{v}_{i}}]\\\
\quad\quad=(\boldsymbol{I}_{ji}-\boldsymbol{\mathcal{H}}_{ji})^{\top}\boldsymbol{R}_{j}^{-\frac{1}{2}}\end{cases}$

Then
$\widetilde{\boldsymbol{C}}^{-1}(i,j)=\sum_{k}\boldsymbol{U}_{ik}\boldsymbol{U}_{jk}^{\top}$
and, as in (11), each block of the precision matrix is:

$\displaystyle\widetilde{\boldsymbol{C}}^{-1}(i,j)$
$\displaystyle=\sum\limits_{\boldsymbol{v}_{k}\in\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})}(\boldsymbol{I}_{ki}-\boldsymbol{H}_{i\to
k})^{\top}\boldsymbol{R}_{k}^{-1}(\boldsymbol{I}_{kj}-\boldsymbol{H}_{j\to
k})$ (23)
$\displaystyle=\sum\limits_{\boldsymbol{v}_{k}\in\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})}(\boldsymbol{I}_{ki}-\boldsymbol{\mathcal{H}}_{ki})^{\top}\boldsymbol{R}_{k}^{-1}(\boldsymbol{I}_{kj}-\boldsymbol{\mathcal{H}}_{kj})$

where $\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=\emptyset$ implies
$\widetilde{\boldsymbol{C}}^{-1}(i,j)=\boldsymbol{O}$ and
$\boldsymbol{I}_{ij}$ a zero matrix unless $i=j$ as it is the $(i,j)$ block of
an identity matrix of dimension $n\times n$.

####  B.2.2 Computation of large matrix inverses

One important aspect in building $\widetilde{\boldsymbol{C}}^{-1}$ is that it
requires the computation of the inverse $\boldsymbol{C}^{-1}_{[i]}$ of
dimension $J_{i}\times J_{i}$ for all nodes with parents, i.e. at $r>0$.
Unlike models which achieve scalable computations by limiting the size of the
parent set (e.g. NNGPs and their blocked variant, or tessellated MGPs), this
inverse is increasingly costlier when $\delta>1$ for nodes at a higher-level
of the tree as those nodes have more parents and hence larger sets of parent
locations (the same conclusion holds for non-reference nodes). However, the
treed structure in ${\cal G}$ allows one to avoid computing the inverse in
$O(J_{i}^{3})$. In fact, suppose we have a symmetric, positive-definite block-
matrix $\boldsymbol{A}$ and we wish to compute its inverse. We write

$\displaystyle\boldsymbol{A}=\begin{bmatrix}C&B\\\ B^{\top}&D\end{bmatrix}$
$\displaystyle\quad\boldsymbol{A}^{-1}=\begin{bmatrix}C^{-1}+C^{-1}BS^{-1}B^{\top}C^{-1}&-C^{-1}BS^{-1}\\\
-S^{-1}B^{\top}C^{-1}&S^{-1}\end{bmatrix},$

where $S=C-BD^{-1}B^{\top}$ is the Schur complement of $D$ in
$\boldsymbol{A}$. If $C^{-1}$ was available, the only necessary inversion is
that of $S$. In SpamTrees with $\delta>1$, suppose
$\boldsymbol{v}_{i},\boldsymbol{v}_{j}$ are two nodes such that
$\text{Pa}[{\boldsymbol{v}_{j}}]=\\{\boldsymbol{v}_{i}\\}\cup\text{Pa}[{\boldsymbol{v}_{i}}]$
– this arises for nodes $\boldsymbol{v}_{j}\in\boldsymbol{V}_{r},r\geq
M_{\delta}$. Regardless of whether $\boldsymbol{v}_{j}$ is a reference node or
not, $\eta^{-1}(\text{Pa}[{\boldsymbol{v}_{j}}])=\\{{\cal S}_{i},{\cal
S}_{[i]}\\}$ and

$\displaystyle\boldsymbol{C}_{[j]}=\begin{bmatrix}\boldsymbol{C}_{[i]}&\boldsymbol{C}_{[i],i}\\\
\boldsymbol{C}_{i,[i]}&\boldsymbol{C}_{i}\end{bmatrix},$
$\displaystyle\quad\boldsymbol{C}_{[j]}^{-1}=\begin{bmatrix}\boldsymbol{C}_{[i]}^{-1}+\boldsymbol{C}_{[i]}^{-1}\boldsymbol{C}_{[i],i}S^{-1}\boldsymbol{C}_{i,[i]}\boldsymbol{C}_{[i]}^{-1}&-\boldsymbol{C}_{[i]}^{-1}\boldsymbol{C}_{[i],i}S^{-1}\\\
-S^{-1}\boldsymbol{C}_{i,[i]}\boldsymbol{C}_{[i]}^{-1}&S^{-1}\end{bmatrix},$

where the Schur complement of $\boldsymbol{C}_{i}$ is
$S=\boldsymbol{C}_{i}-\boldsymbol{C}_{i,[i]}\boldsymbol{C}_{[i]}^{-1}\boldsymbol{C}_{[i],i}=\boldsymbol{R}_{i}$.
Noting that
$\boldsymbol{H}_{i}=\boldsymbol{C}_{i,[i]}\boldsymbol{C}^{-1}_{[i]}$ we write

$\displaystyle\boldsymbol{C}_{[j]}^{-1}=\begin{bmatrix}\boldsymbol{C}_{[i]}^{-1}+\boldsymbol{H}_{i}^{\top}\boldsymbol{R}_{i}^{-1}\boldsymbol{H}_{i}&-\boldsymbol{H}_{i}^{\top}\boldsymbol{R}_{i}^{-1}\\\
-\boldsymbol{R}_{i}^{-1}\boldsymbol{H}_{i}&\boldsymbol{R}_{i}^{-1}\end{bmatrix}.$
(24)

####  B.2.3 Computing
$(\widetilde{\boldsymbol{C}}^{-1}+\boldsymbol{\Sigma})^{-1}$ and its
determinant without sparse Cholesky

Bayesian estimation of regression models requiring the computation of
$(\boldsymbol{Z}\widetilde{\boldsymbol{C}}\boldsymbol{Z}^{\top}+\boldsymbol{D})^{-1}$
and its determinant use the Sherman-Morrison-Woodbury matrix identity to find
$(\boldsymbol{Z}\widetilde{\boldsymbol{C}}\boldsymbol{Z}^{\top}+\boldsymbol{D})^{-1}=\boldsymbol{D}^{-1}-\boldsymbol{D}^{-1}\boldsymbol{Z}(\widetilde{\boldsymbol{C}}^{-1}+\boldsymbol{\Sigma})^{-1}\boldsymbol{Z}^{\top}\boldsymbol{D}^{-1}$,
where
$\boldsymbol{\Sigma}=\boldsymbol{Z}^{\top}\boldsymbol{D}^{-1}\boldsymbol{Z}$.
A sparse Cholesky factorization of
$\widetilde{\boldsymbol{C}}^{-1}+\boldsymbol{\Sigma}$ can be used as typically
$\boldsymbol{\Sigma}$ is diagonal or block-diagonal, thus maintaining the
sparsity structure of $\widetilde{\boldsymbol{C}}^{-1}$. Sparse Cholesky
libraries (e.g. Cholmod, Chen et al., , 2008), which are embedded in software
or high-level languages such as Matlab™ or the Matrix package for R, scale to
large sparse matrices but are either too flexible or too restrictive in our
use cases: (1) we know ${\cal G}$ and its properties in advance; (2) SpamTrees
take advantage of block structures and grouped data. In fact, sparse matrix
libraries typically are agnostic of ${\cal G}$ and heuristically attempt to
infer a sparse ${\cal G}$ given its moralized counterpart. While this
operation is typically performed once, a priori knowledge of ${\cal G}$
implies that reliance on such libraries is in principle unnecessary.

We thus take advantage of the known structure in ${\cal G}$ to derive direct
algorithms for computing
$(\widetilde{\boldsymbol{C}}^{-1}+\boldsymbol{\Sigma})^{-1}$ and its
determinant. In the discussion below we consider $\delta=M$, noting here that
choosing $\delta=1$ simplifies the treatment as
$\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=\\{\boldsymbol{v}_{i}\\}$ if
$\boldsymbol{v}_{i}=\boldsymbol{v}_{j}$, and it is empty otherwise. We now
show how (23) leads to Algorithm 4 for the decomposition of any precision
matrix $\boldsymbol{\Lambda}$ which conforms to ${\cal G}$ – i.e. it has the
same block-sparse structure as a precision matrix built as in Section B.1.
Suppose from $\boldsymbol{\Lambda}$ we seek a block lower-triangular matrix
$\boldsymbol{L}$ and a block diagonal $\boldsymbol{D}$ such that

$\boldsymbol{\Lambda}_{ij}=\sum\limits_{\boldsymbol{v}_{k}\in\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})}(\boldsymbol{I}_{ki}-\boldsymbol{L}_{ki})^{\top}\boldsymbol{D}_{k}(\boldsymbol{I}_{kj}-\boldsymbol{L}_{kj}).$

Start with $\boldsymbol{v}_{i},\boldsymbol{v}_{j}$ taken from the leaf nodes,
i.e. $\boldsymbol{v}_{i},\boldsymbol{v}_{j}\in\boldsymbol{V}_{M}$. Then
$\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=\emptyset$ and we set
$\boldsymbol{L}_{ij}=\boldsymbol{L}(j,i)^{\top}=\boldsymbol{O}=\boldsymbol{\Lambda}_{ij}$.
If $i=j$ then
$\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{i})=\\{\boldsymbol{v}_{i}\\}$
and

$\displaystyle\sum_{\boldsymbol{v}_{k}\in\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{i})}$
$\displaystyle(\boldsymbol{I}_{ki}-\boldsymbol{L}_{ki})^{\top}\boldsymbol{D}_{k}(\boldsymbol{I}_{ki}-\boldsymbol{L}_{ki})=(\boldsymbol{I}_{ii}-\boldsymbol{L}_{ii})^{\top}\boldsymbol{D}_{i}(\boldsymbol{I}_{ii}-\boldsymbol{L}(i,i))$
$\displaystyle=\boldsymbol{D}_{i}-\boldsymbol{D}_{i}\boldsymbol{L}_{ii}-\boldsymbol{L}_{ii}^{\top}\boldsymbol{D}_{i}+\boldsymbol{L}_{ii}^{\top}\boldsymbol{D}_{i}\boldsymbol{L}_{ii};$

we then set $\boldsymbol{L}_{ii}=\boldsymbol{O}$ and get the $i$-th block of
$\boldsymbol{D}_{i}$ simply setting
$\boldsymbol{D}_{i}=\boldsymbol{\Lambda}_{ii}$. Proceeding downwards along
${\cal G}$, if
$\boldsymbol{v}_{j}\in\boldsymbol{V}_{M-1}\cap\text{Pa}[{\boldsymbol{v}_{i}}]$
we have
$\boldsymbol{\Lambda}_{ij}=\boldsymbol{D}_{i}(\boldsymbol{I}_{ij}-\boldsymbol{L}_{ij})=-\boldsymbol{D}_{i}\boldsymbol{L}_{ij}$
and thus set
$\boldsymbol{L}_{ij}=-\boldsymbol{D}_{i}^{-1}\boldsymbol{\Lambda}_{ij}$. We
then note that
$\text{cd}(\boldsymbol{v}_{j},\boldsymbol{v}_{j})=\\{\boldsymbol{v}_{j},\boldsymbol{v}_{i}\\}$
and obtain
$\boldsymbol{\Lambda}_{jj}=\boldsymbol{D}_{j}+\boldsymbol{L}_{ij}^{\top}\boldsymbol{D}_{i}\boldsymbol{L}_{ij}$
where $\boldsymbol{L}_{ij}$ and $\boldsymbol{D}_{i}$ have been fixed at the
previous step; this results in
$\boldsymbol{D}_{j}=\boldsymbol{\Lambda}_{jj}-\boldsymbol{L}_{ij}^{\top}\boldsymbol{D}_{i}\boldsymbol{L}_{ij}$.

Then, the $s$-th (of $M$) step takes
$\boldsymbol{v}_{j}\in\boldsymbol{V}_{M-s}\cap\text{Pa}[{\boldsymbol{v}_{i}}]$
and $\boldsymbol{v}_{i}\in\boldsymbol{V}_{M-s+1}$, implying
$\text{cd}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=\\{\boldsymbol{v}_{i}\\}\cup\text{Ch}[{\boldsymbol{v}_{i}}]$.
Noting that
$\boldsymbol{F}^{*}=\sum_{\boldsymbol{v}_{k}\in\text{Ch}[{\boldsymbol{v}_{i}}]}(\boldsymbol{I}_{ki}-\boldsymbol{L}_{ki})^{\top}\boldsymbol{D}_{k}(\boldsymbol{I}_{kj}-\boldsymbol{L}_{kj})$
has been fixed at previous steps since each $\boldsymbol{v}_{k}$ is at level
$M-s+2$, we split the sum in (23) and get

$\displaystyle\boldsymbol{\Lambda}_{ij}$
$\displaystyle-\boldsymbol{F}^{*}=\boldsymbol{D}_{i}(\boldsymbol{I}_{ij}-\boldsymbol{L}_{ij})=-\boldsymbol{D}_{i}\boldsymbol{L}_{ij},$

where $\boldsymbol{D}_{i}$ has been fixed at step $s-1$, obtaining
$\boldsymbol{L}_{ij}=-\boldsymbol{D}_{i}^{-1}(\boldsymbol{\Lambda}_{ij}-\boldsymbol{F}^{*})$;
$\boldsymbol{D}_{j}$ can be found using the same logic. Proceeding until
$M-s=0$ from the leaves of ${\cal G}$ to the root, we ultimately fill each
non-empty block in $\boldsymbol{L}$ and $\boldsymbol{D}$ resulting in
$\boldsymbol{\Lambda}=(\boldsymbol{I}-\boldsymbol{L})^{\top}\boldsymbol{D}(\boldsymbol{I}-\boldsymbol{L})$.
Algorithm 4 unifies these steps to obtain the block decomposition of any
sparse precision matrix $\boldsymbol{\Lambda}$ conforming to ${\cal G}$
resulting in
$\boldsymbol{\Lambda}=(\boldsymbol{I}-\boldsymbol{L})^{\top}\boldsymbol{D}(\boldsymbol{I}-\boldsymbol{L})$,
where $\boldsymbol{L}$ is block lower triangular and $\boldsymbol{D}$ is block
diagonal. This is akin to a block-LDL decomposition of $\boldsymbol{\Lambda}$
indexed on nodes of ${\cal G}$. Algorithm 5 complements this decomposition by
providing a ${\cal G}$-specific block version of forward substitution for
computing $(\boldsymbol{I}-\boldsymbol{L})^{-1}$ with $\boldsymbol{L}$ as
above.

In practice, a block matrix with $K^{2}$ blocks can be represented as a
$K^{2}$ array with rows and columns indexed by nodes in ${\cal G}$ and matrix
elements which may be zero-dimensional whenever corresponding to blocks of
zeros. The specification of all algorithms in block notation allows us to
never deal with large (sparse) matrices in practice but only with small block
matrices indexed by nodes in ${\cal G}$, bypassing the need for external
sparse matrix libraries. Specifically we use the above algorithms to compute
$\boldsymbol{\Lambda}^{-1}=(\widetilde{\boldsymbol{C}}^{-1}+\boldsymbol{\Sigma})^{-1}$
and its determinant:
$\boldsymbol{\Lambda}^{-1}=(\boldsymbol{I}-\boldsymbol{L})^{-1}\boldsymbol{D}^{-1}(\boldsymbol{I}-\boldsymbol{L})^{-\top}$
and
$|\boldsymbol{\Lambda}^{-1}|=\prod_{i=1}^{M_{S}}1/|\widetilde{\boldsymbol{D}}_{i}|$.
We have not distinguished non-reference and reference nodes in this
discussion. In cases in which the non-reference set is large, we note that the
conditional independence of all non-reference locations, given their parents,
results in $\widetilde{\boldsymbol{C}}^{-1}(i,i)$ being diagonal for all
$\boldsymbol{\ell}\in{\cal U}$ (i.e.
$\eta(\boldsymbol{\ell})=\boldsymbol{v}_{i}\in\boldsymbol{B}$). This portion
of the precision matrix can just be stored as a column vector.

Input : $\boldsymbol{\Lambda}$ $n\times n$ precision matrix conforming to
${\cal G}$

Initialize $\boldsymbol{L}=\boldsymbol{O}_{n\times
n},\boldsymbol{D}=\boldsymbol{O}_{n\times n}$;

for  _$r\in\\{M,\dots,0\\}$_ do // top down from last level

for  _$j:\\{\boldsymbol{v}_{j}\in\boldsymbol{V}_{r}\\}$_ do // [parallel for]

$\boldsymbol{D}_{jj}=\boldsymbol{\Lambda}_{jj}$;

for  _$p:\\{\boldsymbol{v}_{p}\in\text{Pa}[{\boldsymbol{v}_{j}}]\\}$_ do

$\boldsymbol{L}_{jp}=-\boldsymbol{D}^{-1}_{jj}\boldsymbol{\Lambda}_{jp}$;

for  _$g:\\{\boldsymbol{v}_{g}\in\text{Pa}[{\boldsymbol{v}_{j}}]\\}$ _ do

$\boldsymbol{\Lambda}_{pg}=\boldsymbol{\Lambda}_{pg}-\boldsymbol{\Lambda}_{jp}^{\top}\boldsymbol{L}_{jg}$;

$\boldsymbol{\Lambda}_{gp}=\boldsymbol{\Lambda}_{pg}^{\top}$;

Result: Block-lower-triangular $\boldsymbol{L}$ with
$\boldsymbol{L}_{ij}\neq\boldsymbol{O}$ if
$\boldsymbol{v}_{i}\in\text{Pa}[{\boldsymbol{v}_{j}}]$, and block-diagonal
$\boldsymbol{D}$ such that
$(\boldsymbol{I}-\boldsymbol{L})^{\top}\boldsymbol{D}(\boldsymbol{I}-\boldsymbol{L})=\boldsymbol{\Lambda}$.

Algorithm 4 Precision matrix decomposition given treed graph ${\cal G}$ with
$M$ levels.

Input : $\boldsymbol{\Gamma}=\boldsymbol{I}-\boldsymbol{L}$ where
$\boldsymbol{L}$ is as in Algorithm 4.

Initialize $\boldsymbol{\Delta}_{ij}=\boldsymbol{O}_{n_{i},n_{j}}$ for all
$i,j$ such that $\boldsymbol{v}_{j}\in\text{Pa}[{\boldsymbol{v}_{i}}]$;

for  _$r\in\\{0,\dots,M\\}$_ do // bottom up from root of ${\cal G}$

for  _$j:\\{\boldsymbol{v}_{j}\in\boldsymbol{V}_{r}\\}$_ do // [parallel for]

for  _$p:\\{\boldsymbol{v}_{p}\in\widetilde{{\cal P}}_{0\to[j]}\\}$_ do

Set
$\text{chain}(\boldsymbol{v}_{p},\boldsymbol{v}_{j})=\\{\boldsymbol{v}_{p}\\}\cup\\{\widetilde{{\cal
P}}_{0\to[j]}\cap\widetilde{{\cal P}}_{0\to[p]}\\}$;

for
_$g:\\{\boldsymbol{v}_{g}\in\text{chain}(\boldsymbol{v}_{p},\boldsymbol{v}_{j})\\}$
_ do

$\boldsymbol{\Delta}_{jp}=\boldsymbol{\Delta}_{jp}-\boldsymbol{\Gamma}_{jg}\boldsymbol{\Delta}_{gp}$

Result: $\boldsymbol{\Delta}=\boldsymbol{\Gamma}^{-1}$.

Algorithm 5 Calculating the inverse of $\boldsymbol{I}-\boldsymbol{L}$ with
$\boldsymbol{L}$ output from Algorithm 4.

####  B.2.4 Sparsity of $\widetilde{\boldsymbol{C}}^{-1}$

We calculate the sparsity in the precision matrix; considering an enumeration
of nodes by level in ${\cal G}$, denote
$n_{ij}=|\eta^{-1}(\boldsymbol{v}_{ij})|$, $m_{j}=|\boldsymbol{V}_{j}|$, and
$J_{ij}=|\eta^{-1}(\text{Pa}[{\boldsymbol{v}_{ij}}])|$, and noting that by
symmetry
$(\widetilde{\boldsymbol{C}}^{-1}(i,j))^{\top}=\widetilde{\boldsymbol{C}}^{-1}(j,i)$,
the number of nonzero elements of $\widetilde{\boldsymbol{C}}^{-1}$ is

$\text{nnz}(\widetilde{\boldsymbol{C}}^{-1})=\sum_{j=0}^{M}\sum_{i=1}^{m_{j}}\left(2n_{ij}J_{ij}+n_{ij}^{2}\boldsymbol{1}\\{j<M\\}+n_{ij}\boldsymbol{1}\\{j=M\\}\right),$

where $n_{ij}\boldsymbol{1}\\{j=M\\}$ refers to the diagonal elements of the
precision matrix at non-reference locations.

###  B.3 Properties of SpamTrees with $\delta=M$

We outline recursive properties of $\widetilde{\boldsymbol{C}}$ induced by
${\cal G}$ when $\delta=M$. In the case $1<\delta<M$, these properties hold
for nodes at or above level $M_{\delta}$, using $\boldsymbol{A}_{M_{\delta}}$
as root. We focus on paths in ${\cal G}$. These can be represented as
sequences of nodes $\\{\boldsymbol{v}_{i_{1}},\dots,\boldsymbol{v}_{i_{r}}\\}$
such that
$\\{\boldsymbol{v}_{i_{j}},\dots,\boldsymbol{v}_{i_{k}}\\}\subset\text{Pa}[{\boldsymbol{v}_{i_{k+1}}}]$
for $1<j<k<r$. Take two successive elements of such a sequence, i.e.
$\boldsymbol{v}_{i},\boldsymbol{v}_{j}$ such that
$\boldsymbol{v}_{i}\rightarrow\boldsymbol{v}_{j}$ in ${\cal G}$. Consider
$\mbox{E}[\boldsymbol{w}_{j}\,|\,\boldsymbol{w}_{[j]}]=\boldsymbol{H}_{j}\boldsymbol{w}_{[j]}=\boldsymbol{C}_{j,[j]}\boldsymbol{C}_{[j]}^{-1}\boldsymbol{w}_{[j]}$
and
$\boldsymbol{R}_{j}=\mbox{var}\\{\boldsymbol{w}_{j}\,|\,\boldsymbol{w}_{[j]}\\}=\boldsymbol{C}_{j,j}-\boldsymbol{C}_{j,[j]}\boldsymbol{C}_{[j]}^{-1}\boldsymbol{C}_{[j],j}$.
By (24) we can write

$\displaystyle\boldsymbol{H}_{j}$
$\displaystyle\boldsymbol{w}_{[{j}]}=\begin{bmatrix}\boldsymbol{C}_{j,[i]}&\boldsymbol{C}_{j,i}\end{bmatrix}\begin{bmatrix}\boldsymbol{C}_{[i]}^{-1}+\boldsymbol{H}_{i}^{\top}\boldsymbol{R}_{i}^{-1}\boldsymbol{H}_{i}&-\boldsymbol{H}_{i}^{\top}\boldsymbol{R}_{i}^{-1}\\\
-\boldsymbol{R}_{i}^{-1}\boldsymbol{H}_{i}&\boldsymbol{R}_{i}^{-1}\end{bmatrix}\begin{bmatrix}\boldsymbol{w}_{[{i}]}\\\
\boldsymbol{w}_{i}\end{bmatrix}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\boldsymbol{C}_{j,[i]}&\boldsymbol{C}_{j,i}-\boldsymbol{C}_{j,[i]}\boldsymbol{C}^{-1}_{[i]}\boldsymbol{C}_{[i],i}\end{bmatrix}\begin{bmatrix}\boldsymbol{C}^{-1}_{[i]}&\boldsymbol{O}\\\
\boldsymbol{O}&(\boldsymbol{C}_{i,i}-\boldsymbol{C}_{i,[i]}\boldsymbol{C}^{-1}_{[i]}\boldsymbol{C}_{[i],i})^{-1}\end{bmatrix}\begin{bmatrix}\boldsymbol{w}_{[{i}]}\\\
\boldsymbol{w}_{i}-\boldsymbol{C}_{i,[i]}\boldsymbol{C}^{-1}_{[i]}\boldsymbol{w}_{[i]}\end{bmatrix}$
$\displaystyle=$
$\displaystyle\begin{bmatrix}\boldsymbol{C}_{j,[i]}\boldsymbol{C}^{-1}_{[i]}&(\boldsymbol{C}_{j,i}-\boldsymbol{C}_{j,[i]}\boldsymbol{C}^{-1}_{[i]}\boldsymbol{C}_{[i],i})(\boldsymbol{C}_{i,i}-\boldsymbol{C}_{i,[i]}\boldsymbol{C}^{-1}_{[i]}\boldsymbol{C}_{[i],i})^{-1}\end{bmatrix}\begin{bmatrix}\boldsymbol{w}_{[{i}]}\\\
\boldsymbol{w}_{i}-\boldsymbol{C}_{i,[i]}\boldsymbol{C}^{-1}_{[i]}\boldsymbol{w}_{[i]}\end{bmatrix}.$

Now define the covariance function
$\boldsymbol{K}_{i}(\boldsymbol{\ell},\boldsymbol{\ell}^{\prime})=\boldsymbol{C}_{\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}}-\boldsymbol{C}_{\boldsymbol{\ell},[i]}\boldsymbol{C}^{-1}_{[i]}\boldsymbol{C}_{[i],\boldsymbol{\ell}^{\prime}}$;
recalling that the reference set is ${\cal
S}=\cup_{i=0}^{M-1}\cup_{j=1}^{m_{j}}S_{j}$ we use a shorthand notation for
these subsets: $\boldsymbol{K}_{i}(S_{h},S_{k})=\boldsymbol{K}_{i}(h,k)$. Also
denote
$\boldsymbol{e}_{i}=\boldsymbol{w}_{i}-\boldsymbol{C}_{i,[i]}\boldsymbol{C}^{-1}_{[i]}\boldsymbol{w}_{[{i}]}$
for all $i$. The above expression becomes

$\displaystyle\boldsymbol{H}_{j}\boldsymbol{w}_{[{j}]}$
$\displaystyle=\begin{bmatrix}\boldsymbol{C}_{j,[i]}\boldsymbol{C}^{-1}_{[i]}&\boldsymbol{K}_{i}(j,i)\boldsymbol{K}_{i}^{-1}(i,i)\end{bmatrix}\begin{bmatrix}\boldsymbol{w}_{[{i}]}\\\
\boldsymbol{e}_{i}\end{bmatrix}$
$\displaystyle=\boldsymbol{H}_{i}\boldsymbol{w}_{[{i}]}+\boldsymbol{K}_{i}(j,i)\boldsymbol{K}_{i}^{-1}(i,i)\boldsymbol{e}_{i};$

we can use this recursively on
$\\{\boldsymbol{v}_{i_{0}},\boldsymbol{v}_{i_{1}},\dots,\boldsymbol{v}_{i_{r}}\\}$
where $\boldsymbol{v}_{i_{0}}\in\boldsymbol{A}_{0}$ and
$\boldsymbol{v}_{i_{r}}=\boldsymbol{v}_{j}$ and get

$\displaystyle\boldsymbol{H}_{j}\boldsymbol{w}_{[{j}]}$
$\displaystyle=\sum\limits_{s=i_{1}}^{i_{r-1}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}$
$\displaystyle\mbox{E}[\boldsymbol{w}_{j}\,|\,\boldsymbol{w}_{[{j}]}]$
$\displaystyle=\sum\limits_{s=i_{1}}^{i_{r-1}}\mbox{E}_{\boldsymbol{e}_{s}}[\boldsymbol{w}_{j}\,|\,\boldsymbol{e}_{s}],$

where the expectations on the r.h.s. are taken with respect to the
distributions of $\boldsymbol{e}_{s}$ which are Gaussian with mean zero and
$\mbox{var}\\{\boldsymbol{e}_{h}\\}=\boldsymbol{K}_{h}(h,h)$ – this is a
compact expression of the conditionals governing the process as prescribed by
${\cal G}$. We can also write the above as
$\mbox{E}(\boldsymbol{w}_{j}\,|\,\boldsymbol{w}_{[{j}]})=\sum_{s=i_{0}}^{i_{r}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)(\boldsymbol{w}_{s}-\mbox{E}[\boldsymbol{w}_{s}\,|\,\boldsymbol{w}_{[{s}]}])$;
using
$\mbox{E}[\boldsymbol{e}_{h}\,|\,\boldsymbol{w}_{h},\boldsymbol{w}_{[h]}]=0$,
for $h<k$ we find

$\displaystyle\mbox{cov}\\{\boldsymbol{e}_{h},\boldsymbol{e}_{k}\\}$
$\displaystyle=\mbox{E}[\mbox{cov}\\{\boldsymbol{e}_{h},\boldsymbol{e}_{k}\,|\,\boldsymbol{w}_{h},\boldsymbol{w}_{[h]}\\}]+\mbox{cov}\\{\mbox{E}[\boldsymbol{e}_{h}\,|\,\boldsymbol{w}_{h},\boldsymbol{w}_{[h]}],\mbox{E}[\boldsymbol{e}_{k}\,|\,\boldsymbol{w}_{h},\boldsymbol{w}_{[h]}]\\}$
$\displaystyle=\mbox{cov}\\{\mbox{E}[\boldsymbol{e}_{h}\,|\,\boldsymbol{w}_{h},\boldsymbol{w}_{[h]}],\mbox{E}[\boldsymbol{e}_{k}\,|\,\boldsymbol{w}_{h},\boldsymbol{w}_{[h]}]\\}=0.$

The above results also imply
$\boldsymbol{C}_{j,[j]}\boldsymbol{C}^{-1}_{[j]}\boldsymbol{C}_{[j],j}=\sum\limits_{s=i_{0}}^{i_{r}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{K}_{s}(s,j)$
and suggest an additive representation via orthogonal basis functions:

$\boldsymbol{w}_{j}=\sum\limits_{s=i_{0}}^{i_{r-1}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{e}_{j}$
(25)

Finally, considering the same sequence of nodes, recursively introduce the
covariance functions $\boldsymbol{F}_{0}(r,s)=\boldsymbol{C}_{r,s}$ and for
$j>1$,
$\boldsymbol{F}_{j}(r,s)=\boldsymbol{F}_{j-1}(r,s)-\boldsymbol{F}_{j\mbox{-}1}(r,j\mbox{-}1)\boldsymbol{F}_{j\mbox{-}1}^{-1}(j\mbox{-}1,j\mbox{-}1)\boldsymbol{F}_{j\mbox{-}1}(j\mbox{-}1,s)$.
We get

$\displaystyle\boldsymbol{F}_{j+1}(r,s)$
$\displaystyle=\boldsymbol{F}_{j}(r,s)-\boldsymbol{F}_{j}(r,j)\boldsymbol{F}_{j}^{-1}(j,j)\boldsymbol{F}_{j}(j,s)$
using (24)
$\displaystyle=\boldsymbol{F}_{j-1}(r,s)-\boldsymbol{F}_{j-1}(r,[j\mbox{-1:}j])\boldsymbol{F}_{j-1}^{-1}([j\mbox{-1:}j],[j\mbox{-1:}j])\boldsymbol{F}_{j-1}^{-1}([j\mbox{-1:}j],s)$
$\displaystyle=\boldsymbol{C}(r,s)-\boldsymbol{C}(r,[0\mbox{:}j])\boldsymbol{C}^{-1}([0\mbox{:}j],[0\mbox{:}j])\boldsymbol{C}([0\mbox{:}j],s)$
$\displaystyle=\boldsymbol{C}(r,s)-\boldsymbol{C}_{r,[j+1]}\boldsymbol{C}^{-1}_{[j+1]}\boldsymbol{C}_{[j+1],s}$

which can be iterated forward and results in an additional recursive way to
compute covariances in SpamTrees. Notice that while $\boldsymbol{K}_{j}$ is
formulated using the inverse of $J_{j}\times J_{j}$ matrix
$\boldsymbol{C}_{[j]}$, the $\boldsymbol{F}_{j}$’s require inversion of
smaller $n_{j}\times n_{j}$ matrices
$\boldsymbol{F}_{j-1}(j\mbox{-}1,j\mbox{-}1)$.

###  B.4 Properties of $\widetilde{\boldsymbol{C}}$

####  B.4.1 $\delta=1$

Choosing depth $\delta=1$ results in each node having exactly $1$ parent. In
this case the path ${\cal P}_{k\to
j}=\\{\boldsymbol{v}_{i_{1}},\dots,\boldsymbol{v}_{i_{r}}\\}$ from
$\boldsymbol{v}_{k}$ to $\boldsymbol{v}_{j}$, where
$\boldsymbol{v}_{i_{1}}=\boldsymbol{v}_{k}$,
$\boldsymbol{v}_{i_{r}}=\boldsymbol{v}_{j}$ and
$\\{\boldsymbol{v}_{i_{h}}\\}=\text{Pa}[{\boldsymbol{v}_{i_{h+1}}}]$, is
unique, and there is thus no distinction between shortest and longest paths:
${\cal P}_{k\to j}=\bar{{\cal P}}_{k\to j}=\widetilde{{\cal P}}_{k\to j}$.
Then denote $\dddot{\boldsymbol{H}}_{k\to
j}=\boldsymbol{H}_{i_{r}}\cdot\boldsymbol{H}_{i_{r-1}}\cdots\boldsymbol{H}_{i_{1}}$.
Let $\boldsymbol{v}_{z}$ be the concestor between $\boldsymbol{v}_{i}$ and
$\boldsymbol{v}_{j}$ i.e.
$\boldsymbol{v}_{z}=\text{con}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})=\arg\max_{\boldsymbol{v}_{k}\in\boldsymbol{V}}\\{k:{\cal
P}_{k\to i}\cap{\cal P}_{k\to j}\neq\emptyset\\}$ and the associated paths
${\cal P}_{z\to
i}=\\{\boldsymbol{v}_{i_{1}},\dots,\boldsymbol{v}_{i_{r_{i}}}\\}$ and ${\cal
P}_{z\to j}=\\{\boldsymbol{v}_{j_{1}},\dots,\boldsymbol{v}_{j_{r_{j}}}\\}$
where $\boldsymbol{v}_{i_{1}}=\boldsymbol{v}_{j_{1}}=\boldsymbol{v}_{z}$,
$\boldsymbol{v}_{i_{r_{i}}}=\boldsymbol{v}_{i}$ and
$\boldsymbol{v}_{j_{r_{j}}}=\boldsymbol{v}_{j}$. Then we can write
$\boldsymbol{w}_{i}=\boldsymbol{w}_{i_{r_{i}}}=\boldsymbol{H}_{i_{r_{i}}}\boldsymbol{w}_{i_{r_{i}-1}}+\boldsymbol{\nu}_{i_{r_{i}}}$
where $\boldsymbol{\nu}_{i_{r_{i}}}\sim
N(\mathbf{0},\boldsymbol{R}_{i_{r_{i}}})$ and proceed expanding
$\boldsymbol{w}_{i_{r_{i}-1}}$ to get
$\boldsymbol{w}_{i_{r_{i}}}=\boldsymbol{H}_{i_{r_{i}}}(\boldsymbol{H}_{i_{r_{i}-1}}\boldsymbol{w}_{i_{r_{i}-2}}+\boldsymbol{\nu}_{i_{r_{i}-1}})+\boldsymbol{\nu}_{i_{r_{i}}}=\boldsymbol{H}_{i_{r_{i}}}\boldsymbol{H}_{i_{r_{i}-1}}\boldsymbol{w}_{i_{r_{i}-2}}+(\boldsymbol{H}_{i_{r_{i}}}\boldsymbol{\nu}_{i_{r_{i}-1}}+\boldsymbol{\nu}_{i_{r_{i}}})$;
continuing downwardly along the tree we eventually find
$\boldsymbol{w}_{i}=\boldsymbol{H}_{i_{r_{i}}}\cdots\boldsymbol{H}_{i_{1}}\boldsymbol{w}_{i_{1}}+\widetilde{\boldsymbol{v}}_{i}=\dddot{\boldsymbol{H}}_{z\to
i}\boldsymbol{w}_{z}+\widetilde{\boldsymbol{v}}_{i}$ where
$\widetilde{\boldsymbol{v}}_{i}$ is independent of $\boldsymbol{w}_{z}$. After
proceeding analogously with $\boldsymbol{w}_{j}$, take
$\boldsymbol{\ell}_{i},\boldsymbol{\ell}_{j}$ such that
$\eta(\boldsymbol{\ell}_{i})=\boldsymbol{v}_{i}$ and
$\eta(\boldsymbol{\ell}_{j})=\boldsymbol{v}_{j}$. Then

$\text{Cov}_{\widetilde{p}}(w(\boldsymbol{\ell}_{i}),w(\boldsymbol{\ell}_{j}))=\dddot{\boldsymbol{H}}_{z\to
i}(\boldsymbol{\ell}_{i})\boldsymbol{C}_{z}\dddot{\boldsymbol{H}}_{z\to
j}(\boldsymbol{\ell}_{j})^{\top},$ (26)

where $\dddot{\boldsymbol{H}}_{z\to
i}(\boldsymbol{\ell}_{i})=\boldsymbol{C}(\boldsymbol{\ell}_{i},S_{i})\boldsymbol{C}^{-1}_{i}\dddot{\boldsymbol{H}}_{z\to{[i]}}$
and similarly for $\dddot{\boldsymbol{H}}_{z\to j}(\boldsymbol{\ell}_{j})$.

####  B.4.2 $1<\delta<M$

Take two nodes $\boldsymbol{v}_{i},\boldsymbol{v}_{j}\in\boldsymbol{V}$. If
$\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]\neq\emptyset$
then we apply the same logic as in B.4.3 using
$\boldsymbol{v}_{z}=\text{con}(\boldsymbol{v}_{i},\boldsymbol{v}_{j})$ as
root. If
$\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]=\emptyset$
and both nodes are at levels below $M_{\delta}$ then we use B.4.1. The
remaining scenario is thus one in which
$\boldsymbol{v}_{i}\in\boldsymbol{A}_{r}$, $r>M_{\delta}$ and
$\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]=\emptyset$.
We take $\boldsymbol{v}_{j}\in\boldsymbol{A}_{s}$, $s<M_{\delta}$ for
simplicity in exposition and without loss of generality. By (25)

$\displaystyle\boldsymbol{w}_{i}$
$\displaystyle=\sum\limits_{s=i_{M_{\delta}}}^{i_{r-1}}\boldsymbol{K}_{s}(i,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{e}_{i},$
(27)
$\displaystyle=\sum\limits_{s=i_{M_{\delta}+1}}^{i_{r-1}}\boldsymbol{K}_{s}(i,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{C}_{ix}\boldsymbol{C}^{-1}_{x}\boldsymbol{w}_{x}+\boldsymbol{e}_{i},$

where $\boldsymbol{v}_{x}\in\boldsymbol{A}_{M_{\delta}}$ is the parent node of
$\boldsymbol{v}_{i}$ at level $M_{\delta}$. The final result of (16) is then
achieved by noting that the relevant subgraph linking $\boldsymbol{v}_{x}$ and
$\boldsymbol{v}_{j}$ has depth $\delta_{x}=1$ and thus
$\text{Cov}(\boldsymbol{w}_{x},\boldsymbol{w}_{j})$ can be found via B.4.1,
then
$\text{Cov}(\boldsymbol{w}_{i},\boldsymbol{w}_{j})=\boldsymbol{C}_{ix}\boldsymbol{C}^{-1}_{x}\text{Cov}(\boldsymbol{w}_{x},\boldsymbol{w}_{j})=\boldsymbol{F}_{i}\text{Cov}(\boldsymbol{w}_{x},\boldsymbol{w}_{j})$.
Notice that $\boldsymbol{F}_{i}$ directly uses the directed edge
$\boldsymbol{v}_{x}\to\boldsymbol{v}_{i}$ in ${\cal G}$; for this reason the
path between $\boldsymbol{w}_{i}$ and
$\boldsymbol{w}_{z}=\text{con}(\boldsymbol{w}_{x},\boldsymbol{w}_{j})$ is the
actually the shortest path and we have
$\boldsymbol{v}_{z}\to\cdots\to\boldsymbol{v}_{x}\longrightarrow\boldsymbol{v}_{i}$.

####  B.4.3 $\delta=M$

Take $\boldsymbol{v}_{i},\boldsymbol{v}_{j}\in\boldsymbol{V}$ and the full
paths from the root $\widetilde{{\cal P}}_{0\to
i}=\\{i_{0},\dots,i_{r_{i}}\\}$ and $\widetilde{{\cal P}}_{0\to
j}=\\{j_{0},\dots,j_{r_{j}}\\}$, respectively. Then using (25) we have

$\displaystyle\boldsymbol{w}_{i}$
$\displaystyle=\sum\limits_{s\in\widetilde{{\cal P}}_{0\to
i}}\boldsymbol{K}_{s}(i,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{e}_{i}$
(28) $\displaystyle=\sum\limits_{s\in\widetilde{{\cal P}}_{0\to
i}\cap\widetilde{{\cal P}}_{0\to
j}}\boldsymbol{K}_{s}(i,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\sum\limits_{s\in\widetilde{{\cal
P}}_{0\to i}\setminus\widetilde{{\cal P}}_{0\to
j}}\boldsymbol{K}_{s}(i,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{e}_{i}$
$\displaystyle=\sum\limits_{s\in\widetilde{{\cal P}}_{0\to
i}\cap\widetilde{{\cal P}}_{0\to
j}}\boldsymbol{K}_{s}(i,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\widetilde{\boldsymbol{e}}_{i}$
$\displaystyle\boldsymbol{w}_{j}$
$\displaystyle=\sum\limits_{s\in\widetilde{{\cal P}}_{0\to
j}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{e}_{j}$
$\displaystyle=\sum\limits_{s\in\widetilde{{\cal P}}_{0\to
i}\cap\widetilde{{\cal P}}_{0\to
j}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\sum\limits_{s\in\widetilde{{\cal
P}}_{0\to j}\setminus\widetilde{{\cal P}}_{0\to
i}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{e}_{j}$
$\displaystyle=\sum\limits_{s\in\widetilde{{\cal P}}_{0\to
i}\cap\widetilde{{\cal P}}_{0\to
j}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\widetilde{\boldsymbol{e}}_{j},$

where
$\text{Cov}(\widetilde{\boldsymbol{e}}_{i},\widetilde{\boldsymbol{e}}_{j})=0$.
Then since $\boldsymbol{e}_{s}$ are independent and $\boldsymbol{e}_{s}\sim
N(\mathbf{0},\boldsymbol{K}_{s}(s,s))$ we find

$\displaystyle\text{Cov}_{\widetilde{p}}(\boldsymbol{w}_{i},\boldsymbol{w}_{j})$
$\displaystyle=\text{Cov}\left(\sum\limits_{s\in\widetilde{{\cal P}}_{0\to
i}\cap\widetilde{{\cal P}}_{0\to
j}}\boldsymbol{K}_{s}(i,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{e}_{i},\right.$
(29) $\displaystyle\left.\qquad\qquad\qquad\sum\limits_{s\in\widetilde{{\cal
P}}_{0\to i}\cap\widetilde{{\cal P}}_{0\to
j}}\boldsymbol{K}_{s}(j,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{e}_{s}+\boldsymbol{e}_{j}\right)$
$\displaystyle=\sum\limits_{s\in\widetilde{{\cal P}}_{0\to
i}\cap\widetilde{{\cal P}}_{0\to
j}}\boldsymbol{K}_{s}(i,s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{K}_{s}(s,j)+\boldsymbol{1}_{i=j}\\{\boldsymbol{K}_{i}(i,i)\\}.$

We conclude by noting that $\delta=M$ implies $\widetilde{{\cal P}}_{0\to
i}\cap\widetilde{{\cal P}}_{0\to
j}=\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]$;
considering two locations $\boldsymbol{\ell}_{i},\boldsymbol{\ell}_{j}\in{\cal
D}^{*}$ such that $\eta(\boldsymbol{\ell}_{i})=\boldsymbol{v}_{i}$ and
$\eta(\boldsymbol{\ell}_{j})=\boldsymbol{v}_{j}$ we obtain

$\displaystyle\text{Cov}_{\widetilde{p}}(w(\boldsymbol{\ell}_{i}),w(\boldsymbol{\ell}_{j}))$
$\displaystyle=\sum\limits_{s:\\{\boldsymbol{v}_{s}\in\text{Pa}[{\boldsymbol{v}_{i}}]\cap\text{Pa}[{\boldsymbol{v}_{j}}]\\}}\boldsymbol{K}_{s}(\boldsymbol{\ell}_{i},s)\boldsymbol{K}_{s}^{-1}(s,s)\boldsymbol{K}_{s}(s,\boldsymbol{\ell}_{j})+\boldsymbol{1}_{\boldsymbol{\ell}_{i}=\boldsymbol{\ell}_{j}}\\{\boldsymbol{K}_{i}(\boldsymbol{\ell}_{i},\boldsymbol{\ell}_{j})\\}.$
(30)

###  B.5 Computational cost

We make some assumptions here to simplify the calculation of overall cost:
first, we assume that reference locations are all observed ${\cal
S}\subset{\cal T}$, and consequently ${\cal U}={\cal T}\setminus{\cal S}$.
Second, we assume that all reference subsets have the same size i.e.
$|S_{i}|=N_{s}$ for all $i$. Third, we assume all nodes have the same number
of children at the next level in ${\cal G}$, i.e. if
$\boldsymbol{v}_{i}\in\boldsymbol{A}_{r}$ with $r<M-1$, then
$|\text{Ch}[{\boldsymbol{v}_{i}}]\cap\boldsymbol{A}_{r+1}|=C$, whereas if
$r=M-1$ then $|\text{Ch}[{\boldsymbol{v}_{i}}]|=N_{u}$. Fourth, we assume that
all non-reference subsets are singletons i.e. if
$\boldsymbol{v}_{i}\in\boldsymbol{B}$ then $|U_{i}|=1$. The latter two
assumptions imply (7). We also fix $CN_{s}=N_{u}$. As a result, the number of
nodes at level $r=0,\dots,M-1$ is $C^{r}$, therefore
$|\boldsymbol{A}|+|\boldsymbol{B}|=\sum_{r=0}^{M-1}C^{r}+N_{u}C^{M-1}=\frac{C^{M}-1}{C-1}+N_{s}C^{M}$.
Then the sample size is $n=|{\cal T}|=|{\cal S}|+|{\cal
U}|=N_{s}\frac{C^{M+1}-1}{C-1}$ hence $M\approx\log_{C}(n/N_{s})$. Starting
with $\delta=M$, the parent set sizes $J_{i}$ for a node
$\boldsymbol{v}_{i}\in\boldsymbol{A}_{r}$ grow with $r$ as $J_{i}=rN_{s}$ and
if $\boldsymbol{v}_{i}\in\boldsymbol{B}$ then $J_{i}=MN_{s}$. The cost of
computing $p(\boldsymbol{w}\,|\,\boldsymbol{\theta})$ is driven by the
calculation of $\boldsymbol{H}_{j}$, which is $O(r^{2}N_{s}^{3})$ for
reference nodes at level $r$, for a total of
$O(N_{s}^{3}\sum_{r=0}^{M-1}C^{r}r^{2})$. Since for common choices of $C$ and
$M$ we have
$\sum_{r=0}^{M-1}C^{r}r^{2}N_{s}^{3}\leq\sum_{r=0}^{M-1}C^{2r}N_{s}^{3}=\frac{C^{2M}-1}{C^{2}-1}N_{s}^{3}\approx
C^{M}N_{s}^{3}\approx\frac{n}{N_{s}}N_{s}^{3}=nN_{s}^{2}$ then the cost for
reference sets is $O(nN_{s}^{2})$. Analogously for non reference nodes we get
$O(C^{M}M^{2}N_{s}^{3})$ which leads to a cost of $O(nN_{s}^{2})$. The cost of
sampling $\boldsymbol{w}$ is mainly driven by the computation of the Cholesky
factor of a $N_{s}\times N_{s}$ matrix at each of $\frac{C^{M}-1}{C-1}$
reference nodes, which amounts to $O(nN_{s}^{2})$. For the $N_{s}C^{M}$ non-
reference nodes the main cost is in computing
$\boldsymbol{H}_{i}\boldsymbol{w}_{[i]}$ which is $M^{2}N_{s}^{2}$ for overall
cost $O(C^{M}M^{2}N_{s}^{3})$ which again is $O(nN_{s}^{2})$. Obtaining
$\boldsymbol{F}_{i}^{(c)}$ at the root of ${\cal G}$ is associated to a cost
$O(N_{s}^{2}\frac{C^{M}-C}{C-1})$ which is $O(nN_{s})$ but constitutes a
bottleneck if such operation is performed simultaneously to sampling; however
this bottleneck is eliminated in Algorithm 3.

If $\delta=1$ then the parent set sizes $J_{i}$ for all nodes
$\boldsymbol{v}_{i}\in\boldsymbol{V}$ are constant $J_{i}=N_{s}$; since the
nodes at levels $0$ to $M-1$ have $C$ children, the asymptotic cost of
computing $p(\boldsymbol{w}\,|\,\boldsymbol{\theta})$ is