# Error Reduction from Stacked Regressions††thanks:  The authors would like to
thank Matias Cattaneo, Jianqing Fan, and Giles Hooker for insightful
discussions. Klusowski gratefully acknowledges financial support from the
National Science Foundation through CAREER DMS-2239448, DMS-2054808, and HDR
TRIPODS CCF-1934924. Tan gratefully acknowledges support from NUS Start-up
Grant A-8000448-00-00.

Xin Chen Department of Operations Research and Financial Engineering,
Princeton University Jason Klusowski22footnotemark: 2 Yan Shuo Tan Department
of Statistics and Data Science, National University of Singapore

###### Abstract

Stacking regressions is an ensemble technique that forms linear combinations
of different regression estimators to enhance predictive accuracy. The
conventional approach uses cross-validation data to generate predictions from
the constituent estimators, and least-squares with nonnegativity constraints
to learn the combination weights. In this paper, we learn these weights
analogously by minimizing an estimate of the population risk subject to a
nonnegativity constraint. When the constituent estimators are linear least-
squares projections onto nested subspaces separated by at least three
dimensions, we show that thanks to a shrinkage effect, the resulting stacked
estimator has strictly smaller population risk than best single estimator
among them. Here “best” refers to an estimator that minimizes a model
selection criterion such as AIC or BIC. In other words, in this setting, the
best single estimator is inadmissible. Because the optimization problem can be
reformulated as isotonic regression, the stacked estimator requires the same
order of computation as the best single estimator, making it an attractive
alternative in terms of both performance and implementation.

##  1 Introduction

When performing regression, an analyst rarely knows a priori what the true
model is, and instead starts by deriving a collection of candidate models
$\hat{\mu}_{1},\hat{\mu}_{2},\dots,\hat{\mu}_{M}$. Classically, the next step
in the process is to select the best model based on criteria such as
complexity (e.g., AIC or BIC) or out-of-sample error (e.g., cross-validation).
Such procedures come under the umbrella of model selection, and have been
extensively studied (Hastie et al., 2009).

If predictive performance is the primary consideration, Wolpert (1992)
realized that an alternate approach may be more fruitful. Instead of selecting
the best model, one may use predictions from these estimators as inputs for
another (combined) model, a scheme he called _stacked generalizations_.
Breiman (1996) was able to operationalize Wolpert’s stacking idea by
restricting the combined models to be of the form

$\hat{f}_{\text{stack}}(\mathbf{x})=\sum_{k=1}^{M}\hat{\alpha}_{k}\hat{\mu}_{k}(\mathbf{x}).$
(1)

and learning the weights
$\hat{\alpha}_{1},\hat{\alpha}_{2},\ldots,\hat{\alpha}_{M}$ using cross-
validation. Through extensive experiments on real-world and simulated data, he
consistently observed the following:

  1. (i)

The stacked model (1) has lower mean squared error than the single model
having lowest cross-validation error.

  2. (ii)

The optimal weights sum to approximately one, despite no explicit constraint
to do so.

These initial results have proved to be quite robust and have been reproduced
in many different settings. Indeed, stacking has found widespread applications
in industry (often going by the name _blending_), and has also been a
component of several successful solutions to data science competitions,
including some on the platform Kaggle and, perhaps most famously, the Netflix
Prize (Koren, 2009). It is therefore unfortunate that the beneficial impact of
stacking has yet to be properly understood theoretically. In the statistical
community, variations of stacking have been studied under the name of
aggregation. However, the current state of the art is only able to bound the
excess risk that the stacked model has with respect to the best model selected
by an oracle. This does not quite address the empirical observations.

In our paper, we focus on a special setting in order to better illustrate
these stylized features of stacking. We will focus on the case when the base
estimators comprise a nested sequence of regression models. Such a sequence
arises naturally in the context of sieve estimation (series, spline,
polynomial, wavelet), stepwise regression, and pruning decision trees (e.g.,
CART). Hence, it is fairly general. In such a setting, we can rigorously prove
both of Breiman’s empirical observations for a variant of stacking that
estimates the stacking weights via a regularized mean squared error. Our proof
relies on a novel connection between the optimization program and isotonic
regression. This connection also allows our version of the stacked regression
problem to be solved efficiently.

##  2 Preliminaries

Throughout the paper, we consider a nonparametric regression model. Suppose we
have collected data
$(\mathbf{x}_{1},y_{1}),(\mathbf{x}_{2},y_{2}),\dots,(\mathbf{x}_{n},y_{n})$
satisfying

$\displaystyle y_{i}=f(\mathbf{x}_{i})+\sigma\varepsilon_{i},\quad
i=1,2,\dots,n,$ (2)

where $y_{i}\in\mathbb{R}$ is the label of the $i$-th observation,
$\mathbf{x}_{i}\in\mathbb{R}^{d}$ is the $d$-dimensional feature vector of the
$i$-th observation,
$\varepsilon_{i}\overset{\mathrm{iid}}{\sim}\mathcal{N}(0,1)$ is the $i$-th
unobserved noise variable following the standard normal distribution, and $f$
is the unknown regression function to estimate. We consider the fixed design
problem where the $\mathbf{x}_{i}$’s are deterministic. We also assume that
the noise level $\sigma^{2}$ is known a priori.

For a real-valued function $f$, we define
$\mathbf{f}=(f(\mathbf{x}_{1}),\dots,f(\mathbf{x}_{n}))^{\mathrm{T}}$ to be
the $n\times 1$ vector of $f$ evaluated at the design points
$\mathbf{X}=(\mathbf{x}_{1},\dots,\mathbf{x}_{n})^{\mathrm{T}}\in\mathbb{R}^{n\times
d}$. For real-valued functions $f$ and $g$, we let

$\displaystyle\left\|{\mathbf{f}}\right\|^{2}=\frac{1}{n}\sum_{i=1}^{n}\left(f(\mathbf{x}_{i})\right)^{2},\qquad\left\langle\mathbf{f},\mathbf{g}\right\rangle=\frac{1}{n}\sum_{i=1}^{n}f(\mathbf{x}_{i})g(\mathbf{x}_{i})$

denote the squared empirical norm and empirical inner product, respectively.
The response vector $\mathbf{y}=(y_{1},y_{2},\dots,y_{n})^{\mathrm{T}}$ is
viewed as a relation, defined on the design matrix $\mathbf{X}$, that
associates $\mathbf{x}_{i}$ with $y(\mathbf{x}_{i})=y_{i}$. Thus, we write
$\left\|{\mathbf{y}-\mathbf{f}}\right\|^{2}=\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-f(\mathbf{x}_{i})\right)^{2}$
and
$\left\langle\mathbf{y},\mathbf{f}\right\rangle=\frac{1}{n}\sum_{i=1}^{n}y_{i}f(\mathbf{x}_{i})$.
The $\ell_{0}$ norm of a vector, denoted by $\|\cdot\|_{\ell_{0}}$, is defined
as the number of its nonzero components. We write $(z)_{+}=\max\\{0,z\\}$ for
the positive part of a real number $z$.

###  2.1 Nested Regression Models

In this paper, we investigate the problem of approximating a target variable
$y$ by stacking a sequence of least squares projections onto nested subspaces
corresponding to a different level of approximation. Specifically, we consider
a fixed sequence of orthonormal basis functions
$\\{\psi_{l}\\}_{l\in\mathbb{N}}$, that is, $\|\boldsymbol{\psi}_{l}\|^{2}=1$
and $\langle\boldsymbol{\psi}_{l},\boldsymbol{\psi}_{l^{\prime}}\rangle=0$ for
all $l\neq l^{\prime}$. Let $A_{1}\subseteq A_{2}\subseteq\cdots\subseteq
A_{M}$ be a nested sequence of index sets and define the $k$-th linear
subspace $\mathcal{A}_{k}$ as $\text{span}(\\{\psi_{l}:l\in A_{k}\\})$. The
nested structure implies that each subsequent subspace of dimension
$d_{k}=|A_{k}|$ expands the representation space of the previous one. Let
$\hat{\mu}_{k}(\mathbf{x})$ represent the projection of the response values
$\mathbf{y}=(y_{1},y_{2},\dots,y_{n})^{\mathrm{T}}$ onto the $k$-th linear
subspace $\mathcal{A}_{k}$, given explicitly by

$\hat{\mu}_{k}(\mathbf{x})=\sum_{l\in
A_{k}}\langle\mathbf{y},\boldsymbol{\psi}_{l}\rangle\psi_{l}(\mathbf{x}).$ (3)

We define the mean of $\hat{\mu}_{k}$ by $f_{k}(\mathbf{x})=\sum_{l\in
A_{k}}\langle\mathbf{f},\boldsymbol{\psi}_{l}\rangle\psi_{l}(\mathbf{x})$ and
the empirical risk of $\hat{\mu}_{k}$ by
$R_{k}=\|\mathbf{y}-\boldsymbol{\hat{\mu}}_{k}\|^{2}$. For notational
convenience in some of the upcoming expressions, we define $d_{0}=0$,
$R_{0}=\|\mathbf{y}\|^{2}$, and $\hat{\mu}_{0}(\mathbf{x})\equiv 0$,
corresponding to the _null model_. Without loss of generality, we shall assume
that the models are distinct, so that the set inclusions above are strict. In
this case, because the representation spaces are nested, we have $R_{0}\geq
R_{1}\geq\cdots\geq R_{M}$, with strict inequality holding almost
surely.111This is due to the fact that
$(n/\sigma^{2})(R_{0}-R_{k})\sim\chi^{2}(d_{k},\theta)$ follows a (continuous)
noncentral chi-squared distribution with $d_{k}$ degrees of freedom and
noncentrality parameter
$\theta=n\big{(}\|\mathbf{f}\|^{2}-\|\mathbf{f}-\mathbf{f}_{k}\|^{2}\big{)}/\sigma^{2}$.
We therefore assume throughout the paper that $R_{k}<R_{k-1}$ for all $k$, an
event that holds with probability one.

###  2.2 Examples of Nested Regressions

While there are many nested families of regression models that one could
potentially stack (e.g., sieve estimators), Breiman focused on two canonical
forms that we now describe.

####  2.2.1 Subset Regressions

Using validation data, perform stepwise deletion (elimination) on a linear
model with $d$ variables, that is, start with all $d$ variables and then
discard the one that reduces the empirical risk the least, one by one, until
one is left with a simple linear model. This produces $d$ (nested) linear
models such that the $k$-th model has $d-k+1$ nonzero coefficients. These
nested regressions can then be stacked together using the training data.

####  2.2.2 Decision Trees

A decision tree model can be viewed as a least squares projection onto a set
of piecewise constant functions. The set of piecewise-constant functions is
defined in terms of a (possibly data-dependent) partition of the input space
induced by the splits in the tree. The tree output can therefore be written in
the form (3), where the basis elements $\psi_{l}$ are indexed by the internal
nodes of the tree and depend only on the induced tree partition; see for
example, (Cattaneo et al., 2022, Lemma 2.1).

Nested decision trees serve as a canonical example of nested regressions,
since a refinement of a partition naturally defines an ordering. Each tree
$T_{k}$ is then a subtree of its successor: $T_{1}\preceq
T_{2}\preceq\cdots\preceq T_{M}$, where $T^{\prime}\preceq T$ means that
$T^{\prime}$ can be obtained from $T$ by pruning. As each tree model is
linear, $d_{k}$ is equal to the number of internal nodes in its corresponding
tree $T_{k}$.

To form the nested sequence of trees in practice, Breiman (1996) recommended
growing a large tree $T_{M}$ with $M$ terminal nodes ($M-1$ internal nodes)
and then pruning upwards so that $T_{k}$ is the subtree of $T_{M}$ having the
smallest empirical risk among all subtrees of $T_{M}$ with $k$ terminal nodes
($k-1$ internal nodes). This construction implies that the number of models in
the stack could potentially be quite large, as $M$ can grow polynomially with
$n$.

###  2.3 Accounting for Data-adaptive Subspaces

To preserve the integrity of our theory, we require the sequence of subspaces
to be statistically independent of the responses. On the other hand, subset
regression and decision trees as implemented in practice both select these
subspaces data-adaptively, which in the case of decision trees corresponds to
selecting splits based on a data-dependent criterion, such as impurity
decrease for CART. In order to maintain independence, one may perform data
splitting, using a portion of the data to select the subspaces and the rest of
the data to fit the stacked model. This assumption, sometimes known as the
honesty condition, is often made in the literature for theoretical
tractability (Athey and Imbens, 2016).

##  3 Learning the Stacking Weights

For a sequence of weights
$\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},\ldots,\alpha_{M})$, the empirical
risk of the corresponding stacked model is given by

$R(\boldsymbol{\alpha})=\Bigg{\|}\mathbf{y}-\sum_{k=1}^{M}\alpha_{k}\boldsymbol{\hat{\mu}}_{k}\Bigg{\|}^{2}.$
(4)

###  3.1 Breiman’s Stacking

It is easy to see that directly optimizing (4) leads to a weight vector that
is a point mass on the most complex model $\hat{\mu}_{M}$, and so Breiman
(1996) advocated for optimizing the weights with respect to cross-validation
error. Even so, optimizing over $\boldsymbol{\alpha}$ without constraints
often leads to stacked models that overfit. Recognizing this, Breiman
advocated for some form of regularization. Because the $\hat{\mu}_{k}$ values
are highly correlated—since they aim to predict the same outcome—he initially
experimented with ridge constraints. Ultimately, however, he settled on
nonnegative constraints, i.e., $\alpha_{k}\geq 0$ for $k=1,2,\dots,M$, for
their seemingly best performance. Breiman argued that the reason nonnegative
weights (that also sum to one) work well is because they make the stacked
estimator “interpolating” in the sense that
$\min_{k}\hat{\mu}_{k}(\mathbf{x})\leq\hat{f}_{\text{stack}}(\mathbf{x})\leq\max_{k}\hat{\mu}_{k}(\mathbf{x})$.
On the other hand, he did not impose the constraint that the weights had to
sum to unity as this seemed to have little effect on performance.

While Breiman used cross-validation to learn the weights of combination, we
can quickly see that such an approach is not appropriate for the fixed design
setting because the cross-validation error does not give an unbiased estimate
of the in-sample error. As such, we will instead study a variant of stacking
that optimizes a penalized version of (4) based on model degrees of freedom
and dimension. We believe this to be a small departure from Breiman’s
formulation of stacking, as the asymptotic model selection properties of
cross-validation and penalized empirical risk have been shown to be similar
under certain conditions (Stone, 1977; Shao, 1997). Because of this, our
theory can and indeed does explain some of the empirical findings of Breiman.

###  3.2 Dimension of a Model

We define the dimension of a function $\mu$ by

$\text{dim}(\mu)=\min_{A\subseteq\mathbb{N}}\big{\\{}|A|:\mu\in\text{span}(\\{\psi_{l}:l\in
A\\})\big{\\}},$

that is, the fewest number of basis elements needed to represent the function.
If no such representation exists, then $\text{dim}(\mu)=\infty$, and if
$\mu\equiv 0$, then $\text{dim}(\mu)=0$. Using this definition and setting
$\alpha_{0}=1$, it is easy to see that a stacked model of nested regressions
with fixed weight vector $\boldsymbol{\alpha}$ has dimension
$\max_{k=0,1,\dots,M}\big{\\{}d_{k}:\alpha_{k}\neq 0\big{\\}}.$

###  3.3 Degrees of Freedom of a Stacked Model

In fixed design regression, the degrees of freedom of an estimator $\hat{f}$
is defined as

$\text{df}(\hat{f})=\sum_{i=1}^{n}\frac{\text{cov}(y_{i},\hat{f}(x_{i}))}{\sigma^{2}}.$

It quantifies the optimism bias of the empirical risk with respect to the mean
squared error (Hastie et al., 2009, Section 7.5 & 7.7):

$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}\|^{2}\big{]}=\mathbb{E}\big{[}\|\mathbf{y}-\mathbf{\hat{f}}\|^{2}\big{]}+\frac{2\sigma^{2}\text{df}(\hat{f})}{n}-\sigma^{2}.$
(5)

Although it is generally a population level quantity, in the case of linear
regression, it is equivalent to the number of regressors which is known a
priori; in particular, $\text{df}(\hat{\mu}_{k})=d_{k}$ for each
$k=1,2,\ldots,M$. Moreover, because covariance is linear in each of its
arguments, the degrees of freedom for a stacked model with fixed weight vector
$\boldsymbol{\alpha}$ likewise has a known formula. It is equal to

$\text{df}(\boldsymbol{\alpha})=\sum_{k=1}^{M}\alpha_{k}d_{k}.$ (6)

###  3.4 Stacking via Complexity Penalization

The identity (5) immediately suggests that an unbiased estimator of the
population risk of the stacked model is
$R(\boldsymbol{\alpha})+(2\sigma^{2}/n)\text{df}(\boldsymbol{\alpha})$, which
we may then optimize to learn the stacking weights. However, a major problem
with this estimator is that unbiasedness holds as long as the stacking weights
are nonadaptive. So there is no guarantee that
$R(\hat{\boldsymbol{\alpha}})+(2\sigma^{2}/n)\text{df}(\hat{\boldsymbol{\alpha}})$
is an unbiased estimator of the expected population risk of the stacked model
with adaptive weights $\hat{\boldsymbol{\alpha}}$, where
$\hat{\boldsymbol{\alpha}}$ minimizes
$R(\boldsymbol{\alpha})+(2\sigma^{2}/n)\text{df}(\boldsymbol{\alpha})$ subject
to $\boldsymbol{\alpha}\geq\mathbf{0}$. In fact, it turns out that

$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}=\mathbb{E}\Bigg{[}R(\hat{\boldsymbol{\alpha}})+\frac{2\sigma^{2}}{n}\text{df}(\hat{\boldsymbol{\alpha}})+\frac{4\sigma^{2}}{n}\|\hat{\boldsymbol{\alpha}}\|_{\ell_{0}}-\frac{4\sigma^{2}}{n}\hat{\boldsymbol{\alpha}}^{\mathrm{T}}\hat{\boldsymbol{\alpha}}_{0}-\sigma^{2}\Bigg{]},$
(7)

where
$\hat{\boldsymbol{\alpha}}_{0}(k)=\sum_{j=1}^{k}\mathbf{1}(\hat{\alpha}_{j}\neq
0)$. It therefore becomes clear that an accurate estimator of the population
risk of a stacked model with adaptive weights should involve more than just
the empirical risk $R(\boldsymbol{\alpha})$ and degrees of freedom
$\text{df}(\boldsymbol{\alpha})$, and even be discontinuous in the weights.
Inspired by this observation, we propose learning the weights via solving the
following program:

$\displaystyle\text{minimize}\quad
R(\boldsymbol{\alpha})+\frac{2\tau\sigma^{2}}{n}\text{df}(\boldsymbol{\alpha})+\frac{(\lambda-\tau)^{2}_{+}}{\lambda}\frac{\sigma^{2}}{n}\text{dim}(\boldsymbol{\alpha})$
(8) $\displaystyle\text{subject to}\quad\boldsymbol{\alpha}\geq\mathbf{0},$

