# The Computational Curse of Big Data for Bayesian Additive Regression Trees:
A Hitting Time Analysis

Yan Shuo Tan  Department of Statistics and Data Science, National University
of Singapore  Omer Ronen Equal contribution, alphabetical ordering Department
of Statistics, UC Berkeley  Theo Saarinen††footnotemark:  Department of
Statistics, UC Berkeley  Bin Yu  Department of Statistics, UC Berkeley
Department of Electrical Engineering and Computer Sciences, UC Berkeley
Center for Computational Biology, UC Berkeley

###### Abstract

Bayesian Additive Regression Trees (BART) is a popular Bayesian non-parametric
regression model that is commonly used in causal inference and beyond. Its
strong predictive performance is supported by theoretical guarantees that its
posterior distribution concentrates around the true regression function at
optimal rates under various data generative settings and for appropriate prior
choices. In this paper, we show that the BART sampler often converges slowly,
confirming empirical observations by other researchers. Assuming discrete
covariates, we show that, while the BART posterior concentrates on a set
comprising all optimal tree structures (smallest bias and complexity), the
Markov chain’s hitting time for this set increases with $n$ (training sample
size), under several common data generative settings. As $n$ increases, the
approximate BART posterior thus becomes increasingly different from the exact
posterior (for the same number of MCMC samples), contrasting with earlier
concentration results on the exact posterior. This contrast is highlighted by
our simulations showing worsening frequentist undercoverage for approximate
posterior intervals and a growing ratio between the MSE of the approximate
posterior and that obtainable by artificially improving convergence via
averaging multiple sampler chains. Finally, based on our theoretical insights,
possibilities are discussed to improve the BART sampler convergence
performance.

##  1 Introduction

###  1.1 The Rise of BART

Decision tree models such as CART (Breiman et al.,, 1984) are piecewise
constant regression models obtained by recursively partitioning the covariate
space along coordinate axes. They and their ensembles such as Random Forests
(RFs) (Breiman,, 2001) and Gradient Boosted Trees (GBTs) (Friedman,, 2001;
Chen and Guestrin,, 2016) have proved to be enormously successful because of
their strong predictive performance (Caruana and Niculescu-Mizil,, 2006;
Caruana et al.,, 2008; Fernández-Delgado et al.,, 2014). Indeed, RFs and GBTs
regularly outperform even deep learning on medium-sized tabular datasets
(Grinsztajn et al.,, 2022). Nonetheless, these tree-based methods still suffer
from several notable problems: They are defined via algorithms rather than via
statistical models, so it is often difficult to quantify the uncertainty of
their predictions; they use greedy splitting criteria, so there is no
guarantee for the optimality of the fitted model; RFs in particular grow their
trees independently of each other, therefore making them statistically
inefficient when fitted to data with additive structure (Tan et al., 2022a, ).

To address these issues, Chipman et al., (1998) proposed a Bayesian adaptation
of CART (BCART) and later an ensemble of Bayesian CART trees, which they
called Bayesian Additive Regression Trees (BART) (Chipman et al.,, 2010).
These are Bayesian non-parametric regression models, which put a prior on the
space of regression functions, assume a likelihood for the observed data, and
combine these to obtain a posterior. In the case of Bayesian CART and BART,
priors and posteriors are supported on the subspace of functions that can be
realized by decision trees (or their ensembles). Similar to Gaussian process
(GP) regression, the posterior distribution can be used to provide posterior
predictive credible intervals. On the other hand, unlike GP regression, there
is no closed form formula for the BART posterior and one has to sample from it
approximately via a Markov chain Monte Carlo (MCMC) algorithm.

BART has been shown empirically to enjoy strong predictive performance that is
sometimes even superior to that of RFs and GBTs, especially after
hyperparameter optimization (Hill et al.,, 2020). Naturally, it has become
increasingly popular in diverse fields ranging from the social sciences (Green
and Kern,, 2010; Yeager et al.,, 2019) to biostatistics (Wendling et al.,,
2018; Starling et al.,, 2020) and has been particularly enthused causal
inference researchers (Hill,, 2011; Green and Kern,, 2012; Kern et al.,, 2016;
Dorie et al.,, 2019; Hahn et al.,, 2019). Extending and improving BART
methodology remains a highly active area of research, with many variants of
the algorithm proposed over the last few years (see for instance Linero,
(2018); Pratola, (2016); Pratola et al., (2020); Luo and Pratola, (2023), as
well as the survey Hill et al., (2020) and the references therein.)

The strong predictive performance of BART is supported by a burgeoning body of
theoretical evidence regarding the BART posterior. Most significantly,
researchers have shown that the BART posterior concentrates around the true
regression function used to generate the response data as $n$, the number of
training samples, increases, with this concentration happening at optimal
rates under various assumptions on the smoothness and sparsity of the
regression function and for appropriate prior choices (Ročková and Saha,,
2019; Ročková and van der Pas,, 2020; Linero and Yang,, 2018; Jeong and
Rockova,, 2020; Rockova and Rousseau,, 2021; Castillo and Ročková,, 2021).
Achieving these optimal rates does not require any oracle knowledge or
hyperparameter tuning—instead BART automatically adapts to the level of
smoothness and sparsity, with the former even happening at a local level
(Rockova and Rousseau,, 2021).

###  1.2 Observed Poor Mixing of BART and its Significance

While there is evidence that the BART _posterior_ enjoys favorable properties
under a variety of settings, the fact that we can only sample approximately
from the posterior via MCMC creates a gap in our understanding of how and why
BART works. Specifically, if the sampler chain does not converge efficiently
to the posterior distribution, that is, if it does not mix well, the output of
BART _algorithm_ may not enjoy the same desirable inferential properties as
the BART posterior. Most popular BART implementations use remarkably similar
samplers based on the original design of Chipman et al., (2010), which uses a
Bayesian backfitting approach to update one tree at a time via proposed local
changes to the tree nodes coupled with a Metropolis-Hastings filter.
Unfortunately, as described by Hill et al., (2020), “while this algorithm is
often effective, it does not always mix well.” Indeed, poor mixing for this
sampler has been empirically documented by multiple sources (Chipman et al.,,
1998; Carnegie,, 2019).

The literature contains various suggestions on how to improve the mixing time
for the BART sampler. These include parallelization (Pratola et al.,, 2014),
modifying the MCMC proposal moves (Wu et al.,, 2007; Pratola,, 2016; Kim and
Rockova,, 2023), warm starts from greedily constructed tree ensembles (He and
Hahn,, 2021), or running multiple chains (Carnegie,, 2019). Despite this
interest, there has been minimal theoretical work done to quantify the mixing
time and to understand why and under what settings slow mixing occurs.

###  1.3 Main Contributions

In this paper, we show theoretically, assuming discrete covariates, that the
BART sampler often converges slowly to its posterior, confirming the empirical
observations of Hill et al., (2020) and other researchers. In fact, the
convergence unexpectedly becomes _worse_ as $n$ increases, in contrast with
the posterior’s concentration to the true regression function becoming better.

To state our results more formally, note that a regression tree is
parameterized via its _tree structure_ (which features are split on and at
which thresholds) and its leaf parameters (the function value on each leaf).
BART combines these parameters over multiple trees to parameterize a tree
ensemble (see Section 2.2.). Since the leaf parameters can be sampled in
closed form conditionally on tree structures, the original BART sampler, under
a slight modification, can be thought of as a Markov chain on the space of
tree structures.

We show that the BART posterior concentrates on a set comprising all optimal
tree structures (i.e. those with the smallest bias and complexity), which also
forms a highest posterior density region (HPDR). On the other hand, the BART
sampler’s _hitting time_ for this set increases with $n$ under four common
data generative settings. In other words, the sampler requires more and more
steps to reach any optimal tree structure. Note that this is a frequentist
analysis and requires the assumption of a generative model for the data that
can and will be different from the Bayesian parameterization. Our hitting time
lower bounds are summarized in the Table 1.

Data Generating Process | Allowed Moveset | Multiple Trees | Lower Bound
---|---|---|---
Additive | Full | Yes | Square root
Additive | “Grow” and “Prune” | Yes | Polynomial
Contains Pure Interaction | “Grow” and “Prune” | Yes | Square root
Root Dependence | “Grow” and “Prune” | No | Exponential
Table 1: A summary of the HPDR hitting time lower bounds provided by our paper
and their dependence on the training sample size $n$ (last column). The first
two lower bounds apply to additive generative models. The third applies to a
setting where the generative regression function contains a pure interaction
(defined in Section 5.2.) The fourth applies to the setting where the
generative function has “root dependence”, which represents a form of
asymmetric dependence on the features (defined in Section 5.3.) All lower
bounds apart from the fourth allow BART to use multiple trees.

We complement our theoretical analysis with a comprehensive simulation study
of BART involving a wide range of data-generating processes with continuous
covariates. From now on, we use the term approximate posterior to refer to the
distribution obtained from 1000 MCMC samples after a generous burn-in (5000
iterations as opposed to the default of 100). To create a proxy for the exact
posterior, we combine samples from multiple (5) sampler chains, which is known
to improve mixing (Carnegie,, 2019).111Carnegie, (2019) investigated the use
of 1, 4, and 10 chains. Taking reference from this, we used 5 chains but did
not investigate the effect of increasing the number of chains beyond this
value. We believe that the optimal number could possibly vary with the sample
size and DGP. We compare the performance of the approximate and multi-chain
approximate posteriors via two metrics: (i) the RMSE of their posterior mean
function from the true regression function on a held-out test set, (ii) the
empirical coverage of their pointwise credible intervals for the true
regression function. In both cases, the relative performance gap between the
approximate and multi-chain approximate posteriors (measured as a ratio)
increases as $n$ increases. These two findings provide further evidence that
the approximate BART posterior becomes increasingly and meaningfully different
from the exact posterior as $n$ increases. We also perform two other
experiments to validate this claim in the setting of our theoretical lower
bounds. Our theory and simulations thus echo existing advice that BART users
should not blindly take its posterior credible intervals at face value and
should run multiple chains whenever feasible.

Lastly, our theoretical results and their proof strategies yield insights on
why the BART sampler has trouble converging, which leads us to suggest
possible ways to improve its performance. Most importantly, our proof will
show that a major reason why hitting times grow with training sample size is
because the “temperature” of the BART sampler is inversely proportional to the
training sample size. There is no need for temperature and sample size to be
intrinsically tied in this manner, and we conjecture that building in more
flexible temperature control, such as via simulated tempering, may help to
accelerate mixing.

###  1.4 Prior Work on Mixing for BART

Kim and Rockova, (2023) and our prior work (Ronen et al.,, 2022) sought to
analyze mixing times for the Bayesian CART sampler. Both made the surprising
discovery that its mixing time can grow exponentially in the training sample
size, when the only allowed moves are growing new leaves and pruning existing
ones. Kim and Rockova, (2023) studied this in a one-dimensional setting and
further showed that, with a more aggressive move set, Bayesian CART
constrained to dyadic splits has a mixing time upper bound that is linear in
the sample size. Ronen et al., (2022)’s proof strategy was to show that this
Markov chain has a bottleneck state—the trivial tree comprising a single node.
This bottleneck arises because when Bayesian CART makes a wrong first split
followed by other informative splits, the only way to reverse the wrong first
split involves pruning the informative splits, which becomes increasingly
difficult as the training sample size increases.

When we submitted our prior work to a peer review, some reviewers rightly
pointed out that our mixing time lower bounds may be misleading to
practitioners. This is because different tree structures could realize the
same partition of the covariate space and hence implement the same regression
function. In other words, the Bayesian CART model is not identifiable at the
level of tree structures, which means that failure to mix in this space of
tree structures may be benign (Redner and Walker,, 1984). It does not reflect
a failure to mix at the level of regression functions, let alone any
degradation in the inferential properties of the algorithm’s output.

Figure 1: The Bayesian CART model is not identifiable at the level of tree
structures. The two tree structures shown on the left both realize the same
partition of the covariate space and even the same regression function, which
is shown on the right.

A further limitation of our prior analysis is that mixing time is defined as
the worst case time to convergence as we vary over all possible initial
states. Since Bayesian CART is always initialized at the trivial tree, a lower
bound on this worst case quantity may be unreasonably pessimistic.

These issues with relevance are resolved by considering hitting times.

##  2 Data Generation Models for BART and for Frequentist Analysis

###  2.1 Generative Model for Frequentist Analysis

We assume a $d$-dimensional discrete covariate space
$\mathcal{X}=\\{1,2,\ldots,b\\}^{d}$ and a regression function
$f^{*}\colon\mathcal{X}\to\mathbb{R}$. Suppose that we observe a training data
set $\mathcal{D}_{n}$ comprising $n$ independent and identically distributed
tuples $(\mathbf{x}_{i},y_{i})$ with $\mathbf{x}_{i}\sim\nu$ and
$y_{i}=f^{*}(\mathbf{x}_{i})+\epsilon_{i}$ for $i=1,2\ldots,n$. Here, $\nu$ is
a measure on $\mathcal{X}$ with full support, while $\epsilon$ is a sub-
Gaussian random variable. For notational convenience, we will use $\mathbf{X}$
to denote the $n\times d$ matrix formed by stacking the covariate vectors as
rows. We will also use $\mathbf{y}$, $\mathbf{f}^{*}$, and
$\boldsymbol{\epsilon}$ to denote the $n$-dimensional vectors formed by
stacking the responses $y_{i}$, the function values $f^{*}(\mathbf{x}_{i})$
and the noise components $\epsilon_{i}$ respectively. Similarly, we will
denote all vectors and matrices with boldface notation, with subscripts
referencing the index of the vector. Vector coordinates will be denoted using
regular font. We will denote probabilities and expectations with respect to
$\mathcal{D}_{n}$ using $\mathbb{P}_{n}$ and $\mathbb{E}_{n}$ respectively.

###  2.2 Bayesian Model Specification for BART

We first describe the version of BART that we analyze in our paper, before
discussing its differences with the version described by Chipman et al.,
(2010) and which is still most commonly used in practice.

#### Regression trees.

A binary axis-aligned regression tree is parameterized by a tuple
$(\mathfrak{T},\boldsymbol{\mu})$. Here, $\mathfrak{T}$ refers to the _tree
structure_ , which specifies the topology of the tree as a rooted binary tree
planar graph and, given an ordering of the graph’s vertices (e.g. via breadth-
first search), specifies the splitting rule for each internal node $j$. Note
that the splitting rule comprises a feature $v_{j}$ and a threshold $t_{j}$.
The leaves $\mathcal{L}_{1},\mathcal{L}_{2},\ldots,\mathcal{L}_{b}$ of
$\mathfrak{T}$ thus correspond to rectangular regions of the covariate space
that together form a partition of the space. We let
$\boldsymbol{\mu}\in\mathbb{R}^{b}$ be a vector of _leaf parameters_ , one for
each leaf of $\mathfrak{T}$. Together, $(\mathfrak{T},\boldsymbol{\mu})$
specify a piecewise constant function $g$ that outputs

$g(\mathbf{x};\mathfrak{T},\boldsymbol{\mu})=\mu_{l(\mathbf{x})},$

where $l(\mathbf{x})$ is the index of the leaf containing $\mathbf{x}$.

#### Sum-of-trees model.

Given observed data $\mathcal{D}_{n}$, the BART model posits
$y_{i}=f(\mathbf{x}_{i})+e_{i}$ for $i=1,2\ldots,n$, where
$e_{1},e_{2},\ldots,e_{n}\sim_{i.i.d.}\mathcal{N}(0,\sigma^{2})$, and $f$ is a
sum of the outputs of $m$ trees:

$f(\mathbf{x})=g(\mathbf{x};\mathfrak{T}_{1},\boldsymbol{\mu}_{1})+g(\mathbf{x};\mathfrak{T}_{2},\boldsymbol{\mu}_{2})+\cdots+g(\mathbf{x};\mathfrak{T}_{m},\boldsymbol{\mu}_{m}).$

We denote the ordered tuple
$(\mathfrak{T}_{1},\mathfrak{T}_{2},\ldots,\mathfrak{T}_{m})$ by
$\mathfrak{E}$ and call it a _tree structure ensemble_ (TSE). We shall abuse
notation and use $\boldsymbol{\mu}$ to refer to the concatenation of
$\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2},\ldots,\boldsymbol{\mu}_{m}$. Note
that, when conditioned on $\mathfrak{E}$ and $\mathbf{X}$, this is just a
Bayesian linear regression model. To see this, let $\boldsymbol{\Psi}$ denote
the $n\times b$ matrix whose columns are the indicator vectors over the
training set of each leaf in $\mathfrak{E}$. We then have

$\mathbf{y}=\boldsymbol{\Psi}\boldsymbol{\mu}+\mathbf{e}.$ (1)

#### Priors.

We assume a fixed prior distribution $p$ on $\Omega_{\text{TSE},m}$, the space
of TSEs with $m$ trees.222Because the emphasis in our theory is on the
dependence of hitting times on training sample size in large samples, the
specific form of the prior holds no bearing on our results. Conditioned on a
TSE $\mathfrak{E}$, the conditional prior distribution on the leaf parameters
is an isotropic Gaussian, i.e.
$p(\boldsymbol{\mu}|\mathfrak{E})\sim\mathcal{N}(0,(\sigma^{2}/\lambda)\mathbf{I}_{b})$,
where $b$ is the total number of leaves in all trees in $\mathfrak{E}$. Both
$\sigma^{2}$ and $\lambda$ are assumed to be fixed hyperparameters, with
$\sigma^{2}$ taking the same value as that used in the variance of the
additive noise $e_{i}$, while $\lambda$ is a modulation parameter that should
be set to approximately the reciprocal of the signal-to-noise ratio.

#### Differences with in-practice BART.

