Posterior Concentration for Bayesian Regression Trees and Forests

Non-parametric Regression Setup

The remarkable empirical success of Bayesian tree-based regression has raised considerable interest in understanding why and when these methods produce good results. Despite their extensive use in practice, theoretical justifications have, thus far, been unavailable. To narrow this yawning gap, we consider the fundamental problem of making inference about an unknown regression function.

Our setup consists of the nonparametric regression model

where output variables $\boldsymbol{Y}^{(n)}=(Y_{1},\dots,Y_{n})^{\prime}$ are related in a stochastic fashion to a set of $p$ potential covariates $\boldsymbol{x}_{i}=(x_{i1},\dots,x_{ip})^{\prime},1\leq i\leq n$. We assume that the covariates $\boldsymbol{x}_{i}=(x_{i1},\dots,x_{ip})^{\prime}$ are fixed and have been rescaled so that $\boldsymbol{x}_{i}\in^{p}$. The true unknown regression surface $f_{0}(\boldsymbol{x})$ will be assumed to be smooth, possibly involving only a small fraction $q_{0}$ of the $p$ potential covariates.

In the absence of a parametric model, a natural strategy to estimate the unknown regression function is by partitioning the covariate space into cells and then estimating the function locally within each cell from available observations. Such strategies yield histogram reconstructions of the regression surface and have been analyzed theoretically by multiple authors . Regression trees are among the most popular data-dependent histogram methods, where the partitioning scheme is obtained through nested parallel axis splitting. Trees are an integral constituent of ensemble methods that aggregate single tree learners into forests to boost prediction . Tree-based regression, either single or ensemble, is arguably one of the most popular machine learning tools today. In particular, Bayesian variants of these methods (Bayesian CART and BART) have earned a prominent role as one of the top machine learners. While consistency results for classical trees and random forests have been available , theory on the also very widely used Bayesian counterparts is non-existent. Our goal in this paper is to provide the first frequentist optimality results for Bayesian trees, and their ensembles.

Most of the work on Bayesian nonparametric regression has revolved around Gaussian processes . While there are multiple results on recursive partitioning (or histogram) priors for Bayesian density estimation (or non-linear autoregression ), the literature on Bayesian regression histograms is far more deserted. One fundamental contribution is due to Coram and Lalley , who showed consistency of Bayesian binary regression with uniform mixture priors on step functions and with one predictor. More recently, van der Pas and Ročková considered a similar setup for estimating step mean functions in Gaussian regression, again with a single predictor. This paper goes far beyond that framework, addressing (a) the full-fledged high-dimensional setup with a diverging number of potential covariates, (b) tree ensembles for additive regression.

The purpose of this paper is to study the rate of convergence of posterior distributions induced by step function priors on the regression surface when $p>n$. The speed of convergence is measured by the size of the smallest shrinking ball around $f_{0}$ that contains most of the posterior probability. In pioneering works, Ghosal, Ghosh and Van der Vaart and Shen and Wasserman obtained rates of convergence for infinite-dimensional parametric models with iid observations in terms of the size of the model (measured by the metric entropy) and concentration rate of the prior around $f_{0}$. These results were later extended to infinite-dimensional models that are not iid by Ghosal and Van der Vaart . Their general conceptual framework serves as an umbrella for our development.

We initially assume that $f_{0}$ is Hölder continuous (with smoothness not exceeding one) and may depend only on a small fraction of $q_{0}$ predictors. The optimal rate of estimation of a $q_{0}$-variable function, which is known to be $\alpha$-smooth, is $n^{-\alpha/(2\alpha+q_{0})}$ . Our first result shows that, with suitable regularization priors, single Bayesian regression trees achieve this minimax rate (up to a log factor). In other words, the posterior behaves nearly as well as if we knew $\alpha$ and the number of active covariates $q_{0}$, concentrating at a rate that only depends on the number of active predictors. This is the first optimality result for Bayesian regression trees, demonstrating their adaptability and reluctance to overfit in high-dimensional scenarios with $p>n$. The regularization is achieved through our proposed spike-and-tree prior, a new variant of the Bayesian CART prior for dimension reduction and model-free variable selection. Going further, we show that Bayesian additive regression trees also achieve (near) minimax-rate optimal performance when approximating a single smooth function. Finally, the tree ensembles are also shown to be certifiably optimal when the true function is an actual sum of smooth functions, again concentrating at a near minimax rate.

2 Notation

3 Outline

We outline our goals and strategy in Section 2. We then review several useful concepts for analyzing recursive partitioning schemes in Section 3. In Section 4, we state our first main result on the posterior concentration for Bayesian CART. In Section 5, we develop tools for analyzing Bayesian additive regression trees and show their optimal posterior concentration in non-additive regression. Section 6 presents the final development of our theory concerning the recovery of an additive regression function with additive trees. We conclude with a discussion in Section 7. The proofs of our main theorems are presented in Sections 8, 9 and 10.

Background

The key to our approach will be drawing upon the foundational posterior concentration theory for non-iid observations, laid down in the seminal paper by Ghosal and Van der Vaart .

Our results are obtained under a unifying hat of a general result which requires three conditions to hold. Namely, suppose that for a sequence $\varepsilon_{n}^{2}\rightarrow 0$ such that $n\,\varepsilon_{n}^{2}$ is bounded away from zero and sets $\mathcal{F}_{n}\subset\mathcal{F}$ we have

for some $d>2$. Then it follows from Theorem 4 in that the posterior distribution concentrates at the rate $\varepsilon_{n}^{2}$, i.e.

in $P_{f_{0}}^{n}$-probability, as $n\rightarrow\infty$, for any $M_{n}\rightarrow\infty$.

The conditions of Ghosal and Van der Vaart provide a very general recipe for showing posterior concentration in infinite-dimensional models. Our goal in this paper is to obtain tailored statements for Bayesian regression trees and forests. The major challenge will be (a) designing a sequence of approximating spaces (a sieve) $\mathcal{F}_{n}\subset\mathcal{F}$ and (b) endowing $\mathcal{F}$ with a prior distribution such that the three conditions hold simultaneously for $\varepsilon_{n}$ as small as possible. To this end, we will build on, and develop, tools for an analysis of recursive partitioning schemes.

Throughout this work, we assume that the true regression function $f_{0}$ is Hölder continuous and the smoothness parameter does not exceed one, (or an additive composition thereof), in the sense made precise below.

With $\mathcal{H}^{\alpha}_{p}$ we denote the space of uniformly $\alpha$-Hölder continuous functions, i.e.

where $\alpha\in(0,1]$ and where $||f||_{\mathcal{H}^{\alpha}}$ is the Hölder coefficient.

The assumption $\alpha\leq 1$ is standard in the study of piecewise constant estimators and priors, see for example . The reason for this limitation is that step functions are relatively rough; e.g. the approximation error of histograms for functions that are smoother than Lipschitz is at least of the order $1/K$, where $K$ is the number of bins. The number of steps required to approximate a smooth function well is thus too large, creating a costly bias-variance tradeoff.

In some applications, it is reasonable to expect that the regression function $f_{0}$ depends only on a small fraction of input covariates. For a set of indices $\mathcal{S}\subseteq\{1,\dots,p\}$, we define

Regime 1: $f_{0}$ is $\alpha$-Hölder continuous and depends on an unknown subset $\mathcal{S}_{0}$ of $|S_{0}|=q_{0}$ covariates, i.e. $f_{0}\in\mathcal{H}^{\alpha}_{p}\,\cap\,\mathcal{C}(\mathcal{S}_{0})$.

Regime 2: $f_{0}$ is an aggregate of $T_{0}$ $\alpha^{t}$-Hölder continuous functions $f_{0}^{t}$, $1\leq t\leq T_{0}$, each depending on an unknown subset $\mathcal{S}_{0}^{t}$ of $|S_{0}^{t}|=q_{0}^{t}$ covariates, i.e. $f_{0}(\boldsymbol{x})=\sum_{t=1}^{T_{0}}f_{0}^{t}(\boldsymbol{x})$, where $f_{0}^{t}\in\mathcal{H}_{p}^{\alpha^{t}}\,\cap\,\mathcal{C}(\mathcal{S}_{0}^{t})$.