where
$\text{dim}(\boldsymbol{\alpha})=\max_{k=0,1,\dots,M}\big{\\{}d_{k}:\alpha_{k}\neq
0\big{\\}}$ with $\alpha_{0}=1$, and $\tau,\lambda>0$ are tuning parameters
chosen by the user. An equivalent way of formulating the program is to
minimize
$R(\boldsymbol{\alpha})+(2\tau\sigma^{2}/n)\text{df}(\boldsymbol{\alpha})$
subject to the constraints $\boldsymbol{\alpha}\geq\mathbf{0}$ and
$\text{dim}(\boldsymbol{\alpha})\leq\text{dim}(\hat{f}_{\text{best}})=d_{\hat{m}}$.

The objective function (8) modulates the empirical risk of the stacked
ensemble in two distinct ways. Firstly, similar to $\ell_{1}$ regularization
in nonnegative Lasso, it penalizes the weights according to the size of the
models, as reflected by the term $\text{df}(\boldsymbol{\alpha})$. Secondly,
similar to $\ell_{0}$ regularization, the objective function also takes into
account the number of basis functions needed to represent the model through
the term $\text{dim}(\boldsymbol{\alpha})$. This promotes model parsimony by
encouraging a sparse selection of models.

While program (8) is nonconvex due to $\text{dim}(\boldsymbol{\alpha})$, we
can leverage the nested structure of the estimators and reduce it to an
isotonic regression problem, which can be solved in $O(M)$ time, i.e., linear
time in the number of models. Note that this is orderwise the same complexity
as finding the best single model (9). This means that we could potentially
stack thousands of models and not worry about computational issues.

It is worth noting that, like conventional stacking, we do not impose any sum
constraint on the weights. Indeed, as we shall see, such a constraint is
essentially superfluous as the unconstrained solution always satisfies
$\sum_{k=1}^{M}\hat{\alpha}_{k}<1$ with near equality in most cases. This fact
also corroborates with Breiman’s experiments, in particular, when stacking
nested decision trees or linear regressions determined by stepwise deletion.
We have also found both empirically and theoretically that enforcing the
equality constraint can lead to inferior performance, as otherwise one is
limiting the potentially beneficial effects of shrinkage.

###  3.5 Best Single Model

In our analysis, we will compare the performance of the stacked model (8) to
that of the _best single model_ , defined as
$\hat{f}_{\text{best}}(\mathbf{x})=\hat{\mu}_{\hat{m}}(\mathbf{x})$, where

$\hat{m}\in\operatorname*{argmin}_{k=0,1,\dots,M}\;R_{k}+\lambda\frac{\sigma^{2}d_{k}}{n},$
(9)

where $\lambda>0$ is a tuning parameter whose value we will take to be the
same in (8) and (9). If ties exist in (9), we choose the smallest $\hat{m}$.
Note that we have also included the null model (i.e., $\hat{\mu}_{0}\equiv 0$)
as a candidate model, which will simplify some of the forthcoming theory.

The criterion $R_{k}+\lambda\sigma^{2}d_{k}/n$ holds significance in various
model selection contexts. Specifically, for $\lambda=2$, it matches Mallows’s
$C_{p}$. In the present Gaussian linear regression setting, if $\lambda=2$, it
is equivalent to the Akaike Information Criterion (AIC) and Stein’s Unbiased
Risk Estimate (SURE) (Stein, 1981), and if $\lambda=\log(n)$, it corresponds
to the Bayesian Information Criterion (BIC). These model selection criteria
operate on a single model at a time, seeking the one that optimally balances
goodness of fit and complexity. In contrast, stacking leverages the diverse
strengths of multiple models through an optimal linear combination.

###  3.6 Stacking Versus Pruning Decision Trees

When the nested subspaces arise from a sequence of decision trees obtained via
Breiman’s approach (see Section 2.2.2), model selection via (9) corresponds
exactly to cost-complexity pruning (Hastie et al., 2009, Section 9.2.2).
Meanwhile, the stacked model $\hat{f}_{\text{stack}}(\mathbf{x})$ also takes
the form of a decision tree with the same structure as the largest tree in the
stack. It consists of $\text{dim}(\hat{f}_{\text{stack}})$ internal nodes and
outputs
$\hat{f}_{\text{stack}}(\mathbf{x})=\sum_{k=1}^{M}\hat{\alpha}_{k}\overline{y}_{k}$,
where $\overline{y}_{k}$ represents the output of tree $T_{k}$ at
$\mathbf{x}$, which is the sample mean of observations in the cell containing
$\mathbf{x}$. Since for $k=1,2,\ldots,M$, the leaf of $T_{k}$ containing
$\mathbf{x}$ is an ancestor of that containing $\mathbf{x}$ in $T_{M}$, this
has the effect of shrinking the predictions over each leaf in $T_{M}$ to its
ancestors’ values, similar to the method of Agarwal et al. (2022).

##  4 Main Results

Our main result is that stacked model with weights from (8) _strictly_
outperforms the best single model (9) provided the dimensions of the
individual models differ by a constant. Put another way, in the parlance of
statistical decision theory, the best single estimator $\hat{f}_{\text{best}}$
(and by implication, the estimator selected by AIC and BIC) is _inadmissible_.
As far as we are aware, this is the first result showing superior performance
of the stacked model compared to the best single model in any context. We
emphasize again that both $\hat{f}_{\text{best}}$ and $\hat{f}_{\text{stack}}$
have the same computational complexity, that is, $O(M)$.

######  Theorem 4.1.

Suppose $0<\tau<2$ and $d_{k}\geq d_{k-1}+4/(2-\tau)$ for all $k$. The
population risk of the stacked model with weights from (8) is strictly less
than the population risk of the best single model (9); that is,

$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}<\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}.$

######  Remark 1.

Since the $d_{k}$ are integers, the hypothesis of Theorem 4.1 becomes
$d_{k}\geq d_{k-1}+3$ when $0<\tau\leq 2/3$. Thus, a sufficient condition for
the stacked model to be provably superior to the best single model is that the
constituent models differ by at least three dimensions. We presently do not
know if the condition is also necessary. If we knew the theoretically optimal
$\tau$, which involves unknown population level quantities, then the condition
can be improved to $d_{k}\geq d_{k-1}+2$.

The proof hinges on a few novel ingredients that we briefly introduce here,
and will explain further in the next section. First, the nestedness of the
base estimators and the nonnegative weight constraints allow us to recast the
optimization program (8) as isotonic regression (see Lemma 5.1). We then build
on known properties of solutions to isotonic regression to obtain an explicit
representation of $\hat{f}_{\text{stack}}$ and compare it to that of
$\hat{f}_{\text{best}}$ (see Theorem 5.2). Finally, to compare the population
risks of the two representations, we apply an extension of Stein’s Lemma for
discontinuous functions (Tibshirani, 2015), which allows us to decompose the
model degrees of freedom into two portions, one of which is attributable to
the model selection mechanism (Tibshirani called this contribution the _search
degrees of freedom_).

This proof strategy in fact yields the following lower bound for the risk gap
$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}-\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}$:

$\frac{\sigma^{2}\tau(2-\tau)}{n}\mathbb{E}\Bigg{[}\min_{1\leq k\leq
M}\frac{(d_{k}-4k/(2-\tau))^{2}}{(n/\sigma^{2})(R_{0}-R_{k})}\Bigg{]}+2\min\Big{\\{}1,\frac{\tau}{\lambda}\Big{\\}}\frac{\sigma^{2}}{n}\text{sdf}(\hat{f}_{\text{best}}).$
(10)

Here, the second term is in terms of the search degrees of freedom of
$\hat{f}_{\text{best}}$. Meanwhile, the first term is reminiscent of the
improvement that would arise from applying James-Stein shrinkage (James and
Stein, 1961) to an individual estimator, viz.,

$\hat{\mu}_{k}^{\,\text{JS}}(\mathbf{x})=\Bigg{(}1-\frac{d_{k}-2}{(n/\sigma^{2})(R_{0}-R_{k})}\Bigg{)}\;\hat{\mu}_{k}(\mathbf{x}),$

with risk gap
$\mathbb{E}\big{[}\|\mathbf{f}-\boldsymbol{\hat{\mu}}_{k}\|^{2}\big{]}-\mathbb{E}\big{[}\|\mathbf{f}-\boldsymbol{\hat{\mu}}_{k}^{\,\text{JS}}\|^{2}\big{]}$
given by

$\frac{\sigma^{2}}{n}\mathbb{E}\Bigg{[}\frac{(d_{k}-2)^{2}}{(n/\sigma^{2})(R_{0}-R_{k})}\Bigg{]}.$
(11)

The similarity between (10) and (11) is not a coincidence. Our proof strategy
shows that one may think of stacking, in the nested regressions setting, as
essentially performing shrinkage to an adaptively selected subset of models.
We will discuss this connection further in the next two sections. As with
usual James-Stein shrinkage (11), the risk gap (10) tends to be larger when
the signal-to-noise ratio $\|\mathbf{f}\|/\sigma$ or sample size $n$ is
small—agreeing with the adage that ensemble methods tend to work particularly
well in these settings. This is suggested analytically from the fact that, for
each $k$, the quantity $(n/\sigma^{2})(R_{0}-R_{k})\sim\chi^{2}(d_{k},\theta)$
is stochastically increasing in
$\theta=n\big{(}\|\mathbf{f}\|^{2}-\|\mathbf{f}-\mathbf{f}_{k}\|^{2}\big{)}/\sigma^{2}$.

######  Remark 2.

We conjecture that the risk gap is uniformly bounded away from zero as the
number of models increases, i.e.,

$\liminf_{M\rightarrow\infty}\Big{[}\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}-\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}\Big{]}>0,$

which would imply that stacking offers a consistent advantage in predictive
performance over the best single model.

Up until now, we have focused on the virtues of stacking in a relative sense;
that is, how well is does relative to another estimator. Here we show how well
it does relative to an oracle, i.e., the stacked model with oracle weights.

######  Theorem 4.2.

The risk of the stacked model with weights from (8) when $\tau=\lambda=1$
satisfies the following oracle inequality:

$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}\leq\mathbb{E}\Bigg{[}\Bigg{\|}\mathbf{f}-\sum_{k=1}^{M}\alpha^{\star}_{k}\boldsymbol{\hat{\mu}}_{k}\Bigg{\|}^{2}\Bigg{]}+\frac{4\sigma^{2}\mathbb{E}[\|\hat{\boldsymbol{\alpha}}\|_{\ell_{0}}]}{n},$

where
$\boldsymbol{\alpha}^{\star}=(\alpha^{\star}_{1},\alpha^{\star}_{2},\dots,\alpha^{\star}_{M})^{\mathrm{T}}$
is the minimizer of the population risk.

Theorem 4.2 says that the expected number of stacked models, or
$\mathbb{E}[\|\hat{\boldsymbol{\alpha}}\|_{\ell_{0}}]$, controls the error
from the oracle stacked model. Breiman observed that this quantity was
surprisingly small. For example, in his experiments with stacked subset
regressions ($M=40$), he found
$\mathbb{E}[\|\hat{\boldsymbol{\alpha}}\|_{\ell_{0}}]\approx 3.1$. The size of
$\|\hat{\boldsymbol{\alpha}}\|_{\ell_{0}}$ is equal to the number of unique
elements of the isotonic sequence $\hat{\boldsymbol{\gamma}}$ that are less
than one, defined in (13) below. Finding useful bounds (other than the number
of models $M$) appears to be nontrivial. While previous shape-constrained
estimation literature has shown that similar isotonic minimax sequences can
have, on average, $o(M)$ distinct elements (Meyer and Woodroofe, 2000), we
have yet to establish this for $\hat{\boldsymbol{\gamma}}$.

##  5 Reduction to Isotonic Regression

Due to the nested structure of the models and the nonnegative weight
constraint, the problem of determining the stacking weights from problem (8)
can be recast as isotonic regression. To see this, let $\Delta
d_{k}=d_{k}-d_{k-1}$ and $\Delta R_{k}=R_{k-1}-R_{k}$. Recall that
$R_{k}<R_{k-1}$ for all $k$, implying that $\Delta R_{k}>0$ for all $k$.
Making the change of variables $\alpha_{k}=\beta_{k+1}-\beta_{k}$, for
$k=1,2,\dots,M$, with $\beta_{M+1}=1$, the nonnegativity constraint on the
weights $\alpha_{k}$ becomes an isotonic constraint on the $\beta_{k}$. Some
further algebraic manipulation gives the following equivalence of optimization
programs.

######  Lemma 5.1.

Program (8) is equivalent to:

$\displaystyle\text{minimize}\quad\sum_{k=1}^{M}w_{k}(z_{k}-\beta_{k})^{2}+\xi\sum_{k=1}^{M}\Delta
d_{k}\mathbf{1}(\beta_{k}\neq 1)$ (12) $\displaystyle\text{subject
to}\quad\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{M}\leq 1,$

where $w_{k}=\Delta R_{k}>0$, $z_{k}=(\tau\sigma^{2}/n)(\Delta d_{k}/\Delta
R_{k})>0$, and $\xi=(\sigma^{2}/n)((\lambda-\tau)^{2}_{+}/\lambda)$.

We shall henceforth refer to the stacking weights $\hat{\boldsymbol{\alpha}}$
as _the_ solution to program (8), since it turns out to be almost surely
unique. Crucially, using known properties of isotonic regression, the solution
admits a closed form expression in terms of the following nondecreasing
minimax sequence (13), which allows us to connect $\hat{f}_{\text{stack}}$ and
$\hat{f}_{\text{best}}$:

$\hat{\gamma}_{k}=\frac{\sigma^{2}}{n}\min_{k\leq i\leq M}\max_{0\leq
j<k}\frac{d_{i}-d_{j}}{R_{j}-R_{i}},\quad k=1,2,\dots,M.$ (13)

The vector of this sequence is denoted as
$\hat{\boldsymbol{\gamma}}=(\hat{\gamma}_{1},\hat{\gamma}_{2},\dots,\hat{\gamma}_{M})^{\mathrm{T}}$.
In what follows, it will be convenient to define $\hat{\gamma}_{0}=0$ and
$\hat{\gamma}_{M+1}=\infty$.

######  Theorem 5.2.

The (almost surely) unique solution to program (8) is

$\hat{\alpha}_{k}=(1-\tau\hat{\gamma}_{k})\mathbf{1}(\hat{\gamma}_{k}<\gamma)-(1-\tau\hat{\gamma}_{k+1})\mathbf{1}(\hat{\gamma}_{k+1}<\gamma),\quad
k=1,2,\dots,M,$ (14)

where $\gamma=\min\\{1/\tau,1/\lambda\\}$, and, consequently,
$\sum_{k=1}^{M}\hat{\alpha}_{k}=(1-\tau\hat{\gamma}_{1})\mathbf{1}(\hat{\gamma}_{1}<\gamma)<1$.
Furthermore, the stacked model can be written as

$\hat{f}_{\text{stack}}(\mathbf{x})=\sum_{k=1}^{M}(\hat{\mu}_{k}(\mathbf{x})-\hat{\mu}_{k-1}(\mathbf{x}))(1-\tau\hat{\gamma}_{k})\mathbf{1}(\hat{\gamma}_{k}<\gamma),$
(15)

and the best single model can be written as

$\hat{f}_{\text{best}}(\mathbf{x})=\sum_{k=1}^{M}(\hat{\mu}_{k}(\mathbf{x})-\hat{\mu}_{k-1}(\mathbf{x}))\mathbf{1}(\hat{\gamma}_{k}<1/\lambda).$
(16)

According to (16), the best single estimator can be written as a telescoping
sum of predictive differences
$\hat{\mu}_{k}(\mathbf{x})-\hat{\mu}_{k-1}(\mathbf{x})$ across successive
submodels, up to the selected model $\hat{\mu}_{\hat{m}}$. In contrast,
according to (15), stacking additionally shrinks these predictive differences
towards zero, and because $\\{\hat{\gamma}_{k}\\}$ is an increasing sequence,
larger models are shrunk more than smaller models. The adaptive shrinkage
factor $1-\tau\hat{\gamma}_{k}$ in (15) is the driving force behind its
superior performance. So in summary, in the setting of nested regressions,
_stacking performs both model selection and shrinkage simultaneously._

###  5.1 Implementation

Because we know the general form of the solution (14) in terms of
$\hat{\boldsymbol{\gamma}}$, a quantity free of $\lambda$ and $\tau$, we
simply need solve for $\hat{\boldsymbol{\gamma}}$ using the related program

$\displaystyle\text{minimize}\quad\sum_{k=1}^{M}w_{k}(z_{k}-\beta_{k})^{2}$
(17) $\displaystyle\text{subject
to}\quad\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{M},$

which we recognize as a weighted isotonic regression problem. The solution
admits a closed form expression
$\hat{\boldsymbol{\beta}}=\hat{\boldsymbol{\gamma}}$, the minimax sequence
(13), from which the general solution (14) can be directly obtained. Program
(17) can be solved in $O(M)$ time using the Pooled Adjacent Violators
Algorithm (PAVA) (Barlow et al., 1972), thereby permitting $M$ to be in the
thousands without incurring worrisome computational burden. The isoreg()
function from base R will implement program (17).

##  6 Model Selection and Shrinkage in Stacking

According to Theorem 5.2, the version of stacking we analyzed performs both
model selection and shrinkage simultaneously. As discussed after Theorem 4.1
and as will be proved in the appendix, the improvement in prediction
performance over the best single model derives almost entirely from shrinking
the weights in the telescoping sum (16). Since the shrinkage factors are
optimized, they naturally resemble James-Stein shrinkage factors, leading to a
similar formula for the risk gap (10). In this section, we compare shrinkage
from stacking with prior work on shrinkage and discuss the generality of this
interpretation of stacking.

Agarwal et al. (2022) introduced a hierarchical shrinkage procedure to
regularize decision tree models. When viewing the tree model as a linear
regression onto orthogonal basis elements $\psi_{l}$ indexed by the internal
nodes of the tree, they showed that their procedure was equivalent to
performing ridge regression instead of linear regression on this basis. Since
the amount of shrinkage in ridge regression is controlled by a single
parameter, theirs is therefore a more constrained form of shrinkage compared
to what arises from stacking nested subtrees. The latter optimizes the amount
of shrinkage over $M$ different parameters, each controlling the amount of
shrinkage over a separate block of the regression coefficients. On the other
hand, hierarchical shrinkage allows the basis elements $\psi_{l}$ to be
unnormalized, with $\|\psi_{l}\|^{2}=N_{l}/n$, where $N_{l}$ is the number of
samples contained in node $l$. The ridge regression shrinkage factor for the
coefficient of $\psi_{l}$ is then $\frac{N_{l}}{N_{l}+\lambda}$, which means
that splits on nodes with fewer samples are penalized more strongly. A
theoretical or empirical performance comparison of the two forms of shrinkage
is left for future work.

Agarwal et al. (2022) also showed empirically that the best performing tree
model after shrinkage is usually much deeper than the best performing tree
model prior to shrinkage. In other words, by allowing a more nuanced bias-
variance tradeoff, shrinkage allows for a “larger model” and makes use of
features whose inclusion would not be justified otherwise. Breiman (1996)’s
empirical findings for stacked subtrees present a similar story. In one
experiment he found that, among a collection of $M=50$ subtrees, the stacked
tree had $33$ internal nodes (the _dimension of the model_ in our
terminology), whereas the (underperforming) best single tree had only $22$.
Breiman was surprised by this, and explained it as the prediction in the
stacked tree “gathering strength” from the data among ancestor terminal nodes
in the smaller stacked subtrees.