Chipman et al., (1998) proposed a prior on tree structures defined in terms of
a stochastic process. Starting from a single root node, the process
recursively splits each node at depth $d$ with probability
$\alpha(1+d)^{-\beta}$, where $\alpha$ and $\beta$ are hyperparameters with
default values $\alpha=0.95$ and $\beta=2$ respectively. Features and
thresholds for splits are selected uniformly at random. Chipman et al., (2010)
extended this to a prior on TSEs by independence, i.e.
$p(\mathfrak{E})=\prod_{j=1}^{m}p(\mathfrak{T}_{j})$. After rescaling the
response variable to lie between $-0.5$ and $0.5$, the leaf parameter standard
deviation is set to be $\sigma_{\mu}=0.5/k\sqrt{m}$, where $k$ is a further
hyperparameter with default value $k=2$. Finally, an inverse-$\chi^{2}$
hyperprior is placed on the noise variance $\sigma^{2}$ and is calibrated to
the observed data. This last assumption is the only way in which the BART
model we study in this paper departs from that in Chipman et al., (2010). We
make this change for analytical tractability and believe it to be minor, since
simulations show that the posterior on $\sigma^{2}$ quickly converges to a
fixed value and our theoretical guarantees hold for any fixed choice of
$\sigma^{2}$.

##  3 Sampling from BART via MCMC

###  3.1 The In-Practice BART Sampler

The sampler proposed by Chipman et al., (2010) can be described as a
“Metropolis-within-Gibbs MCMC sampler” (Hill et al.,, 2020). More precisely,
for each outer loop of the algorithm, it iterates over the tree indices
$j=1,2,\ldots,m$ and updates the $j$-th pair
$(\mathfrak{T}_{j},\boldsymbol{\mu}_{j})$ using an approximate draw from the
conditional distribution
$p(\mathfrak{T}_{j},\boldsymbol{\mu}_{j}|\mathfrak{E}_{-j},\boldsymbol{\mu}_{-j},\mathbf{y},\sigma^{2})$,
where $\mathfrak{E}_{-j}$ and $\boldsymbol{\mu}_{-j}$ refer to the
concantenation of current tree structures and leaf parameter vectors
respectively, each with the $j$-th index omitted.333Note that the conditional
distribution is of course also conditional on the observed covariate data
$\mathbf{X}$. However, since this is always conditioned upon, we omit it from
our notation to avoid clutter. As a final step in the loop, it updates
$\sigma^{2}$ using a draw from its full conditional distribution
$p(\sigma^{2}|\mathfrak{E},\boldsymbol{\mu},\mathbf{y})$.

To describe how to sample (approximately) from
$p(\mathfrak{T}_{j},\boldsymbol{\mu}_{j}|\mathfrak{E}_{-j},\boldsymbol{\mu}_{-j},\mathbf{y},\sigma^{2})$,
we first describe Chipman et al., (1998)’s algorithm for Bayesian CART, i.e.
when the ensemble comprises a single tree. In this case, we first factorize
the posterior into a conditional posterior on leaf parameters and a marginal
posterior on tree structures:444Chipman et al., (1998)’s formulation of the
Bayesian CART sampler marginalizes out $\sigma^{2}$ instead of conditioning on
it. As this is no longer done for BART, we omit discussing it to avoid
confusing readers.

$p(\mathfrak{T},\boldsymbol{\mu}|\mathbf{y},\sigma^{2})=p(\boldsymbol{\mu}|\mathfrak{T},\mathbf{y},\sigma^{2})p(\mathfrak{T}|\mathbf{y},\sigma^{2}).$
(2)

The first multiplicand on the right,
$p(\boldsymbol{\mu}|\mathfrak{T},\mathbf{y},\sigma^{2})$, is a multivariate
Gaussian (with diagonal covariance) and can be sampled from directly. The
second multiplicand is proportional to
$p(\mathbf{y}|\mathfrak{T},\sigma^{2})p(\mathfrak{T})$, which is the product
between a marginal likelihood of a Bayesian linear regression model and the
prior. The marginal likelihood can be computed using standard techniques, but
the posterior cannot be sampled directly. As such, a Metropolis-Hastings
sampler is used with the following types of proposed moves:

  1. 1.

Pick a leaf in the tree and split it (grow);

  2. 2.

Pick two adjacent leaves and collapse them back into a single leaf (prune);

  3. 3.

Pick an interior node and change the splitting rule (change);

  4. 4.

Pick a pair of parent-child nodes that are both internal and swap their
splitting rules, unless both children of the parent have the same splitting
rules, in which case, swap the splitting rule of the parent with that of both
children (swap).

Note that all selections in these proposed moves (of nodes, splitting rules,
etc.) are made uniformly at random from all available choices.555Splits that
result in empty leaves are not allowed. The proposed move types are chosen
with probabilities $\pi_{g}$, $\pi_{p}$, $\pi_{c}$, and $\pi_{s}$
respectively. Let $Q(-,-)$ denote the transition kernel of the proposal, i.e.
$Q(\mathfrak{T},\mathfrak{T}^{*})$ is the probability of tree structure
$\mathfrak{T}^{*}$ being proposed given current tree structure $\mathfrak{T}$.
With $\mathfrak{T}$ and $\mathfrak{T}^{*}$ thus defined, the Metropolis-
Hastings algorithm accepts the proposal with probability

$\alpha(\mathfrak{T},\mathfrak{T}^{*})\coloneqq\min\left\\{\frac{Q(\mathfrak{T}^{*},\mathfrak{T})p(\mathfrak{T}^{*}|\mathbf{y},\sigma^{2})}{Q(\mathfrak{T},\mathfrak{T}^{*})p(\mathfrak{T}|\mathbf{y},\sigma^{2})},1\right\\}.$

A simple reparameterization trick is used to adapt this sampler to the case
when the ensemble has multiple trees. Because of the independence of the
priors on different trees and the Gaussian likelihood, the conditional
posterior
$p(\mathfrak{T}_{j},\boldsymbol{\mu}_{j}|\mathfrak{E}_{-j},\boldsymbol{\mu}_{-j},\sigma^{2},\mathbf{y})$
can be rewritten in terms of the residual vector

$\mathbf{r}_{-j}\coloneqq\mathbf{y}-\sum_{k\neq
j}g(\mathbf{X};\mathfrak{T}_{k},\boldsymbol{\mu}_{k}).$

Specifically, we have

$p(\mathfrak{T}_{j},\boldsymbol{\mu}_{j}|\mathfrak{E}_{-j},\boldsymbol{\mu}_{-j},\sigma^{2},\mathbf{y})=p(\mathfrak{T}_{j},\boldsymbol{\mu}_{j}|\mathbf{r}_{-j},\sigma^{2}),$

where the right-hand side is the single-tree posterior. A single Metropolis-
Hastings update step as described above is performed to draw an approximate
sample $(\mathfrak{T}_{j},\boldsymbol{\mu}_{j})$ from the conditional
posterior.

###  3.2 The Analyzed BART Sampler

The BART sampler described above is difficult to analyze because the
deterministic Gibbs outer loop makes it a time-varying Markov chain. More
significantly, it is convenient in analyzing Bayesian CART to collapse the
Markov chain state space by marginalizing out the leaf parameters. The
collapsed state space is simply the space of tree structures, which is
discrete and finite. However, we are unable to do this for BART in general
because of the conditioning on the residuals from other trees in the inner
loop. Both of these difficulties make it impossible to apply standard
techniques in Markov chain theory.

To overcome this, we propose an adaptation of the sampler that brings it
closer to Bayesian CART. First, we imitate (2) and factorize the posterior
into a conditional posterior on leaf parameters and a marginal posterior on
tree _ensemble_ structures:

$p(\mathfrak{E},\boldsymbol{\mu}|\mathbf{y})=p(\boldsymbol{\mu}|\mathfrak{E},\mathbf{y})p(\mathfrak{E}|\mathbf{y}).$

The conditional posterior on leaf parameters is still a multivariate Gaussian
and can be sampled from directly, while the marginal posterior on tree
ensemble structures remains proportional to the product of the marginal
likelihood of a Bayesian linear regression model and the prior:
$p(\mathfrak{E}|\mathbf{y})\propto p(\mathbf{y}|\mathfrak{E})p(\mathfrak{E})$
(see (1).) To sample from this marginal posterior, we run Metropolis-Hastings
MCMC similarly to before. However, instead of cycling deterministically
through the trees in an inner loop as before, we pick a tree index uniformly
at random. We propose an updated tree using the same transition kernel
$Q(-,-)$, but write the acceptance probability in terms of the full marginal
posterior instead of conditioning on the residuals from other trees. For
further clarity, the algorithm is summarized in pseudocode as Algorithm 1.

We denote the transition kernel of the sampler using $P(-,-)$, and, to avoid
confusion with randomness arising from sampling the training set, we will
denote all probabilities and expectations with respect to the algorithmic
randomness using $P$ and $E$ respectively.

Algorithm 1 BART sampler.

1:BART($\mathcal{D}_{n}$: data, $m$: no. of trees, $\sigma^{2}$: guess for
noise variance, $\lambda$: guess for reciprocal SNR, $\boldsymbol{\pi}$:
proposal probabilities, $p_{TSE}$: TSE prior, $t_{max}$: no. of sampler
iterations)

2:Initialize $\mathfrak{T}_{1},\mathfrak{T}_{2},\ldots,\mathfrak{T}_{m}$ as
trivial trees.

3:for $t=1,2,\ldots,t_{max}$

4: Sample $k\sim\text{Unif}(\\{1,2,\ldots,m\\})$.

5: Propose $\mathfrak{T}^{*}\sim Q(\mathfrak{T}_{k},\mathfrak{T}^{*})$.

6: Set
$\alpha(\mathfrak{T}_{k},\mathfrak{T}^{*})=\min\left\\{\frac{Q(\mathfrak{T}^{*},\mathfrak{T}_{k})p(\mathfrak{T}_{1},\ldots,\mathfrak{T}_{k-1},\mathfrak{T}^{*},\mathfrak{T}_{k+1},\ldots,\mathfrak{T}_{m}|\mathbf{y})}{Q(\mathfrak{T}_{k},\mathfrak{T}^{*})p(\mathfrak{T}_{1},\ldots,\mathfrak{T}_{k-1},\mathfrak{T}_{k},\mathfrak{T}_{k+1},\ldots,\mathfrak{T}_{m}|\mathbf{y})},1\right\\}$.

7: Set $\mathfrak{T}_{k}=\mathfrak{T}^{*}$ with probability
$\alpha(\mathfrak{T},\mathfrak{T}^{*})$.

##  4 BIC for BART and Posterior Concentration on Optimal TSEs

The goal of this section is to first show how to quantify the bias and
complexity of a TSE separately and then jointly via BIC. We will then show
that, as a function of TSEs, the posterior probability
$p(\mathfrak{E}|\mathbf{y})$ concentrates on the set of TSEs with zero bias
and the lowest possible complexity, and are therefore minimizers of BIC. As
such, as argued in the introduction, the highest posterior density region
contains all of the most desirable TSEs. This implies that lower bounds on the
hitting times of this region reflect computational drawbacks of practical
consequence.

###  4.1 Measuring Bias and Complexity for TSEs

We first discuss how to quantify the bias and complexity of a TSE.

#### Partitions.

A _cell_ $\mathcal{C}$ is a rectangular region of $\mathcal{X}$, i.e.

$\mathcal{C}=\left\\{\mathbf{x}\in\mathcal{X}~{}\colon~{}a_{i}<x_{i}\leq
b_{i}~{}\text{for}~{}i=1,\ldots,d\right\\},$

with lower and upper limits $a_{i}$ and $b_{i}$ respectively in coordinate $i$
for $i=1,2,\ldots,d$. A _partition_ is a collection of disjoint cells
$\mathcal{C}_{1},\ldots,\mathcal{C}_{b}$ whose union is the whole space
$\mathcal{X}$. Every tree structure $\mathfrak{T}$ induces a partition
$\mathcal{P}$ via its leaves. Not only is $\mathcal{P}$ a sufficient statistic
for $\mathfrak{T}$, it also completely characterizes the bias and complexity
of the resulting data model conditioned on $\mathfrak{T}$. Indeed, this data
model is just Bayesian linear regression on the indicator functions on the
leaves of $\mathfrak{T}$. Since the functions are orthogonal, the degrees of
freedom of the regression is equal to the size of the partition.

#### Partition ensemble models (PEMs).

A TSE $\mathfrak{E}$ induces an ensemble of partitions
$\mathcal{P}_{1},\mathcal{P}_{2},\ldots,\mathcal{P}_{m}$. Indeed the data
model conditioned on $\mathfrak{E}$ is still a Bayesian linear regression on
the indicator functions of all leaves in $\mathfrak{E}$. However, these
indicators are no longer orthogonal, which means that different ensembles,
with possibly different numbers of leaves, can give rise to the same subspace
of regressors, making $\mathfrak{E}$ not identifiable from data. To avoid this
issue, we directly consider the subspace of the regression function. Formally,
let $\mathbb{V}\subset L^{2}(\mathcal{X}^{d},\nu)$ be the subspace spanned by
indicators of the cells in $\mathcal{P}_{j}$ for $j=1,2,\ldots,m$. We call
this the _partition ensemble model_ (PEM) associated to the TSE $\mathfrak{E}$
and indicate this association via the mapping
$\mathbb{V}=\mathcal{F}(\mathfrak{E})$.

#### Measuring bias.

Let $\Pi_{\mathfrak{E}}$ denote orthogonal projection onto
$\mathcal{F}(\mathfrak{E})$ in $L^{2}(\nu)$. We define the squared (mangitude
of the) bias of $\mathfrak{E}$ with respect to a regression function $f$ as

$\textnormal{Bias}^{2}(\mathfrak{E};f)\coloneqq\int(f-\Pi_{\mathfrak{E}}[f])^{2}d\nu.$
(3)

This is precisely the squared bias of Bayesian linear regression on
$\mathfrak{E}$ if we ignore the regularization effect from the leaf parameter
priors, which is inconsequential in large sample sizes.

#### Measuring complexity.

When conditioned on a TSE $\mathfrak{E}$ and ignoring regularization from leaf
parameter priors, the degrees of freedom of the resulting Bayesian linear
regression model is just the dimension of $\mathcal{F}(\mathfrak{E})$. We
denote this by $\textnormal{df}(\mathfrak{E})$ and use it as a measure of
complexity of $\mathfrak{E}$. Note that this definition does not depend on the
covariate distribution $\nu$ (see Lemma J.2 in the appendix.)

#### Function dimension and optimal sets.

Excessive complexity leads to overfitting and is hence undesirable. To
quantify the excess, we first define the _$m$ -ensemble dimension_ of a
regression function $f$ as

$\dim_{m}(f)\coloneqq\min\left\\{\textnormal{df}(\mathfrak{E})\colon
f\in\mathcal{F}(\mathfrak{E})~{}\text{and}~{}\mathfrak{E}\in\Omega_{\text{TSE},m}\right\\}.$
(4)

In large sample sizes $n$, which is the setting we are concerned with, the
TSEs that result in the smallest MSE must be bias-free. We hence define the
set of optimal TSEs in $\Omega_{\text{TSE},m}$ to be the minimizers of (4).
More generally, we define a series of nested sets with increasing levels of
suboptimality tolerance via:

$\textnormal{OPT}_{m}(f,k)\coloneqq\left\\{\mathfrak{E}\in\Omega_{\text{TSE},m}\colon
f\in\mathcal{F}(\mathfrak{E})~{}\text{and}~{}\textnormal{df}(\mathfrak{E})\leq\dim_{m}(f)+k\right\\}.$

###  4.2 BIC and BART Posterior Concentration

The Bayesian information criterion (BIC) (Schwarz,, 1978) of a TSE
$\mathfrak{E}$ is given by

$\textnormal{BIC}(\mathfrak{E})=\frac{\mathbf{y}^{T}\left(\mathbf{I}-\mathbf{P}_{\mathfrak{E}}\right)\mathbf{y}}{\sigma^{2}}+\textnormal{df}(\mathfrak{E})\log
n+\log\left(2\pi\sigma^{2}\right)n.$

Here, $\mathbf{P}_{\mathfrak{E}}$ refers to projection onto
$\mathcal{F}(\mathfrak{E})$ with respect to the empirical norm $\|\cdot\|_{n}$
(realized as a matrix.) Ignoring the effect of the noise vector $\mathbb{e}$
for now, we see that the first term, divided by the sample size $n$, is an
estimate for the squared bias. Meanwhile, the second term directly measures
the model complexity. Hence, BIC quantifies the quality of a TSE by accounting
for both bias and complexity. Indeed, under our data generative model (Section
2.1) we have the following concentration lemma:

######  Proposition 4.1 (Concentration of BIC differences).

Consider two TSEs $\mathfrak{E}$ and $\mathfrak{E}^{\prime}$ and denote the
difference in their BIC values as
$\Delta\textnormal{BIC}(\mathfrak{E},\mathfrak{E}^{\prime})=\textnormal{BIC}(\mathfrak{E})-\textnormal{BIC}(\mathfrak{E}^{\prime})$.
Then for any $0<\delta<1$, with probability at least $1-\delta$ with respect
to $\mathbb{P}_{n}$, we have

$\Delta\textnormal{BIC}(\mathfrak{E},\mathfrak{E}^{\prime})=\frac{n}{\sigma^{2}}\left(\textnormal{Bias}^{2}(\mathfrak{E};f^{*})-\textnormal{Bias}^{2}(\mathfrak{E}^{\prime};f^{*})\right)+O\left(\sqrt{n\log(1/\delta)}+\log(1/\delta)\right).$
(5)

If furthermore, both TSEs have the same bias, i.e.
$\Pi_{\mathfrak{E}}[f^{*}]=\Pi_{\mathfrak{E}^{\prime}}[f^{*}]$, then we have

$\Delta\textnormal{BIC}(\mathfrak{E},\mathfrak{E}^{\prime})=\log
n\left(\textnormal{df}(\mathfrak{E})-\textnormal{df}(\mathfrak{E}^{\prime})\right)+O\left(\log(1/\delta)\right).$
(6)