For an estimation procedure to be successful in Regime 1, it needs to be doubly adaptive (with respect to $\alpha$ and $q_{0}$). We will show that both single trees and tree ensembles adapt accordingly, performing at a near-minimax rate. Regime 2 makes the performance discrepancies more apparent, where the additive structure of $f_{0}$ is appreciated by tree ensembles, which can achieve a faster convergence rate than single trees. A variant of Regime 2 was studied by who derived minimax rates for estimating additive smooth functions and showed optimal concentration of additive Gaussian processes. Here, we approximate $f_{0}$ with step functions and their aggregates. We limit considerations to step functions that are underpinned by recursive partitioning schemes.

On Recursive Partitions

Tree-based regression consists of first finding an underlying partitioning scheme that hierarchically subdivides a dataset into more homogeneous subsets, and then learning a piecewise constant function on that partition. This section serves to review several useful concepts for analyzing such nested partitioning rules that will be instrumental in our analysis.

For the first requirement, let us formalize the notion of the cell size in terms of the empirical measure induced by observations $\mathcal{X}=\{\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}\}$. For each cell $\Omega_{k}$, we define by

the proportion of observations falling inside $\Omega_{k}$. Throughout this work, we focus on partitions whose boxes can adaptively stretch and shrink, allowing the splits to arrange themselves in a data-dependent way . The simplest data-adaptive partition is based on statistically equivalent blocks , where all cells have approximately the same number of points, i.e. $\mu(\Omega_{k})\asymp 1/K$. We deviate from such a strict rule by allowing for imbalance and define the so called valid partitions.

(Valid Partitions) Denote by $\boldsymbol{\Omega}=\{\Omega_{k}\}_{k=1}^{K}$ a partition of $^{p}$. We say that $\boldsymbol{\Omega}$ is valid if

Valid partitions have non-empty cells, where the allocation does not need to be balanced. In balanced partitions (introduced in van der Pas and Ročková ), the cells satisfy $\frac{C_{min}^{2}}{K}\leq\mu(\Omega_{k})\leq\frac{C_{max}^{2}}{K}$ for some $C_{min}<1<C_{max}$. One prominent example of such balanced partitions is the median tree (or a $k$-$d$ tree) , which will be discussed in the next section and will be used as a benchmark tree approximation towards establishing Condition 2.2.

For the second requirement, let us formalize the notion of the cell size in terms of the local spread of the data. To this end, we introduce the partition diameter .

(Diameter) Denote by $\boldsymbol{\Omega}=\{\Omega_{k}\}_{k=1}^{K}$ a partition of $^{p}$ and by $\mathcal{X}=\{\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}\}$ a collection of data points in $^{p}$. For an index set $\mathcal{S}\subseteq\{1,\dots,p\}$, we define a diameter of $\Omega_{k}$ w.r.t. $\mathcal{S}$ as

The diameter of $\Omega_{k}$ corresponds to the largest $||\cdot||_{2}$ distance between $\mathcal{S}$-coordinate projections of two points that fall inside $\Omega_{k}$. This is one of the more flexible notions of a cell diameter, which takes into account the data itself, not just the physical cell size. Traditionally, bias and convergence rate analysis of tree-based estimators have been characterized in terms of the cell diameters. The rate depends on how fast the diameters shrink as we move down the tree: the more rapidly, the better. As will be seen later in Section 3.3, controlling the diameter will be essential for obtaining tight bounds on the approximation error.

2 Tree Partitions

We are ultimately interested in partitions that can be obtained with nested axis-parallel splits. Such partitions can be represented by a tree diagram, a hierarchical arrangement of nodes. We will focus on binary tree partitions, where each split yields two children nodes. Namely, starting from a parent node $^{p}$, a binary tree partition is grown by successively applying a splitting rule on a chosen internal node, say ${\Omega}_{k}^{\star}\subset^{p}$. The node ${\Omega}_{k}^{\star}$ is split into two cells by a perpendicular bisection along one of the $p$ coordinates, say $j$, at a chosen observed value $\tau\in\mathcal{X}=\{\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}\}$. The newborn cells $\{\boldsymbol{x}\in{\Omega}_{k}^{\star}:x_{j}\leq\tau\}$ and $\{\boldsymbol{x}\in{\Omega}_{k}^{\star}:x_{j}>\tau\}$ constitute two rectangular regions of $^{p}$, which can be split further (to become internal nodes) or end up being terminal nodes . The terminal cells after $K-1$ splits then yield a box-shaped (tree) partition $\{\Omega_{k}\}_{k=1}^{K}$.

Unlike dyadic trees, where the threshold $\tau$ is preset at a midpoint of the rectangle ${\Omega}_{k}$ along the $j^{th}$ direction, we allow for splits at available observations $\mathcal{X}$. Such data-dependent splits are integral to Bayesian CART and BART implementations . With more opportunities for each split, many more tree topologies can be generated. However, since the tree partitions are arranged in a nested fashion, their combinatorial complexity is not too large (as shown in Lemma 3.1 below).

We will denote each tree-structured $K$-partition by $\mathcal{T}^{K}=\{\Omega_{k}\}_{k=1}^{K}$. With $\mathcal{V}^{K}_{\mathcal{S}}$ we denote a family of valid tree partitions of $^{p}$, obtained by splitting $K-1$ times at observed values in $\mathcal{X}$ along each coordinate direction inside $\mathcal{S}\subseteq\{1,\dots,p\}$ at least once. In particular, each tree $\mathcal{T}\in\mathcal{V}^{K}_{\mathcal{S}}$ is valid according to Definition 3.1 and uses up all covariates in $\mathcal{S}$ for splits, where $\mathcal{S}$ can be regarded as an index set of active predictors. We will refer to the partitioning number $\Delta(\mathcal{V}^{K}_{\mathcal{S}})$ (in a similar vein as in ) as the overall number of distinct partitions of $\mathcal{X}$ that can be induced by members of the partitioning family $\mathcal{V}^{K}_{\mathcal{S}}$.

Denote by $\mathcal{S}\subseteq\{1,\dots,p\}$ an index set of active covariates. Let $\mathcal{V}^{K}_{\mathcal{S}}$ denote the set of valid tree partitions obtained with $K-1$ splits. Then

This follows from the recursive formula $\Delta(\mathcal{V}^{K}_{\mathcal{S}})=\Delta(\mathcal{V}^{K-1}_{\mathcal{S}})(n-K+2)|\mathcal{S}|$, where we use the fact that there are $\Delta(\mathcal{V}^{K-1}_{\mathcal{S}})$ possible trees with $K-1$ cells which have altogether $n-K+2$ potential next splits along one of the $|\mathcal{S}|$ coordinates. ∎

Lemma 3.1 will be fundamental for understanding the combinatorial richness of trees and for obtaining bounds on the covering numbers towards establishing Condition (2.1).

3 On Tree-structured Step Functions

The second critical ingredient in building a tree regressor is learning the piecewise function on a given partition. In this section, we describe some facts about the approximating properties of such tree-structured step functions (further referred to as trees). The understanding of how well we can approximate a smooth regression surface will be vital for establishing Condition (2.2).

For a family of valid tree partitions $\mathcal{V}^{K}_{\mathcal{S}}$, we denote by

The cell diameters oversee how closely one can approximate $f\in\mathcal{H}^{\alpha}_{p}\,\cap\,\mathcal{C}(\mathcal{S})$ with $f_{\mathcal{T},\boldsymbol{\beta}}\in\mathcal{F}(\mathcal{V}^{K}_{\mathcal{S}})$ and it is desirable that they decay fast with $K$. Ideally, the approximation error should be no slower than $q^{\gamma}/K^{1/q}$ for some $\gamma>0$, where $q=|\mathcal{S}|$. To get an intuitive insight into this requirement, consider the following perfectly regular partition: a $q$-dimensional “chess-board” that splits $^{q}$ into $K=s^{q}$ cubes of length $1/s=1/K^{1/q}$. The maximal interpoint distance inside each cube will be at most $\sqrt{q}/K^{1/q}$. This partition is, however, not adaptive and thereby less practical. It turns out that, under a mild requirement on the spread of the data points $\mathcal{X}$, the fast diameter decay is also guaranteed by the data-adaptive $k$-$d$ trees mentioned in Remark 3.1. The “mild requirement” is formalized in our definition below.

