# Fast Linear Model Trees by PILOT

Jakob Raymaekers
Department of Quantitative Economics, Maastricht University, The Netherlands

Peter J. Rousseeuw
Section of Statistics and Data Science, University of Leuven, Belgium

Tim Verdonck
Department of Mathematics, University of Antwerp - imec, Belgium

Ruicong Yao
Section of Statistics and Data Science, University of Leuven, Belgium

(February 7, 2023)

###### Abstract

Linear model trees are regression trees that incorporate linear models in the
leaf nodes. This preserves the intuitive interpretation of decision trees and
at the same time enables them to better capture linear relationships, which is
hard for standard decision trees. But most existing methods for fitting linear
model trees are time consuming and therefore not scalable to large data sets.
In addition, they are more prone to overfitting and extrapolation issues than
standard regression trees. In this paper we introduce PILOT, a new algorithm
for linear model trees that is fast, regularized, stable and interpretable.
PILOT trains in a greedy fashion like classic regression trees, but
incorporates an $L^{2}$ boosting approach and a model selection rule for
fitting linear models in the nodes. The abbreviation PILOT stands for
PIecewise Linear Organic Tree, where ‘organic’ refers to the fact that no
pruning is carried out. PILOT has the same low time and space complexity as
CART without its pruning. An empirical study indicates that PILOT tends to
outperform standard decision trees and other linear model trees on a variety
of data sets. Moreover, we prove its consistency in an additive model setting
under weak assumptions. When the data is generated by a linear model, the
convergence rate is polynomial.

Keywords: Consistency, Piecewise linear model, Regression trees, Scalable
algorithms.

##  1 Introduction

Despite their long history, decision trees such as CART (Breiman et al., 1984)
and C4.5 (Quinlan, 1993) remain popular machine learning tools. Decision trees
can be trained quickly with few parameters to tune, and the final model can be
easily interpreted and visualized which is an appealing advantage in practice.
As a result, these methods are widely applied in a variety of disciplines
including engineering (Shamshirband et al., 2020), bioinformatics (Stiglic et
al., 2012), agriculture (Boryan et al., 2011), and business analysis (Delen et
al., 2013; Noulas et al., 2012). In addition to their standalone use, decision
trees have seen widespread adoption in ensemble methods, often as the best
available “weak learner”. Prime examples are random forests (Breiman, 2001)
and gradient boosting methods such as XGBoost (Chen and Guestrin, 2016) and
LightGBM (Ke et al., 2017). In this work we assume that the target variable is
continuous, so we focus on regression trees rather than on classification.

A limitation of classical regression trees is that their piecewise constant
nature makes them ill-suited to capture continuous relationships. They require
many splits in order to approximate linear functions, which is undesirable.
There are two main approaches to overcome this issue. The first is to use
ensembles of decision trees such as random forests or gradient boosted trees.
These ensembles smooth the prediction and can therefore model a more
continuous relation between predictors and response. A drawback of these
ensembles is the loss of interpretability. Combining multiple regression trees
no longer allows for a simple visualization of the model, or for an
explainable stepwise path from predictors to prediction. The second approach
to capture continuous relationships is to use model trees. Model trees have a
tree-like structure, but allow for non-constant fits in the leaf nodes of the
tree. Model trees thus retain the intuitive interpretability of classical
regression trees while being more flexible. Arguably the most intuitive and
common model tree is the linear model tree, which allows for linear models in
the leaf nodes.

To the best of our knowledge, the first linear model tree algorithm was
introduced by Friedman (1979). We will abbreviate it as FRIED. In each node it
uses univariate piecewise linear fits, to replace the piecewise constant fits
in CART. Once a model is fit on a node, its residuals are passed on to its
child nodes for further fitting. The final prediction is given by the sum of
the linear models along the path. Our experiments suggest that this method may
suffer from overfitting and extrapolation issues. The FRIED algorithm received
much less attention than its more involved successor MARS (Friedman, 1991).

The M5 algorithm (Quinlan, 1992) is by far the most popular linear model tree,
and is still commonly used today (Agrawal and Choudhary, 2016; Golbeck et al.,
2011; Noulas et al., 2012; Quercia et al., 2011; Shamshirband et al., 2020).
It starts by fitting CART. Once this tree has been built, linear regressions
are introduced at the leaf nodes. Pruning and smoothing are then applied to
reduce its generalization error. One potential objection against this
algorithm is that the tree structure is built completely oblivious of the fact
that linear models will be used in the leaf nodes.

The GUIDE algorithm (Loh, 2002) fits multiple linear models to numerical
predictors in each node and then applies a $\chi^{2}$ test comparing positive
and negative residuals to decide on which predictor to split. The algorithm
can also be applied to detect interactions between predictors. However, no
theoretical guarantee is provided to justify the splitting procedure. There
are also algorithms that use ideas from clustering in their splitting rule.
For example, SECRET (Dobra and Gehrke, 2002) uses the EM algorithm to fit two
Gaussian clusters to the data, and locally transforms the regression problem
into a classification problem based on the closeness to these clusters.
Experimental results do not favor this method over GUIDE, and its
computational cost is high.

The SMOTI method (Malerba et al., 2004) uses two types of nodes: regression
nodes and splitting nodes. In each leaf, the final model is the multiple
regression fit to its ‘active’ variables. This is achieved by ‘regressing out’
the variable of a regression node from both the response and the other
variables. The resulting fit is quite attractive, but the algorithm has a
complexity of $\mathcal{O}(n^{2}p^{3})$, which can be prohibitive for large
data sets. This makes it less suitable for an extension to random forests or
gradient boosting.

Most of the above algorithms need a pruning procedure to ensure their
generalization ability, which is typically time consuming. An exception is
LLRT (Vogel et al., 2007) which uses stepwise regression and evaluates the
models in each node via $k$-fold cross validation. To alleviate the
computational cost, the algorithm uses quantiles of the predictors as the
potential splitting points. It also maintains the data matrices for fitting
linear models on both child nodes so that they can be updated. Unfortunately,
the time complexity of LLRT is quite high at $\mathcal{O}(knp^{3}+np^{5})$
where $k$ is the depth of the tree.

It is also interesting to compare ensembles of linear model trees with
ensembles of classical decision trees. Recently, Shi et al. (2019) combined
piecewise linear trees with gradient boosting. For each tree, they used
additive fitting as in Friedman (1979) or half additive fitting where the past
prediction was multiplied by a factor. To control the complexity of the tree,
they set a tuning parameter for the maximal number of predictors that can be
used along each path. Empirical results suggest that this procedure
outperformed classical decision trees in XGBoost (Chen and Guestrin, 2016) and
in LightGBM (Ke et al., 2017) on a variety of data sets, and required fewer
training iterations. This points to a potential strength of linear model trees
for ensembles.

To conclude, the main issue with existing linear model trees is their high
computational cost. Methods that apply multiple regression fits to leaf and/or
internal nodes introduce a factor $p^{2}$ or $p^{3}$ in their time complexity.
And the methods that use simple linear fits, such as FRIED, still require
pruning which is also costly on large data sets. In addition, these methods
can have large extrapolation errors (Loh et al., 2007). Finally, to the best
of our knowledge there is no theoretical support for linear model trees in the
literature.

In response to this challenging combination of issues we propose a novel
linear model tree algorithm. Its acronym PILOT stands for PIecewise Linear
Organic Tree, where ‘organic’ refers to the fact that no pruning is carried
out. The main features of PILOT are:

  * •

Speed: It has the same low time complexity as CART without its pruning.

  * •

Regularized: In each node, a model selection procedure is applied to the
potential linear models. This requires no extra computational complexity.

  * •

Explainable: Thanks to the simplicity of linear models in the leaf nodes, the
final tree remains highly interpretable. Also a measure of feature importance
can be computed.

  * •

Stable extrapolation: Two truncation procedures are applied to avoid extreme
fits on the training data and large extrapolation errors on the test data,
which are common issues in linear model trees.

  * •

Theoretically supported: PILOT has proven consistency in additive models. When
the data is generated by a linear model, PILOT attains a polynomial
convergence rate, in contrast with CART.

The paper is organized as follows. Section 2 describes the PILOT algorithm and
discusses its properties, including a derivation of its time and space
complexity. Section 3 presents the two main theoretical results. First, the
consistency for general additive models is discussed and proven. Second, when
the true underlying function is indeed linear an improved rate of convergence
is demonstrated. We refer to the Appendix for proofs of the theorems,
propositions and lemmas. Empirical evaluations are provided in Section 4,
where PILOT is compared with several alternatives on a variety of benchmark
data sets. It outperformed other tree-based methods on data sets where linear
models are known to fit well, and outperformed other linear model trees on
data sets where CART typically performs well. Section 5 concludes.

##  2 Methodology

In this section we describe the workings of the PILOT learning algorithm. We
begin by explaining how its tree is built and motivate the choices made, and
then derive its computational cost.

We will denote the $n\times p$ design matrix as
$\boldsymbol{X}=({X}_{1},\dots,{X}_{n})^{\top}$ and the $n\times 1$ response
vector as $Y=(y_{1},\dots,y_{n})^{\top}$. We consider the standard regression
model $y=f(X)+\epsilon$, where $f:\mathbb{R}^{p}\rightarrow\mathbb{R}$ is the
unknown regression function and $\epsilon$ has mean zero.

As a side remark, the linear models and trees we consider here are equivariant
to adding a constant to the response, that is, the predictive models would
keep the same fitted parameters except for the intercept terms. Therefore, in
the presentation we will assume that the responses are centered around zero,
that is, $Y_{min}=-Y_{max}$ without loss of generality.

###  2.1 Main structure of PILOT

A typical regression tree has four ingredients: a construction rule, the
evaluation of the models/splits, a stopping rule, and the prediction. Most
regression tree algorithms are built greedily from top to bottom. That is,
they split the original space along a predictor and repeat the procedure on
the subsets. This approach has some major advantages. The first is its speed.
The second is that it starts by investigating the data from a global point of
view, which is helpful when detecting linear relationships.

Algorithm 1 presents a high-level description of PILOT. It starts by sorting
each predictor and storing its ranks. At each node, PILOT then selects a
predictor and a univariate piecewise linear model via a selection procedure
that we will describe in Sections 2.2 and 2.3. Then it fits the model and
passes the residuals on to the child nodes for further fitting. This recursion
is applied until the algorithm stops. There are three stopping triggers. In
particular, the recursion in a node stops when:

  * •

the number of cases in the node is below a threshold value $n_{\text{fit}}$;

  * •

the number of cases in one of the child nodes would be less than a threshold
$n_{\text{leaf}}$;

  * •

none of the candidate models do substantially better than a constant
prediction, as we will describe in Section 2.3.

When all nodes have reached a stopping trigger, the algorithm is done. The
final prediction is given by aggregating the predictions of the linear models
from the root to the leaves, as in the example in Figure 1.

Input: Training set $(\boldsymbol{X},Y)$, $n_{\text{fit}}$: minimal number of
cases to fit a model; $n_{\text{leaf}}$: minimal number of cases in a leaf;
$K_{\text{max}}$: maximal tree depth.

Output: $T$, a decision tree object.

function
_PILOT(_$\boldsymbol{X},Y,n_{\text{fit}},n_{\text{leaf}},K_{\text{max}}$_)_:

$I\leftarrow\\{1,\dots,\mbox{length}(Y)\\}$

$Y_{0}\leftarrow Y$; $p\leftarrow$ number of predictors

$B\leftarrow Y_{max}$

Sort the values $X^{(j)}_{i}$ in each predictor variable $X^{(j)}$, and store
the indices of the sorted cases in an $n\times p$ integer matrix $idx$.

Recursively build the tree using
$(\boldsymbol{X},Y),n_{\text{fit}},n_{\text{leaf}},K_{\text{max}},I,Y_{0},p,idx,B$,
and function Model_Selection. Fit the selected model, update the residuals
$Y$, and the index set $I$ for cases in each child node.

return _T_

end

Algorithm 1 Sketch of the PILOT algorithm

Figure 1: An example of a PILOT tree.

Note that PILOT retains good interpretability. The tree structure visualizes
the decision process, while the simple linear models in each step reveal which
predictor was used and how the prediction changes with respect to a small
change in this predictor. Moreover, one can define a measure of feature
importance similar to that of CART, based on the impurity gain of each node.
This is because PILOT selects only one predictor in each node, which makes the
gain fully dependent on it. This differs from methods such as M5, GUIDE, and
LLRT that use multiple predictors in each node, making it harder to fairly
distribute the gain over several predictors in order to derive an overall
measure of each predictor’s importance.

###  2.2 Models used in the nodes

As stated in the high-level summary in Algorithm 1, at each node PILOT selects
a fit from a number of linear and piecewise linear models. In particular,
PILOT considers the following regression models on each predictor:

  * •

con: A CONstant fit. We stop the recursion in a node after a con model is
fitted.

  * •

lin: Simple LINear regression.

  * •

pcon: A Piecewise CONstant fit, as in CART.

  * •

blin: A Broken LINear fit: a continuous function consisting of two linear
pieces.

  * •

plin: A two-Piece LINear fit that need not be continuous.

To guard against unstable fitting, the lin and blin models are only considered
when the number of unique values in the predictor is at least 5. Similarly,
the plin model is only considered if both potential child nodes have at least
5 unique values of the predictor.

These models extend the piecewise constant pcon fits of CART to linear fits.
In particular, plin is a direct extension of pcon while blin can be regarded
as a regularized and smoothed version of plin. For variables that are
categorical or only have a few unique values, only pcon and con are used.
Finally, if con is selected in a node, the recursion in that node stops.

PILOT reduces to CART (without its pruning) if only pcon models are selected,
since the least squares constant fit in a child node equals its average
response.

Note that a node will not be split when lin is selected. Therefore, when a lin
fit is carried out in a node, PILOT does not increase its reported depth. This
affects the reported depth of the final tree, which is a tuning parameter of
the method.

It is possible that a consecutive series of lin fits is made in the same node,
and this did happen in our empirical studies. In a node where this occurs,
PILOT is in fact executing $L^{2}$ boosting (Bühlmann, 2006), which fits
multiple regression using repeated simple linear regressions. It has been
shown that $L^{2}$ boosting is consistent for high dimensional linear
regression and produces results comparable to the Lasso. It was also shown in
Freund et al. (2017) that its convergence rate can be relatively fast under
certain assumptions. PILOT does not increment the depth value of the node for
lin models to avoid interrupting this boosting procedure.

###  2.3 Model selection rule

In each node, PILOT employs a selection procedure to choose between the five
model types con, lin, pcon, blin and plin. It would be inappropriate to select
the model based on the largest reduction in the residual sum of squares
($\mathrm{RSS}$), since this would always choose plin as that model
generalizes all the others. Therefore, we need some kind of regularization to
select a simpler model when the extra $\mathrm{RSS}$ gain of going all the way
to plin is not substantial enough.

After many experiments, the following regularization scheme was adopted. In
each node, PILOT chooses the combination of a predictor and a regression model
that leads to the lowest BIC value

$n\log\Big{(}\frac{\mathrm{RSS}}{n}\Big{)}+v\log(n)$ (1)

which is a function of the residual sum of squares $\mathrm{RSS}$, the number
of cases $n$ in the node, and the degrees of freedom $v$ of the model. With
this selection rule PILOT mitigates overfitting locally, and applying it
throughout the tree leads to a global regularization effect.

It remains to be addressed what the value of $v$ in (1) should be for each of
the five models. The degrees of freedom are determined by aggregating the
number of model parameters excluding the splits, and the degrees of freedom of
a split. The model types con, lin, pcon, blin and plin contain 1, 2, 2, 3, and
4 coefficients apart from any split points. Based on empirical and theoretical
arguments, Hall et al. (2017) found that a discontinuity point deserves 3
degrees of freedom. We follow this recommendation for pcon and plin, which
each receive 3 additional degrees of freedom. Also blin contains a split
point, but here the fitted function is continuous. To reflect this
intermediate complexity we add 2 degrees of freedom to blin. In conclusion, we
end up with $v$ = 1, 2, 5, 5, and 7 for the model types con, lin, pcon, blin
and plin.

The BIC in (1) is one of several model selection criteria that we could have
used used. We also tried other measures, such as the Akaike information
criterion (AIC) and the adjusted AIC. It turned out that the AIC tended to
choose plin too often, which reduced the regularization effect. The adjusted
AIC required a pruning procedure to perform comparably to the BIC criterion.
As pruning comes at a substantial additional computational cost, we decided
against the adjusted AIC.

Alternatively, one could follow Bühlmann (2006) to compute the degrees of
freedom for the aggregated model in each node, and then compute the adjusted
AIC to decide when to stop. But in this approach hat matrices have to be
maintained for each path, which requires more computation time and memory
space. Moreover, as the number of cases in the nodes changes, the evaluations
of the hat matrices become more complicated. For these reasons, we preferred
to stay with the selection rule (1) in each node.

Algorithm 2 shows the pseudocode of our model selection procedure in a given
node. We iterate through each predictor separately. If the predictor is
numerical we find the univariate model with lowest BIC among the 5 candidate
models. For con and lin this evaluation is immediate. For the other three
models, we need to consider all possible split points. We consider these in
increasing order (recall that the ordering of each predictor was already
obtained at the start of the algorithm), which allows the Gram and moment
matrices for the linear models to be maintained and updated efficiently in
each step. On the other hand, if the predictor is categorical, we follow the
approach of CART as in Section 9.2.4 of Hastie et al. (2009). That is, we
first compute the mean $m_{c}$ of the response of each level $c$ and then sort
the cases according to $m_{c}$. We then fit the model pcon in the usual way.

Once the best model type has been determined for each predictor separately,
Algorithm 2 returns the combination of predictor and model with the lowest BIC
criterion.

Globals : $(\boldsymbol{X},Y)$: predictors and residual vector in the current
node; $p$: number of predictors; $idx$: indices matrix of the $p$ presorted
predictors; $n_{\text{leaf}}$: minimal number of cases in the leaf.

Input: $I$, indices of cases in the current node.

Output: $best\\_pivot$; $best\\_model$; $best\\_pred$; $lm\\_l$, $lm\\_r$,
linear model for the left and right child; $range$, range of the selected
predictors in the current node.

function  _Model_Selection(_$I$_)_:

init:  _best_pivot, best_model, best_bic, best_pred, lm_l, lm_r, range_

for  _each predictor $X^{(j)}$ with $j=1,\ldots,p$_ do

Select the cases sorted according to $X^{(j)}$ using $I$ and $idx$.

Fit con and initialize the parameters of the best model using its result.

if  _$X^{(j)}$ is a numerical predictor_ then

Compute $(D^{\textsc{lin}})^{\top}D^{\textsc{lin}}$ and
$(D^{\textsc{lin}})^{\top}Y$, where $D^{\textsc{lin}}$ is the design matrix of
the linear model. It is $|I|\times 1$ for pcon, and $|I|\times 2$ otherwise.

Fit lin using the above matrices. Store the parameters of the best model.