Our work formalizes Breiman’s notion of “gathering strength” as shrinkage, and
rigorously proves its benefits. On the other hand, we are not yet able to show
theoretically how the prediction performance of stacking benefits from
allowing for a larger model. The optimization program (8) was explicitly
designed so that the additive representations (15) and (16) for
$\hat{f}_{\text{stack}}$ and $\hat{f}_{\text{best}}$ truncate at the same
level. This property is crucial in our proof of Theorem 4.1 because it allows
us to compare the search degrees of freedom of the two estimators,
$\text{sdf}(\hat{f}_{\text{stack}})$ and $\text{sdf}(\hat{f}_{\text{best}})$,
which would otherwise be challenging. Recall that these quantities arise
because we use an extension of Stein’s Lemma for discontinuous functions
(Tibshirani, 2015). We need this to cope with the fact that both
$\hat{f}_{\text{stack}}$ and $\hat{f}_{\text{best}}$ are discontinuous with
respect to the $y$-data—they perform hard thresholding on a data-dependent
quantity, arising from the model selection mechanism. In general, bounding the
search degrees of freedom (e.g., from best subset selection) is known to be
quite challenging (Mikkelsen and Hansen, 2018; Tibshirani and Rosset, 2019).
Truncating the additive representations for $\hat{f}_{\text{stack}}$ and
$\hat{f}_{\text{best}}$ at the same level places their discontinuities (with
respect to the $y$-data) at exactly the same locations, which causes their
search degrees of freedom to be positive and proportional to each other, that
is,
$\text{sdf}(\hat{f}_{\text{stack}})=(1-\tau/\lambda)_{+}\text{sdf}(\hat{f}_{\text{best}})\geq
0$. Nonetheless, we conjecture that stacking nested regressions without this
artificial constraint does lead to a larger and more predictive model. Some
partial results in this direction are presented in the next section.

Furthermore, while we have analyzed stacking in the context of nested
regressions, the benefits of stacking have been empirically observed in much
more general settings without any nested structure. The shrinkage effect of
stacking likely continues to play some role in improving prediction
performance, but it is surely not the only factor at play. We hypothesize that
stacking in these settings benefits from combining different models, some more
complex than $\hat{f}_{\text{best}}$, where each captures separate nuances in
the data.

##  7 Size of the Stacked Model

In this section, we will provide theoretical evidence of the aforementioned
phenomenon observed by Breiman and show that, in addition to having a sum less
than one, for certain stacking weights, the stacked model puts nonzero weight
on some constituent models that have higher complexity than the best single
model when $\lambda=2$ (e.g., the model selected by AIC, SURE, or Mallows’s
$C_{p}$). To this end, note that the quantity inside the expectation in the
right hand side of (7) is an unbiased estimator of the expected prediction
error of the stacked model with adaptive weights $\hat{\boldsymbol{\alpha}}$,
where $\hat{\boldsymbol{\alpha}}$ minimizes
$R(\boldsymbol{\alpha})+(2\sigma^{2}/n)\text{df}(\boldsymbol{\alpha})$ subject
to $\boldsymbol{\alpha}\geq\mathbf{0}$. Dropping the negative term involving
$0\leq\hat{\boldsymbol{\alpha}}^{\mathrm{T}}\hat{\boldsymbol{\alpha}}_{0}\ll\|\hat{\boldsymbol{\alpha}}\|_{\ell_{0}}$,
we may thus use

$R(\boldsymbol{\alpha})+\frac{2\sigma^{2}}{n}\text{df}(\boldsymbol{\alpha})+\frac{4\sigma^{2}}{n}\|\boldsymbol{\alpha}\|_{\ell_{0}},$

to hopefully better estimate the population risk of an (adaptive) stacked
model, and solve the program

$\displaystyle\text{minimize}\quad
R(\boldsymbol{\alpha})+\frac{2\sigma^{2}}{n}\text{df}(\boldsymbol{\alpha})+\frac{4\sigma^{2}}{n}\|\boldsymbol{\alpha}\|_{\ell_{0}}$
(18) $\displaystyle\text{subject to}\quad\boldsymbol{\alpha}\geq\mathbf{0}.$

We see that the objective function penalizes the number of included models
through $\|\boldsymbol{\alpha}\|_{\ell_{0}}$ as well as the complexity of the
included models through $\text{df}(\boldsymbol{\alpha})$. The optimization
problem (18) is equivalent to (weighted) reduced isotonic regression (Gao et
al., 2020; Haiminen et al., 2008) in which one fits an isotonic sequence
subject to a constraint on the number of distinct values it can assume, i.e.,

$\displaystyle\text{minimize}\quad\sum_{k=1}^{M}w_{k}(z_{k}-\beta_{k})^{2}+\frac{4\sigma^{2}}{n}\sum_{k=1}^{M}\mathbf{1}(\beta_{k}\neq\beta_{k+1})$
(19) $\displaystyle\text{subject
to}\quad\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{M}\leq 1,$

where $w_{k}=\Delta R_{k}>0$ and $z_{k}=(\sigma^{2}/n)(\Delta d_{k}/\Delta
R_{k})>0$.

Reduced isotonic regression can be implemented in $O(M^{3})$ time using a two-
stage approach consisting of Pooled Adjacent Violators Algorithm (PAVA)
(Barlow et al., 1972) followed by dynamic programming for the $K$-segmentation
problem (Bellman, 1961), e.g., see the procedure of (Haiminen et al., 2008).
The scar package in R provides a function for reduced isotonic regression.

It turns out that the stacked model with weights from solving (18) will always
assign nonzero weight to models with dimension at least as large as that of
the best single model, agreeing with what Breiman observed using cross-
validated data.

######  Theorem 7.1.

The dimension of the stacked model with weights from (18) is no less than the
dimension of the best single model (9) with $\lambda=2$; that is,

$\text{dim}(\hat{f}_{\text{stack}})\geq\text{dim}(\hat{f}_{\text{best}}).$

Furthermore, the weights sum to less than one, i.e.,
$\sum_{k=1}^{M}\hat{\alpha}_{k}<1$.

##  8 Other Connections

In this section, we discuss connections between stacking and other ways of
combining estimators.

###  8.1 Aggregation

Another line of research known as _aggregation_ covers similar ground, albeit
without an examination of shrinkage and model complexity. The goal of this
research also differs in that it benchmarks against the error of the best
single oracle model, namely, $\min_{1\leq k\leq
M}\mathbb{E}\big{[}\|\mathbf{f}-\boldsymbol{\hat{\mu}}_{k}\|^{2}]$, rather
than what is achieved by conventional model selection criteria, namely,
$\hat{f}_{\text{best}}$.

In this case, the weights assigned to the model (15) are additionally
constrained to sum to one, and the risk bounds that are sought take the form

$\min_{1\leq k\leq
M}\mathbb{E}\big{[}\|\mathbf{f}-\boldsymbol{\hat{\mu}}_{k}\|^{2}]+\frac{\Lambda\sigma^{2}}{n},$
(20)

where $\Lambda$ is a constant that possibly depends on $M$. The term
$\Lambda\sigma^{2}/n$ is the price to pay for aggregating the various models
and from not knowing the oracle model that minimizes the population risk.
Evaluating the relative performance of $\hat{f}_{\text{stack}}$ and
$\hat{f}_{\text{best}}$ using this bound presents challenges, as existing
inequalities in the literature only allow for a comparison of the risks of
$\hat{f}_{\text{best}}$ and the best single oracle model up to constants
(Cavalier et al., 2002). For an unstructured collection of models, the
constant $\Lambda$ in the bound can be taken to be of order $\log(M)$ (Leung
and Barron, 2006; Dai et al., 2012; Bellec, 2018). However, Bellec and Yang
(2020) have recently demonstrated that if the constituent models are ordered
linear smoothers (Kneip, 1994), which includes the nested linear least squares
models examined in this paper, then the constant is universal. Bellec (2018)
provided a best combination oracle inequality for aggregation over the simplex
similar to Theorem 4.2, though paying a price of $\sigma^{2}\sqrt{\log(M)/n}$,
uniformly over all models, instead of the distribution-dependent
$4\sigma^{2}\mathbb{E}[\|\hat{\boldsymbol{\alpha}}\|_{\ell_{0}}]/n$.

Inspired by the $Q$-aggregation procedure (Dai et al., 2012; Rigollet, 2012),
recent work (Bellec, 2018; Bellec and Yang, 2020) adds an additional term to
the unbiased estimator of the population risk
$R(\boldsymbol{\alpha})+(2\sigma^{2}/n)\text{df}(\boldsymbol{\alpha})$,
namely,
$\eta\sum_{m=1}^{M}\alpha_{m}\big{\|}\boldsymbol{\hat{\mu}}_{m}-\sum_{k=1}^{M}\alpha_{k}\boldsymbol{\hat{\mu}}_{k}\big{\|}^{2}$
with $\eta=1/2$, which encourages the solution to be closer to the vertices of
the probability simplex (the search space for $\boldsymbol{\alpha}$). While
this does lead to sparser solutions, the mechanism by which it is done is
quite different from (19), which in contrast explicitly penalizes the number
of nonzero weights. To see why, suppose the individual models are linear least
squares projections onto nested subspaces. Simple algebra (see the appendix)
shows that the program from (Bellec, 2018; Bellec and Yang, 2020) is
equivalent to

$\displaystyle\text{minimize}\quad\sum_{k=1}^{M}w_{k}(z_{k}-\beta_{k})^{2}$
(21) $\displaystyle\text{subject to}\quad
0=\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{M}\leq 1,$

where $w_{k}=\Delta R_{k}$ and
$z_{k}=(1/(1-\eta))(\eta/2+(\sigma^{2}/n)(\Delta d_{k}/\Delta R_{k}))$. The
solution has a closed form expression that we derive in the appendix:

$\check{\beta}_{1}=0,\quad\check{\beta}_{k}=1-\phi\Bigg{(}\frac{1}{1-\eta}\Bigg{(}1-\eta/2-\frac{\sigma^{2}}{n}\min_{k\leq
i\leq M}\max_{1\leq j<k}\frac{d_{i}-d_{j}}{R_{j}-R_{i}}\Bigg{)}\Bigg{)},\quad
k=2,3,\dots,M,$ (22)

where $\phi(z)=\min\\{1,\max\\{0,z\\}\\}$ is the clip function for $z$ at 0
and 1. As can be seen from (22) and the fact that
$\phi\big{(}\pm\frac{\text{const.}}{1-\eta}\big{)}\in\\{0,1\\}$ for
sufficiently large $\eta$, the additional term tends to push
$\check{\beta}_{k}$ into the corners $\\{0,1\\}$, making it more likely that
$\check{\alpha}_{k}=\check{\beta}_{k+1}-\check{\beta}_{k}=0$. However, there
are other ways to achieve sparsity beyond
$\check{\beta}_{k}=\check{\beta}_{k+1}=0$ or
$\check{\beta}_{k}=\check{\beta}_{k+1}=1$. For example, whenever
$\check{\beta}_{k}=\check{\beta}_{k+1}$, i.e., the isotonic solution is
constant across two successive values, $\check{\alpha}_{k}$ will be zero.
Thus, sparsity can also be modulated by controlling the number of constant
pieces of the isotonic solution $\check{\beta}_{k}$—precicely what the
objective function in program (19) seeks to do.

###  8.2 Randomized Ensembles

A hallmark of random forests is that the optimization process for choosing
good constituent tree models (i.e., splitting variables) is randomized. That
is, instead of selecting the best splitting variable among all variables at
each node, only a randomly generated subset of mtry of them are considered
each time. We can use a similar idea here in the context of nested regression
models. To this end, let $m$ be an integer less than or equal to $M$, which
can be thought of as analogous to the mtry parameter from random forests.
Generate a random subset $I$ from $\\{1,2,\dots,M\\}$ with cardinality $m-1$
and set $J=\\{0\\}\cup I$. Define

$\hat{m}\in\operatorname*{argmin}_{k\in
J}\;R_{k}+\lambda\frac{\sigma^{2}d_{k}}{n},$

so that $\hat{\mu}_{\hat{m}}$ is the best single model among
$\\{\hat{\mu}_{k}:k\in J\\}$. Repeating this process independently $B$ times
and averaging the resulting model predictions yields a randomized ensemble
estimator

$\hat{f}_{\text{rand}}(\mathbf{x})=\frac{1}{B}\sum_{b=1}^{B}\hat{\mu}_{\hat{m}_{b}}(\mathbf{x}).$

Using the form (16) of the best single model but specialized to a subset of
the models, we may write

$\hat{\mu}_{\hat{m}}(\mathbf{x})=\sum_{k=1}^{M}(\hat{\mu}_{k}(\mathbf{x})-\hat{\mu}_{k-1}(\mathbf{x}))\mathbf{1}(\hat{\gamma}_{k}(I)<1/\lambda),$

where $\hat{\gamma}_{k}(I)$ is the isotonic sequence (13) with knots in $I$,
that is, for each $k\in I$,

$\hat{\gamma}_{k}(I)=\frac{\sigma^{2}}{n}\min_{i\in I,\;i\geq k}\;\max_{j\in
J,\;j<k}\frac{d_{i}-d_{j}}{R_{j}-R_{i}},$

and $\hat{\gamma}_{k}(I)$ is constant for $k\notin I$ between successive
indices in $I$, taking on the smallest of the two adjacent
$\hat{\gamma}_{l}(I)$ values. When $B$ is large, we have the approximation

$\hat{f}_{\text{rand}}(\mathbf{x})\approx\sum_{k=1}^{M}(\hat{\mu}_{k}(\mathbf{x})-\hat{\mu}_{k-1}(\mathbf{x}))\mathbb{P}_{I}\big{(}\hat{\gamma}_{k}(I)<1/\lambda\big{)},$
(23)

where $\mathbb{P}_{I}$ denotes the probability measure induced from random
sampling. As can be seen from comparing (23) and (15), the randomized ensemble
and stacking both reduce the instability in the best single model, but do so
in different ways. Recall the form of the best single estimator (16) in which
the successive predictive differences
$\hat{\mu}_{k}(\mathbf{x})-\hat{\mu}_{k-1}(\mathbf{x})$ across successive
submodels are multiplied by $\mathbf{1}(\hat{\gamma}_{k}<1/\lambda)$ and then
summed up to the selected model index $\hat{m}$. The randomized ensemble
smooths the jump discontinuity in the indicator function
$\mathbf{1}(\hat{\gamma}_{k}<1/\lambda)$ by instead averaging over the
indicator functions $\mathbf{1}(\hat{\gamma}_{k}(I)<1/\lambda)$ resulting from
many (locally) best single models. On the other hand, stacking shrinks
$\mathbf{1}(\hat{\gamma}_{k}<1/\lambda)$ towards zero, via
$(1-\tau\hat{\gamma}_{k})\mathbf{1}(\hat{\gamma}_{k}<1/\lambda)$, making the
jump discontinuity less pronounced. An interesting question for future work
would be to analytically compare the factors
$\mathbb{P}_{I}\big{(}\hat{\gamma}_{k}(I)<1/\lambda\big{)}$ from (23) and
$(1-\tau\hat{\gamma}_{k})\mathbf{1}(\hat{\gamma}_{k}<1/\lambda)$ from (15) on
predictive accuracy.

We mention in closing that previous work (Bühlmann and Yu, 2002) has also
revealed the smoothing effect of ensembling on indicator functions, but from
subsampling the data (i.e., bagging) as opposed to subsampling the candidate
models.

##  9 Discussion and Conclusion

The findings from this research provide compelling evidence for the enhanced
predictive accuracy of stacked regressions. As we continually seek cheap and
simple ways to improve model performance, ensemble methods like stacking
emerge as valuable tools. Our research reiterates the seminal works of Wolpert
(1992) and Breiman (1996), both of whom touched upon the potential of
stacking.

The fact that the stacked model provably outperforms the best single model for
_any_ regression function $f$, noise level $\sigma^{2}$, and sample size $n$
is of course noteworthy. However, it is essential to consider the conditions
under which this advantage is true. Perhaps most limiting was that we stack
models resulting from least squares projections onto nested subspaces and that
the errors are Gaussian. For tractable theoretical analysis, it is likely the
Gaussian assumption cannot be dropped, but the structure of the models can
most likely be extended. For example, Breiman also considered stacking ridge
regressions and reached similar though not quite as strong conclusions on its
merits. More generally, ridge and nested regression models are both instances
of ordered linear smoothers (Kneip, 1994), i.e., families of models that can
be written as

$\hat{\mu}_{k}(\mathbf{x})=\sum_{l=1}^{n}c_{k,l}\langle\mathbf{y},\boldsymbol{\psi}_{l}\rangle\psi_{l}(\mathbf{x}),\quad
0\leq c_{k,l}\leq c_{k+1,l}\leq 1.$

Interestingly, the representation (16) for the best single model still holds
if $d_{k}$ is replaced by $\text{df}(\hat{\mu}_{k})=\sum_{l=1}^{n}c_{k,l}$.
However, it is not true anymore that the stacking estimator has the form (15),
unless $c_{k,l}\in\\{0,1\\}$. We leave these considerations for future work.

##  Appendix A Proofs

In this section, we provide full proofs of all the main results and other
statements in the main text. We first establish some notation. Given a
sequence of indices $s_{0}=0,s_{1},\dots,s_{u}=M$, we define $\Delta
d_{s_{l}}=d_{s_{l}}-d_{s_{l-1}}$ and $\Delta R_{s_{l}}=R_{s_{l-1}}-R_{s_{l}}$.

###### Proof of Lemma 5.1.

The proof is based on the following lemma.

######  Lemma A.1.

Suppose $h(\mathbf{x})$ has the form

$\displaystyle
h(\mathbf{x})=\sum_{k=1}^{M}(\hat{\mu}_{k}(\mathbf{x})-\hat{\mu}_{k-1}(\mathbf{x}))c_{k}.$

Then,

$\displaystyle\left\|{\mathbf{y}-\mathbf{h}}\right\|^{2}=R_{0}+\sum_{k=1}^{M}\Delta
R_{k}(c_{k}^{2}-2c_{k}).$

###### Proof of Lemma A.1.

Note that for any $k\neq l$,

$\displaystyle\langle\boldsymbol{\hat{\mu}}_{k}-\boldsymbol{\hat{\mu}}_{k-1},\boldsymbol{\hat{\mu}}_{l}-\boldsymbol{\hat{\mu}}_{l-1}\rangle=\Bigg{\langle}{\sum_{i\in
A_{k}\backslash
A_{k-1}}\left\langle\mathbf{y},\boldsymbol{\psi}_{i}\right\rangle\psi_{i}},{\sum_{j\in
A_{l}\backslash
A_{l-1}}\langle\mathbf{y},\boldsymbol{\psi}_{j}\rangle\psi_{j}}\Bigg{\rangle}=0,$

since $\langle\boldsymbol{\psi}_{i},\boldsymbol{\psi}_{j}\rangle=0$ for any
$i\neq j$, and

$\displaystyle\|\boldsymbol{\hat{\mu}}_{k}-\boldsymbol{\hat{\mu}}_{k-1}\|^{2}=\sum_{j\in
A_{k}\backslash
A_{k-1}}\langle\mathbf{y},\boldsymbol{\psi}_{j}\rangle^{2}=\Delta R_{k}.$