Denote by $\smash{\widehat{\mathcal{T}}}=\{\smash{\widehat{\Omega}}_{k}\}_{k=1}^{K}\in\mathcal{V}_{\mathcal{S}}^{K}$ the $k$-$d$ tree where $\mathcal{S}\subseteq\{1,\dots,p\}$ and $K=2^{s\,|\mathcal{S}|}$. We say that a dataset $\mathcal{X}$ is $(M,\mathcal{S})$-regular if

The definition states that in a regular dataset, the maximal diameter in the $k$-$d$ tree partition should not be much larger than a “typical” diameter. This condition ensures that, as more and more data points are collected, the data conform to a structure that does not permit outliers in active directions $\mathcal{S}$. For example, a fixed design on a regular grid (including directions $\mathcal{S}$) will satisfy (3.4). Our notion of regularity goes farther by allowing (a) the predictors to be correlated, (b) the points to scatter unevenly and/or close to a lower-dimensional manifold. The manifold, however, should be varying in active directions $\mathcal{S}$ and do so sufficiently smoothly (or be monotone) so that the cells in the $k$-$d$ tree do not contain isolated clouds of points. For example, data points arising from distributions with atomic marginals (in active directions) would violate regularity.

As will be seen in the following lemma, for regular datasets, $k$-$d$ trees have a faster diameter decay (roughly halving the cell diameters after one round of $q=|\mathcal{S}|$ splits), thereby attaining better approximation error. The following lemma is an important element of our proof skeleton, showing that there exists a tree (a $k$-$d$ tree) that approximates well.

Controlling the approximation error is only one of the critical aspects in our theoretical study. The second will be monitoring the complexity of our approximating function classes. Finding the right balance between the two will be instrumental for obtaining optimal performance.

Now that we have developed the necessary tools, we can dive into the main results of the paper.

Adaptive Dimension Reduction with Trees

In Regime 1, we assume that the target regression surface $f_{0}\in\mathcal{H}^{\alpha}_{p}\cap\mathcal{C}(\mathcal{S}_{0})$, although initially conceived as a function of $\boldsymbol{x}=(x_{1},\dots,x_{p})^{\prime}$, in fact depends only on a small fraction of $q_{0}$ features $\boldsymbol{x}_{\mathcal{S}_{0}}$, where $\mathcal{S}_{0}\subset\{1,\dots,p\}$. If an oracle were able to isolate $\mathcal{S}_{0}$, the $L_{2}$ minimax rate would improve from $n^{-\alpha/(2\alpha+p)}$ to $n^{-\alpha/(2\alpha+q_{0})}$ and it would be the fastest rate possible. When no such oracle information is available, characterized the minimax rate as follows: $\varepsilon_{n}^{2}=\lambda^{2}\left(\sqrt{n}\lambda\right)^{-4\alpha/(2\alpha+q_{0})}+\frac{q_{0}}{n}\log\left(\frac{p}{q_{0}}\right)$ where $\lambda$ is the isotropic Hölder norm. The first term is the classical minimax risk for estimating a $q_{0}$-dimensional smooth function. The second term is the penalty incurred by variable selection uncertainty. While the number of features $p$ is not forbidden from growing to infinity much faster than $n$, we keep in mind that consistent estimation will only be possible in sparse contexts where $q_{0}$ is small relative to $p$ and $n$ (in which case the complexity penalty will be dominated by the first term in the minimax rate).

Bayesian regression tree implementations that do not induce sparsity (when in fact present) are unfit for inference in high-dimensional setups. In particular, the traditional Bayesian CART prior grows trees by splitting each node, indexed by the depth level $d$, with a probability $g(d)=\gamma/(1+d)^{\alpha}$, where $\alpha>1,\gamma\in(0,1)$ are tuning parameters. The splitting variable is picked from $\{x_{1},\dots,x_{p}\}$ uniformly (i.e. with a fixed probability $\theta_{i}=1/p$). Recently, proposed an adaptive variant of this prior by placing a sparsity-inducing Dirichlet prior on the splitting proportions $\theta_{i}$. This prior uses fewer variables in the tree construction and thereby is more reluctant to overfit. In another popular Bayesian CART implementation, suggest directly placing a prior on $K$ and a conditionally uniform prior on tree topologies with $K$ bottom leaves. Again, in its original form, this prior will likely suffer from the curse of dimensionality, failing to harvest the intrinsic lower-dimensional structure. Here, we propose a fix to this problem. To make the Bayesian CART prior of appropriate for high-dimensional setups, we propose a spike-and-tree variant by injecting one more layer: a complexity prior over the active set of predictors.

Bayesian models for feature selection have traditionally involved a hierarchy of priors over subset sizes $q=|\mathcal{S}|$ and subsets $\mathcal{S}\subseteq\{1,\dots,p\}$ . Instead of modeling the mean outcome as a linear functional of active predictors $\boldsymbol{x}_{\mathcal{S}}$, here we grow a tree from $\boldsymbol{x}_{\mathcal{S}}$. We begin by treating $q_{0}$ as unknown with an exponentially decaying prior

Next, given the dimensionality $q$, we assume that all $p\choose q$ subsets $\mathcal{S}$ of $q=|\mathcal{S}|$ covariates are a-priori equally likely, i.e.

Given $q,\mathcal{S}$ and $K$, we assign a uniform prior over valid tree topologies $\mathcal{T}=\{\Omega_{k}\}_{k=1}^{K}\in\mathcal{V}^{K}_{\mathcal{S}}$, i.e.

Similar constraints on trees where each terminal node is assigned a minimal number of data points have been implemented in stochastic search algorithms . At the very least, we can choose $\bar{C}=1$ in (3.2) , merely requiring that the cells be non-empty. Finally, given the partition of size $K$, we assign an iid Gaussian prior on the step heights (similarly as in )

The name spike-and-tree prior deserves a bit of explanation. It follows from the fact that (T1) and (T2) will be satisfied if each covariate has a prior probability $\theta\sim\mathcal{B}(1,p^{u})$ of contributing to the mean regression surface for some $u>1$ . Endowing each covariate $x_{j}$ with a Bernoulli indicator $\gamma_{j}$, where $\Pi(\gamma_{j}=1\,|\>\theta)=\theta$, and building a tree on $\mathcal{S}=\{j:\gamma_{j}=1\}$, one obtains a mixture prior on $f(\boldsymbol{x})$ that pertains to spike-and-slab variable selection. Here, the slab is a tree prior built on active covariates rather than an independent product prior on active regression coefficients. This hierarchical construction has distinct advantages for variable selection. In linear regression, it is customary to select variables by thresholding marginal posterior inclusion probabilities $\Pi(\gamma_{j}=1\,|\>\boldsymbol{Y}^{(n)})$. These will be available also under our spike-and-tree construction. Inspecting the posterior probabilities $\Pi(\gamma_{j}=1\,|\>\boldsymbol{Y}^{(n)})$ obtained with our hierarchical tree prior will be a new avenue for conducting variable selection in Bayesian CART and BART, an alternative to . Thus, our prior has important practical implications for performing principled model-free variable selection.

2 Posterior Concentration for Bayesian CART

The difficulty in properly analyzing Bayesian CART stems from the combinatorial richness of the prior that makes it less tractable analytically. By building on our developments from previous sections, we are now fully equipped to present the first theoretical result concerning this method.

The following theorem solidifies the optimality properties of Bayesian CART by showing that under the hierarchical prior (T1)-(T5), the posterior adapts to both smoothness and sparsity, concentrating at the (near) minimax rate that depends only on the number of strong covariates regardless of how many noise variables are present. The near-minimaxity refers to an additional log factor. The result holds for sparse (high-dimensional) regimes, where $p$ can be potentially much larger than $n$ and where $q_{0}\lesssim\log^{1/2}n$. We will denote by

the collection of all tree-structured step functions (with various tree sizes and split subsets) that can be obtained by partitioning $\mathcal{X}$.