From this proposition, we also see that $\textnormal{OPT}_{m}(f^{*},k)$ for
$k=0,1,2,\ldots$ are just sublevel sets of BIC when $n$ is large enough. We
next show that BIC is closely connected to the log marginal likelihood for
TSEs as follows:

######  Proposition 4.2 (Log marginal likelihood and BIC).

Consider a TSE $\mathfrak{E}$. Then for any $0<\delta<1$, there is a minimal
sample size $N$ so that for all $n\geq N$, with probability at least
$1-\delta$ with respect to $\mathbb{P}_{n}$, the log marginal likelihood
satisfies

$\log
p(\mathbf{y}|\mathfrak{E})=-\frac{\textnormal{BIC}(\mathfrak{E})}{2}+O\left(1\right).$

Consequently, the log marginal posterior also satisfies

$\log
p(\mathfrak{E}|\mathbf{y})=-\frac{\textnormal{BIC}(\mathfrak{E})}{2}-\log
p(\mathbf{y})+O\left(1\right).$

This almost linear relationship implies that $\textnormal{OPT}_{m}(f,k)$ for
$k=0,1,2,\ldots$ are also superlevel sets of the marginal posterior. In other
words, they form highest posterior density regions (HPDR). As advertised,
these will be the target sets of our hitting time analysis. Finally, combining
the previous two propositions gives the following result on posterior
concentration.

######  Proposition 4.3 (BART posterior concentration).

For any $0<\delta,\epsilon<1$, there is a minimal sample size $N$ so that for
all $n\geq N$, with probability at least $1-\delta$ with respect to
$\mathbb{P}_{n}$, the marginal posterior measure on $\Omega_{\text{TSE},m}$
satisfies

$p(\textnormal{OPT}_{m}(f^{*},0))~{}|~{}\mathbf{y})>1-\epsilon.$

##  5 Hitting Time Lower Bounds for BART MCMC

Our primary findings consist of lower bounds on HPDR hitting times for BART,
explored across four distinct settings for BART and the data generating
process (DGP). As discussed in the previous section, these regions also
comprise sublevel sets of BIC. We first define hitting times in a general
setting.

#### Hitting times.

Let $(X_{t})$ be a discrete time Markov chain on a finite state space
$\Omega$. Let $\mathcal{A}\subset\Omega$ be a subset. The _hitting time_ of
$\mathcal{A}$ is defined as:

$\tau_{\mathcal{A}}\coloneqq\min\left\\{t\geq 0\colon
X_{t}\in\mathcal{A}\right\\}.$

Note that this is a random variable and that it, in principle, depends on the
initial state $X_{0}$. In our analysis, the initial state is always chosen to
be an ensemble of trivial trees and so will not be referenced in the notation
to avoid clutter.

###  5.1 Square Root and Polynomial Lower Bounds for Additive Models

Our first two lower bounds are for the setting where the DGP is an additive
model.

######  Theorem 5.1 (Lower bound for additive model).

Let $f^{*}$ be an additive function, i.e.

$f^{*}(\mathbf{x})=f_{1}(x_{1})+f_{2}(x_{2})+\cdots+f_{m^{\prime}}(x_{m^{\prime}})$

with $m^{\prime}\geq 2$. Suppose $x_{1},x_{2},\ldots,x_{m^{\prime}}$ are
independent. Suppose $m\leq m^{\prime}$, and we make arbitrary choices for the
other BART hyperparameters $\sigma^{2}$, $\lambda$, $\boldsymbol{\pi}$,
$p_{TSE}$. Then for any $0<\delta<1$, there is a minimal sample size $N$ so
that for all $n\geq N$, with probability at least $1-\delta$ with respect to
$\mathbb{P}_{n}$, the Markov chain induced by
$\textbf{BART}(\mathcal{D}_{n},m,\sigma^{2},\lambda,\boldsymbol{\pi},p_{\mathfrak{E}},-)$
satisfies

$E\left\\{\tau_{\textnormal{OPT}_{m}(f^{*},(q_{\text{max}}-2)(q_{\text{min}}-2)-2)}\right\\}=\Omega\left(n^{1/2}\right),$

where $q_{\text{max}}=\max_{1\leq i\leq m^{\prime}}\dim_{1}(f_{i})$,
$q_{\text{min}}=\min_{1\leq i\leq m^{\prime}}\dim_{1}(f_{i})$.

If furthermore, $m<m^{\prime}$ and we disallow “change” and “swap” moves, i.e.
$\pi_{c}=\pi_{s}=0$, then we have666Recall that
$\textnormal{OPT}_{m}(f^{*},0)\subset\textnormal{OPT}_{m}(f^{*},k)$ for $k\geq
0$. Since $\textnormal{OPT}_{m}(f^{*},k)=\emptyset$ for negative $k$, these
statements can be interpreted as being meaningless unless $q_{\text{max}}$ and
$q_{\text{min}}$ satisfy the relevant constraints.

$E\left\\{\tau_{\textnormal{OPT}_{m}(f^{*},q_{\text{max}}-q_{\text{min}}-1)}\right\\}=\Omega\left(n^{q_{\text{min}}/2-1}\right).$

Additive models are natural generalizations of linear models that have been
widely studied in statistics and machine learning (Hastie and Tibshirani,,
1986; Hastie et al.,, 2009). When fitted to real world datasets, they often
enjoy good prediction accuracy. Hence, we view an additive generative model to
be a natural class of functions for our study.

Furthermore, several works deriving consistency guarantees for frequentist
greedy decision trees and random forests have used additive models as a
generative function class, as the assumption of additivity helps to circumvent
some of the practical and theoretical difficulties arising from greedy
splitting (Scornet et al.,, 2015; Klusowski,, 2021). On the other hand, Tan et
al., 2022a  showed generalization lower bounds for decision trees in this
setting, with the recommendation that models with multiple trees fit additive
models better than those comprising a single tree (see also Tan et al., 2022b
.) Indeed, given the assumptions of Theorem 5.1, we have
$\dim_{l}(f^{*})>\dim_{m}(f^{*})$ for any $l<m$, while
$\dim_{l}(f^{*})=\dim_{m}(f^{*})$ for any $l\geq m$ (see Proposition F.7.) In
other words, a minimum BIC value is achievable if and only if the number of
trees in the BART model is larger than or equal to the number of components in
the additive model.

However, Theorem 5.1 tells us that even when the BART model is correctly
specified and contains an efficient representation of $f^{*}$, i.e. when
$m^{\prime}=m$, the BART sampler may fail to reach a TSE implementing such a
representation within a reasonable time. This adds another perspective to
recent arguments that there may be value in overparameterization (see for
instance Bartlett et al., (2020).) Allowing for more trees than the number of
additive components may empower the sampler with more freedom of navigation to
avoid potential computational bottlenecks. We confirm this conjecture
empirically in our simulations section.

###  5.2 Square Root Lower Bound for Pure Interactions

Our next lower bound is for the setting when the data generating process has a
pure interaction, which we define as follows. Let $x_{i}$ and $x_{j}$ be two
features, i.e. components of $\mathbf{x}\sim\nu$. We say that they form a
_pure interaction_ with respect to a regression function $f^{*}$ if they are
jointly dependent with the response $y$, but are separately conditionally
independent of $y$ for any conditioning set of indices $I$ unless it includes
the other’s index. Mathematically, we can write this as follows:

  * •

$(x_{i},x_{j})\not\perp\\!\\!\\!\perp y$;

  * •

$x_{i}\perp\\!\\!\\!\perp y~{}|~{}x_{I}$ and $x_{j}\perp\\!\\!\\!\perp
y~{}|~{}x_{I}$ for any $I\subset\left\\{1,2,\ldots,d\right\\}$ such that
$i,j\notin I$.

A canonical example of a pure interaction is the exclusive-or (XOR) function
over binary features with a uniform distribution, i.e.
$f(\mathbf{x})=x_{1}x_{2}$ for $\mathcal{X}=\left\\{-1,1\right\\}^{d}$. This
function is well-known to be difficult to learn using either CART (Syrgkanis
and Zampetakis,, 2020; Mazumder and Wang,, 2024) or neural networks (Abbe et
al.,, 2022). As such, it is perhaps unsurprising that the BART sampler also
experiences difficulties in this setting.

######  Theorem 5.2 (Lower bound for pure interaction).

Let $f^{*}$ contain a pure interaction. Suppose we disallow “change” moves,
i.e. $\pi_{c}=0$, and we make arbitrary choices for all other BART
hyperparameters $m$, $\sigma^{2}$, $\lambda$, $\pi_{g}$, $\pi_{p}$, $\pi_{s}$,
$p_{TSE}$. Then for any $0<\delta<1$, there is a minimal sample size $N$ so
that for all $n\geq N$, with probability at least $1-\delta$ with respect to
$\mathbb{P}_{n}$, the Markov chain induced by
$\textbf{BART}(\mathcal{D}_{n},m,\sigma^{2},\lambda,\boldsymbol{\pi},p_{\mathfrak{E}},-)$
satisfies

$E\left\\{\tau_{\textnormal{OPT}_{m}(f^{*},\infty)}\right\\}=\Omega\left(n^{1/2}\right),$

where
$\textnormal{OPT}_{m}(f^{*},\infty)\coloneqq\cup_{k=0}^{\infty}\textnormal{OPT}_{m}(f^{*},k)$.

Note that the suboptimality gap for Theorem 5.2 is much wider than that for
Theorem 5.1. Indeed, it implies that the only TSEs that are reachable within
$o\left(n^{1/2}\right)$ iterations of the sampler have nonzero bias. As such,
the MSE of the BART sampler output does not even converge to zero with the
training sample size, unless we allow for $\Omega(n^{1/2})$ iterations of the
sampler.

###  5.3 Exponential Lower Bound for Bayesian CART

Our final hitting time lower bound shows that the HPDR hitting time for
Bayesian CART can be exponential in the training sample size. This complements
and improves our results in Ronen et al., (2022), which provided an
exponential lower bound for mixing time for Bayesian CART. While the previous
result relied on extremely weak assumptions, the improved version requires a
new assumption on the asymmetry of the regression function $f^{*}$ in terms of
its dependence on different features. Specifically, we say that a regression
function $f^{*}$ has _root dependence_ if there exists a feature $x_{i}$ and
threshold $t$ such that:

  * •

$\text{Corr}^{2}\left(y,\mathbf{1}\left\\{x_{i}\leq t\right\\}\right)>0$;

  * •

$(i,t)$ does not occur as a root split on any tree structure
$\mathfrak{T}\in\textnormal{OPT}_{1}(f,0)$.

An example of such a function is the “staircase” function
$f^{*}(\mathbf{x})=\sum_{j=1}^{s}\prod_{k=1}^{j}\mathbf{1}\left\\{x_{k}>1\right\\}$.
Any feature $x_{i}$ for $2\leq i\leq s$ satisfies the above two properties. On
the other hand, additive functions on independent features do not satisfy
these properties.

######  Theorem 5.3 (Lower bound for Bayesian CART with root dependence).

Suppose $f^{*}$ has root dependence. Suppose $m=1$ and that we disallow
“change” and “swap” moves, i.e. $\pi_{c}=\pi_{s}=0$. Suppose we make arbitrary
choices for all other BART hyperparameters $\sigma^{2}$, $\lambda$, $\pi_{g}$,
$\pi_{p}$, $p_{TSE}$. Then the Markov chain induced by
$\textbf{BART}(\mathcal{D}_{n},1,\sigma^{2},\lambda,\boldsymbol{\pi},p_{\mathfrak{E}},-)$
satisfies

$\liminf_{n\to\infty}\mathbb{E}_{n}\left\\{\frac{\log
E\left\\{\tau_{\textnormal{OPT}_{1}(f^{*},0)}\right\\}}{n}\right\\}\geq\frac{1}{2\sigma^{2}}\left(\int(f^{*})^{2}d\nu-\left(\int
f^{*}d\nu\right)^{2}\right).$

######  Remark 5.4.

Our hitting time lower bounds directly imply mixing time lower bounds in the
space of PEMs, which, as argued in Section 4, are identifiable from data.
Since the bias and degrees of freedom of a TSE $\mathfrak{E}$ is defined in
terms of its associated PEM $\mathcal{F}(\mathfrak{E})$,
$\textnormal{OPT}_{m}(f^{*},k)$ is the preimage under $\mathcal{F}$ of a set
$\widetilde{\textnormal{OPT}}_{m}(f^{*},k)$ in the space of PEMs. Hence, a
hitting time lower bound for $\textnormal{OPT}_{m}(f^{*},k)$ is simultaneously
a hitting time lower bound for $\widetilde{\textnormal{OPT}}_{m}(f^{*},k)$
when considering the induced Markov chain on the space of PEMs. It is easy to
see that this is a lower bound for the mixing time.

######  Remark 5.5.

For the sake of narrative clarity, we have not tried to optimize the
suboptimality gaps (i.e. $k$ in $\textnormal{OPT}_{m}(f^{*},k)$) in our lower
bounds. Note also that we have not attempted to investigate the dependence of
our lower bounds on other data or algorithmic hyperparameters. These of course
influence the minimum sample size $N$ as well as the hidden constant factor in
the Big-Omega notation of our lower bounds.

##  6 Hitting Time Lower Bounds via Barrier Sets

In this section, we briefly outline the proof strategy we use to derive our
hitting time lower bounds. Our first move is to make use of the standard
interpretation of a symmetric Markov chain as a random walk on a network,
whose vertices comprise the states of the Markov chain and whose edges
comprise pairs of states with positive transition probability. We next notice
that hitting times are closely related to escape probabilities, which can be
interpreted as voltages on the network. Voltages can then be calculated using
standard network simplification techniques. Putting these ingredients together
creates the following recipe for deriving hitting time lower bounds:

######  Proposition 6.1 (Recipe for hitting time lower bounds).

Let $\mathfrak{E}_{0},\mathfrak{E}_{1},\mathfrak{E}_{2},\ldots$ denote the
Markov chain induced by a run of
$\textbf{BART}(\mathcal{D}_{n},m,\sigma^{2},\lambda,\boldsymbol{\pi},p_{\mathfrak{E}},-)$
for any fixed dataset $\mathcal{D}_{n}$ and any choice of hyperparameters. Let
$\mathfrak{E}_{\text{bad}}\in\Omega_{\text{TSE},m}$ be a TSE such that, for
some $0<\delta<1$, with probability at least $1-\delta$ with respect to
$\mathbb{P}_{n}$,

$P\left\\{\tau_{\mathfrak{E}_{\text{bad}}}<\tau_{\textnormal{OPT}_{m}(f^{*},k)}\right\\}=\Omega(1).$

Let $\mathcal{B}\subset\Omega_{\text{TSE},m}$ be a subset such that every path
from $\mathfrak{E}_{\text{bad}}$ to $\textnormal{OPT}_{m}(f^{*},k)$ intersects
$\mathcal{B}$. Then with probability at least $1-2\delta$ with respect to
$\mathbb{P}_{n}$, the hitting time of $\textnormal{OPT}_{m}(f^{*},k)$
satisfies

$E\left\\{\tau_{\textnormal{OPT}_{m}(f^{*},k)}\right\\}=\Omega\left(\exp\left(\frac{1}{2}\min_{\mathfrak{E}\in\mathcal{B}}\Delta\textnormal{BIC}(\mathfrak{E},\mathfrak{E}_{\text{bad}})\right)\right).$

Figure 2:  Visual illustration of Proposition 6.1. The chain is initialized at
$\mathfrak{E}_{0}$ and with positive probability hits a suboptimal TSE
$\mathfrak{E}_{\text{bad}}$ before $\textnormal{OPT}_{m}(f^{*},k)$. This
causes to chain to get stuck, as it can only reach
$\textnormal{OPT}_{m}(f^{*},k)$ by passing through an “impassable” barrier set
$\mathcal{B}$.

In other words, our hitting time lower bounds follow immediately once we are
able to identify (i) a suboptimal TSE $\mathfrak{E}_{\text{bad}}$ that is
reachable by the sampler before it hits $\textnormal{OPT}_{m}(f^{*},k)$ and
(ii) a barrier $\mathcal{B}$ that separates $\mathfrak{E}_{\text{bad}}$ from
$\textnormal{OPT}_{m}(f^{*},k)$, i.e. is a vertex cutset, and which has higher
BIC, therefore acting as an impassable “barrier”. We now briefly describe our
choices for each of the three different lower bound settings we have
considered in this paper.

#### Additive models.

For simplicity, we only discuss the case of $m^{\prime}=m$ (when the number of
components is equal to the number of trees.) For $i=1,2,\ldots,m$, denote
$q_{i}=\dim_{1}(f_{i})$, and let
$0=\xi_{i,0}<\xi_{i,1}<\cdots<\xi_{i,q_{i}}=b$ denote the _knots_ of $f_{j}$,
i.e. the values for which $f_{i}(\xi_{i,j})\neq f_{i}(\xi_{i,j}+1)$, together
with the endpoints.777 $\dim_{1}(f_{i})$ is simply the number of constant
pieces of $f_{i}$ or alternatively, one larger than the number of knots of
$f_{i}$. Without loss of generality, assume that $f_{1},\ldots,f_{m}$ are
ordered in descending order of their 1-ensemble dimension, i.e. $q_{1}\geq
q_{2}\geq\cdots\geq q_{m}$.

We now define $\mathfrak{E}_{\text{bad}}$ and a “bad set” $\mathcal{A}$. To
this end, we define a collection of partition models
$\mathbb{V}_{1},\mathbb{V}_{2},\ldots,\mathbb{V}_{m}$ (spans of indicators in
a single partition) as follows. First, for $i=3,4,\ldots,m$,
$j=1,2,\ldots,q_{i}$, define the cells
$\mathcal{L}_{i,j}\coloneqq\\{\mathbf{x}\colon\xi_{i,{j-1}}<x_{i}\leq\xi_{i,j}\\}$,
and set
$\mathbb{V}_{i}=\textnormal{span}\left(\left\\{\mathbf{1}_{\mathcal{L}_{i,j}}\colon
j=1,2,\ldots,q_{i}\right\\}\right)$, i.e. for each $i$, $\mathbb{V}_{i}$
contains splits only on feature $i$ and only at the knots of $f_{i}$. To
introduce inefficiency, we define each of $\mathbb{V}_{1}$ and
$\mathbb{V}_{2}$ to have splits on both features 1 and 2. This construction is
demonstrated in Figure 3. The formal details are fairly involved and will be
deferred to the appendix. We define $\mathcal{A}$ via

