\begin{document}

\maketitle

\begin{abstract}
\input{abstract}


\end{abstract}

\section{Introduction}
Decision trees are the building blocks of some of the most important and powerful algorithms in statistical learning. For example, ensembles of decision trees are used for some bootstrap aggregated prediction rules (e.g., bagging \citep{breiman1996bagging} and random forests \citep{breiman2001}). In addition, each iteration of gradient tree boosting (e.g., TreeBoost \citep{friedman2001greedy}) fits the pseudo-residuals with decision trees as base learners. From an applied perspective, decision trees have an appealing interpretability and are accompanied by a rich set of analytic and visual diagnostic tools. These attributes make tree-based learning particularly well-suited for applied sciences and related disciplines---which may rely heavily on understanding and interpreting the output from

the training algorithm.

Although, as with many aspects of statistical learning, good empirical performance often comes at the expense of rigor.

Tree-structured learning with decision trees is no exception---statistical guarantees for popular variants, i.e., those that are actually used in practice, are hard to find. Indeed, the recursive manner in which decision trees are constructed makes them unamenable to analysis, especially when the split protocol involves \emph{both} the input and output data. Despite these challenges, we take a step forward in advancing the theory of decision trees and aim to tackle the following fundamental question:
\begin{quotation}
\centering\normalsize\textit{When do decision trees adapt to the sparsity of a predictive model?}
\end{quotation}


To make our work informative to the applied user of decision trees, we strive to make the least departure from practice and therefore focus specifically on Classification and Regression Tree (CART) \citep{breiman1984} methodology---by far the most popular for regression and classification problems. With this methodology, the tree construction importantly depends on both the input and output data and is therefore \emph{data dependent}. This aspect lends itself favorably to the empirical performance of CART, but poses unique mathematical challenges. It is perhaps not surprising then that, despite the widespread use of CART, there have been only a small handful of papers that study its theoretical properties. For example, \citep{scornet2015} study the asymptotic properties of CART in a fixed dimensional regime, en route to establishing consistency of Breiman's random forests for additive regression models. Another notable paper \citep{gey2005model} provides oracle-type inequalities for the CART pruning algorithm proposed by \citep{breiman1984}, though the theory does not imply guarantees for out-out-sample prediction. What the existing literature currently lacks, however, is a more fine-grained analysis that reveals the unique advantages of tree learning with CART over other unstructured regression procedures, like vanilla $k$-NN or other kernel based estimators. Filling this theoretical gap, our main message is that, in certain settings, CART can identify low dimensional, latent structure in the data and adapt accordingly. We illustrate the adaptive properties of CART when the model is \emph{sparse}, namely, when the output depends only on a small, unknown subset of the input variables—thereby circumventing the curse of dimensionality.


Arguably the most difficult technical aspect of studying decision trees (and for that matter, any adaptive partitioning-based predictor) is understanding their approximation error, or at least pinning down conditions on the data that enable such an endeavor. Indeed, most existing convergence results for decision trees or ensembles thereof bound the expected (over the training data) prediction error by the size (i.e., the diameter) of the terminal nodes and show that they vanish with the depth of the tree, ensuring that the approximation error does so also \citep{denil2014, wager2018estimation}.

Others \citep{scornet2015} control the variation of the regression function inside the terminal nodes, without explicitly controlling their diameters, though the effect is similar. While these techniques can be useful to prove consistency statements, they are not generally delicate enough to capture the adaptive properties of the tree. It also often requires making strong assumptions about the tree construction. To address this shortcoming, in contrast, we use the fact that the prediction error is with high-probability (over the training data) bounded by the training error plus a complexity term. One of our crucial insights is that we can avoid using the node diameters as a proxy for the approximation error and, instead, directly bound the training error in terms of data dependent quantities (like the Pearson correlation coefficient) that are more transparent and interpretable, thereby facilitating our analysis and allowing us to prove more fine-grained results. 


\subsection{Learning setting}

Let us now describe the learning setting and framework that we will operate under for the rest of the paper. For clarity and ease of exposition, we focus specifically on \emph{regression trees}, where the target outcome is a continuous real value. 


We assume the training data is $ \calD_n = \{(\bX_1, Y_1), \dots, (\bX_n, Y_n) \} $, where $(\bX_i , Y_i ) $, $1\leq i \leq n $ are i.i.d. with common joint distribution $\mathbb{P}_{\bX, Y} $.

Here, $\bX_i \in [0, 1]^d $ is the input and $Y_i \in \mathbb{R} $ is a continuous response (or output) variable. A generic pair of variables will be denoted as $ (\bX, Y) $. A generic coordinate of $ \bX $ will be denoted by $ X $, unless there is a need to highlight the dependence on the $ j^{\text{th}}$ coordinate index, denoted by $ X_j $, or additionally on the $ i^{\text{th}} $ data point, denoted $ X_{ij} $.


Using squared error loss $ L(Y, Y') = (Y-Y')^2 $ as the performance metric, our goal is to predict $ Y $ at a new point $ \bX = \bx $ via a tree structured prediction rule $ \widehat Y(\bx) = \widehat Y(\bx; \calD_n) $. The training error and mean squared prediction error are, respectively, 
\begin{equation*}
\overline{\text{err}}(\widehat Y) \coloneqq \frac{1}{n}\sum_{i=1}^n(Y_i - \widehat Y(\bX_i))^2 \;\quad \text{and}\; \quad \text{Err}(\widehat Y) \coloneqq \Expect_{(\bX', Y')}[(Y'- \widehat Y(\bX'))^2],
\end{equation*}