Assume $f_{0}\in\mathcal{H}^{\alpha}_{p}\cap\mathcal{C}(\mathcal{S}_{0})$ with $0<\alpha\leq 1$ and $0<q_{0}=|\mathcal{S}_{0}|$ such that $q_{0}\|f_{0}\|_{\mathcal{H}^{\alpha}}\lesssim\log^{1/2}n$ and $\|f_{0}\|_{\infty}\lesssim\log^{1/2}n$. Moreover, we assume that $\log p\lesssim n^{q_{0}/(2\alpha+q_{0})}$ and that $\mathcal{X}$ is $(M,\mathcal{S}_{0})$-regular. We endow $\mathcal{F}_{\mathcal{T}}$ with priors (T1)-(T5). Then with $\varepsilon_{n}=n^{-\alpha/(2\alpha+q_{0})}\log^{1/2}n$ we have

for any $M_{n}\to\infty$ in $P_{0}^{n}$-probability, as $n,p\to\infty$.

It is useful to note that Theorem 4.1 holds also when $p<n$ and when $q_{0}$ is fixed as $n\rightarrow\infty$. When $q_{0}$ is fixed, however, the assumptions $\|f_{0}\|_{\infty}\lesssim\log^{1/2}n$ and $q_{0}\|f_{0}\|_{\mathcal{H}^{\alpha}}\lesssim\log^{1/2}n$ can be omitted. The first assumption is needed here to make sure that the step sizes of an approximating $k$-$d$ tree are well behaved when $q_{0}\rightarrow\infty$. The result holds for a bit slower rate $\varepsilon_{n}=n^{-\alpha/(2\alpha+q_{0})}\log^{\beta}n$ with $\beta>1/2$ under slightly relaxed assumptions $q_{0}\|f_{0}\|_{\mathcal{H}^{\alpha}}\lesssim\log^{\beta}n$ and $\|f_{0}\|_{\infty}\lesssim\log^{\beta}n$.

The assumption of a regular design is an inevitable consequence of treating $\boldsymbol{x}$’s as fixed. As noted by in their study of random forests, point-wise consistency results have been complicated by the difficulty in controlling local (cell) diameters. The regularity assumption guarantees this control and is apt to be satisfied for most realizations of $\boldsymbol{x}$ from reasonable distributions on $^{p}$. The following Corollary certifies that Bayesian CART, under a suitable complexity prior on the number of terminal nodes, is reluctant to overfit. This is seen from the behavior of the posterior, which concentrates on values $K$ that are only a constant multiple larger than the optimal oracle value $n^{q_{0}/(2\alpha+q_{0})}$.

(Bayesian regression trees do not overfit.) Under the assumptions of Theorem 4.1 with $0<\alpha\leq 1$ we have

in $P^{n}_{0}$-probability for a suitable constant $C_{k}>0$.

Corollary 4.1 also reveals a fundamental limitation of trees (step functions) in recovering smoother functions. To see this, consider $f_{0}(x)=x$ which possesses Hölder smoothness $1$ and, by Corollary 4.1, will thus be approximated by trees with at most $n^{1/3}$ leaves (up to multiplicative constants) with high posterior probability. However, $f_{0}$ is also in $\mathcal{C}^{2}$ and the approximation error by a regular histogram with $n^{1/3}$ leaves will be at least of the order $n^{-1/3}$ which is too large to achieve the minimax rate of $n^{-5/2}$ over $\mathcal{C}^{2}$.

Tree Ensembles

Combining multiple trees through additive aggregation has proved to be remarkably effective for enhancing prediction . This section offers new theoretical insights into the mechanics behind the Bayesian variants of such tree ensemble methods. Our approach rests on a detailed analysis of the collective behavior of partitioning cells generated by individual trees. We will see that the overall performance is affected not only by the quality of single trees but also how well they can collaborate .

Additive regression trees grow an ensemble predictor by binding together $T$ tree-shaped regressors. For subsets $\boldsymbol{\mathcal{S}}=\{\mathcal{S}^{1},\dots,\mathcal{S}^{T}\}$ and tree sizes $\boldsymbol{K}=(K^{1},\dots,K^{T})^{\prime}$, we define a sum-of-trees model (forest) as

Sum-of-trees models offer an improved representation flexibility by chopping up the predictor space into more refined segmentations. These segmentations are obtained by superimposing multiple tree partitions, yielding what we define below as a global partition.

For a partition ensemble $\mathcal{E}$, we define a global partition

as the partition obtained by merging all cuts in $\mathcal{T}^{1},\dots,\mathcal{T}^{T}$. We refer to $\widetilde{\Omega}_{k}$’s as global cells in the ensemble.

The concept of the global partition can be better understood from Figure 1, where splits from $T=2$ trees (each having $3$ leaves) are merged to obtain a global partition with ${K}(\mathcal{E})=7$ global cells. Generally, the global partition $\mathcal{E}$ itself is not necessarily a tree and can have as many as $K(\mathcal{E})\leq\prod_{t=1}^{T}K^{t}$ cells. This upper bound can be attained if each trees splits $K^{t}-1$ times on a single variable, where each tree uses a different one. Without loss of generality, the global partition will be assumed to have non-empty cells. This requirement can be met by merging vacuous cells with their nonempty neighbors.

Bayesian additive regression trees were conceived as a collection of weak learners that capture different aspects of the predictor space . To characterize the amount of diversity/correlation between trees in the ensemble, we introduce the so-called stretching matrix.

For a partition ensemble $\mathcal{E}=\{\mathcal{T}^{1},\dots,\mathcal{T}^{T}\}$, we define the stretching matrix $\boldsymbol{A}(\mathcal{E})=(a_{ij})_{i,j}$ as follows: for each $1\leq i\leq K(\mathcal{E})$ and $1\leq j\leq T\bar{K}$ we have

where $1\leq t\leq T$ and $1\leq m\leq K^{t}$ are such that $j=\sum_{l=1}^{t-1}K^{l}+m$ and where $\Omega^{t}_{m}$ is the $m^{th}$ (local) cell in the $t^{th}$ tree $\mathcal{T}^{t}$.

Each row of the stretching matrix $\boldsymbol{A}(\mathcal{E})$ corresponds to one global cell and each column to one local cell. The row entries sum to $T$, indicating which local cells overlap with that global cell (as shown in Figure 2 for partitions from Figure 1). To further characterize the pattern of overlap between trees, we introduce the $(T\bar{K})\times(T\bar{K})$ Gram matrix

The off-diagonal elements measure the “similarity” between local cells, say $\Omega^{t}_{j}$ and $\Omega^{u}_{k}$, in terms of the number of global cells that they share. More formally, let $r[\widetilde{\mathcal{T}}(\mathcal{E}),V]=|\{\widetilde{\Omega}_{k}:\widetilde{\Omega}_{k}\cap V\neq\emptyset\}|$ be the restricted cell count , measuring the number of global cells that intersect with a compact set $V$. For $i=\sum_{s=1}^{t-1}K^{s}+j$ and $l=\sum_{s=1}^{u-1}K^{s}+k$ we can write $\widetilde{a}_{il}=r[\widetilde{\mathcal{T}}(\mathcal{E}),\Omega^{t}_{j}\cap\Omega^{u}_{k}]$. Small off-diagonal entries indicate less overlap, where the individual trees capture more diverse aspects of the predictor space. The diagonal elements, on the other hand, quantify the “persistence” of each local cell, say $\Omega_{k}^{u}$, counting the number of global cells it stretches over. More formally, for $i=\sum_{s=1}^{u-1}K^{s}+k$ we have $\widetilde{a}_{ii}=r[\widetilde{\mathcal{T}}(\mathcal{E}),\Omega_{k}^{u}]$. The amount of diversity (or tree dis-similarity) in the ensemble can be quantified with eigenvalues of $\widetilde{\boldsymbol{A}}(\mathcal{E})$. We denote by $\lambda_{min}(\mathcal{E})$ (resp. $\lambda_{max}(\mathcal{E})$) the minimal (resp. maximal) singular values of ${\boldsymbol{A}}(\mathcal{E})$ (i.e. square roots of extremal nonzero eigenvalues of $\widetilde{\boldsymbol{A}}(\mathcal{E})$). If some trees in the ensemble are redundant, the conditioning number $\kappa(\mathcal{E})\equiv\lambda_{max}(\mathcal{E})/\lambda_{min}(\mathcal{E})$ will be large. The idea of diversifying trees was originally introduced by Breiman via subsampling. One could, in principle, impose a restriction on $\lambda_{min}(\mathcal{E})$ in the prior to encourage the trees to collaborate and diversify. However, this is not required for our theoretical study. We will focus on the so-called valid ensembles which consist of valid trees.