$\mathcal{A}\coloneqq\left\\{(\mathfrak{T}_{1},\mathfrak{T}_{2},\ldots,\mathfrak{T}_{m})\colon\mathcal{F}(\mathfrak{T}_{i})=\mathbb{V}_{i}~{}\text{for}~{}i=1,2,\ldots,m\right\\},$

and set $\mathcal{B}$ to be the _outer boundary_ of $\mathcal{A}$, i.e.

$\mathcal{B}=\left\\{\mathfrak{E}\in\Omega_{\text{TSE},m}\colon\mathfrak{E}~{}\text{has
an edge
to}~{}\mathcal{A}~{}\text{and}~{}\mathfrak{E}\notin\mathcal{A}\right\\}.$

Finally, we pick $\mathfrak{E}_{\text{bad}}$ to be a particular element of
$\mathcal{A}$ whose precise construction will be detailed in the Appendix H.
Therein, we will also show that

  * •

$\mathfrak{E}$ has zero bias for all $\mathfrak{E}\in\mathcal{A}$;

  * •

$\mathcal{A}\cap\textnormal{OPT}_{m}(f^{*},(q_{\text{max}}-2)(q_{\text{min}}-2)-2)=\emptyset$;

  * •

$\min_{\mathfrak{E}\in\mathcal{B}}\Delta\textnormal{BIC}(\mathfrak{E},\mathfrak{E}_{\text{bad}})\geq\log
n-O(1)$.

Figure 3:  Visual illustration of the construction for
$\mathfrak{E}_{\text{bad}}$ used in the proof of Theorem 5.1. The left panel
displays the function $f_{1}(x_{1})+f_{2}(x_{2})$ together with all knots of
$f_{1}$ and $f_{2}$. We define $\mathfrak{E}_{\text{bad}}$ so that its first
two trees $\mathfrak{T}_{1}$ and $\mathfrak{T}_{2}$ induce the partitions
$\mathbb{V}_{1}$ and $\mathbb{V}_{2}$ respectively. These combine for a total
of 13 leaves. On the other hand, an optimal TSE will instead make use of
$\mathbb{V}_{1}^{*}$ and $\mathbb{V}_{2}^{*}$, which combine for a total of
only 8 leaves. Nonetheless, we still have
$f_{1}+f_{2}\in\mathbb{V}_{1}^{*}+\mathbb{V}_{2}^{*}$.

#### Pure interactions.

We first define the notion of _reachability_ as follows: Given
$\mathfrak{E},\mathfrak{E}^{\prime}\in\Omega_{\text{TSE},m}$, we say that
$\mathfrak{E}\succsim\mathfrak{E}^{\prime}$ if $\mathfrak{E}$ and
$\mathfrak{E}^{\prime}$ are connected by an edge and if either
$\textnormal{Bias}^{2}(\mathfrak{E};f^{*})>\textnormal{Bias}^{2}(\mathfrak{E}^{\prime};f^{*})$
or
$\textnormal{Bias}^{2}(\mathfrak{E};f^{*})=\textnormal{Bias}^{2}(\mathfrak{E}^{\prime};f^{*})$
and $\textnormal{df}(\mathfrak{E})\geq\textnormal{df}(\mathfrak{E}^{\prime})$.
Note that, because we allow only “grow” and “prune” moves, for adjacent
$\mathfrak{E}$ and $\mathfrak{E}^{\prime}$, $\mathcal{F}(\mathfrak{E})$ and
$\mathcal{F}(\mathfrak{E}^{\prime})$ are nested subspaces, so that
$\textnormal{Bias}^{2}(\mathfrak{E};f^{*})=\textnormal{Bias}^{2}(\mathfrak{E}^{\prime};f^{*})$
if and only if $\Pi_{\mathfrak{E}}[f^{*}]=\Pi_{\mathfrak{E}^{\prime}}[f^{*}]$.
We say that $\mathfrak{E}$ is reachable from $\mathfrak{E}^{\prime}$, denoted
$\mathfrak{E}\succeq\mathfrak{E}^{\prime}$, if there is a sequence of TSEs
$\mathfrak{E}=\mathfrak{E}^{0},\mathfrak{E}^{1},\ldots,\mathfrak{E}^{k}=\mathfrak{E}^{\prime}$
such that $\mathfrak{E}^{i}\succsim\mathfrak{E}^{i+1}$ for $i=0,1,\ldots,k-1$.

Without loss of generality, let $(x_{1},x_{2})$ be a pure interaction for
$f^{*}$. Let $\mathfrak{E}_{\text{bad}}$ be any TSE such that

  * •

$\mathfrak{E}_{\text{bad}}$ is reachable from $\mathfrak{E}_{\emptyset}$;

  * •

There does not exist $\mathfrak{E}\in\Omega_{\text{TSE},m}$ such that
$\mathfrak{E}$ is reachable from $\mathfrak{E}_{\text{bad}}$ but
$\mathfrak{E}_{\text{bad}}$ is not reachable from $\mathfrak{E}$.

Note that such a TSE exists because $\Omega_{\text{TSE},m}$ is finite and
$\succeq$ is a partial ordering on this space. We set $\mathcal{A}$ to be the
equivalence class of $\mathfrak{E}_{\text{bad}}$ under $\succeq$ and set
$\mathcal{B}$ to be the outer boundary of $\mathcal{A}$. We will show in the
appendix that no TSE in $\mathcal{A}$ makes a split on either $x_{1}$ or
$x_{2}$, which implies that
$\mathcal{A}\cap\textnormal{OPT}_{m}(f^{*},\infty)=\emptyset$. We will also
show that
$\min_{\mathfrak{E}\in\mathcal{B}}\Delta\textnormal{BIC}(\mathfrak{E},\mathfrak{E}_{\text{bad}})\geq\log
n-O(1)$.

#### Bayesian CART.

Without loss of generality, let $x_{1}$ be the feature that gives $f^{*}$ root
dependence. By assumption, there is a threshold $t$ such that splitting the
trivial tree on $x_{1}$ at $t$ gives a decrease in squared bias. We set

$\mathcal{A}=\left\\{\mathfrak{T}\in\Omega_{\text{TSE},1}\colon\mathfrak{T}~{}\text{has
root split on}~{}x_{1}~{}\text{at}~{}t\right\\},$

and

$\mathfrak{T}_{\text{bad}}=\arg\min\left\\{\textnormal{BIC}(\mathfrak{E})\colon\mathfrak{E}\in\mathcal{A}\right\\}.$

Note that the outer boundary of $\mathcal{A}$ is a singleton set comprising
the trivial tree $\mathfrak{T}_{\emptyset}$. By assumption, we have
$\mathcal{A}\cap\textnormal{OPT}_{1}(f^{*},0)=\emptyset$. We will show in the
appendix that

$\Delta\textnormal{BIC}(\mathfrak{T}_{\emptyset},\mathfrak{T}_{\text{bad}})\geq\frac{n}{\sigma^{2}}\left(\int(f^{*})^{2}d\nu-\left(\int
f^{*}d\nu\right)^{2}\right)-o(n).$

From these constructions, we also see that the reason why hitting times grow
with training sample size is because the barrier sets become increasingly
difficult to pass through. Heuristically, we can say that this is because the
intrinsic “temperature” of the BART sampler is inversely proportional to the
training sample size.

##  7 Theoretical Limitations and Future Work

In this section, we detail the limitations of our theoretical results, which
naturally suggest directions for future work.

#### Difference from in-practice BART.

The three differences between the version of BART we analyzed and BART as used
in practice are that:

  * •

We assume a fixed noise parameter $\sigma^{2}$ instead of putting a prior on
it;

  * •

At each iteration, we pick a random tree to update instead of cycling
deterministically through the trees;

  * •

We change the rejection probability in the Metropolis filter to be in terms of
the marginal posteriors on TSEs instead of being conditional posteriors on the
updated tree.

We believe that these differences to be relatively minor and that they may
even improve the mixing properties of BART, albeit at the cost of higher
computational complexity per iteration.888Computing the marginal posterior
involves solving a $d$-dimensional linear regression whereas computing the
conditional posterior involves solving a univariate linear regression.

#### Failure to address MSE and coverage directly.

Let $\mathfrak{E}_{0},\mathfrak{E}_{1},\mathfrak{E}_{2},\ldots$ denote the
Markov chain induced by a run of BART. For $j=0,1,2,\ldots$, let $h_{j}$ be a
draw from the conditional posterior on regression functions,
$p(f|\mathfrak{E}_{j},\mathbf{y})$. The BART algorithm returns the collection

$\left\\{h_{t_{\text{burn-in}}+1},h_{t_{\text{burn-in}}+2},\ldots
h_{t_{\text{max}}}\right\\},$ (7)

where $t_{\text{burn-in}}$ is the number of burn-in iterations to be
discarded. The final fitted function is the mean of this collection, while
approximate credible prediction intervals are derived from quantiles. While
our hitting time lower bounds imply that the Markov chain on TSEs fails to
converge and hence that the $p(f|\mathfrak{E}_{j},\mathbf{y})$’s are
separately suboptimal, this does not preclude the _mean_ of (7) having good
MSE performance. It also does not address coverage of the approximate credible
intervals.

#### Restriction to discrete covariates.

We assumed that our covariate distribution was discrete, i.e.
$\mathcal{X}=\left\\{1,2,\ldots,b\right\\}^{d}$. Many real datasets, of
course, contain continuous features. In this case, in-practice BART computes a
grid of quantile values for each continuous feature, and selects a split
threshold only from among these values. Although this effectively makes the
covariate space discrete, it also means that the space varies with the
training sample size $n$. Our proof relies heavily on the uniform
concentration of all node and split-based quantities across a finite set and
so does not automatically generalize to this setting.

#### Failure to address dependence on other DGP parameters.

We did not analyze the dependence of our hitting time lower bounds on $d$, the
dimension of the feature space, $b$, the number of categories for each
discrete feature, $s$, the sparsity of the true regression function, and
$\nu$, the covariate distribution.

#### Failure to address dependence on other algorithmic parameters.

We did not analyze the dependence of our hitting time lower bounds on
$p_{\mathfrak{E}}$, the prior on $\Omega_{\text{TSE},m}$, the move
probabilities $\pi_{g},\pi_{p},\pi_{c},\pi_{s}$, and the variance parameters
$\sigma^{2}$ and $\lambda$.

#### Asymptotic nature of results.

The hitting time lower bounds hold only when the training sample size is
larger than a minimum number $N$. This number depends on the DGP and
algorithmic parameters, and if tracked, can be exponentially large in some of
them.

##  8 Simulations

In this section, we describe the results of a simulation study designed to
bridge some of the theoretical limitations raised in the previous section.
Specifically, we directly study the RMSE of the posterior mean, with respect
to the true regression function, and the empirical coverage of pointwise
posterior credible intervals for the true regression function. We do this for
the output of the original BART algorithm, investigating how they vary
according to various DGP and algorithmic parameters. Our experiments show the
folowing:

  * •

When the data is generated from an additive model, RMSE and coverage for BART
improve with the number of trees, even when the number of trees is larger than
the number of additive components;

  * •

When using Bayesian CART, the root split is often chosen suboptimally and yet
is rarely reversed. The probability that this root split is reversed decreases
as the training sample size increases;

  * •

Across a wide variety of real and synthetic DGPs, RMSE and coverage improve
when averaging the results of multiple BART sampler chains. The relative
improvement gap becomes increasingly pronounced as the training sample size
increases.

The first and second results validate our hitting time lower bounds (Theorem
5.1 and 5.3 respectively) and suggest that they hold for the original BART
algorithm and for reasonable training sample sizes. Since averaging results
from multiple chains should not make a difference if each chain is well-mixed,
the third result suggests that the tendency of mixing and hitting times to
increase with the number of training samples is a fairly general phenomenon
that holds across a wide variety of DGPs.

#### Code availability.

All the code necessary to reproduce the experiments in this section is
publicly available at  github.com/theo-s/bart-hitting-time-sims The computing
infrastructure used was a Linux cluster managed by Department of Statistics at
UC Berkeley. Most runs of the simulation used a single 24-core node with 128
GB of RAM, while the larger datasets required a large-memory node with 792 GB
RAM and 96 cores.

#### Algorithm settings and hyperparameters.

We use the dbarts R package (Dorie,, 2022) with almost all hyperparameters
kept at their default values. In particular,
$\pi_{g}=\pi_{p}=0.25,\pi_{c}=0.1,\pi_{s}=0.4$, and the number of posterior
samples is ndpost=1000. The only exception is that we increase the number of
burn-in iterations from 100 to 5000 (nskip=5000), in order to highlight that
mixing does not occur within a reasonable number of iterations. The responses
are centered and scaled to have variance one. The prior for the noise variance
is calibrated using the residuals of a linear model fit on all of the
predictors.999This is standard for most software implementations of the
algorithm. Unless otherwise noted in the experiment description, the number of
trees was kept equal to 200.

#### Data.

For each experiment, we generate a training dataset $\mathcal{D}_{n}$
comprising $n$ i.i.d. tuples $(\mathbf{x}_{i},y_{i})$, where
$y_{i}=f^{*}(\mathbf{x}_{i})+\epsilon_{i}$,
$\epsilon\sim\mathcal{N}(0,\eta^{2})$. $f^{*}$, $\eta^{2}$, and the covariate
distribution will be varied from experiment to experiment.

#### Evaluation.

Both RMSE and empirical coverage are calculated over an independent test set
consisting of 1000 points drawn from the same distribution. The results are
averaged over 100 experimental replicates, with error bars representing $\pm
1.96$SE. To compute empirical coverage, we count the proportion of the test
set points whose credible interval for the function value contains the ground
truth.

###  8.1 Experiment 1: Under an Additive Model, More Trees Improves
Performance

Our first experiment studies how RMSE and empirical coverage depends on the
interaction of the number of trees specified in the BART model and the number
of components in an additive DGP. We chose the regression function to be one
of the two forms described below. We varied the number of trees in the grid
$\left\\{1,2,\ldots,10\right\\}$ and the training sample size in the grid
$\left\\{10K,20K,50K,100K\right\\}$.

#### DGP.

We let $\mathbf{x}_{i}\in\mathbb{R}^{10}$,
$\mathbf{x}_{i}\sim\mathcal{N}(0,\boldsymbol{\Sigma})$, where we set
$\Sigma_{ii}=1$, $\Sigma_{ij}=0.01$ for $i\neq j$. We let $\eta^{2}=2$, and
take $f^{*}$ to be an additive function with 5 components taking one the
following two functional forms:

  1. 1.

Linear: $f^{*}(\mathbf{x})=0.2x_{1}-x_{2}+0.6x_{3}-0.9x_{4}+0.85x_{5}$;

  2. 2.

Smooth:
$f^{*}(\mathbf{x})=0.2x_{1}^{2}-x_{2}+0.6\cos(x_{3})-0.9|x_{4}|^{1/2}+0.85\sin(x_{5})$.

#### Results:

The results are displayed in Figure 4 and show that the RMSE for BART
decreases while the empirical coverage increases as the number of trees
increases, even when the number of trees is larger than the number of additive
components. This trend holds over both functional forms and is consistent
across different choices of training sample sizes. As discussed in previous
sections, the BART posterior already achieves the maximum goodness of fit
(i.e. BIC is minimized) when the number of trees is equal to the number of
additive components, implying that further improvement in RMSE and empirical
coverage as we increase the number of trees arises purely from better mixing.
This thus corroborates Theorem 5.1. Furthermore, the decrease in RMSE (and
increase in empirical coverage) has a larger slope for larger training sample
sizes, indicating that the improvement in mixing is more pronounced in these
settings.

---
Figure 4:  When fitted to an additive model, RMSE for BART decreases and
coverage increases as the number of trees increases, even when the number of
trees is larger than the number of additive components (i.e. $m>5$.) This
suggests that more trees leads to better mixing and thereby corroborates
Theorem 5.1. This trend holds over different functional forms and several
choices of training sample sizes. Both RMSE and coverage are calculated on an
independent test set and are averaged over 100 experimental replicates, with
error bars representing $\pm 1.96$SE.

###  8.2 Experiment 2: Root Split Gets Stuck on Suboptimal Feature

In the second experiment, we fit Bayesian CART (BART with
$\textit{ntree}={1}$) to DGPs in which the optimal root split is known and
study the root split behavior of the sampler across its iterations. We chose
the regression function to be one of the two forms described below and varied
the training sample size in $\left\\{100,10K\right\\}$. For each of these
choices, we ran the sampler for 500 iterations and counted (i) the percentage
of iterations for which the root split was correct as well as (ii) the number
of total changes to the root split across the 500 iterations.101010We did not
discard any burn-in samples for this investigation. We repeated this process
100 times.

#### DGP.

We let $\mathbf{x}_{i}\in\mathbb{R}^{10}$,
$\mathbf{x}_{i}\sim\mathcal{N}(0,\boldsymbol{\Sigma})$, where we set
$\Sigma_{ii}=1$, $\Sigma_{ij}=0.01$ for $i\neq j$. We set $\eta=1$ and take
$f^{\star}$ to be a piecewise constant function with a single optimal tree
structure, taking one of the following forms (see also Figure 5):

  1. 1.

DGP 1: $f^{*}(\mathbf{x})=\begin{cases}12,&\text{ if }x_{1}>.1\land
x_{2}>.7\\\ 9,&\text{ if }x_{1}>.1\land x_{2}\leq.7\\\ 3,&\text{ if
}x_{1}\leq.1\land x_{3}>-.2\\\ 6,&\text{ if }x_{1}\leq.1\land
x_{3}\leq-.2\end{cases}$

  2. 2.