In addition,

$\displaystyle\langle\mathbf{y},\boldsymbol{\hat{\mu}}_{k}-\boldsymbol{\hat{\mu}}_{k-1}\rangle=\Bigg{\langle}\mathbf{y},\sum_{j\in
A_{k}\backslash
A_{k-1}}\langle\mathbf{y},\boldsymbol{\psi}_{j}\rangle\psi_{j}\Bigg{\rangle}=\sum_{j\in
A_{k}\backslash A_{k-1}}\langle y,\boldsymbol{\psi}_{j}\rangle^{2}=\Delta
R_{k}.$

Thus,

$\displaystyle\left\|{\mathbf{y}-\mathbf{h}}\right\|^{2}$
$\displaystyle=\left\|{\mathbf{y}-\sum_{k=1}^{M}(\boldsymbol{\hat{\mu}}_{k}-\boldsymbol{\hat{\mu}}_{k-1})c_{k}}\right\|^{2}$
$\displaystyle=\|\mathbf{y}\|^{2}+\sum_{k=1}^{M}c_{k}^{2}\|\boldsymbol{\hat{\mu}}_{k}-\boldsymbol{\hat{\mu}}_{k-1}\|^{2}-2\sum_{k=1}^{M}c_{k}\langle\mathbf{y},\boldsymbol{\hat{\mu}}_{k}-\boldsymbol{\hat{\mu}}_{k-1}\rangle$
$\displaystyle=R_{0}+\sum_{k=1}^{M}\Delta R_{k}(c_{k}^{2}-2c_{k}).\qed$

Let $c_{k}=\sum_{i=k}^{M}\alpha_{i}$, for $k=1,2,\dots,M$. Define
$\beta_{k}=1-c_{k}$, with $\beta_{M+1}=1$. Note that
$\alpha_{k}=\beta_{k+1}-\beta_{k}$. Using summation by parts, we know

$\displaystyle\sum_{k=1}^{M}\alpha_{k}\hat{\mu}_{k}(\mathbf{x})=\sum_{k=1}^{M}(\hat{\mu}_{k}(\mathbf{x})-\hat{\mu}_{k-1}(\mathbf{x}))c_{k},\quad\sum_{k=1}^{M}\alpha_{k}d_{k}=\sum_{k=1}^{M}c_{k}\Delta
d_{k},$

where $\hat{\mu}_{0}\equiv 0$. In addition, note that if there exists some
$1\leq a\leq M$ such that $\alpha_{a}>0=\alpha_{a+1}=\cdots\alpha_{M}$, then
we have $c_{a}>0=c_{a+1}=\cdots=c_{M}$. Thus,

$\displaystyle\text{dim}(\boldsymbol{\alpha})=\sum_{k=1}^{M}\Delta
d_{k}\textbf{1}(c_{k}\neq 0)=\sum_{k=1}^{M}\Delta
d_{k}\textbf{1}(\beta_{k}\neq 1).$

Therefore, with Lemma A.1 we can establish the equivalence between the
solution of (8) and the solution of the following programs:

$\displaystyle\text{minimize}\;\;R_{0}-\sum_{k=1}^{M}\Delta
R_{k}(2c_{k}-c_{k}^{2})+\frac{2\tau\sigma^{2}}{n}\sum_{k=1}^{M}c_{k}\Delta
d_{k}+\xi\sum_{k=1}^{M}\Delta d_{k}\textbf{1}(c_{k}\neq 0)$
$\displaystyle\text{subject to}\;\;c_{1}\geq c_{2}\geq\cdots\geq c_{M}\geq 0$
$\displaystyle\Leftrightarrow~{}$
$\displaystyle\text{minimize}\;\;\sum_{k=1}^{M}\Delta
R_{k}\Bigg{(}c_{k}^{2}-2c_{k}+c_{k}\frac{2\tau\sigma^{2}}{n}\frac{\Delta
d_{k}}{\Delta R_{k}}\Bigg{)}+\xi\sum_{k=1}^{M}\Delta d_{k}\textbf{1}(c_{k}\neq
0)$ $\displaystyle\text{subject to}\;\;c_{1}\geq c_{2}\geq\cdots\geq c_{M}\geq
0$ $\displaystyle\Leftrightarrow~{}$
$\displaystyle\text{minimize}\;\;\sum_{k=1}^{M}\Delta
R_{k}\Bigg{(}c_{k}-\Bigg{(}1-\frac{\tau\sigma^{2}}{n}\frac{\Delta
d_{k}}{\Delta R_{k}}\Bigg{)}\Bigg{)}^{2}+\xi\sum_{k=1}^{M}\Delta
d_{k}\textbf{1}(c_{k}\neq 0)$ $\displaystyle\text{subject to}\;\;c_{1}\geq
c_{2}\geq\cdots\geq c_{M}\geq 0$ $\displaystyle\Leftrightarrow~{}$
$\displaystyle\text{minimize}\;\;\sum_{k=1}^{M}w_{k}(z_{k}-\beta_{k})^{2}+\xi\sum_{k=1}^{M}\Delta
d_{k}\textbf{1}(\beta_{k}\neq 1)$ $\displaystyle\text{subject
to}\;\;\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{M}\leq 1,$

where $w_{k}=\Delta R_{k}>0$, $z_{k}=(\tau\sigma^{2}/n)(\Delta d_{k}/\Delta
R_{k})>0$, and $\xi=(\sigma^{2}/n)((\lambda-\tau)^{2}_{+}/\lambda)$. ∎

###### Proof of Theorem 5.2.

We first show the formula for the best single model (16). It can be directly
obtained from the following lemma, which characterizes the selected index
$\hat{m}$.

######  Lemma A.2.

Let
$\tilde{m}=\max\big{\\{}k\in\\{0,1,\dots,M\\}:\hat{\gamma}_{k}<1/\lambda\big{\\}}$
where $\\{\hat{\gamma}_{k}\\}$ is defined in (13). Then, for any $0\leq
j<\tilde{m}$,

$\displaystyle
R_{\tilde{m}}+\frac{\lambda\sigma^{2}}{n}d_{\tilde{m}}<R_{j}+\frac{\lambda\sigma^{2}}{n}d_{j}.$

and for any $\tilde{m}<j\leq M$

$\displaystyle R_{\tilde{m}}+\frac{\lambda\sigma^{2}}{n}d_{\tilde{m}}\leq
R_{j}+\frac{\lambda\sigma^{2}}{n}d_{j}.$

This shows $\hat{m}=\tilde{m}$.

###### Proof of Lemma A.2.

A well-known geometric property of isotonic regression is that the solution is
the left slope of the _greatest convex minorant_ of the _cumulative sum
diagram_ (Barlow et al., 1972). Let $0=s_{0}<s_{1}<\dots<s_{u}=M$ be the
values of $k$ for which $\hat{\gamma}_{s_{k}+1}>\hat{\gamma}_{s_{k}}$ ($u\geq
k\geq 0$), and $\hat{\gamma}_{0}=0$ and $\hat{\gamma}_{M+1}=\infty$, implying
that $s_{0},s_{1},\dots,s_{u}$ are the change points of the isotonic solution.
Then, according to this property,

$\displaystyle\hat{\gamma}_{s_{l}}=\frac{\sigma^{2}}{n}\frac{d_{s_{l}}-d_{s_{l-1}}}{R_{s_{l-1}}-R_{s_{l}}},\quad
l=1,2,\dots,u.$ (24)

The definition of $\tilde{m}$ implies it should be one of the change points.
We assume $\tilde{m}=s_{a}$. According to (24) and the definition of
$\hat{\gamma}_{s_{a}}$, we know for $0\leq j<\tilde{m}$,

$\displaystyle\frac{\sigma^{2}}{n}\frac{d_{s_{a}}-d_{s_{a-1}}}{R_{s_{a-1}}-R_{s_{a}}}=\hat{\gamma}_{s_{a}}=\max_{1\leq
i\leq
s_{a}-1}\frac{\sigma^{2}}{n}\frac{d_{s_{a}}-d_{i}}{R_{i}-R_{s_{a}}}\geq\frac{\sigma^{2}}{n}\frac{d_{s_{a}}-d_{j}}{R_{j}-R_{s_{a}}}.$

Thus,

$\displaystyle R_{s_{a}}+\frac{\lambda\sigma^{2}}{n}d_{s_{a}}$
$\displaystyle\leq
R_{j}+R_{s_{a}}-R_{j}+\frac{\lambda\sigma^{2}}{n}d_{j}+\lambda\hat{\gamma}_{s_{a}}(R_{j}-R_{s_{a}})$
$\displaystyle=R_{j}+\frac{\lambda\sigma^{2}}{n}d_{j}+(R_{j}-R_{s_{a}})(\lambda\hat{\gamma}_{s_{a}}-1)$
$\displaystyle<R_{j}+\frac{\lambda\sigma^{2}}{n}d_{j}.$

where the last inequality holds since
$\hat{\gamma}_{s_{a}}=\hat{\gamma}_{\tilde{m}}<1/\lambda$. Similarly, for any
$\tilde{m}<j\leq M$,

$\displaystyle\frac{\sigma^{2}}{n}\frac{d_{s_{a+1}}-d_{s_{a}}}{R_{s_{a}}-R_{s_{a+1}}}=\hat{\gamma}_{s_{a}+1}=\min_{s_{a}+1\leq
i\leq
M}\frac{\sigma^{2}}{n}\frac{d_{i}-d_{s_{a}}}{R_{s_{a}}-R_{i}}\leq\frac{\sigma^{2}}{n}\frac{d_{j}-d_{s_{a}}}{R_{s_{a}}-R_{j}}.$

Thus,

$\displaystyle R_{s_{a}}+\frac{\lambda\sigma^{2}}{n}d_{s_{a}}$
$\displaystyle\leq
R_{j}+R_{s_{a}}-R_{j}+\frac{\lambda\sigma^{2}}{n}d_{j}-\lambda\hat{\gamma}_{s_{a}+1}(R_{s_{a}}-R_{j})$
$\displaystyle=R_{j}+\frac{\lambda\sigma^{2}}{n}d_{j}+(R_{s_{a}}-R_{j})(1-\lambda\hat{\gamma}_{s_{a}+1})$
$\displaystyle<R_{j}+\frac{\lambda\sigma^{2}}{n}d_{j}.$

where the last inequality holds since
$\hat{\gamma}_{s_{a}+1}=\hat{\gamma}_{\tilde{m}+1}\geq 1/\lambda$ and
$R_{s_{a}}-R_{j}>0$. ∎

Next we show (15). According to Lemma 5.1, we can focus on program (12). We
need the following lemmas: Lemma A.3 solves the weighted isotonic regression
problem with bounded constraints and Lemma A.4 paves the way for establishing
almost sure uniqueness.

######  Lemma A.3.

Let $w_{k}>0$ ($k=1,2,\dots,M$). Then the solution to the program

$\displaystyle\text{minimize}\quad\sum_{k=1}^{M}w_{k}\left(z_{k}-\mu_{k}\right)^{2}$
(25) $\displaystyle\text{subject to}\quad
a\leq\mu_{1}\leq\cdots\leq\mu_{M}\leq b,$

is $\tilde{\mu}_{k}=\phi_{a,b}\Big{(}\max\limits_{1\leq i\leq
k}\min\limits_{k\leq j\leq
M}\frac{\sum_{l=i}^{j}w_{l}z_{l}}{\sum_{l=i}^{j}w_{l}}\Big{)}$ for
$k=1,2,\dots,M$, where $\phi_{a,b}(z)=\min\\{b,\max\\{a,z\\}\\}$ is the clip
function for $z$ at $a$ and $b$. In particular, when $a=0$ and $b=\infty$, we
have $\tilde{\mu}_{k}=\Big{(}\max\limits_{1\leq i\leq k}\min\limits_{k\leq
j\leq M}\frac{\sum_{l=i}^{j}w_{l}z_{l}}{\sum_{l=i}^{j}w_{l}}\Big{)}_{+}$ for
$k=1,2,\dots,M$.

###### Proof of Lemma A.3.

We first consider another program

$\displaystyle\text{minimize}\quad\sum_{k=1}^{M}w_{k}\left(z_{k}-\mu_{k}\right)^{2}$
(26) $\displaystyle\text{subject
to}\quad\mu_{1}\leq\mu_{2}\leq\cdots\leq\mu_{M},$

A well-known result (Barlow and Brunk, 1972, Equation 2.8) shows the solution
of (LABEL:optim:original), denoted by $\boldsymbol{\mu}^{*}$, is
$\mu_{k}^{*}=\max\limits_{1\leq i\leq k}\min\limits_{k\leq j\leq
M}\frac{\sum_{l=i}^{j}w_{l}z_{l}}{\sum_{l=i}^{j}w_{l}}$ for $1\leq k\leq M$.
Let $\boldsymbol{\eta}=(\eta_{0},\eta_{1},\dots,\eta_{M})$ and
$\boldsymbol{\lambda}=(\lambda_{1},\lambda_{2},\dots,\lambda_{M-1})$ with each
entry of $\boldsymbol{\eta}$ and $\boldsymbol{\lambda}$ nonnegative. We denote
the Lagrangians of (LABEL:optim:truncate) and (LABEL:optim:original)
respectively by

$\displaystyle
L_{1}(\boldsymbol{\mu},\boldsymbol{\eta})=\sum_{k=1}^{M}w_{k}(z_{k}-\mu_{k})^{2}-\eta_{0}(\mu_{1}-a)+\sum_{k=1}^{M-1}\eta_{k}(\mu_{k}-\mu_{k+1})+\eta_{M}(\mu_{M}-b)$

and

$\displaystyle
L_{2}(\boldsymbol{\mu},\boldsymbol{\lambda})=\sum_{k=1}^{M}w_{k}(z_{k}-\mu_{k})^{2}+\sum_{k=1}^{M-1}\lambda_{k}(\mu_{k}-\mu_{k+1}).$

Because the program (LABEL:optim:original) is convex, there exists
$\boldsymbol{\lambda}^{*}=(\lambda_{1}^{*},\lambda_{2}^{*},\dots,\lambda_{M-1}^{*})$
such that $(\boldsymbol{\mu}^{*},\boldsymbol{\lambda}^{*})$ satisfies the
Karush–Kuhn–Tucker (KKT) conditions, i.e.,

$\displaystyle\nabla L_{2}(\boldsymbol{\mu}^{*},\boldsymbol{\lambda}^{*})=0$
$\displaystyle\Rightarrow-2w_{k}(z_{k}-\mu_{k}^{*})-\lambda_{k-1}^{*}+\lambda_{k}^{*}=0,\quad
k=1,2,\dots,M$ (27) $\displaystyle\boldsymbol{\lambda}^{*}$ $\displaystyle\geq
0,$ (28) $\displaystyle\mu_{1}^{*}\leq\mu_{2}^{*}\leq$
$\displaystyle\cdots\leq\mu_{M}^{*},$ (29)
$\displaystyle\lambda_{k}^{*}(\mu_{k}^{*}-\mu_{k+1}^{*})$
$\displaystyle=0,\quad k=1,2,\dots,M-1,$ (30)

where we let $\lambda_{0}^{*}=\lambda_{M}^{*}=0$. Let
$m_{1}=\max\\{k\in[M]:\mu_{k}^{*}<a\\}$ and
$m_{2}=\max\\{k\in[M]:\mu_{k}^{*}\leq b\\}$. If $\mu_{k}^{*}\geq a$ for any
$k=1,2,\dots,M$, then let $m_{1}=0$. Since $\\{\mu_{k}^{*}\\}$ is not
decreasing and $a\leq b$, we have $m_{1}\leq m_{2}$. By (30) and the facts
that $\mu_{m_{1}+1}^{*}-\mu_{m_{1}}^{*}>0$ and
$\mu_{m_{2}+1}^{*}-\mu_{m_{2}}^{*}>0$, we know
$\lambda_{m_{1}}^{*}=\lambda_{m_{2}}^{*}=0$. [Here we define
$\mu_{M+1}^{*}=\infty$.] Now we construct an
$\boldsymbol{\eta}^{*}=(\eta_{0}^{*},\eta_{1}^{*},\dots,\eta_{M}^{*})$
satisfying

$\displaystyle\eta_{k}^{*}=\begin{cases}2w_{k}z_{k}-2w_{k}b+\eta_{k-1}^{*}&\quad
k=m_{2}+1,m_{2}+2,\dots,M\\\ \lambda_{k}^{*}&\quad
k=m_{1},m_{1}+1,\dots,m_{2}\\\ -2w_{k+1}z_{k+1}+2w_{k+1}a+\eta_{k+1}^{*}&\quad
k=0,1,\dots,m_{1}-1.\end{cases}$

In what follows, we verify that
$(\boldsymbol{\tilde{\mu}},\boldsymbol{\eta}^{*})$ satisfies the KKT
conditions for (LABEL:optim:truncate). Since $\mu_{k}^{*}$ is increasing, it
is easy to see $a\leq\tilde{\mu}_{1}\leq\cdots\tilde{\mu}_{M}\leq b$, which
confirms the primal feasibility. Note that
$\eta_{k}^{*}(\tilde{\mu}_{k}-\tilde{\mu}_{k+1})=\lambda_{k}^{*}(\mu_{k}^{*}-\mu_{k+1}^{*})=0$
for $m_{2}>k>m_{1}$. Since
$\eta_{m_{1}}^{*}=\lambda_{m_{1}}^{*}=\eta_{m_{2}}^{*}=\lambda_{m_{2}}^{*}=0$
and $\tilde{\mu}_{k}=a$ ($k\leq m_{1}$), and $\tilde{\mu}_{k}=b$ ($k>m_{2}$),
we also know $\eta_{k}^{*}(\tilde{\mu}_{k}-\tilde{\mu}_{k+1})=0$ for $k\leq
m_{1}$ and $k\geq m_{2}$. These confirm the complementary slackness.
Furthermore, for $m_{2}\geq k>m_{1}$,

$\displaystyle\frac{\partial}{\partial\mu_{k}}L_{1}(\boldsymbol{\tilde{\mu}},\boldsymbol{\eta}^{*})$
$\displaystyle=-2w_{k}(z_{k}-\tilde{\mu}_{k})-\eta_{k-1}^{*}+\eta_{k}^{*}=-2w_{k}(z_{k}-\mu_{k}^{*})-\lambda_{k-1}^{*}+\lambda_{k}^{*}=0,$

and for $k\leq m_{1}$,

$\displaystyle\frac{\partial}{\partial\mu_{k}}L_{1}(\boldsymbol{\tilde{\mu}},\boldsymbol{\eta}^{*})$
$\displaystyle=-2w_{k}(z_{k}-\tilde{\mu}_{k})-\eta_{k-1}^{*}+\eta_{k}^{*}=-2w_{k}z_{k}+2w_{k}a-\eta_{k-1}^{*}+\eta_{k}^{*}=0,$

where we use the definition of $\eta^{*}_{k}$ and $\tilde{\mu}_{k}=a$ ($k\leq
m_{1}$) for the above equality. Similarly, for $k>m_{2}$, we also have
$\frac{\partial}{\partial\mu_{k}}L_{1}(\boldsymbol{\tilde{\mu}},\boldsymbol{\eta}^{*})=0$.
These lead to $\nabla
L_{1}(\boldsymbol{\tilde{\mu}},\boldsymbol{\eta}^{*})=0$. Lastly, to show
$\boldsymbol{\eta}^{*}\geq 0$, it suffices to show $\eta^{*}_{k}\geq 0$ for
$k\leq m_{1}-1$ and $k\geq m_{2}+1$ since $\eta^{*}_{k}=\lambda^{*}_{k}\geq 0$
($m_{2}\geq k\geq m_{1}$). Note that

