The Implicit Regularization of Ordinary Least Squares Ensembles

Daniel LeJeune Hamid Javadi Richard G. Baraniuk Rice University Rice
University Rice University

###### Abstract

Ensemble methods that average over a collection of independent predictors that
are each limited to a subsampling of both the examples and features of the
training data command a significant presence in machine learning, such as the
ever-popular random forest, yet the nature of the subsampling effect,
particularly of the features, is not well understood. We study the case of an
ensemble of linear predictors, where each individual predictor is fit using
ordinary least squares on a random submatrix of the data matrix. We show that,
under standard Gaussianity assumptions, when the number of features selected
for each predictor is optimally tuned, the asymptotic risk of a large ensemble
is equal to the asymptotic ridge regression risk, which is known to be optimal
among linear predictors in this setting. In addition to eliciting this
implicit regularization that results from subsampling, we also connect this
ensemble to the dropout technique used in training deep (neural) networks,
another strategy that has been shown to have a ridge-like regularizing effect.

##  1 INTRODUCTION

_Ensemble methods_ (Breiman, 1996; Amit and Geman, 1997; Josse and Wager,
2016) are an oft-used strategy employed successfully in a broad range of
problems in machine learning and statistics, in which one combines a number of
_weak predictors_ together to obtain one powerful predictor. This is
accomplished by giving each weak learner a different _view_ of the training
data. Various strategies for changing this training data view exist, among
which many are simple sampling-based techniques in which each predictor is
(independently) given access to a subsampling of the rows (examples) and
columns (features) of the training data matrix, such as _bagging_ (Breiman,
1996; Bühlmann and Yu, 2002). Another noteworthy technique is _boosting_
(Freund and Schapire, 1997; Breiman, 1998), in which the training data
examples are reweighted adaptively according to how badly they have been
misclassified while building the ensemble. In this work, we consider the
former class of techniques—those that train each weak predictor using an
independent subsampling of the training data.

Figure 1:  Example (rows) and feature (columns) subsampling of the training
data used in the ordinary least squares fit for one member of the ensemble.
The $i$-th member of the ensemble is only allowed to predict using its subset
of the features (green). It must learn its parameters
$\widehat{{\bm{\beta}}}{}^{(i)}$ by performing ordinary least squares using
the subsampled examples of $\mathbf{y}$ (red) and the subsampled examples
(rows) and features (columns) of the data matrix $\mathbf{X}$ (blue,
crosshatched).

Ensemble methods based on independent example and feature subsampling are
attractive for two reasons. First, they are computationally appealing in that
they are massively parallelizable, and since each member of the ensemble uses
only part of the data, they are able to overcome memory limitations faced by
other methods (Louppe and Geurts, 2012). Second, ensemble methods are known to
achieve lower risk due to the fact that combining several different predictors
reduces variance (Bühlmann and Yu, 2002; Wager et al., 2014; Scornet et al.,
2015), and empirically they have been found to perform very well. _Random
forests_ (Breiman, 2001; Athey et al., 2019; Friedberg et al., 2018), for
example, ensemble methods that combine example and feature subsampling with
decision trees by choosing the most useful feature from a random subset of the
features at each branch of the tree, remain among the best-performing off-the-
shelf machine learning methods available (Cutler and Zhao, 2001; Fernández-
Delgado et al., 2014; Wyner et al., 2017).

Let $\mathbf{X}\in\mathbb{R}^{n\times p}$ be the training data matrix
consisting of $n$ examples of data points each having $p$ features. While
there exist theoretical results on the benefits of _example (row) subsampling_
(Bühlmann and Yu, 2002), the exact nature of the effect of _feature (column)
subsampling_ on ensemble performance remains poorly understood. In this paper,
we study the prototypical form of this problem in the context of linear
regression. That is, given the data matrix $\mathbf{X}$ and target variables
$\mathbf{y}\in\mathbb{R}^{n}$, we study the ensemble
$\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}=\frac{1}{k}\sum_{i=1}^{k}\widehat{{\bm{\beta}}}{}^{(i)}$,
where each $\widehat{{\bm{\beta}}}{}^{(i)}$ is learned using ordinary least
squares on an independent random subsampling of both the examples and features
of the training data. This subsampling is illustrated in Figure 1. We show
that under such a scheme, the resulting predictor of this ensemble performs as
well as the _ridge regression_ (Hoerl and Kennard, 1970; Friedman et al.,
2001) predictor fit using the entire training data, which is known to be the
optimal linear predictor under the data assumptions that we consider. Further,
the asymptotic risk of the ensemble depends _only on the amount of feature
subsampling_ and not on the amount of example subsampling, provided each
individual ordinary least squares problem is underdetermined. Our main result
in Theorem 3.6 can be summarized as follows:

######  Theorem 3.6 (informal statement).

When the features and underlying model weights both follow i.i.d. Gaussian
distributions, the optimal asymptotic risk for an ensemble of ordinary least
squares predictors is equal to the optimal asymptotic ridge regression risk.

We can interpret this result as an example of _implicit regularization_
(Mahoney, 2012; Neyshabur et al., 2014; Gunasekar et al., 2017; Arora et al.,
2019). That is, while the individual ordinary least squares subproblems are
completely unregularized, the ensemble behaves as if it had been regularized
using a ridge regression penalty. Recently, there has been much interest in
investigating the implicit regularization effects of commonly used heuristic
methods, particularly in cases where they enable the training of highly
_overparameterized_ models that generalize well to test data despite having
the capacity to overfit the training data (Zhang et al., 2017; Belkin et al.,
2018). Examples of heuristic techniques that have been shown to have implicit
regularization effects include stochastic gradient descent (Hardt et al.,
2016) and _dropout_ (Srivastava et al., 2014). Incidentally, we show a strong
connection between the ensemble of ordinary least squares predictors and
dropout, which is known to have a ridge-like regularizing effect (Wager et
al., 2013), and we make this link via stochastic gradient descent.

##### Contributions

We summarize our contributions as follows: [C1] We prove that when the amount
of feature subsampling is optimized to minimize risk, an ensemble of ordinary
least squares predictors achieves the same risk as the optimal ridge
regression predictor asymptotically as $n,p\to\infty$ (see Section 3). [C2] We
demonstrate the converge of the ensemble risk to the optimal ridge regression
risk via simulation (see Section 4.1). [C3] We reveal a connection between the
ordinary least squares ensemble and the popular _dropout_ technique used in
deep (neural) network training (see Section 4.3) and from the insight gained
from this connection develop a recipe for mitigating excess risk under
suboptimal feature subsampling via simple output scaling (see Section 4.4).

##  2 ENSEMBLES OF ORDINARY LEAST SQUARES PREDICTORS

We consider the familiar setting of _linear regression_ , where there exists a
linear relationship between the target variable $y\in\mathbb{R}$ and the
feature variables $\mathbf{x}\in\mathbb{R}^{p}$—i.e.,
$y=\langle\mathbf{x},{\bm{\beta}}\rangle$, where
${\bm{\beta}}\in\mathbb{R}^{p}$ is the model parameter vector. The goal of a
machine learning algorithm is to estimate these parameters given $n$ i.i.d.
noisy samples $\left\\{\mathbf{x}^{(i)},y^{(i)}\right\\}_{i=1}^{n}$. The noise
relationship is given by

$\displaystyle\mathbf{y}=\mathbf{X}{\bm{\beta}}+\sigma\mathbf{z},$ (1)

where $[\mathbf{X}]_{ij}=[\mathbf{x}^{(i)}]_{j}$, $[\mathbf{y}]_{i}=y^{(i)}$,
and $[\mathbf{z}]_{i}=z^{(i)}$, where $z^{(i)}$ are i.i.d. zero-mean random
variables with unit variance independent of $\mathbf{X}$. We assume a Gaussian
$\mathcal{N}({\bm{0}},{\bm{\Sigma}})$ distribution on $\mathbf{x}$, and for
the results in this paper, we assume ${\bm{\Sigma}}=\mathbf{I}_{p}$.

Our ensemble consists of $k$ linear predictors each fit using ordinary least
squares on a submatrix of $\mathbf{X}$, and the resulting prediction is the
average of the outputs. Equivalently, our ensemble is defined by its estimate
of the parameters

$\displaystyle\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}\triangleq\frac{1}{k}\sum_{i=1}^{k}\widehat{{\bm{\beta}}}{}^{(i)},$
(2)

where $\widehat{{\bm{\beta}}}{}^{(i)}$ is the parameter estimate of the $i$-th
member of the ensemble. To characterize the estimates
$\widehat{{\bm{\beta}}}{}^{(i)}$, we first introduce some notation. Let the
_selection matrix_ $\mathbf{S}$ corresponding to a subset of indices
$S\subseteq[p]$, where $[p]=\left\\{1,\ldots,p\right\\}$, denote the the
$p\times|S|$ matrix obtained by selecting from $\mathbf{I}_{p}$ the columns
corresponding to the indices in $S$, where $\mathbf{I}_{p}$ denotes the
$p\times p$ identity matrix. With this definition of selection matrices, for
$S\subseteq[p]$ and $T\subseteq[n]$, we have that
$\mathbf{T}^{\top}\mathbf{X}\mathbf{S}$ is the matrix of size $|T|\times|S|$
obtained from $\mathbf{X}$ by selecting (subsampling) the rows and columns
indicated by sets $T$ and $S$. Returning to the ensemble, let
$\mathcal{S}\triangleq(S_{i})_{i=1}^{k}$ and
$\mathcal{T}\triangleq(T_{i})_{i=1}^{k}$ denote the collection of _feature
subsets_ and _example subsets_ , respectively, where each $S_{i}\subseteq[p]$
and each $T_{i}\subseteq[n]$. Then, assuming $|S_{i}|<|T_{i}|$, for each
member of the ensemble we let

$\displaystyle\widehat{{\bm{\beta}}}{}^{(i)}_{S_{i}}$
$\displaystyle=\operatorname*{arg\,min}_{{\bm{\beta}}^{\prime}{}}{\left\|\mathbf{T}_{i}^{\top}\left(\mathbf{X}\mathbf{S}_{i}{\bm{\beta}}^{\prime}-\mathbf{y}\right)\right\|}_{2},$
(3) $\displaystyle\widehat{{\bm{\beta}}}{}^{(i)}_{S_{i}^{c}}$
$\displaystyle={\bm{0}},$ (4)

where $S_{i}^{c}=[p]\setminus S_{i}$ denotes the complement of the set
$S_{i}$. This can alternatively be written in closed form as

$\displaystyle\widehat{{\bm{\beta}}}{}^{(i)}=\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{y},$
(5)

where $(\cdot)^{\dagger}$ denotes the Moore–Penrose pseudoinverse. Thus, the
closed-form expression for the ensemble parameter estimate is given by

$\displaystyle\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}=\frac{1}{k}\sum_{i=1}^{k}\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{y}.$
(6)

##  3 ENSEMBLE RISK

We define the _risk_ of a linear predictor as the expected squared error of a
prediction of the target variable on an independent data point $\mathbf{x}$:

$\displaystyle R({\bm{\beta}}^{\prime})$
$\displaystyle\triangleq{\mathbb{E}}_{\mathbf{x}}\left[\left\langle\mathbf{x},{\bm{\beta}}-{\bm{\beta}}^{\prime}\right\rangle^{2}\right]$
$\displaystyle=\left\langle{\bm{\beta}}-{\bm{\beta}}^{\prime},{\bm{\Sigma}}\left({\bm{\beta}}-{\bm{\beta}}^{\prime}\right)\right\rangle.$
(7)

For any predictor of the form
${\bm{\beta}}^{\prime}{}=f(\mathbf{X})\mathbf{y}$, for some
$f:\mathbb{R}^{n\times p}\to\mathbb{R}^{p\times n}$, we can rewrite parameter
estimation error as

$\displaystyle{\bm{\beta}}-{\bm{\beta}}^{\prime}{}=(\mathbf{I}_{p}-f(\mathbf{X})\mathbf{X}){\bm{\beta}}-\sigma
f(\mathbf{X})\mathbf{z}.$ (8)

Then by the independence of $\mathbf{X}$ and $\mathbf{z}$ and some algebra, we
can decompose the risk into the so-called “bias” and “variance” components

$\displaystyle{\mathbb{E}}_{\mathbf{z}}\left[R({\bm{\beta}}^{\prime}{})\right]={}$
$\displaystyle\underbrace{\left\langle{\bm{\beta}}{\bm{\beta}}^{\top},(\mathbf{I}_{p}-f(\mathbf{X})\mathbf{X})^{\top}{\bm{\Sigma}}(\mathbf{I}_{p}-f(\mathbf{X})\mathbf{X})\right\rangle}_{\mathrm{bias}({\bm{\beta}}^{\prime}{})}$
$\displaystyle+\underbrace{\sigma^{2}\left\langle
f(\mathbf{X}),{\bm{\Sigma}}f(\mathbf{X})\right\rangle}_{\mathrm{variance}({\bm{\beta}}^{\prime}{})}.$
(9)

