Theoretical Limitations of Ensembles in the Age of Overparameterization

## Abstract

Abstract Classic tree-based ensembles generalize better than any single decision tree. In contrast, recent empirical
studies find that modern ensembles of (overparameterized) neural networks may
not provide any inherent generalization advantage over single but larger neural networks. This paper
clarifies how modern overparameterized ensembles differ from their classic
underparameterized counterparts, using ensembles of random feature (RF)
regressors as a basis for developing theory. In contrast to the
underparameterized regime, where ensembling typically induces regularization and
increases generalization, we prove that infinite ensembles of overparameterized RF
regressors become pointwise equivalent to (single) infinite-width RF regressors. This
equivalence, which is exact for ridgeless models and approximate for small
ridge penalties, implies that overparameterized ensembles and single large
models exhibit nearly identical generalization. As a consequence, we can characterize
the predictive variance amongst ensemble members, and demonstrate that it
quantifies the expected effects of increasing capacity rather than capturing
any conventional notion of uncertainty. Our results challenge common
assumptions about the advantages of ensembles in overparameterized settings,
prompting a reconsideration of how well intuitions from underparameterized
ensembles transfer to deep ensembles and the overparameterized regime.

## 1 Introduction

Ensembling is one of the most well-established techniques in machine learning
(e.g. Schapire, ; Hansen & Salamon, ; Opitz & Maclin, ; Dietterich, ).
Historically, most ensembles aggregated component models that are simple by today’s standards.
Common techniques like bagging (Breiman, ), feature selection
(Breiman, ), random projections (Kabán, ; Thanei et al., ), and boosting (Freund, ; Chen & Guestrin, )
were developed and analyzed assuming decision trees, least-squares regressors,
and other *underparameterized* component models incapable of achieving near-zero training error.
Crucially, the resulting ensembles achieve better generalization
than what could be achieved by any individual component model.

Recently, researchers and practitioners have turned to ensembling
large *overparameterized* models, such as neural networks,
which have more than enough capacity to memorize training datasets
and are typically trained with little to no regularization.
Like ensembles of underparameterized models, ensembles of large neural networks
are often used to reduce generalization error (Lee et al., ; Fort et al., ).
Motivated by practical effectiveness and heuristics from classic ensembles (Mentch & Hooker, ), some have further argued that the predictive variance amongst component models in these so-called *deep ensembles*
is a well-calibrated notion of uncertainty (Lakshminarayanan et al., ; Ovadia et al., ; Gustafsson et al., )
that can be used on downstream decision-making tasks (Gal et al., ; Yu et al., ).

While there are few theoretical works analyzing these modern overparameterized ensembles,
recent empirical evidence suggests that
intuitions from their underparameterized counterparts do not hold in this new regime.
For example,
classic methods to increase diversity amongst component models, such as bagging,
are harmful for deep ensembles (Nixon et al., ; Jeffares et al., ; Abe et al., )
despite being nearly universally beneficial for underparameterized ensembles.
Moreover, several recent studies question
whether deep ensembles offer significant improvements in robustness and
uncertainty quantification over what can be achieved by a single (but larger) neural network
(Abe et al., ; Theisen et al., ; Chen et al., ).
These results suggest that an ensemble of (large) overparameterized networks may
not differ fundamentally from a single (extremely large) neural network,
in contrast to the underparameterized regime where ensembles
are a fundamentally different class of predictors Schapire (); Breiman (); Kabán ().

To address this divergence and verify recent empirical findings on deep ensembles,
we develop a theoretical characterization of ensembles in the
overparameterized regime, with the goal of contrasting against
(traditional) underparameterized ensembles.
We answer the following questions:

- 1.
Do large ensembles of overparameterized models differ from single (very large)
models trained on the same data?
Does the capacity of the component models affect this difference?
- 2.
Under a fixed parameter/computation budget, does an ensemble of overparameterized
models provide additional generalization or robustness benefits over a single (larger) model?
- 3.
What does the predictive variance of overparameterized ensembles measure, and does it relate to different notions of uncertainty?

To answer these questions, we analyze ensembles of overparameterized random feature (RF) linear regressors,
a widely used theoretically-tractable approximation of neural networks (e.g. Belkin et al., ; Bartlett et al., ; Mei & Montanari, ).
These models can be interpreted as neural networks where only the last layer is trained (e.g. Rudi & Rosasco, ; Belkin et al., )
or as first-order Taylor approximations of neural networks (e.g. Jacot et al., ).
By averaging models that differ solely in their random features, we emulate the common practice of ensembling neural networks that differ only by random initialization (Lakshminarayanan et al., ).
Our analysis focuses on the practically relevant regime where regressors are trained with little to no regularization.

### 1.1 Related Work

#### Random feature models.

RF models perform regression on a random subset or projection of
a high- (or infinite-) dimensional feature representation.
Originally introduced as a scalable approximation to kernel machines (Rahimi & Recht, ; ; ),
RF regressors have seen growing theoretical interest as simplified models of neural networks
(e.g. Belkin et al., ; Jacot et al., ; Bartlett et al., ; Mei & Montanari, ; Simon et al., ).
This approximation of neural networks becomes exact in the limit of infinite width (e.g. Jacot et al., ; Lee et al., ).

#### Underparameterized random feature models and ensembles.

There are many works theoretically characterizing ensembles of tree-based models (e.g. Schapire & Singer, ; Sexton & Laake, ; Wager et al., ; Mentch & Hooker, ).
Here, we restrict our discussion to analyses of (ensembles of) RF regressors.
Most works of this nature analyze *underparameterized* models,
where the number of random features (i.e., the width) is assumed to be far fewer than the number of data points.
In the underparameterized fixed design setting,
the infinite ensemble of unregularized RF regressors
achieves the same generalization error as ridge regression on the original (unprojected) inputs
(Kabán, ; Thanei et al., ; Bach, ).
We emphasize the distinction between underparameterized component models and their aggregated prediction:
i.e., the ensemble of unregularized regressors is equivalent to a regularized predictor.
(We provide theoretical analysis in Appx. D that further demonstrates
ridge-like behaviour of underpameterized RF ensembles.)

#### Overparameterized random feature models.

Recent works on RF models have focused on the *overparameterized regime*,
often using high-dimensional asymptotics to characterize generalization error
(Adlam & Pennington, ; Hastie et al., ; Mei & Montanari, ; Bach, ).
Implicit in many works is an assumption of *Gaussian universality*,
in which the marginal distributions over the random features are replaced by moment-matched Gaussians.
While such assumptions are common throughout asymptotic random matrix theory (e.g. Tao, ),
our work aims to establish more general results that hold for more general random features.
We demonstrate—both theoretically and empirically—that Gaussianity may be an inappropriate approximation for neural network features
when comparing the pointwise behaviour of ensembles versus single models.
Most related to our work is Jacot et al. (), who characterize the pointwise
expectation and variance of ridge-regularized RF models where the features are draws from Gaussian processes (GP).
Our results for ensembles similarly characterize the pointwise expectation and variance of overparameterized RF models and ensembles;
however, we significantly weaken the assumptions on the random features.

### 1.2 Contributions

Figure: Figure 1: An infinite ensemble of overparameterized RF models is equivalent to a single infinite-width RF model. (Left) We show a sample of $100$ finite-width RF models (blue) with ReLU activations trained on the same $N=6$ data points. Additionally, we show the single infinite-width RF model (pink). The finite-width predictions concentrate around the infinite-width model. (Right) We again show the single infinite-width RF model (pink) and the “infinite” ensemble of $M=10,000$ RF models (blue). We note no perceptible difference between the two.
Refer to caption: x1.png

We consider ensembles of *overparameterized* RF regressors in both
the ridgeless and small ridge regimes.
Unlike prior work, we make minimal assumptions about the distribution of the random features.
Therefore, our results can be assumed to hold for most RF ensembles
rather than only those that are compositions of GP-random features.
Concretely, we make the following contributions:

To answer Question 1: we show that the average ridgeless RF regressor is pointwise
equivalent to its corresponding ridgeless kernel regressor
(LABEL:thm:identity-infinite-model-infinite-ensemble-subexp),
implying that an infinite ensemble of overparameterized RF models
is *exactly* equivalent to a single infinite-width RF model.
We further show that this equivalence approximately holds in the small ridge regime
(LABEL:thm:smoothness-difference-ridgeless-ridge).

To answer Question 2: we use rates established in prior work to
demonstrate that the variance reduction from ensembling overparameterized
RF regressors is very similar to increasing the number of features in a
single model. This shows that ensembles do not offer additional
generalization or robustness advantages over single models under fixed
parameter budgets (see Sec. 3.2).

To answer Question 3: we show that
the predictive variance in an overparameterized ensemble is
the expected squared difference between the predictions from a
(finite-width) RF regressor and its corresponding kernel regressor (i.e., the infinite-width model).
With this finding, we demonstrate that ensemble variance differs from conventional uncertainty quantifications,
except in practically unrealistic cases where the random features are sampled from
a GP (see Sec. 3.2).

Altogether, these results support recent empirical findings
that deep ensembles offer few generalization and uncertainty quantification benefits over single models
(Abe et al., ; Theisen et al., ).
Our theory and experiments demonstrate that these phenomena are not specific to neural networks
but are more general properties of ensembles in the overparameterized regime.

## 2 Setup

We work in a regression setting with a fixed training dataset $\mathcal{D}=\{(x_{i},y_{i})\}_{i=1}^{N}\in({\mathcal{X}}\times\mathbb{R})^{N}$ of size $N$. The vector $y\in\mathbb{R}^{N}$ represents the concatenation of all training responses.

We consider *random feature (RF) models* of the form:

$$ $h_{\mathcal{W}}(x)=\tfrac{1}{\sqrt{D}}\operatorname{{\textstyle\sum}}_{j=1}^{D }\phi(\omega_{j},x)\theta_{j},$ $$

where $\theta_{j}$ are learned parameters, ${\mathcal{W}}=\{\omega_{j}\}_{j=1}^{D}\in\Omega^{D}$ are i.i.d. samples from a distribution $\pi(\cdot)$, and $\phi:\Omega\times{\mathcal{X}}\to\mathbb{R}$ is the *feature extraction function*. For a ReLU-based RF model with $p$-dimensional inputs, we have ${\mathcal{X}}=\Omega=\mathbb{R}^{p}$ and $\phi(\omega_{j},x)=\max(0,\omega_{j}^{\top}x)$. Although RF models cannot fully explain the behavior of neural networks (e.g. Ghorbani et al., ; Li et al., ; Pleiss & Cunningham, ), they serve as useful proxies for studying the effects of overparameterization and model capacity on generalization (e.g. Belkin et al., ; Adlam & Pennington, ; Mallinar et al., ).

#### Notation.

For any $x,x^{\prime}\in{\mathcal{X}}$, let $k(x,x^{\prime})=\mathbb{E}_{\omega}[\phi(\omega,x)\phi(\omega,x^{\prime})]$
denote the second moment of the feature extraction function $\phi(\omega,\cdot)$.
We note that the function $k$ is a positive definite kernel function,
and we will refer to it as such.
We will use the matrices
$K\coloneqq[k(x_{i},x_{j})]_{ij}\in\mathbb{R}^{N\times N}$
and
$\Phi_{\mathcal{W}}\coloneqq[\phi(\omega_{j},x_{i})]_{ij}\in\mathbb{R}^{N\times
D}$
to denote the kernel function applied to all training data pairs and the
feature extraction function applied to all data/feature combinations, respectively.
We will drop the subscript ${\mathcal{W}}$ when the set of random features is clear from context.
We assume that $K$ is invertible.

Throughout our analysis, it will be useful to consider the “whitened” feature matrix $W=R^{-\top}\Phi\in\mathbb{R}^{N\times D}$
where $R^{\top}R=K$ is the Cholesky decomposition of the kernel matrix $K$.
When considering a test point $x^{*}\in{\mathcal{X}}$ (or equivalently a set of test points), we extend the $K$, $R$, $\Phi$, $W$ notation by

$$ $\begin{bmatrix}K&[k(x_{i},x^{*})]_{i}\\ [k(x^{*}\!,x_{j})]_{j}&k(x^{*}\!,x^{*})\end{bmatrix}\!=\!\begin{bmatrix}R&c\\ 0&r_{\perp}\end{bmatrix}^{\top}\!\!\begin{bmatrix}R&c\\ 0&r_{\perp}\end{bmatrix},\quad\begin{bmatrix}W\\ w_{\perp}^{\top}\end{bmatrix}=\begin{bmatrix}R&c\\ 0&r_{\perp}\end{bmatrix}^{-\top}\!\!\begin{bmatrix}\Phi\\ [\phi(\omega_{j},x^{*})]_{j}\end{bmatrix}.$ (1) $$

For fixed training/test points,
$\mathbb{E}_{W}[WW^{\top}]=D\cdot I$, $\mathbb{E}_{w_{\perp}}[w_{\perp}^{\top}w_{\perp}]=D$ and $\mathbb{E}_{W,w_{\perp}}[w_{\perp}^{\top}W^{\top}]=0$ which can be directly derived from $\mathbb{E}_{\Phi}[\Phi\Phi^{\top}]=D\cdot K$ (and similar properties for $\phi^{*}$).
Moreover, the columns $[w_{i};w_{\bot i}]$ of $[W;w_{\bot}]$ are i.i.d. since they are transformations of the i.i.d. columns of $\Phi$.

#### Overparameterized ridge/ridgeless regressors and ensembles.