$\displaystyle\eta^{*}_{m_{1}-1}=-2w_{m_{1}}z_{m_{1}}+2w_{m_{1}}a+\eta^{*}_{m_{1}}\geq-2w_{m_{1}}z_{m_{1}}+\eta^{*}_{m_{1}}+2w_{m_{1}}\mu^{*}_{m_{1}}=\lambda^{*}_{{m_{1}}-1}\geq
0,$ (31)

where we use facts that $\mu^{*}_{m_{1}}<a$ and
$\eta_{m_{1}}^{*}=\lambda_{m_{1}}^{*}$, and the stationary condition (27) for
program (LABEL:optim:original). In addition,

$\displaystyle\eta_{{m_{1}}-2}^{*}=-2w_{{m_{1}}-1}z_{{m_{1}}-1}+2w_{m_{1}-1}a+\eta^{*}_{m_{1}-1}\geq-2w_{m_{1}-1}z_{m_{1}-1}+\lambda^{*}_{m_{1}-1}+2w_{m_{1}-1}\mu^{*}_{m_{1}-1}=\lambda^{*}_{m_{1}-2}\geq
0,$

where for the first inequality we use (31) and the fact that
$\mu^{*}_{m_{1}-1}\leq\mu^{*}_{m_{1}}<0$. For the last equality we use (27)
again. Therefore, with an induction argument we can establish
$\eta_{k}^{*}\geq 0$ for $k\leq m_{1}-1$. Similarly, we know $\eta_{k}^{*}\geq
0$ for $k\geq m_{2}+1$, which implies $\boldsymbol{\eta}^{*}\geq 0$. In
conclusion, we show that $(\boldsymbol{\tilde{\mu}},\boldsymbol{\eta}^{*})$
satisfies the KKT conditions for program (LABEL:optim:truncate). Since
(LABEL:optim:truncate) is a convex program, we know $\boldsymbol{\tilde{\mu}}$
and $\boldsymbol{\eta}^{*}$ should be its primal and dual optimizers, which
completes the proof. ∎

######  Lemma A.4.

Given any constant $c>0$. The probability that there exists an index $1\leq
k\leq M$ such that $\hat{\gamma}_{k}=c$ is zero.

###### Proof of Lemma A.4.

We show that given any constant $c>0$, the probability of there existing two
indices $0\leq a<b\leq M$ such that
$\frac{\sigma^{2}}{n}\frac{d_{b}-d_{a}}{R_{a}-R_{b}}=c$ is zero. To this end,
given indices $a$ and $b$, note that $(n/\sigma^{2})(R_{a}-R_{b})$ follows a
noncentral chi-squared distribution with v degrees of freedom and
noncentrality parameter
$(n/\sigma^{2})(\|\mathbf{f}-\mathbf{f}_{a}\|^{2}-\|\mathbf{f}-\mathbf{f}_{b}\|^{2})$.
Therefore, we have that
$\mathbb{P}(\frac{\sigma^{2}}{n}\frac{d_{b}-d_{a}}{R_{a}-R_{b}}=c)=0$. The
lemma can thus be obtained by noticing that the countable union of measure-
zero sets has measure zero. ∎

Now we start to prove (15). We begin by considering the case that
$\lambda\leq\tau$. In this case, program (12) is reduced to the bounded
isotonic regression program

$\displaystyle\text{minimize}\quad\sum_{k=1}^{M}\Delta
R_{k}\bigg{(}\frac{\tau\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}-\beta_{k}\bigg{)}^{2}$ (32) $\displaystyle\text{subject
to}\quad\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{M}\leq 1.$

From Lemma A.3, we know its solution is
$\\{\tau\hat{\gamma}_{k}\mathbf{1}(\hat{\gamma}_{k}\leq
1/\tau)+\mathbf{1}(\hat{\gamma}_{k}>1/\tau),\;k=1,2,\dots,M\\}$, where
$\hat{\gamma}_{k}$ is defined in (13). Then the solution to program (8) is

$\displaystyle\hat{\alpha}_{k}$
$\displaystyle=(1-\tau\hat{\gamma}_{k})\mathbf{1}(\hat{\gamma}_{k}\leq
1/\tau)-(1-\tau\hat{\gamma}_{k+1})\mathbf{1}(\hat{\gamma}_{k+1}\leq 1/\tau)$
$\displaystyle=(1-\tau\hat{\gamma}_{k})_{+}-(1-\tau\hat{\gamma}_{k+1})_{+},\quad
k=1,2,\dots,M.$

We next consider the case that $\lambda>\tau$. Recall that in the proof of
Lemma A.2, we assume $s_{0}=0,s_{1},\dots,s_{u}=M$ correspond to the change
points of $\hat{\boldsymbol{\gamma}}$ and also assume
$\hat{m}=\tilde{m}=s_{a}$ where $0\leq a\leq M$. Let
$\hat{\boldsymbol{\beta}}$ be the solution to program (12) and suppose
$r_{q}=\max\big{\\{}k\in\\{0,1,\dots,M\\}:\hat{\beta}_{k}<1\big{\\}}$. Note
that we want to show almost surely
$\hat{\alpha}_{k}=(1-\tau\hat{\gamma}_{k})\mathbf{1}(\hat{\gamma}_{k}<1/\lambda)-(1-\tau\hat{\gamma}_{k+1})\mathbf{1}(\hat{\gamma}_{k+1}<1/\lambda)$,
which implies that $\hat{\alpha}_{k}=0$ and $\hat{\beta}_{k}=1$ when
$k>s_{a}$. Hence, we first show that almost surely, $r_{q}=\hat{m}=s_{a}$. The
general idea of the proof is that if $r_{q}\neq s_{a}$, we can always
construct an alternative valid solution with a smaller objective value,
thereby establishing a contradiction.

Suppose $r_{q}>s_{a}$. Note that when given $r_{q}$, program (12) is reduced
to the problem

$\displaystyle\text{minimize}\quad\sum_{k=1}^{r_{q}}\Delta
R_{k}\bigg{(}\frac{\tau\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}-\beta_{k}\bigg{)}^{2}+\frac{(\lambda-\tau)^{2}}{\lambda}\frac{\sigma^{2}}{n}\sum_{k=1}^{r_{q}}\Delta
d_{k}\mathbf{1}(\beta_{k}\neq 1)$ (33) $\displaystyle\text{subject
to}\quad\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{r_{q}}\leq 1.$

Suppose the change points that correspond to its solution
$\hat{\boldsymbol{\beta}}$ are $r_{0}=0,r_{1},\dots,r_{q}$ and also suppose
$s_{b-1}<r_{q}\leq s_{b}$ where $b>a$. The solution $\hat{\boldsymbol{\beta}}$
has the form
$\hat{\beta}_{k}=\frac{\tau\sigma^{2}}{n}\frac{d_{r_{l}}-d_{r_{l-1}}}{R_{r_{l-1}}-R_{r_{l}}}$
when $r_{l-1}<k\leq r_{l}$ $(l\leq q)$. Note that for $s_{b-1}<k\leq r_{q}$,

$\displaystyle\hat{\beta}_{k}=\frac{\tau\sigma^{2}}{n}\min_{k\leq j\leq
r_{q}}\max_{1\leq i\leq
k}\frac{d_{j}-d_{i-1}}{R_{i-1}-R_{j}}\geq\tau\hat{\gamma}_{k}\geq\tau\hat{\gamma}_{s_{a}+1}\geq\frac{\tau}{\lambda}$

where we use the fact that $k\geq s_{b-1}+1\geq s_{a}+1$. Thus, leveraging
Lemma A.4, we conclude almost surely $\tau/\lambda<\hat{\beta}_{r_{q}}<1$. Let
$\Delta R_{r_{l}}=R_{r_{l-1}}-R_{r_{l}}$ and $\Delta
d_{r_{l}}=d_{r_{l}}-d_{r_{l-1}}$ for any $l\leq q$. Note that

$\displaystyle\sum_{k=r_{q-1}+1}^{r_{q}}\Delta
R_{k}\bigg{(}\frac{\tau\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}-\hat{\beta}_{k}\bigg{)}^{2}+\frac{(\lambda-\tau)^{2}}{\lambda}\frac{\sigma^{2}}{n}(d_{r_{q}}-d_{r_{q-1}})-\sum_{k=r_{q-1}+1}^{r_{q}}\Delta
R_{k}\bigg{(}\frac{\tau\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}-1\bigg{)}^{2}$ (34) $\displaystyle=~{}\Delta
R_{r_{q}}\bigg{(}-\hat{\beta}_{r_{q}}^{2}+\frac{(\lambda-\tau)^{2}}{\tau\lambda}\hat{\beta}_{r_{q}}-1+2\hat{\beta}_{r_{q}}\bigg{)}$
$\displaystyle=~{}-\frac{\Delta
R_{r_{q}}}{\tau\lambda}(\tau\hat{\beta}_{r_{q}}-\lambda)(\lambda\hat{\beta}_{r_{q}}-\tau)>0.$

This means if we change $\hat{\beta}_{k}$ to 1 for any $r_{q-1}<k\leq r_{q}$,
then we will get a smaller objective value of program (12). This establishes a
contradiction, leading to the conclusion that $r_{q}\leq s_{a}$.

Suppose $r_{q}<s_{a}$. Again, we assume $s_{b-1}<r_{q}\leq s_{b}$ where $b\leq
a$. If there exists $b<a$ such that $r_{q}=s_{b}$, then for any $s_{b}+1\leq
k\leq s_{b+1}$, we consider changing $\hat{\beta}_{k}$ from 1 to
$\check{\beta}=\frac{\tau\sigma^{2}}{n}\frac{d_{s_{b+1}}-d_{s_{b}}}{R_{s_{b}}-R_{s_{b+1}}}$.
Following a similar inequality as (34), we know

$\displaystyle\sum_{k=s_{b}+1}^{s_{b+1}}\Delta
R_{k}\bigg{(}\frac{\tau\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}-\check{\beta}\bigg{)}^{2}+\frac{(\lambda-\tau)^{2}}{\lambda}\frac{\sigma^{2}}{n}(d_{s_{b+1}}-d_{s_{b}})-\sum_{k=s_{b}+1}^{s_{b+1}}\Delta
R_{k}\bigg{(}\frac{\tau\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}-1\bigg{)}^{2}$
$\displaystyle=~{}-\frac{R_{s_{b}}-R_{s_{b+1}}}{\tau\lambda}(\tau\check{\beta}-\lambda)(\lambda\check{\beta}-\tau)<0,$

where in the last inequality we use the fact that $\check{\beta}<\tau/\lambda$
since $b<a$. This means changing $\hat{\beta}_{k}$ $(s_{b}+1\leq k\leq
s_{b+1})$ from 1 to $\check{\beta}$ can get a smaller objective value, which
leads to a contradiction.

Thus, we can assume $s_{b-1}<r_{q}<s_{b}$ where $b\leq a$. When
$\hat{\beta}_{r_{q}}>\tau/\lambda$, with a similar argument as above, we can
change $\hat{\beta}_{k}$ to 1 for any $r_{q-1}<k\leq r_{q}$ and get a smaller
objective value of program (12). We thus only need to focus on the case of
$\hat{\beta}_{r_{q}}\leq\tau/\lambda$, and almost surely we can assume
$\hat{\beta}_{r_{q}}<\tau/\lambda$. For any $k\leq s_{b}$ we let
$s_{c-1}<k\leq s_{c}$ where $c\leq b$. Then, the solution $\hat{\beta}_{k}$ to
program (LABEL:prog:reduced) satisfies

$\tau\hat{\gamma}_{k}=\frac{\tau\sigma^{2}}{n}\min_{k\leq j\leq M}\max_{1\leq
i\leq
k}\frac{d_{j}-d_{i-1}}{R_{i-1}-R_{j}}\leq\hat{\beta}_{k}=\frac{\tau\sigma^{2}}{n}\min_{k\leq
j\leq r_{q}}\max_{1\leq i\leq
k}\frac{d_{j}-d_{i-1}}{R_{i-1}-R_{j}}\leq\frac{\tau\sigma^{2}}{n}\max_{1\leq
i\leq k}\frac{d_{s_{c}}-d_{i-1}}{R_{i-1}-R_{s_{c}}}=\tau\hat{\gamma}_{k},$
(35)

which means $\hat{\beta}_{k}=\tau\hat{\gamma}_{k}$ for any $k\leq s_{b-1}$ and
also $r_{l}=s_{l}$ for any $l\leq b-1$. Let
$\tilde{\beta}=\frac{\tau\sigma^{2}}{n}\frac{d_{s_{b}}-d_{s_{b-1}}}{R_{s_{b-1}}-R_{s_{b}}}$.
We next show that if we change $\hat{\beta}_{k}$ to $\tilde{\beta}$ for any
$s_{b-1}<k\leq s_{b}$, then we obtain a smaller objective value of program
(12). It is worth noting that $\tilde{\beta}\geq\hat{\beta}_{k}$ for any
$k\leq s_{b-1}$ since they are the solutions of isotonic regression, and
$\tilde{\beta}<\tau/\lambda<1$. Hence, the replacement of $\hat{\beta}_{k}$
$(s_{b-1}<k\leq s_{b})$ with $\tilde{\beta}$ is indeed valid. Thus, it
suffices to show

$\displaystyle\sum_{k=r_{b-1}+1=s_{b-1}+1}^{r_{q}}w_{k}(z_{k}-\hat{\beta}_{k})^{2}+\frac{(\lambda-\tau)^{2}}{\lambda}\frac{\sigma^{2}}{n}(d_{r_{q}}-d_{r_{b-1}})+\sum_{k=r_{q}+1}^{s_{b}}w_{k}(z_{k}-1)^{2}$
(36)
$\displaystyle>~{}\sum_{k=r_{b-1}+1=s_{b-1}+1}^{s_{b}}w_{k}(z_{k}-\tilde{\beta})^{2}+\frac{(\lambda-\tau)^{2}}{\lambda}\frac{\sigma^{2}}{n}(d_{s_{b}}-d_{s_{b-1}})$

where $w_{k}=\Delta R_{k}>0$ and $z_{k}=(\tau\sigma^{2}/n)(\Delta d_{k}/\Delta
R_{k})>0$, and we have constraints
$\tilde{\beta}<\hat{\beta}_{r_{b}}<\cdots<\hat{\beta}_{r_{q}}<\tau/\lambda$.
By taking
$\hat{\beta}_{k}=\frac{\tau\sigma^{2}}{n}\frac{d_{r_{l}}-d_{r_{l-1}}}{R_{r_{l-1}}-R_{r_{l}}}$
when $r_{l-1}<k\leq r_{l}$ $(l\leq q)$ and $\tilde{\beta}$ into (36), we know
(36) is equivalent to

$\sum_{u=b}^{q}\frac{(\sum_{k=r_{u-1}+1}^{r_{u}}w_{k}z_{k})^{2}}{\sum_{k=r_{u-1}+1}^{r_{u}}w_{k}}\leq\frac{(\sum_{k=s_{b-1}+1}^{s_{b}}w_{k}z_{k})^{2}}{\sum_{k=s_{b-1}+1}^{s_{b}}w_{k}}+\sum_{k=r_{q}+1}^{s_{b}}w_{k}-\frac{\lambda^{2}+\tau^{2}}{\lambda\tau}\sum_{k=r_{q}+1}^{s_{b}}w_{k}z_{k}.$
(37)

For simplicity of notation, let $p=q-b$,
$t_{i}=\frac{\sum_{k=r_{b+i-1}+1}^{r_{b+i}}w_{k}z_{k}}{\sum_{k=r_{b+i-1}+1}^{r_{b+i}}w_{k}}$
and
$c_{i}=\frac{\sum_{k=r_{b+i-1}+1}^{r_{b+i}}w_{k}}{\sum_{k=r_{b-1}+1}^{s_{b}}w_{k}}$
where $i=0,1,\dots,p$. In addition, let
$\tilde{b}=\frac{\sum_{k=r_{q}+1}^{s_{b}}w_{k}z_{k}}{\sum_{k=r_{q}+1}^{s_{b}}w_{k}}$
and
$\tilde{c}=\frac{\sum_{k=r_{q}+1}^{s_{b}}w_{k}}{\sum_{k=r_{b-1}+1}^{s_{b}}w_{k}}$.
Then (37) becomes equivalent to

$\displaystyle
c_{0}t_{0}^{2}+c_{1}t_{1}^{2}+\cdots+c_{p}t_{p}^{2}\leq(c_{0}t_{0}+c_{1}t_{1}+\cdots+c_{p}t_{p}+\tilde{c}\tilde{b})^{2}+\tilde{c}-\frac{\lambda^{2}+\tau^{2}}{\lambda\tau}\tilde{b}\tilde{c},$
(38)

and we have constraints
$\tilde{\beta}<\hat{\beta}_{r_{b}}<\cdots<\hat{\beta}_{r_{q}}<\tau/\lambda$
equivalent to $c_{0}+c_{1}+\dots+c_{p}+\tilde{c}=1$ and

$c_{0}t_{0}+c_{1}t_{1}+\cdots+c_{p}t_{p}+\tilde{c}\tilde{b}<t_{0}<\dots<t_{p}\leq\tau/\lambda.$
(39)

Note that since $r_{q}<s_{b}$, we have $0<\tilde{c}\leq 1$. Then (39) can be
rewritten as

$\tilde{b}<\frac{t_{0}-(c_{0}t_{0}+c_{1}t_{1}+\dots+c_{p}t_{p})}{\tilde{c}}\quad\text{and}\quad
t_{0}<t_{1}<\dots<t_{p}<\tau/\lambda.$ (40)

Note that we can write the right hand side of (38) as a quadratic function of
$\tilde{b}$, which we denote by $g(\tilde{b})$. Then,

$g(\tilde{b})=\tilde{c}^{2}\tilde{b}^{2}+2\tilde{c}(c_{0}t_{0}+c_{1}t_{1}+\dots+c_{p}t_{p}-\frac{\lambda^{2}+\tau^{2}}{2\lambda\tau})\tilde{b}+(c_{0}t_{0}+c_{1}t_{1}+\dots+c_{p}t_{p})^{2}+\tilde{c}.$
(41)

Since $t_{0}<\tau/\lambda<1\leq\frac{\lambda^{2}+\tau^{2}}{2\lambda\tau}$ and
$\tilde{b}\leq\frac{t_{0}-(c_{0}t_{0}+c_{1}t_{1}+\dots+c_{p}t_{p})}{\tilde{c}}$,
we know $g(\tilde{b})\geq
g\Big{(}\frac{t_{0}-(c_{0}t_{0}+c_{1}t_{1}+\dots+c_{p}t_{p})}{\tilde{c}}\Big{)}$.
Then to establish (38), it suffices to show