For the ensemble, we obtain for the bias and variance

$\displaystyle\mathrm{bias}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})$
$\displaystyle=\frac{1}{k^{2}}\sum_{i,j=1}^{k}\mathrm{bias}_{ij}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})$
(10) $\displaystyle\mathrm{variance}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})$
$\displaystyle=\frac{1}{k^{2}}\sum_{i,j=1}^{k}\mathrm{variance}_{ij}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}),$
(11)

where

$\displaystyle\mathrm{bias}_{ij}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})$
$\displaystyle={}\Big{\langle}{\bm{\beta}}{\bm{\beta}}^{\top},\left(\mathbf{I}_{p}-\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\right)^{\top}$
$\displaystyle\times{\bm{\Sigma}}\left(\mathbf{I}_{p}-\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\mathbf{X}\right)\Big{\rangle},$
(12)
$\displaystyle\mathrm{variance}_{ij}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})={}\sigma^{2}\Big{\langle}$
$\displaystyle\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top},$
$\displaystyle{\bm{\Sigma}}\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\Big{\rangle}.$
(13)

Thus, evaluating the risk of the ensemble is a matter of evaluating these
pairwise interaction terms.

To begin evaluating the above terms, we need to introduce additional
assumptions. Specifically, we assume that the subsets are independent and that
all indices are equally likely to be included in each subset.

######  Assumption 3.1 (finite subsampling).

The subsets in the collections $\mathcal{S}$ and $\mathcal{T}$ are selected at
random such that $|S_{i}|<|T_{i}|-1$ and that the following hold:

  * •

$\Pr(j\in S_{i})=|S_{i}|/p$ for all $j\in[p]$,

  * •

$\Pr(m\in T_{i})=|T_{i}|/n$ for all $m\in[n]$,

  * •

The subsets $S_{1},S_{2},\ldots,S_{k},T_{1},T_{2},\ldots,T_{k}$ are
conditionally independent given the example subset sizes
$(|T_{i}|)_{i=1}^{k}$.

A simple sampling strategy that satisfies these assumptions is to fix
$|S_{i}|$ and $|T_{i}|$ such that $|S_{i}|<|T_{i}|-1$ and select subsets
uniformly at random of the given sizes. Another strategy is to construct the
subsets by flipping a coin for each index, rejecting any resulting subsets
that fail to satisfy $|S_{i}|<|T_{i}|-1$.

With Assumption 3.1, we are now equipped to evaluate the pairwise interaction
terms. The following two lemmas enable us to characterize the bias and
variance components of the risk in the finite-dimensional setting. The proofs
of these lemmas are exercises in linear algebra and conditional expectations
and can be found in the Appendix.

With some slight abuse of notation, we allow
${\mathbb{E}}_{\mathcal{S},\mathcal{T}}$ to denote the expectation taken with
respect to the choice of indices in the subsets, but not their sizes. In other
words, ${\mathbb{E}}_{\mathcal{S},\mathcal{T}}$ indicates the conditional
expectation over $\mathcal{S}$ and $\mathcal{T}$, conditioned on the subset
sizes indicated by the context.

######  Lemma 3.2 (bias).

Assume that ${\bm{\Sigma}}=\mathbf{I}_{p}$ and that Assumption 3.1 holds. Then

${\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathrm{bias}_{ij}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})\right]\\\
=\begin{cases}\frac{|S_{i}^{c}\cap S_{j}^{c}|}{p}\left(1+\frac{|S_{i}\cap
S_{j}|}{n-|S_{i}\cap
S_{j}|-1}\right){\left\|{\bm{\beta}}\right\|}_{2}^{2}&{\rm{if}}\;i\neq j,\\\
\frac{|S_{i}^{c}|}{p}\left(1+\frac{|S_{i}|}{|T_{i}|-|S_{i}|-1}\right){\left\|{\bm{\beta}}\right\|}_{2}^{2}&{\rm{if}}\;i=j.\end{cases}$
(14)

######  Lemma 3.3 (variance).

Assume that ${\bm{\Sigma}}=\mathbf{I}_{p}$ and that Assumption 3.1 holds. Then

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathrm{variance}_{ij}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})\right]=\begin{cases}\frac{\sigma^{2}|S_{i}\cap
S_{j}|}{n-|S_{i}\cap S_{j}|-1}&{\rm{if}}\;i\neq j,\\\
\frac{\sigma^{2}|S_{i}|}{|T_{i}|-|S_{i}|-1}&{\rm{if}}\;i=j.\end{cases}$ (15)

One observation that we can make already from these results is that the
example subsampling only affects the terms where $i=j$. Assuming that the
subsampling procedure is the same for each $i$, so that for large $k$ the
$i\neq j$ terms are sure to dominate the sum, this means that in the limit as
$k\to\infty$, the effects of example subsampling are non-existent. We note
that this is a result of the assumption that $|S_{i}|<|T_{i}|$, and that if we
were to have $|S_{i}|>|T_{i}|$, then we would observe effects of example
subsampling when $i\neq j$, which we discuss further in Section 5.2.

We now turn our attention to the setting where $n,p\to\infty$ in order to
better reason about the results contained in these lemmas. We introduce the
following additional assumption.

######  Assumption 3.4 (asymptotic subsampling).

For some $\alpha,\eta\in[0,1]$, the subsets in the collections $\mathcal{S}$
and $\mathcal{T}$ are selected randomly such that
$|S_{i}|/p\xrightarrow{{\rm{a.s.}}}\alpha$ as $p\to\infty$ and
$|T_{i}|/n\xrightarrow{{\rm{a.s.}}}\eta$ as $n\to\infty$ for all $i\in[k]$.

This assumption is easily satisfied. For example, in the sampling strategy
where we fix $|S_{i}|$ and $|T_{i}|$, we can choose $|S_{i}|=\lfloor\alpha
p\rfloor$ and $|T_{i}|=\lfloor\eta n\rfloor$. For the coin-flipping strategy,
we can select feature subsets with a coin of probability $\alpha$ and example
subsets with a coin of probability $\eta$.

Under this assumption, and additionally assuming without loss of generality
that ${\left\|{\bm{\beta}}\right\|}_{2}=1$, if $n,p\to\infty$ such that
$p/n\to\gamma$ and $\eta>\alpha\gamma$, the quantities in (14) and (15)
converge almost surely as follows:

${\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathrm{bias}_{ij}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})\right]\\\
\xrightarrow{a.s.}\begin{cases}(1-\alpha)^{2}\left(1+\frac{\alpha^{2}\gamma}{1-\alpha^{2}\gamma}\right)&\text{if
}i\neq j,\\\
(1-\alpha)\left(1+\frac{\alpha\gamma}{\eta-\alpha\gamma}\right)&\text{if
}i=j,\end{cases}$ (16)

and

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathrm{variance}_{ij}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})\right]\xrightarrow{a.s.}\begin{cases}\frac{\sigma^{2}\alpha^{2}\gamma}{1-\alpha^{2}\gamma}&\text{if
}i\neq j,\\\ \frac{\sigma^{2}\alpha\gamma}{\eta-\alpha\gamma}&\text{if
}i=j.\end{cases}$ (17)

We are now equipped to state our asymptotic risk result for the ensemble of
ordinary least squares predictors. Denote for an ensemble satisfying
Assumptions 3.1 and 3.4 with parameters $\alpha$, $\eta$, and $k$ the
_limiting risk_

$\displaystyle
R_{\alpha,\eta,k}^{\mathrm{ens}}\triangleq\lim_{n,p\to\infty}{\mathbb{E}}_{\mathbf{X},\mathbf{z},\mathcal{S},\mathcal{T}}\left[R(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})\right].$
(18)

From (10) and (11), we know that both the bias and variance components of the
limiting risk are the averages of $k^{2}$ terms, and from (16) and (17), we
know that the $k(k-1)$ terms where $i\neq j$ will take one value and the
remaining $k$ terms where $i=j$ will take another. Thus we have the _limiting
bias_

$\lim_{n,p\to\infty}{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathrm{bias}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})\right]\\\
=\frac{k-1}{k}\left(\frac{(1-\alpha)^{2}}{1-\alpha^{2}\gamma}\right)+\frac{1}{k}\left(\frac{\eta(1-\alpha)}{\eta-\alpha\gamma}\right)$
(19)

and _limiting variance_

$\lim_{n,p\to\infty}{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathrm{variance}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})\right]\\\
=\frac{k-1}{k}\left(\frac{\sigma^{2}\alpha^{2}\gamma}{1-\alpha^{2}\gamma}\right)+\frac{1}{k}\left(\frac{\sigma^{2}\alpha\gamma}{\eta-\alpha\gamma}\right).$
(20)

Upon careful examination of these quantities, we observe that in fact both the
limiting bias _and_ the limiting variance are _decreasing_ in $k$, and thus
the ensemble serves not only as a means to reduce variance (as is well
understood), _but also_ to reduce bias. We defer further discussion to Section
4.2. Adding the limiting bias and variance together yields the following
result.

######  Theorem 3.5 (limiting risk).

Assume that ${\bm{\Sigma}}=\mathbf{I}_{p}$ and
${{\left\|{\bm{\beta}}\right\|}_{2}=1}$ and that Assumptions 3.1 and 3.4 hold.
Then in the limit as $n,p\to\infty$ with $p/n\to\gamma$, for
$\eta>\alpha\gamma$, we have almost surely that

$R_{\alpha,\eta,k}^{\mathrm{ens}}=\frac{k-1}{k}\left(\frac{(1-\alpha)^{2}+\sigma^{2}\alpha^{2}\gamma}{1-\alpha^{2}\gamma}\right)\\\
+\frac{1}{k}\left(\frac{\eta(1-\alpha)+\sigma^{2}\alpha\gamma}{\eta-\alpha\gamma}\right).$
(21)

Here we see again more explicitly that for large $k$, the effect of example
subsampling vanishes. This leaves us with the _large-ensemble_ risk

$\displaystyle R_{\alpha}^{\mathrm{ens}}{}$
$\displaystyle\triangleq\lim_{k\to\infty}R_{\alpha,\eta,k}^{\mathrm{ens}}$
$\displaystyle=\frac{(1-\alpha)^{2}+\sigma^{2}\alpha^{2}\gamma}{1-\alpha^{2}\gamma}.$
(22)

We note that while the large-ensemble risk depends only upon $\alpha$, we
cannot realize this risk with an ensemble if $\eta\leq\alpha\gamma$. Our
remaining results concern the large-ensemble risk and therefore assume that
$\eta=1$ for simplicity, but we caution the reader that some of these results
may not be valid for some smaller values of $\eta$, depending on $\sigma$ and
$\gamma$.

Because $\alpha$ is an algorithmic hyperparameter, it can be tuned to minimize
the risk. If we do so, then what we obtain is the perhaps surprising result
that the optimal large-ensemble risk of the ordinary least squares predictor
is equal to the limiting risk of the _ridge regression_ predictor under our
assumptions. The ridge regression predictor with parameter $\lambda$ is
defined as

$\displaystyle\widehat{{\bm{\beta}}}{}_{\lambda}^{\mathrm{ridge}}{}$
$\displaystyle\triangleq\operatorname*{arg\,min}_{{\bm{\beta}}^{\prime}{}}{\left\|\mathbf{X}{\bm{\beta}}^{\prime}{}-\mathbf{y}\right\|}_{2}^{2}+\lambda{\left\|{\bm{\beta}}^{\prime}{}\right\|}_{2}^{2}$
$\displaystyle=\left(\mathbf{X}^{\top}\mathbf{X}+\lambda\mathbf{I}_{p}\right)^{-1}\mathbf{X}^{\top}\mathbf{y}.$
(23)

We formally state this result in the following theorem, which leverages the
recent analysis of the limiting risk of ridge regression by Dobriban and Wager
(2018).111We note that results on MMSE estimation error from the wireless
communication community (see, e.g., Tulino and Verdú, 2004) predate the more
general result of Dobriban and Wager (2018), and that these apply to the
${\bm{\Sigma}}=\mathbf{I}_{p}$ setting we consider, where risk is equal to
estimation error. The proof is found in the Appendix.

######  Theorem 3.6.