As our focus is the overparameterized regime,
we assume a computational budget of $D>N$ features
(${\mathcal{W}}=\{\omega_{1},\ldots,\omega_{D}\}\sim\pi^{D}$)
to construct an RF regressor $h_{\mathcal{W}}(x)=\tfrac{1}{\sqrt{D}}\phi_{\mathcal{W}}(x)^{\top}\theta$.
We train the regressor parameters $\theta$ to minimize the loss $\|\frac{1}{\sqrt{D}}\Phi_{\mathcal{W}}\theta-y\|_{2}^{2}+\lambda\|\theta\|_{2}
^{2}$ for some ridge parameter $\lambda\geq 0$.
When $\lambda>0$, this optimization problem admits the closed-form solution
$\theta_{{\mathcal{W}},\lambda}^{(\mathrm{RR})}=\tfrac{1}{\sqrt{D}}\Phi_{
\mathcal{W}}^{\top}\left(\tfrac{1}{D}\cdot\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}
^{\top}+\lambda I\right)^{-1}y.$
Although the learning problem is underspecified when $\lambda=0$ (i.e., in the ridgeless case),
the implicit bias of (stochastic) gradient descent initialized at zero leads to the minimum norm interpolating solution
$\theta^{(\mathrm{LN})}_{\mathcal{W}}=\tfrac{1}{\sqrt{D}}(\Phi)^{\top}\left(
\tfrac{1}{D}\cdot\Phi\Phi^{\top}\right)^{-1}y.$
We denote the resulting ridge(less) regressors as

$$ $h^{(\mathrm{LN})}_{\mathcal{W}}(\cdot):=\tfrac{1}{\sqrt{D}}\begin{bmatrix}\phi (\omega_{j},\cdot)\end{bmatrix}_{j}\theta^{(\mathrm{LN})}_{\mathcal{W}},\qquad h ^{(\mathrm{RR})}_{{\mathcal{W}},\lambda}(\cdot):=\tfrac{1}{\sqrt{D}}\begin{ bmatrix}\phi(\omega_{j},\cdot)\end{bmatrix}_{j}\theta_{{\mathcal{W}},\lambda}^ {(\mathrm{RR})}.$ $$

We also consider ensembles of $M$ ridge(less) regressors.
We assume that each is trained on the training same dataset but using different sets of i.i.d. random features ${\mathcal{W}}_{1},\ldots,{\mathcal{W}}_{M}\sim\pi^{D}$, with $D>N$.
Thus, the only source of randomness in these ensembles comes from the random selection of features ${\mathcal{W}}_{i}$,
analogous to the standard training procedure of deep ensembles (Lakshminarayanan et al., ).
The ensemble prediction is given by the arithmetic average of the individual models

$$ $\bar{h}_{{\mathcal{W}}_{1:M}}(\cdot)=\tfrac{1}{M}\operatorname{{\textstyle\sum }}_{m=1}^{M}h_{{\mathcal{W}}_{m},\lambda}(\cdot)=\tfrac{1}{M}\operatorname{{ \textstyle\sum}}_{m=1}^{M}\left[\tfrac{1}{D}\begin{bmatrix}\phi(\omega_{mj}, \cdot)\end{bmatrix}_{j}\Phi_{{\mathcal{W}}_{m}}^{\top}\left(\tfrac{1}{D}\cdot \Phi_{{\mathcal{W}}_{m}}\Phi_{{\mathcal{W}}_{m}}^{\top}+\lambda I\right)^{-1} \right]y.$ $$

#### Assumptions.

A key difference between this paper and prior literature is the set of assumptions about the random feature distribution $\pi(\cdot)$.
It is commonly assumed that entries in the extended whitened feature matrix $[W;w_{\bot}]$ are i.i.d. draws from a zero-mean sub-Gaussian distribution (e.g. Bartlett et al., ; Bach, ),
which implicitly places constraints on $\phi(\cdot,\cdot)$ and $\pi(\cdot)$.
Many works further assume *Gaussian universality*—i.e. that the distribution of $W,w_{\bot}$ can be modeled by i.i.d. standard Gaussian random variables (Adlam & Pennington, ; Mei & Montanari, ; Simon et al., )—implying that $\phi(\omega_{i},\cdot)$ are draws from a Gaussian process with covariance $k$.(^1^11 Assume the entries of $W$ and $w_{\perp}$ are i.i.d. Gaussian. Then the $i^{\mathrm{th}}$ feature applied to training/test inputs ($[R^{\top}w_{i};c^{\top}w_{i}\!+\!r_{\perp}w_{\perp i}]$) is multivariate Gaussian. This fact holds for any train/test data; thus, the $i^{\mathrm{th}}$ feature is a GP by definition (e.g. Rasmussen & Williams, , Ch. 2).)
We argue this assumption is unrealistic when considering features that resemble those from neural networks.
For example, ReLU features are always non-negative; thus, the mean of $W=R^{-\top}\Phi$ is almost surely non-zero.
Moreover, if $\mathcal{X}\subseteq\mathbb{R}^{p}$ with $p<N$, then feature extraction functions of the form $\sigma(\omega^{\top}x)$
are fully specified by a $p$-dimensional random variable.
Thus, knowing $N$ evaluations of $\omega_{j}^{\top}x_{i}$ allows one to
infer $\omega_{j}$, making $w_{\perp}$ deterministic given $W$.
We instead consider the following less restrictive assumptions about
$W$, $w_{\bot}$, which implicitly specify properties of $\pi(\cdot)$:
{restatable}assumptionass:subexponential-assumptionAssumption of subexponentiality
We have that

- 1.
$w_{i}w_{\bot i}$ (where $w_{i}$ is the $i^{\mathrm{th}}$ column of $W$) is sub-exponential $\forall i\in\{1,...,D\}$ and
- 2.
$\sum_{i=1}^{D}w_{i}w_{i}^{\top}$ is almost surely positive definite for any $D\geq N$.

The first condition is in particular
fulfilled whenever $w_{\bot i}$ and $w_{i}^{\top}$ are sub-Gaussian but potentially dependent,
a condition that is true for bounded activation functions or piecewise linear activation functions with standard normal weights.
The second condition is equivalent to $\Phi$ having almost surely full rank.
We note that the second condition will be violated by ReLUs and leaky-ReLUs;
however, it holds for arbitrarily precise approximations.(^2^22 E.g., $\phi_{\alpha}(\omega,x)=\frac{1}{\alpha}\log(1+e^{\alpha\omega^{\top}x})$, $\alpha>0$ yields an a.s. full-rank $\Phi$, and $\phi_{\alpha}(\omega,x)\overset{\alpha\to\infty}{\to}\mathrm{ReLU}(\omega^{ \top}x)$.)
Note that we make no assumptions about the mean, independence, or even the uncorrelation
of entries in a given column of $[W,w_{\bot}]$.

## 3 Main results

Figure: Figure 2: Empirically, the term $\mathbb{E}[w_{\bot}^{\top}W^{\top}(WW^{\top})^{-1}]$ is consistently zero. We plot the distribution of the first index of $w_{\bot}^{\top}W^{\top}(WW^{\top})^{-1}$, which captures the difference between the infinite-width single model and a smaller overparameterized RF model (see Eq. 2). (Left) We use ReLU as activation function, $x_{i}\in\mathbb{R}$, and $N=6,D=200$. (Right) We use the Gaussian Error activation function, the California Housing dataset (Kelley Pace & Barry, ), and $N=12,D=200$.
Refer to caption: x3.png

### 3.1 Equivalence of Infinite Ensembles and the Infinite-Width Single Models

We at first assume an infinite computational budget and consider the following two limiting predictors, for which we will show pointwise equivalence in predictions:

- 1.
An infinite-width least norm predictor, $h^{(\mathrm{LN})}_{\infty}$, the a.s. limit of $h^{(\mathrm{LN})}_{\mathcal{W}}$ as $|{\mathcal{W}}|=D\to\infty$
- 2.
An infinite ensemble of finite-width least norm predictors, $\bar{h}^{(LN)}_{\infty}$, which is the almost sure limit of $\bar{h}^{(LN)}_{{\mathcal{W}}_{1:M}}$ as $M\to\infty$, with $N<D<\infty$ remaining constant.

These limiting predictors do not only serve as approximations to large ensembles and very large single models but will also prove useful in characterizing the variance and generalization error of finite overparameterized ensembles, as discussed in Sec. 3.2.

Define $k_{N}(\cdot):{\mathcal{X}}\to\mathbb{R}^{N}$ as the vector of kernel evaluations with the training data $k_{N}(\cdot)=\begin{bmatrix}k(x_{1},\cdot)&\cdots&k(x_{N},\cdot)\end{bmatrix}^
{\top}\in\mathbb{R}^{N}$.
As $D\to\infty$, the
minimum norm interpolating model converges pointwise almost surely to the
ridgeless kernel regressor by the Strong Law of Large Numbers:

$$ $h^{(\mathrm{LN})}_{\mathcal{W}}(\cdot)\overset{\mathrm{a.s.}}{\longrightarrow} h^{(\mathrm{LN})}_{\infty}(\cdot),\qquad h^{(\mathrm{LN})}_{\infty}(\cdot):=k_ {N}(\cdot)^{\top}K^{-1}y.$ $$

On the other hand, using $W$ and $w_{\perp}$ as introduced in Sec. 2 we can rewrite the infinite ensemble prediction $\bar{h}^{(LN)}_{\infty}(x^{*})$ as (for a derivation of this, see Sec. B.1)

$$ $\bar{h}^{(LN)}_{\infty}(x^{*})=h^{(\mathrm{LN})}_{\infty}(x^{*})\>\>+\>\>{r_{ \bot}\mathbb{E}_{W,w_{\bot}}\left[w\bot^{\top}W^{\top}\left(WW^{\top}\right)^{ -1}\right]R^{-\top}y}$ (2) $$

Figure: Figure 3: Infinite overparameterized ensembles are equivalent to a single infinite-width model regardless of width, while underparameterized ensembles behave fundamentally differently. We present the average absolute difference between the infinite ensemble and the single infinite model for different feature counts $D$. (Left) ReLU activations, $N=6$, data are from the setting in Fig. 1. (Right) softplus activations, $N=12$, California Housing dataset. Both exhibit a “hockey stick” pattern: there is a substantial difference between the underparameterized ensemble and the infinite-width model; however, this difference vanishes in the overparameterized regime.
Refer to caption: x5.png

To prove the pointwise equivalence
of the infinite ensemble and infinite-width single model, we need to show that
$\mathbb{E}_{W,w_{\bot}}[w_{\bot}^{\top}W^{\top}(WW^{\top})^{-1}]$ term in Eq. 2 is zero.
Note that this result trivially holds when
the entries of $W$ and $w_{\perp}$ are i.i.d., as assumed in prior work (e.g. Jacot et al., ).
Here, we show that this term is zero even when $w_{\perp}$ and $W$ are dependent,
which—as described in Sec. 2—is a more realistic assumption for neural network features.
Empirically, in Fig. 2 we observe that the entries of the random variable $w_{\bot}^{\top}W^{\top}(WW^{\top})^{-1}\in\mathbb{R}^{N}$ have a mean of zero for both ReLU and Gaussian error
function features, both of which violate independence assumptions between $w_{\bot}$ and $W$ (as noted in Sec. 2).
We formalize this observation in the following lemma:
{restatable}lemmalmm:identity-infinite-model-infinite-ensemble-subexp
Under LABEL:ass:subexponential-assumption, it holds that $\mathbb{E}_{W,w_{\bot}}[w_{\bot}^{\top}W^{\top}(WW^{\top})^{-1}]=0$.

###### Proof sketch.

(See Sec. B.1 for a full proof.)
We start by applying the Woodbury formula to express the matrix inverse $(WW^{\top})^{-1}$ as a decomposition involving the matrix $A_{-i}=(WW^{\top}-w_{i}w_{i}^{\top})$ and the individual column $w_{i}$ of $W$. This decomposition yields the expression:

$$ $w_{\bot}^{\top}W^{\top}(WW^{\top})^{-1}=\operatorname{{\textstyle\sum}}_{i=1}^ {D}(w_{\bot i}w_{i}^{\top})/(1+w_{i}^{\top}A_{-i}^{-1}w_{i})A_{-i}$ $$

Next, using sub-exponential concentration inequalities in conjunction with the Weak Law of Large Numbers,
we show that the conditional expectation
$\mathbb{E}_{w_{\bot i},w_{i}}\left[(w_{\bot i}w_{i}^{\top})/(1+w_{i}^{\top}A_{
-i}^{-1}w_{i})\mid A_{-i}\right]$
exists and is zero for all invertible $A_{-i}$.
This result implies that:
$\mathbb{E}[w_{\bot}W^{\top}(WW^{\top})^{-1}]=\operatorname{{\textstyle\sum}}_{
i=1}^{D}\mathbb{E}_{A_{-i}}[\mathbb{E}_{w_{\bot i},w_{i}}[(w_{\bot i}w_{i}^{
\top})/(1+w_{i}^{\top}A_{-i}^{-1}w_{i})]A_{-i}^{-1}]=0.$
∎

Combining LABEL:lmm:identity-infinite-model-infinite-ensemble-subexp and Eq. 2 yields the pointwise equivalence of $\bar{h}^{(LN)}_{\infty}$ and $h^{(\mathrm{LN})}_{\infty}$:

theoremthm:identity-infinite-model-infinite-ensemble-subexpEquivalence of infinite-width single model and infinite ensembles
Under LABEL:ass:subexponential-assumption, the infinite ensemble of finite-width (but overparameterized) RF regressors $\bar{h}^{(LN)}_{\infty}$ is pointwise almost surely equivalent to the (single) infinite-width RF regressor $h^{(\mathrm{LN})}_{\infty}$.