An ensemble $\mathcal{E}=\{\mathcal{T}^{1},\dots,\mathcal{T}^{T}\}$ is valid if each $\mathcal{T}^{t}$ is valid according to Definition 3.1. For tree sizes ${\boldsymbol{K}}=(K^{1},\dots,K^{T})^{\prime}$ and subsets $\boldsymbol{\mathcal{S}}=\{\mathcal{S}^{1},\dots,\mathcal{S}^{T}\}$, we denote the set of all valid ensembles by $\mathcal{V}\mathcal{E}_{\boldsymbol{\mathcal{S}}}^{\boldsymbol{K}}$.

The representation flexibility of additive trees also pertains to jump sizes. The global step size coefficients under additive trees are intertwined due to the tree overlap. This can be seen from the following, more compact, representation of (5.1):

where $\boldsymbol{A}(\mathcal{E})$ is the stretching matrix defined in (5.2). This link unfolds the theoretical analysis of tree ensembles using tools that we have already developed for single trees. Note that the condition number $\kappa(\mathcal{E})$ determines how much the relative change in $\bar{\boldsymbol{\beta}}$ influences the relative change in $\mathcal{B}$.

The mapping (5.4) can be in principle many-to-one in the sense that many tree-structured step functions can sum towards the same target (5.1). Such over-parametrization occurs, for instance, when ${K}(\mathcal{E})<T\bar{K}$ or, more generally, when ${\boldsymbol{A}}(\mathcal{E})$ has zero eigenvalues. This redundancy is not entirely unwanted and, in fact, it endows sum-of-trees models with a lot of modeling freedom.

We now formally define the space of approximating additive trees. For variable sets $\boldsymbol{\mathcal{S}}=\{\mathcal{S}^{1},\dots,\mathcal{S}^{T}\}$ and a vector of tree sizes $\boldsymbol{K}=(K^{1},\dots,K^{T})^{\prime}$, we denote by

the set of all additive tree step functions supported on valid ensembles $\mathcal{V}\mathcal{E}_{\boldsymbol{\mathcal{S}}}^{\boldsymbol{K}}$. The union of these over the number of trees $T$, all possible sets $\boldsymbol{\mathcal{S}}$ of sizes $\boldsymbol{q}=(q^{1},\dots,q^{T})^{\prime}$ and tree sizes $\boldsymbol{K}$ gives rise to

our approximating space of additive regression tree functions.

2 Additive Regression Trees are Adaptive

This section provides an interesting initial perspective on the behavior of Bayesian additive regression trees in Regime 1. We will continue with the more general Regime 2 in the next section. We will focus on a variant of the popular Bayesian Additive Regression Trees (BART) model , modified in three ways. First, the tree prior will be according to rather than . The second modification is that the trees are built on the same set of variables, endowed with a subset selection prior construction. In the next section, we allow for the fully general case where each tree builds on a potentially different set of variables. Third, rather than fixing the number of trees, we endow $T$ with a prior distribution.

We will see that having a good control of the regression function variation inside each global cell together with a good choice of the prior on the total number of leaves $\sum K^{t}$ will be sufficient to ensure optimal behavior. The approximation ability of tree ensembles hinges on the diameter of the global partition. Each tree partition does not need to have a small diameter (i.e. can be a weak learner), as long as the global one does. An important building block in our proof will be the construction of a single tree ensemble that can approximate well. As will be shown in Lemma 10.1, we can construct such ensemble by first finding a single $k$-$d$ tree from Lemma 3.2 (a strong learner) and then redistributing the cuts among small trees (weak learners) in a way that the global partition is exactly equal to the $k$-$d$ tree. An example of this deconstruction is depicted in Figure 3, where a full symmetric tree from Figure 4, say $\widehat{\mathcal{T}}$, has been trimmed into many smaller imbalanced trees which add up towards $\widehat{\mathcal{T}}$. More details on this decomposition are in the proof of Lemma 10.1.

The following theorem is an ensemble variant of Theorem 4.1 which will serve as a useful stepping stone towards the full-fledged result presented in the next section. Instead of approximating $f_{0}\in\mathcal{H}^{\alpha}_{p}\cap\mathcal{C}(\mathcal{S}_{0})$ with one large tree in Regime 1, we build a forest made up of many smaller trees (weak learners). The first two layers of the prior (T1) and (T2) are the same. Next, we assign a prior on the number of trees

which is sufficiently diffuse so as to promote ensembles with many trees. Next, conditionally on $T$, we assign a prior on $\boldsymbol{K}=(K^{1},\dots,K^{T})^{\prime}$ as an independent product of Poisson distributions (T3). An important distinction between the Poisson prior deployed for single trees in Section 4 and the one deployed here for ensembles is that the hyper-parameter depends on $T$. Namely, our prior on leaves satisfies

As a prior on $\mathcal{E}$, given $(\mathcal{S},\boldsymbol{K})$, we use the uniform prior over valid ensembles

where $\Delta(\mathcal{V}\mathcal{E}_{\mathcal{S}}^{\boldsymbol{K}})$ is the overall number of ensembles consisting of $T$ valid trees that can be obtained by splitting on data points $\mathcal{X}$. Recall that throughout this section, all trees are constrained to use split variables in the same set $\mathcal{S}$. To mark this difference, we have denoted the partition ensembles with $\mathcal{V}\mathcal{E}_{\mathcal{S}}^{\boldsymbol{K}}$ instead of $\mathcal{V}\mathcal{E}_{\boldsymbol{\mathcal{S}}}^{\boldsymbol{K}}$.

The following theorem shows that additive regression trees, in combination with a subset selection prior, can nicely adapt to the ambient dimensionality and smoothness, also achieving the optimal concentration rate in Regime $1$.

Assume $f_{0}\in\mathcal{H}^{\alpha}_{p}\cap\mathcal{C}(\mathcal{S}_{0})$ with $0<\alpha\leq 1$ and $0<q_{0}=|\mathcal{S}_{0}|$ such that $q_{0}\|f_{0}\|_{\mathcal{H}^{\alpha}}\lesssim\log n$ and $\|f_{0}\|_{\infty}\lesssim\log n$. Moreover, we assume $\log p\lesssim n^{q_{0}/(2\alpha+q_{0})}$ and that $\mathcal{X}$ is $(M,\mathcal{S}_{0})$-regular. We endow $\mathcal{F}_{\mathcal{E}}$ with priors (T1),(T2),(T3*)-(T6*). With $\varepsilon_{n}=n^{-\alpha/(2\alpha+q_{0})}\log n$ we have

for any $M_{n}\to\infty$ in $P_{0}^{n}$-probability, as $n,p\to\infty$.

Similarly as for single trees, we obtain the following corollary which states that the posterior concentrates on ensembles whose overall number of leaves is not much larger than the optimal value $n^{q_{0}/(2\alpha+q_{0})}\log n$.

Under the assumptions of Theorem 5.1 with $0<\alpha\leq 1$ we have

in $P_{0}^{n}$-probability, as $n,p\to\infty$, for a suitable constant $C_{k}>0$.

This follows from Lemma 1 of and Section 10.3. ∎

Corollary 5.1 shows that the posterior distribution rewards either many weak learners or a few strong ones. The compromise between the two is regulated by the prior (T3*), where stronger shrinkage (i.e. larger $C_{T}$) will result in fewer trees. This corollary provides an important theoretical justification for why Bayesian additive regression trees have been so resilient to overfitting in practice.

The dependence on $T$ in the Poisson prior (T4⋆) works nicely in tandem with the exponential prior (T3*). Removing $T$ from (T4⋆) would have to be balanced with a bit stronger prior $\pi(T)$. Theorem 5.1 also holds with $T$ fixed (with slight modifications of the proof). The prior $\pi(T)$ will be instrumental in the additive case (next section). The dependence on $T$ in (T6*), though recommended in practice , is not needed in Theorem 5.1.