Assume that ${\bm{\Sigma}}=\mathbf{I}_{p}$ and
${\bm{\beta}}\sim\mathcal{N}({\bm{0}},p^{-1}\mathbf{I}_{p})$ and that
Assumptions 3.1 and 3.4 hold with $\eta=1$. Then in the limit as
$n,p\to\infty$ with $p/n\to\gamma$, we have almost surely that

$\displaystyle\inf_{\alpha<\gamma^{-1}}R_{\alpha}^{\mathrm{ens}}=\inf_{\lambda}R\left(\widehat{{\bm{\beta}}}^{\mathrm{ridge}}_{\lambda}\right).$
(24)

The implication of Theorem 3.6 is quite strong. Under the assumption of the
theorem that true parameters ${\bm{\beta}}$ have a Gaussian distribution with
covariance $p^{-1}\mathbf{I}_{p}$, the ridge regression predictor (the maximum
a posteriori estimator for this setting) is the predictor with the lowest
expected risk of all predictors of the form
${\bm{\beta}}^{\prime}{}=f(\mathbf{X})\mathbf{y}$. To see this, note that if
we take the expectation of (3) with respect to ${\bm{\beta}}$, we find that
the optimal $f(\mathbf{X})$ must satisfy the first order optimality condition

$\displaystyle{\bm{\Sigma}}f(\mathbf{X})(\mathbf{X}\mathbf{X}^{\top}+p\sigma^{2}\mathbf{I}_{p})={\bm{\Sigma}}\mathbf{X}^{\top},$
(25)

which for invertible ${\bm{\Sigma}}$ yields the optimally tuned ridge
regression predictor. Thus, in the ${\bm{\Sigma}}=\mathbf{I}_{p}$ setting, the
optimally tuned ensemble achieves the optimal risk for any linear predictor.

A curious result obtained during the proof of this theorem is the following
corollary relating the optimal large ensemble risk to the optimal choice of
the hyperparameter $\alpha$.

######  Corollary 3.7.

Assume that ${\bm{\Sigma}}=\mathbf{I}_{p}$ and
${{\left\|{\bm{\beta}}\right\|}_{2}=1}$ and that Assumptions 3.1 and 3.4 hold
with $\eta=1$. Then in the limit as $n,p\to\infty$ with $p/n\to\gamma$, we
have almost surely that

$\displaystyle R_{\alpha_{*}}^{\mathrm{ens}}=1-\alpha_{*},$ (26)

where
$\alpha_{*}=\operatorname*{arg\,min}_{\alpha<\gamma^{-1}}R_{\alpha}^{\mathrm{ens}}$.

##  4 DISCUSSION

###  4.1 Convergence

Figure 2: Approximate limiting risk (averaged over 50 trials with
$n=200,\sigma=1$) when using $\eta=1$ (solid) and $\eta=1.1\times\alpha\gamma$
(dotted). For each value of $\gamma$, both ensembles converge to the
theoretical optimal ridge regression risk (dashed). Figure 3: Approximate
limiting risk (averaged over 100 trials with $n=200,p=400$) when using
$\alpha=\alpha_{*}$ (top) and
$\alpha=\operatorname*{arg\,min}_{\alpha^{\prime}}R_{\alpha^{\prime},\eta,k}^{\mathrm{ens}}$
(bottom). As $k$ increases, in both cases the risk converges to the
theoretical optimal ridge regression risk (black dashed).

In practice, any ensemble will have only a finite number of members.
Therefore, it is important to understand the rates at which the risk of the
ensemble converges to large-ensemble risk in (3). From Theorem 3.5, it is
clear that as a function of $k$, the limiting risk converges to the large-
ensemble risk at a rate $O(1/k)$. However, as the choice of $\eta$ approaches
$\alpha\gamma$, this rate becomes slower. In Figure 2, we plot222See
https://github.com/dlej/ensemble-ols. the convergence in $k$ of the limiting
risk to the large-ensemble risk for $\eta=1$ (using all examples) and for
$\eta=1.1\times(\alpha\gamma)$ (near to as small as possible while still
having $|S_{i}|<|T_{i}|$). We plot these curves for $\sigma=1$ and for three
different values of $\gamma$, using $n=200$, which is sufficient to realize
the convergence in $n$ and $p$. We choose $\alpha=\alpha_{*}$, the minimizer
of the large-ensemble risk. What we observe is that, indeed, for both choices
of $\eta$, the risks converge to the optimal ridge risk. As expected, however,
with the smaller choice of $\eta$ the risk converges nearly an order of
magnitude more slowly.

While the choice of $\alpha=\alpha_{*}$ will result in optimal risk for large
enough ensembles, for finite $k$ this choice can in some cases be undesirable.
For instance, consider the setting where $\eta=1$ and $\gamma>1$. Then as
$\sigma\to 0$, $\alpha_{*}\to\gamma^{-1}$ (see expressions for $\alpha_{*}$ in
the Appendix). This obviously yields the optimal large-ensemble risk, by
definition, but for any finite $k$, the limiting risk tends to infinity for
this choice of $\alpha$. However, if we know what the size of our ensemble
will be, we can tune $\alpha$ to the limiting risk for finite $k$ instead of
the large ensemble risk. In general, this means choosing an $\alpha$ smaller
than $\alpha_{*}$. In Figure 3, we demonstrate the convergence in $k$ to the
large-ensemble risk as a function of $\sigma$ for $\alpha=\alpha_{*}$ and for
$\alpha=\operatorname*{arg\,min}_{\alpha^{\prime}}R_{\alpha^{\prime},\eta,k}^{\mathrm{ens}}$.
We plot these curves for $\gamma=2$ and $\sigma\in[0.1,10]$, using $n=200$.
While for both choices of $\alpha$ we see convergence in $k$ for each
$\sigma$, as $\sigma\to 0$, the risk is very large for $\alpha=\alpha_{*}$.
For $\alpha$ adapted to the choice of $k$, however, this effect is mitigated.

###  4.2 Bias and Variance Decrease with Ensemble Size

We return here to the observation made in Section 3 that the limiting bias and
variance are both _decreasing_ in $k$. Thus, although there is a bias–variance
tradeoff in $\alpha$, there is no such tradeoff with $k$. This can be seen by
comparing the $i=j$ and $i\neq j$ terms in each case. In the case of bias, for
the bias to be decreasing, it must be that

$\displaystyle\frac{(1-\alpha)^{2}}{1-\alpha^{2}\gamma}<\frac{\eta(1-\alpha)}{\eta-\alpha\gamma}.$
(27)

Since $\alpha^{2}\gamma<1$ and $\eta>\alpha\gamma$, after some algebra, this
reduces to

$\displaystyle\gamma(\alpha-1)<\eta(1-\alpha\gamma).$ (28)

Because $\alpha\leq 1$, the left-hand side is non-positive, and since
$\alpha<\gamma^{-1}$, the right-hand side is strictly positive. Thus this
inequality always holds, and the bias is decreasing.

In the case of variance, for the variance to be decreasing, we must have

$\displaystyle\frac{\alpha^{2}\gamma}{1-\alpha^{2}\gamma}<\frac{\alpha\gamma}{\eta-\alpha\gamma}.$
(29)

Again since $\alpha^{2}\gamma<1$ and $\eta>\alpha\gamma$, this reduces to

$\displaystyle\alpha\eta<1.$ (30)

So, unless both $\alpha=1$ and $\eta=1$, in which case every member of the
ensemble is the ordinary least squares predictor fit using the entire training
data, the variance is decreasing.

###  4.3 Dropout and Ridge Regression

There is an interesting connection between the ordinary least squares ensemble
with $\eta=1$ and the popular _dropout_ technique (Srivastava et al., 2014)
used in deep (neural) network training, which consists of randomly masking the
features at each iteration of (stochastic) gradient descent. To draw this
connection, define

$\displaystyle\ell_{i}({\bm{\beta}}^{\prime}{})={\left\|\mathbf{X}\mathbf{S}_{i}\mathbf{S}_{i}^{\top}{\bm{\beta}}^{\prime}-\mathbf{y}\right\|}_{2}^{2}.$
(31)

Then our ensemble member parameter estimates are minimizers of this loss
function.

$\displaystyle\widehat{{\bm{\beta}}}{}^{(i)}=\operatorname*{arg\,min}_{{\bm{\beta}}^{\prime}{}}\ell_{i}({\bm{\beta}}^{\prime}{})\text{
s.t. }{\bm{\beta}}_{S_{i}^{c}}^{\prime}$ $\displaystyle={\bm{0}}.$ (32)

For each $i$, the $i$-th member of the ensemble is able to solve its
subproblem independently of the other members. As a result, we can consider
the ensemble to be a model with $\sum_{i=1}^{k}|S_{i}|$ parameters that are
eventually averaged to reduce them down to $p$ parameters. If we were to
instead constrain ourselves so that we were allowed to use _only_ $p$
parameters, such that we could not optimize each member of the ensemble
independently, we might try to optimize them jointly by minimizing the average
loss. That is,

$\displaystyle\widehat{{\bm{\beta}}}{}=\operatorname*{arg\,min}_{{\bm{\beta}}^{\prime}{}}\frac{1}{k}\sum_{i=1}^{k}\ell_{i}({\bm{\beta}}^{\prime}{}).$
(33)

If we go a step further and let $k\to\infty$ and optimize this loss using
stochastic gradient descent where at each iteration we use the gradient of an
individual $\ell_{i}$ selected at random, then our ensemble becomes equivalent
to the predictor learned using dropout. It is well-known that dropout with
linear regression has a very strong connection to ridge regression (Srivastava
et al., 2014); specifically, we find that

$\displaystyle\widehat{{\bm{\beta}}}{}=\frac{1}{\alpha}\left(\mathbf{X}^{\top}\mathbf{X}+\frac{1-\alpha}{\alpha}\mathrm{diag}(\mathbf{X}^{\top}\mathbf{X})\right)^{-1}\mathbf{X}^{\top}\mathbf{y}.$
(34)

In the case of ${\bm{\Sigma}}=\mathbf{I}_{p}$,
$n^{-1}\mathrm{diag}(\mathbf{X}^{\top}\mathbf{X})$ will converge to
$\mathbf{I}_{p}$ as $n,p\to\infty$, in which case dropout and ridge regression
are equivalent up to a rescaling. We discuss the case where
${\bm{\Sigma}}\neq\mathbf{I}_{p}$ in Section 5.1.

###  4.4 Scaled Ensembles

Figure 4: $\mu$-scaled large-ensemble risk (theoretical, $\gamma=0.5$) when
using $\mu=1$ (solid) and $\mu=\mu_{*}$ (dotted). For both the setting where
we use fewer features than optimal with $\alpha=\alpha_{*}/2$ (blue) and the
fixed $\alpha=1/2$ setting (red), we see significantly improved risk by
scaling.

Our ensemble combines the individual predictors by simple averaging. However,
in light of the fact that dropout is only equivalent to ridge regression up to
a rescaling of the output, it is worth considering the effect of using an
equally-weighted linear combination but using different weights from $1/k$ in
constructing the ensemble predictor. That is, we consider the risk of the
_$\mu$ -scaled_ predictor
$\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}_{\mu}=(\mu/k)\sum_{i=1}^{k}\widehat{{\bm{\beta}}}{}^{(i)}$.
A simple calculation, proved in the Appendix, shows that under the assumptions
of Theorem 3.5 the large-ensemble risk of the $\mu$-scaled predictor is given
by

$\displaystyle
R_{\alpha,\mu}^{\mathrm{ens}}{}=\mu^{2}R_{\alpha}^{\mathrm{ens}}{}+(1-\mu)^{2}+2\mu(1-\mu)(1-\alpha).$
(35)

Hence, it is possible to minimize the risk of
$\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}_{\mu}$ over the choice of parameter
$\mu$. This results in

$\displaystyle\mu_{*}$
$\displaystyle=\frac{\alpha}{R_{\alpha}^{\mathrm{ens}}{}+2\alpha-1}$ (36)

as the optimal choice for $\mu$ and

$\displaystyle R_{\alpha,\mu_{*}}^{\mathrm{ens}}{}$
$\displaystyle=1-\frac{\alpha^{2}}{2\alpha-1+R_{\alpha}^{\mathrm{ens}}{}}$
(37)