This result shows that, in the case of overparameterized RF regression,
ensembling yields *exactly* the same predictions as what could be
achieved by simply increasing the capacity of a single model (see Fig. 1 for
a visualization).
As a result, we should not expect substantial differences in generalization between large single models
and overparameterized ensembles.
We emphasize a contrast with the underparameterized regime,
where RF ensembles match the generalization error of kernel ridge regression
(see Appx. D or Bach, , Sec. 10.2.2).
While width controls the implicit ridge parameter in the underparameterized regime (see Sec. 1.1),
width does not affect the ensemble predictor in the overparameterized regime.
We confirm this difference in Fig. 3
which shows that RF ensembles are equivalent to the ridgeless kernel regressor when $D>N$
but not when $D<N$.

###### Proof sketch.

### 3.2 Variance of Ensemble Predictions

We now analyze the predictive variance amongst component models in an overparameterized RF ensemble, a quantity used to quantify predictive uncertainty and provide insights about the generalization error.
Using LABEL:thm:identity-infinite-model-infinite-ensemble-subexp,
the variance of the predictions of a single RF model with respect to its random features can be expressed as (see Sec. B.2 for a derivation)

$$ $\displaystyle\mathrm{Var}_{{\mathcal{W}}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x^{* })]$ $\displaystyle=r_{\perp}^{2}\left(y^{\top}R^{-1}\>\mathbb{E}_{W,w_{\bot}}[(WW^{ \top})^{-T}Ww_{\perp}w_{\perp}^{\top}W^{\top}(WW^{\top})^{-1}]\>R^{-\top}y \right).$ (3) $$

In the special case where $W$ and $w_{\perp}$ are i.i.d. standard normal,
this expression simplifies to

$$ $\mathrm{Var}_{\mathcal{W}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x^{*})]=r_{\perp}^{ 2}\left(\tfrac{\|h^{(\mathrm{LN})}_{\infty}\|_{k}^{2}}{D-N-1}\right),$ (4) $$

where $\|h^{(\mathrm{LN})}_{\infty}\|_{K}^{2}$ represents the squared norm of $h^{(\mathrm{LN})}_{\infty}$ in the RKHS defined by the kernel $k(\cdot,\cdot)$.
From this equation, we note that the variance decreases with RF regressor width as $\sim 1/D$,
scales with the complexity of $h^{(\mathrm{LN})}_{\infty}$,
and only depends on $x^{*}$ through the quantity $r_{\perp}^{2}$ (a term we will analyze later).

Unfortunately, $\mathbb{E}_{W,w_{\bot}}[(WW^{\top})^{-T}Ww_{\perp}w_{\perp}^{\top}W^{\top}(WW^
{\top})^{-1}]$ generally does not have simple expression for arbitrary $W,w_{\bot}$ satisfying LABEL:ass:subexponential-assumption.
Without assuming *Gaussian universality*, the variance depends on $x^{*}$ through both $r_{\perp}^{2}$ as well as through the expectation from Eq. 3 involving $w_{\bot}$.
Still, prior works and empirical results (see Fig. 5 and Sec. A.3) suggest that the variance of RF models decays with $\sim 1/D$
under a variety of distributions (e.g. Adlam & Pennington, ).

Figure: Figure 4: Ensemble variance (left) and Bayesian notions of uncertainty (right) can differ significantly. For $N=6$ and $D=200$ with ReLU activations, we show the overparameterized ensemble variance (left) and the posterior variance of a Gaussian process with prior covariance $k(\cdot,\cdot)$ (right) across the input range. Empirically, we observe substantial differences between the two quantities.
Refer to caption: x7.png

#### Implications for uncertainty quantification.

A common approach to uncertainty quantification with ensembles is to examine
the predictive variance of their members at a specific test point $x^{*}$
(Lakshminarayanan et al., ).
Before diving into an analysis of Eqs. 3 and 4,
it is worth reflecting on the implications that LABEL:thm:identity-infinite-model-infinite-ensemble-subexp
has for ensemble variance as uncertainty quantification.
Because the expected overparameterized RF model is the infinite-width RF model,
we can exactly characterize the ensemble variance as the expected squared difference between the predictions
of a large (i.e., infinite-width) model versus a smaller (finite-width but still overparameterized) model.
This reveals that ensemble variance provides a non-standard notion of
uncertainty, differing from both conventional frequentist and Bayesian
interpretations.

A notable exception is when $W$ and $w_{\bot}$ are i.i.d. standard normal.
Recall by Eq. 4 that the variance under the *Gaussian universality* assumption
only depends on $x^{*}$ through the quantity $r_{\perp}^{2}$.
By Eq. 1 we see that $r_{\perp}^{2}$ is equal to $k(x^{*},x^{*})-k_{N}(x^{*})^{\top}K^{-1}k_{N}(x^{*})$,
which is exactly the Gaussian process posterior variance with prior covariance $k(\cdot,\cdot)$ (e.g. Rasmussen & Williams, ).
In this case, ensembles provide a scaled version of a classic Bayesian estimate of uncertainty.

However, relaxing from Gaussianity to LABEL:ass:subexponential-assumption makes
the relationship between ensemble variance and $r_{\perp}^{2}$ more complex, as the
independence between $W$ and $w_{\perp}$ is no longer guaranteed.
As can be seen in Eq. 3,
the variance generally depends on $x^{*}$ through both $r_{\perp}^{2}$ and a complicated expectation involving $W$ and $w_{\perp}$.
In Sec. B.2 we demonstrate with a simple example that this expectation
can indeed depend on $x^{*}$, implying that ensemble variance does not exactly correspond to a scalar multiple of $r_{\perp}^{2}$.
In our numerical experiments using realistic (i.e., non-Gaussian) random feature distributions,
(Fig. 4 and Sec. A.3),
we observe significant deviations between the ensemble variance and the Gaussian process posterior variance,
further implying that one cannot view ensembles through a classic framework of uncertainty.
These discrepancies are particularly important for uncertainty estimation in
safety-critical applications or active learning
(e.g. Gal et al., ; Beluch et al., ).

Figure: Figure 5: The variance and generalization error of overparameterized ensembles and single large models scale similarly with the total number of features. (Left) We show that the variance of a single model with $MD$ features decays as $\sim\tfrac{1}{MD}$, consistent with the scaling behavior of the variance in an ensemble. (Right) We present the generalization error of an ensemble of $M$ models, each with $D=200$ features, compared to a single model with $MD$ total features. Both exhibit nearly identical dependence on the total feature budget. The results use a ReLU activation function, the California Housing dataset, and $N=12$.
Refer to caption: x9.png

#### Ensembles versus larger single models under a finite feature budget.

Our characterization of ensemble variance also holds implications for the generalization error of ensembles versus single models under a finite computational budget.
We compare ensembles of $M$ models with $D$ features each ($\bar{h}^{(LN)}_{{\mathcal{W}}_{1:M}}=\tfrac{1}{M}\operatorname{{\textstyle\sum
}}_{m=1}^{M}h^{(\mathrm{LN})}_{{\mathcal{W}}_{m}}$) to single models with $MD$ features $h^{(\mathrm{LN})}_{{\mathcal{W}}^{*}}(\cdot)$ (i.e., here $|{\mathcal{W}}_{m}|=D$ for all $m$ and $|{\mathcal{W}}^{*}|=MD$).
The expected generalization error of either predictor
can be decomposed into standard bias and variance terms:

$$ $\mathbb{E}_{h}\left[L\left(h\right)\right]:=\mathbb{E}_{h}[\>\mathbb{E}_{x}[ \left(h(x)-\mathbb{E}[y\mid x]\right)^{2}]=L\left(\mathbb{E}_{h}\left[h\right] \right)+\mathbb{E}_{x}\left[\operatorname{Var}_{h}\left(h(x)\right)\right].$ $$

Since $h^{(\mathrm{LN})}_{{\mathcal{W}}^{*}}$ and $h^{(\mathrm{LN})}_{{\mathcal{W}}_{1}},\ldots,h^{(\mathrm{LN})}_{{\mathcal{W}}_
{M}}$ share the same expected predictor
(as established in
LABEL:thm:identity-infinite-model-infinite-ensemble-subexp), the only
difference in the generalization of $h^{(\mathrm{LN})}_{{\mathcal{W}}^{*}}$ and $\bar{h}^{(LN)}_{{\mathcal{W}}_{1:M}}$ arises from their variances.
Due to the independence between ensemble members,
we have that $\mathrm{Var}_{{\mathcal{W}}_{1:M}}[\bar{h}^{(LN)}_{{\mathcal{W}}_{1:M}}(x)]=
\tfrac{1}{M}\mathrm{Var}_{{\mathcal{W}}_{m}}[h^{(\mathrm{LN})}_{{\mathcal{W}}_
{m}}(x)]$.
Moreover, since the variance of a single RF model is inversely proportional to the number of features
(exactly in the case of Gaussian features and approximately in the general case, as discussed above),
we have that $\mathrm{Var}_{{\mathcal{W}}^{*}}[h^{(\mathrm{LN})}_{{\mathcal{W}}^{*}}(x)]\>/
\>\mathrm{Var}_{{\mathcal{W}}_{m}}[h^{(\mathrm{LN})}_{{\mathcal{W}}_{m}}]
\asymp 1/M$.
Altogether, this suggests that the generalization error of finite ensembles and finite-width single models
decay at similar rates.
We confirm this similar rate of decay in Fig. 5 and Sec. A.3,
which compare ensembles versus single models under various feature budgets.
These results further demonstrate that ensembles do not provide any meaningful generalization advantage over
single models.

### 3.3 Equivalence of the limiting predictors in the small ridge regime

Having established the pointwise equivalence between infinite ensembles and infinite-width single models in the ridgeless regime, we now extend our analysis to the case where a small ridge parameter $\lambda>0$ is introduced.
While $h^{(RR)}_{\infty,\lambda}$, the infinite-width limit of $h_{{\mathcal{W}},\lambda}^{(RR)}$ as $|{\mathcal{W}}|=D\to\infty$, almost surely converges to the kernel ridge regressor with ridge $\lambda$, the infinite ensemble $\bar{h}^{(RR)}_{\infty,\lambda}\coloneqq\mathbb{E}_{{\mathcal{W}}}[h_{{
\mathcal{W}},\lambda}^{(RR)}(x)]$ does not generally maintain pointwise equivalence with $h^{(RR)}_{\infty,\lambda}$.
This divergence occurs even under the *Gaussian universality* assumption (Jacot et al., ).
However, we hypothesize that the difference between these limiting predictors is small when $\lambda$ is close to zero, which is common in practical applications. To analyze this regime, we introduce a minor additional assumption:
{restatable}assumptionass:inverse-third-moment-constrained
We assume that $\mathbb{E}_{\mathcal{W}}[\left(\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{\top}
\right)^{-1}]$ is finite for all $|{\mathcal{W}}|=D>N$.

Figure: Figure 6: Lipschitz continuity of predictions for an infinite ensemble and kernel regressor with respect to the ridge parameter. (Left) We plot the $|\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})-\bar{h}^{(LS)}_{\infty}(x^{*})|$ as a function of $\lambda$ for 500 test points. (Right) We show the evolution of $|h^{(RR)}_{\infty,\lambda}(x^{*})-h^{(LS)}_{\infty}(x^{*})|$ for the same test points. Both plots use the ReLU activation function and the California Housing Dataset with $N=12$ and $D=200$. While the direct difference $|\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})-h^{(RR)}_{\infty,\lambda}(x^{*})|$ is not shown due to reasons outlined in Sec. A.4, it can be bounded by a sum of the shown differences (see Sec. C.1).
Refer to caption: x11.png

Under Assumptions LABEL:ass:subexponential-assumption and LABEL:ass:inverse-third-moment-constrained, we show that the difference is Lipschitz-continuous with respect to $\lambda$:

theoremthm:smoothness-difference-ridgeless-ridgeThe difference between ensembles and large single models is smooth with respect to $\lambda$.
Under Assumptions LABEL:ass:subexponential-assumption and LABEL:ass:inverse-third-moment-constrained, the difference
$|\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})-h^{(RR)}_{\infty,\lambda}(x^{*})|$
between the infinite ensemble and the single infinite-width model trained
with ridge $\lambda$ is Lipschitz-continuous in $\lambda$ for $\lambda\geq 0$. The Lipschitz constant is independent of $x^{*}$ for compact $\mathcal{X}$.

###### Proof sketch.

We first prove a lemma that shows the predictions of infinite-width RF
regressors $h^{(RR)}_{\infty,\lambda}(x^{*})$ are Lipschitz-continuous in
$\lambda$ (see LABEL:lmm:bound-difference-ridge-regressors). Using a
similar proof strategy and noting that an equivalent statement to
LABEL:lmm:identity-infinite-model-infinite-ensemble-subexp
holds in the ridge regime, we also prove that the predictions of infinite
ensembles $\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})$ are Lipschitz-continuous in
$\lambda$ (see
LABEL:lmm:bound-difference-random-feature-ridge-regressors). In
Fig. 6, we show how the differences
between these predictions and their ridgeless counterparts evolve with
respect to $\lambda$ for various test points. Combining these two results
and using a triangle inequality yields our theorem. For the full proof, see
Sec. C.1.
∎

Since LABEL:thm:identity-infinite-model-infinite-ensemble-subexp
ensures that $|\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})-h^{(RR)}_{\infty,\lambda}(x^{*})|=0$
for $\lambda=0$, we can conclude that the pointwise difference grows at most linearly with $\lambda$. Specifically, we have the
following bound

$$ $\left|\bar{h}^{(RR)}_{\infty,\lambda}-h^{(RR)}_{\infty,\lambda}(x)\right|\leq C \cdot\lambda,$ $$

for some constant $C$ independent of $x^{*}$, provided that $\mathcal{X}$ is compact.
In practical terms, this result indicates that for sufficiently small values of $\lambda$, the predictions of large ensembles and large single models remain nearly indistinguishable,
reinforcing our findings from the ridgeless regime.

###### Proof sketch.

## 4 Conclusion