$\displaystyle c_{0}t_{0}^{2}+c_{1}t_{1}^{2}+\cdots+c_{p}t_{p}^{2}$
$\displaystyle<g\Big{(}\frac{t_{0}-(c_{0}t_{0}+c_{1}t_{1}+\dots+c_{p}t_{p})}{\tilde{c}}\Big{)}$
(42)
$\displaystyle=t_{0}^{2}+\tilde{c}-\frac{\lambda^{2}+\tau^{2}}{\lambda\tau}(t_{0}-c_{0}t_{0}-c_{1}t_{1}-\dots-
c_{p}t_{p}),$ (43)

which is equivalent to

$\displaystyle\sum_{j=1}^{p}c_{j}h(t_{j})<(1-c_{0})t_{0}^{2}-\frac{\lambda^{2}+\tau^{2}}{\lambda\tau}(1-c_{0})t_{0}+\tilde{c},$
(44)

where we let $h(z)=z^{2}-\frac{\lambda^{2}+\tau^{2}}{\lambda\tau}z$. Note that
$h(z)$ is also a quadratic function and $t_{0}<t_{j}<\tau/\lambda$ for any
$j=1,2,\dots,p$. Then

$\displaystyle\sum_{j=1}^{p}c_{j}h(t_{j})<\sum_{j=1}^{p}c_{j}h(t_{1}).$

In addition, note that

$\displaystyle\sum_{j=1}^{p}c_{j}h(t_{0})<(1-c_{0})t_{0}^{2}-\frac{\lambda^{2}+\tau^{2}}{\lambda\tau}(1-c_{0})t_{0}+\tilde{c}$

is equivalent to

$\displaystyle\frac{\tilde{c}}{\lambda\tau}(\tau t_{0}-\lambda)(\lambda
t_{0}-\tau)>0.$

This is satisfied by the constraint that $t_{0}<\tau/\lambda$, and thus we
prove (38). Overall, we know changing $\hat{\beta}_{k}$ to $\tilde{\beta}$ for
any $s_{b-1}<k\leq s_{b}$ can lead to a smaller objective value of program
(12), which gives a contradiction. Therefore, we proved almost surely
$r_{q}=s_{a}$. In this case, with a similar argument as (35), we know
$\hat{\beta}_{k}=\tau\hat{\gamma}_{k}$ for any $k\leq r_{q}$. Thus,

$\hat{\alpha}_{k}=(1-\tau\hat{\gamma}_{k})\mathbf{1}(\hat{\gamma}_{k}<1/\lambda)-(1-\tau\hat{\gamma}_{k+1})\mathbf{1}(\hat{\gamma}_{k+1}<1/\lambda),\quad
k=1,2,\dots,M,$

and this is the almost surely unique solution. ∎

###### Proof of Theorem 4.1.

To prove this theorem, we need the following lemma, which is an extension of
Stein’s Lemma for discontinuous functions.

######  Lemma A.5 (Corollary 1 in Tibshirani (2015)).

Let $X\sim N(\mu,\sigma^{2})$, where $\mu\in\mathbb{R}$ and $\sigma>0$. Let
$h$ be a piecewise absolutely continuous function, with discontinuity set
$\\{\delta_{1},\delta_{2},\dots,\delta_{m}\\}$, and derivative $h^{\prime}$
satisfying $\mathbb{E}[|h^{\prime}(X)|]<\infty$. Then,

$\displaystyle\frac{1}{\sigma^{2}}\mathbb{E}[(X-\mu)h(X)]=\mathbb{E}[h^{\prime}(X)]+\frac{1}{\sigma}\sum_{k=1}^{m}\phi\bigg{(}\frac{\delta_{k}-\mu}{\sigma}\bigg{)}[h(\delta_{k})_{+}-h(\delta_{k})_{-}],$

where $\phi$ is the standard normal density function, and
$h(x)_{+}=\lim_{t\downarrow x}h(t)$ and $h(x)_{-}=\lim_{t\uparrow x}h(t)$.

We first characterize
$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}$.
By the usual population risk decomposition formula,

$\displaystyle\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}=\mathbb{E}\big{[}\|\mathbf{y}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}-\sigma^{2}+\frac{2\sigma^{2}}{n}\text{df}(\hat{f}_{\text{best}}).$

For simplicity of notation, we let
$\tilde{\beta}_{k}=\mathbf{1}(\hat{\gamma}_{k}\geq 1/\lambda)$, where
$k=1,2,\dots,M$. Then
$\hat{f}_{\text{best}}(\mathbf{x})=\sum_{k=1}^{M}(\hat{\mu}_{k}(\mathbf{x})-\hat{\mu}_{k-1}(\mathbf{x}))(1-\tilde{\beta}_{k})$.
To compute $\text{df}(\hat{f}_{\text{best}})$, note that

$\displaystyle\text{df}(\hat{f}_{\text{best}})=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\text{cov}\left(\hat{f}_{\text{best}}(\mathbf{x}_{i}),y_{i}\right)=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\sum_{k=1}^{M}\sum_{l\in
A_{k}\backslash
A_{k-1}}\text{cov}\left(\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\psi_{l}(\mathbf{x}_{i})(1-\tilde{\beta}_{k}),y_{i}\right).$
(45)

In addition,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\text{cov}\left(\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\psi_{l}(\mathbf{x}_{i})(1-\tilde{\beta}_{k}),y_{i}\right)$
$\displaystyle=~{}\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\psi_{l}(\mathbf{x}_{i})y_{i}(1-\tilde{\beta}_{k})\right]-\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\psi_{l}(\mathbf{x}_{i})f(\mathbf{x}_{i})(1-\tilde{\beta}_{k})\right]$
$\displaystyle=~{}\mathbb{E}\left[\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle^{2}(1-\tilde{\beta}_{k})\right]-\mathbb{E}\left[\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\left\langle
f,\psi_{l}\right\rangle(1-\tilde{\beta}_{k})\right]$
$\displaystyle=~{}\text{cov}\left(\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle,\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle(1-\tilde{\beta}_{k})\right).$

Since $\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\sim
N(\left\langle\mathbf{f},\boldsymbol{\psi}_{l}\right\rangle,\sigma^{2}/n)$ and
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle(1-\tilde{\beta}_{k})$
is a piecewise absolutely continuous function in
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle$, we can leverage
Lemma A.5. Specifically, we have

$\displaystyle\text{df}(\hat{f}_{\text{best}})$
$\displaystyle=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\sum_{k=1}^{M}\sum_{l\in
A_{k}\backslash
A_{k-1}}\text{cov}\left(\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\psi_{l}(\mathbf{x}_{i})(1-\tilde{\beta}_{k}),y_{i}\right)$
(46) $\displaystyle=\frac{n}{\sigma^{2}}\sum_{k=1}^{M}\sum_{l\in
A_{k}\backslash
A_{k-1}}\text{cov}\left(\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle,\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle(1-\tilde{\beta}_{k})\right)$
$\displaystyle=\sum_{k=1}^{M}\sum_{l\in A_{k}\backslash
A_{k-1}}\mathbb{E}\left[\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle(1-\tilde{\beta}_{k})\right]+\text{sdf}(\hat{f}_{\text{best}})$
$\displaystyle=\sum_{k=1}^{M}\sum_{l\in A_{k}\backslash
A_{k-1}}\mathbb{E}\left[1-\hat{\beta}_{k}+\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}(1-\tilde{\beta}_{k})\right]+\text{sdf}(\hat{f}_{\text{best}})$
$\displaystyle=\sum_{k=1}^{M}\Delta
d_{k}(1-\tilde{\beta}_{k})+\text{sdf}(\hat{f}_{\text{best}}),$

where we denote $\text{sdf}(\hat{f}_{\text{best}})$ as the risk arising from
the discontinuous points of $\hat{f}_{\text{best}}$, also called the _search
degrees of freedom_ (Tibshirani, 2015). Note that in the third equality above,
it is necessary to first condition on
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l^{\prime}}\right\rangle$ for all
$l^{\prime}\neq l$, and subsequently absorb the conditional expectation. For
the sake of simplicity in notation, we will omit this step in the following
discussion. The last equality holds since in each piece, we have
$\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}(1-\hat{\beta}_{k})=0$.
Therefore, with Lemma A.1 and (46), we know

$\displaystyle\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}$
$\displaystyle=\mathbb{E}\big{[}\|\mathbf{y}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}-\sigma^{2}+\frac{2\sigma^{2}}{n}\text{df}(\hat{f}_{\text{best}})$
$\displaystyle=\mathbb{E}\left[R_{0}+\sum_{k=1}^{M}\Delta
R_{k}\left((1-\tilde{\beta}_{k})^{2}-2(1-\tilde{\beta}_{k})\right)+\sum_{k=1}^{M}\frac{2\sigma^{2}}{n}\Delta
d_{k}(1-\tilde{\beta}_{k})\right]+\frac{2\sigma^{2}}{n}\text{sdf}(\hat{f}_{\text{best}})-\sigma^{2}$
$\displaystyle=\mathbb{E}\left[R_{0}+\sum_{k=1}^{M}\Delta
R_{k}\left(\tilde{\beta}_{k}^{2}-1\right)+\sum_{k=1}^{M}\frac{2\sigma^{2}}{n}\Delta
d_{k}(1-\tilde{\beta}_{k})\right]+\frac{2\sigma^{2}}{n}\text{sdf}(\hat{f}_{\text{best}})-\sigma^{2}.$

We defer the characterization of $\text{sdf}(\hat{f}_{\text{best}})$ to a
later point. To compute
$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}$
we need the following lemma.

######  Lemma A.6.

The following results hold.

  1. (i)

For any $1\leq k\leq M$ and $l\in A_{M}$, given
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l^{\prime}}\right\rangle$ for all
$l^{\prime}\neq l$, we have $\hat{\gamma}_{k}$ is nonincreasing with
$|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|$.

  2. (ii)

For any $1\leq k\leq M$ and $l\in A_{M}$, given
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l^{\prime}}\right\rangle$ for all
$l^{\prime}\neq l$, we have
$(1-\tau\hat{\gamma}_{k})\mathbf{1}(\hat{\gamma}_{k}<\gamma)$ is piecewise
absolutely continuous with
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle$.

###### Proof of Lemma A.6.

To prove part (i), suppose that $l\in A_{a}\backslash A_{a-1}$ where $1\leq
a\leq M$. If $\hat{\gamma}_{a}<\hat{\gamma}_{k}$, then as
$|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|$ increases, since
$l\in A_{a}\backslash A_{a-1}$ and
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle^{2}$ appears in the
denominator $R_{j-1}-R_{i}$ of
$\hat{\gamma}_{a}=\frac{\sigma^{2}}{n}\min_{a\leq i\leq M}\max_{1\leq j\leq
a}\frac{d_{i}-d_{j-1}}{R_{j-1}-R_{i}}$, the $\hat{\gamma}_{a}$ will keep
decreasing while $\hat{\gamma}_{k}$ remains unchanged. If
$\hat{\gamma}_{a}\geq\hat{\gamma}_{k}$, then as
$|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|$ increases,
$\hat{\gamma}_{a}$ keeps decreasing and $\hat{\gamma}_{k}$ still remains
unchanged until $\hat{\gamma}_{a}$ reaches $\hat{\gamma}_{k}$. After this
point, they remain equal and continue to decrease.

To prove part (ii), by symmetry it suffices to show the piecewise absolutely
continuity on
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\in[0,\infty)$. Let
$\hat{\tau}_{k}=(1-\tau\hat{\gamma}_{k})\mathbf{1}(\hat{\gamma}_{k}<\gamma)$.
We first consider the case of $\lambda>\tau$. Since $\hat{\gamma}_{k}$ is the
solution of vanilla isotonic regression, it is continuous. By part (i), we
know there is at most one discontinuous point when
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\geq 0$. We denote
it by $t$ and note that it satisfies $\hat{\gamma}_{k}=1/\lambda$ when
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle=t$. Then from part
(i) we know $\hat{\tau}_{k}=0$ when
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\leq t$ and thus it
is absolutely continuous. When
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle>t$, note that we
can write
$\hat{\gamma}_{k}=\frac{\sigma^{2}}{n}\frac{d_{c}-d_{b}}{R_{b}-R_{c}}$, where
$1\leq b<c\leq M$. Here $a,b$ are two change points of
$\hat{\boldsymbol{\gamma}}$, and they are both functions of
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle$. Since $a,b$ have
only finite number of choices, we can thus rewrite $\hat{\gamma}_{k}$ as

$\displaystyle\hat{\gamma}_{k}=\sum_{i=1}^{K}\frac{\sigma^{2}}{n}\frac{d_{c_{i}}-d_{b_{i}}}{R_{b_{i}}-R_{c_{i}}}\mathbf{1}(t_{i-1}<\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\leq
t_{i})+\frac{\sigma^{2}}{n}\frac{d_{c_{K+1}}-d_{b_{K+1}}}{R_{b_{K+1}}-R_{c_{K+1}}}\mathbf{1}(t_{K}<\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle),$
(47)

where $1\leq b_{i}<c_{i}\leq M$ for any $i=1,2,\dots,K+1$ and
$t=t_{0}<t_{1}<\cdots<t_{K}$. Here $b_{i},c_{i}$ are constants independent of
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle$ for any $i$, and
$K$ is also a constant independent of
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle$. To clarify, it is
important to note that $b_{i},c_{i}$, and $K$ should depend on
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l^{\prime}}\right\rangle$ where
$l^{\prime}\neq l$. However, in this particular lemma, we are provided with
$\\{\left\langle\mathbf{y},\boldsymbol{\psi}_{l^{\prime}}\right\rangle\\}_{l^{\prime}\neq
l}$. Therefore, for the purposes of this lemma, we treat $b_{i},c_{i}$, and
$K$ as constants. Note that
$\frac{\sigma^{2}}{n}\frac{d_{c_{i}}-d_{b_{i}}}{R_{b_{i}}-R_{c_{i}}}$ is
absolutely continuous. Hence, we know $\hat{\gamma}_{k}$ is absolutely
continuous when $\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle>t$,
and so is $(1-\tau\hat{\gamma}_{k})$. Therefore, we proved $\hat{\tau}_{k}$ is
piecewise absolutely continuous when $\lambda>\tau$. In the case of
$\lambda\leq\tau$, we have
$\hat{\tau}_{k}=(1-\tau\hat{\gamma}_{k})\mathbf{1}(\tau\hat{\gamma}_{k}<1)$ is
continuous. With a similar argument as before, we also know that
$\hat{\tau}_{k}$ is absolutely continuous. ∎

Now we start to compute
$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}$.
Since

$\displaystyle\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}=\mathbb{E}\big{[}\|\mathbf{y}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}+\frac{2\sigma^{2}}{n}\text{df}(\hat{f}_{\text{stack}})-\sigma^{2},$

we need to compute $\text{df}(\hat{f}_{\text{stack}})$. In the case of
$\lambda>\tau$, recall that in the proof of Theorem 5.2, we know
$\hat{\beta}_{k}=\tau\hat{\gamma}_{k}$ for any $k\leq s_{a}=\tilde{m}$ where
$\tilde{m}=\max\big{\\{}k\in\\{0,1,\dots,M\\}:\hat{\gamma}_{k}<1/\lambda\big{\\}}$
and $\hat{\beta}_{k}=1$ for any $k>\tilde{m}$. Using a similar argument as in
the computation of $\text{df}(\hat{f}_{\text{best}})$, we have

$\displaystyle\text{df}(\hat{f}_{\text{stack}})$
$\displaystyle=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\sum_{k=1}^{M}\sum_{l\in
A_{k}\backslash
A_{k-1}}\text{cov}\left(\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\psi_{l}(\mathbf{x}_{i})(1-\hat{\beta}_{k}),y_{i}\right)$
(48) $\displaystyle=\frac{n}{\sigma^{2}}\sum_{k=1}^{M}\sum_{l\in
A_{k}\backslash
A_{k-1}}\text{cov}\left(\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle,\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle(1-\hat{\beta}_{k})\right)$
$\displaystyle=\sum_{k=1}^{M}\sum_{l\in A_{k}\backslash
A_{k-1}}\mathbb{E}\left[\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle(1-\hat{\beta}_{k})\right]+\text{sdf}(\hat{f}_{\text{stack}})$
$\displaystyle=\sum_{k=1}^{M}\sum_{l\in A_{k}\backslash
A_{k-1}}\mathbb{E}\left[1-\hat{\beta}_{k}+\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}(1-\hat{\beta}_{k})\right]+\text{sdf}(\hat{f}_{\text{stack}}),$

where $\text{sdf}(\hat{f}_{\text{stack}})$ is the search degrees of freedom of
$\hat{f}_{\text{stack}}$. To compute
$\mathbb{E}\left[\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}\hat{\beta}_{k}\right]$,
according to (47) and by the symmetry of the case when
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle<0$, we have

$\displaystyle\mathbb{E}\left[\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}\hat{\beta}_{k}\right]$
(49)
$\displaystyle=~{}\mathbb{E}\left[\Bigg{(}\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}\hat{\beta}_{k}\Bigg{)}\Bigg{(}\sum_{i=1}^{K}\mathbf{1}(t_{i-1}<|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|\leq
t_{i})+\mathbf{1}(t_{i}<|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|)\Bigg{)}\right]$
$\displaystyle=~{}\sum_{i=1}^{K}\mathbb{E}\Bigg{[}\Bigg{(}\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}\frac{\tau\sigma^{2}}{n}\frac{d_{c_{i}}-d_{b_{i}}}{R_{b_{i}}-R_{c_{i}}}\Bigg{)}\mathbf{1}(t_{i-1}<|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|\leq
t_{i})\Bigg{]}$
$\displaystyle\qquad+\mathbb{E}\Bigg{[}\Bigg{(}\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}\frac{\tau\sigma^{2}}{n}\frac{d_{c_{K+1}}-d_{b_{K+1}}}{R_{b_{K+1}}-R_{c_{K+1}}}\Bigg{)}\mathbf{1}(t_{K}<|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|)\Bigg{]}.$

Note that in (48), we know $l\in A_{k}\backslash A_{k-1}$ which means
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle^{2}$ always appears
in the denominator $R_{j-1}-R_{i}$ of
$\hat{\gamma}_{k}=\frac{\sigma^{2}}{n}\min_{k\leq i\leq M}\max_{1\leq j\leq
k}\frac{d_{i}-d_{j-1}}{R_{j-1}-R_{i}}$. Thus, in (49), we know
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle$ should also always
appear in $R_{b_{i}}-R_{c_{i}}$ for any $i=1,2,\dots,K+1$. Recall that
$s_{0}=0,s_{1},\dots,s_{u}=M$ correspond to the change points of
$\hat{\boldsymbol{\gamma}}$. Suppose $a(k)$ is the index for which
$s_{a(k)-1}<k\leq s_{a(k)}$. Here $s_{0},s_{1},\dots,s_{u}$, and $a(k)$ are
all functions of $\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle$.
Then we can continue from (49) and obtain