as the achieved risk for the optimally-scaled ensemble. Note that as a result
of Corollary 3.7, $R_{\alpha_{*}}^{\mathrm{ens}}{}=1-\alpha_{*}$. Therefore,
for ensembles with optimally-tuned $\alpha=\alpha_{*}$ we have $\mu_{*}=1$,
and any scaling in constructing the ensemble predictor will not further
improve the achieved risk. However, it is easy to see that when
$\alpha>\alpha_{*}$ (the ensemble members select _more_ features than is
optimal), $\mu_{*}<1$, and the risk is improved by adding extra shrinkage to
the ensemble predictor. Similarly, if $\alpha<\alpha_{*}$, (the ensemble
members select _less_ features than is optimal), $\mu_{*}>1$, and the risk is
improved by inflating the ensemble predictor. We illustrate the improvement in
risk to be had in Figure 4, where we plot the risk with ($\mu=\mu_{*}$) and
without ($\mu=1$) optimal scaling for two choices of $\alpha$—one where we
always select half as many features as optimal ($\alpha=\alpha_{*}/2)$, and
one where we always use half of the available features $(\alpha=1/2)$.

##  5 FUTURE DIRECTIONS

###  5.1 Non-Identity Covariance

Of course, it is important to understand the behavior of the ordinary least
squares ensemble in the case where ${\bm{\Sigma}}\neq\mathbf{I}_{p}$ when
considering applications of the method to real data. As discussed in Section
3, provided ${\bm{\Sigma}}$ is invertible, ridge regression remains the
optimal linear predictor, and whether the ensemble (or extensions thereto)
still achieves the optimal risk in this setting remains an open question.

By inspection of the closed-form solution of dropout in (34), we see that it
is no longer equivalent (as $n,p\to\infty$) to ridge regression in this
setting and is therefore no longer optimal. We believe that this is likely the
case for the ensemble as well. However, if we extend the coin-flipping
strategy for feature subset selection to one where we have a collection of
coin with probabilities ${\bm{\alpha}}\in[0,1]^{p}$, one for each feature, we
can extend the result in (34) to obtain the closed-form dropout solution

$\displaystyle\widehat{{\bm{\beta}}}{}=\mathbf{A}^{-1}\left(\mathbf{X}^{\top}\mathbf{X}+\left(\mathbf{I}_{p}-\mathbf{A}\right)\mathbf{A}^{-1}\mathrm{diag}(\mathbf{X}^{\top}\mathbf{X})\right)^{-1}\mathbf{X}^{\top}\mathbf{y},$
(38)

where $\mathbf{A}=\mathrm{diag}({\bm{\alpha}})$. We prove this result in the
Appendix. Thus, if ${\bm{\alpha}}$ is chosen such that

$\displaystyle\frac{1-\alpha_{j}}{\alpha_{j}}=\frac{\lambda}{n[{\bm{\Sigma}}]_{jj}},$
(39)

then the corrected dropout estimator

$\displaystyle\widetilde{{\bm{\beta}}}=\mathbf{A}\widehat{{\bm{\beta}}}{}$
(40)

is equivalent to ridge regression with parameter $\lambda$ as $n,p\to\infty$.
This leads us to believe that the optimal ensemble in the
${\bm{\Sigma}}\neq\mathbf{I}_{p}$ setting should also use non-uniform feature
sampling, and extending our analysis to this case is an interesting area for
future work.

###  5.2 Beyond Ordinary Least Squares: Ensembles of Interpolators

Throughout this work we have assumed that the members of the ensemble solve
their subproblems using ordinary least squares, which yields the unique
solution that minimizes the squared error given $|T_{i}|$ observations of
$|S_{i}|$ variables, and this uniqueness requires that $|T_{i}|$ be no less
than $|S_{i}|$. In the case where $|T_{i}|<|S_{i}|$, there are infinitely many
solutions that minimize the squared error. However, we could in this case opt
to _regularize_ the solution to solve this problem. While analysis of the
effect of regularizing the solution of the subproblems in the ensemble is
beyond the scope of this work, we comment briefly on what would happen if we
were to simply use the same solution presented in (5)—i.e., use the
pseudoinverse solution, which has the smallest $\ell^{2}$ norm of all
solutions to the least squares problem. In this case, when $\eta=1$, the
learned predictor would be an _interpolator_ (Belkin et al., 2018; Hastie et
al., 2019) of the training data, and such methods have recently become
increasingly of interest given the ability of deep (neural) network methods to
have extremely good test performance while having (nearly) zero training error
(Zhang et al., 2017; Belkin et al., 2019).

Specifically, it becomes immediately clear that in this setting, the effect of
the choice of $\eta$ does not vanish as $k\to\infty$. Lemma 3.3 can easily be
extended to this setting, since the roles of $S_{i}$ and $T_{i}$ in (3) can
simply be reversed, and as $n,p\to\infty$, we obtain

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathrm{variance}_{ij}(\widehat{{\bm{\beta}}}{}^{\mathrm{ens}})\right]\xrightarrow{a.s.}\begin{cases}\frac{\sigma^{2}\eta^{2}}{\gamma-\eta^{2}}&\text{if
}i\neq j,\\\ \frac{\sigma^{2}\eta}{\alpha\gamma-\eta}&\text{if
}i=j.\end{cases}$ (41)

Thus, the variance component of the large-ensemble risk in this setting is
equal to $\sigma^{2}\eta^{2}/(\gamma-\eta^{2})$ and does not depend upon
$\alpha$. In future work, we plan to extend our analysis for the bias
component of the large-ensemble risk to this setting, and we expect that in
this case the bias will depend on _both_ $\alpha$ and $\eta$.

###  5.3 Optimal Ensemble Mixing

In the ordinary least squares ensemble, we have used equal weighting when
taking the average of our predictors. Instead, we could extend the idea
presented in Section 4.4 to consider unequal weighting parameterized by
${\bm{\mu}}\in\mathbb{R}^{k}$, giving us the ensemble parameter estimate
$\widehat{{\bm{\beta}}}{}_{\bm{\mu}}^{\mathrm{ens}}=\sum_{i=1}^{k}\mu_{i}\widehat{{\bm{\beta}}}{}^{(i)}$.
While equal weighting gives us optimal risk in the setting where
${\bm{\beta}}\sim\mathcal{N}({\bm{0}},p^{-1}\mathbf{I}_{p})$, where ridge
regression is optimal, under other distributional assumptions on
${\bm{\beta}}$, such as sparsity, where ridge regression is not optimal,
unequal weighting has the potential to yield better ensembles.

Using the sparsity example, consider ${\bm{\beta}}$ such that
${\left\|{\bm{\beta}}\right\|}_{0}=s\ll p$, and suppose that for some $i$,
$S_{i}=S_{\bm{\beta}}$, where $S_{\bm{\beta}}=\left\\{j:\beta_{j}\neq
0\right\\}$. For simplicity, assume that $\eta=1$, so that $T_{i}=[n]$ for all
$i$. In this case, any predictor that uses the remaining $p-s$ features
injects noise into its predictions, so the best predictor uses only the $s$
features in $S_{\bm{\beta}}$. Under the i.i.d. Gaussian noise assumption, the
predictor with lowest risk is in fact

$\displaystyle\widehat{{\bm{\beta}}}$
$\displaystyle=\operatorname*{arg\,min}_{{\bm{\beta}}^{\prime}{}:{\bm{\beta}}_{S_{\bm{\beta}}^{c}}^{\prime}={\bm{0}}}{\left\|\mathbf{y}-\mathbf{X}{\bm{\beta}}^{\prime}{}\right\|}_{2}=\widehat{{\bm{\beta}}}{}^{(i)},$
(42)

where $i$ is such that $S_{i}=S_{\bm{\beta}}$. Thus an optimal weighting
${\bm{\mu}}$ is given by

$\displaystyle\mu_{i}=\begin{cases}\frac{1}{C}&\text{if
}S_{i}=S_{\bm{\beta}},\\\ 0&\text{otherwise},\end{cases}$ (43)

where $C=|\left\\{i:S_{i}=S_{\bm{\beta}}\right\\}|$. This optimal weighting is
decidedly non-uniform, and this raises the question of what schemes could be
employed, either adaptively or non-adaptively, to minimize risk, and how they
would fit into this analysis framework.

#### Acknowledgements

We would like to thank Ryan Tibshirani for helpful discussions and the
anonymous reviewers for their helpful feedback. This work was supported by NSF
grants CCF-1911094, IIS-1838177, and IIS-1730574; ONR grants N00014-18-12571
and N00014-17-1-2551; AFOSR grant FA9550-18-1-0478; DARPA grant G001534-7500;
and a Vannevar Bush Faculty Fellowship, ONR grant N00014-18-1-2047.

## References

  * Amit and Geman (1997) Y. Amit and D. Geman.  Shape quantization and recognition with randomized trees.  _Neural Computation_ , 9(7):1545–1588, 1997.
  * Arora et al. (2019) S. Arora, N. Cohen, W. Hu, and Y. Luo.  Implicit regularization in deep matrix factorization.  In _Advances in Neural Information Processing Systems 32_ , pages 7413–7424. 2019.
  * Athey et al. (2019) S. Athey, J. Tibshirani, and S. Wager.  Generalized random forests.  _The Annals of Statistics_ , 47(2):1148–1178, Apr. 2019.
  * Belkin et al. (2018) M. Belkin, D. J. Hsu, and P. Mitra.  Overfitting or perfect fitting? Risk bounds for classification and regression rules that interpolate.  In _Advances in Neural Information Processing Systems 31_ , pages 2300–2311. 2018.
  * Belkin et al. (2019) M. Belkin, D. Hsu, S. Ma, and S. Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  _Proceedings of the National Academy of Sciences_ , 116(32):15849–15854, 2019.
  * Breiman (1996) L. Breiman.  Bagging predictors.  _Machine Learning_ , 24(2):123–140, Aug. 1996\.
  * Breiman (1998) L. Breiman.  Arcing classifier (with discussion and a rejoinder by the author).  _The Annals of Statistics_ , 26(3):801–849, June 1998.
  * Breiman (2001) L. Breiman.  Random forests.  _Machine Learning_ , 45(1):5–32, Oct. 2001.
  * Bühlmann and Yu (2002) P. Bühlmann and B. Yu.  Analyzing bagging.  _The Annals of Statistics_ , 30(4):927–961, Aug. 2002.
  * Cutler and Zhao (2001) A. Cutler and G. Zhao.  PERT - perfect random tree ensembles.  _Computing Science and Statistics_ , page 497, 2001.
  * Dobriban and Wager (2018) E. Dobriban and S. Wager.  High-dimensional asymptotics of prediction: Ridge regression and classification.  _The Annals of Statistics_ , 46(1):247–279, Feb. 2018.
  * Fernández-Delgado et al. (2014) M. Fernández-Delgado, E. Cernadas, S. Barro, and D. Amorim.  Do we need hundreds of classifiers to solve real world classification problems?  _Journal of Machine Learning Research_ , 15:3133–3181, 2014\.
  * Freund and Schapire (1997) Y. Freund and R. E. Schapire.  A decision-theoretic generalization of on-line learning and an application to boosting.  _Journal of Computer and System Sciences_ , 55(1):119–139, 1997.
  * Friedberg et al. (2018) R. Friedberg, J. Tibshirani, S. Athey, and S. Wager.  Local linear forests.  _arXiv preprint arXiv:1807.11408_ , 2018.
  * Friedman et al. (2001) J. Friedman, T. Hastie, and R. Tibshirani.  _The Elements of Statistical Learning_.  Springer Series in Statistics, 2001.
  * Gunasekar et al. (2017) S. Gunasekar, B. E. Woodworth, S. Bhojanapalli, B. Neyshabur, and N. Srebro.  Implicit regularization in matrix factorization.  In _Advances in Neural Information Processing Systems 30_ , pages 6151–6159. 2017.
  * Hardt et al. (2016) M. Hardt, B. Recht, and Y. Singer.  Train faster, generalize better: Stability of stochastic gradient descent.  In _Proceedings of The 33rd International Conference on Machine Learning_ , volume 48, pages 1225–1234, June 2016.
  * Hastie et al. (2019) T. Hastie, A. Montanari, S. Rosset, and R. J. Tibshirani.  Surprises in high-dimensional ridgeless least squares interpolation.  _arXiv preprint arXiv:1903.08560_ , 2019.
  * Hoerl and Kennard (1970) A. E. Hoerl and R. W. Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  _Technometrics_ , 12(1):55–67, 1970.
  * Josse and Wager (2016) J. Josse and S. Wager.  Bootstrap-based regularization for low-rank matrix estimation.  _Journal of Machine Learning Research_ , 17(1):4227–4255, Jan. 2016.
  * Louppe and Geurts (2012) G. Louppe and P. Geurts.  Ensembles on random patches.  In _Machine Learning and Knowledge Discovery in Databases_ , pages 346–361, Berlin, Heidelberg, 2012.
  * Mahoney (2012) M. W. Mahoney.  Approximate computation and implicit regularization for very large-scale data analysis.  In _Proceedings of the 31st ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems_ , PODS ’12, pages 143–154, 2012.
  * Neyshabur et al. (2014) B. Neyshabur, R. Tomioka, and N. Srebro.  In search of the real inductive bias: On the role of implicit regularization in deep learning.  _arXiv preprint arXiv:1412.6614_ , 2014.
  * Scornet et al. (2015) E. Scornet, G. Biau, and J.-P. Vert.  Consistency of random forests.  _The Annals of Statistics_ , 43(4):1716–1741, Aug. 2015.
  * Srivastava et al. (2014) N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _Journal of Machine Learning Research_ , 15(1):1929–1958, Jan. 2014.
  * Tulino and Verdú (2004) A. M. Tulino and S. Verdú.  Random matrix theory and wireless communications.  _Foundations and Trends® in Communications and Information Theory_ , 1(1):1–182, 2004.
  * Wager et al. (2013) S. Wager, S. Wang, and P. S. Liang.  Dropout training as adaptive regularization.  In _Advances in Neural Information Processing Systems 26_ , pages 351–359. 2013.
  * Wager et al. (2014) S. Wager, T. Hastie, and B. Efron.  Confidence intervals for random forests: The jackknife and the infinitesimal jackknife.  _Journal of Machine Learning Research_ , 15:1625–1651, 2014\.
  * Wyner et al. (2017) A. J. Wyner, M. Olson, J. Bleich, and D. Mease.  Explaining the success of AdaBoost and random forests as interpolating classifiers.  _Journal of Machine Learning Research_ , 18(48):1–33, 2017.
  * Zhang et al. (2017) C. Zhang, S. Bengio, M. Hardt, B. Recht, and O. Vinyals.  Understanding deep learning requires rethinking generalization.  In _5th International Conference on Learning Representations, ICLR 2017_ , 2017.