For *Question 1*, we demonstrated that under weak conditions, infinite ensembles, and single infinite-width models are pointwise equivalent in the ridgeless regime and nearly identical with a small ridge, significantly expanding on prior results
(e.g. Jacot et al., ).
These results verify recent empirical findings (e.g. Abe et al., )
that much of the benefit attributed to overparameterized ensembles, such as improved predictive performance and robustness, can be explained by their similarity to larger single models.
We contrast these findings to the underparameterized regime,
where ensembling typically induces regularization and improves generalization.
Similarly, for *Question 2*, we argued that the variance reduction from ensembling is asymptotically equivalent to increasing the number of features of a single model.
This result further strengthens our findings on Question 1 and demonstrates functional similarities under relatively small computational budgets.
For *Question 3*, we found that the ensemble variance measures the
expected difference to a single larger model and is thus a non-standard measure
for uncertainty. Significant deviations from the Gaussian process posterior variance indicate that caution is needed when using ensemble variance for
uncertainty quantification, especially in safety-critical settings.
Again, these results reinforce empirical findings from (Abe et al., ) about overparameterized neural network ensembles.

Overall, while our results do not contradict the utility of overparameterized ensembles,
they suggest that their benefits may often be explained by their similarity to larger models and that further research is needed to improve uncertainty quantification methods.

#### Limitations.

The practical implications of our work are limited by the theoretical abstractions we employ.
While these abstractions provide valuable theoretical insights, they may not always hold in real-world, finite settings.
Most notably, we approximate neural networks using RF models and focus on infinite single models and infinite ensembles as approximations for large models and large ensembles.
Nevertheless, we emphasize that our theoretical results on RF regressors align with recent empirical observations on deep ensembles (Abe et al., ; Theisen et al., ),
further supporting the growing body of work that uses RF models to provide insights into deep learning phenomena (e.g. Belkin et al., ; Hastie et al., ; Simon et al., ).

In addition to these theoretical assumptions, our empirical results are
constrained in terms of scale and complexity. Due to numerical stability issues (see Sec. A.2), we
primarily considered a small number of samples, a relatively large number of
random features, and simple data-generating functions.
Again, we refer the readers to the afformentioned empirical work for larger-scale experiments.

## Reproducibility and Ethics Statements

#### Reproducibility.

The primary contribution of this paper is a theoretical analysis to explain empirical phenomena studied in the recent works of Abe et al. (; ); Theisen et al. ().
All proofs and derivations are largely self-contained, either in the main text or the appendix.
We supply references to all background material where applicable.

Empirical results are not the main focus of this work.
Nevertheless, we provide the simulation code used to generate all figures in the text,
and a complete description of the experiments can be found in Sec. A.1.
We also include a discussion on the numerical stability of our experiments in Sec. A.2.

#### Ethics.

We believe there are no significant ethical concerns stemming from this work,
as it is largely a theoretical analysis of previous empirical results.
However, we do note that this work studies ensembles of neural networks and their uncertainty estimates,
which have the potential to be used in safety-critical applications (Lakshminarayanan et al., ; Ovadia et al., ).

## Acknowledgements

Resources used in this research were provided, in part, by
the Province of Ontario, the Government of Canada, and
companies sponsoring the Vector Institute.
JPC is supported by the Gatsby Charitable Foundation (GAT3708), the Simons Foundation (542963), the NSF AI Institute for Artificial and Natural Intelligence (ARNI: NSF DBI 2229929) and the Kavli Foundation.
GP acknowledges support from NSERC and the Canada CIFAR AI Chair program.

## Appendix A Experimental setup and additional results

The code to run all our experiments can be found on GitHub: [https://github.com/nic-dern/theoretical-limitations-overparameterized-ensembles](https://github.com/nic-dern/theoretical-limitations-overparameterized-ensembles). It contains a README.md file that explains how to set up and run the experiments.

### A.1 Experimental setup

We had two setups using which we performed most of our experiments:

Figure: Figure 7: True function $f(x)=\sin(5\cdot b^{\top}x)$ with different random seeds. The blue line shows the true function, while red dots represent training samples for two distinct random seeds.
Refer to caption: x13.png

- 1.
We generate training and test points uniformly at random from $[-5,5]^{d}$ using the function $f(x)=\sin(5\cdot b^{\top}x)$, where $b$ is a vector (depending on the random seed) and the noise parameter is $\sigma=0.05$ (we assume Gaussian noise with mean 0 0). In this setting, we use $N=6,D=200$, and data from $\mathbb{R}$ (i.e., $d=1$) if not specified otherwise. You can find a plot of an example true function in Fig. 7.
- 2.
We use the California Housing (Kelley Pace & Barry, ) dataset and sample distinct training and test points from it (randomly permutating the dataset initially). In this setting, we use $N=12,D=200$ if not differently specified. The data dimension is $\mathbb{R}^{8}$ here. In contrast to the first setting, we employ a data normalization using a max-min normalization *on the entire dataset* since we experimentally found this makes our methods more stable.

We calculate the generalization error using $N=1000$ test points in both settings. In the first setting, we calculate the variance of the predictions of a single model using $M=20,000$ models, while in the second setting, we use $M=4,000$ models. Apart from Fig. 2 where we use $100,000$ samples, “infinite” ensembles consist of $M=10,000$ models.

As distribution $\tau(\cdot)$ of the elements $\omega_{i}\in{\mathcal{W}}$ we always use $\mathcal{N}(0,I)$. As activation functions, we use ReLU, the Gaussian error function, and the softplus function $\frac{1}{\beta}\cdot\log(1+\exp(\beta\cdot\omega^{\top}x))$ with $\beta=1$. For the first two activation functions, there exist analytically calculatable limiting kernels, the arc-cosine kernel (Cho & Saul, ) and the erf-kernel (Williams, ). The closed forms for these are

$$ $k_{\text{arc-cosine}}(x,x^{\prime})=\frac{1}{2\pi}\|x\|\|x^{\prime}\|\left( \sin\theta+(\pi-\theta)\cos\theta\right),$ $$

where $\theta=\cos^{-1}\left(\frac{x^{\top}x^{\prime}}{\|x\|\|x^{\prime}\|}\right)$ and

$$ $k_{\text{erf}}(x,x^{\prime})=\frac{2}{\pi}\sin^{-1}\left(\frac{2x^{\top}x^{ \prime}}{\sqrt{(1+2\|x\|^{2})(1+2\|x^{\prime}\|^{2})}}\right).$ $$

For the softplus function, we approximate the kernel by estimating the second moment $k(x,x^{\prime})=\mathbb{E}[\phi(\omega,x)\phi(\omega,x^{\prime})\mid x,x^{
\prime}]$ of the feature extraction using $10^{7}$ samples from $\tau(\cdot)$.
For sampling Gaussian features (i.e., testing under the assumption *Gaussian universality*), we use the same approach as described by Jacot et al. ().

Before training on data, we always append a $1$ in the zeroeth-dimension of the data before calculating the dot product with $\omega$ (correspondingly, the dimension of $\omega$ is $d+1$) and applying the activation function. In the ridgeless case, we use $\lambda=10^{-8}$ to avoid numerical issues.

### A.2 Notes on stability

During our experiments, we encountered challenges related to both mathematical stability (i.e., matrices being truly singular rather than nearly singular) and numerical stability. This section outlines these issues and describes the steps we took to mitigate them.

Most importantly, the matrix $\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{\top}$ is not almost surely invertible
when using the ReLU activation function, meaning that technically, the second
condition of our LABEL:ass:subexponential-assumption is not fulfilled. In
numerical experiments, this results in cases where $(\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{\top})^{-1}$ is
nearly singular (though stabilized with $\lambda=10^{-8}$).

Figure: Figure 8: Visualization of hyperplanes separating training points. We illustrate how a series of hyperplanes can separate a growing subset of the training points, leading to a triangular, invertible matrix structure as a subset of $\Phi$.
Refer to caption: x15.png

On the other hand, when $D$ is sufficiently large relative to $N$, $\Phi_{\mathcal{W}}$ is full rank with high probability, which implies that $\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{\top}$ is invertible with high probability. Given our data transformation of appending a 1 in the zeroeth dimension, one can see this as there exists a series of (non-zero probability sets of)
hyperplanes separating an increasing subset of the training points, leading to a subset of $\Phi_{\mathcal{W}}$’s columns that form a triangular, invertible matrix (see Fig. 8 for a visualization). Intuitively, higher data dimensionality and better separability of the points increase the probability of $\Phi_{\mathcal{W}}$ having full rank.

Figure: Figure 9: An adversarial example where the infinite ensemble of overparameterized RF models is numerically not equivalent to a single infinite-width RF model. (Left) We show a sample of $100$ RF models (blue) with ReLU activations trained on the same $N=15$ densely clustered data points. Additionally, we show the single infinite-width RF model (pink). (Right) We again show the single infinite-width RF model (blue) and the “infinite” ensemble of $M=10,000$ RF models (pink). A significant difference between the two models is observed in this adversarial case, indicating instability.
Refer to caption: x16.png

As an example of the discussed instabilities, see the adversarial scenario shown in Fig. 9, where $N=15$ and many training points are placed very close to each other. In this case, individual RF regressors exhibit relatively high variance output values (due to numerical instabilities), which are not averaged out in the “infinite” ensemble. Similar issues were also observed when using the Gaussian error function as the activation function, although they were generally less pronounced.

Figure: Figure 10: Using softplus activations instead of ReLU activations reduces instabilities in overparameterized RF ensembles. The plots show the average absolute difference between the predictions of an infinite ensemble and a single infinite-width model for varying feature counts $D$, using $N=12$ training samples from the California Housing dataset. (Left) ReLU activations exhibit significant instability, especially for $D>N,D\approx N$, and do not consistently show the expected pointwise equivalence between the infinite ensemble and the single infinite-width model. (Right) Softplus activations — as equivalently shown in Fig. 3 — smooth out these instabilities and more consistently show the expected pointwise equivalence.
Refer to caption: x18.png

To alleviate these issues, we used the following approaches:

- •
We used a relatively low number of samples, $N=6$ or $N=12$, compared to $D=200$. As shown in Fig. 1, even with $D=200$, there is still a considerable amount of variance in the RF regressors (i.e., the individual RF regressors are not yet closely approximating the limiting kernel ridge regressor).
- •
We appended a $1$ in the zeroeth dimension of the data before calculating the dot product with $\omega$.
- •
We performed additional experiments using the softplus function with $\beta=1$ as a smooth approximation of the ReLU activation function. This often helped stabilize the numerical computations, as seen in Fig. 10, where we repeated a part of the experiment from Fig. 3 using the ReLU function as activation function which increased the numerical instability for low $D$ values.
- •
We used a ridge term $\lambda=10^{-8}$ in the ridgeless case to stabilize the inversion of $\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{\top}$.
- •
We used *double precision* for all computations and used the torch.linalg.lstsq function with the driver gelsd (for not-well-conditioned matrices) to solve linear systems.
- •
We applied max-min normalization to the entire California Housing dataset to improve stability.

### A.3 Additional experiments for the ridgeless case

#### Additional experiments on the ensemble variance.

We observed a different behavior of the RF regressor variance and $r_{\perp}^{2}$ as shown in Fig. 4 consistently across different random seeds and dimensions for both ReLU and the Gaussian error function activations as activation functions. In Fig. 11, we present additional examples for the Gaussian error function in one dimension and the ReLU activation in two dimensions.

Figure: Figure 11: Variance and $r_{\perp}^{2}$ for different activations and dimensions. (Top left) Variance of RF model predictions across the input range for $D=200$ and $N=6$, using the erf activation function. (Top right) Corresponding $r_{\perp}^{2}$ values across the input range using the erf kernel. (Bottom left) Variance of RF model predictions across the input range for $D=200$, $p=2$, and $N=12$, using the ReLU activation function. (Bottom right) Corresponding $r_{\perp}^{2}$ values across the input range using the arc-cosine kernel.
Refer to caption: x19.png

#### Additional experiment on generalization error and variance scaling.

In Fig. 5, we demonstrated variance and generalization error decay for the ReLU activation function. To verify the consistency of these trends, we repeated the experiment using the Gaussian error function and the corresponding erf-kernel. The results are very similar, shown in Fig. 12.

Figure: Figure 12: Variance and generalization error scale similarly with the number of features, consistent with Fig. 5. In (a), the variance of a single model with $MD$ features decays as $\sim\tfrac{1}{MD}$, matching the ensemble’s behavior. In (b), the generalization error of an ensemble with $M$ models and $D=200$ features shows a similar decay to that of a single model with $MD$ features. Results use the Gaussian error function, California Housing dataset, and $N=12$.
Refer to caption: x23.png

### A.4 More experiments for the ridge case

#### Additional experiments for the convergence of the expected value term.

In Appx. C, we show that a variant of LABEL:lmm:identity-infinite-model-infinite-ensemble-subexp also holds in the ridge case. More precisely, we show that

$$ $\mathbb{E}_{W,w_{\bot}}\left[w_{\bot}^{\top}W^{\top}\left(WW^{\top}+D\cdot \lambda\cdot R^{-\top}R^{-1}\right)^{-1}\right]=0$ $$

under LABEL:ass:subexponential-assumption. We repeated the experiment from Fig. 2 for the ridge case to verify this experimentally. The results are shown in Fig. 13.

Figure: Figure 13: Empirically, the term $\mathbb{E}_{W,w_{\bot}}\left[w_{\bot}^{\top}W^{\top}\left(WW^{\top}+D\cdot \lambda\cdot R^{-\top}R^{-1}\right)^{-1}\right]$ is consistently zero. We show the empirical distribution of the third—just because it looks more interesting—index of $w_{\bot}^{\top}W^{\top}\left(WW^{\top}+D\cdot\lambda\cdot R^{-\top}R^{-1} \right)^{-1}\in\mathbb{R}^{N}$, which captures the difference in predictions between $c^{\top}\mathbb{E}_{W,w_{\bot}}\left[WW^{\top}\left(WW^{\top}+D\cdot\lambda \cdot R^{-\top}R^{-1}\right)^{-1}\right]R^{-\top}y$ and a finite-sized overparameterized RF model (see Sec. B.1). We use $\lambda=1.0$ in both plots. (Left) We use a ReLU activation function, $x_{i}\in\mathbb{R}$, and $N=6,D=200$. (Right) We use the Gaussian Error Function as activation function, the California Housing dataset, and $N=12,D=200$.
Refer to caption: x25.png

#### Additional notes.

In Fig. 6, we illustrate the Lipschitz continuity of the predictions for an infinite ensemble and a kernel regressor with respect to the ridge parameter. Rather than directly presenting the difference $\left|\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})-h^{(RR)}_{\infty,\lambda}(x^{*})\right|$, we show the evolution of $\left|\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})-\bar{h}^{(LS)}_{\infty}(x^{*})\right|$ and $\left|h^{(RR)}_{\infty,\lambda}(x^{*})-h^{(LS)}_{\infty}(x^{*})\right|$. This choice was made because the upper bound we obtained was not consistently tight for settings with large $D$. In particular, the pointwise predictions of the infinite ensemble $\bar{h}^{(RR)}_{\infty,\lambda}$ and the single infinite-width model $h^{(RR)}_{\infty,\lambda}$ trained with ridge $\lambda$ were already very close for non-zero $\lambda$. We opted to display the upper bounds rather than the direct difference to avoid cherry-picking favorable settings.