For pcon and plin initialize the Gram matrices $G^{\textsc{plin}}_{l,r}$ and
moment matrices $M^{\textsc{plin}}_{l,r}$ of the simple linear models in the
left and right child such that
$G^{\textsc{plin}}_{l}=M^{\textsc{plin}}_{l}=\boldsymbol{0}$,
$G^{\textsc{plin}}_{r}=(D^{\textsc{lin}})^{\top}D^{\textsc{lin}}$,
$M^{\textsc{plin}}_{r}=(D^{\textsc{lin}})^{\top}Y$.

For blin initialize the Gram matrix
$G^{\textsc{blin}}=(D^{\textsc{blin}})^{\top}D^{\textsc{blin}}$ and the moment
matrix $M^{\textsc{blin}}=(D^{\textsc{blin}})^{\top}Y$ where
$D^{\textsc{blin}}$ is the design matrix of blin with knot equal to $\min
X^{(j)}$.

for  _each possible splitting pivot_ do

Update $G^{\textsc{blin}},M^{\textsc{blin}}$,
$G^{\textsc{plin}}_{l,r},M^{\textsc{plin}}_{l,r}$.

If $n_{\text{leaf}}$ is satisfied, fit plin, pcon and blin using the Gram and
moment matrices, yielding the left and right models $lm\\_l$ and $lm\\_r$.

Store the parameters of the best model.

end for

else

Compute the mean $m_{c}$ of the response for each categorical level $c$.

Re-sort the cases according to their $m_{c}$.

Greedily search the best $value$ such that pcon with partition
{$c|m_{c}\leqslant value$} and {$c|m_{c}>value$} has the lowest BIC.

Store the parameters of the best model. Pivot is the set {$c|m_{c}\leqslant
value$}.

end for

return:  _best_pivot, best_model, best_bic, best_pred, lm_l, lm_r, range_

end

Algorithm 2 Model selection and split finding in a node of the training set

###  2.4 Truncation of predictions

PILOT relies on two truncation procedures to avoid unstable predictions on the
training data and the test data.

Figure 2: Left: An example of the first truncation method in one node of the
tree. Right: An example of the second truncation method.

The first truncation procedure is motivated as follows. The left panel of
Figure 2 shows the data in a node of PILOT. The selected predictor is on the
horizontal axis. The plot illustrates a situation where most of the data
points in a node are concentrated in a part of the predictor range, causing
the linear model to put a lot of weight on these points and little weight on
the other(s). This can result in extreme predictions, for instance for the red
point on the right. Although such unstable predictions might be corrected by
deeper nodes, this could induce unnecessarily complicated models with extra
splits. Moreover, if this effect occurs at a leaf node, the prediction cannot
be corrected. To avoid this issue we will truncate the prediction function.
Note that if we add a constant to the response, PILOT will yield the same
estimated parameters except for the intercept terms that change the way they
should. Therefore, we can assume without loss of generality that
$Y_{max}=-Y_{min}$ . Calling this quantity $B$, the original responses have
the range $[-B,B]$. We then clip the prediction so it belongs to $[-3B,3B]$.
Then we compute the response $Y$ minus the truncated prediction to get the new
residuals, and proceed with building the tree for the next depth. The whole
process is summarized in Algorithm 3. Truncating at $\pm 3B$ may seem rather
loose, and one could truncate at a smaller multiple of $B$ if so desired. But
we found empirically that truncating at $\pm 3B$ works fine for our purpose.

Globals : $p$; $idx$; $n_{\text{fit}}$: minimal number of cases required to
fit a model; $n_{\text{leaf}}$: minimal number of cases in a leaf;
$K_{\text{max}}$: maximal depth; $Y_{0}$: original response; $B$: maximum of
centered response.

Input: $(\boldsymbol{X},Y)$: Training set; $K$: current depth; $I$: indices of
cases.

Output: A tree object $T$ with attributes $m$: the model; $l,r$: the left and
the right child node; $lm_{l},lm_{r}$: model coefficients in the left and
right child nodes.

function  _Build_Tree(_$(\boldsymbol{X},Y),K,I$_)_:

$T\leftarrow$ a None node

if  _$K <K_{max}\text{ and length($I$)}\geqslant n_{\text{fit}}$_ then

Do Model_Selection($I$) and record the parameters for the best model

$K\leftarrow K+1$ if the $best\\_model$ is not lin or con

$T\leftarrow$ A tree node with the recorded parameters

if  _best_model is con_ then

return _$T$_

else if  _best_model is lin_ then

$raw\\_pred[I]\;\leftarrow\;$ the prediction of $lm_{l}$ on
$\boldsymbol{X}[I]$

Truncate the actual prediction $Y_{0}[I]-Y[I]+raw\\_pred[I]$ at $\pm 3B$.
Update the residuals $Y[I]$ with this new prediction.

$T.l\leftarrow\texttt{Build\\_Tree}((\boldsymbol{X},Y),K,I)$

else

$I_{l}$, $I_{r}\leftarrow$ the indices of cases in the left and right child

$raw\\_pred[I_{l}]\;\leftarrow\;$ the prediction of $lm_{l}$ on
$\boldsymbol{X}[I_{l}]$

$raw\\_pred[I_{r}]\;\leftarrow\;$ the prediction of $lm_{r}$ on
$\boldsymbol{X}[I_{r}]$

Truncate the actual predictions $Y_{0}[I_{l}]-Y[I_{l}]+raw\\_pred[I_{l}]$ and
$Y_{0}[I_{r}]-Y[I_{r}]+raw\\_pred[I_{r}]$ at $\pm 3B$, and update the
residuals $Y[I_{l}]$ and $Y[I_{r}]$ in the child nodes.

$T.l\leftarrow\texttt{Build\\_Tree}((\boldsymbol{X},Y),K,I_{l})$

$T.r\leftarrow\texttt{Build\\_Tree}((\boldsymbol{X},Y),K,I_{r})$

end if

else

return _$T$_

end if

return $T$

Algorithm 3 Truncation during tree building on the training data

The first truncation method is not only applied during training, but also to
the prediction on new data. The range of the response of the training data is
stored, and when a prediction for new data would be outside the proposed range
$[-3B,3B]$ it is truncated in the same way. This procedure works as a basic
safeguard for the prediction. However, this is not enough because a new issue
can arise: it may happen that values of a predictor in the test data fall
outside the range of the same predictor in the training data. This is
illustrated in the right panel of Figure 2. The vertical dashed lines indicate
the range of the predictor on the training data. However, the predictor value
of the test case (the red point) lies quite far from that range. Predicting
its response by the straight line that was fitted during training would be an
extrapolation, which could be unrealistic. This is a known problem of linear
model trees. For instance, in our experiments in Section 4 the predictions of
FRIED and M5 exploded on some data sets, which turned out to be due to this
problem.

Therefore we add a second truncation procedure, which is only applied to new
data (because in the training data, all predictor values are inside their
range). Therefore, we propose to truncate the prediction value as done by Loh
et al. (2007). More precisely, during training we record the range
$[x_{\mbox{\tiny min}},x_{\mbox{\tiny max}}]$ of the predictor selected in
this node, and store the range of the corresponding predictions
$[\hat{f}(x_{\mbox{\tiny min}}),\hat{f}(x_{\mbox{\tiny max}})]$. When
predicting on new data, we truncate the predictions so they stay in the
training range. More precisely, we replace the original predictions
$\hat{f}(x_{\mbox{\tiny test}})$ on this node by
$\max(\min(\hat{f}(x_{\mbox{\tiny test}}),\hat{f}(x_{\mbox{\tiny
max}})),\hat{f}(x_{\mbox{\tiny min}}))$.

The two truncation methods complement each other. For instance, the first
approach would not suffice in the right panel of Figure 2 since the linear
prediction at the new case would still lie in $[-3B,3B]$. Also, the second
approach would not work in the left panel of Figure 2 as the unusual case is
included in the range of the training data. Our empirical studies indicate
that combining both truncation approaches helps eliminate extreme predictions
and therefore improve the stability of our model. Moreover, this did not
detract from the overall predictive power. Algorithm 4 describes how both
truncation methods are applied when PILOT makes predictions on new data.

Input: $T$: a decision tree object; $x$: a new observation; $B$: maximum value
of the centered response in the training data.

Output: $pred$, the prediction in $x$.

function  _Predict(_$T,X$_)_:

$pred\leftarrow 0$

while  _T is not None_ do

if  _$x$ is in the left child of $T$_ then

$X_{min},X_{max}\leftarrow$ The lower and upper bound of $T.range$ on the
selected predictor $j$ in the training data.

$X^{*}\leftarrow\min(\max(X^{(j)},X_{min}),X_{max})$

$pred\\_add$ $\leftarrow$ The prediction of $T.lm\\_l$ on $X^{*}$

$pred$ $\leftarrow$ $\max(\min(pred+pred\\_add,3B)-3B)$

$T\leftarrow T.l$

else

$X_{min},X_{max}\leftarrow$ The lower and upper bound of $T.range$ on the
selected predictor $j$ in the training data.

$X^{*}\leftarrow\max(\min(X^{(j)},X_{min}),X_{max})$

$pred\\_add$ $\leftarrow$ The prediction of $T.lm\\_r$ on $X^{*}$

$pred$ $\leftarrow$ $\max(\min(pred+pred\\_add,3B),-3B)$

$T\leftarrow T.r$

end if

end while

return _$pred$_

end

Algorithm 4 Truncation during prediction on new data

###  2.5 Stopping rules versus pruning

It is well known that decision trees have a tendency to overfit the data if
the tree is allowed to become very deep. A first step toward avoiding this is
to require a minimal number of cases in a node before it can be split, and a
minimal number of cases in each child node. In addition, PILOT also stops
splitting a node if the BIC model selection criterion selects con.

Several decision trees, including CART, take a different approach. They let
the tree grow, and afterward prune it. This indeed helps to achieve a lower
generalization error. However, pruning can be very time consuming. One often
needs a cross-validation procedure to select pruning parameters. It would be
possible for PILOT to also incorporate cost-complexity pruning as in CART. (If
a series of lin fits occurs in the same node, PILOT would either keep all of
them or remove them together.) However, in our empirical study we found that
adding this pruning to PILOT did not outperform the existing con stopping rule
on a variety of datasets. Also, we will see in the next section that the
generalization error vanishes asymptotically. Moreover, not pruning is much
more efficient computationally. For these reasons, PILOT does not employ
pruning. The computational gain of this choice makes it more feasible to use
PILOT in ensemble methods.

###  2.6 Time and space complexity

Finally, we demonstrate that the time complexity and space complexity of PILOT
are the same as for CART without its pruning. The main point of the proof is
that for a single predictor, each of the five models can be evaluated in a
single pass through the predictor values. Evaluating the next split point
requires only a rank-one update of the Gram and moment matrices, which can be
done in constant time. Naturally, the implicit proportionality factors in
front of the $\mathcal{O}(.)$ complexities are higher for PILOT than for CART,
but the algorithm is quite fast in practice.

######  Proposition 1.

PILOT has the same time and space complexities as CART without its pruning.

###### Proof.

For both PILOT and CART we assume that the $p$ predictors have been presorted,
which only takes $\mathcal{O}(np\log(n))$ time once.

We first check the complexity of the model selection in Algorithm 2. It is
known that the time complexity of CART for split finding is $\mathcal{O}(np)$.
For the pcon model, PILOT uses the same algorithm as CART. For con we only
need to compute one average in $\mathcal{O}(n)$ time. For lin, the evaluation
of $(D^{\textsc{lin}})^{\top}D^{\textsc{lin}}$ and
$(D^{\textsc{lin}})^{\top}Y$ also takes $\mathcal{O}(n)$. In the evaluation of
the plin model, the Gram and moment matrices always satisfy
$G^{\textsc{plin}}_{l}+G^{\textsc{plin}}_{r}=G^{\textsc{lin}}=(D^{\textsc{lin}})^{\top}D^{\textsc{lin}}$
and $M^{\textsc{plin}}_{l}+M^{\textsc{plin}}_{r}=(D^{\textsc{lin}})^{\top}Y$
by definition. These matrices have at most 4 entries, and can be incrementally
updated in $\mathcal{O}(1)$ time as we evaluate the next split point. For blin
the reasoning is analogous. In each model, inverting the Gram matrix only
takes $\mathcal{O}(1)$ because its size is fixed. Therefore, we can evaluate
all options on one presorted predictor in one pass through the data, which
remains $\mathcal{O}(n)$, so for all predictors this takes $\mathcal{O}(np)$
time, the same as for CART. For the space complexity, the Gram and moment
matrices only require $\mathcal{O}(1)$ of storage. Since CART also has to
store the average response in each child node, which takes $\mathcal{O}(1)$
storage as well, the space complexity of both methods is the same.

For the tree building in Algorithm 3, the time complexity of computing the
indices $I_{l,r}$ and the residuals is $\mathcal{O}(n)$ at each depth
$1\leqslant k\leqslant K_{max}$. CART also computes the indices $I_{l,r}$
which requires the same complexity $\mathcal{O}(n)$ for each step. Therefore,
the overall time for both methods is $\mathcal{O}(K_{max}n)$. Since the PILOT
tree only has $\mathcal{O}(1)$ more attributes in each node than CART, the
space complexity of PILOT remains the same as that of CART.

For the initialization and prediction parts, the results are straightforward.
∎

##  3 Theoretical results

###  3.1 Universal consistency

In this section we prove the consistency of PILOT. We follow Klusowski (2021)
and assume the underlying function $f\in\mathcal{F}\subset L^{2}([0,1]^{p})$
admits an additive form

$f(X):=f_{1}(X^{(1)})+\dots+f_{p}(X^{(p)})$ (2)

where $f_{i}$ has bounded variation and $X^{(j)}$ is the $j$-th predictor. We
define the total variation norm $||f||_{TV}$ of $f\in\mathcal{F}$ as the
infimum of $\sum_{i=1}^{p}||f_{i}||_{TV}$ over all possible representations of
$f$, and assume that the representation in (2) attains this infimum.

For all $f,g\in L^{2}([0,1]^{p})$ and a dataset $X_{1},\ldots,X_{n}$ of size
$n$, we define the empirical norm and the empirical inner product as

$||f||_{n}^{2}\,:=\frac{1}{n}\sum_{i=1}^{n}|f(X_{i})|^{2}\quad\text{and}\quad\langle
f,g\rangle_{n}:=\frac{1}{n}\sum_{i=1}^{n}f(X_{i})g(X_{i}).$ (3)

For the response vector $Y=(Y_{1},\ldots,Y_{n})$ we denote

$||Y-f||_{n}^{2}\,:=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-f(X_{i}))^{2}\quad\text{and}\quad\langle
Y,f\rangle_{n}:=\frac{1}{n}\sum_{i=1}^{n}Y_{i}f(X_{i}).$ (4)

The general $L^{2}$ norm on vectors and functions will be denoted as
$||\cdot||$ without subscript.

To indicate the norm and the inner product of a function on a node $T$ with t
observations we replace the subscript $n$ by $t$. We denote by
$\mathcal{T}_{k}$ the set of tree nodes at depth $k$, plus the leaf nodes of
depth lower than $k$. In particular, $\mathcal{T}_{K}$ contains all the leaf
nodes of a $K$-depth tree.

Now we can state the main theorem:

######  Theorem 2.

Let $f\in\mathcal{F}$ with $||f||_{TV}\,\leqslant A$ and denote by
$\hat{f}(\mathcal{T}_{K})$ the prediction of a $K$-depth PILOT tree. Suppose
$X\sim P$ on $[0,1]^{p_{n}}$, the responses are a.s. bounded in $[-B,B]$, and
the depth $K_{n}$ and the number of predictors $p_{n}$ satisfy
$K_{n}\rightarrow\infty$ and $2^{K_{n}}p_{n}\log(np_{n})/n\rightarrow 0$. Then
PILOT is consistent, that is

$\lim_{n\rightarrow\infty}\mathbb{E}[||f-\hat{f}(\mathcal{T}_{K_{n}})||^{2}]=0.$
(5)

######  Remark 3.

Note that the conditions on the depth and the dimension can easily be
satisfied if we let $K_{n}=\log_{2}(n)/r$ for some $r>1$ and $p_{n}=n^{s}$
such that $0<s<1-1/r$. The resulting convergence rate is
$\mathcal{O}(1/\log(n))$. This is the same rate as obtained by Klusowski
(2021) for CART under very similar assumptions.

The key of the proof of Theorem 2 is to establish a recursive formula for the
training error $R_{k}:=||Y-\hat{f}(\mathcal{T}_{k})||^{2}_{n}-||Y-f||^{2}_{n}$
for each depth $k$. Then we can leverage results from empirical process theory
(Györfi et al., 2002) to derive an oracle inequality for PILOT, and finally
prove consistency.

For the recursive formula we first note that
$R_{k}=\sum_{T\in\mathcal{T}_{K}}w(T)R_{k}(T)$ where $w(T):=t/n$ is the weight
of node $T$ and
$R_{k}(T):=||Y-\hat{f}(\mathcal{T}_{k})||^{2}_{t}-||Y-f||^{2}_{t}$. Then we
immediately have $R_{k+1}=R_{k}-\sum_{T\in\mathcal{T}_{k}}w(T)\Delta^{k+1}(T)$
where

$\Delta^{k+1}(T):=||Y-\hat{f}(\mathcal{T}_{k})||^{2}_{t}\;-\;t_{l}||Y-\hat{f}(\mathcal{T}_{k+1})||^{2}_{t_{l}}/t\;-\;t_{r}||Y-\hat{f}(\mathcal{T}_{k+1})||^{2}_{t_{r}}/t$

is the impurity gain of the model on $T$. Here $t_{l}$ and $t_{r}$ denote the
number of cases in the left and right child node. If no split occurs, one of
these numbers is set to zero. For PILOT we need to remove con nodes that do
not produce an impurity gain from the recursive formula. To do this, we define
$C_{k}^{+}=\\{T|T=\textsc{con},T\in\mathcal{T}_{k-1},R_{k}(T)>0\\}$, i.e., the
set of ‘bad’ nodes on which con is fitted, and similarly $C_{k}^{-}$ for those
with $R_{k}(T)\leqslant 0$. Then we can define $R_{C_{k}^{+}}:=\sum_{T\in
C_{k}^{+}}w(T)R_{k}(T)$, the positive training error corresponding to the
cases in these con nodes. We can then consider
$\widetilde{R}_{k}:=R_{k}-R_{C_{k}^{+}}$ and show that asymptotically both
$R_{k}$ and $R_{C_{k}^{+}}$ become small.

Now the problem is reduced to relating the impurity gain of the selected model
to the training error before this step. Recently, Klusowski (2021) introduced
a novel approach for this kind of estimation in CART. However, his proof used
the fact that in each step the prediction is piecewise constant, which is not
the case here. Therefore, we derived a generalized result for pcon in PILOT
(Lemma 12 in Appendix B.1). For controlling the gain of the other models, we
can use the following proposition:

######  Proposition 4.

Let $\Delta_{1}$, $\Delta_{2}$ and $v_{1}$, $v_{2}$ be the impurity gains and
degrees of freedom of two regression models on some node $T$ with $t$ cases.
Let $R_{0}$ be the initial residual sum of squares in $T$. We have that

  * •