##  Appendix A USEFUL LEMMAS

The following two lemmas will be useful in deriving the bias and variance
terms of the ensemble risk. Their proofs can be found in Section F.

######  Lemma A.1.

Let $S\subseteq[p]$ be a subset with corresponding selection matrix
$\mathbf{S}$, and let $\mathbf{S}^{c}$ be the selection matrix corresponding
to $S^{c}$. Then for a random matrix $\mathbf{X}\in\mathbb{R}^{n\times p}$
with rows independently drawn from $\mathcal{N}(\bm{0},\mathbf{I}_{p})$ such
that $n>|S|$, and for any random function
$f:\mathbb{R}^{n\times|S|}\to\mathbb{R}^{n\times|S|}$ that
$f(\mathbf{X}\mathbf{S})$ and $\mathbf{X}\mathbf{S}^{c}$ are independent,

$\displaystyle{\mathbb{E}}_{\mathbf{X}\mathbf{S}^{c}}\left[\mathbf{S}^{\top}\mathbf{X}^{\dagger}\right]=\left(\mathbf{X}\mathbf{S}\right)^{\dagger}$
(44)

and

$\displaystyle{\mathbb{E}}_{\mathbf{X}\mathbf{S}^{c}}\left[{\mathbf{S}^{c}}^{\top}\mathbf{X}^{\top}f(\mathbf{X}\mathbf{S})\mathbf{S}^{\top}\mathbf{X}^{\dagger}\right]=\bf{0}.$
(45)

######  Lemma A.2.

Let $T_{1},T_{2}\subseteq[n]$ be independent random subsets with corresponding
selection matrices $\mathbf{T}_{1},\mathbf{T}_{2}$ such that
${\mathbb{E}}\left[\mathbf{T}_{j}\mathbf{T}_{j}^{\top}\right]=\frac{|T_{j}|}{n}\mathbf{I}_{n}$.
Then for random matrices $\mathbf{X}\in\mathbb{R}^{n\times
p_{X}},\mathbf{Y}\in\mathbb{R}^{n\times p_{Y}}$ independent of $T_{1}$ and
$T_{2}$ with independent and identically distributed rows such that
$\mathbf{X}^{\top}\mathbf{T}_{j}\mathbf{T}_{j}^{\top}\mathbf{X}$ and
$\mathbf{Y}^{\top}\mathbf{T}_{j}\mathbf{T}_{j}^{\top}\mathbf{Y}$ are
invertible, and for any matrix $\mathbf{A}\in\mathbb{R}^{p_{X}\times p_{Y}}$,

$\displaystyle{\mathbb{E}}_{T_{1},T_{2}}\left[\left(\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{\dagger}\mathbf{T}_{1}^{\top}\left(\left(\mathbf{T}_{2}^{\top}\mathbf{X}\right)^{\dagger}\mathbf{T}_{2}^{\top}\right)^{\top}\right]=\left(\mathbf{X}^{\top}\mathbf{X}\right)^{\dagger}$
(46)

and

$\displaystyle{\mathbb{E}}_{T_{1},T_{2}}\left[\left(\left(\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{\dagger}\mathbf{T}_{1}^{\top}\right)^{\top}\mathbf{A}\left(\mathbf{T}_{2}^{\top}\mathbf{Y}\right)^{\dagger}\mathbf{T}_{2}^{\top}\right]=\left(\mathbf{X}^{\dagger}\right)^{\top}\mathbf{A}\mathbf{Y}^{\dagger}.$
(47)

##  Appendix B PROOF OF LEMMA 3.2 (BIAS)

To compute the bias, we need to evaluate terms of the form

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left\langle{\bm{\beta}}{\bm{\beta}}^{\top},\left(\mathbf{I}_{p}-\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\right)^{\top}\left(\mathbf{I}_{p}-\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\mathbf{X}\right)\right\rangle.$
(48)

First, we note that since
$\mathbf{S}_{i}\mathbf{S}_{i}^{\top}+\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}=\mathbf{I}_{p}$,

$\displaystyle\mathbf{I}_{p}-\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}$
$\displaystyle=\mathbf{I}_{p}-\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\left(\mathbf{S}_{i}\mathbf{S}_{i}^{\top}+\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}\right)$
(49)
$\displaystyle=\mathbf{I}_{p}-\mathbf{S}_{i}\mathbf{S}_{i}^{\top}-\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}$
(50)
$\displaystyle=\left(\mathbf{I}_{p}-\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\right)\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}.$
(51)

So, we can equivalently evaluate

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left\langle{\bm{\beta}}{\bm{\beta}}^{\top},\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}\left[\mathbf{I}_{p}-\mathbf{X}^{\top}\mathbf{T}_{i}\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\mathbf{T}_{i}\right)^{\dagger}\mathbf{S}_{i}^{\top}\right]\left[\mathbf{I}_{p}-\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\mathbf{X}\right]\mathbf{S}_{j}^{c}{\mathbf{S}_{j}^{c}}^{\top}\right\rangle.$
(52)

It suffices to evaluate the expectation of the second argument of the inner
product:

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\bigg{[}$
$\displaystyle\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}\left[\mathbf{I}_{p}-\mathbf{X}^{\top}\mathbf{T}_{i}\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\mathbf{T}_{i}\right)^{\dagger}\mathbf{S}_{i}^{\top}\right]\left[\mathbf{I}_{p}-\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\mathbf{X}\right]\mathbf{S}_{j}^{c}{\mathbf{S}_{j}^{c}}^{\top}\bigg{]}$
$\displaystyle={}$
$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\bigg{[}\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}\mathbf{X}^{\top}\mathbf{T}_{i}\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\mathbf{T}_{i}\right)^{\dagger}\mathbf{S}_{i}^{\top}\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{i}^{\top}\mathbf{S}_{j}^{c}{\mathbf{S}_{j}^{c}}^{\top}$
$\displaystyle-\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}\mathbf{X}^{\top}\mathbf{T}_{i}\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\mathbf{T}_{i}\right)^{\dagger}\mathbf{S}_{i}^{\top}-\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}^{c}{\mathbf{S}_{j}^{c}}^{\top}+\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}\mathbf{S}_{j}^{c}{\mathbf{S}_{j}^{c}}^{\top}\bigg{]}.$
(53)

The second and third terms are zero in expectation. To see this for the second
term, observe that ${\mathbf{S}_{i}^{c}}^{\top}\mathbf{X}^{\top}$ and
$\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}$ are independent and each zero-mean.
An analogous argument applies to the third term. The fourth term is equal to

$\displaystyle\frac{|S_{i}^{c}\cap S_{j}^{c}|}{p}\mathbf{I}_{p}.$ (54)

We now consider the case where $i\neq j$. To evaluate the first term, we first
apply Lemma A.2. This simplifies the expression to

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{S}}\left[\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}\mathbf{X}^{\top}\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\right)^{\dagger}\mathbf{S}_{i}^{\top}\mathbf{S}_{j}\left(\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{X}\mathbf{S}_{i}^{\top}\mathbf{S}_{j}^{c}{\mathbf{S}_{j}^{c}}^{\top}\right].$
(55)

Now let $\mathbf{S}_{i\cap j}$, $\mathbf{S}_{i\setminus j}$,
$\mathbf{S}_{j\setminus i}$, and $\mathbf{S}_{i\cup j}^{c}$ denote the
selection matrices corresponding to the sets $S_{i}\cap S_{j}$,
$S_{i}\setminus S_{j}$, $S_{j}\setminus S_{i}$, and $S_{i}^{c}\cap S_{j}^{c}$,
respectively. Without loss of generality, consider when
$\mathbf{S}_{i}^{c}=\begin{bmatrix}\mathbf{S}_{j\setminus
i}\;\mathbf{S}_{i\cup j}^{c}\end{bmatrix}$ and
$\mathbf{S}_{j}^{c}=\begin{bmatrix}\mathbf{S}_{i\setminus
j}\;\mathbf{S}_{i\cup j}^{c}\end{bmatrix}$. Then the matrix inside this
expectation is of the form

$\displaystyle\mathbf{S}_{i}^{c}\begin{bmatrix}\mathbf{A}&\mathbf{B}\\\
\mathbf{C}&\mathbf{D}\end{bmatrix}{\mathbf{S}_{j}^{c}}^{\top},$ (56)

where