Our best explanation for this phenomenon is that infinite ensembles under LABEL:ass:subexponential-assumption in the ridge regime often behave similarly to the single infinite-width model $h^{(RR)}_{\infty,\tilde{\lambda}}$ with an *implicit ridge* parameter $\tilde{\lambda}$, which solves the equation

$$ $\tilde{\lambda}=\lambda+\frac{\tilde{\lambda}}{D}\sum_{i=1}^{N}\frac{d_{i}}{ \tilde{\lambda}+d_{i}}$ $$

where $d_{i}$ are the eigenvalues of the kernel matrix $K$, as shown by Jacot et al. () under *Gaussian universality*. Intuitively and empirically, for large $D$, the implicit ridge $\tilde{\lambda}$ tends to be very close to the true ridge $\lambda$. Using LABEL:lmm:bound-difference-ridge-regressors, this suggests that for small values of $\lambda$, the difference between the infinite ensemble and the infinite-width single model $h^{(RR)}_{\infty,\lambda}$ with ridge $\lambda$ is already minimal before $\lambda$ approaches zero.

Interestingly, our findings (see Fig. 3) suggest that in the ridgeless case, the similarity to the ridge regressor with the implicit ridge only holds in the overparameterized regime. Note that this does not violate the results from Jacot et al. () since the constants in their bounds blow up as $\lambda\rightarrow 0$ in both the underparameterized and overparameterized regimes.

## Appendix B Proofs for Overparameterized Ridgeless Regression

### B.1 Equivalence of Infinite Ensemble and Infinite Single Model.

We start by proving the equivalent formulation of the infinite ensemble prediction stated in Eq. 2 using the terms $W$ and $w_{\perp}$ as introduced in Sec. 2:

###### Proof.

Defining $\phi_{\mathcal{W}}^{*}=[\phi(\omega_{i},x^{*})]_{i}\in\mathbb{R}^{D}$,
we have

$$ $\displaystyle\bar{h}_{\infty}(x^{*})$ $\displaystyle=\mathbb{E}_{{\mathcal{W}}}\left[\tfrac{1}{D}\phi_{\mathcal{W}}^{ *}\Phi_{\mathcal{W}}^{\top}\left(\tfrac{1}{D}\cdot\Phi_{\mathcal{W}}\Phi_{ \mathcal{W}}^{\top}\right)^{-1}\right]y$ $\displaystyle=\mathbb{E}_{W,w_{\bot}}\left[\left(c^{\top}W+r_{\bot}w_{\bot}^{ \top}\right)W^{\top}R\left(R^{\top}WW^{\top}R\right)^{-1}\right]y$ $\displaystyle=\mathbb{E}_{W,w_{\bot}}\left[\left(c^{\top}W+r_{\bot}w_{\bot}^{ \top}\right)W^{\top}\left(WW^{\top}\right)^{-1}\right]R^{-\top}y$ $\displaystyle=c^{\top}R^{-\top}y\>\>+\>\>r_{\bot}\mathbb{E}_{W,w_{\bot}}\left[ w_{\bot}^{\top}W^{\top}\left(WW^{\top}\right)^{-1}\right]R^{-\top}y,$ (5) $$

where $c,R,r_{\bot}$ are as defined in Eq. 1.
The left term in Eq. 5
is equal to $h^{(\mathrm{LN})}_{\infty}(x)$:

$$ $h^{(\mathrm{LN})}_{\infty}(x^{*})=\begin{bmatrix}k(x_{i},x^{*})\end{bmatrix}_{ i=1}^{N}K^{-1}y=c^{\top}RR^{-1}R^{-\top}y=c^{\top}R^{-\top}y.$ $$

∎

In the case of $\lambda>0$, we can similarly see that

$$ $\displaystyle\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})$ $\displaystyle=\mathbb{E}_{{\mathcal{W}}}\left[\tfrac{1}{D}\phi_{\mathcal{W}}^{ *}\Phi_{\mathcal{W}}^{\top}\left(\tfrac{1}{D}\cdot\Phi_{\mathcal{W}}\Phi_{ \mathcal{W}}^{\top}+\lambda I\right)^{-1}\right]y$ $\displaystyle=\mathbb{E}_{W,w_{\bot}}\left[\left(c^{\top}W+r_{\bot}w_{\bot}^{ \top}\right)W^{\top}R\left(R^{\top}WW^{\top}R+D\cdot\lambda\cdot R^{\top}R^{- \top}R^{-1}R\right)^{-1}\right]y$ $\displaystyle=\mathbb{E}_{W,w_{\bot}}\left[\left(c^{\top}W+r_{\bot}w_{\bot}^{ \top}\right)W^{\top}\left(WW^{\top}+D\cdot\lambda\cdot R^{-\top}R^{-1}\right)^ {-1}\right]R^{-\top}y$ $\displaystyle=c^{\top}\mathbb{E}_{W,w_{\bot}}\left[WW^{\top}\left(WW^{\top}+D \cdot\lambda\cdot R^{-\top}R^{-1}\right)^{-1}\right]R^{-\top}y$ $\displaystyle+r_{\bot}\mathbb{E}_{W,w_{\bot}}\left[w_{\bot}^{\top}W^{\top} \left(WW^{\top}+D\cdot\lambda\cdot R^{-\top}R^{-1}\right)^{-1}\right]R^{-\top}y.$ (6) $$

Next up, we show that the expected value $\mathbb{E}_{W,w_{\bot}}\left[w_{\bot}^{\top}W^{\top}\left(WW^{\top}\right)^{-1
}\right]$ is zero under LABEL:ass:subexponential-assumption. This directly implies the pointwise equivalence of the infinite ensemble and the single infinite-width model (see LABEL:thm:identity-infinite-model-infinite-ensemble-subexp).

###### Proof.

Define $A_{-i}=(WW^{\top}-w_{i}w_{i}^{\top})$. Note that $A_{-1}$ is almost surely invertible and positive definite by assumption LABEL:ass:subexponential-assumption.

By the Woodbury formula, for almost every $WW^{\top}$ we have that

$$ $(WW^{\top})^{-1}=(A_{-i}+w_{i}w_{i}^{\top})^{-1}=A_{-i}^{-1}-\tfrac{A_{-i}^{-1 }w_{i}w_{i}^{\top}A_{-i}^{-1}}{1+w_{i}^{\top}A_{-i}^{-1}w_{i}},$ $$

which implies that

$$ $\displaystyle w_{\bot}^{\top}W^{\top}(WW^{\top})^{-1}$ $\displaystyle=\operatorname{{\textstyle\sum}}_{i=1}^{D}w_{\bot i}w_{i}^{\top} \left(A_{-i}^{-1}-\tfrac{A_{-i}^{-1}w_{i}w_{i}^{\top}A_{-i}^{-1}}{1+w_{i}^{ \top}A_{-i}^{-1}w_{i}}\right)$ $\displaystyle=\operatorname{{\textstyle\sum}}_{i=1}^{D}w_{\bot i}\left(w_{i}^{ \top}\tfrac{A_{-i}^{-1}+w_{i}^{\top}A_{-i}^{-1}w_{i}^{\top}A_{-i}^{-1}w_{i}}{1 +w_{i}^{\top}A_{-i}^{-1}w_{i}}-\tfrac{w_{i}^{\top}A_{-i}^{-1}w_{i}w_{i}^{\top} A_{-i}^{-1}}{1+w_{i}^{\top}A_{-i}^{-1}w_{i}}\right)$ $\displaystyle=\operatorname{{\textstyle\sum}}_{i=1}^{D}\tfrac{w_{\bot i}w_{i}^ {\top}}{1+w_{i}^{\top}A_{-i}^{-1}w_{i}}A_{-i}^{-1}.$ $$

For any positive definite matrix $B\in\mathbb{R}^{N\times N}$ and any vector $v\in\mathbb{R}^{N};\|v\|=1$ and any $i\in\{1,...,D\}$, we have

$$ $\displaystyle\left|\mathbb{E}_{w_{\bot i},w_{i}}\left[\tfrac{w_{\bot i}w_{i}^{ \top}}{1+w_{i}^{\top}Bw_{i}}\right]v\right|$ $\displaystyle\leq\mathbb{E}_{w_{\bot i},w_{i}}\left[\left|\tfrac{w_{\bot i}w_{ i}^{\top}v}{1+w_{i}^{\top}Bw_{i}}\right|\right]$ $\displaystyle=\int_{0}^{\infty}\mathbb{P}\left[\left|\tfrac{w_{\bot i}w_{i}^{ \top}v}{1+w_{i}^{\top}Bw_{i}}\right|\geq t\right]dt$ $\displaystyle=\int_{0}^{\infty}\mathbb{P}\left[\left|w_{\bot i}w_{i}^{\top}v \right|\geq\left(1+w_{i}^{\top}Bw_{i}\right)t\right]dt$ $\displaystyle\leq\int_{0}^{\infty}\mathbb{P}\left[\left|w_{\bot i}w_{i}^{\top} \right|>t\right]dt$ $\displaystyle\leq\int_{0}^{\nu^{2}/\alpha}2\exp\left(-\tfrac{t^{2}}{2\nu} \right)dt+\int_{\nu^{2}/\alpha}^{\infty}2\exp\left(-\tfrac{t}{2\alpha}\right)dt,$ (7) $$

where the last inequality is a standard sub-exponential bound applied to
$w_{\bot i}w_{i}$. Note that we here use the fact that $\mathbb{E}[w_{\bot i}w_{i}^{\top}]=0$ and the $(\nu^{2},\alpha)$-sub-exponentiality of $\left|w_{\bot i}w_{i}^{\top}\right|$.

Since the last two integrals in Sec. B.1 are finite, the expectation
$\mathbb{E}_{W,w_{\bot}}\left[(w_{\bot i}w_{i}^{\top})/(1+w_{i}^{\top}Bw_{i})
\right]v$
is finite.
By the weak law of large numbers, for i.i.d. random variables $w_{i}^{(j)}$ and $w_{\bot i}^{(j)}$ across different $j$’s, we have

$$ $\mathbb{P}\left[\left|\tfrac{1}{M}\operatorname{{\textstyle\sum}}_{j=1}^{M} \tfrac{w_{\bot i}^{(j)}(w_{i}^{(j)})^{\top}v}{1+(w_{i}^{(j)})^{\top}Bw_{i}^{(j )}}-\mathbb{E}_{W,w_{\perp}}\left[\tfrac{w_{\bot i}w_{i}^{\top}}{1+w_{i}^{\top }Bw_{i}}\right]v\right|>t\right]\to 0,$ $$

for any $t>0$ and $v\in\mathbb{R}^{N}$ such that $\|v\|=1$ as $M\to\infty$.
At the same time, repeating the sub-exponential argument above,
we have that

$$ $\displaystyle\mathbb{P}\left[\left|\tfrac{1}{M}\operatorname{{\textstyle\sum}} _{j=1}^{M}\tfrac{w_{\bot i}^{(j)}(w_{i}^{(j)})^{\top}v}{1+(w_{i}^{(j)})^{\top} Bw_{i}^{(j)}}\right|>t\right]$ $\displaystyle\leq\mathbb{P}\left[\left|\tfrac{1}{M}\operatorname{{\textstyle \sum}}_{i=1}^{M}w_{\bot i}^{(j)}(w_{i}^{(j)})^{\top}v\right|>t\right]$ $\displaystyle\leq\begin{cases}2\exp\left(-\tfrac{Mt^{2}}{2\nu}\right)&0<t\leq \nu^{2}/\alpha\\ 2\exp\left(-\tfrac{Mt}{2\alpha}\right)&t>\nu^{2}/\alpha\end{cases}$ $\displaystyle\to 0$ $$

as $M\to\infty$. Here we use the property that the sum of $M$ $(\nu^{2},\alpha)$-sub-exponential random variables is $(M\nu^{2},\alpha)$-sub-exponential.

Together, these results imply that
$\mathbb{E}_{W,w_{\bot}}\left[(w_{\bot i}w_{i}^{\top})/(1+w_{i}^{\top}Bw_{i})
\right]=0$
for every positive definite $B$.
Since the random matrix $A_{-i}$ is positive semidefinite, almost surely invertible (by the second half of LABEL:ass:subexponential-assumption), and independent of $w_{i},w_{\bot i}$, we have that