$\displaystyle\mathbb{E}\left[\frac{\partial}{\partial\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle}\hat{\beta}_{k}\right]$
$\displaystyle=\sum_{i=1}^{K}\mathbb{E}\Bigg{[}-2\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\frac{\tau\sigma^{2}}{n}\frac{d_{c_{i}}-d_{b_{i}}}{(R_{b_{i}}-R_{c_{i}})^{2}}\mathbf{1}(t_{i-1}<|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|\leq
t_{i})\Bigg{]}$ (50)
$\displaystyle\quad+\mathbb{E}\Bigg{[}-2\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\frac{\tau\sigma^{2}}{n}\frac{d_{c_{K+1}}-d_{b_{K+1}}}{(R_{b_{K+1}}-R_{c_{K+1}})^{2}}\mathbf{1}(t_{K}<|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|)\Bigg{]}$
(51)
$\displaystyle=\mathbb{E}\Bigg{[}-2\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\frac{\tau\sigma^{2}}{n}\frac{d_{s_{a(k)}}-d_{s_{a(k)-1}}}{(R_{s_{a(k)-1}}-R_{s_{a(k)}})^{2}}\mathbf{1}(t_{0}<|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|)\Bigg{]}$
(52)
$\displaystyle=\mathbb{E}\Bigg{[}-2\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle\frac{\tau\sigma^{2}}{n}\frac{d_{s_{a(k)}}-d_{s_{a(k)-1}}}{(R_{s_{a(k)-1}}-R_{s_{a(k)}})^{2}}\mathbf{1}(k\leq
s_{a})\Bigg{]}.$ (53)

Continuing from (48), we have

$\displaystyle\text{df}(\hat{f}_{\text{stack}})$ (54)
$\displaystyle=~{}\mathbb{E}\left[\sum_{k=1}^{M}\Delta
d_{k}(1-\hat{\beta}_{k})+2\sum_{k=1}^{s_{a}}\sum_{l\in A_{k}\backslash
A_{k-1}}\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle^{2}\frac{\tau\sigma^{2}}{n}\frac{d_{s_{a(k)}}-d_{s_{a(k)-1}}}{\left(R_{s_{a(k)-1}}-R_{s_{a(k)}}\right)^{2}}\right]+\text{sdf}(\hat{f}_{\text{stack}})$
(55) $\displaystyle=~{}\mathbb{E}\left[\sum_{k=1}^{M}\Delta
d_{k}(1-\hat{\beta}_{k})+2\sum_{l=1}^{a}\sum_{l\in A_{s_{l}}\backslash
A_{s_{l-1}}}\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle^{2}\frac{\tau\sigma^{2}}{n}\frac{d_{s_{l}}-d_{s_{l-1}}}{\left(R_{s_{l-1}}-R_{s_{l}}\right)^{2}}\right]+\text{sdf}(\hat{f}_{\text{stack}})$
(56) $\displaystyle=~{}\mathbb{E}\left[\sum_{k=1}^{M}\Delta
d_{k}(1-\hat{\beta}_{k})+2\sum_{l=1}^{a}\hat{\beta}_{s_{l}}\right]+\text{sdf}(\hat{f}_{\text{stack}}).$
(57)

In the case of $\lambda\leq\tau$, let $\tilde{a}$ satisfy
$s_{\tilde{a}}=\max\big{\\{}k\in\\{0,1,\dots,M\\}:\hat{\gamma}_{k}<1/\tau\big{\\}}$.
Since $\hat{\beta}_{k}$ is continuous in this case, it should hold that
$\text{sdf}(\hat{f}_{\text{stack}})=0$. Thus, with a similar argument as
above, we have

$\displaystyle\text{df}(\hat{f}_{\text{stack}})=\mathbb{E}\left[\sum_{k=1}^{M}\Delta
d_{k}(1-\hat{\beta}_{k})+2\sum_{l=1}^{\tilde{a}}\hat{\beta}_{s_{l}}\right].$
(58)

Next we start to characterize $\text{sdf}(\hat{f}_{\text{best}})$ and
$\text{sdf}(\hat{f}_{\text{stack}})$. Given $1\leq k\leq M$ and $l\in
A_{k}\backslash A_{k-1}$, we know both $\hat{\beta}_{k}$ and
$\tilde{\beta}_{k}$ should have the same discontinuous points with respect to
$\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle$. In addition,
suppose $t\geq 0$ is the discontinuous point then $-t$ is also the
discontinuous point. By part (i) of Lemma A.6, we know $\hat{\gamma}_{k}$ is
nonincreasing with
$|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|$ and thus, both
$\hat{\beta}_{k}$ and $\tilde{\beta}_{k}$ have at most two discontinuous
points. To indicate the dependence with $l$ and $k$, we denote the two
discontinuous points by $t_{l,k}$ and $-t_{l,k}$. Again since
$\hat{\gamma}_{k}$ is nonincreasing with
$|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|$, we know
$\tilde{\beta}_{k}=\mathbf{1}(|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|\leq
t_{l,k})$ and
$\hat{\beta}_{k}=\mathbf{1}(|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|\leq
t_{l,k})+\tau\hat{\gamma}_{k}\mathbf{1}(|\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle|>t_{l,k})$.
Therefore, applying Lemma A.5 to
$h(\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle)=\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle(1-\tilde{\beta}_{k})$,
we have

$\displaystyle\text{sdf}(\hat{f}_{\text{best}})=\frac{\sqrt{n}}{\sigma}\sum_{k=1}^{M}\sum_{l\in
A_{k}\backslash
A_{k-1}}\mathbb{E}\Bigg{[}t_{l,k}\Bigg{(}\phi\Bigg{(}\frac{\sqrt{n}(t_{l,k}-\mu_{l})}{\sigma}\Bigg{)}+\phi\Bigg{(}\frac{\sqrt{n}(-t_{l,k}-\mu_{l})}{\sigma}\Bigg{)}\Bigg{)}\Bigg{]}>0,$

where $\mu_{l}=\left\langle\mathbf{f},\boldsymbol{\psi}_{l}\right\rangle$,
c.f., (Mikkelsen and Hansen, 2018; Tibshirani, 2015; Tibshirani and Rosset,
2019). Also, when $\lambda>\tau$, we take
$h(\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle)=\left\langle\mathbf{y},\boldsymbol{\psi}_{l}\right\rangle(1-\hat{\beta}_{k})$
and obtain

$\displaystyle\text{sdf}(\hat{f}_{\text{stack}})$
$\displaystyle=\frac{\sqrt{n}}{\sigma}\sum_{k=1}^{M}\sum_{l\in A_{k}\backslash
A_{k-1}}\mathbb{E}\Bigg{[}t_{l,k}(1-\tau/\lambda)\Bigg{(}\phi\Bigg{(}\frac{\sqrt{n}(t_{l,k}-\mu_{l})}{\sigma}\Bigg{)}+\phi\Bigg{(}\frac{\sqrt{n}(-t_{l,k}-\mu_{l})}{\sigma}\Bigg{)}\Bigg{)}\Bigg{]}$
$\displaystyle=(1-\tau/\lambda)\text{sdf}(\hat{f}_{\text{best}}).$

When $\lambda\leq\tau$, according to (58), we already know that
$\text{sdf}(\hat{f}_{\text{stack}})=0$. Now we are ready to show
$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}<\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}$.
In the case of $\lambda>\tau$, recall that $0=s_{0},s_{1},\dots,s_{u}=M$
correspond to the change points of $\hat{\boldsymbol{\gamma}}$ and
$\tilde{\beta}_{k}=0$ when $k\leq s_{a}$ and 1 otherwise. Thus, we know

$\displaystyle\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}$
$\displaystyle=\mathbb{E}\left[R_{0}+\sum_{k=1}^{M}\Delta
R_{k}\left(\tilde{\beta}_{k}^{2}-1\right)+\sum_{k=1}^{M}\frac{2\sigma^{2}}{n}\Delta
d_{k}(1-\tilde{\beta}_{k})\right]+\frac{2\sigma^{2}}{n}\text{sdf}(\hat{f}_{\text{best}})-\sigma^{2}$
$\displaystyle=\mathbb{E}\left[R_{0}+\sum_{l=1}^{a}\Delta
R_{s_{l}}\left(\tilde{\beta}_{s_{l}}^{2}-1\right)+\sum_{l=1}^{a}\frac{2\sigma^{2}}{n}\Delta
d_{s_{l}}(1-\tilde{\beta}_{s_{l}})\right]+\frac{2\sigma^{2}}{n}\text{sdf}(\hat{f}_{\text{best}})-\sigma^{2},$

where $\tilde{\beta}_{s_{l}}=0$ for any $1\leq l\leq a$. With a similar
argument, we have

$\displaystyle\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}$
$\displaystyle=~{}\mathbb{E}\left[R_{0}+\sum_{l=1}^{a}\Delta
R_{s_{l}}\left(\hat{\beta}_{s_{l}}^{2}-1\right)+\sum_{l=1}^{a}\frac{2\sigma^{2}}{n}\Delta
d_{s_{l}}(1-\hat{\beta}_{s_{l}})+\frac{4\sigma^{2}}{n}\sum_{l=1}^{a}\hat{\beta}_{s_{l}}\right]+\frac{2\sigma^{2}}{n}\text{sdf}(\hat{f}_{\text{stack}})-\sigma^{2},$

where
$\hat{\beta}_{s_{l}}=\tau\hat{\gamma}_{s_{l}}=\frac{\tau\sigma^{2}}{n}\frac{d_{s_{l}}-d_{s_{l-1}}}{R_{s_{l-1}}-R_{s_{l}}}$
for any $1\leq l\leq a$. Define the quadratic function $h_{l}(z)=\Delta
R_{s_{l}}z^{2}-\frac{2\sigma^{2}}{n}(\Delta d_{s_{l}}-2)z$. It is evident that
if $\Delta d_{s_{l}}\geq\min_{k}\Delta d_{k}\geq 4/(2-\tau)$, then
$h_{l}(\hat{\beta}_{s_{l}})\leq h_{l}(0)$ holds for all $1\leq l\leq a$. Also,
we have proved
$\text{sdf}(\hat{f}_{\text{stack}})=(1-\tau/\lambda)\text{sdf}(\hat{f}_{\text{best}})<\text{sdf}(\hat{f}_{\text{best}})$,
since $\text{sdf}(\hat{f}_{\text{best}})>0$. Thus,
$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}<\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}$
when $\lambda>\tau$. In addition, the risk gap
$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}-\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}$
is equal to

$\frac{\sigma^{4}\tau(2-\tau)}{n^{2}}\mathbb{E}\Bigg{[}\sum_{l=1}^{a}\frac{\Delta
d_{s_{l}}(\Delta d_{s_{l}}-4/(2-\tau))}{\Delta
R_{s_{l}}}\Bigg{]}+\frac{2\tau}{\lambda}\frac{\sigma^{2}}{n}\text{sdf}(\hat{f}_{\text{best}}).$
(59)

We can further lower bound the first term in (59) by

$\displaystyle\frac{\sigma^{4}\tau(2-\tau)}{n^{2}}\mathbb{E}\Bigg{[}\sum_{l=1}^{a}\frac{(\Delta
d_{s_{l}}-4/(2-\tau))^{2}}{\Delta R_{s_{l}}}\Bigg{]}$
$\displaystyle\geq~{}\frac{\sigma^{4}\tau(2-\tau)}{n^{2}}\mathbb{E}\Bigg{[}\frac{(d_{s_{a}}-4a/(2-\tau))^{2}}{R_{0}-R_{s_{a}}}\Bigg{]}$
(60)
$\displaystyle\geq~{}\frac{\sigma^{4}\tau(2-\tau)}{n^{2}}\mathbb{E}\Bigg{[}\frac{(d_{s_{a}}-4s_{a}/(2-\tau))^{2}}{R_{0}-R_{s_{a}}}\Bigg{]}$
$\displaystyle\geq~{}\frac{\sigma^{4}\tau(2-\tau)}{n^{2}}\mathbb{E}\Bigg{[}\min_{1\leq
k\leq M}\frac{(d_{k}-4k/(2-\tau))^{2}}{R_{0}-R_{k}}\Bigg{]},$

where we apply the Cauchy-Schwarz inequality to obtain the first inequality.
This establishes (10) when $\lambda>\tau$. In the case of $\lambda\leq\tau$,
by (58) we have

$\displaystyle\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}$
(61) $\displaystyle=~{}\mathbb{E}\left[R_{0}+\sum_{l=1}^{\tilde{a}}\Delta
R_{s_{l}}\left(\hat{\beta}_{s_{l}}^{2}-1\right)+\sum_{l=1}^{\tilde{a}}\frac{2\sigma^{2}}{n}\Delta
d_{s_{l}}(1-\hat{\beta}_{s_{l}})+\frac{4\sigma^{2}}{n}\sum_{l=1}^{\tilde{a}}\hat{\beta}_{s_{l}}\right]-\sigma^{2},$

where we recall $\tilde{a}$ satisfies
$s_{\tilde{a}}=\max\big{\\{}k\in\\{0,1,\dots,M\\}:\hat{\gamma}_{k}<1/\tau\big{\\}}$.
Again, we know $h_{l}(\hat{\beta}_{s_{l}})\leq h_{l}(0)$ for any $1\leq
l\leq\tilde{a}$ when $\Delta d_{s_{l}}\geq 4/(2-\tau)$. Note that
$\tilde{a}\leq a$ and for any $\tilde{a}<l\leq a$, according to the definition
of $\tilde{a}$, we have $\hat{\gamma}_{s_{l}}=\frac{\sigma^{2}}{n}\frac{\Delta
d_{s_{l}}}{\Delta R_{s_{l}}}\geq 1/\tau>1/2$. Thus for any $\tilde{a}<l\leq
a$,

$\Delta
R_{s_{l}}\left(\tilde{\beta}_{s_{l}}^{2}-1\right)+\frac{2\sigma^{2}}{n}\Delta
d_{s_{l}}(1-\tilde{\beta}_{s_{l}})=-\Delta
R_{s_{l}}+\frac{2\sigma^{2}}{n}\Delta d_{s_{l}}>0.$

Therefore, we obtain
$\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}<\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{best}}\|^{2}\big{]}$,
and with a similar argument as (60), we have can lower bound the risk gap by

$\displaystyle\frac{\sigma^{4}\tau(2-\tau)}{n^{2}}\mathbb{E}\Bigg{[}\sum_{l=1}^{a}\frac{\Delta
d_{s_{l}}(\Delta d_{s_{l}}-4/(2-\tau))}{\Delta
R_{s_{l}}}\Bigg{]}+\frac{2\sigma^{2}}{n}\text{sdf}(\hat{f}_{\text{best}})$
$\displaystyle>~{}\frac{\sigma^{4}\tau(2-\tau)}{n^{2}}\mathbb{E}\Bigg{[}\min_{1\leq
k\leq
M}\frac{(d_{k}-4k/(2-\tau))^{2}}{R_{0}-R_{k}}\Bigg{]}+\frac{2\sigma^{2}}{n}\text{sdf}(\hat{f}_{\text{best}}),$

which establishes (10) when $\lambda\leq\tau$. ∎

###### Proof of Theorem 4.2.

With equality (61) and the assumption that $\tau=\lambda=1$, we know

$\displaystyle\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}$
(62) $\displaystyle=~{}\mathbb{E}\left[R_{0}+\sum_{l=1}^{\tilde{a}}\Delta
R_{s_{l}}\left(\hat{\beta}_{s_{l}}^{2}-1\right)+\sum_{l=1}^{\tilde{a}}\frac{2\sigma^{2}}{n}\Delta
d_{s_{l}}(1-\hat{\beta}_{s_{l}})+\frac{4\sigma^{2}}{n}\sum_{l=1}^{\tilde{a}}\hat{\beta}_{s_{l}}\right]-\sigma^{2}$
(63) $\displaystyle=~{}\mathbb{E}\left[R_{0}+\sum_{k=1}^{M}\Delta
R_{k}\left(\hat{\beta}_{k}^{2}-1\right)+\sum_{k=1}^{M}\frac{2\sigma^{2}}{n}\Delta
d_{k}(1-\hat{\beta}_{k})+\frac{4\sigma^{2}}{n}\sum_{k=1}^{M}\hat{\beta}_{k}\mathbf{1}(\hat{\beta}_{k}\neq\hat{\beta}_{k+1})\right]-\sigma^{2}$
(64) $\displaystyle\leq~{}\mathbb{E}\left[R_{0}+\sum_{k=1}^{M}\Delta
R_{k}\left(\beta_{k}^{2}-1\right)+\sum_{k=1}^{M}\frac{2\sigma^{2}}{n}\Delta
d_{k}(1-\beta_{k})\right]+\mathbb{E}\left[\frac{4\sigma^{2}}{n}\sum_{k=1}^{M}\hat{\beta}_{k}\mathbf{1}(\hat{\beta}_{k}\neq\hat{\beta}_{k+1})\right]-\sigma^{2},$
(65)

where the last inequality holds for any deterministic or random
$\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{M}\leq 1$, since
$\hat{\boldsymbol{\beta}}$ is the solution of program (12) and
$\tau=\lambda=1$. For the second equality, we use the fact that
$\hat{\beta}_{k}=1$ when $k>s_{a}$. Recall that
$\alpha_{k}=\beta_{k+1}-\beta_{k}\geq 0$. Note that with deterministic
$\\{\beta_{k}\\}$ (and $\\{\alpha_{k}\\}$),

$\displaystyle\mathbb{E}\Bigg{[}\Bigg{\|}\mathbf{f}-\sum_{k=1}^{M}\alpha_{k}\boldsymbol{\hat{\mu}}_{k}\Bigg{\|}^{2}\Bigg{]}$
$\displaystyle=\mathbb{E}\Bigg{[}\Bigg{\|}\mathbf{y}-\sum_{k=1}^{M}\alpha_{k}\boldsymbol{\hat{\mu}}_{k}\Bigg{\|}^{2}\Bigg{]}-\sigma^{2}+\sum_{k=1}^{M}\frac{2\sigma^{2}}{n}d_{k}\alpha_{k}$
$\displaystyle=\mathbb{E}\left[R_{0}+\sum_{k=1}^{M}\Delta
R_{k}\left(\beta_{k}^{2}-1\right)+\sum_{k=1}^{M}\frac{2\sigma^{2}}{n}\Delta
d_{k}(1-\beta_{k})\right]-\sigma^{2}.$

Since $\boldsymbol{\alpha}^{\star}$ is deterministic, we thus have

$\displaystyle\mathbb{E}\big{[}\|\mathbf{f}-\mathbf{\hat{f}}_{\text{stack}}\|^{2}\big{]}\leq\mathbb{E}\Bigg{[}\Bigg{\|}\mathbf{f}-\sum_{k=1}^{M}\alpha^{\star}_{k}\boldsymbol{\hat{\mu}}_{k}\Bigg{\|}^{2}\Bigg{]}+\mathbb{E}\left[\frac{4\sigma^{2}}{n}\sum_{k=1}^{M}\hat{\beta}_{k}\mathbf{1}(\hat{\beta}_{k}\neq\hat{\beta}_{k+1})\right].$

Theorem 4.2 follows from the fact that
$\sum_{k=1}^{M}\hat{\beta}_{k}\mathbf{1}(\hat{\beta}_{k}\neq\hat{\beta}_{k+1})=\sum_{k=1}^{M}\Big{(}1-\sum_{m=k}^{M}\hat{\alpha}_{m}\Big{)}\mathbf{1}(\hat{\alpha}_{k}\neq
0)$, and the inequality

$\mathbb{E}\left[\sum_{k=1}^{M}\Bigg{(}1-\sum_{m=k}^{M}\hat{\alpha}_{m}\Bigg{)}\mathbf{1}(\hat{\alpha}_{k}\neq
0)\right]\leq\mathbb{E}\left[\sum_{k=1}^{M}\mathbf{1}(\hat{\alpha}_{k}\neq
0)\right]=\mathbb{E}\left[\|\hat{\boldsymbol{\alpha}}\|_{\ell_{0}}\right].\qed$