$\displaystyle\mathbf{A}$ $\displaystyle=\mathbf{S}_{j\setminus
i}^{\top}\mathbf{X}^{\top}\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\right)^{\dagger}\mathbf{S}_{i}^{\top}\mathbf{S}_{j}\left(\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{X}\mathbf{S}_{i\setminus
j}$ (57) $\displaystyle\mathbf{B}$ $\displaystyle=\mathbf{S}_{j\setminus
i}^{\top}\mathbf{X}^{\top}\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\right)^{\dagger}\mathbf{S}_{i}^{\top}\mathbf{S}_{j}\left(\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{X}\mathbf{S}_{i\cup
j}^{c}$ (58) $\displaystyle\mathbf{C}$ $\displaystyle={\mathbf{S}_{i\cup
j}^{c}}^{\top}\mathbf{X}^{\top}\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\right)^{\dagger}\mathbf{S}_{i}^{\top}\mathbf{S}_{j}\left(\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{X}\mathbf{S}_{i\setminus
j}$ (59) $\displaystyle\mathbf{D}$ $\displaystyle={\mathbf{S}_{i\cup
j}^{c}}^{\top}\mathbf{X}^{\top}\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\right)^{\dagger}\mathbf{S}_{i}^{\top}\mathbf{S}_{j}\left(\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{X}\mathbf{S}_{i\cup
j}^{c}.$ (60)

In the case of $\mathbf{B}$ and $\mathbf{C}$, because
$\mathbf{X}\mathbf{S}_{i\cup j}^{c}$ is independent of the remainder of the
factors, ${\mathbb{E}}_{\mathbf{X}}\left[\mathbf{B}\right]$ and
${\mathbb{E}}_{\mathbf{X}}\left[\mathbf{C}\right]$ are equal to $\bm{0}$. By
applying the second claim of Lemma A.1, we observe that
${\mathbb{E}}_{\mathbf{X}}\left[\mathbf{A}\right]$ is also equal to $\bm{0}$.
This leaves

$\displaystyle{\mathbb{E}}_{\mathbf{X}}\left[\mathbf{D}\right]$
$\displaystyle={\mathbb{E}}_{\mathbf{X}}\left[{\mathbf{S}_{i\cup
j}^{c}}^{\top}\mathbf{X}^{\top}{\mathbb{E}}_{\mathbf{X}\mathbf{S}_{j\setminus
i}}\left[\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\right)^{\dagger}\mathbf{S}_{i}^{\top}\mathbf{S}_{i\cap
j}\right]{\mathbb{E}}_{\mathbf{X}\mathbf{S}_{j\setminus
i}}\left[\mathbf{S}_{i\cap
j}^{\top}\mathbf{S}_{j}\left(\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\right]\mathbf{X}\mathbf{S}_{i\cup
j}^{c}\right]$ (61)
$\displaystyle={\mathbb{E}}_{\mathbf{X}}\left[{\mathbf{S}_{i\cup
j}^{c}}^{\top}\mathbf{X}^{\top}\left(\mathbf{X}\mathbf{S}_{i\cap
j}\mathbf{S}_{i\cap
j}^{\top}\mathbf{X}^{\top}\right)^{\dagger}\mathbf{X}\mathbf{S}_{i\cup
j}^{c}\right].$ (62)

We can evaluate the expectation of the pseudoinverse on its own since
$\mathbf{X}\mathbf{S}_{i\cap j}$ and $\mathbf{X}\mathbf{S}_{i\cup j}^{c}$ are
independent. This matrix has a generalized inverse Wishart distribution with
scale matrix $\mathbf{I}_{n}$ and $|S_{i}\cap S_{j}|$ degrees of freedom,
which yields

$\displaystyle{\mathbb{E}}_{\mathbf{X}}\left[\left(\mathbf{X}\mathbf{S}_{i\cap
j}\mathbf{S}_{i\cap j}^{\top}\mathbf{X}^{\top}\right)^{\dagger}\right]$
$\displaystyle=\frac{|S_{i}\cap S_{j}|}{n(n-|S_{i}\cap
S_{j}|-1)}\mathbf{I}_{n}.$ (63)

This leaves

$\displaystyle{\mathbb{E}}_{\mathbf{X}}\left[{\mathbf{S}_{i\cup
j}^{c}}^{\top}\mathbf{X}^{\top}\left(\frac{|S_{i}\cap S_{j}|}{n(n-|S_{i}\cap
S_{j}|-1)}\mathbf{I}_{n}\right)\mathbf{X}\mathbf{S}_{i\cup j}^{c}\right]$
$\displaystyle=\frac{|S_{i}\cap S_{j}|}{n-|S_{i}\cap
S_{j}|-1}\mathbf{I}_{|S_{i}^{c}\cap S_{j}^{c}|}.$ (64)

Then the expectation in (55) becomes

$\displaystyle{\mathbb{E}}_{\mathcal{S}}\left[\frac{|S_{i}\cap
S_{j}|}{(n-|S_{i}\cap
S_{j}|-1)}\mathbf{S}_{i}^{c}\begin{bmatrix}\bm{0}&\bm{0}\\\
\bm{0}&\mathbf{I}_{|S_{i}^{c}\cap
S_{j}^{c}|}\end{bmatrix}{\mathbf{S}_{j}^{c}}^{\top}\right]$
$\displaystyle=\frac{|S_{i}\cap S_{j}||S_{i}^{c}\cap
S_{j}^{c}|}{p(n-|S_{i}\cap S_{j}|-1)}\mathbf{I}_{p},$ (65)

and combing with (54), we have that the bias is equal to

$\displaystyle\frac{|S_{i}^{c}\cap S_{j}^{c}|}{p}\left(1+\frac{|S_{i}\cap
S_{j}|}{n-|S_{i}\cap S_{j}|-1}\right){\left\|{\bm{\beta}}\right\|}_{2}^{2}.$
(66)

When $i=j$, by a similar argument, without the need to apply Lemma A.2, it
follows that the bias is equal to

$\displaystyle\frac{|S_{i}^{c}|}{p}\left(1+\frac{|S_{i}|}{|T_{i}|-|S_{i}|-1}\right){\left\|{\bm{\beta}}\right\|}_{2}^{2}.$
(67)

##  Appendix C PROOF OF LEMMA 3.3 (VARIANCE)

To compute the variance, we need to evaluate the terms of the form

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{T}}\left\langle\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top},\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\right\rangle.$
(68)

Let $\mathbf{S}$ be the selection matrix corresponding to the set $S_{i}\cap
S_{j}$. Then

$\displaystyle{\mathbb{E}}\left\langle\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top},\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\right\rangle$
$\displaystyle={\mathbb{E}}\left\langle\mathbf{S}^{\top}\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top},\mathbf{S}^{\top}\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\mathbf{T}_{j}^{\top}\right\rangle$
(69)
$\displaystyle={\mathbb{E}}\left\langle{\mathbb{E}}_{\mathbf{X}^{S_{i}\setminus
S_{j}}}\left[\mathbf{S}^{\top}\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\right]\mathbf{T}_{i}^{\top},{\mathbb{E}}_{\mathbf{X}^{S_{j}\setminus
S_{i}}}\left[\mathbf{S}^{\top}\mathbf{S}_{j}\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}_{j}\right)^{\dagger}\right]\mathbf{T}_{j}^{\top}\right\rangle$
(70)
$\displaystyle={\mathbb{E}}\left\langle\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}\right)^{\dagger}\mathbf{T}_{i}^{\top},\left(\mathbf{T}_{j}^{\top}\mathbf{X}\mathbf{S}\right)^{\dagger}\mathbf{T}_{j}^{\top}\right\rangle.$
(71)

The equality (71) is the result of two applications of Lemma A.1.

In the case that $i\neq j$, an application of Lemma A.2 simplifies the above
to

$\displaystyle{\rm
tr}\left({\mathbb{E}}_{\mathbf{X}}\left[\left(\mathbf{S}^{\top}\mathbf{X}^{\top}\mathbf{X}\mathbf{S}\right)^{-1}\right]\right)=\frac{|S_{i}\cap
S_{j}|}{n-|S_{i}\cap S_{j}|-1}.$ (72)

The equality comes from
$\left(\mathbf{S}^{\top}\mathbf{X}^{\top}\mathbf{X}\mathbf{S}\right)^{-1}$
having an inverse Wishart distribution with scale matrix
$\mathbf{I}_{|S_{i}\cap S_{j}|}$ and $n$ degrees of freedom.

When $i=j$, we obtain a similar result directly without needing Lemma A.2. The
above simplifies to

$\displaystyle{\rm
tr}\left({\mathbb{E}}_{\mathbf{X}}\left[\left(\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\mathbf{T}_{i}\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{-1}\right]\right)=\frac{|S_{i}|}{|T_{i}|-|S_{i}|-1}.$
(73)

##  Appendix D PROOF OF THEOREM 3.6

We first introduce the result due to Dobriban and Wager (2018). We note again,
as we noted in the main text, that in the setting of
${\bm{\Sigma}}=\mathbf{I}_{p}$, where the optimal ridge regression risk is
equal to the estimation error of the minimum mean squared error (MMSE)
estimator, results on the value of this quantity predate the result of
Dobriban and Wager (2018). We refer the reader, for example, to the wireless
communication literature (see, e.g., Tulino and Verdú, 2004). However,
Dobriban and Wager (2018) have developed the first results on ridge regression
risk for general ${\bm{\Sigma}}$, and their clean theorem statement is simple
and straightforward to use, even in the ${\bm{\Sigma}}=\mathbf{I}_{p}$ case.

######  Proposition D.1 (from Dobriban and Wager, 2018, Theorem 2.1).

Assume that $\Sigma=\mathbf{I}_{p}$ and
${\bm{\beta}}\sim\mathcal{N}({\bm{0}},p^{-1}\mathbf{I}_{p})$. Then in the
limit as $n,p\to\infty$ with $p/n\to\gamma$, we have almost surely that

$\displaystyle\inf_{\lambda}R(\widehat{{\bm{\beta}}}{}_{\lambda}^{\mathrm{ridge}})=\frac{1}{2}\left(\frac{\gamma-1}{\gamma}-\sigma^{2}+\sqrt{\left(\sigma^{2}-\frac{\gamma-1}{\gamma}\right)^{2}+4\sigma^{2}}\right).$
(74)

We note that this expression is equal to
$\sigma^{2}(R^{*}(1/\sigma^{2},\gamma)-1)$ in the notation of Dobriban and
Wager (2018), where this transformation is necessary because we assume
${\left\|{\bm{\beta}}\right\|}_{2}=1$ rather than $\sigma=1$ and because we
evaluate the noise-free risk.

The minimizer of the large ensemble risk should satisfy the first-order
optimality condition, so we begin by taking its derivative.

$\displaystyle\frac{dR_{\alpha}^{\mathrm{ens}}}{d\alpha}$
$\displaystyle=\frac{(-2(1-\alpha)+2\sigma^{2}\alpha\gamma)(1-\alpha^{2}\gamma)-((1-\alpha)^{2}+\sigma^{2}\alpha^{2}\gamma)(-2\alpha\gamma)}{(1-\alpha^{2}\gamma)^{2}}$
(75)
$\displaystyle=\frac{-\alpha^{2}\gamma+(\gamma(\sigma^{2}+1)+1)-1}{(1-\alpha^{2}\gamma)^{2}}.$
(76)

Thus the minimizer $\alpha_{*}$ should satisfy

$\displaystyle\alpha_{*}^{2}\gamma-\alpha_{*}(\gamma(\sigma^{2}+1)+1)+1=0.$
(77)

From here, it is simply a matter of cumbersome algebra to show that the choice

$\displaystyle\alpha_{*}=\frac{\gamma(\sigma^{2}+1)+1-\sqrt{(\gamma(\sigma^{2}+1)+1)^{2}-4\gamma}}{2\gamma}$
(78)

is the valid root of this quadratic expression and is such that
$R_{\alpha_{*}}^{\mathrm{ens}}=\inf_{\lambda}R(\widehat{{\bm{\beta}}}{}_{\lambda}^{\mathrm{ridge}})$.
We here show a slightly more interesting approach, leading to Corollary 3.7.
First, we start from (77) and add a root of $\alpha_{*}=0$, and then we
proceed to manipulate the resulting equation.

$\displaystyle\alpha_{*}(\alpha_{*}^{2}\gamma-\alpha_{*}(\gamma(\sigma^{2}+1)+1)+1)$
$\displaystyle=0$ (79)
$\displaystyle\alpha_{*}-\alpha_{*}^{2}(\gamma(\sigma^{2}+1)+1)$
$\displaystyle=-\alpha_{*}^{3}\gamma$ (80) $\displaystyle
2\alpha_{*}-\alpha_{*}^{2}(\gamma(\sigma^{2}+1)+1)$
$\displaystyle=\alpha_{*}(1-\alpha_{*}^{2}\gamma)$ (81)
$\displaystyle\frac{2\alpha_{*}-\alpha_{*}^{2}(\gamma(\sigma^{2}+1)+1)}{1-\alpha_{*}^{2}\gamma}$
$\displaystyle=\alpha_{*}.$ (82)

Continuing from this last equation,

$\displaystyle\alpha_{*}$
$\displaystyle=\frac{2\alpha_{*}-\alpha_{*}^{2}(\gamma(\sigma^{2}+1)+1)}{1-\alpha_{*}^{2}\gamma}$
(83)
$\displaystyle=\frac{2\alpha_{*}-\sigma^{2}\alpha_{*}^{2}\gamma-\alpha^{2}\gamma-\alpha^{2}+1-1}{1-\alpha_{*}^{2}\gamma}$
(84)
$\displaystyle=\frac{1-\alpha^{2}\gamma-(1-2\alpha_{*}+\alpha_{*}^{2})-\sigma^{2}\alpha_{*}^{2}\gamma}{1-\alpha_{*}^{2}\gamma}$
(85)
$\displaystyle=1-\frac{(1-\alpha_{*})^{2}+\sigma^{2}\alpha_{*}^{2}\gamma}{1-\alpha_{*}^{2}\gamma}$
(86) $\displaystyle=1-R_{\alpha_{*}}^{\mathrm{ens}}.$ (87)

Thus, if $\alpha_{*}$ is a root of (77) or $\alpha_{*}=0$, then
$R_{\alpha_{*}}^{\mathrm{ens}}=1-\alpha_{*}$. We proceed by checking the
larger root of (77), but before doing so, we derive the following equality:

$\displaystyle(\gamma(\sigma^{2}+1)+1)^{2}-4\gamma$
$\displaystyle=(\gamma(\sigma^{2}+1)+1)^{2}-(4\gamma^{2}(\sigma^{2}+1)+4\gamma)+4\gamma^{2}+4\sigma^{2}\gamma^{2}$
(88)
$\displaystyle=(\gamma(\sigma^{2}+1)+1-2\gamma)^{2}+4\sigma^{2}\gamma^{2}$
(89) $\displaystyle=(\gamma(\sigma^{2}-1)+1)^{2}+4\sigma^{2}\gamma^{2}.$ (90)

Now, we observe for the larger root (which we denote as $\alpha_{*}^{\prime}$)
that

$\displaystyle\alpha_{*}^{\prime}$
$\displaystyle=\frac{\gamma(\sigma^{2}+1)+1+\sqrt{(\gamma(\sigma^{2}-1)+1)^{2}+4\sigma^{2}\gamma^{2}}}{2\gamma}$
(91)
$\displaystyle\geq\frac{1}{2}\left(\sigma^{2}+1+\frac{1}{\gamma}+\left|\sigma^{2}-1+\frac{1}{\gamma}\right|\right)$
(92) $\displaystyle=\begin{cases}\sigma^{2}+\frac{1}{\gamma}&\text{if
}\frac{1}{\gamma}>1-\sigma^{2}\\\ 1&\text{if }\frac{1}{\gamma}\leq
1-\sigma^{2}.\end{cases}$ (93)

Thus the only case where $\alpha_{*}^{\prime}$ is a valid hyperparameter
choice (that is, $\alpha_{*}\leq\min\left\\{1,\gamma^{-1}\right\\}$) is when
$\sigma^{2}=0$ and $\gamma=1$, in which case $\alpha_{*}=1$ is a double root
of (77). So it suffices to evaluate the smaller root even in that case. Now
that we know that $\alpha_{*}^{\prime}$ is not conatined in
$[0,\min\left\\{1,\gamma^{-1}\right\\}]$ (except in the aforementioned special
case) and that by inspection of $R_{\alpha}^{\mathrm{ens}}$ the asymptote at
$\alpha=\gamma^{-1/2}$ is not contained in this interval, if we can show that
the smaller root (which we denote simply as $\alpha_{*}$) of (77) is contained
in this interval, then we know that it is the minimizer of
$R_{\alpha}^{\mathrm{ens}}$.

For the smaller root, it is clear from (78) that $\alpha_{*}\geq 0$. We show
by a series of equivalences that $\alpha_{*}\leq 1/\gamma$:

$\displaystyle\alpha_{*}=\frac{1}{2}\left(\sigma^{2}+1+\frac{1}{\gamma}-\sqrt{\left(\sigma^{2}-1+\frac{1}{\gamma}\right)^{2}+4\sigma^{2}}\right)\leq\frac{1}{\gamma}$
(94) $\displaystyle\iff{}$
$\displaystyle\sigma^{2}+1-\frac{1}{\gamma}\leq\sqrt{\left(\sigma^{2}-1+\frac{1}{\gamma}\right)^{2}+4\sigma^{2}}$
(95) $\displaystyle\iff{}$
$\displaystyle\left(\sigma^{2}+1-\frac{1}{\gamma}\right)^{2}\leq\left(\sigma^{2}-1+\frac{1}{\gamma}\right)^{2}+4\sigma^{2}$
(96) $\displaystyle\iff{}$
$\displaystyle\left(\sigma^{2}-1+\frac{1}{\gamma}\right)^{2}+4\sigma^{2}-4\frac{\sigma^{2}}{\gamma}\leq\left(\sigma^{2}-1+\frac{1}{\gamma}\right)^{2}+4\sigma^{2}$
(97) $\displaystyle\iff{}$
$\displaystyle\left(\sigma^{2}-1+\frac{1}{\gamma}\right)^{2}+4\sigma^{2}-4\frac{\sigma^{2}}{\gamma}\leq\left(\sigma^{2}-1+\frac{1}{\gamma}\right)^{2}+4\sigma^{2}$
(98) $\displaystyle\iff{}$ $\displaystyle 0\leq\frac{\sigma^{2}}{\gamma}.$
(99)

The last inequality is always true. Further, we note that every equivalence
here still holds under strict inequalities, so for $\sigma>0$, we have that
$\alpha_{*}<\gamma^{-1}$. By a similar argument, we can show that
$\alpha_{*}\leq 1$ and that $\alpha_{*}<1$ if and only if $\sigma>0$. By the
form of the derivative in (76), we know that $\alpha_{*}$, as the smaller
root, is a local minimum, and therefore it must be the minimum of
$R_{\alpha}^{\mathrm{ens}}$ on $[0,\min\left\\{1,\gamma^{-1}\right\\}]$.
Evaluating the risk at $\alpha_{*}$, we have

$\displaystyle R_{\alpha_{*}}^{\mathrm{ens}}$ $\displaystyle=1-\alpha_{*}$
(100)
$\displaystyle=1-\frac{1}{2}\left(\sigma^{2}+1+\frac{1}{\gamma}-\sqrt{\left(\sigma^{2}-1+\frac{1}{\gamma}\right)^{2}+4\sigma^{2}}\right)$
(101)
$\displaystyle=\frac{1}{2}\left(1-\sigma^{2}-\frac{1}{\gamma}+\sqrt{\left(\sigma^{2}-1+\frac{1}{\gamma}\right)^{2}+4\sigma^{2}}\right)$
(102)
$\displaystyle=\frac{1}{2}\left(\frac{\gamma-1}{\gamma}-\sigma^{2}+\sqrt{\left(\sigma^{2}-\frac{\gamma-1}{\gamma}\right)^{2}+4\sigma^{2}}\right)$
(103)
$\displaystyle=\inf_{\lambda}R(\widehat{{\bm{\beta}}}{}_{\lambda}^{\mathrm{ridge}}).$
(104)

##  Appendix E PROOFS OF DISCUSSION RESULTS

###  E.1 Proof of Equation (35) ($\mu$-scaled Risk)

Under the assumption that ${\bm{\Sigma}}=\mathbf{I}_{p}$, the $\mu$-scaled
risk is given by

$\displaystyle R(\mu\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{})$
$\displaystyle={\left\|{\bm{\beta}}-\mu\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{}\right\|}_{2}^{2}$
(105)
$\displaystyle={\left\|(1-\mu){\bm{\beta}}+\mu({\bm{\beta}}-\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{})\right\|}_{2}^{2}$
(106)
$\displaystyle=(1-\mu)^{2}{\left\|{\bm{\beta}}\right\|}_{2}^{2}+2(1-\mu)\mu\left\langle{\bm{\beta}},{\bm{\beta}}-\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{}\right\rangle+\mu^{2}{\left\|{\bm{\beta}}-\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{}\right\|}_{2}^{2}$
(107)

Examining the inner product, we find that

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathbf{z},\mathcal{S},\mathcal{T}}\left[\left\langle{\bm{\beta}},{\bm{\beta}}-\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{}\right\rangle\right]$
$\displaystyle=\left\langle{\bm{\beta}},{\mathbb{E}}_{\mathbf{X},\mathbf{z},\mathcal{S},\mathcal{T}}\left[{\bm{\beta}}-\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{}\right]\right\rangle$
(108)
$\displaystyle=\left\langle{\bm{\beta}},{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathbf{I}_{p}-\frac{1}{k}\sum_{i=1}^{k}\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\right]{\bm{\beta}}\right\rangle,$
(109)