If model 1 does better than con, i.e. $BIC_{con}>BIC_{1}$, we have
$t\Delta_{1}/R_{0}>C(v_{1},t)>0$ for some positive function $C$ depending on
$v_{1}$ and $t$.

  * •

If $BIC_{con}>BIC_{1}$ and $BIC_{2}>BIC_{1}$ with $v_{2}\geqslant v_{1}$ we
have $\frac{\Delta_{1}}{\Delta_{2}}\geqslant\frac{v_{1}-1}{v_{2}-1}.$

Moreover, if con is chosen at a node, it will also be chosen in subsequent
models.

We thus know that the gain of the selected model is always comparable to that
of pcon, up to a constant factor.

The proposition also justifies the con stopping rule from a model selection
point of view: when con is selected for a node, all the subsequent models for
this node will still be con. A con model is selected when the gain of other
regression models do not make up for the complexity they introduce. Since con
only regresses out the mean, we can stop the algorithm within the node once
the first con fit is encountered.

We can show the following recursive formula:

######  Theorem 5.

Under the same assumptions as Theorem 2 we have for any $K\geqslant 1$ and any
$f\in\mathcal{F}$ that

$||Y-\hat{f}(\mathcal{T}_{K})||^{2}_{n}\,\leqslant||Y-f||^{2}_{n}\,+\frac{(||f||_{TV}+6B)^{2}}{6(K+3)}+R_{C_{K}^{+}}\;.$
(6)

Moreover, if we let $K=\log_{2}(n)/r$ with $r>1$ we have
$R_{C_{K}^{+}}\sim\mathcal{O}(\sqrt{\log(n)/n^{(r-1)/r}})$.

The BIC criterion ensures that the training errors in all con nodes vanish as
$n\rightarrow\infty$. This allows us to prove the main theorem by using the
preceding results and Theorems 11.4 and 9.4 and Lemma 13.1 of Györfi et al.
(2002). All intermediate lemmas and proofs can be found in the Appendix B.

###  3.2 Convergence rates on linear models

The convergence rate resulting from Theorem 2 is the same as that of CART
derived by Klusowski (2021). This is because we make no specific assumptions
on the properties of the true underlying additive function. Therefore this
rate of convergence holds for a wide variety of functions, including poorly
behaved ones. In practice however, the convergence could be much faster if the
function is somehow well behaved. In particular, given that PILOT incorporates
linear models in its nodes, it is likely that it will perform better when the
true underlying function is linear, which is a special case of an additive
model. We will show that PILOT indeed has an improved rate of convergence in
that setting. This result does not hold for CART, and to the best of our
knowledge was not proved for any other linear model tree algorithm.

In order to prove the convergence rate of PILOT in linear models, we apply a
recent result of Freund et al. (2017). For the convergence rate of the $L^{2}$
boosting algorithm, Freund et al. (2017) showed in their Theorem 2.1 that the
mean squared error (MSE) decays exponentially with respect to the depth $K$,
for a fixed design matrix. The proof is based on the connection between
$L^{2}$ boosting and quadratic programming (QP) on the squared loss
$||Y-\boldsymbol{X}\beta||^{2}_{n}$. In fact, the $L^{\infty}$ norm of the
gradient of the squared loss is equal to the largest impurity gain among all
simple linear models, assuming each predictor has been centered and
standardized. Therefore, the gain of the $k$-th step depends linearly on the
mean squared error in step $k-1$, which results in a faster convergence rate.
By the nature of quadratic programming, the rate depends on the smallest
eigenvalue $\lambda_{min}$ of $\boldsymbol{X}^{\top}\boldsymbol{X}$.

In our framework we can consider a local QP problem in each node to estimate
the gain of lin, which gives a lower bound for the improvement of the selected
model (except when con is selected). To ensure that we have a constant ratio
for the error decay, we have to make some assumptions on the smallest
eigenvalue of the local correlation matrix of the predictors. This eigenvalue
can be regarded as a measure of multicollinearity between the predictors,
which is a key factor that affects the convergence rate of least squares
regression.

To be precise, we require that the smallest eigenvalue of the theoretical
correlation matrix (the $jl$-th entry of which is
$\mathbb{E}[(X-\bar{X})(X-\bar{X})^{\top}]_{jl}/(\sigma_{j}\sigma_{l})$ where
$\sigma_{j}$ is the standard deviation of the $j$-th predictor) is lower
bounded by some constant $\lambda_{0}$ in any cubic region that is
sufficiently small. Then we apply two concentration inequalities to show that
the smallest eigenvalue of a covariance matrix of data in such a cube is
larger than $\lambda_{0}$ with high probability. Finally, we can show that the
expectation of the error decays exponentially, which leads to a fast
convergence rate.

Our conditions are:

  * •

Condition 1: The PILOT algorithm stops splitting a node whenever

    * –

the number of cases in the node is less than $n_{\mbox{\tiny
min}}=n^{\alpha}$, $0<\alpha<1$.

    * –

the variance of some predictor is less than $2\sigma_{0}^{2}$ where
$0<\sigma_{0}<1$.

    * –

the volume of the node is less than a threshold $\eta$, $0<\eta<1$.

  * •

Condition 2: We assume that $X\sim P$ on $[0,1]^{p}$ and the error $\epsilon$
has a finite fourth moment. Moreover, for any cube $C$ with volume smaller
than $\eta$ we assume that$\lambda_{min}(Cor(X|X\in C))\geqslant
2\lambda_{0}>0$, where $Cor(X)$ is the correlation matrix.

The first condition is natural for tree-based methods. For the second one, an
immediate example would be a setting when all the predictors are independent.
Another example is when they follow a multivariate Gaussian distribution,
which becomes a truncated Gaussian when restricted to $[0,1]^{p}$ or to some
other cube $C$. It can be shown that the Gaussian distribution has relatively
small correlation in cubes $C$ of small volume; see, e.g., Muthen (1990) for
the two-dimensional case. This is not surprising since on small cubes the
density tends to a constant. As a result, the smallest eigenvalue of the
correlation matrix on such cubes will be bounded from below.

Under these conditions we obtain the following result, whose proof can be
found in Appendix C.

######  Theorem 6.

Assume the data are generated by the linear model
$Y\sim\boldsymbol{X}\beta+\epsilon$ with $n$ cases. Under conditions 1 and 2,
the difference between the training loss $L^{k}_{n}$ of PILOT at depth $k$ and
the training loss $L^{*}_{n}$ of least squares linear regression satisfies

$\mathbb{E}[L^{k}_{n}-L^{*}_{n}]\leqslant\gamma^{k}\sqrt{\mathbb{E}[(L^{0}_{n}-L^{*}_{n})^{2}]}+\mathcal{O}(N_{\text{leaves}}\log(n)/n)$
(7)

where

$\gamma:=1-\frac{\lambda_{0}}{4p}\;.$

Combining this result with the techniques in the proof of Theorem 2, we can
show polynomial rate convergence of our method on linear functions.

######  Corollary 7 (Fast convergence on linear models).

Assume the conditions of Theorem 6 hold and that $|Y|$ is a.s. bounded by some
constant $B$. Let $K_{n}=log_{\gamma}(n)$. Then we have for any $0<\alpha<1$
and corresponding $n_{\mbox{\tiny min}}=n^{\alpha}$ that

$\mathbb{E}[||\hat{f}(\mathcal{T}_{K})-\boldsymbol{X}\beta||^{2}]\leqslant\mathcal{O}\Big{(}\frac{\log(n)}{n^{\alpha}}\Big{)}\;.$
(8)

The choice of the tuning parameter $\alpha$ thus determines the convergence
rate. When we use a low $\alpha$ we get smaller nodes hence more nodes,
yielding looser bounds in the oracle inequality and for the errors in con
nodes, leading to a slower convergence. This is natural, as we have fewer
cases for the estimation in each node. If we choose $\alpha$ relatively high
we obtain a faster rate, which is intuitive since for $\alpha\rightarrow 1$ we
would have only a single node.

Cattaneo et al. (2022) prove a polynomial convergence rate for oblique
decision trees with a more general function class, and it may be possible to
adapt this work to PILOT. However, we note that their result makes the strong
assumption that the training cases should be approximately evenly distributed
over the leaf nodes. Our conditions here focus on the rules of the algorithm
and the distribution of the data itself.

##  4 Empirical evaluation

In this section we evaluate the proposed PILOT algorithm empirically and
compare its results with those of some popular competitors. We start by
describing the data and methods under comparison.

###  4.1 Data sets and methods