DGP 2: $f^{*}(\mathbf{x})=\begin{cases}12,&\text{ if }x_{1}>.1\land
x_{2}>.7\land x_{4}>.1\\\ 9,&\text{ if }x_{1}>.1\land x_{2}>.7\land
x_{4}\leq.1\\\ 3,&\text{ if }x_{1}>.1\land x_{2}\leq.7\land x_{5}>.1\text{ or
}x_{1}\leq.1\land x_{3}>-.2\land x_{6}>.2\\\ 4,&\text{ if }x_{1}>.1\land
x_{2}\leq.7\land x_{5}\leq.1\text{ or }x_{1}\leq.1\land x_{3}>-.2\land
x_{6}\leq.2\\\ 2,&\text{ if }x_{1}\leq.1\land x_{3}\leq-.2\land x_{7}>.15\\\
1,&\text{ if }x_{1}\leq.1\land x_{3}\leq-.2\land x_{7}\leq.15\end{cases}$.

In both cases the optimal first split is on $x_{1}$.

Figure 5:  The two regression functions used in Experiment 2 to explore root
split behavior of Bayesian CART.

#### Results.

The distributions of both quantities across the 100 replicates are displayed
in Figure 6. They show that, as the training sample size increases, the
percentage of correct first splits does not converge to 1, but instead seems
to develop a bimodal distribution. Moreover, the number of changes to the root
split decreases. This suggests that the sampler becomes more likely to get
stuck in a set comprising trees with an incorrect root split, which agrees
with Theorem 5.3. What is especially notable is that while Theorem 5.3
requires a restriction to only “grow” and “prune” moves, these results do not
have such a restriction.

Figure 6:  When the DGP has a unique optimal root split, as the training
sample size increases, the Bayesian CART sampler becomes more likely to get
stuck in a set comprising trees with incorrect root splits. For each of 100
experimental replicates, we initialize a sampler at the trivial tree and run
it for 500 iterations. We then calculate (i) the percentage of sampler
iterations with a correct root split and (ii) the total number of changes to
the root split. The top and bottom panels plot the distribution of the first
and second quantities respectively. As the training sample size increases, the
distribution for (i) seems to become bimodal around 0 and 1. Meanwhile, (ii)
becomes more concentrated around 0.

###  8.3 Experiment 3: Multiple BART Chains Improves Mixing

In the final experiment, we investigate the effect of training sample size on
the mixing performance of BART as we fit it to a variety of DGPs, described
below. We varied the number of training samples in the grid
$\left\\{1K,10K,20K,50K,100K\right\\}$ and, for each dataset, varied the
number of chains for the BART sampler in $\\{1,2,5,10\\}$. For each number of
chains, we divided a fixed budget of 1000 posterior samples evenly amongst the
chains and obtained RMSE and empirical coverage values by averaging the
posterior samples thus obtained. Note that we still use the same number of
burn-in iterations for each chain.

#### DGP.

We study the following DGPs:

  1. 1.

Low (Lei and Candès,, 2021): $f^{*}(\mathbf{x})=g(x_{1})g(x_{2}),\quad
g(x)=\frac{2}{1+\exp\\{-12(x-0.5)\\}}$. We let
$\mathbf{x}_{i}\in\mathbb{R}^{10}$,
$\mathbf{x}_{i}\sim\mathcal{N}(0,\boldsymbol{\Sigma})$, where we set
$\Sigma_{ii}=1$, $\Sigma_{ij}=0.01$ for $i\neq j$.

  2. 2.

High (Lei and Candès,, 2021): $f^{*}(\mathbf{x})=g(x_{1})g(x_{2}),\quad
g(x)=\frac{2}{1+\exp\\{-12(x-0.5)\\}}$. We let
$\mathbf{x}_{i}\in\mathbb{R}^{100}$,
$\mathbf{x}_{i}\sim\mathcal{N}(0,\boldsymbol{\Sigma})$, where we set
$\Sigma_{ii}=1$, $\Sigma_{ij}=0.01$ for $i\neq j$.

  3. 3.

Local Sparse Spiky (Behr et al.,, 2021):
$f^{*}(\mathbf{x})=2\cdot\mathbf{1}\left\\{x_{1}<0,x_{3}>0)\right\\}-3\cdot
1_{(x_{5}>0,x_{6}>1}+.8\cdot\mathbf{1}\left\\{x_{3}<1.5,x_{5}<1\right\\}.$ We
let $\mathbf{x}_{i}\in\mathbb{R}^{10}$,
$\mathbf{x}_{i}\sim\mathcal{N}(0,\boldsymbol{\Sigma})$, where we set
$\Sigma_{ii}=1$, $\Sigma_{ij}=0.01$ for $i\neq j$

  4. 4.

Piecewise Linear (Künzel et al.,, 2019):
$f^{*}(\mathbf{x})=\begin{cases}\boldsymbol{\beta}_{1}^{T}\mathbf{x},&\text{
if }x_{20}<-.4\\\ \boldsymbol{\beta}_{2}^{T}\mathbf{x},&\text{ if
}x_{20}<.4\\\ \boldsymbol{\beta}_{3}^{T}\mathbf{x},&\text{ otherwise
}\end{cases}$, where
$\boldsymbol{\beta}_{1},\boldsymbol{\beta}_{2},\boldsymbol{\beta}_{3}\sim\text{Unif}\left([-15,15]^{20}\right)$.
We let $\mathbf{x}_{i}\in\mathbb{R}^{20}$,
$\mathbf{x}_{i}\sim\mathcal{N}(0,\boldsymbol{\Sigma})$, where we set
$\Sigma_{ii}=1$, $\Sigma_{ij}=0.01$ for $i\neq j$.

  5. 5.

Sum: $f^{*}(\mathbf{x})=\sum_{j=1}^{10}x_{j}$. We let
$\mathbf{x}_{i}\in\mathbb{R}^{20}$,
$\mathbf{x}_{i}\sim\mathcal{N}(0,\boldsymbol{\Sigma})$, where we set
$\Sigma_{ij}=\begin{cases}1&\text{if}~{}i=j\\\ 0.01&\text{if}~{}i=j+10\\\
0&\text{otherwise}\end{cases}.$

  6. 6.

Tree: $f^{*}(\mathbf{x})=\mathfrak{T}(\mathbf{x})$, where $\mathfrak{T}$ is a
decision tree function fitted to a standard Gaussian response vector, using
the CART algorithm with maximal depth of 7. We let
$\mathbf{x}_{i}\in\mathbb{R}^{10}$,
$\mathbf{x}_{i}\sim\mathcal{N}(0,\boldsymbol{\Sigma})$, where we set
$\Sigma_{ii}=1$, $\Sigma_{ij}=0.01$ for $i\neq j$

We set $\eta$ such that Tree and Piecewise have signal to noise ratios of 1,
High and Low have signal to noise ratios of .5, and Local Sparse Spiky and Sum
have signal to noise ratios of .75.

#### Results.

The results are displayed in Figure 7. Note that instead of plotting RMSE, we
have chosen to plot relative RMSE, which measures the ratio between the RMSE
obtained from multiple chains and that obtained for a single chain for a given
data setting. The results show that both RMSE and empirical coverage improves,
often quite significantly, as we add multiple chains. This means that
different chains give rise to significantly different distributions, implying
that the BART sampler has not mixed even after the large number of burn-in
iterations. Furthermore, as the number of training samples increases, the
relative performance gap between a single chain and multiple chains increases.
Our results therefore provide evidence that the tendency of HPDR hitting times
to grow with training sample size is consistent across a wider range of DGPs
than was studied theoretically.

#### Real data simulations.

We also performed a similar experiment with a number of benchmark datasets.
Results and further description for this experiment are provided in Appendix
K.

---
Figure 7:  The RMSE (top panel) and empirical coverage (bottom panel) of BART
improve if we average posterior samples from multiple sampler chains, given a
fixed total budget of posterior samples. This improvement means that the BART
sampler has not mixed ever after 5000 burn-in iterations. Furthermore, the
relative performance gaps increases with the number of training samples,
providing evidence that the tendency of HPDR hitting time to grow with
training sample size is consistent across a wide range of DGPs. Both RMSE and
coverage are calculated on an independent test set and are averaged over 100
experimental replicates, with error bars representing $\pm 1.96$SE.

##  9 Discussion

It is widely accepted that Chipman et al., (2010)’s BART sampler often has
issues with mixing and that there is much room for computational improvement.
In this paper, we vastly improve upon our earlier work (Ronen et al.,, 2022)
to provide theoretical computational lower bounds for a slightly modified
version of the BART sampler. Ours is the first work to analyze BART rather
than Bayesian CART. Furthermore, we create a new framework for analysis by
studying hitting times of HPDRs instead of mixing times, which resolves issues
of identifiability and leads to more meaningful computational lower bounds. We
derive these lower bounds under four different, fairly realistic algorithmic
and data generative settings and show in all cases that they grow with the
training sample size. We complement our theoretical results with a simulation
study that validates our results and also suggests that our central thesis,
that BART mixing and hitting times increase with the number of training
samples, is a fairly general phenomenon that holds across a wide variety of
DGPs. We also argue that this is due to the unnecessarily coarse way in which
the BART sampler relates temperature and training sample size.

Our results give theoretical and empirical support to some of the choices that
BART practitioners often already make. Specifically, they support the use of
more trees in the BART ensemble, and they support the use of multiple BART
sampler chains. How to select the optimal number of trees and chains is an
intriguing and important question and will be left to future work. In
addition, our results advocate for the design of better BART samplers and
suggest possible approaches for improvement. First, we believe that there is
great potential in exploring various forms of temperature control. This could
take the form of simulated annealing (Van Laarhoven et al.,, 1987), or
simulated tempering (Marinari and Parisi,, 1992), which has previously been
explored for Bayesian CART (Angelopoulos and Cussens,, 2005), but has yet to
be adapted to BART. Second, we believe that the proposal distribution should
favor more promising split directions, instead of being uniform at random.
There is now a vast literature on how using gradient and even Hessian
information can help to accelerate MCMC in continuous state spaces (see for
instance Neal et al., (2011)), and there is recent work in extending this to
discrete spaces (Zanella,, 2020). Third, we believe instead of constraining
the proposal distribution to “local” moves, it could benefit from
incorporating moves that alter tree structures more drastically. This has been
explored somewhat by Kim and Rockova, (2023).

### Acknowledgements

We would like to thank Antonio Linero and Jungguem Kim for insightful
discussions. We would also like to thank the Statistical Computing Facility
(SCF) of the Department of Statistics at University of California, Berkeley
for computing support, and in particular, Jacob Steinhardt for access to his
group’s large-memory nodes within the SCF Linux cluster.
We gratefully acknowledge partial support from NSF TRIPODS Grant 1740855,
DMS-1613002, 1953191, 2015341, 2209975, 20241842, IIS 1741340, ONR grant
N00014-17-1-2176, NSF grant 2023505 on Collaborative Research: Foundations of
Data Science Institute (FODSI), NSF grant MC2378 to the Institute for
Artificial CyberThreat Intelligence and OperatioN (ACTION), the NSF and the
Simons Foundation for the Collaboration on the Theoretical Foundations of Deep
Learning through awards DMS-2031883 and 814639, and a Weill Neurohub grant. YT
was partially supported by NUS Start-up Grant A-8000448-00-00.

## References

  * Abbe et al., (2022) Abbe, E., Adsera, E. B., and Misiakiewicz, T. (2022).  The merged-staircase property: a necessary and nearly sufficient condition for sgd learning of sparse functions on two-layer neural networks.  In Conference on Learning Theory, pages 4782–4887. PMLR.
  * Agarwal et al., (2022) Agarwal, A., Tan, Y. S., Ronen, O., Singh, C., and Yu, B. (2022).  Hierarchical shrinkage: Improving the accuracy and interpretability of tree-based models.  In International Conference on Machine Learning, pages 111–135. PMLR.
  * Angelopoulos and Cussens, (2005) Angelopoulos, N. and Cussens, J. (2005).  Tempering for Bayesian CART.  In Proceedings of the 22nd International Conference on Machine Learning, pages 17–24.
  * Bartlett et al., (2020) Bartlett, P. L., Long, P. M., Lugosi, G., and Tsigler, A. (2020).  Benign overfitting in linear regression.  Proceedings of the National Academy of Sciences, 117(48):30063–30070.
  * Behr et al., (2021) Behr, M., Wang, Y., Li, X., and Yu, B. (2021).  Provable boolean interaction recovery from tree ensemble obtained via random forests.  arXiv preprint arXiv:2102.11800.
  * Breiman, (2001) Breiman, L. (2001).  Random forests.  Machine Learning, 45(1):5–32.
  * Breiman et al., (1984) Breiman, L., Friedman, J., Olshen, R., and Stone, C. J. (1984).  Classification and regression trees.  Chapman and Hall/CRC.
  * Carnegie, (2019) Carnegie, N. B. (2019).  Comment: Contributions of model features to BART causal inference performance using ACIC 2016 competition data.  Statistical Science, 34(1):90 – 93.
  * Caruana et al., (2008) Caruana, R., Karampatziakis, N., and Yessenalina, A. (2008).  An empirical evaluation of supervised learning in high dimensions.  In Proceedings of the 25th International Conference on Machine Learning, pages 96–103.
  * Caruana and Niculescu-Mizil, (2006) Caruana, R. and Niculescu-Mizil, A. (2006).  An empirical comparison of supervised learning algorithms.  In Proceedings of the 23rd International Conference on Machine Learning, pages 161–168.
  * Castillo and Ročková, (2021) Castillo, I. and Ročková, V. (2021).  Uncertainty quantification for Bayesian CART.  The Annals of Statistics, 49(6):3482–3509.  Publisher: Institute of Mathematical Statistics.
  * Chen and Guestrin, (2016) Chen, T. and Guestrin, C. (2016).  XGBoost: A scalable tree boosting system.  In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pages 785–794.
  * Chipman et al., (1998) Chipman, H. A., George, E. I., and McCulloch, R. E. (1998).  Bayesian CART model search.  Journal of the American Statistical Association, 93(443):935–948.
  * Chipman et al., (2010) Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).  BART: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1):266–298.
  * Dorie, (2022) Dorie, V. (2022).  dbarts: Discrete Bayesian additive regression trees sampler.  R package version 0.9-22.
  * Dorie et al., (2019) Dorie, V., Hill, J., Shalit, U., Scott, M., and Cervone, D. (2019).  Automated versus do-it-yourself methods for causal inference: Lessons learned from a data analysis competition.  Statistical Science, 34(1):43–68.
  * Efron et al., (2004) Efron, B., Hastie, T., Johnstone, I., and Tibshirani, R. (2004).  Least angle regression.  The Annals of Statistics, 32(2):407–499.
  * Fernández-Delgado et al., (2014) Fernández-Delgado, M., Cernadas, E., Barro, S., and Amorim, D. (2014).  Do we need hundreds of classifiers to solve real world classification problems?  The Journal of Machine Learning Research, 15(1):3133–3181.
  * Friedman, (2001) Friedman, J. H. (2001).  Greedy function approximation: a gradient boosting machine.  Annals of Statistics, 29(5):1189–1232.
  * Green and Kern, (2010) Green, D. P. and Kern, H. L. (2010).  Modeling heterogeneous treatment effects in large-scale experiments using Bayesian additive regression trees.  In The Annual Summer Meeting of the Society of Political Methodology, pages 100–110.
  * Green and Kern, (2012) Green, D. P. and Kern, H. L. (2012).  Modeling heterogeneous treatment effects in survey experiments with Bayesian additive regression trees.  Public Opinion Quarterly, 76(3):491–511.
  * Grinsztajn et al., (2022) Grinsztajn, L., Oyallon, E., and Varoquaux, G. (2022).  Why do tree-based models still outperform deep learning on typical tabular data?  Advances in Neural Information Processing Systems, 35:507–520.
  * Hahn et al., (2019) Hahn, P. R., Dorie, V., and Murray, J. S. (2019).  Atlantic causal inference conference (acic) data analysis challenge 2017.  arXiv preprint arXiv:1905.09515.
  * Hastie and Tibshirani, (1986) Hastie, T. and Tibshirani, R. (1986).  Generalized additive models.  Statistical Science, 1(3):297–318.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., Friedman, J. H., and Friedman, J. H. (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction, volume 2.  Springer.
  * He and Hahn, (2021) He, J. and Hahn, P. R. (2021).  Stochastic tree ensembles for regularized nonlinear regression.  Journal of the American Statistical Association, pages 1–20.
  * Hill et al., (2020) Hill, J., Linero, A., and Murray, J. (2020).  Bayesian additive regression trees: A review and look forward.  Annual Review of Statistics and Its Application, 7(1).
  * Hill, (2011) Hill, J. L. (2011).  Bayesian nonparametric modeling for causal inference.  Journal of Computational and Graphical Statistics, 20(1):217–240.
  * Jeong and Rockova, (2020) Jeong, S. and Rockova, V. (2020).  The art of BART: On flexibility of Bayesian forests.  arXiv preprint arXiv:2008.06620, 3(69):146.
  * Kern et al., (2016) Kern, H. L., Stuart, E. A., Hill, J., and Green, D. P. (2016).  Assessing methods for generalizing experimental impact estimates to target populations.  Journal of Research on Educational Effectiveness, 9(1):103–127.
  * Kim and Rockova, (2023) Kim, J. and Rockova, V. (2023).  On mixing rates for Bayesian CART.  arXiv preprint arXiv:2306.00126.
  * Klusowski, (2021) Klusowski, J. M. (2021).  Universal consistency of decision trees in high dimensions.  arXiv preprint arXiv:2104.13881.
  * Künzel et al., (2019) Künzel, S. R., Sekhon, J. S., Bickel, P. J., and Yu, B. (2019).  Metalearners for estimating heterogeneous treatment effects using machine learning.  Proceedings of the National Academy of Sciences, 116(10):4156–4165.
  * Lei and Candès, (2021) Lei, L. and Candès, E. J. (2021).  Conformal inference of counterfactuals and individual treatment effects.  Journal of the Royal Statistical Society: Series B (Statistical Methodology).
  * Levin et al., (2006) Levin, D. A., Peres, Y., and Wilmer, E. L. (2006).  Markov chains and mixing times.  American Mathematical Society.
  * Linero, (2018) Linero, A. R. (2018).  Bayesian regression trees for high-dimensional prediction and variable selection.  Journal of the American Statistical Association, 113(522):626–636.
  * Linero and Yang, (2018) Linero, A. R. and Yang, Y. (2018).  Bayesian regression tree ensembles that adapt to smoothness and sparsity.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(5):1087–1110.
  * Luo and Pratola, (2023) Luo, H. and Pratola, M. T. (2023).  Sharded Bayesian additive regression trees.  arXiv preprint arXiv:2306.00361.
  * Marinari and Parisi, (1992) Marinari, E. and Parisi, G. (1992).  Simulated tempering: a new monte carlo scheme.  Europhysics Letters, 19(6):451.
  * Mazumder and Wang, (2024) Mazumder, R. and Wang, H. (2024).  On the convergence of CART under sufficient impurity decrease condition.  Advances in Neural Information Processing Systems, 36.
  * Nash et al., (1994) Nash, W. J., Sellers, T. L., Talbot, S. R., Cawthorn, A. J., and Ford, W. B. (1994).  The population biology of abalone (haliotis species) in tasmania. i. blacklip abalone (h. rubra) from the north coast and islands of bass strait.  Sea Fisheries Division, Technical Report, 48:p411.
  * Neal et al., (2011) Neal, R. M. et al. (2011).  MCMC using Hamiltonian dynamics.  Handbook of Markov Chain Monte Carlo, 2(11):2.
  * Pace and Barry, (1997) Pace, K. and Barry, R. (1997).  Sparse spatial autoregressions.  Statistics & Probability Letters, 33(3):291–297.
  * Pratola, (2016) Pratola, M. T. (2016).  Efficient metropolis–hastings proposal mechanisms for Bayesian regression tree models.  Bayesian analysis, 11(3):885–911.
  * Pratola et al., (2014) Pratola, M. T., Chipman, H. A., Gattiker, J. R., Higdon, D. M., McCulloch, R., and Rust, W. N. (2014).  Parallel Bayesian Additive Regression Trees.  Journal of Computational and Graphical Statistics, 23(3):830–852.  Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/10618600.2013.841584.
  * Pratola et al., (2020) Pratola, M. T., Chipman, H. A., George, E. I., and McCulloch, R. E. (2020).  Heteroscedastic BART via multiplicative regression trees.  Journal of Computational and Graphical Statistics, 29(2):405–417.
  * Redner and Walker, (1984) Redner, R. A. and Walker, H. F. (1984).  Mixture densities, maximum likelihood and the EM algorithm.  SIAM Review, 26(2):195–239.
  * Rockova and Rousseau, (2021) Rockova, V. and Rousseau, J. (2021).  Ideal Bayesian spatial adaptation.  arXiv preprint arXiv:2105.12793.
  * Ročková and Saha, (2019) Ročková, V. and Saha, E. (2019).  On theory for BART.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2839–2848. PMLR.
  * Ročková and van der Pas, (2020) Ročková, V. and van der Pas, S. (2020).  Posterior concentration for Bayesian regression trees and forests.  The Annals of Statistics, 48(4):2108–2131.
  * Romano et al., (2021) Romano, J. D., Le, T. T., La Cava, W., Gregg, J. T., Goldberg, D. J., Chakraborty, P., Ray, N. L., Himmelstein, D., Fu, W., and Moore, J. H. (2021).  PMLB v1.0: an open-source dataset collection for benchmarking machine learning methods.  Bioinformatics, 38(3):878–880.
  * Ronen et al., (2022) Ronen, O., Saarinen, T., Tan, Y. S., Duncan, J., and Yu, B. (2022).  A mixing time lower bound for a simplified version of BART.  arXiv preprint arXiv:2210.09352.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of Statistics, 6(2):461–464.
  * Scornet et al., (2015) Scornet, E., Biau, G., and Vert, J.-P. (2015).  Consistency of random forests.  Annals of Statistics, 43(4):1716–1741.
  * Starling et al., (2020) Starling, J. E., Murray, J. S., Carvalho, C. M., Bukowski, R. K., and Scott, J. G. (2020).  BART with targeted smoothing: An analysis of patient-specific stillbirth risk.  The Annals of Applied Statistics, 14(1):28–50.
  * Syrgkanis and Zampetakis, (2020) Syrgkanis, V. and Zampetakis, M. (2020).  Estimation and inference with trees and forests in high dimensions.  In Conference on Learning Theory, pages 3453–3454. PMLR.
  * (57) Tan, Y. S., Agarwal, A., and Yu, B. (2022a).  A cautionary tale on fitting decision trees to data from additive models: generalization lower bounds.  In International Conference on Artificial Intelligence and Statistics, pages 9663–9685. PMLR.
  * (58) Tan, Y. S., Singh, C., Nasseri, K., Agarwal, A., and Yu, B. (2022b).  Fast interpretable greedy-tree sums (FIGS).  arXiv preprint arXiv:2201.11931.
  * Van Laarhoven et al., (1987) Van Laarhoven, P. J., Aarts, E. H., van Laarhoven, P. J., and Aarts, E. H. (1987).  Simulated Annealing.  Springer.
  * Vershynin, (2018) Vershynin, R. (2018).  High-Dimensional Probability: An Introduction with Applications in Data Science, volume 47.  Cambridge University Press.
  * Wainwright, (2019) Wainwright, M. J. (2019).  High-Dimensional Statistics: A Non-Asymptotic Viewpoint, volume 48.  Cambridge University Press.
  * Wendling et al., (2018) Wendling, T., Jung, K., Callahan, A., Schuler, A., Shah, N. H., and Gallego, B. (2018).  Comparing methods for estimation of heterogeneous treatment effects using observational data from health care databases.  Statistics in medicine, 37(23):3309–3324.
  * Wu et al., (2007) Wu, Y., Tjelmeland, H., and West, M. (2007).  Bayesian CART: Prior specification and posterior simulation.  Journal of Computational and Graphical Statistics, 16(1):44–66.
  * Yeager et al., (2019) Yeager, D. S., Hanselman, P., Walton, G. M., Murray, J. S., Crosnoe, R., Muller, C., Tipton, E., Schneider, B., Hulleman, C. S., Hinojosa, C. P., et al. (2019).  A national experiment reveals where a growth mindset improves achievement.  Nature, 573(7774):364–369.
  * Zanella, (2020) Zanella, G. (2020).  Informed proposals for local MCMC in discrete spaces.  Journal of the American Statistical Association, 115(530):852–865.

##  Appendix A Literature summary for BART posterior concentration results

We list here all results on posterior concentration that we are aware of:

  * •

Posterior concentration at the optimal $n^{-\alpha/(2\alpha+p)}$ rate for
BCART and BART when the regression function is Hölder $\alpha$-smooth,
$0<\alpha\leq 1$ (Ročková and Saha,, 2019).

  * •

Posterior concentration at the optimal $n^{-\alpha/(2\alpha+s)}$ rate for
BCART and BART when the regression function is Hölder $\alpha$-smooth,
$0<\alpha\leq 1$, and $s$-sparse (Ročková and van der Pas,, 2020).

  * •

Posterior concentration at the almost optimal $n^{-\alpha/(2\alpha+1)}\log n$
rate for BART when the regression function is an additive sum of univariate
components, each of which is Hölder $\alpha$-smooth, $0<\alpha\leq 1$ (Ročková
and van der Pas,, 2020).

  * •

Extensions of the above results to $\alpha>1$, provided BART is modified to
allow for “soft” splits (Linero and Yang,, 2018).

  * •

Posterior concentration at the optimal rate (up to log factors) for BART when
the regression function is piecewise heterogeneous anisotropic Hölder smooth
(Jeong and Rockova,, 2020; Rockova and Rousseau,, 2021).

##  Appendix B Proofs for Section 4.2

We first introduce some additional notation that will be used throughout the
rest of the appendix. We denote $D_{f}=\|f^{*}\|_{\infty}$ and let
$D_{\epsilon}$ denote the sub-Gaussian parameter of the noise random variable
$\epsilon$ (Wainwright,, 2019). $C$ will be used to denote a universal
constant (i.e. not depending on any parameters) that may change from line to
line. For simplicity, many of our results are written using big-$O$ notation,
which will indicate leading order dependence on the sample size $n$ as well as
the error probability $\delta$.

###  B.1 Main Proofs

######  Lemma B.1 (Concentration of empirical risk difference).

Let $\mathfrak{E}$ and $\mathfrak{E}^{\prime}$ be two partition ensemble
models. Then for any $0<\delta<1$, with probability at least $1-\delta$, the
risk difference between them satisfies

$\left|\mathbf{y}^{T}\left(\mathbf{P}_{\mathfrak{E}}-\mathbf{P}_{\mathfrak{E}^{\prime}}\right)\mathbf{y}-n\int\left(\Pi_{\mathfrak{E}}f^{*}\right)^{2}-\left(\Pi_{\mathfrak{E}^{\prime}}f^{*}\right)^{2}d\nu\right|=O\left(\sqrt{n\log(1/\delta)}\right).$
(8)

If furthermore, $\Pi_{\mathfrak{E}}f^{*}=\Pi_{\mathfrak{E}^{\prime}}f^{*}$,
then the above bound can be improved to

$\left|\mathbf{y}^{T}\left(\mathbf{P}_{\mathfrak{E}}-\mathbf{P}_{\mathfrak{E}^{\prime}}\right)\mathbf{y}\right|=O\left(\log(1/\delta)\right).$
(9)

#### Proof of Proposition 4.1

Since

$\Delta\textnormal{BIC}(\mathfrak{E},\mathfrak{E}^{\prime})=\frac{\mathbf{y}^{T}\left(\mathbf{P}_{\mathcal{\mathfrak{E}}}-\mathbf{P}_{\mathfrak{E}^{\prime}}\right)\mathbf{y}}{\sigma^{2}}+(\textnormal{df}(\mathfrak{E})-\textnormal{df}(\mathfrak{E}^{\prime}))\log
n,$

the desired concentration follows immediately from Lemma B.1.∎

######  Lemma B.2 (Log marginal likelihood formula).

Let $\mathfrak{E}$ be a tree ensemble structure (TSE), and let
$\boldsymbol{\Psi}$ be an $n\times b$ matrix whose columns comprise the
indicators of the leaves in $\mathfrak{E}$. The log marginal likelihood
satisfies

$\displaystyle-2\log p(\mathbf{y}|\mathbf{X},\mathfrak{E})$
$\displaystyle=n\log\left(2\pi\sigma^{2}\right)+\log\det\left(\lambda^{-1}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\mathbf{I}\right)$
$\displaystyle\quad\quad+\frac{1}{\sigma^{2}}\left(\left\|\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}-\mathbf{y}\right\|_{2}^{2}+\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\left(\mathbf{I}-\boldsymbol{\Psi}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)^{-1}\boldsymbol{\Psi}^{T}\right)\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}\right),$
(10)