where the equation (109) holds because ${\mathbb{E}}[\mathbf{z}]={\bm{0}}$.
Because the subsamplings are identically distributed, we have

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathbf{I}_{p}-\frac{1}{k}\sum_{i=1}^{k}\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\right]$
$\displaystyle=\mathbf{I}_{p}-{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\right]$
(110)
$\displaystyle=\mathbf{I}_{p}-{\mathbb{E}}_{\mathbf{X},\mathcal{S},\mathcal{T}}\left[\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\left(\mathbf{S}_{i}\mathbf{S}_{i}^{\top}+\mathbf{S}_{i}^{c}{\mathbf{S}_{i}^{c}}^{\top}\right)\right]$
(111)
$\displaystyle=\mathbf{I}_{p}-{\mathbb{E}}_{\mathcal{S}}\left[\mathbf{S}_{i}\left(\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\right)^{\dagger}\mathbf{T}_{i}^{\top}\mathbf{X}\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\right]$
(112)
$\displaystyle=\mathbf{I}_{p}-{\mathbb{E}}_{\mathcal{S}}\left[\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\right]$
(113) $\displaystyle=(1-\alpha)\mathbf{I}_{p},$ (114)

where the equation (112) holds because
${\mathbb{E}}[\mathbf{X}\mathbf{S}_{i}^{c}]={\bm{0}}$. Thus

$\displaystyle{\mathbb{E}}_{\mathbf{X},\mathbf{z},\mathcal{S},\mathcal{T}}[R(\mu\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{})]$
$\displaystyle=(1-\mu)^{2}{\left\|{\bm{\beta}}\right\|}_{2}^{2}+2(1-\mu)\mu{\mathbb{E}}_{\mathbf{X},\mathbf{z},\mathcal{S},\mathcal{T}}\left[\left\langle{\bm{\beta}},{\bm{\beta}}-\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{}\right\rangle\right]+\mu^{2}{\mathbb{E}}_{\mathbf{X},\mathbf{z},\mathcal{S},\mathcal{T}}\left[{\left\|{\bm{\beta}}-\widehat{{\bm{\beta}}}{}^{\mathrm{ens}}{}\right\|}_{2}^{2}\right]$
(115)
$\displaystyle=(1-\mu)^{2}+2(1-\mu)\mu(1-\alpha)+\mu^{2}R_{\alpha}^{\mathrm{ens}},$
(116)

where the last equality holds because
$\langle{\bm{\beta}},{\bm{\beta}}\rangle={\left\|{\bm{\beta}}\right\|}_{2}^{2}=1$.

###  E.2 Proof of Equation (38) (Generalized Dropout)

For $k\to\infty$, dropout minimizes the expected loss:

$\displaystyle{\mathbb{E}}_{S_{i}}\left[\ell_{i}({\bm{\beta}}^{\prime})\right]$
$\displaystyle={\mathbb{E}}_{S_{i}}\left[{\left\|\mathbf{X}\mathbf{S}_{i}\mathbf{S}_{i}^{\top}{\bm{\beta}}^{\prime}-\mathbf{y}\right\|}_{2}^{2}\right].$
(117)

The expected loss is convex in ${\bm{\beta}}^{\prime}$, so we can find its
minimizer by the first order optimality condition:

$\displaystyle\nabla_{{\bm{\beta}}^{\prime}}{\mathbb{E}}_{S_{i}}\left[\ell_{i}({\bm{\beta}}^{\prime})\right]={\mathbb{E}}_{S_{i}}\left[\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\left(\mathbf{X}\mathbf{S}_{i}\mathbf{S}_{i}^{\top}{\bm{\beta}}^{\prime}-\mathbf{y}\right)\right]=0$
(118)

Thus,

$\displaystyle\widehat{{\bm{\beta}}}{}=\left({\mathbb{E}}_{S_{i}}\left[\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\mathbf{X}\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\right]\right)^{-1}{\mathbb{E}}_{S_{i}}\left[\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\right]\mathbf{X}^{\top}\mathbf{y}.$
(119)

Turning first to the inverse, consider that

$\displaystyle\left[{\mathbb{E}}_{S_{i}}\left[\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\mathbf{X}\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\right]\right]_{j\ell}$
$\displaystyle=\left[\mathbf{X}^{\top}\mathbf{X}\right]_{j\ell}\Pr(j\in
S_{i},\ell\in S_{i}),$ (120)

and that

$\displaystyle\Pr(j\in S_{i},\ell\in S_{i})$
$\displaystyle=\begin{cases}\alpha_{j}&\text{if }j=\ell,\\\
\alpha_{j}\alpha_{\ell}&\text{otherwise}.\end{cases}$ (121)

This gives us

$\displaystyle{\mathbb{E}}_{S_{i}}\left[\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\mathbf{X}^{\top}\mathbf{X}\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\right]$
$\displaystyle=\mathbf{A}\mathbf{X}^{\top}\mathbf{X}\mathbf{A}+\mathrm{diag}(\mathbf{X}^{\top}\mathbf{X})(\mathbf{A}-\mathbf{A}^{2}),$
(122)

where $\mathbf{A}=\mathrm{diag}({\bm{\alpha}})$. By a similar and simpler
argument,

$\displaystyle{\mathbb{E}}_{S_{i}}\left[\mathbf{S}_{i}\mathbf{S}_{i}^{\top}\right]=\mathbf{A},$
(123)

which all together yields

$\displaystyle\widehat{{\bm{\beta}}}{}$
$\displaystyle=\left(\mathbf{A}\mathbf{X}^{\top}\mathbf{X}\mathbf{A}+\mathrm{diag}(\mathbf{X}^{\top}\mathbf{X})(\mathbf{I}_{p}-\mathbf{A})\mathbf{A}\right)^{-1}\mathbf{A}\mathbf{X}^{\top}\mathbf{y}$
(124)
$\displaystyle=\mathbf{A}^{-1}\left(\mathbf{X}^{\top}\mathbf{X}+\mathbf{A}^{-1}\mathrm{diag}(\mathbf{X}^{\top}\mathbf{X})(\mathbf{I}_{p}-\mathbf{A})\right)^{-1}\mathbf{X}^{\top}\mathbf{y}.$
(125)