We analyzed 20 benchmark data sets, with the number of cases ranging from 71
to 21263 and the number of predictors ranging from 6 to 4088. From the UCI
repository (Dua and Graff, 2017) we used the data sets Abalone, Airfoil, Bike,
Communities, Concrete, Diabetes, Electricity, Energy, Residential, Skills,
Superconductor, Temperature and Wine. From Kaggle we obtained Graduate
Admission, Bodyfat, Boston Housing and Walmart. The California Housing data
came from the StatLib repository (http://lib.stat.cmu.edu/datasets/). The
Ozone data is in the R-package hdi (Dezeure et al., 2015), and Riboflavin came
from the R-package missMDA (Josse and Husson, 2016). Table 1 gives an
alphabetical list of the data sets with their sizes.

Table 1: The data sets used in the empirical study. Source | Dataset | $n$ | $p$
---|---|---|---
UCI repository | Abalone | 4177 | 8
Kaggle | Graduate Admission | 400 | 9
UCI repository | Airfoil | 1503 | 6
UCI repository | Bike | 17389 | 16
Kaggle | Bodyfat | 252 | 14
Kaggle | Boston Housing | 506 | 13
StatLib repository | California Housing | 20640 | 8
UCI repository | Communities | 1994 | 128
UCI repository | Concrete | 1030 | 9
UCI repository | Diabetes | 442 | 10
UCI repository | Electricity | 10000 | 14
UCI repository | Energy | 768 | 8
hdi | Ozone | 112 | 11
UCI repository | Residential | 372 | 105
missMDA | Riboflavin | 71 | 4088
UCI repository | Skills | 3395 | 20
UCI repository | Superconductor | 21263 | 81
UCI repository | Temperature | 7750 | 25
Kaggle | Walmart | 6435 | 8
UCI repository | Wine | 4898 | 12

In the comparative study we ran the following tree-based methods: the proposed
PILOT algorithm we implemented in Python with NUMBA acceleration in the split
function, FRIED (Friedman, 1979) implemented by us, M5 (Quinlan, 1992) from
the R package Rweka (Hornik et al., 2009), and CART. By default we let
$K_{\text{max}}=12$, $n_{\text{fit}}=10$, and $n_{\text{leaf}}=5$ which are
common choices for regression trees. We also ran ridge linear regression from
the Scikit-Learn package (Pedregosa et al., 2011) and the lasso (Tibshirani,
1996). When running methods that were not devised to deal with categorical
variables, we one-hot encoded such variables first.

We randomly split each dataset into 5 folds and computed the cross-validated
mean square error, that is, the average test MSE of each method on each of the
5 folds. The final score of each method is its MSE divided by the lowest MSE
on that dataset. A score of 1 thus corresponds with the method that performed
best on that dataset, and the scores of the other methods say how much larger
their MSE was.

###  4.2 Results

The results are summarized in Table 2. The blue entries correspond to scores
of at most 1.05, that is, methods whose MSE was no more than 1.05 times that
of the best MSE on that dataset. PILOT outperformed FRIED 16 times out of 19.
It also did better than CART on 17 datasets, and outperformed M5 on 15. In the
bottom two lines of the table we see that PILOT achieved the best average
score with the lowest standard deviation.

Table 2: MSE ratios relative to best, for the datasets of Table 1.

| PILOT | Ridge | LASSO | FRIED | CART | M5 |
---|---|---|---|---|---|---|---
Abalone | 1.00 | 1.02 | 1.01 | 1.11 | 1.12 | 1.56 |
Admission | 1.19 | 1.00 | 1.00 | 1.11 | 1.41 | (1.21) |
Airfoil | 1.98 | 4.40 | 4.40 | 1.00 | 1.75 | 9.52 |
Bike | 1.00 | 3.03 | 3.05 | 1.84 | 1.05 | 1.45 |
Bodyfat | 1.15 | 1.00 | 1.03 | 1.47 | 1.55 | (1.16) |
Boston | 1.02 | 1.00 | 1.06 | 1.17 | 1.16 | 3.99 |
California | 1.00 | 1.95 | 1.95 | 1.08 | 1.42 | (11.74) |
Communities | 1.08 | 1.01 | 1.00 | 1.17 | 1.32 | 1.04 |
Concrete | 1.00 | 2.61 | 2.62 | 3.42 | 1.38 | 2.48 |
Diabetes | 1.07 | 1.00 | 1.00 | 1.24 | 1.31 | 1.12 |
Electricity | 1.00 | 2.75 | 2.75 | 1.28 | 1.88 | 2.19 |
Energy | 1.17 | 3.60 | 3.57 | 1.00 | 1.29 | 1.12 |
Ozone | 1.36 | 1.05 | 1.00 | 1.37 | 1.18 | 1.22 |
Residential | 1.20 | 1.08 | 1.01 | 116.31 | 1.99 | 1.00 |
Riboflavin | 2.21 | 1.00 | 1.14 | 3.43 | 2.89 | 2.44 |
Skills | 1.00 | 1.02 | 1.03 | 1.06 | 1.17 | (1.06) |
Superconductor | 1.00 | 2.11 | 2.12 | ** | ** | (9.74) |
Temperature | 1.07 | 1.48 | 1.54 | 1.36 | 2.19 | 1.00 |
Walmart | 1.00 | 1.36 | 1.37 | 11.96 | 12.61 | 1.30 |
Wine | 1.00 | 1.11 | 1.11 | 1.07 | 1.08 | 1.35 |
Mean | 1.18 | 1.73 | 1.74 | 8.08 | 2.09 | (2.88) |
Std | 0.33 | 1.03 | 1.02 | 26.33 | 2.59 | (3.32) |

Each entry is the MSE of that method on that dataset, divided by the

lowest MSE of that row. Entries of at most 1.05 are shown in blue.

**: Exceeded preset maximal wall time.

(xxx): M5 results based on the rescaled predictors.

Table 2 has some unusual entries that require an explanation. The column of
FRIED has an extremely high test error for the Residential dataset. A careful
examination revealed that this was caused by extrapolation. At one node, the
value of the predictor of a test case was outside the range of the training
cases, so lin extrapolated and gave an unreasonable prediction. In the last
column, the M5 results were partially adjusted because its predictions are not
invariant to scaling the predictors. Therefore we ran M5 both with and without
predictor scaling, and present the best result here. Those based on scaling
are shown in parentheses. Unfortunately, the performance of M5 often remained
poor. On the Superconductor data the training time for FRIED and CART exceeded
the preset wall time, despite applying optimized functions (from Scikit-Learn)
for cross validation and pruning.

A closer look at the results reveals that PILOT struggled on the Airfoil data.
A potential explanation is the fact that the Airfoil data has a rather
complicated underlying model, which is often modeled using ensemble methods
such as gradient boosting and random forests (Patri and Patnaik, 2015). The
tree-based methods including PILOT gave relatively poor predictions on the
Riboflavin data, which is high-dimensional with only 78 observations but more
than 7000 predictors, and whose response has a strong linear relationship with
the regressors. This explains why ridge regression and the lasso performed so
well on these data. The tree-based methods suffered from the relatively small
number of cases in the training set, which resulted in a higher test error.
Still, PILOT did better than CART, FRIED and M5 on these data, suggesting that
it is indeed less prone to overfitting and better at capturing linear
relationships.

There are several other datasets on which Ridge and/or Lasso do very well,
indicating that the data have linear patterns. The results show that PILOT
tended to do quite well on those datasets and to outperform CART on them. This
confirms that the linear components of PILOT make it well suited for such
data, whereas we saw that it usually did at least as well as the other tree-
based methods on the other datasets.

###  4.3 Results after transforming predictors

It has been observed that linear regression methods tend to perform better if
the numerical predictors are not too skewed. In order to address potentially
skewed predictors, we reran the empirical experiments after preprocessing all
numerical predictors by the Yeo-Johnson (YJ) transformation (Yeo and Johnson,
2000). The parameter of the YJ transformation was fit using maximum likelihood
on the training set, after which the fitted transformation parameter was
applied to the test set. One way in which such a transformation could help is
by mitigating the extrapolation issue as this is expected to occur mainly on
the long tailed side of a skewed predictor, which becomes shorter by the
transformation.

Table 3: MSE ratios relative to best, after transforming the predictors.

| PILOT | Ridge | LASSO | FRIED | CART | M5
---|---|---|---|---|---|---
Abalone | 0.98 | 1.02 | 1.00 | 1.21 | 1.12 | 1.08
Admission | 1.19 | 1.00 | 1.00 | 1.12 | 1.39 | 1.20
Airfoil | 1.66 | 4.70 | 4.71 | 1.11 | 1.78 | 9.98
Bike | 1.01 | 3.04 | 3.05 | 1.29 | 1.05 | 2.20
Bodyfat | 1.13 | 0.96 | 0.94 | 1.49 | 1.58 | 1.16
Boston | 1.06 | 0.93 | 0.96 | 1.09 | 1.12 | 5.97
California | 1.03 | 1.87 | 1.87 | 1.15 | 1.43 | 24.45
Communities | 1.05 | 1.04 | 1.03 | 1.16 | 1.30 | 1.01
Concrete | 0.95 | 1.24 | 1.24 | 1.23 | 1.39 | 2.48
Diabetes | 1.09 | 1.01 | 1.01 | 1.22 | 1.32 | 1.20
Electricity | 1.00 | 2.75 | 2.75 | 1.29 | 1.87 | 2.71
Energy | 1.34 | 3.84 | 3.86 | 1.01 | 1.29 | 1.11
Ozone | 1.25 | 1.07 | 1.00 | 1.16 | 1.01 | 1.22
Residential | 0.78 | 0.98 | 0.81 | 116.31 | 2.01 | 0.81
Riboflavin | 1.53 | 1.01 | 1.28 | 3.57 | 3.00 | 2.65
Skills | 0.96 | 0.93 | 0.93 | 1.03 | 1.17 | 1.06
Superconductor | 1.05 | 2.04 | 2.05 | *** | *** | 9.74
Temperature | 1.10 | 1.53 | 1.61 | 1.42 | 2.17 | 1.60
Walmart | 1.00 | 1.36 | 1.38 | 11.89 | 11.82 | 1.30
Wine | 1.01 | 1.09 | 1.09 | 1.07 | 1.08 | 1.35
Mean | 1.11 | 1.67 | 1.68 | 7.94 | 2.05 | 3.71
Std | 0.20 | 1.08 | 1.09 | 26.36 | 2.42 | 5.59

Each entry is the MSE of that method on that dataset,

divided by the lowest MSE of the same row in Table 2.

**: Exceeded preset maximal wall time.

Table 3 shows the results after transforming the predictors. The MSE ratios
are relative to the lowest MSE of each row in Table 2, so the entries in Table
3 can be compared directly with those in Table 2. We see that the
transformation often enhanced the performance of PILOT as well as that of some
other methods. On the high dimensional Riboflavin data set the score of PILOT
came closer to Ridge and Lasso, whereas the other decision trees kept
struggling. Also on the other data sets with linear structure PILOT was
comparable with both linear methods. The average and the standard deviation of
its score were reduced by the transformation as well. The overall performance
of PILOT remained the best.

##  5 Conclusion

In this paper we presented a new linear model tree algorithm called PILOT. It
is computationally efficient as it has the same time and space complexity as
CART. PILOT is regularized by a model selection rule to avoid overfitting.
This regularization adds no computational complexity, and no pruning is
required once the tree has been built. This makes PILOT faster than most
existing linear model trees. The prediction can easily be interpreted due to
its tree structure and its additive piecewise linear nature. To guide variable
selection, a measure of feature importance can be defined in the same way as
in CART. PILOT applies two truncation procedures to the predictions, in order
to avoid the extreme extrapolation errors that were sometimes observed with
other linear model trees. An empirical study found that PILOT often
outperformed the CART, M5, and FRIED decision trees on a variety of datasets.
When applied to roughly linear data, PILOT behaved more similarly to high-
dimensional linear methods than other tree-based approaches did, indicating a
better ability to discover linear structures. We proved a theoretical
guarantee for its consistency on a general class of functions. When applied to
data generated by a linear model, the convergence rate of PILOT is polynomial.
To the best of our knowledge, this is the first linear model tree with proven
theoretical guarantees.

We feel that PILOT is particularly well suited for fields that require both
performance and explainability, such as healthcare (Ahmad et al., 2018),
business analytics (Bohanec et al., 2017; Delen et al., 2013) and public
policy (Brennan and Oliver, 2013), where it could support decision making.

A future research direction is to integrate PILOT trees as base learners in
ensemble methods, as it is well known that the accuracy and diversity of base
learners benefit the performance of the ensemble. On the one hand we have seen
that PILOT gave accurate predictions on a number of datasets. On the other
hand, the wide choice of models available at each node allows for greater
diversity of the base learner. For these reasons and the fact that PILOT
requires little computational cost, we expect it to be a promising base
learner for random forests and gradient boosting.

## References

  * Agrawal and Choudhary (2016) Ankit Agrawal and Alok Choudhary.  Perspective: Materials informatics and big data: Realization of the “fourth paradigm” of science in materials science.  _APL Materials_ , 4(5):053208, 2016.
  * Ahmad et al. (2018) Muhammad Aurangzeb Ahmad, Carly Eckert, and Ankur Teredesai.  Interpretable machine learning in healthcare.  In _Proceedings of the 2018 ACM International Conference on Bioinformatics, Computational Biology, and Health Informatics_ , pages 559–560, 2018.
  * Bohanec et al. (2017) Marko Bohanec, Mirjana Kljajić Borštnar, and Marko Robnik-Šikonja.  Explaining machine learning models in sales predictions.  _Expert Systems with Applications_ , 71:416–428, 2017.
  * Boryan et al. (2011) Claire Boryan, Zhengwei Yang, Rick Mueller, and Mike Craig.  Monitoring US agriculture: the US Department of Agriculture, National Agricultural Statistics Service, Cropland Data Layer Program.  _Geocarto International_ , 26(5):341–358, 2011\.
  * Breiman (2001) Leo Breiman.  Random forests.  _Machine learning_ , 45(1):5–32, 2001.
  * Breiman et al. (1984) Leo Breiman, Jerome H Friedman, Richard A Olshen, and Charles J Stone.  _Classification and regression trees_.  Routledge, 1984.
  * Brennan and Oliver (2013) Tim Brennan and William L Oliver.  Emergence of machine learning techniques in criminology: implications of complexity in our data and in research questions.  _Criminology & Public Policy_, 12:551, 2013.
  * Bühlmann (2006) Peter Bühlmann.  Boosting for high-dimensional linear models.  _The Annals of Statistics_ , 34(2):559–583, 2006\.
  * Cattaneo et al. (2022) Matias D Cattaneo, Rajita Chandak, and Jason M Klusowski.  Convergence rates of oblique regression trees for flexible function libraries.  _arXiv preprint arXiv:2210.14429_ , 2022.
  * Chen and Guestrin (2016) Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 785–794, 2016.
  * Delen et al. (2013) Dursun Delen, Cemil Kuzey, and Ali Uyar.  Measuring firm performance using financial ratios: A decision tree approach.  _Expert Systems with Applications_ , 40(10):3970–3983, 2013.
  * Dezeure et al. (2015) Ruben Dezeure, Peter Bühlmann, Lukas Meier, and Nicolai Meinshausen.  High-dimensional inference: Confidence intervals, p-values and R-software hdi.  _Statistical Science_ , 30(4):533–558, 2015.
  * Dobra and Gehrke (2002) Alin Dobra and Johannes Gehrke.  SECRET: A scalable linear regression tree algorithm.  In _Proceedings of the 8th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 481–487, 2002.
  * Dua and Graff (2017) Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Freund et al. (2017) Robert M Freund, Paul Grigas, and Rahul Mazumder.  A new perspective on boosting in linear regression via subgradient optimization and relatives.  _The Annals of Statistics_ , 45(6):2328–2364, 2017.
  * Friedman (1979) Jerome H Friedman.  A tree-structured approach to nonparametric multiple regression.  In Th. Gasser and M. Rosenblatt, editors, _Smoothing Techniques for Curve Estimation_ , pages 5–22. Springer, 1979.
  * Friedman (1991) Jerome H Friedman.  Multivariate adaptive regression splines.  _The Annals of Statistics_ , 19(1):1–67, 1991\.
  * Golbeck et al. (2011) Jennifer Golbeck, Cristina Robles, and Karen Turner.  Predicting personality with social media.  In _CHI’11 Extended Abstracts on Human Factors in Computing Systems_ , pages 253–262. ACM Digital Library, 2011.
  * Györfi et al. (2002) László Györfi, Michael Kohler, Adam Krzyzak, and Harro Walk.  _A distribution-free theory of nonparametric regression_.  Springer Series in Statistics, 2002.
  * Hall et al. (2017) Alastair R Hall, Denise R Osborn, and Nikolaos Sakkas.  The asymptotic behaviour of the residual sum of squares in models with multiple break points.  _Econometric Reviews_ , 36:667–698, 2017.
  * Hastie et al. (2009) Trevor Hastie, Robert Tibshirani, and Jerome H Friedman.  _The Elements of Statistical Learning: Data Mining, Inference, and Prediction_.  Springer, 2009.
  * Hornik et al. (2009) Kurt Hornik, Christian Buchta, and Achim Zeileis.  Open-source machine learning: R meets Weka.  _Computational Statistics_ , 24(2):225–232, 2009\.
  * Josse and Husson (2016) Julie Josse and François Husson.  missMDA: A package for handling missing values in multivariate data analysis.  _Journal of Statistical Software_ , 70(1):1–31, 2016.
  * Ke et al. (2017) Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu.  LightGBM: A highly efficient gradient boosting decision tree.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Klusowski (2021) Jason M Klusowski.  Universal consistency of decision trees in high dimensions.  _arXiv preprint arXiv:2104.13881_ , 2021.
  * Loh (2002) Wei-Yin Loh.  Regression tress with unbiased variable selection and interaction detection.  _Statistica Sinica_ , pages 361–386, 2002.
  * Loh et al. (2007) Wei-Yin Loh, Chien-Wei Chen, and Wei Zheng.  Extrapolation errors in linear model trees.  _ACM Transactions on Knowledge Discovery from Data (TKDD)_ , 1(2):1–17, 2007.
  * Malerba et al. (2004) Donato Malerba, Floriana Esposito, Michelangelo Ceci, and Annalisa Appice.  Top-down induction of model trees with regression and splitting nodes.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 26(5):612–625, 2004.
  * Maurer and Pontil (2009) Andreas Maurer and Massimiliano Pontil.  Empirical Bernstein bounds and sample variance penalization.  _arXiv preprint arXiv:0907.3740_ , 2009.
  * Muthen (1990) Bengt Muthen.  Moments of the censored and truncated bivariate normal distribution.  _British Journal of Mathematical and Statistical Psychology_ , 43(1):131–143, 1990.
  * Noulas et al. (2012) Anastasios Noulas, Salvatore Scellato, Neal Lathia, and Cecilia Mascolo.  Mining user mobility features for next place prediction in location-based services.  In _2012 IEEE 12th International Conference on Data Mining_ , pages 1038–1043. IEEE, 2012.
  * Patri and Patnaik (2015) Ashutosh Patri and Yugesh Patnaik.  Random forest and stochastic gradient tree boosting based approach for the prediction of airfoil self-noise.  _Procedia Computer Science_ , 46:109–121, 2015.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Quercia et al. (2011) Daniele Quercia, Michal Kosinski, David Stillwell, and Jon Crowcroft.  Our twitter profiles, our selves: Predicting personality with twitter.  In _2011 IEEE 3rd International Conference on Privacy, Security, Risk and Trust and 2011 IEEE 3rd International Conference on Social Computing_ , pages 180–185. IEEE, 2011.
  * Quinlan (1992) J R Quinlan.  Learning with continuous classes.  In _5th Australian Joint Conference on Artificial Intelligence_ , volume 92, pages 343–348. World Scientific, 1992.
  * Quinlan (1993) J R Quinlan.  C 4.5: Programs for machine learning.  _The Morgan Kaufmann Series in Machine Learning_ , 1993.
  * Scornet et al. (2015) Erwan Scornet, Gérard Biau, and Jean-Philippe Vert.  Consistency of random forests.  _The Annals of Statistics_ , 43(4):1716–1741, 2015.
  * Shamshirband et al. (2020) Shahabbodin Shamshirband, Sajjad Hashemi, Hana Salimi, Saeed Samadianfard, Esmaeil Asadi, Sadra Shadkani, Katayoun Kargar, Amir Mosavi, Narjes Nabipour, and Kwok-Wing Chau.  Predicting standardized streamflow index for hydrological drought using machine learning models.  _Engineering Applications of Computational Fluid Mechanics_ , 14(1):339–350, 2020.
  * Shi et al. (2019) Yu Shi, Jian Li, and Zhize Li.  Gradient boosting with piece-wise linear regression trees.  In _Proceedings of the 28th International Joint Conference on Artificial Intelligence_ , IJCAI’19, page 3432–3438. AAAI Press, 2019.  ISBN 9780999241141.
  * Stiglic et al. (2012) Gregor Stiglic, Simon Kocbek, Igor Pernek, and Peter Kokol.  Comprehensive decision tree models in bioinformatics.  _PloS One_ , 7(3):e33812, 2012.
  * Tibshirani (1996) Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Tropp et al. (2015) Joel A Tropp et al.  An introduction to matrix concentration inequalities.  _Foundations and Trends® in Machine Learning_ , 8(1-2):1–230, 2015.
  * Vogel et al. (2007) David S Vogel, Ognian Asparouhov, and Tobias Scheffer.  Scalable look-ahead linear regression trees.  In _Proceedings of the 13th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 757–764, 2007.
  * Yeo and Johnson (2000) In-Kwon Yeo and Richard A Johnson.  A new family of power transformations to improve normality or symmetry.  _Biometrika_ , 87(4):954–959, 2000.

APPENDIX

##  Appendix A Preliminary results

In this section we provide notations and preliminary results for the theorems
in Section 3.

###  A.1 Notation

We follow the notations in Sections 2 and 3. The $n$ response values form the
vector $Y$ with mean $\overline{Y}$. The values of the predictors of one case
are combined in $X\in\mathbb{R}^{p}$. The design matrix of all $n$ cases is
denoted as $\boldsymbol{X}\in\mathbb{R}^{n\times p}$. Given some tree node $T$
with $t$ cases, we denote by
$\boldsymbol{X}_{T}:=(X_{T_{1}},\dots,X_{T_{t}})^{\top}\in\mathbb{R}^{t\times
p}$ the restricted data matrix and by $X^{(j)}_{T}$ its $j$-th column. The
variance of the $j$-th column is given by
$\hat{\sigma}_{j,T}^{2}:=\sum_{k=1}^{t}(X_{T_{k}}^{(j)}-\overline{X^{(j)}_{T}})^{2}/t$.
We also let $(\hat{\sigma}_{j,T}^{U})^{2}$ be the classical unbiased variance
estimates with denominator $t-1$. $T$ is omitted in these symbols when $T$ is
the root node or if there is no ambiguity.

We denote by $\mathcal{T}_{k}$ the set of nodes of depth $k$ together with the
leaf nodes of depth less than $k$. The PILOT prediction on these nodes is
denoted as $\hat{f}(\mathcal{T}_{k})\in\mathbb{R}^{n}$, and $\hat{f}^{k}_{T}$
denotes the prediction of the selected model (e.g. plin) on some
$T\in\mathcal{T}_{k}$ , obtained by fitting the residuals
$Y-\hat{f}(\mathcal{T}_{k-1})$ of the previous level. The impurity gain of a
model on a node $T$ with $t$ cases can be written as

$\widehat{\Delta}^{k+1}(T)=||Y-\hat{f}(\mathcal{T}_{k})||_{t}^{2}-P(t_{l})||Y-\hat{f}(\mathcal{T}_{k})-\hat{f}^{k+1}_{T_{l}}||_{t_{l}}^{2}-P(t_{r})||Y-\hat{f}(\mathcal{T}_{k})-\hat{f}^{k+1}_{T_{r}}||_{t_{r}}^{2},$

where $T_{l}$, $T_{r}$ are the left and right child nodes, containing $t_{l}$
and $t_{r}$ cases, and $P(t_{l})=t_{l}/t$, $P(t_{r})=t_{r}/t$.

Since $\hat{f}(\mathcal{T}_{k})$ is additive and we assume that all
$f\in\mathcal{F}$ are too, we can write the functions on the $j$-th predictor
as $\hat{f}_{j}(\mathcal{T}_{k})$ and $f_{j}$ . The total variation of a
function $f$ on $T$ is denoted as $||f||_{TV(T)}$. For nonzero values $A_{n}$
and $B_{n}$ which depend on $n\rightarrow\infty$, we write $A_{n}\precsim
B_{n}$ if and only if $A_{n}/B_{n}\leq\mathcal{O}(1)$, and $\succsim$ is
analogous. We write $A_{n}\asymp B_{n}$ if $A_{n}/B_{n}=\mathcal{O}(1)$. The
difference of the sets $C_{1}$ and $C_{2}$ is denoted as $C_{1}\backslash
C_{2}$.

###  A.2 Representation of the gain and residuals

In this section we show that the impurity gain of lin and pcon can be written
as the square of an inner product. It is also shown that plin can be regarded
as a combination of these two models. These formulas form the starting point
of the proof of the consistency of PILOT and its convergence rate on linear
models.

######  Proposition 8.

Let $\hat{\beta}$ be the least squares estimator of a linear model fitting
$Y\in\mathbb{R}^{p}$ in function of a design matrix
$\boldsymbol{X}\in\mathbb{R}^{n\times p}$. Then we have

$\langle
Y,\boldsymbol{X}\hat{\beta}\rangle_{n}=||\boldsymbol{X}\hat{\beta}||^{2}_{n}\;.$
(9)

###### Proof.

This follows from
$\hat{\beta}=(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1}\boldsymbol{X}^{\top}Y$
and $\langle
Y,\boldsymbol{X}\hat{\beta}\rangle_{n}=Y^{\top}\boldsymbol{X}\hat{\beta}/n$. ∎

######  Lemma 9 (Representation of LIN).

Let $T$ be a node with $t$ cases. For a lin model on the variable $X^{(j)}$ it
holds that:

$\widehat{\Delta}_{lin}^{k+1}=\left\lvert\Bigg{\langle}Y-\hat{f}(\mathcal{T}_{k}),\frac{X^{(j)}_{T}-\overline{X^{(j)}_{T}}}{\hat{\sigma}_{j,T}}\Bigg{\rangle}_{t}\right\rvert^{2}+\left(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})}\right)^{2}.$
(10)

Moreover, if we normalize the second term in the last expression and denote

$\widetilde{f}^{T}_{lin}:=\frac{X^{(j)}-\overline{X^{(j)}_{T}}}{\sqrt{w(t)}\hat{\sigma}_{j,T}}\mathbbm{1}\\{X\in
T\\}$

where $w(t)=t/n$, we have

$\hat{f}^{k+1}_{T}=\langle
Y-\hat{f}(\mathcal{T}_{k}),\widetilde{f}^{T}_{lin}\rangle_{n}\;\widetilde{f}^{T}_{lin}+\overline{Y-\hat{f}(\mathcal{T}_{k})}\;.$

###### Proof.

By Proposition 8 we have

$\displaystyle\widehat{\Delta}_{lin}^{k+1}$
$\displaystyle=||Y-\hat{f}(\mathcal{T}_{k})||_{t}^{2}-||(Y-\hat{f}(\mathcal{T}_{k}))-(\hat{\alpha}+\hat{\beta}X^{(j)}_{T})||_{t}^{2}$
$\displaystyle=2\langle
Y-\hat{f}(\mathcal{T}_{k}),\hat{\alpha}+\hat{\beta}X^{(j)}_{T}\rangle_{t}-||\hat{\alpha}+\hat{\beta}X^{(j)}_{T}||^{2}_{t}$
$\displaystyle=\langle
Y-\hat{f}(\mathcal{T}_{k}),\hat{\alpha}+\hat{\beta}X^{(j)}_{T}\rangle_{t}\;.$

For $\hat{\beta}$ we don’t need to take the mean of the residuals into
account, hence

$\hat{\beta}=\frac{(X^{(j)}_{T}-\overline{X^{(j)}_{T}})^{\top}(Y-\hat{f}(\mathcal{T}_{k}))}{(X^{(j)}_{T}-\overline{X^{(j)}_{T}})^{\top}(X^{(j)}_{T}-\overline{X^{(j)}_{T}})}=\frac{\langle
Y-\hat{f}(\mathcal{T}_{k}),X^{(j)}_{T}-\overline{X^{(j)}_{T}}\rangle_{t}}{\hat{\sigma}_{j,T}^{2}}\;.$

On the other hand
$\hat{\alpha}=\overline{Y-\hat{f}(\mathcal{T}_{k})}-\hat{\beta}\overline{X^{(j)}_{T}}$,
therefore

$\displaystyle\widehat{\Delta}_{lin}^{k+1}$ $\displaystyle=\hat{\beta}\langle
Y-\hat{f}(\mathcal{T}_{k}),X^{(j)}-\overline{X^{(j)}}\rangle_{t}+\left(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})}\right)^{2}$
$\displaystyle=\left\lvert\Bigg{\langle}Y-\hat{f}(\mathcal{T}_{k}),\frac{X^{(j)}-\overline{X^{(j)}}}{\hat{\sigma}_{j,T}}\Bigg{\rangle}_{t}\right\rvert^{2}+\left(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})}\right)^{2}\;.$