$$ $\displaystyle\mathbb{E}_{w_{\bot},W}\left[w_{\bot}W^{\top}\left(WW^{\top} \right)^{-1}\right]$ $\displaystyle=\operatorname{{\textstyle\sum}}_{i=1}^{D}\mathbb{E}_{w_{\bot i}, w_{i},A_{-i}}\left[\tfrac{w_{\bot i}w_{i}^{\top}}{1+w_{i}^{\top}A_{-i}^{-1}w_{ i}}A_{-i}^{-1}\right]$ $\displaystyle=\operatorname{{\textstyle\sum}}_{i=1}^{D}\mathbb{E}_{A_{-i}} \left[\mathbb{E}_{w_{\bot i},w_{i}}\left[\tfrac{w_{\bot i}w_{i}^{\top}}{1+w_{i }^{\top}A_{-i}^{-1}w_{i}}\right]A_{-i}^{-1}\right]=0.$ $$

∎

We remark that this proof equivalently holds for the ridge-regression case, i.e., $\mathbb{E}_{W,w_{\bot}}\left[w_{\bot}^{\top}W^{\top}\left(WW^{\top}+D\cdot
\lambda\cdot R^{-\top}R^{-1}\right)^{-1}\right]=0$ since the proof does not rely on the specific form of the matrix $A_{-i}$ other than it being positive definite. Thus by Sec. B.1 we directly get that under LABEL:ass:subexponential-assumption it holds that

$$ $\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})=c^{\top}\mathbb{E}_{W,w_{\bot}}\left[WW ^{\top}\left(WW^{\top}+D\cdot\lambda\cdot R^{-\top}R^{-1}\right)^{-1}\right]R^ {-\top}y.$ (8) $$

###### Proof.

###### Proof.

### B.2 Variance of Ensemble Predictions.

In the next step, we show the formula for the variance of a single model prediction under *Gaussian universality*. Note that one could also get this result by slightly extending proofs by Jacot et al. ().

lemmathm:variance-gaussian-universalityVariance of single model predictions
Under *Gaussian universality* and assuming $D>N+1$, the variance of single model prediction at a test point $x^{*}$ is given by

$$ $\mathrm{Var}_{{\mathcal{W}}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x^{*})]=r_{\perp} ^{2}\frac{\|h^{(\mathrm{LN})}_{\infty}\|_{\mathcal{H}}^{2}}{D-N-1},$ (9) $$

where $\|\cdot\|_{\mathcal{H}}$ is norm defined by the RKHS associated with kernel $k(\cdot,\cdot)$.

###### Proof.

We start by writing down the variance of the prediction of a single model:

$$ $\mathrm{Var}_{{\mathcal{W}}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x^{*})]=\mathbb{E }_{{\mathcal{W}}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x^{*})^{2}]-\mathbb{E}_{{ \mathcal{W}}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x^{*})]^{2}$ $$

Using LABEL:thm:identity-infinite-model-infinite-ensemble-subexp, the definition of the prediction of a single model and the definition of $W$ and $w_{\perp}$, we can expand this expression:

$$ $\displaystyle=\mathbb{E}_{{\mathcal{W}}}[\phi_{\mathcal{W}}^{*}\Phi_{\mathcal{ W}}^{\top}(\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{\top})^{-1}yy^{\top}(\Phi_{ \mathcal{W}}\Phi_{\mathcal{W}}^{\top})^{-\top}\Phi_{\mathcal{W}}\phi_{\mathcal {W}}^{*\top}]-(h^{(\mathrm{LN})}_{\infty}(x^{*}))^{2}$ $\displaystyle=\mathbb{E}_{W,w_{\bot}}[(r_{\perp}w_{\perp}^{\top}+c^{\top}W)W^{ \top}R\left(R^{\top}WW^{\top}R\right)^{-1}yy^{\top}\left(R^{\top}WW^{\top}R \right)^{-\top}R^{\top}W(r_{\perp}w_{\perp}^{\top}+c^{\top}W)^{\top}]$ $\displaystyle-(h^{(\mathrm{LN})}_{\infty}(x))^{2}$ $\displaystyle=\mathbb{E}_{W,w_{\bot}}[(r_{\perp}w_{\perp}^{\top}+c^{\top}W)W^{ \top}\left(WW^{\top}\right)^{-1}R^{-\top}yy^{\top}R^{-1}\left(WW^{\top}\right) ^{-\top}W(r_{\perp}w_{\perp}^{\top}+c^{\top}W)^{\top}]$ $\displaystyle-(h^{(\mathrm{LN})}_{\infty}(x))^{2}$ $\displaystyle=(c^{\top}R^{-\top}y)^{2}-(h^{(\mathrm{LN})}_{\infty}(x))^{2}$ $\displaystyle+2\cdot r_{\bot}^{\top}\mathbb{E}_{W,w_{\bot}}[w_{\perp}^{\top}W^ {\top}(WW^{\top})^{-1}]R^{-\top}yy^{\top}R^{-1}c$ $\displaystyle+r_{\perp}^{2}\mathbb{E}_{W,w_{\bot}}[w_{\perp}^{\top}W^{\top}(WW ^{\top})^{-1}R^{-\top}yy^{\top}R^{-1}(WW^{\top})^{-T}Ww_{\perp}]$ $$

Now we can see that the first two terms cancel out (since $h^{(\mathrm{LN})}_{\infty}(x)=c^{\top}R^{-\top}y$) and the third term is zero by LABEL:lmm:identity-infinite-model-infinite-ensemble-subexp. We are left with the fourth term, which we can slightly rewrite:

$$ $\displaystyle\mathrm{Var}_{{\mathcal{W}}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x^{* })]$ $\displaystyle=r_{\perp}^{2}\mathbb{E}_{W,w_{\bot}}[w_{\perp}^{\top}W^{\top}(WW ^{\top})^{-1}R^{-\top}yy^{\top}R^{-1}(WW^{\top})^{-T}Ww_{\perp}]$ $\displaystyle=r_{\perp}^{2}y^{\top}R^{-1}\mathbb{E}_{W,w_{\bot}}[(WW^{\top})^{ -T}Ww_{\perp}w_{\perp}^{\top}W^{\top}(WW^{\top})^{-1}]R^{-\top}y$ (10) $$

Using the tower rule for conditional expectations, we have:

$$ $\displaystyle\mathrm{Var}_{{\mathcal{W}}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x)]$ $\displaystyle=r_{\perp}^{2}y^{\top}R^{-1}\mathbb{E}_{W,w_{\perp}}[(WW^{\top})^ {-T}Ww_{\perp}w_{\perp}^{\top}W^{\top}(WW^{\top})^{-1}]R^{-\top}y$ $\displaystyle=r_{\perp}^{2}y^{\top}R^{-1}\mathbb{E}_{W}[(WW^{\top})^{-T}W \mathbb{E}_{w_{\perp}|W}[w_{\perp}w_{\perp}^{\top}|W]W^{\top}(WW^{\top})^{-1}] R^{-\top}y$ $$

Since the *Gaussian universality* assumption implies $W$ and $w_{\perp}$ are independent, we get:

$$ $\displaystyle\mathrm{Var}_{{\mathcal{W}}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x)]$ $\displaystyle=r_{\perp}^{2}y^{\top}R^{-1}\mathbb{E}_{W}[(WW^{\top})^{-T}W \mathbb{E}_{w_{\perp}}[w_{\perp}w_{\perp}^{\top}]W^{\top}(WW^{\top})^{-1}]R^{- \top}y$ $$

Moreover, since by *Gaussian universality* $w_{\perp}$ and $W$ are multivariate Gaussians with the identity matrix as covariance, we get (via the expected value of a Wishart and an inverse Wishart distribution; note that for getting this expected value, we need to assume that $D>N+1$):

$$ $\displaystyle\mathrm{Var}_{{\mathcal{W}}}[h^{(\mathrm{LN})}_{\mathcal{W}}(x)]$ $\displaystyle=r_{\perp}^{2}y^{\top}R^{-1}\mathbb{E}_{W}[(WW^{\top})^{-T}(WW^{ \top})(WW^{\top})^{-1}]R^{-\top}y$ $\displaystyle=r_{\perp}^{2}y^{\top}R^{-1}\mathbb{E}_{W}[(WW^{\top})^{-T}]R^{- \top}y$ $\displaystyle=r_{\perp}^{2}\tfrac{y^{\top}R^{-1}R^{-\top}y}{D-N-1}$ $\displaystyle=r_{\perp}^{2}\tfrac{y^{\top}K^{-1}y}{D-N-1}.$ $$

Recognizing that $y^{\top}K^{-1}y=\|h^{(\mathrm{LN})}_{\infty}\|_{\mathcal{H}}^{2}$ (e.g. Wainwright, , Ch. 12) completes the proof.
∎

An equivalent argument does not work under the more general LABEL:ass:subexponential-assumption since $w_{\perp}$ and $W$ are not necessarily independent.
Even in the case of independence, $\mathbb{E}_{W}[(WW^{\top})^{-1}]$ might not be known.

###### Proof.

#### Counterexample for subexponential case.

We now give an explicit counterexample showing that when only assuming uncorrelatedness between $W$ and $w_{\perp}$ the term

$$ $E\coloneqq\mathbb{E}_{W,w_{\perp}}[(WW^{\top})^{-T}Ww_{\perp}w_{\perp}^{\top}W ^{\top}(WW^{\top})^{-1}]$ $$

from Sec. B.2 depends on $x^{*}$ implying that the variance does not only depend on $x^{*}$ via $r_{\perp}^{2}$.

Let us assume $N=D=1$ and let $W$ be uniformly distributed across the set $\left\{-\tfrac{4}{\sqrt{12.5}},-\tfrac{3}{\sqrt{12.5}},\tfrac{3}{\sqrt{12.5}},
\tfrac{4}{\sqrt{12.5}}\right\}$. Then we have $\mathbb{E}[W]=0$ and $\mathbb{E}[W^{2}]=\tfrac{1}{2}\cdot\tfrac{16}{12.5}+\tfrac{1}{2}\cdot\tfrac{9}
{12.5}=1$.

Now consider an $x^{*}$ that produces a $w_{\perp}$ so that $w_{\perp}=\sqrt{2}$ when $W=\left\{-\tfrac{3}{\sqrt{12.5}},\tfrac{3}{\sqrt{12.5}}\right\}$ and $w_{\perp}=0$ otherwise. Then we have $\mathbb{E}[w_{\perp}^{\top}W]=0$ and $\mathbb{E}[w_{\perp}^{2}]=1$. The value of $E$ is now $\tfrac{12.5}{9}$.

Furthermore, consider an $x^{*}$ that produces a $w_{\perp}$ so that $w_{\perp}=\sqrt{2}$ when $W=\left\{-\tfrac{4}{\sqrt{12.5}},\tfrac{4}{\sqrt{12.5}}\right\}$ and $w_{\perp}=0$ otherwise. Then we have $\mathbb{E}[w_{\perp}^{\top}W]=0$ and $\mathbb{E}[w_{\perp}^{2}]=1$. The value of $E$ is now $\tfrac{12.5}{16}$.

## Appendix C Proofs for Overparameterized Ridge Regression

### C.1 Difference between the infinite ensemble and infinite single model.

We begin with a lemma, which shows that the prediction of kernel regressors is Lipschitz-continuous in $\lambda$ for any $x^{*}$ and $\lambda\geq 0$. We will denote the kernel ridge regressor with regularization parameter $\lambda$ as $h^{(RR)}_{\infty,\lambda}$, as introduced in Sec. 3.3.

lemmalmm:bound-difference-ridge-regressorsBound on the difference between the kernel ridge regressors

Let $\lambda,\lambda^{\prime}\geq 0$ be two regularization parameters. Then, for any $x^{*}\in{\mathcal{X}}$ it holds that:

$$ $|h^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})-h^{(RR)}_{\infty,\lambda}(x^{*})| \leq\sqrt{n}\cdot C_{1}\cdot|\lambda^{\prime}-\lambda|\cdot\sqrt{y^{T}K^{-4}y}$ $$

where we assume $k(x_{i},x^{*})\leq C_{1}$ for all $i\in[N]$.

###### Proof.

We can write the kernel ridge regressors as $h^{(RR)}_{\infty,\lambda}(x^{*})=\operatorname{{\textstyle\sum}}_{i=1}^{n}
\alpha_{1,i}k(x_{i},x^{*})$ and $h^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})=\operatorname{{\textstyle\sum}}_{i=1
}^{n}\alpha_{2,i}k(x_{i},x^{*})$ with coefficients $\alpha_{1}$ and $\alpha_{2}$ given by:

$$ $\displaystyle\alpha_{1}$ $\displaystyle=(K+\lambda I)^{-1}y$ $\displaystyle\alpha_{2}$ $\displaystyle=(K+\lambda^{\prime}I)^{-1}y$ $$

We now write $y$ in the orthonormal basis of the eigenvectors of $K$, i.e. $y=\operatorname{{\textstyle\sum}}_{i=1}^{n}a_{i}v_{i}$. We call the corresponding eigenvalues of $K$ $d_{1},\ldots,d_{n}>0$.

The matrix $(K+\lambda I)^{-1}$ has the same eigenvectors as $K$ and the eigenvalues are $0<\tilde{d}_{i}=\tfrac{1}{d_{i}+\lambda}\leq\tfrac{1}{\lambda}$. Thus, we can write $\alpha_{1}=\operatorname{{\textstyle\sum}}_{i=1}^{n}a_{i}\tfrac{1}{d_{i}+
\lambda}v_{i}$ and $\alpha_{2}=\operatorname{{\textstyle\sum}}_{i=1}^{n}a_{i}\tfrac{1}{d_{i}+
\lambda^{\prime}}v_{i}$.