Tree Ensembles in Additive Regression

In Section 4.2 we have shown that the posterior distribution under the Bayesian CART prior has optimal properties. However, it is now well known that the practical deployments of Bayesian CART suffer from poor MCMC mixing. Additive aggregations of single small trees have proven to have far more superior mixing properties. One may wonder whether the benefits of additive trees are purely computational or whether there are some aspects that make them more attractive also theoretically. We will address this fundamental question.

For estimating a single smooth function, we were not able to tell apart single trees from tree ensemble in terms of their convergence rate (besides perhaps a small difference in the log factor). They are both optimal in Regime 1. Tree ensembles are inherently additive and, as such, are well-equipped for approximating additive $f_{0}$ (Regime 2). Throughout this section, we assume

where $f_{0}^{t}\in\mathcal{H}^{\alpha^{t}}_{p}\cap\mathcal{C}(\mathcal{S}_{0}^{t})$. Note that each component $f_{0}^{t}$ depends only on a potentially very small subset $\mathcal{S}_{0}^{t}$ of covariates, where $|\mathcal{S}_{0}^{t}|=q_{0}^{t}$. However, the additive structure allows $f_{0}$ to depend on a larger number of variables, say $q_{0}$, where $\max_{1\leq t\leq{T_{0}}}q_{0}^{t}\leq q_{0}\leq\sum_{t=1}^{T_{0}}q_{0}^{t}$. The minimax rate $r_{n}^{2}$ for estimating $f_{0}$ in Regime 2 satisfies $C_{1}\underline{\varepsilon}_{n}^{2}\leq r_{n}^{2}\leq C_{2}\bar{\varepsilon}_{n}^{2}$ , where

and where $\lambda^{t}$ is the Hölder norm of $f_{0}^{t}$. The sparsity constraint in Regime 2 is less strict than in Regime 1, where $q_{0}$ can be potentially larger than $\log n$ while still allowing for consistent estimation. For the isotropic case ($\alpha^{t}=\alpha$ and $q_{0}^{t}=\widetilde{q}_{0}$), single trees can achieve the slower rate $n^{-2\alpha/(2\alpha+q_{0})}$ where $\widetilde{q}_{0}\leq q_{0}\leq T_{0}\widetilde{q}_{0}$. As will be shown below, tree ensembles can achieve a faster rate $\varepsilon_{n}^{2}=\sum_{t=1}^{T_{0}}(\varepsilon_{n}^{t})^{2}$ where $\underline{\varepsilon}_{n}^{2}\lesssim\varepsilon_{n}^{2}\lesssim\bar{\varepsilon}_{n}^{2}$ and where $(\varepsilon_{n}^{t})^{2}=\lambda^{t2}\left(\sqrt{n}\lambda^{t}\right)^{-4\alpha^{t}/(2\alpha^{t}+q_{0}^{t})}+\frac{q_{0}^{t}}{n}\log\left(\frac{p}{q_{0}^{t}}\right)$.

We will approximate $f_{0}$ with tree ensembles $f_{\mathcal{E},\mathcal{B}}\in\mathcal{F}_{\mathcal{E}}$. The ensembles here differ from the ones considered in Section 5.2. The crucial difference is that now we allow each of the trees $\mathcal{T}^{t}$ to depend on a different set of variables $\mathcal{S}^{t}$. Now we have a vector of subset sizes $\boldsymbol{q}=(q^{1},\dots,q^{T})^{\prime}$ and a set of subsets $\boldsymbol{\mathcal{S}}=\{\mathcal{S}^{1},\dots,\mathcal{S}^{T}\}$, one for each tree. We consider the following independent product variant of the complexity prior (T1)

and a product prior variant of (T2), given $(T,\boldsymbol{q})$,

The prior on the number of trees, the number of leaves, ensembles and step sizes is the same as in (T3*), (T4*), (T5*).

We are now ready to present our final result showing that the posterior concentration for Bayesian additive regression trees is near-minimax rate optimal when $f_{0}$ has an additive structure.

Assume that $f_{0}$ is as in (6.1) where $f_{0}^{t}\in\mathcal{H}^{\alpha^{t}}_{p}\,\cap\,\mathcal{C}(\mathcal{S}_{0}^{t})$ with $0<\alpha^{t}\leq 1$ and $q_{0}^{t}=|\mathcal{S}_{0}^{t}|$ such that $0<q_{0}^{t}\|f_{0}^{t}\|_{\mathcal{H}_{\alpha^{t}}}\lesssim\log^{1/2}n$ and $\|f_{0}^{t}\|_{\infty}\lesssim\log^{1/2}n$. Moreover, assume $\log p\lesssim\min\limits_{1\leq t\leq T_{0}}n^{q_{0}^{t}/(2\alpha^{t}+q_{0}^{t})}$ and that $\mathcal{X}$ is $(M,\mathcal{S}_{0}^{t})$-regular for $1\leq t\leq T_{0}$. We endow $\mathcal{F}_{\mathcal{E}}$ with priors (T1*)-(T6*). With $\varepsilon_{n}^{2}=\sum_{t=1}^{T_{0}}n^{-2\alpha^{t}/(2\alpha^{t}+q_{0}^{t})}\log n$ and $T_{0}\lesssim n$, we have

for any $M_{n}\to\infty$ in $P_{0}^{n}$-probability, as $n,p\to\infty$ .

Under the assumptions $\log p\lesssim n^{q_{0}^{t}/(2\alpha^{t}+q_{0}^{t})}$ and $q_{0}^{t}\|f_{0}^{t}\|_{\mathcal{H}_{\alpha^{t}}}\lesssim\log^{1/2}n$, the second term $\frac{q_{0}^{t}}{n}\log\left(\frac{p}{q_{0}^{t}}\right)$ in $(\varepsilon_{n}^{t})^{2}$ (defined earlier) is dominated by the first term $n^{-2\alpha^{t}/(2\alpha^{t}+q_{0}^{t})}\log n$. This is why the second term does not appear in the rate $\varepsilon_{n}^{2}$ in Theorem 6.1.

Failing to recognize the additive structure in $f_{0}$, single regression trees achieve the slower rate $n^{-\alpha/(2\alpha+q_{0})}\log^{1/2}n$, according to Theorem 4.1. Theorem 6.1 thus provides an additional theoretical justification for Bayesian additive tree models suggesting their performance superiority over single trees when $f_{0}$ is additive.

Our priors differ from the widely used BART implementations in three ways: (1) we focus on the uniform prior of Denison et al. , (2) we assign a prior distribution on the number of trees and (c) we deploy the spike-and-slab wrapper. Implementations of our priors are feasible with some modifications of the existing software. For the Bayesian CART prior that we analyze, Denison et al. propose a reversible jump MCMC implementation. While this algorithm is different from BART, the acceptance ratios in the Metropolis-Hastings step differ only very slightly. Liu, Ročková and Wang extended their sampler to the spike-and-tree (spike-and-forest) versions in two ways. The first one is a Metropolis-Hasting strategy that consists of joint sampling from variable subsets as well as trees (forests). As a faster alternative, they proposed an approximate ABC sampling strategy based on data splitting (called ABC Bayesian Forests).

Discussion

In this work, we have laid down foundations for the theoretical study of Bayesian regression trees and their additive variants. We have shown an optimal behavior of Bayesian CART, the first theoretical result on this method. We have developed several useful tools for analyzing additive regression trees (variants of the BART method), showing their optimal performance in both additive and non-additive regression. The smoothness order of studied functions is restricted to values not exceeding one, a main limitation of our approach due to the fact that our approximations are piecewise constants . While in the one-dimensional case, step functions are not appealing estimators of a regression function that is thought to be smooth, methods like CART and BART are attractive and feasible solutions in complex high-dimensional data. The limitation $\alpha\leq 1$ could be overcome by extending our approach to piecewise polynomials or kernels, an elaboration that we leave for future investigation. One such extension was recently proposed in a related paper by Linero and Yang . These authors obtained concentration results for a kernel method that can be regarded as a smooth variant of BART. The results of do not apply for single trees, only aggregates of kernels. In contrast, we study actual posteriors of single trees, as well as forests, and analyze sieves of step functions which are the essence of the actual BART method.