The second result follows from the definition of $\hat{\alpha}$, $\hat{\beta}$
and the above formula. ∎

Next we show a similar result for pcon. Note that here our
$\hat{f}(\mathcal{T}_{k})$ is not a constant on each node, so we need to adapt
the result for CART in Klusowski (2021) to our case.

######  Lemma 10 (Representation of PCON).

Let $T$ be a node with $t$ cases and $T_{l}$, $T_{r}$ be its left and right
children. We then have

$\widehat{\Delta}_{pcon}^{k+1}=\left\lvert\Bigg{\langle}Y-\hat{f}(\mathcal{T}_{k}),\frac{\mathbbm{1}\\{X_{T}\in
T_{l}\\}t_{r}-\mathbbm{1}\\{X_{T}\in
T_{r}\\}t_{l}}{\sqrt{t_{r}t_{l}}}\Bigg{\rangle}_{t}\right\rvert^{2}+\left(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})}\right)^{2}.$

Moreover, if we normalize the second term in the inner product and denote

$\widetilde{f}^{T}_{pcon}:=\frac{\mathbbm{1}\\{X_{T}\in
T_{l}\\}P(t_{r})-\mathbbm{1}\\{X_{T}\in
T_{r}\\}P(t_{l})}{\sqrt{w(t)P(t_{r})P(t_{l})}},$

where $w(t)=t/n$, we have

$\hat{f}^{k+1}_{T}=\langle
Y-\hat{f}(\mathcal{T}_{k}),\widetilde{f}^{T}_{pcon}\rangle_{n}\;\widetilde{f}^{T}_{pcon}+\overline{Y-\hat{f}(\mathcal{T}_{k})}\;.$

###### Proof.

By the definition of pcon, we have

$\hat{f}^{k+1}=(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l})\mathbbm{1}(X_{T}\in
T_{l})+(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r})\mathbbm{1}(X_{T}\in
T_{r})$

where $\overline{Y}_{l}$ is the mean of $Y$ in $T_{l}$,
$\overline{\hat{f}(\mathcal{T}_{k})}_{l}$ is the mean of
$\hat{f}(\mathcal{T}_{k})$ in $T_{l}$, and similarly for the cases in the
right child node. In the following, we also denote the constant predictions
for the left and right node as
$\hat{f}^{k+1}_{l}:=\hat{f}^{k+1}\mathbbm{1}(X_{T}\in T_{l})$ and
$\hat{f}^{k+1}_{r}:=\hat{f}^{k+1}\mathbbm{1}(X_{T}\in T_{r})$. Thus we have

$\displaystyle\widehat{\Delta}_{pcon}^{k+1}$
$\displaystyle=\frac{1}{t}\Big{[}2\sum_{X_{T}\in
T_{l}}(Y-\hat{f}(\mathcal{T}_{k}))\hat{f}^{k+1}-\sum_{X_{T}\in
T_{l}}(\hat{f}^{k+1})^{2}$ $\displaystyle+2\sum_{X_{T}\in
T_{r}}(Y-\hat{f}(\mathcal{T}_{k}))\hat{f}^{k+1}-\sum_{X_{T}\in
T_{r}}(\hat{f}^{k+1})^{2}\Big{]}$
$\displaystyle=\frac{1}{t}[2t_{l}(\hat{f}^{k+1}_{l})^{2}-t_{l}(\hat{f}^{k+1}_{l})^{2}+2t_{r}(\hat{f}^{k+1}_{r})^{2}-t_{r}(\hat{f}^{k+1}_{r})^{2}]$
$\displaystyle=\frac{1}{t}[t_{l}\left(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l}\right)^{2}+t_{r}\left(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r}\right)^{2}]$
$\displaystyle=\frac{t_{l}t_{r}}{t^{2}}\Big{[}\left(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l}\right)^{2}+\left(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r}\right)^{2}+\frac{t_{l}}{t_{r}}\left(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l}\right)^{2}+\frac{t_{r}}{t_{l}}\left(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r}\right)^{2}\Big{]}$
$\displaystyle\overset{(i)}{=}\frac{t_{l}t_{r}}{t^{2}}\Big{[}\left(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l}\right)^{2}+\left(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r}\right)^{2}-2(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l})(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r})\Big{]}$
$\displaystyle\quad\quad+\left(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})}\right)^{2}$
$\displaystyle=\frac{t_{l}t_{r}}{t^{2}}\Big{[}(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l})-(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r})\Big{]}^{2}+\left(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})}\right)^{2}$

where (i) follows from the fact that

$\displaystyle\frac{t_{l}}{t_{r}}\left(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l}\right)^{2}$
$\displaystyle+\frac{t_{r}}{t_{l}}\left(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r}\right)^{2}+2(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l})(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r})$
$\displaystyle=\frac{1}{t_{l}t_{r}}\Big{(}\sum_{X_{T}\in
T_{l}}(Y-\hat{f}(\mathcal{T}_{k}))+\sum_{X_{T}\in
T_{r}}(Y-\hat{f}(\mathcal{T}_{k}))\Big{)}^{2}$
$\displaystyle=\frac{t^{2}}{t_{l}t_{r}}\left(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})}\right)^{2}.$

Therefore,

$\displaystyle\widehat{\Delta}_{pcon}^{k+1}$
$\displaystyle=\left\lvert\frac{1}{t}\frac{t_{l}t_{r}(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l})-t_{l}t_{r}(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r})}{\sqrt{t_{l}t_{r}}}\right\rvert^{2}+\left(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})}\right)^{2}$
$\displaystyle=\left\lvert\Bigg{\langle}Y-\hat{f}(\mathcal{T}_{k}),\frac{\mathbbm{1}\\{X_{T}\in
T_{l}\\}t_{r}-\mathbbm{1}\\{X_{T}\in
T_{r}\\}t_{l}}{\sqrt{t_{r}t_{l}}}\Bigg{\rangle}_{t}\right\rvert^{2}+\left(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})}\right)^{2}$

Finally we have

$\displaystyle\quad\langle
Y-\hat{f}(\mathcal{T}_{k}),\widetilde{f}^{T}_{pcon}\rangle_{n}\widetilde{f}^{T}_{pcon}+\overline{Y-\hat{f}(\mathcal{T}_{k})}$
$\displaystyle=\frac{1}{t}[(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l})-(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r})](\mathbbm{1}\\{X_{T}\in
T_{l}\\}t_{r}-\mathbbm{1}\\{X_{T}\in
T_{r}\\}t_{l})+\overline{Y-\hat{f}(\mathcal{T}_{k})}$
$\displaystyle=\frac{t_{l}+t_{r}}{t}[(\overline{Y}_{l}-\overline{\hat{f}(\mathcal{T}_{k})}_{l})\mathbbm{1}\\{X_{T}\in
T_{l}\\}+(\overline{Y}_{r}-\overline{\hat{f}(\mathcal{T}_{k})}_{r})\mathbbm{1}\\{X_{T}\in
T_{r}\\}]$
$\displaystyle\quad-\overline{Y-\hat{f}(\mathcal{T}_{k})}(\mathbbm{1}\\{X_{T}\in
T_{l}\\}+\mathbbm{1}\\{X_{T}\in
T_{r}\\})+\overline{Y-\hat{f}(\mathcal{T}_{k})}$
$\displaystyle=\hat{f}^{k}_{T}$

∎

Note that plin is equivalent to a combination of the two preceding models
(pcon followed by two lin models). As a result, we have a similar expansion
for its prediction.

######  Proposition 11 (Representation of PLIN).

For plin, there exists a function $\tilde{f}^{T}_{\text{plin}}$ depending on
$X$, $T$, $T_{l}$ and $T_{r}$ such that

$\hat{f}^{k+1}_{T}=\langle
Y-\hat{f}(\mathcal{T}_{k}),\widetilde{f}^{T}_{plin}\rangle_{n}\;\widetilde{f}^{T}_{plin}+\overline{Y-\hat{f}(\mathcal{T}_{k})}$

###### Proof.

As plin can be regarded as a pcon followed by two lin fits on the same
predictor, we have

$\displaystyle\hat{f}^{k+1}_{T}$ $\displaystyle=\langle
Y-\hat{f}(\mathcal{T}_{k}),\widetilde{f}^{T}_{pcon}\rangle_{n}\;\widetilde{f}^{T}_{pcon}+\langle
Y-\hat{f}(\mathcal{T}_{k})-\hat{f}_{pcon},\widetilde{f}^{T_{l}}_{lin}\rangle_{n}\;\widetilde{f}^{T_{l}}_{lin}+$
$\displaystyle\quad+\langle
Y-\hat{f}(\mathcal{T}_{k})-\hat{f}_{pcon},\widetilde{f}^{T_{r}}_{lin}\rangle_{n}\;\widetilde{f}^{T_{r}}_{lin}+\overline{Y-\hat{f}(\mathcal{T}_{k})}$
$\displaystyle=\langle
Y-\hat{f}(\mathcal{T}_{k}),\widetilde{f}^{T}_{pcon}\rangle_{n}\;\widetilde{f}^{T}_{pcon}+\langle
Y-\hat{f}(\mathcal{T}_{k}),\widetilde{f}^{T_{l}}_{lin}\rangle_{n}\;\widetilde{f}^{T_{l}}_{lin}+\langle
Y-\hat{f}(\mathcal{T}_{k}),\widetilde{f}^{T_{r}}_{lin}\rangle_{n}\;\widetilde{f}^{T_{r}}_{lin}$
$\displaystyle\quad+\overline{Y-\hat{f}(\mathcal{T}_{k})}$

where the last equation follows from $\langle
1,\widetilde{f}_{lin}\rangle_{n}=0$ and the fact that the prediction of pcon
is constant on two child nodes. Moreover, since $\widetilde{f}^{T_{r}}_{lin}$,
$\widetilde{f}^{T_{l}}_{lin}$, $\widetilde{f}^{T}_{pcon}$ are orthogonal to
each other, we can deduce that

$\displaystyle\hat{f}^{k+1}$ $\displaystyle=\langle
Y-\hat{f}(\mathcal{T}_{k}),A\widetilde{f}^{T}_{pcon}+B\widetilde{f}^{T_{l}}_{lin}+C\widetilde{f}^{T_{r}}_{lin}\rangle_{n}\;(A\widetilde{f}^{T}_{pcon}+B\widetilde{f}^{T_{l}}_{lin}+C\widetilde{f}^{T_{r}}_{lin})$
$\displaystyle\quad+\overline{Y-\hat{f}(\mathcal{T}_{k})}$

for some $A,B$ and $C$ depending on $X$ and the nodes $T$, $T_{l}$ and
$T_{r}$. ∎

##  Appendix B Proof of the universal consistency of PILOT

In this section we show the consistency of PILOT for general additive models.
As discussed in Section 3, we begin with the estimation of the impurity gains
of PILOT. Then we derive a recursion formula which allows us to develop an
oracle inequality and prove the consistency at the end. For all the
intermediate results, we assume the conditions in Theorem 2 holds.

###  B.1 A lower bound for $\Delta_{pcon}$

The key to the proof of the consistency is to relate the impurity gain
$\Delta^{k+1}$ to the training error
$R_{k}:=||Y-\hat{f}(\mathcal{T}_{k})||^{2}_{n}-||Y-f||^{2}_{n}$ so that a
recursion formula on $R_{k}$ can be derived. Recently, Klusowski (2021)
developed such a recursion formula for CART. He showed in Lemma 7.1 that for
CART, $\Delta^{k+1}$ is lower bounded by $R_{k}^{2}$, up to some factors.
However, the proof used the fact that $\langle
Y-\overline{Y}_{t},\overline{Y}_{t}\rangle=0$, and this formula no longer
holds when $\overline{Y}_{t}$ is replaced by $\hat{f}(\mathcal{T}_{k})$. The
reason is that we cannot assume that the predictions in the leaf nodes of a
PILOT tree are given by the mean response in that node, due to earlier nodes.
To overcome this issue, we provide a generalized result for pcon in the
following.

######  Lemma 12.

Assuming $R_{k}(T)>0$ in some node $T$, then the impurity gain of pcon on this
node satisfies:

$\hat{\Delta}^{k+1}_{pcon}(\hat{s},\hat{j},T)\geqslant\frac{R_{k}^{2}(T)}{(||f||_{TV}+6B)^{2}}$

where $\hat{s}$, $\hat{j}$ is the optimal splitting point and prediction.

###### Proof.

Throughout this proof, the computations are related to a single node $T$ so
that we sometimes drop the $T$ in this proof for notational convenience. We
first consider the case where the mean of the previous residuals in the node
is zero. We will show that

$\hat{\Delta}^{k+1}_{pcon}(\hat{s},\hat{j},t)\geqslant\frac{|\langle
Y-\hat{f}(\mathcal{T}_{k}),f-\hat{f}(\mathcal{T}_{k})\rangle_{t}|^{2}}{4||f-\hat{f}(\mathcal{T}_{k})||^{2}_{TV}}\;.$

Note that $\hat{f}(\mathcal{T}_{k})$ is the sum of $p_{n}$ linear functions
defined on $p_{n}$ predictors. Let $\hat{\beta}_{j}$ be the slope of the
linear function on $j$-th predictor. Then we define a probability measure
$\Pi(s,j)$ on the space $\mathbb{R}\times\\{1,\dots,p_{n}\\}$ by the following
Radon-Nikodym derivative:

$\frac{d\Pi(s,j)}{d(s,j)}=\frac{|f^{\prime}_{j}(s)-a_{j}|\sqrt{P(t_{l})P(t_{r})}}{\sum_{j=1}^{p_{n}}\int|f^{\prime}_{j}(s^{\prime})-a_{j}|\sqrt{P(t^{\prime}_{l})P(t^{\prime}_{r})}ds^{\prime}}\;.$

Here, $(t_{l},t_{r})$ and $(t^{\prime}_{l},t^{\prime}_{r})$ are the child
nodes due to the split point $s$ and predictor $j$, and $a_{j}$ and
$f^{\prime}_{j}(x)$ are defined as

$a_{j}=\hat{\beta}_{j},\quad
f^{\prime}_{j}(s)=\begin{dcases}\frac{f_{j}(X_{(i)}^{(j)})-f_{j}(X_{(i-1)}^{(j)})}{X_{(i)}^{(j)}-X_{(i-1)}^{(j)}}&X_{(i-1)}^{(j)}<s<X_{(i)}^{(j)}\\\
0&\text{otherwise}\end{dcases}$

where $X_{(1)}^{(j)},\dots X_{(t)}^{(j)}$ are the ordered data points along
the $j$-th predictor.

The idea is that the optimal impurity gain
$\hat{\Delta}_{pcon}(\hat{s},\hat{j},t)$ is always larger than the average
impurity gain with respect to this probability measure:

$\displaystyle\hat{\Delta}_{pcon}(\hat{s},\hat{j},t)\geqslant$
$\displaystyle\int\hat{\Delta}(s,j,t)d\Pi(s,j)$ $\displaystyle=$
$\displaystyle\int|\langle
Y-\hat{f}(\mathcal{T}_{k}),\sqrt{w(t)}\widetilde{f}^{T}_{pcon}\rangle_{t}|^{2}d\Pi(s,j)$
$\displaystyle\geqslant$ $\displaystyle\left(\int|\langle
Y-\hat{f}(\mathcal{T}_{k}),\hat{f}_{t}\rangle_{t}|d\Pi(s,j)\right)^{2}$

where the last inequality follows from Jensen’s inequality, and the
representation in Lemma 10 is used. We now focus on the term inside the
brackets, for which we obtain an analog of (29) in Klusowski (2021):

$\int|\langle
Y-\hat{f}(\mathcal{T}_{k}),\hat{f}_{t}\rangle_{t}|d\Pi(s,j)\geqslant\frac{|\langle
Y-\hat{f}(\mathcal{T}_{k}),\sum_{j=1}^{p_{n}}\int(f^{\prime}_{j}(s)-a_{j})\mathbbm{1}_{\\{X^{(j)}>s\\}}ds\rangle_{t}|}{\sum_{j=1}^{p}\int|f^{\prime}_{j}(s^{\prime})-a_{j}|\sqrt{P_{t^{\prime}_{l}}P_{t^{\prime}_{r}}}ds^{\prime}}\,.$

For the numerator we have

$\displaystyle\Big{\langle}Y-\hat{f}(\mathcal{T}_{k}),\sum_{j=1}^{p_{n}}\int(f^{\prime}_{j}(s)-a_{j})\mathbbm{1}_{\\{X^{(j)}>s\\}}ds\Big{\rangle}_{t}$
$\displaystyle\overset{(i)}{=}$
$\displaystyle\Big{\langle}Y-\hat{f}(\mathcal{T}_{k}),\sum_{j=1}^{p_{n}}\int(f^{\prime}_{j}(s)-a_{j})\mathbbm{1}_{\\{X^{(j)}>s\\}}ds+\sum_{j=1}^{p_{n}}(f_{j}(X_{(1)}^{(j)})-\hat{f}_{j}(\mathcal{T}_{k})|_{x=X_{(1)}^{(j)}})\Big{\rangle}_{t}$
$\displaystyle=$
$\displaystyle\Big{\langle}Y-\hat{f}(\mathcal{T}_{k}),\sum_{j=1}^{p_{n}}\Big{(}f_{j}(X^{(j)}_{(1)})+\sum_{X^{(j)}_{(2)}\leqslant
X^{(j)}_{(i)}\leqslant
X^{(j)}}\frac{f_{j}(X^{(j)}_{(i)})-f_{j}(X^{(j)}_{(i-1)})}{X^{(j)}_{(i)}-X^{(j)}_{(i-1)}}(X^{(j)}_{(i)}-X^{(j)}_{(i-1)})$
$\displaystyle-\hat{f}_{j}(\mathcal{T}_{k})|_{x=X_{(1)}^{(j)}}-a_{j}(X^{(j)}-X^{(j)}_{(1)})\Big{)}\Big{\rangle}_{t}$
$\displaystyle=$
$\displaystyle\Big{\langle}Y-\hat{f}(\mathcal{T}_{k}),\sum_{j=1}^{p_{n}}(f_{j}-\hat{f}_{j}(\mathcal{T}_{k}))\Big{\rangle}_{t}$
$\displaystyle=$
$\displaystyle\Big{\langle}Y-\hat{f}(\mathcal{T}_{k}),f-\hat{f}(\mathcal{T}_{k})\Big{\rangle}_{t}$

where (i) again used $\langle 1,Y-\hat{f}(\mathcal{T}_{k})\rangle_{t}=0$.

For the denominator restricted to one predictor, we have