In the next step, we bound $\|\alpha_{1}-\alpha_{2}\|_{2}^{2}$: Using the orthonormality of the eigenvectors, we get:

$$ $\|\alpha_{1}-\alpha_{2}\|_{2}^{2}=\operatorname{{\textstyle\sum}}_{i=1}^{n} \left(a_{i}\left(\tfrac{1}{d_{i}+\lambda}-\tfrac{1}{d_{i}+\lambda^{\prime}} \right)\right)^{2}\\$ $$

Now we bound $\left|\tfrac{1}{\lambda+d_{i}}-\tfrac{1}{\lambda^{\prime}+d_{i}}\right|\leq
\left|\tfrac{\lambda^{\prime}-\lambda}{\lambda\lambda^{\prime}+(\lambda+
\lambda^{\prime})d_{i}+d_{i}^{2}}\right|\leq\tfrac{|\lambda^{\prime}-\lambda|}
{d_{i}^{2}}$ which gives us:

$$ $\|\alpha_{1}-\alpha_{2}\|_{2}^{2}\leq\operatorname{{\textstyle\sum}}_{i=1}^{n} \left(\tfrac{a_{i}|\lambda^{\prime}-\lambda|}{d_{i}^{2}}\right)^{2}\leq| \lambda^{\prime}-\lambda|^{2}y^{T}K^{-4}y$ $$

Using this result, we can bound the difference between the predictions of the two kernel regressors at a single point $x^{*}$:

$$ $|h^{(RR)}_{\infty,\lambda}(x^{*})-h^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})|= \left|\operatorname{{\textstyle\sum}}_{i=1}^{n}(\alpha_{1,i}-\alpha_{2,i})k(x_ {i},x^{*})\right|\leq\operatorname{{\textstyle\sum}}_{i=1}^{n}|\alpha_{1,i}- \alpha_{2,i}|k(x_{i},x^{*})$ $$

Since $k(x_{i},x^{*})\leq C_{1}$, we get (using the relation between the 1-norm and the 2-norm):

$$ $|f_{\lambda}(x^{*})-f_{\lambda^{\prime}}(x^{*})|\leq C_{1}\operatorname{{ \textstyle\sum}}_{i=1}^{n}|\alpha_{1,i}-\alpha_{2,i}|\leq C_{1}\|\alpha_{1}- \alpha_{2}\|_{2}\sqrt{n}\leq\sqrt{n}\cdot C_{1}\cdot|\lambda^{\prime}-\lambda| \cdot\sqrt{y^{\top}K^{-4}y}$ $$

∎

Using similar arguments, we now show that the expected prediction of RF regressors, i.e., the prediction of the infinite ensemble of RF regressors, is Lipschitz-continuous for any $x^{*}$ and $\lambda\geq 0$:

lemmalmm:bound-difference-random-feature-ridge-regressorsBound on the difference between expected RF Regressors
Under LABEL:ass:subexponential-assumption and LABEL:ass:inverse-third-moment-constrained, the expected value of the prediction of RF regressors is Lipschitz-continuous in $\lambda$ for any $x^{*}$ and $\lambda\geq 0$, i.e., for any $x^{*}$ it holds that:

$$ $|\bar{h}^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})-\bar{h}^{(RR)}_{\infty, \lambda}(x^{*})|\leq\|c^{\top}R^{-\top}\|\|y\|DC_{2}\left|\lambda^{\prime}- \lambda\right|$ $$

where $C_{2}$ is a constant depending on the distribution of $\Phi$.

###### Proof.

We use the characterization of $\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})$ from Eq. 8, which gives us the difference as

$$ $\left|c^{\top}\mathbb{E}_{W,w_{\bot}}\left[WW^{\top}\left(\left(WW^{\top}+D \cdot\lambda^{\prime}\cdot R^{-\top}R^{-1}\right)^{-1}-\left(WW^{\top}+D\cdot \lambda^{\prime}\cdot R^{-\top}R^{-1}\right)^{-1}\right)\right]R^{-\top}y \right|.$ $$

We can now reverse some steps we made to get this characterization and write it in terms of $\Phi$ again:

$$ $\left|c^{\top}R^{-\top}\mathbb{E}_{{\mathcal{W}}}\left[\Phi_{\mathcal{W}}\Phi_ {\mathcal{W}}^{\top}\left(\left(\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{\top}+D \cdot\lambda^{\prime}\cdot I\right)^{-1}-\left(\Phi_{\mathcal{W}}\Phi_{ \mathcal{W}}^{\top}+D\cdot\lambda\cdot I\right)^{-1}\right)\right]y\right|.$ $$

And now, using Jensen’s inequality and the convexity of the two-norm, we can pull out the expected value to the outside of the difference:

$$ $\|c^{\top}R^{-\top}\|\cdot\mathbb{E}_{{\mathcal{W}}}\left[\|\Phi_{\mathcal{W}} \Phi_{\mathcal{W}}^{\top}\left(\left(\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{ \top}+D\cdot\lambda^{\prime}\cdot I\right)^{-1}-\left(\Phi_{\mathcal{W}}\Phi_{ \mathcal{W}}^{\top}+D\cdot\lambda\cdot I\right)^{-1}\right)y\|\right].$ $$

Similarly to the proof of LABEL:lmm:bound-difference-ridge-regressors, we can write $y$ in the orthonormal basis of the eigenvectors of $\Phi\Phi^{\top}$
(note that we drop the subscript ${\mathcal{W}}$ for notational simplicity),
i.e. $y=\operatorname{{\textstyle\sum}}_{i=1}^{n}a_{i}v_{i}$. Furthermore we define the eigenvalues of $\Phi\Phi^{\top}$ as $d_{1},\ldots,d_{n}>0$. The matrix $(\Phi\Phi^{\top}+D\cdot\lambda I)^{-1}$ again has the same eigenvectors as $\Phi\Phi^{\top}$ and the eigenvalues are $0<\tfrac{1}{d_{i}+D\cdot\lambda}\leq\tfrac{1}{D\cdot\lambda}$.

Multiplying $y$ with $\Phi\Phi^{\top}(\Phi\Phi^{\top}+D\cdot\lambda I)^{-1}$ and $\Phi\Phi^{\top}(\Phi\Phi^{\top}+D\cdot\lambda^{\prime}I)^{-1}$ then gives us:

$$ $\displaystyle\Phi\Phi^{\top}(\Phi\Phi^{\top}+D\cdot\lambda I)^{-1}y$ $\displaystyle=\operatorname{{\textstyle\sum}}_{i=1}^{n}a_{i}\tfrac{d_{i}}{d_{i }+D\cdot\lambda}v_{i}$ $\displaystyle\Phi\Phi^{\top}(\Phi\Phi^{\top}+D\cdot\lambda^{\prime}I)^{-1}y$ $\displaystyle=\operatorname{{\textstyle\sum}}_{i=1}^{n}a_{i}\tfrac{d_{i}}{d_{i }+D\cdot\lambda^{\prime}}v_{i}$ $$

We can now calculate the difference of these two vectors using the orthonormality of the eigenvectors:

$$ $\|\Phi\Phi^{\top}(\Phi\Phi^{\top}+D\cdot\lambda^{\prime}I)^{-1}y-\Phi\Phi^{ \top}(\Phi\Phi^{\top}+D\cdot\lambda I)^{-1}y\|_{2}^{2}=\operatorname{{ \textstyle\sum}}_{i=1}^{n}\left(a_{i}\left(\tfrac{d_{i}}{d_{i}+D\cdot\lambda}- \tfrac{d_{i}}{d_{i}+D\cdot\lambda^{\prime}}\right)\right)^{2}$ $$

Now we look at the difference between the two coefficients and see that for each $i$, we have:

$$ $\left|\tfrac{d_{i}}{d_{i}+D\cdot\lambda}-\tfrac{d_{i}}{d_{i}+D\cdot\lambda^{ \prime}}\right|\leq\tfrac{D\cdot|\lambda^{\prime}-\lambda|}{d_{i}}$ $$

Thus, we have that the difference is bounded by:

$$ $\|\Phi\Phi^{\top}(\Phi\Phi^{\top}+D\cdot\lambda^{\prime}I)^{-1}y-\Phi\Phi^{ \top}(\Phi\Phi^{\top}+D\cdot\lambda I)^{-1}y\|_{2}^{2}\leq\tfrac{D^{2}\cdot| \lambda-\lambda^{\prime}|^{2}}{d_{N}^{2}}\|y\|_{2}^{2}.$ $$

All together, we can now bound the difference of the expected values of the predictions of RF regressors via:

$$ $|\bar{h}^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})-\bar{h}^{(RR)}_{\infty, \lambda}(x^{*})|\leq\|c^{\top}R^{-\top}\|\|y\|D|\lambda^{\prime}-\lambda| \mathbb{E}_{d_{N}}\left[\tfrac{1}{d_{N}}\right]$ $$

Since $\operatorname{tr}((\Phi\Phi^{\top})^{-1})=\operatorname{{\textstyle\sum}}_{i=1
}^{n}\tfrac{1}{d_{i}}$, and the trace is a linear operator, we can write:

$$ $\mathbb{E}_{d_{N}}\left[\tfrac{1}{d_{N}}\right]\leq\mathbb{E}_{{\mathcal{W}}} \left[(\operatorname{tr}(\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{\top})^{-1}) \right]=\operatorname{tr}(\mathbb{E}_{{\mathcal{W}}}\left[(\Phi_{\mathcal{W}} \Phi_{\mathcal{W}}^{\top})^{-1}\right])\eqqcolon C_{2}$ $$

which is finite whenever $\mathbb{E}_{{\mathcal{W}}}\left[(\Phi_{\mathcal{W}}\Phi_{\mathcal{W}}^{\top})^
{-1}\right]$ is finite, i.e. LABEL:ass:inverse-third-moment-constrained holds.
∎

Using LABEL:lmm:bound-difference-ridge-regressors and LABEL:lmm:bound-difference-random-feature-ridge-regressors we can now show that the difference between the infinite ensemble where each model has ridge $\lambda$ and the infinite single model with ridge $\lambda$ is Lipschtiz-continuous in $\lambda$ for $\lambda\geq 0$:

###### Proof.

We bound difference $\left||\bar{h}^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})-h^{(RR)}_{\infty,
\lambda^{\prime}}(x^{*})|-|\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})-h^{(RR)}_{
\infty,\lambda}(x^{*})|\right|$ by using first the inverse, then the normal triangle inequality:

$$ $\displaystyle\left||\bar{h}^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})-h^{(RR)}_{ \infty,\lambda^{\prime}}(x^{*})|-|\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})-h^{( RR)}_{\infty,\lambda}(x^{*})|\right|$ $\displaystyle\leq|\bar{h}^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})-\bar{h}^{(RR )}_{\infty,\lambda}(x^{*})+h^{(RR)}_{\infty,\lambda}(x^{*})-h^{(RR)}_{\infty, \lambda^{\prime}}(x^{*})|$ $\displaystyle\leq|\bar{h}^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})-\bar{h}^{(RR )}_{\infty,\lambda}(x^{*})|+|h^{(RR)}_{\infty,\lambda}(x^{*})-h^{(RR)}_{\infty ,\lambda^{\prime}}(x^{*})|$ $$

Using the bound from LABEL:lmm:bound-difference-ridge-regressors and LABEL:lmm:bound-difference-random-feature-ridge-regressors (and summarizing the the corresponding constants as $c_{1}$ and $c_{2}$) we can bound this by:

$$ $|\bar{h}^{(RR)}_{\infty,\lambda^{\prime}}(x^{*})-h^{(RR)}_{\infty,\lambda^{ \prime}}(x^{*})|-|\bar{h}^{(RR)}_{\infty,\lambda}(x^{*})-h^{(RR)}_{\infty, \lambda}(x^{*})|\leq c_{1}|\lambda^{\prime}-\lambda|+c_{2}|\lambda^{\prime}-\lambda|$ $$

Thus we have Lipschitz-continuity in $\lambda$ for $\lambda\geq 0$.

The Lipschitz constant is independent of $x^{*}$ for ${\mathcal{X}}$ compact since the Lipschitz constants from LABEL:lmm:bound-difference-ridge-regressors and LABEL:lmm:bound-difference-random-feature-ridge-regressors depend on $x^{*}$ in a continuous fashion.
∎

Note that an equivalent argument in combination with Jacot et al. ()[Proposition 4.2], i.e. $\tilde{\lambda}\leq\tfrac{\gamma}{\gamma-1}\lambda$, directly gives the Lipschitz-continuity in $\lambda$ for $\lambda\geq 0$ for the difference between the infinite ensemble and the infinite-width single model with effective ridge in the overparameterized regime.

###### Proof.

###### Proof.

###### Proof.

## Appendix D Underparameterized Ensembles

Here, we offer a proof that infinite, unregularized, underparameterized RF ensembles
are equivalent to kernel ridge regression under a transformed kernel function.
We emphasize the difference from the overparameterized case—the central focus of our paper—in which the infinite ensemble is equivalent to a ridgeless kernel regressor.
Thus, underparameterized ensembles induce regularization, while overparameterized ensembles do not.

Other works have explored the ridge behavior of underparameterized RF ensembles
(Kabán, ; Thanei et al., ; Bach, );
however, these works often focus on an equivalence in generalization error
whereas we establish a pointwise equivalence.
To the best of our knowledge, the following result is novel:

###### Lemma 1 .

If the expected orthogonal projection matrix
$\mathbb{E}_{\tilde{W}}\left[R^{\top}\tilde{W}\left(W^{\top}RR^{\top}W\right)^{
-1}\tilde{W}^{\top}R\right]$
is well defined, and a contraction (i.e., singular values strictly less than 1),
then the infinite underparameterized RF ensemble $\bar{h}^{(LN)}_{\infty}(x^{*})$
is equivalent to kernel ridge regression under some kernel function $\tilde{k}(\cdot,\cdot)$.