###### Proof of Theorem 7.1.

By setting $\lambda=\tau=1$ in (8), and noting that
$\|\boldsymbol{\alpha}\|_{\ell_{0}}=\sum_{k=1}^{M}\mathbf{1}(\hat{\beta}_{k}\neq\hat{\beta}_{k+1})$,
we can establish that the solution to (18) is equivalent to the solution of
the following program:

$\displaystyle\text{minimize}\quad\sum_{k=1}^{M}w_{k}(z_{k}-\beta_{k})^{2}+\frac{4\sigma^{2}}{n}\sum_{k=1}^{M}\mathbf{1}(\beta_{k}\neq\beta_{k+1})$
(66) $\displaystyle\text{subject
to}\quad\beta_{1}\leq\beta_{2}\leq\cdots\leq\beta_{M}\leq 1,$

where $w_{k}=\Delta R_{k}>0$ and $z_{k}=(\sigma^{2}/n)(\Delta d_{k}/\Delta
R_{k})>0$. Recall that we assume
$r_{q}=\max\big{\\{}k\in\\{0,1,\dots,M\\}:\hat{\beta}_{k}<1\big{\\}}$.
According to Lemma A.2, to show
$\text{dim}(\hat{f}_{\text{stack}})\geq\text{dim}(\hat{f}_{\text{best}})$, it
suffices to show $r_{q}\geq\tilde{m}$. If $r_{q}=M$, then
$M=r_{q}\geq\tilde{m}$ is already satisfied. If $r_{q}<M$, we first show for
any $r_{q}<j\leq M$,

$\displaystyle\frac{\hat{\beta}_{r_{q}}+1}{2}\leq\frac{\sigma^{2}}{n}\frac{d_{j}-d_{r_{q}}}{R_{r_{q}}-R_{j}}.$
(67)

Consider another sequence $\boldsymbol{\tilde{\beta}}$ where
$\tilde{\beta}_{k}=\hat{\beta}_{k}\mathbf{1}(1\leq k\leq
r_{q})+\hat{\beta}_{r_{q}}\mathbf{1}(r_{q}<k\leq j)+\mathbf{1}(j<k\leq M)$. By
the optimality of $\hat{\boldsymbol{\beta}}$, we have

$\displaystyle\sum_{k=1}^{M}\Delta
R_{k}\Bigg{(}-\hat{\beta}_{k}+\frac{\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}\Bigg{)}^{2}+\frac{4\sigma^{2}}{n}\sum_{k=1}^{M}\hat{\beta}_{k}\mathbf{1}(\hat{\beta}_{k}\neq\hat{\beta}_{k+1})$
$\displaystyle\leq~{}\sum_{k=1}^{M}\Delta
R_{k}\Bigg{(}-\tilde{\beta}_{k}+\frac{\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}\Bigg{)}^{2}+\frac{4\sigma^{2}}{n}\sum_{k=1}^{M}\tilde{\beta}_{k}\mathbf{1}(\tilde{\beta}_{k}\neq\tilde{\beta}_{k+1}).$

This implies

$\displaystyle\sum_{k=r_{q}+1}^{j}\Delta
R_{k}\Bigg{(}1-\frac{\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}\Bigg{)}^{2}+\sum_{k=j+1}^{M}\Delta
R_{k}\Bigg{(}1-\frac{\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}\Bigg{)}^{2}$ $\displaystyle\leq~{}\sum_{k=r_{q}+1}^{j}\Delta
R_{k}\Bigg{(}\hat{\beta}_{r_{q}}-\frac{\sigma^{2}}{n}\frac{\Delta
d_{k}}{\Delta R_{k}}\Bigg{)}^{2}+\sum_{k=j+1}^{M}\Delta
R_{k}\Bigg{(}1-\frac{\sigma^{2}}{n}\frac{\Delta d_{k}}{\Delta
R_{k}}\Bigg{)}^{2}.$

After simplification, we have

$\displaystyle\left(R_{r_{q}}-R_{j}\right)\left(1-\frac{\sigma^{2}}{n}\frac{d_{j}-d_{r_{q}}}{R_{r_{q}}-R_{j}}\right)^{2}$
$\displaystyle\leq\left(R_{r_{q}}-R_{j}\right)\left(\hat{\beta}_{r_{q}}-\frac{\sigma^{2}}{n}\frac{d_{j}-d_{r_{q}}}{R_{r_{q}}-R_{j}}\right)^{2}$
$\displaystyle\Leftrightarrow\left|1-\frac{\sigma^{2}}{n}\frac{d_{j}-d_{r_{q}}}{R_{r_{q}}-R_{j}}\right|$
$\displaystyle\leq\left|\hat{\beta}_{r_{q}}-\frac{\sigma^{2}}{n}\frac{d_{j}-d_{r_{q}}}{R_{r_{q}}-R_{j}}\right|,$

which implies (67) since $\hat{\beta}_{r_{q}}<1$. Suppose $r_{q}<\tilde{m}$.
Then

$\displaystyle\frac{1}{2}\geq\tilde{\gamma}_{\tilde{m}}=\max_{1\leq
i\leq\tilde{m}}\frac{\sigma^{2}}{n}\frac{d_{\tilde{m}}-d_{i-1}}{R_{i-1}-R_{\tilde{m}}}\geq\frac{\sigma^{2}}{n}\frac{d_{\tilde{m}}-d_{r_{q}}}{R_{r_{q}}-R_{\tilde{m}}}.$

Noting that $\hat{\beta}_{r_{q}}>0$, we obtain a contradiction when taking
$j=\tilde{m}$ in (67). Therefore $r_{q}\geq\tilde{m}.$

To demonstrate $\sum_{k=1}^{M}\hat{\alpha}_{k}<1$, it suffices to show
$\hat{\beta}_{1}>0$. Using an argument similar to (Gao et al., 2020, Lemma
5.1) to characterize the form of the solution to (weighted) reduced isotonic
regression, we can show that there exists an index $1\leq k\leq M$ for which
$\hat{\beta}_{1}=\frac{\sigma^{2}}{n}\frac{d_{k}}{R_{0}-R_{k}}>0$. ∎

###### Proof of the solution of program (LABEL:loss:q-agg).

The program proposed in (Bellec, 2018; Bellec and Yang, 2020) is

$\displaystyle\text{minimize}\quad
R(\boldsymbol{\alpha})+\frac{2\sigma^{2}}{n}\sum_{k=1}^{M}\alpha_{k}d_{k}+\eta\sum_{m=1}^{M}\alpha_{m}\frac{1}{n}\sum_{i=1}^{n}\Bigg{(}\hat{\mu}_{m}(\mathbf{x}_{i})-\sum_{k=1}^{M}\alpha_{k}\hat{\mu}_{k}(\mathbf{x}_{i})\Bigg{)}^{2}$
(68) $\displaystyle\text{subject to}\quad\alpha_{k}\geq
0,\;\;k=1,2,\dots,M,\;\;\text{and}\;\;\sum_{k=1}^{M}\alpha_{k}=1,$

with $\eta=1/2$. Thanks to a bias-variance decomposition, when
$\sum_{k=1}^{M}\alpha_{k}=1$,

$\displaystyle
R(\boldsymbol{\alpha})+\frac{2\sigma^{2}}{n}\sum_{k=1}^{M}\alpha_{k}d_{k}+\eta\sum_{m=1}^{M}\alpha_{m}\frac{1}{n}\sum_{i=1}^{n}\Bigg{(}\hat{\mu}_{m}(\mathbf{x}_{i})-\sum_{k=1}^{M}\alpha_{k}\hat{\mu}_{k}(\mathbf{x}_{i})\Bigg{)}^{2}$
(69)
$\displaystyle=~{}(1-\eta)R(\boldsymbol{\alpha})+\frac{2\sigma^{2}}{n}\sum_{k=1}^{M}\alpha_{k}d_{k}+\eta\sum_{k=1}^{M}\alpha_{k}R_{k}.$
(70)

Thus, program (LABEL:prog:bellec_yang) corresponds to Lemma A.7, which also
demonstrates its equivalence to program (LABEL:loss:q-agg). ∎

######  Lemma A.7.

Consider the program

$\displaystyle\text{minimize}\quad
R(\boldsymbol{\alpha})+\sum_{k=1}^{M}\alpha_{k}\Bigg{(}\frac{\frac{2\sigma^{2}}{n}d_{k}+\eta
R_{k}}{1-\eta}\Bigg{)}$ (71) $\displaystyle\text{subject
to}\quad\alpha_{k}\geq 0,\;\;k=1,2,\dots,M.$

If $0<\eta<1$ and $\sum_{k=1}^{M}\alpha_{k}=1$, then the solution to
(LABEL:prog:general) is

$\check{\alpha}_{k}=\phi\Bigg{(}\frac{1}{2}\Bigg{(}1+\frac{1-2\check{\gamma}_{k}}{1-\eta}\Bigg{)}\Bigg{)}-\phi\Bigg{(}\frac{1}{2}\Bigg{(}1+\frac{1-2\check{\gamma}_{k+1}}{1-\eta}\Bigg{)}\Bigg{)},$
(72)

where $\phi(z)=\min\\{1,\max\\{0,z\\}\\}$ is the clip function for $z$ at 0
and 1, and

$\check{\gamma}_{1}=0,\quad\check{\gamma}_{k}=\frac{\sigma^{2}}{n}\min_{k\leq
i\leq M}\max_{1\leq j<k}\frac{d_{i}-d_{j}}{R_{j}-R_{i}},\quad k=2,3,\dots,M.$
(73)

###### Proof.

Recall that $\alpha_{k}=\beta_{k+1}-\beta_{k}$ and $c_{k}=1-\beta_{k}$. Using
summation by parts and Lemma A.1, we can establish the equivalence between the
solution of (LABEL:prog:general) and the solution of the following program:

$\displaystyle\text{minimize}\quad\frac{\eta}{1-\eta}R_{0}c_{1}+\sum_{k=1}^{M}\Delta
R_{k}\Bigg{(}c_{k}-\frac{1}{2}\Bigg{(}\frac{2-\eta}{1-\eta}-\frac{2}{1-\eta}\frac{\sigma^{2}}{n}\frac{\Delta
d_{k}}{\Delta R_{k}}\Bigg{)}\Bigg{)}^{2}$ (74) $\displaystyle\text{subject
to}\quad c_{1}\geq c_{2}\geq\cdots\geq c_{M}\geq 0.$

When $0<\eta<1$ and $\sum_{k=1}^{M}\alpha_{k}=1$, we note that
$c_{1}=\sum_{k=1}^{M}\alpha_{k}=1$ and $\beta_{1}=1-c_{1}=0$. Thus,
(LABEL:prog:simplify) can be further reduced to (LABEL:loss:q-agg). By Lemma
A.3, we know the solution is $\check{\beta}_{1}=0$ and
$\check{\beta}_{k}=1-\phi\Big{(}\frac{1}{2}\Big{(}1+\frac{1-2\check{\gamma}_{k}}{1-\eta}\Big{)}\Big{)}=1-\phi\Big{(}\frac{1}{1-\eta}\Big{(}1-\eta/2-\check{\gamma}_{k}\Big{)}\Big{)}$
for $k=2,3,\dots,M$. Consequently,

$\check{\alpha}_{k}=\phi\Bigg{(}\frac{1}{2}\Bigg{(}1+\frac{1-2\check{\gamma}_{k}}{1-\eta}\Bigg{)}\Bigg{)}-\phi\Bigg{(}\frac{1}{2}\Bigg{(}1+\frac{1-2\check{\gamma}_{k+1}}{1-\eta}\Bigg{)}\Bigg{)},\quad
k=1,2,\dots,M.$ (75)

∎

## References

  * Agarwal et al. [2022] Abhineet Agarwal, Yan Shuo Tan, Omer Ronen, Chandan Singh, and Bin Yu.  Hierarchical shrinkage: Improving the accuracy and interpretability of tree-based models.  In _International Conference on Machine Learning_ , pages 111–135. PMLR, 2022.
  * Athey and Imbens [2016] Susan Athey and Guido Imbens.  Recursive partitioning for heterogeneous causal effects.  _Proceedings of the National Academy of Sciences_ , 113(27):7353–7360, 2016.  doi: 10.1073/pnas.1510489113.  URL https://www.pnas.org/doi/abs/10.1073/pnas.1510489113.
  * Barlow and Brunk [1972] R. E. Barlow and H. D. Brunk.  The isotonic regression problem and its dual.  _Journal of the American Statistical Association_ , 67(337):140–147, 1972.  doi: 10.1080/01621459.1972.10481216.  URL https://www.tandfonline.com/doi/abs/10.1080/01621459.1972.10481216.
  * Barlow et al. [1972] Richard E Barlow, David J Bartholomew, J Malcolm Bremner, and Hugh D Brunk.  Statistical inference under order restrictions: The theory and application of isotonic regression.  _Wiley_ , 47(3):375–382, 1972.
  * Bellec and Yang [2020] Pierre Bellec and Dana Yang.  The cost-free nature of optimally tuning tikhonov regularizers and other ordered smoothers.  In Hal Daumé III and Aarti Singh, editors, _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pages 746–755. PMLR, 13–18 Jul 2020.  URL https://proceedings.mlr.press/v119/bellec20a.html.
  * Bellec [2018] Pierre C. Bellec.  Optimal bounds for aggregation of affine estimators.  _The Annals of Statistics_ , 46(1):30 – 59, 2018\.  doi: 10.1214/17-AOS1540.  URL https://doi.org/10.1214/17-AOS1540.
  * Bellman [1961] Richard Bellman.  On the approximation of curves by line segments using dynamic programming.  _Commun. ACM_ , 4(6):284, jun 1961.  ISSN 0001-0782.  doi: 10.1145/366573.366611.  URL https://doi.org/10.1145/366573.366611.
  * Breiman [1996] Leo Breiman.  Stacked regressions.  _Machine Learning_ , 24(1):49–64, 1996.  doi: 10.1007/BF00117832.  URL https://doi.org/10.1007/BF00117832.
  * Bühlmann and Yu [2002] Peter Bühlmann and Bin Yu.  Analyzing bagging.  _The Annals of Statistics_ , 30(4):927 – 961, 2002.  doi: 10.1214/aos/1031689014.  URL https://doi.org/10.1214/aos/1031689014.
  * 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.
  * Cavalier et al. [2002] L. Cavalier, G. K. Golubev, D. Picard, and A. B. Tsybakov.  Oracle inequalities for inverse problems.  _The Annals of Statistics_ , 30(3):843 – 874, 2002.  doi: 10.1214/aos/1028674843.  URL https://doi.org/10.1214/aos/1028674843.
  * Dai et al. [2012] Dong Dai, Philippe Rigollet, and Tong Zhang.  Deviation optimal learning using greedy $Q$-aggregation.  _The Annals of Statistics_ , 40(3):1878 – 1905, 2012.  doi: 10.1214/12-AOS1025.  URL https://doi.org/10.1214/12-AOS1025.
  * Gao et al. [2020] Chao Gao, Fang Han, and Cun-Hui Zhang.  On estimation of isotonic piecewise constant signals.  _The Annals of Statistics_ , 48(2):629 – 654, 2020.  doi: 10.1214/18-AOS1792.  URL https://doi.org/10.1214/18-AOS1792.
  * Haiminen et al. [2008] Niina Haiminen, Aristides Gionis, and Kari Laasonen.  Algorithms for unimodal segmentation with applications to unimodality detection.  _Knowledge and Information Systems_ , 14(1):39–57, 2008.  doi: 10.1007/s10115-006-0053-3.  URL https://doi.org/10.1007/s10115-006-0053-3.
  * Hastie et al. [2009] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  _The elements of statistical learning_.  Springer Series in Statistics. Springer-Verlag, New York, 2009.
  * James and Stein [1961] W James and Charles Stein.  Estimation with quadratic loss.  In _Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics_ , volume 4, pages 361–380. University of California Press, 1961.
  * Kneip [1994] Alois Kneip.  Ordered Linear Smoothers.  _The Annals of Statistics_ , 22(2):835 – 866, 1994.  doi: 10.1214/aos/1176325498.  URL https://doi.org/10.1214/aos/1176325498.
  * Koren [2009] Yehuda Koren.  The bellkor solution to the netflix grand prize.  _Netflix prize documentation_ , 81(2009):1–10, 2009.
  * Leung and Barron [2006] G. Leung and A.R. Barron.  Information theory and mixing least-squares regressions.  _IEEE Transactions on Information Theory_ , 52(8):3396–3410, 2006.  doi: 10.1109/TIT.2006.878172.
  * Meyer and Woodroofe [2000] Mary Meyer and Michael Woodroofe.  On the degrees of freedom in shape-restricted regression.  _The Annals of Statistics_ , 28(4):1083 – 1104, 2000.  doi: 10.1214/aos/1015956708.  URL https://doi.org/10.1214/aos/1015956708.
  * Mikkelsen and Hansen [2018] Frederik Riis Mikkelsen and Niels Richard Hansen.  Degrees of freedom for piecewise Lipschitz estimators.  _Annales de l’Institut Henri Poincaré, Probabilités et Statistiques_ , 54(2):819 – 841, 2018.  doi: 10.1214/17-AIHP822.  URL https://doi.org/10.1214/17-AIHP822.
  * Rigollet [2012] Philippe Rigollet.  Kullback–Leibler aggregation and misspecified generalized linear models.  _The Annals of Statistics_ , 40(2):639 – 665, 2012.  doi: 10.1214/11-AOS961.  URL https://doi.org/10.1214/11-AOS961.
  * Shao [1997] Jun Shao.  An asymptotic theory for linear model selection.  _Statistica Sinica_ , 7(2):221–242, 1997.  ISSN 10170405, 19968507.  URL http://www.jstor.org/stable/24306073.
  * Stein [1981] Charles M. Stein.  Estimation of the Mean of a Multivariate Normal Distribution.  _The Annals of Statistics_ , 9(6):1135 – 1151, 1981.  doi: 10.1214/aos/1176345632.  URL https://doi.org/10.1214/aos/1176345632.
  * Stone [1977] M. Stone.  An asymptotic equivalence of choice of model by cross-validation and akaike’s criterion.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):44–47, 1977.  doi: https://doi.org/10.1111/j.2517-6161.1977.tb01603.x.  URL https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/j.2517-6161.1977.tb01603.x.
  * Tibshirani [2015] Ryan J. Tibshirani.  Degrees of freedom and model search.  _Statistica Sinica_ , 25(3):1265–1296, 2015.  ISSN 10170405, 19968507.  URL http://www.jstor.org/stable/24721231.
  * Tibshirani and Rosset [2019] Ryan J. Tibshirani and Saharon Rosset.  Excess optimism: How biased is the apparent error of an estimator tuned by sure?  _Journal of the American Statistical Association_ , 114(526):697–712, 2019.  doi: 10.1080/01621459.2018.1429276.  URL https://doi.org/10.1080/01621459.2018.1429276.
  * Wolpert [1992] David H. Wolpert.  Stacked generalization.  _Neural Networks_ , 5(2):241–259, 1992.  ISSN 0893-6080.  doi: https://doi.org/10.1016/S0893-6080(05)80023-1.  URL https://www.sciencedirect.com/science/article/pii/S0893608005800231.