$\displaystyle\int|f^{\prime}_{j}(s^{\prime})-a_{j}|\sqrt{P_{t^{\prime}_{l}}P_{t^{\prime}_{r}}}ds^{\prime}$
$\displaystyle=$
$\displaystyle\sum_{i=0}^{N(t)}\int_{N(t)P^{\prime}(t^{\prime}_{l})=i}|f^{\prime}_{j}(s^{\prime})-a_{j}|\sqrt{(i/N(t))(1-i/N(t))}ds^{\prime}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N(t)-1}\int_{X_{(i)}^{(j)}}^{X_{(i+1)}^{(j)}}|f^{\prime}_{j}(s^{\prime})-a_{j}|ds^{\prime}\sqrt{(i/N(t))(1-i/N(t))}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N(t)-1}|f_{j}(X_{(i+1)}^{(j)})-a_{j}X_{(i+1)}^{(j)}-f_{j}(X_{(i)}^{(j)})+a_{j}X_{(i)}^{(j)}|\sqrt{(i/N(t))(1-i/N(t))}$
$\displaystyle\leqslant$ $\displaystyle\text{
}\frac{1}{2}\sum_{i=1}^{N(t)-1}|f_{j}(X_{(i+1)}^{(j)})-a_{j}X_{(i+1)}^{(j)}-f_{j}(X_{(i)}^{(j)})+a_{j}X_{(i)}^{(j)}|$
$\displaystyle=$ $\displaystyle\text{
}\frac{1}{2}\sum_{i=1}^{N(t)-1}|(f_{j}(X_{(i+1)}^{(j)})-\hat{f}_{j}(\mathcal{T}_{k})|_{x=X_{(i+1)}^{(j)}})-(f_{j}(X_{(i)}^{(j)})-\hat{f}_{j}(\mathcal{T}_{k})|_{x=X_{(i)}^{(j)}})|$
$\displaystyle=$ $\displaystyle\text{
}\frac{1}{2}||f_{j}-\hat{f}(\mathcal{T}_{k})_{j}||_{TV}(T)\,.$

If we sum over all $j$ predictors we have

$\displaystyle\sum_{j^{\prime}=1}^{p_{n}}\int|f^{\prime}_{j}(s^{\prime})-a_{j}|\sqrt{P_{t^{\prime}_{L}}P_{t^{\prime}_{R}}}ds^{\prime}$
$\displaystyle\leqslant\frac{1}{2}||f-\hat{f}(\mathcal{T}_{k})||_{TV(T)}\;.$

Putting the bounds on the numerator and denominator together, we obtain:

$\displaystyle\int|\langle
Y-\hat{f}(\mathcal{T}_{k}),\hat{f}_{t}\rangle_{t}|d\Pi(s,j)\geqslant$
$\displaystyle\frac{|\langle
Y-\hat{f}(\mathcal{T}_{k}),\sum_{j=1}^{p_{n}}\int(f^{\prime}_{j}(s)-a_{j})\mathbbm{1}_{\\{X^{(j)}>s\\}}ds\rangle_{t}|}{\sum_{j=1}^{p}\int|f^{\prime}_{j}(s^{\prime})-a_{j}|\sqrt{P_{t^{\prime}_{l}}P_{t^{\prime}_{r}}}ds^{\prime}}$
$\displaystyle\geqslant$
$\displaystyle\frac{\Big{\langle}Y-\hat{f}(\mathcal{T}_{k}),f-\hat{f}(\mathcal{T}_{k})\Big{\rangle}_{t}}{\frac{1}{2}||f-\hat{f}(\mathcal{T}_{k})||_{TV(T)}}\;.$

Now note that

$\displaystyle|\langle
Y-\hat{f}(\mathcal{T}_{k}),f-\hat{f}(\mathcal{T}_{k})\rangle_{t}|=$
$\displaystyle|\langle
Y-\hat{f}(\mathcal{T}_{k}),Y-\hat{f}(\mathcal{T}_{k})\rangle_{t}+\langle
Y-\hat{f}(\mathcal{T}_{k}),f-Y\rangle_{t}|$
$\displaystyle\overset{(i)}{\geqslant}$
$\displaystyle||Y-\hat{f}(\mathcal{T}_{k})||^{2}_{t}-|\langle
Y-\hat{f}(\mathcal{T}_{k}),f-Y\rangle_{t}|$
$\displaystyle\overset{}{\geqslant}$
$\displaystyle||Y-\hat{f}(\mathcal{T}_{k})||^{2}_{t}-||Y-f||_{t}||Y-\hat{f}(\mathcal{T}_{k})||_{t}$
$\displaystyle\overset{}{\geqslant}$
$\displaystyle\frac{1}{2}||Y-\hat{f}(\mathcal{T}_{k})||^{2}_{t}-\frac{1}{2}||Y-f||^{2}_{t}$
$\displaystyle=$ $\displaystyle\frac{1}{2}R_{k}(t)\overset{(ii)}{>}0$

where (i) is the triangle inequality and (ii) follows from the assumption that
$R_{k}(T)$ is strictly positive.

Finally, we deal with $\frac{1}{2}||f-\hat{f}(\mathcal{T}_{k})||_{TV(T)}$ for
which we need an upper bound. We have
$\frac{1}{2}||f-\hat{f}(\mathcal{T}_{k})||_{TV(T)}\,\leqslant\frac{1}{2}||f||_{TV}+\frac{1}{2}||\hat{f}(\mathcal{T}_{k})||_{TV(T)}\,$.
It thus remains to bound the total variation of $\hat{f}(\mathcal{T}_{k})$ on
the node $T$. We now show that the sum of the total variations along all the
predictors is bounded by
$(\max\hat{f}(\mathcal{T}_{k})-\min\hat{f}(\mathcal{T}_{k}))\leqslant 6B$.
Since $\hat{f}(\mathcal{T}_{k})$ is linear, the maximum and minimum values are
always attained on vertices of the node (or cube). We can therefore assume
without loss of generality that the cube is $[0,1]^{p_{n}}$ and
$\hat{f}(\mathcal{T}_{k})=\sum_{j=1}^{p_{n}}\hat{\beta}_{j}X^{(j)}$ s.t. for
$\forall j<j^{\prime}$, $\hat{\beta}_{j}<\hat{\beta}_{j^{\prime}}$ and
$\hat{\beta}_{j}<0\Leftrightarrow j<p_{0}$. Now we start with the vertex
$v_{0}\in[0,1]^{p_{n}}$ such that its $j$-th variable is $1$ if and only if
$j<p_{0}$ (otherwise $0$). Then we move to another vertex $v_{1}$ (along the
edge) which is identical to $v_{0}$ except on the first entry, i.e. its first
entry is 0. In addition, we have
$||\hat{f}_{1}(\mathcal{T}_{k})||_{TV}\,=\hat{f}(\mathcal{T}_{k})|_{X=v_{1}}-\hat{f}(\mathcal{T}_{k})|_{X=v_{0}}\,$.
Similarly, for $\forall j$, we let the vertex $v_{j+1}$ be identical to
$v_{j}$ except for the $j$-th variable (by changing $1$ to $0$ or $0$ to $1$)
and we have
$||\hat{f}_{j}(\mathcal{T}_{k})||_{TV}\,=\hat{f}(\mathcal{T}_{k})|_{X=v_{j}}-\hat{f}(\mathcal{T}_{k})|_{X=v_{j-1}}\,$.
Furthermore, it holds that
$\sum_{j=1}^{p_{n}}||\hat{f}(\mathcal{T}_{k})_{j}||_{TV}\,=\max\hat{f}(\mathcal{T}_{k})-\min\hat{f}(\mathcal{T}_{k})\leqslant
6B$. Therefore, we obtain
$\left(\frac{1}{2}||f-\hat{f}(\mathcal{T}_{k})||_{TV(T)}\right)^{2}\leqslant\frac{1}{4}(||f||_{TV}+6B)^{2}$.

Putting everything together yields

$\hat{\Delta}^{k+1}_{pcon}(\hat{s},\hat{j},t)\geqslant\frac{|\langle
Y-\hat{f}(\mathcal{T}_{k}),f-\hat{f}(\mathcal{T}_{k})\rangle_{t}|^{2}}{\frac{1}{4}||f-\hat{f}(\mathcal{T}_{k})||^{2}_{TV}}\geqslant\frac{R_{k}^{2}(T)}{(||f||_{TV}+6B)^{2}}.$

Finally, we treat the general case where the mean of the residuals in the node
is not necessarily zero (which can happen after fitting a blin model). Note
that
$R_{k}(T)=\\{||Y-\hat{f}(\mathcal{T}_{k})-(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})})||^{2}_{t}-||Y-f||^{2}_{t}\\}+(\overline{Y}-\overline{\hat{f}(\mathcal{T}_{k})})^{2}:=R_{k(1)}+R_{k(2)}$,
where $R_{k(1)}$ is the squared training error after regressing out the mean
and $R_{k(2)}$ is the squared mean. Since the first truncation procedure
ensures that $R_{k}(T)\leqslant 16B^{2}$, we have
$(||f||_{TV}+6B)^{2}\geqslant 36B^{2}\geqslant
2R_{k}(T)=2(R_{k(1)}+R_{k(2)})$, so that

$\frac{R_{k(1)}^{2}}{(||f||_{TV}+6B)^{2}}+R_{k(2)}\geqslant\frac{(R_{k(1)}+R_{k(2)})^{2}}{(||f||_{TV}+6B)^{2}}.$

Thus, by the preceding results and the fact that pcon makes the mean of the
residuals zero, we have

$\Delta_{pcon}^{k+1}(\hat{s},\hat{j},t)\geq\frac{R_{k(1)}^{2}}{(||f||_{TV}+6B)^{2}}+R_{k(2)}\geqslant\frac{R_{k}(T)^{2}}{(||f||_{TV}+6B)^{2}}\,.$

∎

###  B.2 Proof of Proposition 4

As plin generalizes pcon, we always have
$\hat{\Delta}^{k}_{plin}\geqslant\hat{\Delta}^{k}_{pcon}$. For other models,
however, it is not possible to develop a lower bound by similarly constructing
a carefully designed probability measure as that in Lemma 12. This is because
we do not have an indicator function which naturally generates $f$ when
associated with $f^{\prime}$ in the integral (see the estimation of the
numerator). Fortunately, we can rely on the BIC criterion to obtain bounds on
the relative impurity gains between the different models. More precisely,
Proposition 4 shows that if blin or lin are chosen, their gain has to be
comparable (i.e., differ by at most a constant factor) to that of pcon. This
ensures no underfitting occurs. Moreover, if con is chosen at one node, all
the subsequent models on that node would also be con which justifies the use
of the con stopping rule. The proof is given in the following.

###### Proof.

Without loss of generality, we may assume no model fits perfectly. Let us
start with the first claim. By the assumption that $BIC_{1}>BIC_{2}$ we have
for any $v_{1},v_{2}$ that

$\displaystyle BIC_{2}$ $\displaystyle<BIC_{1}$ $\displaystyle
t\log\left(\frac{R_{0}-t\Delta_{2}}{t}\right)+v_{2}\log(t)$
$\displaystyle<t\log\left(\frac{R_{0}-t\Delta_{1}}{t}\right)+v_{1}\log(t)$
$\displaystyle\log\Big{(}1-\frac{t\Delta_{2}}{R_{0}}\Big{)}-\log\Big{(}1-\frac{t\Delta_{1}}{R_{0}}\Big{)}$
$\displaystyle<(v_{1}-v_{2})\frac{\log t}{t}$
$\displaystyle\Big{(}1-\frac{t\Delta_{2}}{R_{0}}\Big{)}\Big{/}\Big{(}1-\frac{t\Delta_{1}}{R_{0}}\Big{)}$
$\displaystyle<\exp\Big{(}(v_{1}-v_{2})\frac{\log t}{t}\Big{)}$
$\displaystyle\frac{t\Delta_{2}}{R_{0}}$
$\displaystyle>1+\exp\Big{(}(v_{1}-v_{2})\frac{\log
t}{t}\Big{)}\Big{(}\frac{t\Delta_{1}}{R_{0}}-1\Big{)}$ (11)

Therefore, for any model $i$ that is selected over con by the BIC criterion,
we have

$\frac{t\Delta_{i}}{R_{0}}>1-\exp\Big{(}(1-v_{i})\frac{\log t}{t}\Big{)}$ (12)

Since $t\geqslant 2$, $\log(t)/t$ is always positive, so is the lower bound in
the preceding equations.

We now proceed with the second claim. If (12) holds for model 1, we have by
(B.2),

$\displaystyle\frac{\Delta_{2}}{\Delta_{1}}$
$\displaystyle>\Big{(}1-\exp\Big{(}(v_{1}-v_{2})\frac{\log
t}{t}\Big{)}\Big{)}\frac{R_{0}}{t\Delta_{1}}+\exp\Big{(}(v_{1}-v_{2})\frac{\log
t}{t}\Big{)}$ $\displaystyle>\Big{(}1-\exp\Big{(}(v_{1}-v_{2})\frac{\log
t}{t}\Big{)}\Big{)}\Big{/}\Big{(}1-\exp\Big{(}(1-v_{1})\frac{\log
t}{t}\Big{)}\Big{)}+\exp\Big{(}(v_{1}-v_{2})\frac{\log t}{t}\Big{)}$
$\displaystyle=\frac{1-\exp((1-v_{2})\log t/t)}{1-\exp((1-v_{1})\log t/t)}$
$\displaystyle:=C^{*}\left(\frac{\log t}{t}\right)>0$

On the other hand, if (12) does not hold for model 1, we still get
${\Delta_{2}}>C^{*}{\Delta_{1}}$ by using (B.2) for model 2 and con, its
inverse inequality for model 1 and con, and the fact that function
$f(x)=(1+ax)/(1+bx)$ is monotonically increasing for $x\in[-1,0]$ and
$1>a>b>0$.

Next we find the minimum of $C^{*}$. We compute the derivative of the
numerator and denominator of $C^{*}$ with respect to $\log t/t$ to get
$(v_{2}-1)\exp((1-v_{2})\log t/t)$ and $(v_{1}-1)\exp((1-v_{1})\log t/t)$.
Therefore, for any $\log t/t:=s\in(0,1/2]$ we have

$\displaystyle C^{*}(s)$
$\displaystyle=\frac{\int_{0}^{s}(v_{2}-1)\exp((1-v_{2})r)dr}{\int_{0}^{s}(v_{1}-1)\exp((1-v_{1})r)dr}$
$\displaystyle\geqslant\frac{(v_{2}-1)\int_{0}^{s}\exp((1-v_{1})r)dr}{(v_{1}-1)\int_{0}^{s}\exp((1-v_{1})r)dr}=\frac{v_{2}-1}{v_{1}-1}$

due to $v_{1}\geqslant v_{2}$. Therefore we conclude that
$C^{*}\geqslant(v_{2}-1)/(v_{1}-1)$ for any $t\geqslant 2$.

By previous lemmas we know that the impurity gain of the models can be divided
into two parts. The first is from regressing out the mean and the second is
from the model assumption, which does not depend on the mean of the response
variable. Therefore, we can let $\Delta_{2}=\Delta_{1}+\Delta_{2}^{\prime}$ in
(B.2) when model 1 is con. Here, $\Delta_{2}^{\prime}$ is the gain after
regressing out the mean. Now, if con is better, we have by the inverse
inequality of (B.2) that

$\displaystyle\frac{t\Delta_{2}^{\prime}}{R_{0}}+\frac{t\Delta_{1}}{R_{0}}$
$\displaystyle<1-\exp\Big{(}(1-v_{2})\frac{\log
t}{t}\Big{)}\Big{)}\Big{(}1-\frac{t\Delta_{1}}{R_{0}}\Big{)}$
$\displaystyle\Longleftrightarrow\frac{t\Delta_{2}^{\prime}}{R_{0}}$
$\displaystyle<\Big{(}1-\exp\Big{(}(1-v_{2})\frac{\log
t}{t}\Big{)}\Big{)}\Big{(}1-\frac{t\Delta_{1}}{R_{0}}\Big{)}$
$\displaystyle<1-\exp\Big{(}(1-v_{2})\frac{\log t}{t}\Big{)},$ (13)

which means that the subsequent node still prefers con, even if its gain is 0.
∎

As an example, for our choice of degrees of freedom,
$\Delta_{lin}\geqslant\Delta_{pcon}/4$ if lin is chosen. Similarly, if blin is
the preferred model, we must have $\Delta_{blin}\geqslant\Delta_{pcon}$ since
their degrees of freedom are the same.

###  B.3 Proof of Theorem 5

The remaining issue is that the RSS in con nodes does not get improved as the
depth increases, since subsequent nodes will select con again. Therefore, we
first construct a recursion formula which excludes the terms corresponding to
the RSS in the con nodes. Then we will show that the training error in these
con nodes vanishes asymptotically, which justifies the con stopping rule.