where

$\hat{\boldsymbol{\mu}}_{LS}\coloneqq\operatorname*{\arg\\!\min}_{\boldsymbol{\mu}}\|\boldsymbol{\Psi}\boldsymbol{\mu}-\mathbf{y}\|_{2}^{2}$

is the solution to the least squares problem.

######  Lemma B.3 (Concentration of log-determinent).

With the notation of Lemma B.2, denote
$\hat{\boldsymbol{\Sigma}}\coloneqq\frac{1}{n}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}$
and
$\boldsymbol{\Sigma}\coloneqq\mathbb{E}\left\\{\hat{\boldsymbol{\Sigma}}\right\\}$.
Then for any $0<\delta<1$, for
$n\geq\max\left\\{64m^{3/2}s_{\min}^{-2}\log(2\textnormal{df}(\mathfrak{E})/\delta),4\lambda
s_{\min}^{-1}\right\\}$, with probability at least $1-\delta$, we have

$\left|\log\det\left(\lambda^{-1}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\mathbf{I}\right)-\textnormal{df}(\mathfrak{E})\log(n)-\sum_{i=1}^{\textnormal{df}(\mathfrak{E})}\log(s_{i}/\lambda)\right|=O\left(\sqrt{\frac{\log(1/\delta)}{n}}\right),$
(11)

where $s_{1},s_{2},\ldots,s_{\textnormal{df}(\mathfrak{E})}$ denote the values
of the nonzero eigenvalues of $\boldsymbol{\Sigma}$.

######  Lemma B.4 (Concentration of error term).

With the notation of Lemma B.3, for any $0<\delta<1$, when $n\geq
16m^{3/2}s_{\min}^{-2}\log(2\textnormal{df}(\mathfrak{E})/\delta)$, where
$s_{\min}$ is the minimum nonzero eigenvalue of $\boldsymbol{\Sigma}$, with
probability at least $1-\delta$, we have

$\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\left(\mathbf{I}-\boldsymbol{\Psi}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)^{-1}\boldsymbol{\Psi}^{T}\right)\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}=O\left(1\right).$
(12)

#### Proof of Proposition 4.2.

Starting with equation (B.2) from Lemma B.2, plug in equations (11) and (12)
from Lemma B.2 and Lemma B.3 respectively. Notice that
$\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}=\mathbf{P}_{\mathfrak{E}}\mathbf{y}$,
so that

$\left\|\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}-\mathbf{y}\right\|_{2}^{2}=\left\|\left(\mathbf{I}-\mathbf{P}_{\mathfrak{E}}\right)\mathbf{y}\right\|_{2}^{2}=\mathbf{y}^{T}\left(\mathbf{I}-\mathbf{P}_{\mathfrak{E}}\right)\mathbf{y}.$

This completes the proof.∎

#### Proof of Proposition 4.3.

Let $\mathfrak{E}^{*}\in\textnormal{OPT}_{m}(f^{*},0)$. We will show that for
any $\mathfrak{E}\not\in\textnormal{OPT}_{m}(f^{*},0)$, there is some
$N_{\mathfrak{E}}$ such that

$p(\mathfrak{E}|\mathbf{y})\leq\frac{\epsilon
p(\mathfrak{E}^{*}|\mathbf{y})}{\left|\Omega_{\text{TSE},m}\right|}$ (13)

with probability at least $1-\delta/\left|\Omega_{\text{TSE},m}\right|$ for
all $n\geq N_{\mathfrak{E}}$. If this is true, set
$N=\max_{\mathfrak{E}\notin\textnormal{OPT}_{m}(f^{*},0)}N_{\mathfrak{E}}$ and
take $n\geq N$. On the intersection of all these events, we have

$p(\textnormal{OPT}_{m}(f^{*},0)^{c}~{}|\mathbf{y})=\sum_{\mathfrak{E}\notin\textnormal{OPT}_{m}(f^{*},0)}p(\mathfrak{E}|\mathbf{y})\leq\sum_{\mathfrak{E}\notin\textnormal{OPT}_{m}(f^{*},0)}\frac{\epsilon
p(\mathfrak{E}^{*}|\mathbf{y})}{\left|\Omega_{\text{TSE},m}\right|}\leq\epsilon.$

To prove (13), fix $\mathfrak{E}$. Using Proposition 4.2 and Proposition 4.1,
we get a probability $1-\delta/|\Omega_{\text{TSE},m}|$ event on which

$\displaystyle\log p(\mathfrak{E}^{*}|\mathbf{y})-\log
p(\mathfrak{E}|\mathbf{y})$
$\displaystyle\quad=\begin{cases}\frac{n}{2\sigma^{2}}\left(\textnormal{Bias}^{2}(\mathfrak{E};f^{*})\right)+O\left(\sqrt{n\log(1/\delta)}+\log(1/\delta)\right)&\text{if}~{}\textnormal{Bias}^{2}(\mathfrak{E};f^{*})\neq
0\\\ \frac{\log
n}{2}\left(\textnormal{df}(\mathfrak{E})-\textnormal{df}(\mathfrak{E}^{*})\right)+O\left(\log(1/\delta)\right)&\text{otherwise}.\end{cases}$

In either case, we get
$\frac{p(\mathfrak{E}|\mathbf{y})}{p(\mathfrak{E}^{*}|\mathbf{y})}\to 0$ as
$n\to\infty$.∎

###  B.2 Further Details

#### Proof of Lemma B.1.

Recall that $\Pi_{\mathfrak{E}}$ refers to orthogonal projection onto
$\mathcal{F}(\mathfrak{E})$ in $L^{2}(\nu)$, while $\mathbf{P}_{\mathfrak{E}}$
refers to orthogonal projection onto $\mathcal{F}(\mathfrak{E})$ with respect
to the empirical norm $\|\cdot\|_{n}$. With this in mind, decompose
$y=\epsilon+(f^{*}(x)-\Pi_{\mathfrak{E}}f^{*}(x))+\Pi_{\mathfrak{E}}f^{*}(x)$
and write this in vector form as
$\mathbf{y}=\boldsymbol{\epsilon}+\left(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}}\right)+\mathbf{f}^{*}_{\mathfrak{E}}$.
Since $\Pi_{\mathfrak{E}}f^{*}\in\mathcal{F}(\mathfrak{E})$, we have
$\mathbf{P}_{\mathfrak{E}}\mathbf{f}^{*}_{\mathfrak{E}}=\mathbf{f}^{*}_{\mathfrak{E}}$.
We can then therefore expand the quadratic form
$\mathbf{y}^{T}\mathbf{P}_{\mathfrak{E}}\mathbf{y}$ as follows:

$\displaystyle\mathbf{y}^{T}\mathbf{P}_{\mathfrak{E}}\mathbf{y}$
$\displaystyle=\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\boldsymbol{\epsilon}+(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})+2\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})+(\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}\mathbf{f}^{*}_{\mathfrak{E}}+2\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\mathbf{f}^{*}_{\mathfrak{E}}+2(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}\mathbf{f}^{*}_{\mathfrak{E}}$
$\displaystyle=\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\boldsymbol{\epsilon}+(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})+2\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})+(\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{f}^{*}_{\mathfrak{E}}+2\boldsymbol{\epsilon}^{T}\mathbf{f}^{*}_{\mathfrak{E}}+2(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{f}^{*}_{\mathfrak{E}}.$
(14)

Note that $\epsilon$, $(f^{*}(x)-\Pi_{\mathfrak{E}}f^{*}(x))$, and
$\Pi_{\mathfrak{E}}f^{*}(x)$ are uncorrelated random variables, with
$\epsilon$ being also independent of the other two variables. This implies
that the third, fifth and sixth terms in (B.2) have zero mean. On the other
hand, because of finite sample fluctuations,
$(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})$ and
$\mathbf{f}^{*}_{\mathfrak{E}}$, are not necessarily orthogonal as vectors.

To bound the expectation of (B.2), first observe that
$f^{*}-\Pi_{\mathfrak{E}}f^{*}$ and $\Pi_{\mathfrak{E}}f^{*}$ are bounded
random variables and thus have both standard deviation and sub-Gaussian norm
bounded by $D_{f}$ (Wainwright,, 2019). We then compute

$\displaystyle\mathbb{E}\left\\{\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\boldsymbol{\epsilon}+(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})+2\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})+(\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{f}^{*}_{\mathfrak{E}}+2\boldsymbol{\epsilon}^{T}\mathbf{f}^{*}_{\mathfrak{E}}+2(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{f}^{*}_{\mathfrak{E}}\right\\}$
$\displaystyle=~{}$
$\displaystyle\mathbb{E}\left\\{\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\boldsymbol{\epsilon}\right\\}+\mathbb{E}\left\\{(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})\right\\}+\mathbb{E}\left\\{(\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{f}^{*}_{\mathfrak{E}}\right\\}$
$\displaystyle=~{}$
$\displaystyle\textnormal{tr}\left\\{\mathbf{P}_{\mathfrak{E}}\right\\}\left(\textnormal{Var}\left\\{\epsilon\right\\}+\textnormal{Var}\left\\{f^{*}-f^{*}_{\mathfrak{E}}\right\\}\right)+n\int\left(\Pi_{\mathfrak{E}}f^{*}\right)^{2}d\nu$
$\displaystyle\leq~{}$ $\displaystyle
C\left(D_{\epsilon}+D_{f}\right)\textnormal{df}(\mathfrak{E})+n\int\left(\Pi_{\mathfrak{E}}f^{*}\right)^{2}d\nu.$
(15)

Next, we bound the fluctuations of each term in (B.2) separately. Using the
Hanson-Wright inequality (Wainwright,, 2019; Vershynin,, 2018), we get
$1-\delta$ probability events over which

$\left|\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\boldsymbol{\epsilon}-\mathbb{E}\left\\{\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\boldsymbol{\epsilon}\right\\}\right|\leq
CD_{\epsilon}\max\left\\{\log(1/\delta),\sqrt{\textnormal{df}(\mathfrak{E})\log(1/\delta)}\right\\},$
(16)

and

$\left|(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})-\mathbb{E}\left\\{(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})\right\\}\right|\leq
CD_{f}\max\left\\{\log(1/\delta),\sqrt{\textnormal{df}(\mathfrak{E})\log(1/\delta)}\right\\}.$
(17)

Using Hoeffding’s inequality (Wainwright,, 2019), we have further $1-\delta$
probability events over which

$\left|(\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{f}^{*}_{\mathfrak{E}}-n\int\left(\Pi_{\mathfrak{E}}f^{*}\right)^{2}d\nu\right|\leq
CD_{f}^{2}\sqrt{n\log(1/\delta)},$ (18)
$\left|\boldsymbol{\epsilon}^{T}\mathbf{f}_{\mathfrak{E}}^{*}\right|\leq
CD_{f}D_{\epsilon}\sqrt{n\log(1/\delta)},$ (19)
$\left|\left(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}}\right)^{T}\mathbf{f}_{\mathfrak{E}}^{*}\right|\leq
CD_{f}^{2}\sqrt{n\log(1/\delta)}.$ (20)

For the third term in (B.2), we use Cauchy-Schwarz followed by Young’s
inequality to get

$\displaystyle
2\left|\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})\right|$
$\displaystyle\leq
2\left(\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\boldsymbol{\epsilon}\right)^{1/2}\left((\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})\right)^{1/2}$
$\displaystyle\leq\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\boldsymbol{\epsilon}+(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\mathbf{P}_{\mathfrak{E}}(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}}).$
(21)

Conditioning on all the events guaranteeing (16) to (20) and plugging in these
bounds together with (B.2) into (B.2), we get

$\displaystyle\left|\mathbf{y}^{T}\mathbf{P}_{\mathfrak{E}}\mathbf{y}-n\int\left(\Pi_{\mathfrak{E}}f^{*}\right)^{2}d\nu\right|$
$\displaystyle\leq~{}$ $\displaystyle
C\left(\left(D_{\epsilon}+D_{f}\right)\textnormal{df}(\mathfrak{E})+\left(D_{f}+D_{\epsilon}\right)\log(1/\delta)+D_{f}(D_{f}+D_{\epsilon})\sqrt{n\log(1/\delta)}+(D_{f}+D_{\epsilon})\sqrt{\textnormal{df}(\mathfrak{E})\log(1/\delta)}\right)$
$\displaystyle=~{}$ $\displaystyle
O\left(\sqrt{n\log(1/\delta)}+\log(1/\delta)\right).$

Repeating the same argument for $\mathfrak{E}^{\prime}$ and adjusting $\delta$
so that the intersection of all events conditioned on has probability at least
$1-\delta$ completes the proof of (8).

To prove (9), observe that under the additional assumption, we can cancel
terms in (B.2) to get

$\mathbf{y}^{T}\left(\mathbf{P}_{\mathfrak{E}}-\mathbf{P}_{\mathfrak{E}^{\prime}}\right)\mathbf{y}=\boldsymbol{\epsilon}^{T}\left(\mathbf{P}_{\mathfrak{E}}-\mathbf{P}_{\mathfrak{E}^{\prime}}\right)\boldsymbol{\epsilon}+(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})^{T}\left(\mathbf{P}_{\mathfrak{E}}-\mathbf{P}_{\mathfrak{E}^{\prime}}\right)(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}})+2\boldsymbol{\epsilon}^{T}\left(\mathbf{P}_{\mathfrak{E}}-\mathbf{P}_{\mathfrak{E}^{\prime}}\right)(\mathbf{f}^{*}-\mathbf{f}^{*}_{\mathfrak{E}}).$

Applying (B.2) followed by (16) and (17) completes the proof.

#### Proof of Lemma B.2.

Recall that the full log likelihood satisfies

$p(\mathbf{y}|\mathbf{X},\mathfrak{E},\boldsymbol{\mu})=\left(2\pi\sigma^{2}\right)^{-n/2}\exp\left(-\frac{\left\|\boldsymbol{\Psi}\boldsymbol{\mu}-\mathbf{y}\right\|_{2}^{2}}{2\sigma^{2}}\right),$

while conditioned on $\mathfrak{E}$, the prior on $\boldsymbol{\mu}$ satisfies

$p(\boldsymbol{\mu}|\mathfrak{E})=\left(2\pi\sigma^{2}\lambda^{-1}\right)^{-b/2}\exp\left(-\frac{\lambda\|\boldsymbol{\mu}\|_{2}^{2}}{2\sigma^{2}}\right),$

where $b$ is the number of columns in $\boldsymbol{\Psi}$. Hence

$p(\mathbf{y}|\mathbf{X},\mathfrak{E},\boldsymbol{\mu})p(\boldsymbol{\mu}|\mathfrak{E})=\left(2\pi\sigma^{2}\right)^{-n/2}\left(2\pi\sigma^{2}\lambda^{-1}\right)^{-b/2}\exp\left(-\frac{1}{2\sigma^{2}}\left(\left\|\boldsymbol{\Psi}\boldsymbol{\mu}-\mathbf{y}\right\|_{2}^{2}+\lambda\left\|\boldsymbol{\mu}\right\|_{2}^{2}\right)\right).$
(22)

Consider the orthogonal decomposition

$\left\|\boldsymbol{\Psi}\boldsymbol{\mu}-\mathbf{y}\right\|_{2}^{2}=\left\|\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}-\mathbf{y}\right\|_{2}^{2}+\left\|\boldsymbol{\Psi}(\boldsymbol{\mu}-\hat{\boldsymbol{\mu}}_{LS})\right\|_{2}^{2}.$

We next add the second term on the right to the exponent in the prior and
complete the square:

$\displaystyle\left\|\boldsymbol{\Psi}(\boldsymbol{\mu}-\hat{\boldsymbol{\mu}}_{LS})\right\|_{2}^{2}+\lambda\|\boldsymbol{\mu}\|_{2}^{2}$
$\displaystyle=~{}$
$\displaystyle\boldsymbol{\mu}^{T}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)\boldsymbol{\mu}-2\left(\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}\right)^{T}\boldsymbol{\mu}+\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}$
$\displaystyle=~{}$
$\displaystyle\left(\boldsymbol{\mu}-\boldsymbol{\mu}_{0}\right)^{T}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)\left(\boldsymbol{\mu}-\boldsymbol{\mu}_{0}\right)-\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)^{-1}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}+\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}.$
(23)

where

$\boldsymbol{\mu}_{0}=\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)^{-1}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}.$

The constant term in (B.2) is

$\displaystyle-\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)^{-1}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}+\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}$
$\displaystyle=~{}$
$\displaystyle\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\left(\mathbf{I}-\boldsymbol{\Psi}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)^{-1}\boldsymbol{\Psi}^{T}\right)\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}.$
(24)

Plugging (B.2) back into (22) and integrating, we get

$\displaystyle p(\mathbf{y}|\mathbf{X},\mathfrak{E})$ $\displaystyle=\int
p(\mathbf{y}|\mathbf{X},\mathfrak{E},\boldsymbol{\mu})p(\boldsymbol{\mu}|\mathfrak{E})d\boldsymbol{\mu}$
$\displaystyle=\left(2\pi\sigma^{2}\right)^{-n/2}\exp\left(-\frac{\left\|\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}-\mathbf{y}\right\|_{2}^{2}+\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\left(\mathbf{I}-\boldsymbol{\Psi}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)^{-1}\boldsymbol{\Psi}^{T}\right)\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}}{2\sigma^{2}}\right)$
$\displaystyle\quad\cdot\int\left(2\pi\sigma^{2}\lambda^{-1}\right)^{-b/2}\exp\left(-\frac{\left(\boldsymbol{\mu}^{T}-\boldsymbol{\mu}_{0}\right)^{T}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)\left(\boldsymbol{\mu}^{T}-\boldsymbol{\mu}_{0}\right)}{2\sigma^{2}}\right)d\boldsymbol{\mu}.$
(25)

By a change of variables, the integral can be computed as

$\int\left(2\pi\sigma^{2}\lambda^{-1}\right)^{-b/2}\exp\left(-\frac{\left(\boldsymbol{\mu}^{T}-\boldsymbol{\mu}_{0}\right)^{T}\left(\lambda^{-1}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\mathbf{I}\right)\left(\boldsymbol{\mu}^{T}-\boldsymbol{\mu}_{0}\right)}{2\sigma^{2}\lambda^{-1}}\right)d\boldsymbol{\mu}=\det\left(\lambda^{-1}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\mathbf{I}\right)^{-1/2}.$

Plugging this back into (B.2) and taking logarithms yields (B.2).

#### Proof of Lemma B.3.

Let $\hat{s}_{1},\hat{s}_{2},\ldots,\hat{s}_{b}$ be the eigenvalues of
$\hat{\Sigma}$. Using Lemma B.5, we have $\hat{s}_{i}=0$ for any
$i>\textnormal{df}(\mathfrak{E})$. We may therefore compute

$\displaystyle\log\det\left(\lambda^{-1}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\mathbf{I}\right)$
$\displaystyle=\sum_{i=1}^{\textnormal{df}(\mathfrak{E})}\log(n\lambda^{-1}\hat{s}_{i}+1)$
$\displaystyle=\sum_{i=1}^{\textnormal{df}(\mathfrak{E})}\log\left(\frac{\hat{s}_{i}+\lambda/n}{s_{i}}\right)+\textnormal{df}(\mathfrak{E})\log
n+\sum_{i=1}^{\textnormal{df}(\mathfrak{E})}\log(s_{i}/\lambda).$

It remains to bound the first term. To this end, we first condition on the
$1-\delta$ probability event guaranteed by Lemma B.6. Then, we observe that

$\displaystyle\left|\frac{\hat{s}_{i}+\lambda/n}{s_{i}}-1\right|$
$\displaystyle\leq\frac{1}{s_{i}}\left(\left|\hat{s}_{i}-s_{i}\right|+\frac{\lambda}{n}\right)$
$\displaystyle\leq\frac{1}{s_{\min}}\left(\left\|\hat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}\right\|+\frac{\lambda}{n}\right)$
$\displaystyle\leq\frac{1}{s_{\min}}\left(\max\left\\{\sqrt{\frac{4m^{3/2}\log(2\textnormal{df}(\mathfrak{E})/\delta)}{n}},\frac{4\sqrt{m}\log(2\textnormal{df}(\mathfrak{E})/\delta)}{n}\right\\}+\frac{\lambda}{n}\right).$
(26)

Here, the second inequality makes use of Weyl’s inequality.

Recall the elementary inequality

$\left|\log x\right|\leq 2\left|x-1\right|$

for $0<x<1/2$. Using this together with (B.2), we get

$\sum_{i=1}^{\textnormal{df}(\mathfrak{E})}\log\left(\frac{\hat{s}_{i}+\lambda/n}{s_{i}}\right)=O\left(\sqrt{\frac{\log(1/\delta)}{n}}\right)$

when
$n\geq\max\left\\{64m^{3/2}s_{\min}^{-2}\log(2\textnormal{df}(\mathfrak{E})/\delta),4\lambda
s_{\min}^{-1}\right\\}$.∎

#### Proof of Lemma B.4.

Recall that
$\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}=\mathbf{P}_{\mathfrak{E}}\mathbf{y}$.
We thus rewrite and bound the error term as

$\displaystyle\hat{\boldsymbol{\mu}}_{LS}^{T}\boldsymbol{\Psi}^{T}\left(\mathbf{I}-\boldsymbol{\Psi}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)^{-1}\boldsymbol{\Psi}^{T}\right)\boldsymbol{\Psi}\hat{\boldsymbol{\mu}}_{LS}$
$\displaystyle=\mathbf{y}^{T}\mathbf{P}_{\mathfrak{E}}\mathbf{M}\mathbf{P}_{\mathfrak{E}}\mathbf{y}$
$\displaystyle\leq\left\|\mathbf{P}_{\mathfrak{E}}\mathbf{y}\right\|_{2}^{2}\|\mathbf{M}\|,$
(27)

where

$\mathbf{M}=\mathbf{P}_{\mathfrak{E}}-\boldsymbol{\Psi}\left(\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}+\lambda\mathbf{I}\right)^{-1}\boldsymbol{\Psi}^{T}.$

Using similar arguments as in the proof of Lemma B.1, we bound

$\displaystyle\left\|\mathbf{P}_{\mathfrak{E}}\mathbf{y}\right\|_{2}^{2}$
$\displaystyle=(\mathbf{f}_{\mathfrak{E}}^{*})^{T}\mathbf{f}_{\mathfrak{E}}^{*}+\boldsymbol{\epsilon}^{T}\mathbf{P}_{\mathfrak{E}}\boldsymbol{\epsilon}$
$\displaystyle\leq
D_{f}^{2}\left(n+\sqrt{n\log(1/\delta)}\right)+D_{\epsilon}^{2}\left(n+\sqrt{\textnormal{df}(\mathfrak{E})\log(1/\delta)}\right).$
(28)

Meanwhile, the nonzero eigenvalues of $\mathbf{M}$ are of the form

$1-\frac{\hat{s}_{i}}{\hat{s}_{i}+\lambda/n}=\frac{\lambda}{n\hat{s}_{i}+\lambda}$

for $i=1,2,\ldots,\textnormal{df}(\mathfrak{E})$. These can be further bounded
as

$\displaystyle\frac{\lambda}{n\hat{s}_{i}+\lambda}$
$\displaystyle\leq\frac{\lambda}{n(s_{i}-\left|\hat{s}_{i}-s_{i}\right|)}$
$\displaystyle\leq\frac{\lambda}{n\left(s_{\min}-\left\|\hat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}\right\|\right)}.$
(29)

Taking $n\geq
16m^{3/2}s_{\min}^{-2}\log(2\textnormal{df}(\mathfrak{E})/\delta)$ and
conditioning on the $1-\delta$ probability event guaranteed by Lemma B.6, we
can further bound (B.2) by $\lambda/2ns_{\min}$, which gives
$\|\mathbf{M}\|=O\left(n^{-1}\right)$. Combining this with (B.2) and plugging
them back into (B.2) finishes the proof.∎

######  Lemma B.5 (Rank of empirical covariance matrix).

With the notation of Lemma B.3, we have
$\textnormal{rank}(\hat{\boldsymbol{\Sigma}})\leq\textnormal{rank}\left(\boldsymbol{\Sigma}\right)$.

###### Proof.