##  Appendix F PROOFS OF LEMMAS A.1 and A.2

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

Without loss of generality, let
$\begin{bmatrix}\mathbf{X}_{1}\;\mathbf{X}_{2}\end{bmatrix}=\mathbf{X}$, such
that $\mathbf{X}_{2}=\mathbf{X}\mathbf{S}$. Let
$\begin{bmatrix}\mathbf{Y}_{1}\\\
\mathbf{Y}_{2}\end{bmatrix}=\mathbf{X}^{\dagger}$ be a partitioning of the
pseudo-inverse of $\mathbf{X}$ in the same manner, such that
$\mathbf{Y}_{2}=\mathbf{S}^{\top}\mathbf{X}^{\dagger}=\mathbf{S}^{\top}\left(\mathbf{X}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}$.
Then the Gram matrix can be written as

$\displaystyle\mathbf{X}^{\top}\mathbf{X}=\begin{bmatrix}\mathbf{X}_{1}^{\top}\mathbf{X}_{1}&\mathbf{X}_{1}^{\top}\mathbf{X}_{2}\\\
\mathbf{X}_{2}^{\top}\mathbf{X}_{1}&\mathbf{X}_{2}^{\top}\mathbf{X}_{2}\end{bmatrix},$
(126)

and using block matrix inversion, the inverse admits the form
$\left(\mathbf{X}^{\top}\mathbf{X}\right)^{-1}=\begin{bmatrix}\mathbf{A}&\mathbf{B}\\\
\mathbf{C}&\mathbf{D}\end{bmatrix}$. The relevant quantities are

$\displaystyle\mathbf{C}$
$\displaystyle=-\mathbf{D}\mathbf{X}_{2}^{\top}\mathbf{X}_{1}\left(\mathbf{X}_{1}^{\top}\mathbf{X}_{1}\right)^{-1}$
(127) $\displaystyle\mathbf{D}$
$\displaystyle=\left(\mathbf{X}_{2}^{\top}\mathbf{X}_{2}-\mathbf{X}_{2}^{\top}\mathbf{X}_{1}\left(\mathbf{X}_{1}^{\top}\mathbf{X}_{1}\right)^{-1}\mathbf{X}_{1}^{\top}\mathbf{X}_{2}\right)^{-1}$
(128)
$\displaystyle=\left(\mathbf{X}_{2}^{\top}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{1}^{\top})}\mathbf{X}_{2}\right)^{-1},$
(129)

Where
${\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{1}^{\top})}\triangleq\mathbf{I}_{n}-\left(\mathbf{X}_{1}^{\top}\right)^{\dagger}\mathbf{X}_{1}^{\top}$
denotes the projection onto the column space of $\mathbf{X}_{1}$. This gives

$\displaystyle\mathbf{Y}_{2}$
$\displaystyle=\mathbf{C}\mathbf{X}_{1}^{\top}+\mathbf{D}\mathbf{X}_{2}^{\top}$
(130)
$\displaystyle=\mathbf{D}\mathbf{X}_{2}^{\top}\left(\mathbf{I}_{n}-\mathbf{X}_{1}\left(\mathbf{X}_{1}^{\top}\mathbf{X}_{1}\right)^{-1}\mathbf{X}_{1}^{\top}\right)$
(131)
$\displaystyle=\mathbf{D}\mathbf{X}_{2}^{\top}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{1}^{\top})}\left(\mathbf{X}_{2}\mathbf{X}_{2}^{\dagger}+{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{2}^{\top})}\right)$
(132)
$\displaystyle=\mathbf{X}_{2}^{\dagger}+\mathbf{D}\mathbf{X}_{2}^{\top}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{1}^{\top})}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{2}^{\top})}.$
(133)

Let $\mathbf{U}$, $\mathbf{U}_{*}$, and $\mathbf{V}$ be the matrices
containing the left singular vectors of $\mathbf{X}_{2}$,
${\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{2}^{\top})}$, and
${\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{1}^{\top})}$, respectively. Because the
rows of $\mathbf{X}$ are independently drawn from a spherical Gaussian
distribution, $\mathbf{V}$ has a uniform distribution over orthogonal matrices
in $\mathbb{R}^{n\times|S^{c}|}$. As such,
${\mathbb{E}}_{\mathbf{V}}\left[\mathbf{V}^{\top}\mathbf{U}_{*}|\mathbf{V}^{\top}\mathbf{U}\right]=\bm{0}$.
Then

$\displaystyle{\mathbb{E}}_{\mathbf{X}_{1}}\left[\mathbf{D}\mathbf{X}_{2}^{\top}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{1}^{\top})}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{2}^{\top})}\right]$
$\displaystyle={\mathbb{E}}_{\mathbf{V}}\left[\left(\mathbf{X}_{2}^{\top}\mathbf{V}\mathbf{V}^{\top}\mathbf{X}_{2}\right)^{-1}\mathbf{X}_{2}^{\top}\mathbf{V}\mathbf{V}^{\top}\mathbf{U}_{*}\mathbf{U}_{*}^{\top}\right]$
(134)
$\displaystyle={\mathbb{E}}_{\mathbf{V}}\left[\left(\mathbf{X}_{2}^{\top}\mathbf{V}\mathbf{V}^{\top}\mathbf{X}_{2}\right)^{-1}\mathbf{X}_{2}^{\top}\mathbf{V}{\mathbb{E}}_{\mathbf{V}}\left[\mathbf{V}^{\top}\mathbf{U}_{*}\Big{|}\mathbf{V}^{\top}\mathbf{U}\right]\mathbf{U}_{*}^{\top}\right]$
(135) $\displaystyle=\bm{0},$ (136)

which combined with (133) yields the first claim.

For the second claim, let $\mathbf{V}_{*}$ denote the left singular vectors of
$\mathbf{X}_{1}$, and observe that
${\mathbb{E}}_{\mathbf{V}}\left[\mathbf{V}^{\top}\mathbf{U}_{*}|\mathbf{V}^{\top}\mathbf{U},\mathbf{V}_{*}\right]=\bm{0}$.
Then using similar arguments,

$\displaystyle{\mathbb{E}}_{\mathbf{X}_{1}}\left[\mathbf{X}_{1}^{\top}f(\mathbf{X}_{2})\mathbf{S}^{\top}\mathbf{X}^{\dagger}\right]$
$\displaystyle={\mathbb{E}}_{\mathbf{X}_{1}}\left[\mathbf{X}_{1}^{\top}f(\mathbf{X}_{2})\mathbf{S}^{\top}\left(\mathbf{X}_{2}^{\dagger}+\mathbf{D}\mathbf{X}_{2}^{\top}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{1}^{\top})}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}_{2}^{\top})}\right)\right]$
(137)
$\displaystyle={\mathbb{E}}_{\mathbf{X}_{1}}\left[\mathbf{X}_{1}^{\top}f(\mathbf{X}_{2})\mathbf{S}^{\top}\left(\mathbf{X}_{2}^{\dagger}+\left(\mathbf{X}_{2}^{\top}\mathbf{V}\mathbf{V}^{\top}\mathbf{X}_{2}\right)^{-1}\mathbf{X}_{2}^{\top}\mathbf{V}{\mathbb{E}}_{\mathbf{V}}\left[\mathbf{V}^{\top}\mathbf{U}_{*}\Big{|}\mathbf{V}^{\top}\mathbf{U},\mathbf{V}_{*}\right]\mathbf{U}_{*}^{\top}\right)\right]$
(138)
$\displaystyle={\mathbb{E}}_{\mathbf{X}_{1}}\left[\mathbf{X}_{1}^{\top}f(\mathbf{X}_{2})\mathbf{S}^{\top}\mathbf{X}_{2}^{\dagger}\right]$
(139) $\displaystyle=\bm{0}.$ (140)

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

Define
${\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}^{\top})}\triangleq\mathbf{I}_{n}-\left(\mathbf{X}^{\top}\right)^{\dagger}\mathbf{X}^{\top}$,
the projection operator onto the null space of $\mathbf{X}^{\top}$. Then for
the first claim,

$\displaystyle{\mathbb{E}}_{T_{1},T_{2}}\left[\left(\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{\dagger}\mathbf{T}_{1}^{\top}\left(\left(\mathbf{T}_{2}^{\top}\mathbf{X}\right)^{\dagger}\mathbf{T}_{2}^{\top}\right)^{\top}\right]$
$\displaystyle={\mathbb{E}}_{T_{1},T_{2}}\left[\left(\mathbf{X}^{\top}\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{X}\left(\mathbf{X}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{X}\right)^{-1}\right]$
(141)
$\displaystyle={\mathbb{E}}_{T_{1},T_{2}}\left[\left(\mathbf{X}^{\top}\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\left(\mathbf{X}\left(\mathbf{X}\mathbf{X}^{\top}\right)^{\dagger}\mathbf{X}^{\top}+{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}^{\top})}\right)\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{X}\left(\mathbf{X}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{X}\right)^{-1}\right]$
(142)
$\displaystyle=\left(\mathbf{X}\mathbf{X}^{\top}\right)^{\dagger}+{\mathbb{E}}_{T_{1},T_{2}}\left[\left(\mathbf{X}^{\top}\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}\mathbf{T}_{1}\mathbf{T}_{1}^{\top}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}^{\top})}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{X}\left(\mathbf{X}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{X}\right)^{-1}\right]$
(143)
$\displaystyle=\left(\mathbf{X}\mathbf{X}^{\top}\right)^{\dagger}+\frac{|T_{1}||T_{2}|}{n^{2}}\left(\mathbf{X}^{\top}\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}^{\top})}\mathbf{X}\left(\mathbf{X}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{X}\right)^{-1}$
(144) $\displaystyle=\left(\mathbf{X}^{\top}\mathbf{X}\right)^{\dagger}.$
(145)

The equality (144) follows due the fact that, because of the distributional
assumption on the rows of $\mathbf{X}$,
$\mathbf{X}^{\top}\mathbf{T}_{j}\mathbf{T}_{j}^{\top}\mathbf{X}$ and
$\mathbf{T}_{j}\mathbf{T}_{j}^{\top}$ are conditionally independent given
$|T_{j}|$. The equality (145) follows because
$\mathbf{X}^{\top}{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}^{\top})}=\bf{0}$.

For the second claim,

$\displaystyle{\mathbb{E}}_{T_{1},T_{2}}\left[\left(\left(\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{\dagger}\mathbf{T}_{1}^{\top}\right)^{\top}\mathbf{A}\left(\mathbf{T}_{2}^{\top}\mathbf{Y}\right)^{\dagger}\mathbf{T}_{2}^{\top}\right]$
$\displaystyle={\mathbb{E}}_{T_{1},T_{2}}\left[\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\mathbf{X}\left(\mathbf{X}^{\top}\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{-1}\mathbf{A}\left(\mathbf{Y}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{Y}\right)^{-1}\mathbf{Y}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\right]$
(146)
$\displaystyle={\mathbb{E}}_{T_{1},T_{2}}\left[\left(\left(\mathbf{X}^{\top}\right)^{\dagger}\mathbf{X}^{\top}+{\bm{\Pi}}_{\mathrm{Null}(\mathbf{X}^{\top})}\right)\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\mathbf{X}\left(\mathbf{X}^{\top}\mathbf{T}_{1}\mathbf{T}_{1}^{\top}\mathbf{X}\right)^{-1}\mathbf{A}\left(\mathbf{Y}\top\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{Y}\right)^{-1}\mathbf{Y}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\right]$
(147)
$\displaystyle={\mathbb{E}}_{\Pi(T_{2})}\left[\left(\mathbf{X}^{\dagger}\right)^{\top}\mathbf{A}\left(\mathbf{Y}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{Y}\right)^{-1}\mathbf{Y}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\right]$
(148)
$\displaystyle={\mathbb{E}}_{\Pi(T_{2})}\left[\left(\mathbf{X}^{\dagger}\right)^{\top}\mathbf{A}\left(\mathbf{Y}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\mathbf{Y}\right)^{-1}\mathbf{Y}^{\top}\mathbf{T}_{2}\mathbf{T}_{2}^{\top}\left(\mathbf{Y}\mathbf{Y}^{\dagger}+{\bm{\Pi}}_{\mathrm{Null}(\mathbf{Y}^{\top})}\right)\right]$
(149)
$\displaystyle=\left(\mathbf{X}^{\dagger}\right)^{\top}\mathbf{A}\mathbf{Y}^{\dagger},$
(150)

where the equations (148) and (150) follow by similar arguments to those used
to show the first claim.