Recall the definitions of
$C_{k}^{+}=\\{T|T=\textsc{con},T\in\mathcal{T}_{k-1},R_{k-1}(T)>0\\}$, the set
of nodes on which con is fitted before the $k$-th step, and $C_{k}^{-}$ for
those with $R_{k}(T)\leqslant 0$. Further define
$C_{k}^{*}:=C_{k}^{+}\backslash C_{k-1}^{+}$ and
$C_{k}^{\\#}:=C_{k}^{-}\backslash C_{k-1}^{-}$. Finally, let
$A_{k}^{+}:=\\{T|T\in\mathcal{T}_{k},T\neq\textsc{con},R_{k}(T)>0\\}$ and
$A_{k}^{-}:=\\{T|T\in\mathcal{T}_{k},T\neq\textsc{con},R_{k}(T)\leqslant
0\\}$. Note that with these disjoint sets of nodes we now have

$\mathcal{T}_{k-1}=C_{k-1}^{+}\cup C_{k-1}^{-}\cup A_{k-1}^{+}\cup
A_{k-1}^{-}\cup C_{k}^{*}\cup C_{k}^{\\#}\;.$

Next we can calculate the errors in these sets of nodes. Let
$R_{C_{k}^{+}}:=\sum_{T\in C_{k}^{+}}w(T)R(T)$ and define $R_{C_{k}^{-}}$,
$R_{C_{k}^{*}}$, $R_{C_{k}^{\\#}}$, $R_{A_{k-1}^{+}}$, $R_{A_{k-1}^{-}}$ in
similar fashion. This yields

$R_{k-1}=R_{C_{k-1}^{+}}+R_{C_{k-1}^{-}}+R_{A_{k-1}^{+}}+R_{A_{k-1}^{-}}+R_{C_{k}^{\\#}}+R_{C_{k}^{*}}\;.$

Finally, we let $\widetilde{R}_{k}:=R_{k}-R_{C_{k}^{+}}$.

Throughout the proof we assume without loss of generality that the gain of the
con node fitted at depth k is already included into $R_{k}$, therefore
$R_{C^{+}_{k}}$ can be regarded as the remaining error after regressing out
the mean for con nodes before depth $k$.

###### Proof.

Without loss of generality we may assume that $R_{C_{k}^{+}}\leqslant R_{k}$
so that $R_{k}>0$. Otherwise the claim follows immediately. We have

$\displaystyle\widetilde{R}_{k}$ $\displaystyle=R_{k}-R_{C_{k}^{+}}$
$\displaystyle\leqslant R_{k-1}-R_{C_{k}^{+}}-\sum_{T\in
A_{k-1}^{+}}w(T)\Delta(T)$
$\displaystyle\leqslant\widetilde{R}_{k-1}-R_{C_{k}^{*}}-\sum_{T\in
A_{k-1}^{+}}\frac{w(T)}{F}R^{2}_{k-1}(T)$
$\displaystyle\overset{(i)}{\leqslant}\widetilde{R}_{k-1}-R_{C_{k}^{*}}-\frac{1}{F}\Big{(}\sum_{T\in
A_{k-1}^{+}}{w(T)}R_{k-1}(T)\Big{)}^{2}$
$\displaystyle=\widetilde{R}_{k-1}-R_{C_{k}^{*}}-\frac{1}{F}\Big{(}R_{k-1}-R_{C_{k-1}^{+}}-R_{C_{k-1}^{-}}-\sum_{T\in
A_{k-1}^{-}}w(T)R_{k-1}(T)-R_{C_{k}^{\\#}}-R_{C_{k}^{*}}\Big{)}^{2}$
$\displaystyle=:\widetilde{R}_{k-1}-R_{C_{k}^{*}}-\frac{1}{F}\Big{(}R^{*}_{k-1}-R_{C_{k}^{*}}\Big{)}^{2}$

where the second inequality follows from the preceding Lemma and Proposition,
and (i) is Jensen’s inequality. Since $R_{k-1}-R_{C_{k-1}^{+}}$ is positive
and the other three terms in $R^{*}_{k-1}$ are negative, we may replace
$R^{*}_{k-1}$ by $R_{k-1}-R_{C_{k-1}^{+}}=\widetilde{R}_{k-1}>0$ on the right
hand side of the above inequality. Furthermore, we have by definition,

$-R_{C_{k}^{*}}-(\widetilde{R}_{k-1}-R_{C_{k}^{*}})^{2}/F\leqslant\begin{dcases}-\widetilde{R}_{k-1}+F/4<-\widetilde{R}_{k-1}/2&\text{if
}\widetilde{R}_{k-1}>F/2\\\
-\widetilde{R}_{k-1}^{2}/F&\text{otherwise.}\end{dcases}$ (14)

In fact, by assuming that the first regression model is not con, we already
have $\widetilde{R}_{1}\leqslant F/4$ for any initial $R_{0}$ by either
inequality in (14). Therefore, all the subsequent estimations follow the
second inequality and we have for any $k\geqslant 2$,

$\widetilde{R}_{k}\leqslant\widetilde{R}_{k-1}-\frac{1}{F}\widetilde{R}_{k-1}^{2}.$

By a similar induction argument to Lemma 4.1 of Klusowski (2021), we get the
estimation of $R_{K}$ after moving $R_{C^{+}_{K}}$ to the right hand side.

It remains to control the errors in the con nodes. We can assume that the con
models have regressed out the mean in each $C_{k}^{+}$. The first step is to
use Proposition 4 and Lemma 12 to get the following upper bound for any con
node $T$:

$\frac{R(T)^{2}}{FR_{0}}\leqslant\frac{\Delta_{\text{pcon}}}{R_{0}}\leqslant\frac{1}{t}\Big{(}1-\exp\Big{(}(1-v_{pcon})\frac{\log
t}{t}\Big{)}\Big{)}$ (15)

where we used Lemma 12 for the first inequality and (B.2) for the second.
Since we assume the response is almost surely bounded by $\pm B$, we have

$R(T)\leqslant\sqrt{B^{2}F}\sqrt{1-\exp\Big{(}(1-v_{pcon})\frac{\log
t}{t}\Big{)}}\precsim\sqrt{\frac{(v_{pcon}-1)\log t}{t}}$

Therefore the weighted sum of $R(T)$ on $C_{K}^{+}$ is asymptotically

$\sum_{T\in C_{K}^{+}}w(T)R(T)\precsim\frac{1}{n}\sum_{T\in
C_{K}^{+}}\sqrt{t\log t}\;.$

By the Cauchy-Schwarz inequality we obtain

$\sum_{T\in C_{K}^{+}}\sqrt{t\log t}\leqslant\sqrt{2^{K}}\Big{(}{\sum_{T\in
C_{K}^{+}}t\log t}\Big{)}^{1/2}\leqslant\sqrt{2^{K}}\Big{(}{\log n\sum_{T\in
C_{K}^{+}}t}\Big{)}^{1/2}=\sqrt{2^{K}n\log n}\;.$

Therefore if we let $K=\log_{2}(n)/r$ with $r>1$, $R_{C_{K}^{+}}$ is of order
$\mathcal{O}(\sqrt{\log(n)/n^{(r-1)/r}})$, which tends to zero as the number
of cases goes to infinity. ∎

###  B.4 Proof of Theorem 2

Now we can apply Theorem 5 to derive an oracle inequality for our method and
finally prove its consistency. We will use Theorem 11.4 together with Theorem
9.4 and Lemma 13.1 from Györfi et al. (2002), since these results do not
require the estimator to be an empirical risk minimizer.

###### Proof.

Let $f$ be the true underlying function in $\mathcal{F}$, and
$\mathcal{F}_{n}$ be the class of the linear model tree. We first write the
$L^{2}$ error as $||f-\hat{f}(\mathcal{T}_{K})||^{2}=E_{1}+E_{2}$ to apply
Theorem 11.4, where

$E_{1}:=||f-\hat{f}(\mathcal{T}_{K})||^{2}-2(||Y-\hat{f}(\mathcal{T}_{K})||^{2}_{n}-||Y-f||^{2}_{n})-\alpha-\beta$

and

$E_{2}:=2(||Y-\hat{f}(\mathcal{T}_{K})||_{n}^{2}-||Y-f||_{n}^{2})+\alpha+\beta.$

By our assumption, we can define $B_{0}=\max\\{3B,1\\}$ so that the class of
functions and $Y$ are a.s. bounded by $[-B_{0},B_{0}]$, which fulfills the
condition of Theorem 11.4. Thus by Theorem 11.4 with $\epsilon=1/2$ in their
notation, it holds that

$P(\\{\exists\hat{f}({\mathcal{T}_{K}})\in{\mathcal{F}_{n}}\text{ s.t.
}E_{1}>0\\})\leqslant
14\sup_{\mathbf{x}^{n}}\mathcal{N}\Big{(}\frac{\beta}{40B_{0}},\mathcal{F}_{n},L_{1}(\nu_{\mathbf{x}^{n}})\Big{)}\exp\Big{(}-\frac{\alpha
n}{2568B_{0}^{4}}\Big{)}$

where $\alpha,\beta\rightarrow 0$ as $n\rightarrow\infty$.

Theorem 9.4 gives the estimation for the covering numbers of $\mathcal{G}_{n}$
which denotes the class of functions in the leaf nodes. We note that the
condition of the theorem $\beta/40B_{0}<B_{0}/2$ is automatically satisfied
for sufficiently large $n$ if $\beta\rightarrow 0$ is well-defined.

Combining Theorem 9.4 and Lemma 13.1, we get the estimation for the covering
number of $\mathcal{F}_{n}$. Specifically we have,

$\displaystyle\mathcal{N}\Big{(}\frac{\beta}{40B_{0}},\mathcal{F}_{n},L_{1}(\nu_{\boldsymbol{x}^{n}})\Big{)}$
$\displaystyle\leqslant\Gamma_{n}(\Lambda_{n})\Big{[}3\Big{(}4eB_{0}\frac{40B_{0}}{\beta}\log\Big{(}6eB_{0}\frac{40B_{0}}{\beta}\Big{)}\Big{)}^{V(\mathcal{G}_{n})}\Big{]}^{2^{K}}$
$\displaystyle\leqslant\Gamma_{n}(\Lambda_{n})\Big{[}3\Big{(}\frac{160eB_{0}^{2}}{\beta}\log\frac{240eB_{0}^{2}}{\beta}\Big{)}^{V(\mathcal{G}_{n})}\Big{]}^{2^{K}}$
$\displaystyle\leqslant(np)^{2^{K}}\Big{[}3\Big{(}\frac{160eB_{0}^{2}}{\beta}\log\frac{240eB_{0}^{2}}{\beta}\Big{)}^{p+1}\Big{]}^{2^{K}}.$

Here $\Lambda_{n}$ is the set of all possible binary trees on the training set
of $n$ cases, and $\Gamma_{n}(\Lambda_{n})$ is the upper bound for the number
of different partitions on that set that can be induced by binary trees. The
upper bound on $\Gamma_{n}(\Lambda_{n})$ follows from Klusowski (2021) and
Scornet et al. (2015). The exponent $V(\mathcal{G}_{n})$ is the VC dimension
of $\mathcal{G}_{n}$, and since we have multivariate linear predictions
$V(\mathcal{G}_{n})=p+1$.

If we further let $\beta\asymp\frac{B_{0}^{2}}{n}$ and $\alpha\asymp
3B_{0}^{4}\log(np\log n)2^{K+\log(p+1)}/n$ we have

$\mathcal{N}\Big{(}\frac{\beta}{40B_{0}},\mathcal{F}_{n},L_{1}(\nu_{\boldsymbol{x}^{n}})\Big{)}\precsim(np)^{2^{K}}(n\log
n)^{2^{K+\log(p+1)}}$

and

$\displaystyle P(E_{1}>0)$ $\displaystyle\leqslant
P(\\{\exists\hat{f}({\mathcal{T}_{K}})\in{\mathcal{F}_{n}}\text{ s.t.
}E_{1}>0\\})$ $\displaystyle\precsim(np)^{2^{K}}(n\log
n)^{2^{K+\log(p+1)}}\frac{1}{(np\log n)^{3\times 2^{K+\log(p+1)}}}$
$\displaystyle\precsim\frac{1}{(np\log n)^{2^{K+\log(p+1)}}}$
$\displaystyle\leqslant\mathcal{O}\Big{(}\frac{1}{n}\Big{)}.$

Moreover, $E_{1}\leqslant C_{1}$ for some constant $C_{1}$ by the fact that
$||f||_{TV}\;\leqslant A$ and both $Y$ and $f(\mathcal{T}_{k})$ are bounded in
$[-B_{0},B_{0}]$ for any $k$. Thus $\mathbb{E}[E_{1}]\leqslant C_{1}/n$ for
some $C_{1}>0$. By Theorem 5 we also have for any $f\in\mathcal{F}$ that

$\displaystyle\mathbb{E}[E_{2}]$
$\displaystyle\leqslant\frac{2F}{K+3}+\frac{C_{2}\sqrt{2^{K}n\log
n}}{n}+\alpha+\beta,$

and therefore summing up everything we have,

$\displaystyle\mathbb{E}[||f-\hat{f}(\mathcal{T}_{k})||^{2}]\leqslant$
$\displaystyle\frac{2F}{K+3}+\frac{C_{3}\sqrt{2^{K}n\log
n}}{n}+\frac{C_{4}\log(np\log n)2^{K+\log(p+1)}}{n}$ (16)

where $C_{3}$ and $C_{4}$ only depend on $A$ and $B$.

By our assumption, the expected error tends to zero. Moreover, if we pick
$K_{n}=\log(n)/r$ for some $r>1$ and $p_{n}=n^{s}$ such that $0<s<1-1/r$, the
first term in (16) is $\mathcal{O}(1/\log(n))$ and the second and third terms
turn out to be $o(1/\log(n))$. Therefore, the overall convergence rate becomes
$\mathcal{O}(1/\log(n))$. ∎

##  Appendix C Convergence rate for linear model data

###  C.1 Preliminaries and Ideas

In the previous section we derived the consistency of PILOT in case the
underlying is a general additive function. The convergence rate obtained in
that setting is $\mathcal{O}(1/\log(n))$, the same as CART but rather slow. Of
course this is merely an upper bound, and it is not unlikely that we can
obtain much faster convergence if the true underlying function is somehow
well-behaved.

In this section, we consider such a scenario. In particular, we aim to show a
faster convergence rate of PILOT for linear functions. In order to tackle this
problem, we cannot use the same approach as for the general consistency
result. In particular, a counterpart of Lemma 12 for the gain of lin would
still result in a squared error in the left hand side which would restrain the
convergence rate. Therefore, we take a different approach and turn to a recent
tool proposed by Freund et al. (2017) in the context of boosting models. In
our case, we consider a linear regression of $Y_{T}$ on $\boldsymbol{X}_{T}$
in each node $T$ and estimate the gradient of a related quadratic programming
(QP) problem, thereby estimating the impurity gain.

Throughout the proof we use the same notation as that in Freund et al. (2017)
to indicate a subtle difference in the excess error compared to Klusowski
(2021). We first define the notation of the relevant quantities related to the
loss. Let
$L^{*}_{n}(T):=\min_{\hat{\beta}}||Y-\boldsymbol{X}_{T}\hat{\beta}||^{2}_{T}$
be the least squares loss on the node $T$ and
$L^{*}_{n}:=\min_{\hat{\beta}}||Y-\boldsymbol{X}\hat{\beta}||^{2}_{n}$ be the
least squares loss on the full data. We further denote the loss of a $k$ depth
PILOT by $L^{k}_{n}:=||Y-\hat{f}(\mathcal{T}_{k})||^{2}_{n}$. Finally,
$L^{k}_{n}(T):=||Y_{T}-\boldsymbol{X}_{T}\beta_{T}||^{2}_{t}$ is its loss in
node $T$ for some $\beta_{T}$ (we can write the loss like this, since on $T$
the prediction function is linear). When the depth $k$ is not important, we
omit the superscript.

In the following we use the notation $\boldsymbol{\tilde{X}}_{T}$ for the
standardized predictor matrix obtained by dividing each (non-intercept) column
$j$ in $\boldsymbol{X}$ by $\sqrt{n}\hat{\sigma}^{U}_{j}$ where
$\hat{\sigma}^{U}_{j}$ is the usual estimate of the standard deviation of
column $j$. We then consider
$L_{n}(T):=||Y_{T}-\boldsymbol{\tilde{X}}_{T}\tilde{\beta}_{T}||^{2}_{t}$ and
write its gradient as $\nabla L_{n}(T)|_{\tilde{\beta}=\tilde{\beta}_{T}}\,$.

Theorem 2.1 of Freund et al. (2017) used a fixed design matrix so that a
global eigenvalue condition can be imposed on it. As we are fitting linear
models locally, we need a more flexible result. To this end we will impose a
distributional condition that makes use of the following matrix concentration
inequality:

######  Theorem 13 (Theorem 1.6.2 of Tropp et al. (2015)).

Let $\boldsymbol{S}_{1},\dots,\boldsymbol{S}_{n}$ be independent, centered
random matrices with common dimension $d_{1}\times d_{2}$, and assume that
$\mathbb{E}[\boldsymbol{S}_{i}]=0$ and $||\boldsymbol{S}_{i}||\leqslant L$.
Let $\boldsymbol{Z}=\sum_{k=1}^{n}\boldsymbol{S}_{i}$ and define

$v(\boldsymbol{Z})=\max\\{||\mathbb{E}[\boldsymbol{Z}\boldsymbol{Z}^{\top}]||,||\mathbb{E}[\boldsymbol{Z}^{\top}\boldsymbol{Z}]||\\}.$

Then

$P[||Z||\geqslant
t]\leqslant(d_{1}+d_{2})\exp\Big{(}\frac{-t^{2}/2}{v(\boldsymbol{Z})+Lt/3}\Big{)}\;\;\;\text{
for all }t\geqslant 0.$

###  C.2 A probabilistic lower bound for $\Delta_{lin}$

Our first step is to show a probabilistic bound for the impurity gain of
$\Delta_{lin}$ based on the ideas in the preceding section.

######  Lemma 14.

Let $T$ be a node with $n$ cases. We define $L_{n}(T)$ as the training loss of
PILOT and $L^{*}_{n}(T)$ as that of least squares regression. In addition to
conditions 1 and 2, we assume that the residuals of the response on $T$ have
zero mean. Then there exists a constant $C$ such that the lin model satisfies

$P\Bigg{[}\Delta_{\text{lin}}\geqslant\frac{\lambda_{0}(L_{n}(T)-L^{*}_{n}(T))}{4p}\Bigg{]}\geqslant
1-\exp(C_{\lambda_{0},\sigma_{0},p}n).$

###### Proof.

By the definition of the gradient we have,

$\displaystyle\begin{split}\frac{n||\nabla
L_{n}(T)|_{\tilde{\beta}=\tilde{\beta}_{T}}||_{\infty}}{2}&=||\tilde{\boldsymbol{X}}^{\top}r||_{\infty}=\max_{j\in
1\dots p}\\{|r^{\top}\tilde{X}^{(j)}|\\}\\\ &=\max_{j\in 1\dots
p}\frac{|r^{\top}X^{(j)}|}{\sqrt{n}\hat{\sigma}_{j}^{U}}\\\
&\leqslant\max_{j\in 1\dots
p}\frac{|r^{\top}(X^{(j)}-\overline{X}^{(j)})|}{\sqrt{n}\hat{\sigma}_{j}}\\\
&=\sqrt{n\Delta_{lin}}\end{split}$ (17)

where $r$ denotes the residuals $Y-\tilde{\boldsymbol{X}}\tilde{\beta}_{T}$
and the last equality follows from (10) and the assumption that the residuals
of the response on $T$ have zero mean.

By letting $h(\cdot)=L_{n}(\cdot)$,
$Q=\frac{2}{n}\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}}$ in
Proposition A.1 of Freund et al. (2017) and (17), we have (assuming the
difference in the loss is positive),

$\displaystyle
P\Bigg{[}\Delta_{\text{lin}}\geqslant\frac{\lambda_{0}(L_{n}(T)-L^{*}_{n})}{4p}\Bigg{]}$
$\displaystyle\geqslant P\Bigg{[}\frac{n||\nabla
L_{n}|_{\tilde{\beta}=\tilde{\beta}_{T}}||^{2}_{\infty}}{4}\geqslant\frac{\lambda_{0}(L_{n}(T)-L^{*}_{n})}{4p}\Bigg{]}$
$\displaystyle\geqslant P\Bigg{[}\frac{n}{p}||\nabla
L_{n}|_{\tilde{\beta}=\tilde{\beta}_{T}}||_{2}^{2}\geqslant\frac{\lambda_{0}(L_{n}(\tilde{\beta}_{T})-L^{*}_{n})}{p}\Bigg{]}$
$\displaystyle\geqslant
P[\lambda_{min}(\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}})>\lambda_{0}]\,.$

All that is left to show is that

$P[\lambda_{min}(\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}})>\lambda_{0}]\geqslant
1-\exp(-C_{\lambda_{0},\sigma_{0},p}n).$

In order to do so, we compute the minimum eigenvalue of
$\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}}$. We define the vector
$S_{i}\in\mathbb{R}^{p}$ by its coordinates
$S_{i}^{(j)}=X_{i}^{(j)}/\sigma_{j}$ where $\sigma_{j}$ is the true standard
deviation of the $j$-th predictor. By the properties of eigenvalues and
singular values, we can split the estimation into three parts

$\displaystyle\lambda_{min}(\mathbb{E}[SS^{\top}])$
$\displaystyle=\sigma_{min}(\mathbb{E}[SS^{\top}])$
$\displaystyle=\sigma_{min}\Big{[}\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}}+\underbrace{\frac{1}{n}\sum_{i=1}^{n}S_{i}S_{i}^{\top}-\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}}}_{\tilde{A}_{n}}+\underbrace{\frac{1}{n}\sum_{i=1}^{n}\Big{(}\mathbb{E}[SS^{\top}]-S_{i}S_{i}^{\top}}_{\tilde{B}_{n}}\Big{)}]$
$\displaystyle\leqslant\lambda_{min}(\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}})+\sigma_{max}(\tilde{A})+\sigma_{max}(\tilde{B})$