While our priors do not exactly match the BART prior, BART could be adapted to achieve the same optimality properties. The first modification is the splitting probability (as pointed out in Remark 4.3). The second modification is the spike-and-slab wrapper (as pointed out in Section 4.1) or a modification of the prior on the split variables (as in ). The prior distribution on the number of trees will only be beneficial in the additive model and is not needed when $f_{0}$ has one layer.

The assumption of a known $\sigma_{0}=\mathsf{Var}\,(\varepsilon_{i})=1$ can be relaxed. It has been noted in the literature (e.g. ) that the general result of Ghosal and van der Vaart (which we build upon) can be extended to the unknown $\sigma_{0}$ case. Such an extension was formally proved in Jonge and van Zanten (2013), who assume that $\sigma_{0}$ belongs to a compact interval $[a,b]$ and the prior $\pi(\sigma)$ concentrates on $[a,b]$. We could obtain our results under this restriction as well by verifying suitably adapted conditions (2.1)-(2.3). In related work, Yoo and Ghosal show optimal posterior concentration (in both $L_{2}$ and $L_{\infty}$ sense) in non-parametric regression with unknown variance and B-spline tensor product priors. Next, they show that under an inverse-gamma prior, the posterior for $\sigma$ contracts at $\sigma_{0}$ at the same rate. Moreover, for any prior on $\sigma$ with positive and continuous density, the posterior of $\sigma$ is consistent. These results are obtained under the assumption that $f_{0}$ is uniformly bounded. We anticipate that similar results will hold also for our priors when $f_{0}$ is uniformly bounded.

Acknowledgement We are grateful to the Associate Editor and two referees for helpful comments, and to Johannes Schmidt-Hieber for helpful discussions on the assumption $\alpha\leq 1$.

Proof of Theorem 4.1

where $\mathcal{F}(\mathcal{V}_{\mathcal{S}}^{K})$ was defined in (3.3). The optimal choice of $k_{n}$ and $q_{n}$ will follow from our considerations below.

We start with a useful lemma that characterizes a useful upper bound on the covering number of the smaller sets $\mathcal{F}(\mathcal{V}_{\mathcal{S}}^{K})$.

Let $\mathcal{F}(\mathcal{V}^{K}_{\mathcal{S}})$ be the class of step functions (3.3). Then

where $\Delta(\mathcal{V}^{K}_{\mathcal{S}})$ is the partitioning number of $\mathcal{V}^{K}_{\mathcal{S}}$.

Denote by $f_{\mathcal{T},\smash{\widehat{\boldsymbol{\beta}}}}$ the $\|\cdot\|_{n}$ projection of $f_{0}$ onto $\mathcal{F}(\mathcal{T})\subset\mathcal{F}(\mathcal{V}_{\mathcal{S}}^{K})$, the set of all step functions that live on a given partition $\mathcal{T}$. Then $\{\boldsymbol{\beta}:\|f_{\mathcal{T},\boldsymbol{\beta}}-f_{0}\|_{n}\leq\varepsilon\}\subset\{\boldsymbol{\beta}:\|\boldsymbol{\beta}-\smash{\widehat{\boldsymbol{\beta}}}\|_{2}\leq\varepsilon\sqrt{n}/\bar{C}\}$ and $\{\boldsymbol{\beta}:\|f_{\mathcal{T},\boldsymbol{\beta}}-f_{0}\|_{n}\leq\varepsilon/36\}\supset\{\boldsymbol{\beta}:\|\boldsymbol{\beta}-\smash{\widehat{\boldsymbol{\beta}}}\|_{2}\leq\varepsilon/36)\}$. This relationship shows that the $\varepsilon/36$ covering number of an $||\cdot||_{n}$ ball $\{\boldsymbol{\beta}:\|f_{\mathcal{T},\boldsymbol{\beta}}-f_{0}\|_{n}\leq\varepsilon\}$ can be bounded from above by the $\varepsilon/36$ covering number of an Euclidean ball of a radius $\varepsilon\sqrt{n}/\bar{C}$, which is bounded by $\left(\frac{108}{\bar{C}}\sqrt{n}\right)^{K}$. We can repeat this argument by projecting $f_{0}$ onto any valid tree topology $\mathcal{T}\in\mathcal{V}^{K}_{\mathcal{S}}$. The number of such valid trees is no larger than $\Delta(\mathcal{V}^{K}_{\mathcal{S}})$, which completes the proof.∎

The covering number for the entire sieve $\mathcal{F}_{\mathcal{T}}^{n}$ is then seen to satisfy

From Lemma 8.1 and Lemma 3.1, we obtain the following upper bound

where we used the fact $K!<K^{K}$. Next, using the regularized incomplete beta function representation of the Binomial cdf, we can write

where $C=\left(\frac{108}{\bar{C}}\right)$. Finally, the entropy condition requires that the log-covering number, now upper bounded by

2 Condition (2.2)

We wish to show that the prior assigns enough mass around the truth in the sense that

for some $C_{1}>0$. The proof Lemma 3.2 is in the Supplemental Material (Section 3).

To continue with the proof of Theorem 4.1, We find the smallest $K=2^{s\,q_{0}}$ such that the function $f_{\widetilde{\mathcal{T}},\widetilde{\boldsymbol{\beta}}}\in\mathcal{F}(\mathcal{V}_{\mathcal{S}_{0}}^{K})$ in (8.5) safely approximates $f_{0}$ with an error that is no larger than $\varepsilon_{n}/2$, a constant multiple of the target rate. Such a $K$ will be denoted by $a_{n}$ and is defined as the smallest $K$ such that $\varepsilon_{n}\geq 2C_{0}q_{0}K^{-\alpha/q_{0}}$ for $C_{0}=||f_{0}||_{\mathcal{H}^{\alpha}}C_{1}$. Then we have

Then the statement $\|\boldsymbol{\beta}-\smash{\widehat{\boldsymbol{\beta}}}\|_{2}<\varepsilon_{n}/2$ implies $\|f_{0}-f_{\smash{\widehat{\mathcal{T}}},\boldsymbol{\beta}}\|_{n}<\|f_{0}-f_{\smash{\widehat{\mathcal{T}}},\smash{\widehat{\boldsymbol{\beta}}}}\|_{n}+\varepsilon/2<\varepsilon,$ where the last inequality follows from the definition of $a_{n}$. Thus, we have

From the triangle inequality (and because $\smash{\widehat{\mathcal{T}}}$ is balanced) we have

Taking minus the log of this quantity, Condition (2.2) will be met when

is smaller than a constant multiple of $n\varepsilon_{n}^{2}$. Above, we omitted the small terms ${\varepsilon_{n}^{2}/8}$ (since $\varepsilon_{n}\to 0$) and $\log(a_{n}/2)$. First, we note that

With $q_{0}\lesssim\log^{\beta}n$ and $\log p\lesssim n^{q_{0}/(2\alpha+q_{0})}$, we obtain $-\log c(p,q_{0},a_{n})\lesssim n\varepsilon_{n}^{2}$. Next, focusing on the last term in (8.11), we obtain (from the left inequality in (8.6)) the following bound

and hence $q_{0}/\varepsilon_{n}\lesssim a_{n}$ for $\alpha\leq 1$ and $q_{0}\geq 1$. Moreover, from the right inequality in (8.6) we obtain for $\varepsilon_{n}=n^{-\alpha/(2\alpha+q_{0})}\log^{\beta}n$ and $2C_{0}q_{0}\lesssim\log^{\beta}n$

Under our assumption $\|f_{0}\|_{\infty}\lesssim\log^{\beta}n$, (8.12) immediately yields $a_{n}\|f_{0}\|_{\infty}^{2}\lesssim n\varepsilon_{n}^{2}$. All of these considerations, combined with the fact $-\log{\pi(a_{n})}\lesssim a_{n}\log a_{n}$, yield the following leading term behind the last three summands in (8.11): $a_{n}\log(a_{n}^{3/2}n)$. Using (8.12), we obtain $a_{n}\log(a_{n}^{3/2}n)\lesssim a_{n}\log n\lesssim n\varepsilon_{n}^{2}$ for $\beta\geq 1/2$. Altogether, there exists $d>0$ such that (8.4) is satisfied.