Let $\mathbf{l}_{1},\mathbf{l}_{2},\ldots,\mathbf{l}_{n}$ denote the rows of
$\boldsymbol{\Psi}$, noting that they are i.i.d. random vectors. Note that if
$\boldsymbol{\Sigma}\mathbf{v}=0$ for some $\mathbf{v}$, this implies that
$\textnormal{Cov}\left\\{\mathbf{v}^{T}\mathbf{l}\right\\}=\mathbf{v}^{T}\boldsymbol{\Sigma}\mathbf{v}=0$,
and so $\mathbf{v}^{T}\mathbf{l}\equiv 0$ as a random variable. In particular,
we have $\mathbf{v}^{T}\mathbf{l}_{i}=0$ for $i=1,2\ldots,n$, and we also get
$\hat{\boldsymbol{\Sigma}}\mathbf{v}=0$. As such, the nullspace for
$\boldsymbol{\Sigma}$ is contained within the nullspace for
$\hat{\boldsymbol{\Sigma}}$. The conclusion follows. ∎

######  Lemma B.6 (Concentration of empirical covariance matrix).

With the notation of Lemma B.3, for any $0<\delta<1$, with probability at
least $1-\delta$, we have

$\left\|\frac{1}{n}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}-\boldsymbol{\Sigma}\right\|\leq\max\left\\{\sqrt{\frac{4m^{3/2}\log(2\textnormal{df}(\mathfrak{E})/\delta)}{n}},\frac{4\sqrt{m}\log(2\textnormal{df}(\mathfrak{E})/\delta)}{n}\right\\}.$

###### Proof.

Let $\mathbf{l}_{1},\mathbf{l}_{2},\ldots,\mathbf{l}_{n}$ denote the rows of
$\boldsymbol{\Psi}$ as before. Since each point can only be contained in a
single leaf on each tree, we have
$\left\|\mathbf{l}_{j}\right\|_{2}=\sqrt{m}$, while $\boldsymbol{\Sigma}$ also
satisfies $\left\|\boldsymbol{\Sigma}\right\|\leq m$. Using Corollary
6.20111111While Corollary 6.20 is stated with the assumption that the rows
have mean zero, the proof in Wainwright, (2019) illustrates that this is
unnecessary. in Wainwright, (2019), we therefore have

$\left\|\frac{1}{n}\boldsymbol{\Psi}^{T}\boldsymbol{\Psi}-\boldsymbol{\Sigma}\right\|\leq
2\textnormal{df}(\mathfrak{E})\exp\left(-\frac{nt^{2}}{2\sqrt{m}(m+t)}\right)$

for any $t>0$. Rearranging this equation completes the proof. ∎

##  Appendix C Background on Markov Chains

For the whole of this section, let $X_{0},X_{1},\ldots$ be an irreducible and
aperiodic discrete time Markov chain on a finite state space $\Omega$, with
stationary distribution $\pi$.

###  C.1 Networks and Voltages

#### Harmonic functions.

Let $P$ be the transition kernel of $(X_{t})$. We call a function
$h\colon\Omega\to\mathbb{R}$ harmonic for $P$ at a state $x$ if

$h(x)=\sum_{y\in\Omega}P(x,y)h(y).$ (30)

######  Lemma C.1 (Uniqueness of harmonic extensions, Proposition 9.1. in
Levin et al., (2006)).

Let $\mathcal{A}\subset\Omega$ be a subset of the state space. Let
$h_{\mathcal{A}}\colon\mathcal{A}\to\mathbb{R}$ be a function defined on
$\mathcal{A}$. The function $h\colon\Omega\to\mathbb{R}$ defined by
$h(x)\coloneqq\mathbb{E}\\{h_{\mathcal{A}}(X_{\tau_{\mathcal{A}}})|X_{0}=x\\}$
is the unique extension of $h_{\mathcal{A}}$ such that
$h(x)=h_{\mathcal{A}}(x)$ for all $x\in\mathcal{A}$ and $h$ is harmonic for
$P$ at all $x\in\Omega\backslash\mathcal{A}$.

The values of a harmonic function can be computed by solving the system of
linear equations given by (30) for each $x\in\Omega\backslash\mathcal{A}$.
This is hard to do directly by hand for complicated state spaces, but when the
Markov chain is symmetric, i.e. the stationary distribution satisfies
$\pi(x)P(x,y)=\pi(y)P(y,x)$ for any $x,y$, we can use several operations to
_simplify the state space while preserving the values of the harmonic function
on the remaining state space_. Indeed, harmonic functions are equivalent to
voltages on electrical circuits, and it is well-known how to simplify circuits
in order to calculate voltages:

  1. 1.

(Gluing) Points on the circuit with the same voltage can be joined.

  2. 2.

(Series law) Two resistors in series with resistances $r_{1}$ and $r_{2}$ can
be merged with the new resistor having resistance $r_{1}+r_{2}$.

  3. 3.

(Parallel law) Two resistors in parallel with resistances $r_{1}$ and $r_{2}$
can be merged with the new resistor having resistance $1/(1/r_{1}+1/r_{2})$.

We make this connection rigorous by introducing the following definitions and
via the subsequent lemmas.

#### Networks, conductance, resistance.

A _network_ $(\Omega,c)$ is a tuple comprising a finite state space $\Omega$
and a symmetric function $c\colon\Omega\times\Omega\to\mathbb{R}_{+}$ called
the _conductance_. The _resistance_ function is defined as
$r(x,y)=\frac{1}{c(x,y)}$, and can take the value of positive infinity. We say
that $\\{x,y\\}$ is an _edge_ in the network if $c(x,y)>0$. Any network
$(\Omega,c)$ has an associated Markov chain whose transition probabilities are
defined by $P(x,y)=\frac{c(x,y)}{\sum_{z\in\Omega}c(x,z)}$.

#### Voltage and current flow.

We say that a function is harmonic on the network if it is harmonic with
respect to $P$. Given $a,z\in\Omega$, a _voltage_ $W$ between $a$ and $z$ is a
function that is harmonic on $\Omega\backslash\\{a,z\\}$. The _current flow_
$I\colon\Omega\times\Omega\to\mathbb{R}$ associated with $W$ is defined as
$I(x,y)=c(x,y)\left(W(x)-W(y)\right)$. The strength of the current flow is
defined as

$\|I\|\coloneqq\sum_{y\in\Omega}c(a,y)\left(W(a)-W(y)\right).$

#### Effective resistance and effective conductance.

The _effective resistance_ between $a$ and $z$ is defined as

$R(a\leftrightarrow z)\coloneqq\frac{W(a)-W(z)}{\|I\|},$

noting that this is independent of the choice of $W$ by the uniqueness
property in Lemma C.1. The _effective conductance_ between $a$ and $z$ is
defined as $C(a\leftrightarrow z)\coloneqq 1/R(a\leftrightarrow z)$.

######  Lemma C.2 (Network simplification rules).

Consider a network $(\Omega,c)$. Define the following operations that each
produces a modified network $(\Omega^{\prime},c^{\prime})$.

  1. 1.

(Gluing) Given $u,v\in\Omega$, define
$\Omega^{\prime}\coloneqq\Omega\backslash\\{v\\}$ and

$c^{\prime}(x,y)=\begin{cases}c(x,y)&x,y\neq u,\\\ c(x,u)+c(x,v)&x\neq
u,y=u,\\\ c(u,y)+c(v,y)&x=u,y\neq u\\\ c(u,v)+c(u,y)+c(u,x)&x=y=u.\end{cases}$

  2. 2.

(Parallel and series laws) Given $u,v,w\in\Omega$ with $c(v,x)=0$ for all
$x\notin\\{u,w\\}$, define $\Omega^{\prime}\coloneqq\Omega\backslash\\{v\\}$
and

$c^{\prime}(x,y)=\begin{cases}c(u,w)+\frac{c(u,v)c(v,w)}{c(u,v)+c(v,w)}&(x,y)=(u,w)~{}\text{or}~{}(w,u),\\\
c(x,y)&\text{otherwise.}\end{cases}$

Consider a function $h$ that is harmonic on $\Omega\backslash\mathcal{A}$. The
following hold:

  1. 1.

If we glue states $u,v\in\Omega$ such that $h(u)=h(v)$, then $h$ remains
harmonic on $\Omega^{\prime}\backslash\mathcal{A}$ with respect to the
modified transition matrix $P^{\prime}$.

  2. 2.

If we apply the parallel and series laws to $u,v,w\in\Omega$ with
$v\notin\mathcal{A}$, then $h$ remains harmonic on
$\Omega^{\prime}\backslash\mathcal{A}$ with respect to the modified transition
matrix $P^{\prime}$.

Furthermore, if $h$ is a voltage between two states $a,z\in\Omega$, then
applying the operations does not change the effective conductance and
resistance between them.

###### Proof.

The first statement is obvious as the mean value equation for harmonic
functions can be repeated almost verbatim. For the second statement, we just
have to check the mean value equation for $h$ at $u$. This is equivalent to
the equation

$\sum_{x\in\Omega^{\prime}}c^{\prime}(u,x)\left(h(x)-h(u)\right)=0.$ (31)

First note that under the original network, our assumption on $c(v,x)$ and the
mean value property at $v$ gives

$c(v,w)\left(h(w)-h(v)\right)=c(u,v)\left(h(v)-h(u)\right).$ (32)

Next, we compute

$\displaystyle c^{\prime}(u,w)\left(h(w)-h(u)\right)$
$\displaystyle=\left(c(u,w)+\frac{c(u,v)c(v,w)}{c(u,v)+c(v,w)}\right)\left(h(w)-h(u)\right)$
$\displaystyle=c(u,w)\left(h(w)-h(u)\right)+\frac{c(u,v)c(v,w)}{c(u,v)+c(v,w)}\left((h(w)-h(v))+(h(v)-h(u)\right)$
$\displaystyle=c(u,w)\left(h(w)-h(u)\right)+\frac{c(u,v)^{2}+c(u,v)c(v,w)}{c(u,v)+c(v,w)}\left(h(v)-h(u)\right)$
$\displaystyle=c(u,w)\left(h(w)-h(u)\right)+c(v,w)\left(h(v)-h(u)\right),$
(33)

where the third equality follows from (32). We therefore have

$\sum_{x\in\Omega^{\prime}}c^{\prime}(u,x)\left(h(x)-h(u)\right)=\sum_{x\in\Omega}c(u,x)\left(h(x)-h(u)\right),$

and the mean value equation at $u$ for the original network implies (31).

Finally, to conclude invariance of effective conductance, we observe that it
is is defined in terms of voltages and current flows. We have already shown
that voltages are unchanged, so we just need to argue that the strength of the
current flow is similarly unchanged. This is immediate whenever
$a\notin\\{u,w\\}$. When $a=u$, this follows from (C.1). ∎

######  Lemma C.3 (Rayleigh’s Monotonicity Law, Theorem 9.12 in Levin et al.,
(2006)).

Given a network with two resistance functions
$r,r^{\prime}\colon\Omega\times\Omega\to\mathbb{R}_{+}\cup\\{\infty\\}$,
suppose that $r\leq r^{\prime}$ pointwise. Then we have

$R(a\leftrightarrow z;r)\leq R(a\leftrightarrow z;r^{\prime}).$

###  C.2 Hitting Precedence Probabilities

#### Hitting precedence probabilities.

Let $\mathcal{A},\mathcal{B}\subset\Omega$ be disjoint subsets. The _hitting
precedence probability_ of $\mathcal{A}$ relative to $\mathcal{B}$ is defined
as the following function on $\Omega$:

$\textnormal{HPP}(x;\mathcal{A},\mathcal{B})\coloneqq
P\left\\{\tau_{\mathcal{A}}<\tau_{\mathcal{B}}~{}|~{}X_{0}=x\right\\}.$

######  Lemma C.4.

$\textnormal{HPP}(-;\mathcal{A},\mathcal{B})$ is a harmonic function on
$\Omega\backslash\left(\mathcal{A}\cup\mathcal{B}\right)$.

###### Proof.

Write $h_{\mathcal{A}\cup\mathcal{B}}(z)=\mathbf{1}\\{z\in\mathcal{A}\\}$. It
is easy to see that

$h_{\mathcal{A}\cup\mathcal{B}}(X_{\tau_{\mathcal{A}\cup\mathcal{B}}})=\begin{cases}1&\tau_{\mathcal{A}}<\tau_{\mathcal{B}}\\\
0&\tau_{\mathcal{A}}>\tau_{\mathcal{B}}.\end{cases}$

Hence,

$\mathbb{E}\left\\{h_{\mathcal{A}\cup\mathcal{B}}(X_{\tau_{\mathcal{A}\cup\mathcal{B}}})|X_{0}=x\right\\}=\mathbb{P}\left\\{\tau_{\mathcal{A}}<\tau_{\mathcal{B}}~{}|~{}X_{0}=x\right\\}=\textnormal{HPP}(x;\mathcal{A},\mathcal{B}).$

By Lemma C.1, the left hand side is a harmonic function on
$\Omega\backslash\left(\mathcal{A}\cup\mathcal{B}\right)$. ∎

######  Lemma C.5 (HPP from a bottleneck state).

Given a network $(\Omega,c)$ with states $a,x,z$ and such that $c(u,z)=0$ for
all $u\notin\\{x,z\\}$. Let $x=x_{0},x_{1},\ldots,x_{k}=a$ be any sequence of
states such that there is some $\rho>0$ for which
$c(x_{i-1},x_{i})\geq\rho^{-1}c(x,z)$ for $i=1,\ldots,k$. Then

$\mathbb{P}\\{\tau_{a}<\tau_{z}|X_{0}=x\\}\geq\frac{1}{k\rho+1}.$

###### Proof.

Let $(\Omega,c^{\prime})$ be the modified network in which we set

$c^{\prime}(u,v)=\begin{cases}c(u,v)&\\{u,v\\}\in\\{\\{x_{i-1},x_{i}\\}\colon
i=1,\ldots,k\\}\cup\\{\\{x,z\\}\\}\\\ 0&\text{otherwise}.\end{cases}$

Then we have

$\displaystyle\mathbb{P}\\{\tau_{a}<\tau_{z}|X_{0}=x\\}$
$\displaystyle=\frac{r(x,z)}{R(a\leftrightarrow x;r)+r(x,z)}$
$\displaystyle\geq\frac{r(x,z)}{R(a\leftrightarrow x;r^{\prime})+r(x,z)},$
(34)

where the equality uses Lemma C.6 and the inequality uses Lemma C.3. Next,
using the series law from Lemma C.2, we have

$\displaystyle R(a\leftrightarrow x;r^{\prime})$
$\displaystyle=\sum_{i=1}^{k}r(x_{i-1},x_{i})$ $\displaystyle\leq k\rho
r(x,z).$

Plugging this into (C.2) and cancelling $r(x,z)$ in the numerator and
denominator completes the proof. ∎

######  Lemma C.6.

Given a network $(\Omega,c)$ with states $a,x,z$ and such that $c(u,z)=0$ for
all $u\notin\\{x,z\\}$. Then we have

$P\\{\tau_{a}<\tau_{z}~{}|~{}X_{0}=x\\}=\frac{C(a\leftrightarrow
x)}{C(a\leftrightarrow x)+c(x,z)}.$

###### Proof.

Let $h(y)\coloneqq\mathbb{P}\\{\tau_{a}<\tau_{z}|X_{0}=y\\}$, and note that
$h$ is a voltage between $a$ and $z$. As such, we have

$R(a\leftrightarrow z)=\frac{h(a)-h(z)}{\|I\|}=\frac{1}{\|I\|}.$

On the other hand, $h$ is also a voltage between $a$ and $x$ on the reduced
state space $\Omega\backslash\\{z\\}$, which gives

$R(a\leftrightarrow x)=\frac{h(a)-h(x)}{\|I\|}.$

Finally, by the series law, we have

$R(a\leftrightarrow z)=R(a\leftrightarrow x)+r(x,z).$

Putting everything together, we get

$\displaystyle P\\{\tau_{a}<\tau_{z}|X_{0}=x\\}$
$\displaystyle=1-\left(h(a)-h(x)\right)$ $\displaystyle=1-R(a\leftrightarrow
x)\|I\|$ $\displaystyle=1-\frac{R(a\leftrightarrow x)}{R(a\leftrightarrow z)}$
$\displaystyle=\frac{r(x,z)}{R(a\leftrightarrow x)+r(x,z)}$
$\displaystyle=\frac{C(a\leftrightarrow x)}{C(a\leftrightarrow x)+c(x,z)},$

as we wanted. ∎

##  Appendix D Proof of Theorem 5.2

We will first present the proofs for Theorem 5.2 and Theorem 5.3 because they
are relatively simple compared to that of Theorem 5.1. For convenience, we
repeat the relevant constructions and definitions here.

#### Reachability.

Given $\mathfrak{E},\mathfrak{E}^{\prime}\in\Omega_{\text{TSE},m}$, we say
that $\mathfrak{E}\succsim\mathfrak{E}^{\prime}$ if $\mathfrak{E}$ and
$\mathfrak{E}^{\prime}$ are connected by an edge, and if either
$\textnormal{Bias}^{2}(\mathfrak{E};f^{*})>\textnormal{Bias}^{2}(\mathfrak{E}^{\prime};f^{*})$
or
$\textnormal{Bias}^{2}(\mathfrak{E};f^{*})=\textnormal{Bias}^{2}(\mathfrak{E}^{\prime};f^{*})$
and $\textnormal{df}(\mathfrak{E})\geq\textnormal{df}(\mathfrak{E}^{\prime})$.
Note that, because we allow only “grow” and “prune” moves, for adjacent
$\mathfrak{E}$ and $\mathfrak{E}^{\prime}$, $\mathcal{F}(\mathfrak{E})$ and
$\mathcal{F}(\mathfrak{E}^{\prime})$ are nested subspaces, so that
$\textnormal{Bias}^{2}(\mathfrak{E};f^{*})=\textnormal{Bias}^{2}(\mathfrak{E}^{\prime};f^{*})$
if and only if
$\Pi_{\mathcal{F}(\mathfrak{E})}[f^{*}]=\Pi_{\mathcal{F}(\mathfrak{E}^{\prime})}[f^{*}]$.
We say that $\mathfrak{E}$ is reachable from $\mathfrak{E}^{\prime}$, denoted
$\mathfrak{E}\succeq\mathfrak{E}^{\prime}$, if there is a sequence of TSEs