Therefore, we may lower bound the eigenvalues of the matrix
$\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}}$ using the singular
values of 3 other matrices. The left hand side is already bounded by
$2\lambda_{0}$ by our assumption. For $\sigma_{max}(\tilde{A}_{n})$, we note
that the two matrices inside differ only in the estimation of the standard
deviation. If we can show that the difference in the entries are uniformly
bounded by some constant, then we could also upper bound its largest singular
value. To do this, we use Theorem 10 in Maurer and Pontil (2009), which states
that for i.i.d. random variables $X_{1},\dots,X_{n}$ in [0,1] and $\delta>0$,
the classical unbiased variance estimates $(\hat{\sigma}^{U})^{2}$ and the
true variance $\sigma^{2}$ satisfy

$P\Bigg{[}|\sigma-\hat{\sigma}^{U}|>\sqrt{\frac{2\log(1/\delta)}{n-1}}\Bigg{]}<2\delta\;.$

Due to

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}S_{i}S_{i}^{\top}-\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}[S_{i}S_{i}^{\top}-\tilde{X_{i}}\tilde{X_{i}}^{\top}]$
$\displaystyle=\Bigg{[}\frac{\sum_{i=1}^{n}X^{(k)}_{i}X^{(l)}_{i}}{n}\Big{(}\frac{1}{\sigma_{k}\sigma_{l}}-\frac{1}{\hat{\sigma}^{U}_{k}\hat{\sigma}^{U}_{l}}\Big{)}\Bigg{]}_{kl}$

we only need to control the difference in the variance term, as each
$X_{i}^{(k)}$ is bounded by 1. If we put
$\sqrt{2\log(1/\delta)/(n-1)}=\lambda_{0}\sigma^{4}_{0}/2p$ and let $E_{k}$ be
the event that
$|\sigma_{k}-\hat{\sigma}^{U}_{k}|\leqslant\frac{\lambda_{0}\sigma^{4}_{0}}{2p}$,
we have

$P[E_{k}]\geqslant 1-\exp(-C_{\lambda_{0},\sigma_{0},p}n)\;\;\;\text{ and
}\;\;\;\sigma_{k}>\sigma_{0}\text{ on }E_{k}$

due to $\sigma_{0}<1$, $\lambda_{0}<2\lambda_{0}\leqslant
tr(\mathbb{E}[(S-\overline{S})(S-\overline{S})^{\top}])\leqslant p$ and our
assumption that the estimated variance should be lower bounded by
$2\sigma_{0}^{2}$. On $E:=\bigcup_{k}E_{k}$ we have

$\max_{k,l}(|\tilde{A}_{kl}|)\leqslant\max_{k,l}\frac{|\sigma_{k}\sigma_{l}-\hat{\sigma}_{k}^{U}\hat{\sigma}_{l}^{U}|}{\sigma_{k}\sigma_{l}\hat{\sigma}_{k}^{U}\hat{\sigma}_{l}^{U}}\leqslant\max_{k,l}\frac{\sigma_{k}(|\sigma_{l}-\hat{\sigma}_{l}^{U}|)+\hat{\sigma}_{l}^{U}(|\sigma_{k}-\hat{\sigma}_{k}^{U}|)}{2\sigma_{0}^{4}}\leqslant\frac{\lambda_{0}}{2p}$

where the last inequality is due to $\sigma_{k},\sigma_{l}^{U}\leqslant 1$.
Thus $P[\max_{k,l}(|\tilde{A}_{kl}|)\leqslant\lambda_{0}/(2p)]\geqslant
1-p^{2}\exp(-C_{\lambda_{0},\sigma_{0},p}n)$ by the union bound. As
$\max_{||x||_{2}=1}||\tilde{A}x||_{2}\,\leqslant\max_{k,l}(|\tilde{A}_{kl}|)p$
we have

$\displaystyle P[\sigma_{max}(\tilde{A})\leqslant\lambda_{0}/2]$
$\displaystyle=P[\max_{||x||_{2}=1}||\tilde{A}x||_{2}\leqslant\lambda_{0}/2]$
$\displaystyle\geqslant
P[\max_{k,l}(|\tilde{A}_{kl}|)\leqslant\lambda_{0}/2p]\geqslant
1-p^{2}\exp(-C_{\lambda_{0},\sigma_{0},p}n)$

For $\lambda_{max}(\tilde{B}_{n})$ we can apply Theorem 13 to obtain that it
is lower bounded by $\lambda_{0}/2$ with high probability. In fact, on $E$ the
squared $L^{2}$ norm of $S_{i}$ is bounded by $p/\sigma_{0}^{2}$, and
therefore $||S_{i}S_{i}^{\top}||=\lambda_{max}(S_{i}S_{i}^{\top})\leqslant
p/\sigma_{0}^{2}$ so that
$||\mathbb{E}[SS^{\top}]||\leqslant\mathbb{E}[||S_{i}||^{2}]\leqslant
p/\sigma_{0}^{2}$ by Jensen’s inequality. Moreover, for
$\boldsymbol{S}_{i}:=\mathbb{E}[SS^{\top}]-S_{i}S_{i}^{\top}$ we have
$\mathbb{E}[\boldsymbol{S}_{i}]=0$ and

$||\boldsymbol{S}_{i}||\leqslant\frac{1}{n}(||S_{i}S_{i}^{\top}||+||\mathbb{E}[SS^{\top}]||)\leqslant\frac{2p}{n\sigma_{0}^{2}}\,.$

Using the same argument in Section 1.6.3 of Tropp et al. (2015) we find that
$v(\tilde{B}_{n})\leqslant p^{2}/n\sigma_{0}^{4}$ so that on $E$,

$\displaystyle P(||\tilde{B}_{n}||>\frac{\lambda_{0}}{2})$
$\displaystyle\leqslant P\left(||\tilde{B}_{n}||>\frac{\lambda_{0}}{2}\text{
on }E\right)+P(E^{c})$ $\displaystyle\leqslant
2p\exp\Bigg{(}\frac{-\lambda_{0}^{2}\sigma_{0}^{4}n/8}{p^{2}+p\lambda_{0}\sigma^{2}_{0}/3}\Bigg{)}+p^{2}\exp(-C_{\lambda_{0},\sigma_{0},p}n).$

At the end we note that

$2\lambda_{0}\leqslant\lambda_{min}(Cor(X))=\lambda_{min}(\mathbb{E}[(S-\overline{S})(S-\overline{S})^{\top}])\leqslant\lambda_{min}(\mathbb{E}[SS^{\top}]).$

Summing up all the estimations, we know that by a union bound there exists a
constant $C_{\lambda_{0},\sigma_{0},p}>0$ such that
$P[\lambda_{min}(\tilde{\boldsymbol{X}}^{\top}\tilde{\boldsymbol{X}})>\lambda_{0}]\geqslant
1-\exp(-C_{\lambda_{0},\sigma_{0},p}n)$. ∎

###  C.3 Proof of Theorem 6

We can now prove the recursion formula for the expectation of the excess
error.

###### Proof.

We first deal with the case where con is not included. We know that whatever
non-con model is chosen on $T$, its impurity gain is always larger than that
of lin (because the degrees of freedom of lin is the smallest). Moreover, the
first truncation procedure does not deteriorate the gain of the model, since
it only eliminates bad predictions. Therefore when $L^{k}_{n}-L^{*}_{n}$ is
given and positive, by Lemma 14, with probability at least
$\sum_{T\in\mathcal{T}_{k}}\exp(-C_{\lambda_{0},\sigma_{0},p}t)$, it holds
that

$\displaystyle L_{n}^{k+1}-L^{*}_{n}$
$\displaystyle=L_{n}^{k+1}-L_{n}^{k}+L_{n}^{k}-L^{*}_{n}=L_{n}^{k}-L^{*}_{n}-(L_{n}^{k+1}-L_{n}^{k})$
$\displaystyle\leqslant
L_{n}^{k}-L^{*}_{n}-\sum_{T\in\mathcal{T}_{K}}w(T)\Delta_{lin}(T)$
$\displaystyle=L_{n}^{k}-L^{*}_{n}-\sum_{T\in\mathcal{T}_{K}}w(T)(L^{k+1}_{lin,n}(T)-L^{k}_{n}(T))$
$\displaystyle\overset{(i)}{=}L_{n}^{k}-L^{*}_{n}-\sum_{T\in\mathcal{T}_{K}}w(T)\Big{(}(\overline{Y}-\overline{X_{T}\beta_{T}})^{2}+\Big{\langle}Y-X_{T}\beta_{T},\frac{X^{(j)}-\overline{X^{(j)}}}{\sigma_{j}}\Big{\rangle}_{t}^{2}\Big{)}$
$\displaystyle\overset{(ii)}{=}L_{n}^{k}-L^{*}_{n}-\sum_{T\in\mathcal{T}_{K}}w(T)(\bar{r}^{k}_{T})^{2}-\sum_{T\in\mathcal{T}_{K}}w(T)\tilde{\Delta}_{lin}(T)$
$\displaystyle\overset{(iii)}{\leqslant}L_{n}^{k}-L^{*}_{n}-\sum_{T\in\mathcal{T}_{K}}w(T)(\bar{r}^{k}_{T})^{2}$
$\displaystyle\quad-\sum_{T\in\mathcal{T}_{K}}\frac{\lambda_{0}}{4p}w(T)(||r^{k}_{T}-\bar{r}^{k}_{T}||^{2}_{t}-L^{*}_{n}(T))$
$\displaystyle\overset{(iv)}{\leqslant}L_{n}^{k}-L^{*}_{n}-\sum_{T\in\mathcal{T}_{K}}\frac{\lambda_{0}}{4p}w(T)(L^{k}_{n}(T)-L^{*}_{n}(T))$
$\displaystyle\overset{(v)}{\leqslant}L_{n}^{k}-L^{*}_{n}-\frac{\lambda_{0}}{4p}(L^{k}_{n}-L^{*}_{n})$
$\displaystyle=(L_{n}^{k}-L^{*}_{n})\Big{(}1-\frac{\lambda_{0}}{4p}\Big{)}$

where $r^{k}_{T}:=Y-X_{T}\beta_{T}$ are the residuals on $T$. Here, (i)
follows from equation (10). The $\tilde{\Delta}_{lin}(T)$ in (ii) means the
gain of lin after regressing out the mean. The last equation in (iii) follows
from Lemma 14. (iv) is due to
$||r^{k}_{T}-\bar{r}^{k}_{T}||^{2}_{t}+||\bar{r}^{k}_{T}||^{2}_{t}=||r^{k}_{T}||^{2}_{t}$
and the fact that the coefficient $\frac{\lambda_{0}}{4p}$ is less than 1. (v)
follows from the fact that the weighted sum of the optimal loss in each leaf
node $\sum_{T\in\mathcal{T}_{k}}w(T)L^{*}_{n}(T)$ is less than $L^{*}_{n}$ for
each $k$. As argued in the previous lemma, the factor
$1-\frac{\lambda_{0}}{4p}$ is always positive.

It remains to control the ‘bad’ probability, i.e. the case where
$\lambda_{min}(T)<\lambda_{0}$. We have that

$P(\exists T\in\mathcal{T}_{K}\text{ s.t.
}\lambda_{min}(T)<\lambda_{0})\leqslant\sum_{T\in\cup_{k=1}^{K}\mathcal{T}_{k}}\exp(-C_{\lambda_{0},\sigma_{0},p}t)\leqslant
K\sum_{T\in\mathcal{T}_{K}}\exp(-C_{\lambda_{0},\sigma_{0},p}t).$

By condition 1 it holds that

$\displaystyle
K\sum_{T\in\mathcal{T}_{K}}\exp(-C_{\lambda_{0},\sigma_{0},p}t)$
$\displaystyle\leqslant
Kn^{1-\alpha}\exp(-C_{\lambda_{0},\sigma_{0},p}n^{\alpha})$
$\displaystyle=\exp(-C_{\lambda_{0},\sigma_{0},p}n^{\alpha}+\log
K+(1-\alpha)\log n)$
$\displaystyle\precsim\exp(-C^{\prime}_{\lambda_{0},\sigma_{0},p}n^{\alpha})$

as long as $\log K\precsim n^{\alpha}$ which is satisfied by choosing $K$ of
order $\log n$. To sum up everything, if we let $G:=\\{\forall
T\in\mathcal{T}_{K},\lambda_{pmin}(T)>\lambda_{0}\\}$ and assume
$\mathbb{E}[\epsilon^{4}]<\infty$, we have for sufficiently large $n$,

$\displaystyle\mathbb{E}[L^{k}_{n}-L^{*}_{n}]$
$\displaystyle=\mathbb{E}_{X}[\mathbbm{1}_{G}\mathbb{E}_{\epsilon|X}[L^{k}_{n}-L^{*}_{n}]+\mathbbm{1}_{G^{c}}\mathbb{E}_{\epsilon|X}[L^{k}_{n}-L^{*}_{n}]]$
$\displaystyle\leqslant\sqrt{P(X\in
G)}\sqrt{\mathbb{E}_{X}[(\mathbb{E}_{\epsilon|X\in
G}[L^{k}_{n}-L^{*}_{n}])^{2}]}$ $\displaystyle\quad+\sqrt{P(X\in
G^{c})}\sqrt{\mathbb{E}_{X}[(\mathbb{E}_{\epsilon|X\in
G^{c}}[L^{0}_{n}-L^{*}_{n}])^{2}]}$
$\displaystyle\leqslant(1-\exp(-C^{\prime}_{\lambda_{0},\sigma_{0},p}n^{\alpha}))^{1/2}\Big{(}1-\frac{\lambda_{0}}{4p}\Big{)}^{k}\sqrt{\mathbb{E}[(L^{0}_{n}-L^{*}_{n})^{2}]}$
$\displaystyle\quad+\exp(-C^{\prime\prime}_{\lambda_{0},\sigma_{0},p}n^{\alpha})\sqrt{\mathbb{E}[(L^{0}_{n}-L^{*}_{n})^{2}]}$
$\displaystyle\leqslant\Big{(}1-\frac{\lambda_{0}}{4p}\Big{)}^{k}\sqrt{\mathbb{E}[(L^{0}_{n}-L^{*}_{n})^{2}]}+\mathcal{O}(1/n).$

Since we have an exponential decay in the probability, the error can be
bounded by $\mathcal{O}(1/n)$ if we let $K\asymp\log n$.

When con is included, we can derive an analog of (15) by using Lemma 14 and
Proposition 4 to get for all $T$ with con:

$P_{X,\epsilon}\Bigg{[}\lambda_{0}\frac{L_{n}(T)-L^{*}_{n}(T)}{4pR_{0}}<\frac{1}{t}\Big{(}1-\exp\Big{(}(1-v_{lin})\frac{\log
t}{t}\Big{)}\Big{)}\Bigg{]}\geqslant 1-\exp(C_{\lambda_{0},\sigma_{0},p}n).$

Moreover, since the error term has a finite fourth moment and the true
underlying $f$ is linear, we have with probability at least
$1-\exp(C_{\lambda_{0},\sigma_{0},p}n)$ that

$\mathbb{E}_{\epsilon|X}[L_{n}(T)-L^{*}_{n}(T)]\precsim\mathbb{E}_{\epsilon|X}[Y^{2}]{\Big{(}1-\exp\Big{(}(1-v_{lin})\frac{\log
t}{t}\Big{)}\Big{)}}\precsim{\frac{(v_{lin}-1)\log t}{t}}.$

Let us denote by $\mathcal{T}_{K}^{\textsc{con}}$ all the con nodes in the
tree up to depth $K$. Then the expectation of the weighted sum of the errors
on all of these nodes satisfies

$\displaystyle\mathbb{E}_{X,\epsilon}\Big{[}\sum_{T\in\mathcal{T}_{K}^{\textsc{con}}}w(T)(L_{n}(T)-L^{*}_{n}(T))\Big{]}$
$\displaystyle\leqslant\mathbb{E}_{X}\Big{[}\mathbbm{1}_{G^{c}}\sum_{T\in\mathcal{T}_{K}^{\textsc{con}}}\mathbb{E}_{\epsilon|X}[w(T)(L_{n}(T)-L^{*}_{n}(T))]\Big{]}$
$\displaystyle+\mathbb{E}_{X}\Big{[}\mathbbm{1}_{G}\sum_{T\in\mathcal{T}_{K}^{\textsc{con}}}\frac{(v_{lin}-1)\log(t)}{n}\Big{]}$
$\displaystyle\precsim\mathcal{O}\Big{(}\frac{1}{n}\Big{)}+\mathcal{O}\Big{(}\frac{N_{leaves}\log(n)}{n}\Big{)}.$

Noticing that in our previous estimation for non-con nodes, the estimation of
$L^{k+1}_{n}(T)-L^{*}_{n}(T)$ is always positive if
$L^{k}_{n}(T)-L^{*}_{n}(T)>0$ due to $\lambda_{0}\leqslant p$ (see the proof
of Lemma 14), we were able to first ignore the con nodes to get the recursion
formula, and at the end to add back the errors in all con nodes to get the
final results. ∎

###  C.4 Proof of Corollary 7

Finally, by applying the preceding results and the proof of the consistency
theorem, we get the convergence rate of PILOT on linear model data.

###### Proof.

We know that
$\mathbb{E}[L^{*}_{n}-||\boldsymbol{X}\beta||^{2}_{n}]=\mathcal{O}(1/n)$ for
the true parameter $\beta$. So if we choose $K_{n}=log_{\gamma}(n)$, by the
preceding theorem we have
$\mathbb{E}[L^{K_{n}}_{n}-||\boldsymbol{X}\beta||^{2}_{n}]=\mathcal{O}(N_{leaves}\log(n)/n)$.
Finally we can follow the proof of Theorem 2 in which we replace $2^{K}$ by
$N_{leaves}$ and use Theorem 6 to derive a similar oracle inequality as (16):

$\displaystyle\mathbb{E}[||\boldsymbol{X}\beta-\hat{f}(\mathcal{T}_{K})||^{2}]$
$\displaystyle\precsim\mathcal{O}\Big{(}\frac{N_{leaves}\log(n)}{n}\Big{)}+2\mathbb{E}[|L^{K_{n}}_{n}-||\boldsymbol{X}\beta||^{2}_{n}|]$
$\displaystyle\precsim\mathcal{O}\Big{(}\frac{\log(n)}{n^{\alpha}}\Big{)}.$

∎

As a side remark, in the situation where some predictors are correlated with
each other, in practice plin and pcon can lead to even faster convergence than
pure $L^{2}$ boosting.