###### Proof.

When $D<N$, the infinite ridgeless RF ensemble is given by

$$ $\displaystyle\bar{h}^{(LN)}_{\infty}(x^{*})$ $\displaystyle=\mathbb{E}_{{\mathcal{W}}}\left[\tfrac{1}{D}\operatorname{{ \textstyle\sum}}_{j=1}^{D}\phi(\omega_{j},x^{*})\left(\tfrac{1}{D}\Phi_{ \mathcal{W}}^{\top}\Phi_{\mathcal{W}}\right)^{-1}\Phi_{\mathcal{W}}^{\top} \right]y$ $\displaystyle=\mathbb{E}_{W,w_{\bot}}\left[\left(r_{\bot}w_{\bot}^{\top}+c^{ \top}W\right)\left(W^{\top}RR^{\top}W\right)^{-1}W^{\top}\right]Ry,$ (11) $$

where $W,w_{\bot},r_{\bot},c,R$ are as defined in Sec. 2.
Defining the following block matrices:

$$ $\tilde{W}=\begin{bmatrix}W\\ w_{\bot}^{\top}\end{bmatrix}\in\mathbb{R}^{(N+1)\times D},\qquad\tilde{R}= \begin{bmatrix}R\\ 0\end{bmatrix}\in\mathbb{R}^{(N+1)\times N},\qquad\tilde{c}=\begin{bmatrix}c\\ r_{\bot}\end{bmatrix}\in\mathbb{R}^{(N+1)},$ $$

we can rewrite Eq. 11 as

$$ $\displaystyle\bar{h}^{(LN)}_{\infty}(x^{*})$ $\displaystyle=\tilde{c}^{\top}\left(\mathbb{E}_{\tilde{W}}\left[\tilde{W}\left (W^{\top}RR^{\top}W\right)^{-1}\tilde{W}^{\top}\right]\right)\tilde{R}y.$ $$

By adding and subtracting $\tilde{R}\tilde{R}^{\top}$ inside the outer parenthesis,
we can massage this expression into kernel ridge regression in a transformed coordinate system:

$$ $\displaystyle\bar{h}^{(LN)}_{\infty}(x^{*})$ $\displaystyle=\tilde{c}^{\top}\left(\tilde{R}\tilde{R}^{\top}+\underbracket{ \left(\mathbb{E}_{\tilde{W}}\left[\tilde{W}\left(W^{\top}RR^{\top}W\right)^{-1 }\tilde{W}^{\top}\right]\right)^{-1}-\tilde{R}\tilde{R}^{\top}}_{:=\tilde{A}} \right)^{-1}\tilde{R}y.$ $\displaystyle=\tilde{c}^{\top}\tilde{A}^{-1}\tilde{R}\left(\tilde{R}^{\top} \tilde{A}^{-1}\tilde{R}+I\right)^{-1}y.$ (12) $$

Applying the Woodbury inversion lemma to $\tilde{A}^{-1}$, we have:

$$ $\displaystyle\tilde{A}^{-1}$ $\displaystyle=\mathbb{E}_{\tilde{W}}\left[\tilde{W}\left(W^{\top}RR^{\top}W \right)^{-1}\tilde{W}^{\top}\right]$ (13) $\displaystyle+\mathbb{E}_{\tilde{W}}\left[\tilde{W}\left(W^{\top}RR^{\top}W \right)^{-1}W^{\top}R\right]\left(I-\mathbb{E}_{W}[P_{W}]\right)^{-1}\mathbb{E }_{\tilde{W}}\left[R^{\top}W\left(W^{\top}RR^{\top}W\right)^{-1}\tilde{W}^{ \top}\right],$ $$

where $P_{W}$ is the (random) orthogonal projection matrix onto the span of the columns of $R^{\top}W$:

$$ $P_{W}=R^{\top}W\left(W^{\top}RR^{\top}W\right)^{-1}W^{\top}R.$ $$

Because $P_{W}$ is an orthogonal projection matrix, we have that $\|P_{W}\|_{2}=1$,
and thus (by Jensen’s inequality) $\|\mathbb{E}_{W}[P_{W}]\|_{2}\leq 1$.
If this inequality is strict so that $I-\mathbb{E}_{W}[P_{W}]$ is invertible,
we have by inspection of Eq. 13 that $\tilde{A}$ is positive definite.
Therefore, the block matrix

$$ $\begin{bmatrix}\tilde{R}^{\top}\\ \tilde{c}^{\top}\end{bmatrix}\tilde{A}^{-1}\begin{bmatrix}\tilde{R}&\tilde{c} \end{bmatrix}=\begin{bmatrix}\tilde{R}^{\top}\tilde{A}^{-1}\tilde{R}&\tilde{R} ^{\top}\tilde{A}^{-1}\tilde{c}\\ \tilde{c}^{\top}\tilde{A}^{-1}\tilde{R}&\tilde{c}^{\top}\tilde{A}^{-1}\tilde{c }\end{bmatrix}$ (14) $$

is also positive definite and thus the realization of some kernel function $\tilde{k}(\cdot,\cdot)$;
i.e.

$$ $\begin{bmatrix}\tilde{R}^{\top}\tilde{A}^{-1}\tilde{R}&\tilde{R}^{\top}\tilde{ A}^{-1}\tilde{c}\\ \tilde{c}^{\top}\tilde{A}^{-1}\tilde{R}&\tilde{c}^{\top}\tilde{A}^{-1}\tilde{c }\end{bmatrix}=\begin{bmatrix}\tilde{k}(x_{1},x_{1})&\cdots&\tilde{k}(x_{1},x_ {N})&\tilde{k}(x_{1},x^{*})\\ \vdots&\ddots&\vdots&\vdots\\ \tilde{k}(x_{N},x_{1})&\cdots&\tilde{k}(x_{N},x_{N})&\tilde{k}(x_{N},x^{*})\\ \tilde{k}(x^{*},x_{1})&\cdots&\tilde{k}(x^{*},x_{N})&\tilde{k}(x^{*},x^{*}) \end{bmatrix}.$ $$

Note that if $\tilde{A}=I$ then by Eq. 1 we recover the original kernel matrix

$$ $\begin{bmatrix}\tilde{R}^{\top}\tilde{R}&\tilde{R}^{\top}\tilde{c}\\ \tilde{c}^{\top}\tilde{R}&\tilde{c}^{\top}\tilde{c}\end{bmatrix}=\begin{ bmatrix}k(x_{1},x_{1})&\cdots&k(x_{1},x_{N})&k(x_{1},x^{*})\\ \vdots&\ddots&\vdots&\vdots\\ k(x_{N},x_{1})&\cdots&k(x_{N},x_{N})&k(x_{N},x^{*})\\ k(x^{*},x_{1})&\cdots&k(x^{*},x_{N})&k(x^{*},x^{*})\end{bmatrix}.$ $$

Thus, the underparameterized ensemble in Eq. 12 simplifies to

$$ $\displaystyle\bar{h}^{(LN)}_{\infty}(x^{*})$ $\displaystyle=\begin{bmatrix}\tilde{k}(x^{*},x_{1})&\cdots&\tilde{k}(x^{*},x_{ N})\end{bmatrix}\left(\begin{bmatrix}\tilde{k}(x_{1},x_{1})&\cdots&\tilde{k}(x _{1},x_{N})\\ \vdots&\ddots&\vdots\\ \tilde{k}(x_{N},x_{1})&\cdots&\tilde{k}(x_{N},x_{N})\end{bmatrix}+I\right)^{-1 }y,$ $$

which is kernel ridge regression with respect to the kernel $\tilde{k}(\cdot,\cdot)$.
∎

###### Lemma 1 .

###### Proof.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: Can Theorem 1 be generalized to account for feature learning? For instance, does the result hold true if 1 step of gradient descent is performed on the first layer weights? Is it possible for the analysis in Theorem 2 to be refined so as to place bounds on the implicit regularization $\tilde \lambda$ of the ensembled model? Also, can the discussion in Section 3.2 be used to evaluate robustness of the ensemble to distribution shifts?
soundness: 3
strengths: The paper is clearly written and provides a good overview of recent papers in the area. The strongest result in the paper is Theorem 1. This theorem establishes equivalence results under Assumption 1, which assumes subexponential distributions on the whitened random features and allows for dependence. This is weaker than the gaussian process assumption required by Jacot et al (2020). The key technical result is Lemma 1, which leverages concentration inequalities for subexponential random variables to handle an error term in the pointwise analysis of Theorem 1. The experiments are clear and support Theorem 1 and its corollaries.
summary: This paper studies the effect of ensembling in random features (RF) models. The basic question in this area is whether ensembling provides any advantages over a single larger model, such as implicit regularization or better generalization. Building on the work of Jacot et al. (2020), the authors demonstrate that in the regression setting, an infinite ensemble average of ridgeless RF models is equivalent to a single infinitely-wide ridgeless RF model. Notably, this result does not rely on the gaussian universality assumption common in the literature. The authors then discuss the positive ridge case, showing that the same pointwise result holds up to an error proportional to lambda. As a consequence, they point out that popular ensemble-based uncertainty measures capture deviations from an infinite-width limit, not the true response, and therefore must be interpreted carefully.
weaknesses: One concern is that Theorem 1 does not offer substantially new insights into ensembling compared to Jacot et al. (2020), which also showed the equivalence of ridgeless ensembles and a ridgeless infinite-width limit, albeit under stronger assumptions.
Also, I believe the paper could be strengthened if it provided more theoretical analysis in the positive ridge setting. The authors prove that the error is Lipschitz in the regularization $\lambda$ in Theorem 2. However, in Jacot et al. (2020), the authors characterize the regularization $\tilde \lambda$ of the ensemble and show that $\tilde \lambda$ is strictly greater than $\lambda$. If the results in Theorem 2 are sharpened, can this implicit regularization be recovered?
Finally, I believe the discussion of uncertainty quantification in Section 3.2 could use further development. In particular, can the results in Theorem 1 or Theorem 2 be used to demonstrate a lack of robustness to distribution shifts, as discussed in Abe et al. (2022)?

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: - Assumption 1.2 (a.s. Positive definite): Can you provide some examples of features distributions for which this assumption is true. For the ReLU and Leaky-ReLU approximations, can you provide the precise rate at which the approximations yield the positive definite matrix. For example, in footnote (2), what is the dependence on $\alpha$ in the condition number of the matrix \sum_i w_i w_i^\top?
- Theorem 1: Could you please explain where in the proof of this theorem would you need the assumption that D >> N? What is the rate of convergence in terms of number of elements in the ensembles and in terms of the number of features of the single model?
- Theorem 1: What are the main challenges in obtaining a finite dimensional version of this theorem? 
- Theorem 1: Can you please instantiate the rates of the theorem for a few different classes of feature vectors which satisfy Assumption 1?
- What classes of features and distributions satisfy assumption 2?
- Theorem 2: what is the exact characterization of the Lipschitz constant \lambda? Can you provide a few examples for this constant for different distributions/features?
- What are the implications of these results for actual deep networks? How well do these results transfer empirically to the deep learning regime? It would be useful to provide some experiments with small real-world datasets like Cifar-10 to study the implications of these results.
soundness: 3
strengths: - Builds on existing theoretical framework to address an important question of whether there is any advantage of training ensembles vs a bigger model.
- The assumptions are clearly stated and the proofs of the theorems are complete.
- The notation and setup are clearly explained
summary: This paper explores the generalization and uncertainty quantification aspects of overparameterized ensembles, particularly in neural networks. It examines how modern overparameterized ensembles, such as those constructed from neural networks, differ from traditional, underparameterized ensembles. Using random feature (RF) regressors, the paper theoretically demonstrates that overparameterized ensembles tend to be equivalent to a single large model in terms of generalization and predictive variance. Key contributions include proving that, under certain conditions, overparameterized ensembles provide no inherent generalization advantage over a single large model, and showing that predictive variance in overparameterized ensembles does not equate to traditional uncertainty measures.
weaknesses: (See Questions for details)
- The paper does not explore empirical implications of the theoretical results 
- The results are preliminary in nature - focus mostly on asymptotic regime and corresponding finite-regime rates are not provided
- Very few examples to illustrate the assumptions and the actual implications of the theorems.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: See weaknesses.
soundness: 2
strengths: 1. The story is clear, and the paper is well-organized.
2. The figures verify the theoretical results, and make the statement persuasive.
summary: This paper provided an asymptotic theoretical analysis for model ensemble in over-parameterization regime.
weaknesses: 1. The result in Theorem 1 (question 1) has been proposed in a work titled "On the Benefits of Over-parameterization for Out-of-Distribution Generalization"(see [1]), which has been posted on March 2024. They have analysed such phenomenon under the same setting. Please add the citation, as well as comparing the difference between your result and the result in this work.
2. Assumption 2 might related to the distribution on data x, please characterise it with respect to x, to make the assumption more clear.

[1] Hao Y, Lin Y, Zou D, et al. On the Benefits of Over-parameterization for Out-of-Distribution Generalization[J]. arXiv preprint arXiv:2403.17592, 2024.

Meta-review: This paper studies the theoretical limitations of ensembles under overparameterization. The reviewers raised a number of concerns, including the limited innovation and insight gained by generalizing the results of  Jacot et al (2020) to more general distributions (what does one gain by generalizing to non-Gaussian? can you be specific?), and the unclear scope of some of the assumptions. Further, a closely related work by Hao et al. (2024) has not been discussed in the original submission. While the authors provided a discussion in the rebuttal, the reviewer who raised it was not able to answer, making it unclear if the discussion was complete and adequate. Given these issues, it will be better for the community if the authors revise their paper and resubmit to a future venue.

Decision: Reject