where $ (\bX', Y') $ denotes an independent copy of $ (\bX, Y) $. For data $ \{(\bX_1, U_1, V_1), \dots, (\bX_n, U_n, V_n)\} $, we let $$ \widehat\rho\,(U, V \mid \bX \in A) \coloneqq \frac{\frac{1}{N}\sum_{\bX_i\in A}(U_i-\overline U)(V_i-\overline V)}{\sqrt{\frac{1}{N}\sum_{\bX_i\in A}(U_i-\overline U)^2\times \frac{1}{N}\sum_{\bX_i\in A}(V_i-\overline V)^2}}, $$  ($ A $ is a subset, $ N = \#\{\bX_i \in A \} $, $ \overline U = \frac{1}{N}\sum_{\bX_i \in A} U_i $, and $ \overline V = \frac{1}{N}\sum_{\bX_i \in A} V_i $) denote the empirical Pearson product-moment correlation coefficient, given $ \bX \in A $, and let $ \rho(U, V \mid \bX \in A) $ be its infinite sample counterpart. If $ U_i = g(X_{ij}) $ for a univariate function $ g(\cdot) $ of a coordinate $ X_j $, we write $ \widehat\rho\,(g(X_j), V \mid \bX \in A) $ or $ \rho\,(g(X_j), V \mid \bX \in A) $. For brevity, we let $ \widehat\sigma^2_Y $ denote the sample variance of the response values $ Y_1, Y_2, \dots, Y_n $ in the training data. The $ r^{\text{th}} $ derivative of a real valued function $ g(\cdot) $ is denoted by $ g^{(r)}(\cdot) $. Proofs of all forthcoming results are given in the supplement.


 
 
 
 
 


\section{Preliminaries} \label{sec:prelim}

As mentioned earlier, regression trees are commonly constructed with Classification and Regression Tree (CART) \citep{breiman1984} methodology. The primary objective of CART is to find partitions of the input variables that produce minimal variance of the response values (i.e., minimal sum of squares error with respect to the average response values). Because of the computational infeasibility of choosing the best overall partition, CART decision trees are greedily grown with a procedure in which binary splits recursively partition the tree into near-homogeneous terminal nodes. An effective binary split partitions the data from the parent tree node into two daughter nodes so that the resultant homogeneity of the daughter nodes, as measured through their \emph{impurity}, is improved from the homogeneity of the parent node.


The CART algorithm is comprised of two elements---a growing procedure and a pruning procedure. The growing procedure constructs from the data a maximal binary tree $T_{\text{max}}$ by the recursive partitioning scheme; the pruning procedure selects, among all the subtrees of $T_{\text{max}}$, a sequence of subtrees that greedily optimize a cost function. 

\subsection{Growing the tree}
Let us now describe the tree construction algorithm with additional detail. Consider splitting a regression tree $T$ at a node $\bt$. Let $s$ be a candidate split point for a generic variable $X$ that divides $\bt$ into left and right daughter nodes $\bt_L$ and $\bt_R$ according to whether $ X \leq s $ or $ X > s $. These two nodes will be denoted by $\bt_L = \{\bX \in \bt: X \leq s\} $ and $\bt_R = \{\bX \in \bt: X > s\} $. As mentioned previously, a tree is grown by recursively reducing node impurity. Impurity for regression trees is determined by the within node sample variance 
\begin{equation} \label{eq:impurity} \widehat\Delta(\bt) \coloneqq \widehat{\text{VAR}}(Y \mid \bX \in \bt) = \frac{1}{N(\bt) }\sum_{\bX_i \in \bt}(Y_i - \overline Y_{\bt})^2, \end{equation}
where $\overline Y_{\bt} = \frac{1}{N(\bt)}\sum_{\bX_i \in \bt}Y_i $ is the sample mean for $\bt$ and $N(\bt) =\# \{\bX_i \in \bt \}  $ is the number of data points in $\bt$. Similarly, the within node sample variance for a daughter node is $$ \widehat\Delta(\bt_L) = \frac{1}{N(\bt_L)}\sum_{\bX_i \in \bt_L}(Y_i - \overline Y_{\bt_L})^2, \quad \widehat\Delta(\bt_R) = \frac{1}{N(\bt_R)}\sum_{\bX_i \in \bt_R}(Y_i - \overline Y_{\bt_R})^2, $$ where $\overline Y_{\bt_L}$ is the sample mean for $\bt_L$ and $N(\bt_L)$ is the sample size of $\bt_L$ (similar definitions apply to $\bt_R$). The parent node $ \bt $ is split into two daughter nodes using the variable and split point producing the largest decrease in impurity. For a candidate split $s$ for $X$, this decrease in impurity equals \cite[Definition 8.13]{breiman1984} \begin{equation} \label{eq:data} \widehat\Delta(s, \bt) \coloneqq \widehat\Delta(\bt) - [\widehat P(\bt_L)\widehat\Delta(\bt_L) + \widehat P(\bt_R)\widehat\Delta(\bt_R)], \end{equation} where $ \widehat P(\bt_L) = N(\bt_L) /N(\bt) $ and $\widehat P(\bt_R) = N(\bt_R) /N(\bt) $ are the proportions of data points in $ \bt $ that are contained in $\bt_L$ and $\bt_R$, respectively.

The tree $T$ is grown recursively by finding the variable $ \hat\jmath $ and split point $\hat s = \hat s_{\hat\jmath}$ that maximizes $\widehat\Delta(s, \bt)$. Note that for notational brevity, we suppress the dependence on the input coordinate index $ j $. The output $ \widehat Y(T) $ of the tree at a terminal node $ \bt $ is the least squares predictor, namely, $ \widehat Y(T,\bx) = \overline Y_{\bt} $ for all $ \bx \in \bt $.
\subsection{Pruning the tree}


The CART growing procedure stops once a maximal binary tree $T_{\text{max}}$ is grown (i.e., when the terminal nodes contain at least a single data point). However, $\widehat Y(T_{\text{max}})$ is generally not a good predictor, since it will tend to overfit the data and therefore generalize poorly to unseen data. This effect can be mitigated by complexity regularization. Removing portions of the overly complex tree (i.e., via pruning) is one way of reducing its complexity and improving performance. We will now describe such a procedure.

We say that $ T $ is a pruned subtree of $ T' $, written as $ T\preceq T' $, if $ T $ can be obtained from $ T' $ by iteratively merging

any number of its internal nodes. A pruned subtree of $T_{\text{max}}$ is defined as any binary subtree of $T_{\text{max}}$ having the same root node as $T_{\text{max}}$. The number of terminal nodes in a tree $ T $ is denoted $ |T| $. Given a subtree $ T $ and temperature $ \alpha > 0 $, we define the penalized cost function 
\begin{equation} \label{eq:cost}
R_{\alpha}(\widehat Y(T)) \coloneqq \overline{\text{err}}(\widehat Y(T))+\alpha|T|.
\end{equation}
As shown in \citep[Section 10.2]{breiman1984}, the smallest minimizing subtree for the temperature $\alpha$,
$$
\widehat T \in \argmin_{T\preceq T_{\text{max}}}R_{\alpha}(\widehat Y(T)),
$$
exists and is unique (smallest in the sense that if $ T' \in \argmin_{T\preceq T_{\text{max}}}R_{\alpha}(\widehat Y(T)) $, then $ \widehat T \preceq T' $). For a fixed $ \alpha $, the optimal subtree $ \widehat T $ can be found efficiently by weakest link pruning \citep{breiman1984, hastie2009elements}, i.e., by successively collapsing the internal node that increases $ \overline{\text{err}}(\widehat Y(T)) $ the least, until we arrive at the single-node tree consisting of the root node. Good values of $ \alpha $ can be selected using cross-validation, for example, though analyzing the effect of such a procedure is outside the scope of the present paper.

Our first result shows that, with high probability, the test error of the pruned tree $ \widehat T $ on new data is bounded by a multiple of $ \min_{T\preceq T_{\text{max}}}R_{\alpha}(\widehat Y(T)) $ plus lower order terms.

\begin{theorem} \label{thm:test}
Let $ \widehat T $ be the smallest minimizer of \prettyref{eq:cost}. Suppose $ Y = f(\bX) $, $ B = \sup_{\bx}|f(\bx)| < \infty $, $ n > (d+1)/2 $, 

and $ \alpha > \frac{27B^2(d+1)\log(2en/(d+1))}{n} $.

Then, with probability at least $ 1-\delta $ over the training sample $ \calD_n $,
$$

\text{Err}(\widehat Y(\widehat T)) \leq 4\min_{T\preceq T_{\text{max}}}R_{\alpha}(\widehat Y(T)) + \frac{54B^2\log(2/\delta)}{n}.
$$
\end{theorem}


Similar bounds hold for the Bayes risk for binary classification, i.e., $ Y \in \{ 0, 1 \} $, since, in this case, the squared error impurity \prettyref{eq:impurity} equals one-half of the so-called \emph{Gini} impurity used for classification trees. See also \citep{nobel2002analysis} for results of a similar flavor when the penalty is proportional to $ \sqrt{|T|} $. 


In what follows, we let $T_K \preceq T_{\text{max}}$ denote a fully grown binary tree of depth $ K = \Theta(\log_2(n))  $, i.e., we stop splitting if (1) the node contains a single data point, (2) all response values in the node are the same, or (3) a depth of $ K $ is reached, whichever occurs sooner. We also let $\widehat T $ be the smallest minimizer of the cost function \prettyref{eq:cost} with temperature $ \alpha = \Theta((d/n)\log(n/d)) $.


\section{Bounding the Training Error} \label{sec:training}
In the previous section, \prettyref{thm:test} showed that, with high probability, the test error is bounded by a multiple of the cost function \prettyref{eq:cost} at its minimum (plus lower order terms). Since the cost function is defined as the training error plus penalty term, the next step in our course of study is to understand how the training error of CART behaves.

\subsection{Splitting criterion and Pearson correlation} \label{sec:pearson}

Before we begin our analysis of the training error, we first digress back to the tree construction algorithm and give an alternative characterization of the objective. Now, the use of the sum of squares impurity criterion $ \widehat\Delta(s, \bt) $ with averages in the terminal nodes permits further simplifications of the formula \prettyref{eq:data} above. For example, using the sum of squares decomposition, $ \widehat\Delta(s, \bt) $ can equivalently be expressed as \citep[Section 9.3]{breiman1984}
\begin{equation} \label{eq:simple}
\widehat P(\bt_L)\widehat P(\bt_R)\big(\overline Y_{\bt_L}-\overline Y_{\bt_R}\big)^2,
\end{equation}
which is commonly used for its computational appeal---that is, one can find the best split for a continuous variable with just a single pass over the data, without the need to calculate multiple averages and sums of squared differences for these averages, as required with \prettyref{eq:data}.
Yet another way to view $ \widehat\Delta(s, \bt) $, which does not appear to have been considered in past literature and will prove to be useful for our purposes, is via its equivalent representation as $ \widehat\Delta(\bt) \times \widehat\rho^{\,2}(\widetilde Y, Y \mid \bX \in \bt) $, where 
\begin{equation} \label{eq:correlation} \widehat\rho\,(\widetilde Y, Y \mid \bX \in \bt) \coloneqq \frac{\frac{1}{N(\bt)}\sum_{\bX_i\in\bt}(\widetilde Y_i- \overline Y_{\bt})(Y_i- \overline Y_{\bt})}{\sqrt{\frac{1}{N(\bt) }\sum_{\bX_i \in \bt}(\widetilde Y_i - \overline Y_{\bt})^2\times\frac{1}{N(\bt) }\sum_{\bX_i \in \bt}(Y_i - \overline Y_{\bt})^2}} \geq 0
\end{equation}
is the Pearson product-moment correlation coefficient between the decision stump 
\begin{equation} \label{eq:stump}
 \widetilde Y \coloneqq \overline Y_{\bt_L}\mathbf{1}(X \leq s) + \overline Y_{\bt_R}\mathbf{1}(X > s) 
\end{equation}
and response variable $ Y $ within $ \bt $ (for the proof, see \prettyref{lmm:identity} in the supplement).


Hence, at each node, CART seeks the decision stump most correlated in magnitude with the response variable along a particular variable, i.e.,
\begin{equation} \label{eq:split_cor}
\hat s \in \argmax_{s}\widehat\Delta(s, \bt) = \argmax_{s}\widehat\rho\,(\widetilde Y, Y \mid \bX \in \bt).
\end{equation}
Equivalently, CART splits along variables with decision stumps that are most correlated with the residuals $ Y_i  - \overline Y_{\bt} $ of the current fit $ \overline Y_{\bt} $. As with $ r^2 $ for simple linear regression, we see that the squared correlation $ \widehat\rho^{\,2}(\widetilde Y, Y \mid \bX \in \bt) $ equals the \emph{coefficient of determination} $ R^2 $, in the sense that it describes the fraction of the variance in $ Y $ that is explained by a decision stump $ \widetilde Y $ in $ X $, since $ \widehat\rho^{\,2}(\widetilde Y, Y \mid \bX \in \bt) = 1 - \frac{1}{N(\bt)}\sum_{\bX\in\bt}(Y_i-\widetilde Y_i)^2/\frac{1}{N(\bt)}\sum_{\bX\in\bt}(Y_i-\overline Y_{\bt})^2 $. 


\begin{definition}
We let $ \widehat Y  $ denote a decision stump $ \widetilde Y $ with an optimal direction $ \hat\jmath \in \argmax_{j=1,2,\dots, d}\widehat\Delta(\hat s, \bt) $ and corresponding optimal split $ \hat s $.
\end{definition}
 It should be stressed that the alternative characterization of the splitting criterion \prettyref{eq:data} in terms of a correlation is unique to the squared error impurity with (constant) averages in the terminal nodes.
 


We now introduce a data dependent quantity that will play a central role in determining the rates of convergence of the prediction error. For a univariate function class $ \calG $, we let $ \widehat\rho_{\calG} $ be the largest Pearson correlation between the response data $ Y $ and a function in $ \calG $ of a single input coordinate for a worst-case node, i.e.,
\begin{equation} \label{eq:cor_minimax}
\widehat\rho_{\calG} \coloneqq \adjustlimits\min_{\bt} \sup_{g(\cdot)\in\calG, \; j = 1, 2, \dots, d}|\widehat\rho\,(g(X_j), Y \mid \bX \in \bt)|,
\end{equation}
where the minimum runs over all internal nodes $ \bt $ in $ T_K $. We will specifically focus on classes $ \calG $ that consist of decision stumps, and more generally, monotone functions.

\subsection{Location of splits and Pearson correlation} \label{sec:bias}

Having already revealed the intimate role the correlation between the decision stump and response values \prettyref{eq:correlation} plays in the tree construction, it is instructive to explore this relationship with the location of the splits. In order to study this cleanly, let us for the moment work in an asymptotic data setting to determine the coordinates to split and their split points, i.e.,
\begin{equation} \label{eq:infinite} 
\widehat\Delta(s, \bt) \underset{n\rightarrow\infty}{\rightarrow} \Delta(s, \bt) \coloneqq \Delta(\bt) - [P(\bt_L)\Delta(\bt_L) + P(\bt_R)\Delta(\bt_R)], \end{equation}
where quantities without hats are the population level counterparts of the empirical quantities defined previously in \prettyref{eq:data}. A decision stump \prettyref{eq:stump} with an optimal theoretical direction $ j^* $ and corresponding optimal theoretical split $ s^* = s^*_{j^*} $ is denoted by $ \widehat Y^* $. Now, if the number of data points within $ \bt $ is large and $ \Delta(s, \bt) $ has a unique global maximum, then we can expect $ \hat s  \approx s^*$ (via an empirical process argument) and hence the infinite sample setting is a good approximation to CART with empirical splits, giving us some insights into its dynamics. Indeed, if $ s^* $ is unique, \citep[Theorem 2]{ishwaran2015effect} shows that $ \hat s $ converges in probability to $ s^* $. With additional assumptions, one can go even further and characterize the rate of convergence. For example, \citep[Section 3.4.2]{buhlmann2002analyzing} and \citep{banerjee2007confidence} provide cube root asymptotics for $ \hat s $, i.e., $ n^{1/3}(\hat s - s^*) $ converges in distribution. 


Each node $ \bt $ is a Cartesian product of intervals. As such, the interval along variable $ X $ in $ \bt $ is denoted by $ [a, b] $, where $ a < b $. The next theorem characterizes the relationship between an optimal theoretical split $ s^* $ and infinite sample correlation $ \rho(\widehat Y , Y \mid \bX \in \bt) \overset{\text{a.s.}}{\coloneqq} \lim_n\widehat\rho\,(\widehat Y^* , Y \mid \bX \in \bt) $ for a deterministic node $ \bt $ (the limit exists by the uniform law of large numbers). The proof is based on the first-order optimality condition, namely, $ \frac{\partial}{\partial s}\Delta(s, \bt) \mid_{s=s^*} = 0 $.


\begin{theorem} \label{thm:mainprob} 


Suppose $ \bX $ is uniformly distributed on $ [0, 1]^d $ and $ \Delta(s^*, \bt) > 0 $. For a deterministic parent node $ \bt $, an optimal theoretical split $ s^* \in [a, b] $ along variable $ X $ has the form
\begin{equation} \label{eq:uniformsplit}
\frac{a+b}{2}\;\pm\;\frac{b-a}{2}\sqrt{\frac{v}{ v+\rho^2(\widehat Y^* , Y \mid \bX \in \bt)}},
\end{equation}
where
$
v = \frac{(\expect{Y\mid \bX \in \bt, \; X = s^*}-\expect{Y\mid \bX \in \bt})^2}{\text{VAR}(Y\mid\bX\in\bt)}.
$


\end{theorem}

Expression \prettyref{eq:uniformsplit} in \prettyref{thm:mainprob} reveals that an optimal theoretical split $ s^* $ is a perturbation of the median $ (a+b)/2 $ of the conditional distribution $ X \mid \bX \in \bt $, where the gap is governed by the correlation $ \rho(\widehat Y^* , Y \mid \bX \in \bt) $. These correlations control the local and directional granularity of the partition of the input domain. Splits along input coordinates that contain a strong signal, i.e., $ \rho(\widehat Y^* , Y \mid \bX \in \bt) \gg 0 $, tend to be further away from the parent node edges, thereby producing side lengths $ [a,b] $ that are on average narrower. At the other extreme, the correlation is weakest when there is no signal in the splitting direction or when the response values in the node are not fit well by a decision stump---yielding either 
$ s^*\approx a + (b-a)\rho^2(\widehat Y^* , Y \mid \bX \in \bt)/(4v) $ or $ s^* \approx b - (b-a)\rho^2(\widehat Y^* , Y \mid \bX \in \bt)/(4v)$---and hence the predicted output in one of the daughter nodes does not change by much. 
For example, if $ Y = g(X) $ is a sinusoidal waveform with large frequency $ w $ (not fit well by a single decision stump) and $ \bt $ is the root node $ [0, 1]^d $, then $ v = \Theta(1) $ and $ \rho(\widehat Y^* , Y \mid \bX \in \bt) = \Theta(1/\sqrt{w}) $, and hence by \prettyref{eq:uniformsplit}, either $ s^* = \Theta(1/w) $ or $ s^* = 1-\Theta(1/w)  $ (see \prettyref{lmm:example} in the supplement).  
This phenomenon, where optimal splits concentrate at the endpoints of the node along noisy directions, has been dubbed `end-cut preference' in the literature and has been known empirically since the inception of CART \citep{ishwaran2015effect}, \citep[Section 11.8]{breiman1984}. The theory above is also consistent with empirical studies on the adaptive properties of Breiman's random forests which use CART \citep[Section 4]{lin2006}.


\subsection{Training error and Pearson correlation} 
In addition to determining the location of the splits, the correlation is also directly connected to the training error. Intuitively, the training error should small when CART finds decision stumps that have strong correlation with the response values in each node.

More precisely, the following lemma reveals the importance of the correlation \prettyref{eq:correlation} in controlling the training error. It shows that each time a node $ \bt $ is split, yielding an optimal decision stump $ \widehat Y $, the training error in $ \bt $ is reduced by a constant factor, namely, $ \exp(-\widehat\rho^{\,2}(\widehat Y , Y \mid \bX \in \bt) ) $ or, uniformly, by $ \exp(-\widehat\rho^{\,2}_{\calH}) $, where $ \calH $ is the collection of all decision stumps (i.e., a step function with two constant pieces) and $ \widehat\rho_{\calH} $ is the quantity defined in \prettyref{eq:cor_minimax}. Recursing this contraction inequality over nodes at each level of the tree leads to the conclusion that the training error should be exponentially small in the depth $ K $, provided the correlation at each node is large.


\begin{lemma}\label{lmm:training}
Almost surely,
\begin{equation} \label{eq:training}

\frac{1}{N(\bt)}\sum_{\bX_i\in\bt}(Y_i- \widehat Y_i)^2 \leq \frac{1}{N(\bt)}\sum_{\bX_i\in\bt}(Y_i-\overline Y_{\bt})^2 \times \exp(-\widehat\rho^{\,2}(\widehat Y, Y \mid \bX \in \bt) ),
\end{equation}
and hence
\begin{equation} \label{eq:train_exp}
\overline{\text{err}}(\widehat Y(T_K)) \leq \widehat\sigma^2_Y\exp(-K\times \widehat \rho_{\calH}^{\,2}),
\end{equation}
where $ \calH $ is the collection of all decision stumps.

\end{lemma}


\subsection{Size of Pearson correlation via comparison inequalities} \label{sec:size}

Due to the importance of the correlation in controlling the training error, it is natural to ask when it will be large. We accomplish this by studying its size relative to the correlation between the data and another more flexible model. That is, we fit an arbitrary univariate function $ g(X) $ of a generic coordinate $ X $ to the data in the node and ask how large $ \widehat\rho\,(\widehat Y , Y \mid \bX \in \bt) $ is relative to $ |\widehat\rho\,(g(X), Y \mid \bX \in \bt)| $. Such a relationship will enable us to conclude that if $ Y $ is locally correlated with $ g(X) $ in the node, then so will $ Y $ with an optimal decision stump $ \widehat Y  $.  

Before we continue, let us mention that studying $ \widehat\rho\,(\widehat Y , Y \mid \bX \in \bt) $ directly is hopeless since it not at all straightforward to disentangle the dependence on the data. Even if this could be done and a target population level quantity could be identified, it is difficult to rely on concentration of measure when $ \bt $ contains very few data points; a likely situation among deep nodes. Nevertheless, by definition of $ \widehat Y  $ via \prettyref{eq:split_cor}, we can construct a prior $ \Pi(j, s) $ on coordinates $ j $ and splits $ s $, and lower bound $ \widehat\rho\,(\widehat Y , Y \mid \bX \in \bt) $ by
\begin{equation} \label{eq:integral}
\int \widehat\rho\,(\widetilde Y, Y \mid \bX \in \bt)d\hspace{0.02cm}\Pi(j, s),
\end{equation}
which is much less burdensome to analyze. Importantly, the prior can involve unknown quantities from the distribution of $ (\bX, Y) $. For a special choice of prior $ \Pi $, \prettyref{eq:integral} can be further lower bounded by
\begin{equation} \label{eq:cor_constant}
 \widehat\rho\,(\widehat Y , Y \mid \bX \in \bt)  \geq \text{constant} \times |\widehat\rho\,(g(X), Y \mid \bX \in \bt)|.
\end{equation}
The constant in \prettyref{eq:cor_constant} depends on $ g(\cdot) $, though importantly it is invariant to the scale of $ g(\cdot) $. If $ g(\cdot) $ belongs to a univariate model class $ \calG $, this constant can either be studied directly for the specific $ g(\cdot) $ or minimized over $ g(\cdot) \in \calG $ to yield a more insightful lower bound. For certain model classes $ \calG $, the minimization problem turns out to be equivalent to a quadratic program, and


the solution can be obtained explicitly and used to prove the next set of results.


Our first application of this technique shows that, despite fitting the data to a decision stump with \emph{one} degree of freedom (i.e., the location of the split), 

CART behaves almost as if it fit the data to a monotone function with $ N(\bt)-1$ degrees of freedom, at the expense of a sublogarithmic factor in $ N(\bt) $. For example, the correlation between the response variable and the decision stump is, up to a sub-logarithmic factor, at least as strong as the correlation between the response variable and a linear or

isotonic fit. 


\begin{fact} \label{fact:cor}
Almost surely, uniformly over all monotone functions $ g(\cdot) $ of $ X $ in the node, we have
\begin{equation} \label{eq:cor}
\widehat\rho\,(\widehat Y , Y \mid \bX \in \bt) \geq \frac{1}{\sqrt{1+\log(2N(\bt))}}\times |\widehat\rho\,(g(X), Y \mid \bX \in \bt)|.
\end{equation}
\end{fact}

 \begin{wrapfigure}[37]{r}{0.35\textwidth}
\begin{subfigure}[t]{0.35\textwidth}
\centering
  \includegraphics[width=1\linewidth]{sparse.pdf}
  \caption{Synthetic data. Prediction error of CART vs. $k$-NN as $ d $ varies.}
   \label{fig:example}
 \end{subfigure}
 \begin{subfigure}[t]{0.35\textwidth}
\centering
  \includegraphics[width=1\linewidth]{boston.pdf}
  \caption{Boston housing data. Prediction error of CART vs. $k$-NN as $ d $ varies.}
   \label{fig:boston}
 \end{subfigure}
 \begin{subfigure}[t]{0.35\textwidth}
 \centering
  \includegraphics[width=1\linewidth]{correlation.pdf}
  \caption{Minimum correlation $ \widehat\rho_{\calH} $ (averaged over $ 10 $ independent replications) of all nodes of pruned CART tree as sparsity $ d_0 $ varies.}
   \label{fig:correlation}
 \end{subfigure}
\end{wrapfigure}
The previous fact also suggests that CART is quite good at fitting response values that have a local, low-dimensional, monotone relationship with the input variables. 

Note that because correlation is merely a measure of linear association, $ |\widehat\rho\,(g(X), Y \mid \bX \in \bt)| $ can still be large for some monotone $ g(\cdot) $, even if $ Y $ is not approximately monotone in one coordinate. That is, $ Y $ need only be locally correlated with such a function.  On the other hand, if $ Y $ has no signal in $ X $, then \cite[Lemma 1, Supplement]{li2019debiased} show that, with high probability, $ \widehat\rho\,(\widehat Y, Y \mid \bX \in \bt) $ is $ \calO(\sqrt{(\log N(\bt))/N(\bt)}) $.


\section{Main Results} \label{sec:results}
In this section, we use the training error bound \prettyref{eq:train_exp} and the device \prettyref{eq:integral} for obtaining correlation comparison inequalities (\`a la \prettyref{fact:cor}) to give bounds on the prediction error of CART. We first outline the high-level strategy. By \prettyref{thm:test}, with high probability, the leading behavior of the prediction error $ \text{Err}(\widehat Y(\widehat T)) $ is governed by
$ \inf_{T \preceq T_{\text{max}}}R_{\alpha}(\widehat Y(T)), $
which is smaller than the minimum of $ R_{\alpha}(\widehat Y(T_K)) = \overline{\text{err}}(\widehat Y(T_K)) + \alpha|T_K| $ over all fully grown trees $ T_K $ of depth $ K $ with $ |T_K| \leq 2^K $, i.e.,
\begin{equation} \label{eq:resolve_K_body}
\inf_{K \geq 1}\{\overline{\text{err}}(\widehat Y(T_K)) + \alpha 2^K\}.
\end{equation}
Coupled with an informative bound on $ \overline{\text{err}}(\widehat Y(T_K)) $, \prettyref{eq:resolve_K_body} can then be further bounded and solved. The proofs reveal that a good balance between the tree size and its goodness of fit occurs when $ K $ is logarithmic in the sample size. 

\subsection{Asymptotic consistency rates for sparse additive models}

Applying the training error bound \prettyref{eq:train_exp} to \prettyref{eq:resolve_K_body} with $ K = \lceil(\widehat\rho_{\calH}^{\,2}+\log 2)^{-1}\log(\widehat\sigma^2_Y/\alpha) \rceil $, we have from \prettyref{thm:test} that with probability at least $ 1-\delta $,
\begin{equation} \label{eq:test_poly}

\text{Err}(\widehat Y(\widehat T)) = \calO\Big(\widehat\sigma^2_Y\Big(\frac{d\log(n/d)}{n\widehat\sigma^2_Y}\Big)^{\frac{\widehat\rho_{\calH}^{\,2}}{\widehat\rho_{\calH}^{\,2}+\log 2}} + \frac{\log(1/\delta)}{n} \Big).
\end{equation}


It turns out that if $ \bX $ is uniformly distributed and $ Y $ is a sparse additive model with $ d_0 $ component functions $ g_j(\cdot) $, then $ \widehat\rho_{\calH}^{\,2} $ is asymptotically lower bounded by a constant multiple of $ 1/d_0 $. 


Thus, we find from \prettyref{eq:test_poly} that if $ d_0 $ is fixed, then
$
\lim_n\text{Err}(\widehat Y(\widehat T)) \overset{\text{a.s.}}{=} 0
$
even when the ambient dimension grows as $ d = o(n) $. Note that such a statement is not possible for vanilla $ k $-NN or other kernel based regression methods with nonadaptive weights, unless feature selection is performed beforehand. In fact, we show next that the prediction error rate that CART achieves is the same as what would be achieved by a standard kernel predictor \emph{if} one had a priori knowledge of the locations of the $ d_0 $ relevant input variables that determine the output. A routine computer experiment on synthetic or real data easily confirms this theory. 
In \prettyref{fig:example} and \prettyref{fig:correlation}, we generate $ 1000 $ samples from the model $ Y = \sum_{j=1}^{d_0} g_j(X_j) $, where each $ g_j(X_j) $ equals $ \pm X^2_j $ (alternating signs) and $ \bX \sim \mbox{Uniform}([0, 1]^d) $. In \prettyref{fig:example}, we plot the test error, averaged over $ 10 $ independent replications, of pruned CART vs. $ k$-NN (with cross-validated $ k $) as $ d $ ranges from $ 5 $ to $ 100 $ with $ d_0 = 5 $ fixed. A similar experiment is performed in \prettyref{fig:boston} on the Boston housing dataset \citep[Section 8.2]{breiman1984} ($d_0 = 10$ and $ n = 506 $), where we scale the inputs to be in $ [0, 1] $ and add $ d-d_0 $ noisy $ \text{Uniform}([0, 1]) $ input variables.
 According to \prettyref{thm:dim_limit}, the convergence rate of CART depends primarily on the sparsity $ d_0 $ and therefore its performance should not be adversely affected by growing $ d $. Consistent with our theory, the prediction error of CART remains stable as $ d $ increases, whereas $k$-NN does not adapt to the sparsity. Furthermore, \prettyref{fig:correlation} illustrates how $ \widehat\rho^{\,2}_{\calH} $ decays with $ d_0 $ if $ d = 20 $, thus corroborating with the aforementioned asymptotic behavior of $ \Omega(1/d_0) $.


\begin{theorem} \label{thm:dim_limit}
Suppose $ \bX $ is uniformly distributed on $ [0, 1]^d$ and
$Y = \sum_j g_j(X_j) $ is a sparse additive model with $ d_0 \ll d $ smooth component functions $ g_j(\cdot) $, where each function is not too locally `flat' in the sense that


\begin{equation} \label{eq:hypothesis}
\sup_x \,\inf\{r \geq 1:g^{(r)}_j(\cdot)\; \text{exists, continuous, and nonzero at}\; x \} < \infty.
\end{equation}

Then there exists a constant $ C > 0 $ that is independent of $ d_0 $ such that, almost surely,

$
\liminf_n \widehat\rho_{\calH}^{\,2}  \geq C/d_0,
$

and
\begin{equation} \label{eq:limsup} \limsup_{n} \frac{\text{Err}(\widehat Y(\widehat T))}{((d/n)\log(n/d))^{\Omega(1/d_0)}} \overset{\text{a.s.}}{=} \calO(1). \end{equation}


\end{theorem}


\begin{remark}
For independent, continuous marginal input variables $ X_j $, there is no loss of generality in assuming uniform distributions in \prettyref{thm:dim_limit}. Indeed, CART decision trees are invariant to strictly monotone transformations of $ X_j $. One such transformation is the marginal cumulative distribution function $ F_{X_j}(\cdot) $, for which $ F_{X_j}(X_j) \sim \mbox{Uniform}([0, 1]) $—and so the problem can immediately be reduced to the uniform case.
\end{remark}
Any nonconstant component function $ g_j(\cdot) $ that admits a power series representation satisfies the hypothesis of \prettyref{thm:dim_limit}, though, in general, the condition \prettyref{eq:hypothesis} accommodates functions that are not analytic or infinitely differentiable.
In fact, even differentiability is not necessary—similar results hold if the $ g_j(\cdot) $ are \emph{step functions}, as we now show.


To this end, assume that $ Y = \sum_j g_j(X_j) $ is an additive model, where each component function $ g_j(\cdot) $ is a bounded step function and the total number of constant pieces of $ Y $ is $ V $.


We show in the supplement that each optimal split $ \hat s $ in a node $ \bt $ satisfies
\begin{equation} \label{eq:splitinterval}
\max\{X_{ij} \in I: \bX_i \in \bt\}\leq \hat s \leq \min\{X_{ij} \in I': \bX_i \in \bt\},
\end{equation}
for some direction $ j $ and successive intervals $ I $ and $ I' $ on which $ g_j(\cdot) $ is constant. For example, if $ Y = c_1 \mathbf{1}(X_1< s_1)+ c_2\mathbf{1}(X_2< s_2) $, then the first split separates the data along $ X_1 $ at $ s_1 $ (resp. $ X_2 $ at $ s_2 $), and at the next level down, CART separates the data in both daughter nodes along $ X_2 $ at $ s_2 $ (resp. $ X_1 $ at $ s_1 $). 
Thus, in general, each empirical split always perfectly separates the data in the node between two adjacent constant values of a component function.


A CART decision tree grown in this way will achieve zero training error once it has at least $ V $ terminal nodes, i.e., $ |T| \geq V $. This is in fact the same training error that would result from the global least squares projection onto the space of all step functions with $ V $ constant pieces. From \prettyref{thm:test}, we immediately obtain the following performance bound, which is the optimal $\calO(1/n) $ parametric rate for prediction, up to a logarithmic factor in the sample size \citep{yao1989least}. Notice that we do not make any assumptions on input distribution.
\begin{theorem} \label{thm:pconstant}
Suppose $ Y = \sum_j g_j(X_j) $, where each component function $ g_j(\cdot) $ is a bounded step function and the total number of constant pieces of $ Y $ is $ V $. With probability at least $ 1-\delta $,
\begin{equation} \label{eq:steprisk}
\text{Err}(\widehat Y(\widehat T)) = \calO\Big(\frac{Vd\log(n/d)}{n} +\frac{\log(1/\delta)}{n} \Big).
\end{equation}
\end{theorem}


\subsection{Finite sample consistency rates and general sparse models} \label{sec:finite}


Using \prettyref{fact:cor}, we now provide results of a similar flavor for more general regression models under a mild assumption on the largest number of data points in a node at level $ k $ in $ T_K $, denoted by $ N_k $. 


Importantly, our theory only requires that each $ N(\bt) $ is \emph{upper} bounded at every level of the tree. This condition still allows for nodes that have very few data points, which is typical for trees trained in practice. Contrast this assumption with past work on tree learning (including tree ensembles like random forests) that requires each $ N(\bt) $ to be \emph{lower} bounded \cite[Section 12.2]{breiman1984}, \citep{meinshausen2006quantile, wager2018estimation, denil2014}. 


\begin{assumption} \label{ass:samp_size1}
For some constants $ a \geq 0 $ and $ A > 0 $, the largest number of data points in a node at level $ k $ in $ T_K $ satisfies
$ N_k \leq Ank^a/2^k $, for $ k = 1, 2, \dots, K = \Theta(\log_2(n)) $.
\end{assumption}


Recall the quantity $ \widehat\rho_{\calG} $ defined in \prettyref{eq:cor_minimax}, namely, the largest correlation between the response data $ Y $ and a function in $ \calG $ for a worst-case node. Our next theorem shows that if $ \calM $ is the collection of all monotone (i.e., increasing or decreasing) functions, then 
$$ \widehat\rho_{\calM} = \adjustlimits\min_{\bt} \sup_{g(\cdot) \; \text{monotone}, \; j = 1, 2, \dots, d}|\widehat\rho\,(g(X_j), Y \mid \bX \in \bt)| $$
governs the rate at which the training error and prediction error decrease. Both errors are small if the local monotone dependence between $ \bX $ and $ Y $ is high; that is, if CART partitions the input domain into pieces where the response variable is locally monotone in a few of the input coordinates. 
\begin{theorem} \label{thm:monotone}
Let $ Y = f(\bX) $, where $ f(\cdot) $ is a bounded function. Under \prettyref{ass:samp_size1}, almost surely,
\begin{equation} \label{eq:train1}
\overline{\text{err}}(\widehat Y(T_K)) 

\leq \widehat\sigma^2_Y\Big(1 - \frac{K}{\log_2(4K^aAn)}\Big)^{\widehat\rho^{\,2}_{\calM}}.
\end{equation}


Furthermore, with probability at least $ 1-\delta $,
\begin{equation} \label{eq:test1}
\text{Err}(\widehat Y(\widehat T)) = \calO\Big( \widehat\sigma^2_Y\Big( \frac{\log((d/\widehat\sigma^2_Y)\log^{2+a}(n))}{\log(n)}\Big)^{\widehat\rho^{\,2}_{\calM}} + \frac{\log(1/\delta)}{n} \Big).
\end{equation}


\end{theorem}


We will now argue that $ \widehat \rho_{\calM} $ is an empirical measure of the local dimensionality of $ Y $. More specifically, we argue that if CART effectively partitions the input domain so that,

in each node, $ Y $ is locally correlated with sparse additive models with $ d_0 \ll d $ monotone component functions,

then $ \widehat\rho^{\,2}_{\calM}= \Omega(1/d_0) $.


To see why this assertion is true, suppose $ g_1(X_1), g_2(X_2), \dots, g_d(X_d) $ is an arbitrary collection of $ d $ univariate functions from $ \calM $.

However, suppose that only a subset of $d_0$ of the input variables $ X_1, X_2, \dots, X_d $ locally affect $ Y $ in each node. Then,

it can be shown (see \prettyref{lmm:correlation_dim} in the supplement) that there is some node $ \bt $ and sparse additive model $Y_0 $ with $d_0$ component functions of the form $ \pm g_j(X_j) $, corresponding to the $ d_0 $ input variables that locally affect $Y$, such that, almost surely,
\begin{equation} \label{eq:cor_dim}

 \widehat\rho^{\,2}_{\calM} \geq \frac{\widehat\rho^{\,2}(Y_0, Y \mid \bX \in \bt)}{d_0} = \Omega(1/d_0).
\end{equation}


The above statement is reminiscent of \prettyref{thm:dim_limit} in which $ \widehat\rho^{\,2}_{\calH} = \Omega(1/d_0) $ controls the convergence rate of the prediction error when the underlying regression model is additive. Though, in contrast, note that \prettyref{eq:cor_dim} holds regardless of the dependence structure between the $ d_0 $ input coordinates that matter and the $ d-d_0 $ input coordinates that do not. 


Thus, \prettyref{eq:cor_dim} and \prettyref{thm:monotone} together suggest that it is possible to achieve rates of the form
$ (\log(d)/\log(n))^{\Omega(1/d_0)}$ in fairly general settings.


\section{Extensions to Tree Ensembles}

One of the key insights for our analysis of CART was the ability to connect the training error to the objective function of the growing procedure, as in \prettyref{lmm:training}. 

Establishing similar relationships is not as easy with trees that are constructed from bootstrap samples or random selections of input variables. Nevertheless, we mention a few ideas for future work. Suppose $ \widehat Y = (1/m) \sum_T \widehat Y(T) $ is the output of an ensemble of $ m $ decision trees. By convexity of the squared error loss \citep[Section 11]{breiman2001} or \citep[Section 4.1]{breiman1996bagging}, we have $ \text{Err}(\widehat Y) \leq  (1/m)\sum_T \text{Err}(\widehat Y(T)) $, where the prediction error is averaged with respect to the tree randomization mechanism. Using the law of total expectation by conditioning on each realization of the (random) tree, $ \text{Err}(\widehat Y(T)) $ can be decomposed into quantities that involve the prediction error of a fixed (non-random) tree, for which our previous results can be applied. 


We will leave the exact details of these extensions for future work.

\section{Conclusion} \label{sec:discussion}


A key strength of CART decision trees is that they can exploit local, low dimensionality of the model—via a built-in, automatic dimension reduction mechanism. This is particularly useful since many real-world input/output systems are locally approximated by simple model forms with only a few variables. 

Adaptivity with CART is made possible by the recursive partitioning of the input space, in which optimal splits are increasingly affected by local qualities of the data as the tree is grown. To illustrate this ability, we identified settings where CART adapts to the unknown sparsity of the model. To the best of our knowledge, the consistency rates given here are the first of their kind for CART decision trees.


\newpage


\section*{Broader Impact}

\textbf{Who may benefit from this research.}  There are at least two groups of people who will benefit—either directly or indirectly—from this research.
\begin{enumerate}
\item   Decision makers across a variety of domains, especially those with limited training in statistics. CART has enabled data-driven decision making in multiple high-stakes domains (e.g., business, medicine, and policy) over the past three decades. In particular, those who do not have a formal quantitative background will benefit from the intuitive and interpretable nature of CART and its quick and easy implementation.
\item Members of the society who may be have faced ethical/fairness concerns associated with their data and its use. As this paper has demonstrated, CART forms predictions by emphasizing variables that are more relevant to the output. In a social science context, this suggests that CART may focus more on key diagnostic information (e.g., education, income) without being influenced by potentially non-diagnostic variables that other methods may have falsely deemed relevant (e.g., gender, race).
\end{enumerate}


\textbf{Who may be put at a disadvantage from this research.} There is no foreseeable population who may be put at a disadvantage from this research.


\textbf{What are the consequences of failure of the system.} 

Overreliance on any prediction method can have obvious, negative real-world consequences, particularly when the prediction method is prone to failure. CART suffers from a couple of pitfalls: instability (i.e., small perturbations in the training samples may significantly change the structure of an optimal tree and consequent predictions) and difficulty in accurately approximating certain simple models, such as linear or, more generally, additive, if given insufficient or low quality data. 


\textbf{Whether the task/method leverages biases in the data.} 
While CART is not impervious to all pre-existing biases in the data (e.g., those arising from systematic measurement errors at the data collection stage), as we have shown, it is less susceptible to the presence of additional, non-diagnostic variables in the data. 

Consequently, CART has the potential to mitigate the negative consequences of biasing information that is inevitable with most datasets.


\begin{ack}
The author is indebted to Min Xu, Minge Xie, Samory Kpotufe, Robert McCulloch, Andy Liaw, Richard Baumgartner, Regina Liu, Cun-Hui Zhang, Michael Kosorok, Jianqing Fan, and Matias Cattaneo for helpful discussions and feedback. This research was supported in part by NSF DMS $\#$1915932 and NSF TRIPODS DATA-INSPIRE CCF $\#$1934924.
\end{ack}


\small
\bibliographystyle{plainnat}
\bibliography{cart.bib}

\end{document}
{"config": "arxiv", "file": "2006.04266/main.tex"}