3 Condition (2.3)

for a large enough constant $C_{q}>0$. Indeed, for $p>n$ we have $q_{n}\log p\asymp n\,\varepsilon_{n}^{2}$. For $p<n\,\varepsilon_{n}^{2}/\log n$ we have $q_{n}=p$ and $\Pi(q>q_{n})=0$. Finally, for $n\,\varepsilon_{n}^{2}/\log n\leq p\leq n$, we have $q_{n}\log n\asymp n\,\varepsilon_{n}^{2}$ and $\log p\geq\frac{q_{0}}{2\alpha+q_{0}}\log n+(2\beta-1)\log\log n$. For $c>1$ and $\beta\geq 1/2$, we can write $q_{n}[\log c+a\log p]\leq C_{q}\frac{aq_{0}}{2\alpha+q_{0}}n^{q_{0}/(2\alpha+q_{0})}\log^{2\beta}n$ and (8.13) holds for $C_{q}$ large enough. Next, we apply the Chernoff bound for $\Pi(K>k_{n})$. Namely, for any $t>0$ we can write

With our choice $k_{n}=\lfloor C_{k}n\varepsilon_{n}^{2}/\log n\rfloor\asymp C_{k}n^{q_{0}/(2\alpha+q_{0})}\log^{2\beta-1}n$ (Section 8.1) and with $t=\log k_{n}$ we obtain

Proof of Theorem 6.1

We aim to establish conditions (2.1), (2.2) and (2.3) for $\varepsilon_{n}^{2}=\sum_{t=1}^{T_{0}}({\varepsilon}_{n}^{t})^{2}$, where ${\varepsilon}_{n}^{t}=n^{-\alpha^{t}/(2\alpha^{t}+q_{0}^{t})}\log^{\beta^{t}}n$. Our sieve consists of valid forests with either (a) many trees that are small (weak learners), or (b) a few large trees (strong learners). We impose a joint requirement on $\sum_{t=1}^{T}K^{t}$ so that the overall number of leaves in the ensemble is small. At the same time, we require that $\sum_{t=1}^{T}q^{t}$ (the upper bound on the number of active variables in the ensemble) is small as well. The sieve is constructed as follows:

for some integer values $s_{n}$ and $z_{n}$. Throughout this section we denote $\bar{K}=\frac{1}{T}\sum_{t=1}^{T}K^{t}$.

We first obtain the following upper bound on the log-covering number

Now we find an upper bound on the number of valid ensembles $\mathcal{E}\in\mathcal{V}\mathcal{E}_{\boldsymbol{\mathcal{S}}}^{\boldsymbol{K}}$ inside the sieve $\mathcal{F}_{\mathcal{E}}^{n}$. To start, we note that given $(T,\boldsymbol{q},\boldsymbol{K},\boldsymbol{\mathcal{S}})$, there are at most $\prod_{t=1}^{T}(K^{t}\,q^{t}\,n)^{K^{t}}$ valid ensembles $\mathcal{V}\mathcal{E}_{\boldsymbol{\mathcal{S}}}^{\boldsymbol{K}}$. This bound is obtained from Lemma 3.1 by combining all possible $T$-tuples of trees.The order of trees in $\mathcal{E}$ matters. Given $(T,\boldsymbol{q})$, there are $\prod_{t=1}^{T}{p\choose q^{t}}$ sets of subsets $\boldsymbol{\mathcal{S}}=\{\mathcal{S}^{1},\dots,\mathcal{S}^{T}\}$ satisfying the constraint $|\mathcal{S}^{t}|=q^{t}$. This leads to an overall upper bound

Combining this bound with (9.2), we obtain the following bound

Condition 2.1 will be met when (9.3) is smaller than (a constant multiple of) $n\varepsilon_{n}^{2}=\sum_{t=1}^{T_{0}}n(\varepsilon_{n}^{t})^{2}$. With the choice $z_{n}=\lfloor C_{z}n\varepsilon_{n}^{2}/\log n\rfloor$ and $s_{n}=\lceil C_{s}n^{q_{0}/(2\alpha+q_{0})}\log^{2\beta}n/\log(p\vee n)\rceil$, where $C_{s}$ and $C_{z}$ are large enough constants to be determined later, this condition is satisfied.

2 Condition (2.2)

To establish Condition (2.2) for tree ensembles, we begin by finding a single additive tree that approximates well. We will heavily leverage our findings from Section 8.2, noting that the problem of approximating an additive function $f_{0}$ with a sum of trees can be decomposed into smaller problems of approximating each layer $f_{0}^{t}$ separately.

Denote by $a_{n}^{t}$ the smallest leaf size of a $k$-$d$ tree (defined in Remark 3.1) needed to approximate $f_{0}^{t}$ with an error smaller than $\varepsilon_{n}^{t}/2$, where $\varepsilon_{n}^{t}=n^{-\alpha^{t}/{(2\alpha^{t}+q^{t}_{0})}}\log^{\beta^{t}}n$. Such a tree step function approximation exists according to Lemma 3.2 when $\mathcal{X}$ is $(M,\mathcal{S}_{0}^{t})$-regular. We will denote this approximation with $f_{\smash{\widehat{\mathcal{T}}}^{t},{\smash{\widehat{\boldsymbol{\beta}}}}^{t}}(\boldsymbol{x})$. Moreover, with $\boldsymbol{a}_{n}=(a_{n}^{1},\dots,a_{n}^{T_{0}})^{\prime}$ we denote the vector of such minimal tree sizes, where each $a_{n}^{t}$ satisfies (8.6) with $q_{0}^{t},\alpha^{t}$ and $\varepsilon_{n}^{t}$. Next, we will denote by $\smash{\widehat{\mathcal{E}}}=\{\smash{\widehat{\mathcal{T}}}^{1},\dots,\smash{\widehat{\mathcal{T}}}^{T_{0}}\}$ the approximating partition ensemble with step heights $\smash{\widehat{\mathcal{B}}}=(\smash{\widehat{\boldsymbol{\beta}}}^{1\prime},\dots,\smash{\widehat{\boldsymbol{\beta}}}^{T_{0}\prime})^{\prime}$. The individual tree approximations $f_{\smash{\widehat{\mathcal{T}}}^{t},\smash{\widehat{\boldsymbol{\beta}}}^{t}}(\boldsymbol{x})$ are woven into an approximating forest $f_{\smash{\widehat{\mathcal{E}}},\smash{\widehat{\mathcal{B}}}}(\boldsymbol{x})=\sum_{t=1}^{T_{0}}f_{\smash{\widehat{\mathcal{T}}}^{t},\smash{\widehat{\boldsymbol{\beta}}}^{t}}(\boldsymbol{x})\in\mathcal{F}(\mathcal{V}\mathcal{E}_{\boldsymbol{\mathcal{S}}_{0}}^{\boldsymbol{a}_{n}})$, where $\boldsymbol{\mathcal{S}}_{0}=\{\mathcal{S}^{1}_{0},\dots,\mathcal{S}^{T_{0}}_{0}\}$.

Because we assumed $\beta_{j}^{t}\sim\mathcal{N}(0,1/T)$, given $T$, we can directly use (8.9) to lower-bound the above with

Because each tree $\smash{\widehat{\mathcal{T}}}^{t}$ is a $k$-$d$ tree and is by definition balanced, we have $\|\smash{\widehat{\boldsymbol{\beta}}}^{t}\|_{n}^{2}\lesssim a_{n}^{t}\|f_{0}^{t}\|_{\infty}^{2}$. Now we can directly apply all our calculations from Section 8.2. In particular, using (9.4) and noting that $\Delta(\mathcal{V}\mathcal{E}_{\boldsymbol{\mathcal{S}}_{0}}^{\boldsymbol{a}_{n}})<\prod_{t=1}^{T_{0}}\Delta(\mathcal{V}_{\mathcal{S}_{0}^{t}}^{a_{n}^{t}})$, we obtain

where $L(\cdot)$ was defined in (8.8). It follows from Section (8.2) that

3 Condition (2.3)

for any $\gamma>0$, where we used the fact

With $\gamma=\log(p\vee n)$ and $a>2$, we can write