Learning Curves for Noisy Heterogeneous Feature-Subsampled Ridge Ensembles
Authors: Benjamin Samuel Ruben
Cengiz Pehlevan
Venue: NeurIPS 2023 poster (2023)
Keywords: ridge regression
ensembling methods
TL;DR: We derive learning curves for feature-subsampled ridge ensembles and show that structural heterogeneity acts as an implicit regularizer.

Abstract: Feature bagging is a well-established ensembling method which aims to reduce
prediction variance by combining predictions of many estimators trained on subsets
or projections of features. Here, we develop a theory of feature-bagging in noisy
least-squares ridge ensembles and simplify the resulting learning curves in the special
case of equicorrelated data. Using analytical learning curves, we demonstrate
that subsampling shifts the double-descent peak of a linear predictor. This leads
us to introduce heterogeneous feature ensembling, with estimators built on varying
numbers of feature dimensions, as a computationally efficient method to mitigate
double-descent. Then, we compare the performance of a feature-subsampling
ensemble to a single linear predictor, describing a trade-off between noise amplification
due to subsampling and noise reduction due to ensembling. Our qualitative
insights carry over to linear classifiers applied to image classification tasks with
realistic datasets constructed using a state-of-the-art deep learning feature map.

Learning Curves for Noisy Heterogeneous Feature-Subsampled Ridge Ensembles
Benjamin S. Ruben', Cengiz Pehlevan 2.3,4
Biophysics Graduate Program
Center for Brain Science, JJhhn A. Paulson School of Engineering and Applied Sciences, *Kempner Institute for the Study Natural and Artificial Intelligence, Harvard University
Cambridge, MA 02138
benrubengg. harvard.edu, cpehlevandseas hww.arrr.eee... edu
Abstract
Feature bagging is well-established ensembling method which aims to reduce prediction variance by combining predictions of many estimators trained on subsets projections of features. Here, develop theory of feature-bagging in noisy least-squares ridge ensembles and simplify the resulting learning curves in the spe- cial case of equicorrelated data. Using analytical learning curves, demonstrate that subsampling shifts the double-descent peak of a linear predictor. This leads to introduce heterogeneous feature ensembling, with estimators built varying numbers of feature dimensions, as a computationally efficient method to mitigate double-descent. Then, compare the performance of feature-sufee Seature-subsampling ensemble to a single linear predictor, describing a trade-off between noise ampli- fication due subsampling and noise reduction due ensembling. Our qualitative insights carry over to linear classifiers applied to image classification tasks with realistic datasets constructed using a state-of-the-arr deep learning feature map.
Introduction
Ensembling methods are ubiquitous in machine learning practice [1] A class of ensembling methods (known as attribute bagging (2) or the random subspace method (3] is based on feature subsampling 2 where predictors are independently trained subsets of the features, and their predictions are combined to achieve stronger prediction. The random forest method a popular example B
While commonly used in practice, a theoretical understanding of ensembling via feature subsampling well developed. Here, we provide an analysis of this technique the linear ridge regression setting. Using methods from statistical physics 12] we obtain analytical expressions for typical- case generalization error in linear ridge ensembles (proposition [ and simplify these expressions in the special case of equicorrelated data with isotropic feature noise (proposition [2 The result provides a powerful tool to quickly probe the generalization error of ensembled regression under a rich set of conditions. In section 3 we study the behavior of single feature- feature-subsampling regression model. observe that subsampling shifts the location
predictor's sample-wise double-descent peak 13-15). This motivates section 4 where we study ensembles built on predictors which are heterogencous the number of features they access, as a method to mitigate double-descent. We demonstrate this method's 's efficacy in a realistic image classification task. In section 5 we apply theory to the trade-off between ensembling and subsampling in resource-constraine scttings. We
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

characterize how a variety of factors influence the optimal ensembling strategy, finding particular significance to the level of noise in the predictions made by ensemble members.
In summary, we make the following contributions:
Using the replica trick from statistical physics 8 we derive the generalization error of ensembled least-squares ridge regression a general setting, and simplify the resulting expressions in the tractable special case where features are equicorrelated.
We demonstrate benefits heterogeneous ensembling as robust and computationally efficient regularizer for mitigating descett with analytical theory and in realistic image classification task.
We describe the ensembling-subee suubsbmminggs trade-off in essurce-constrainee settings, and charac- terize the effect of label noise, feature noise, readout noise, regularization, sample size and task structure the optimal ensembling strategy.
Related works: A substantial body of work has elucidated the behavior of linear predictors for a variety of feature maps [14. 30) Several recent works have extended this research to characterize the behavior of ensembled regression using solvable models [24 31 Additional recent works study the performance of ridge ensembles with example-wise subsampling [4 335 and simultane-
ous subsampling of features and examples [2), finding that subsampling behaves as an implicit regularization. Methods from statistical physics have long been used for machine learning theory [10122 26 27 30 36, 37.. Relevant work in this domain include (38) which studied ensembling by
data-subsampling in linear regression.
2 Learning Curves for Ensembled Ridge Regression
We consider noisy ensembled ridge regression in the setting where ensemble members are trained independently on masked versions of the available features. We derive our main analytical formula for generalization error of ensembled linear regression, as well as analytical expressions for generalization error in the special case of equicorrelated features with isotropic noise.
2.1 Problem Setup
Consider a training set D = ("b", Vn of size P. The training examples T5"h E M are drawn from Gaussian distribution with Gaussian feature noise: vhtt r/24 + at, where ap" NO(, Σ.) and T (0, Eo). Data and noise are drawn i.i.d. so that IE' [2blab" mnd Eg. Labels are generated from a noisy teacher function y" VM w obaa + ed where e (O00 (²). Label noises are drawn i.i.d. so that E[e"te] = 5
8nl Eg and E o
We seek analyze the quality of predictions which averaged over an ensemble of ridge regression models, each with access to subset of the features. We consider l: linear predictors with weights w, R.N, r 1...., k. Critically. we allow N, a Np for r * r', which allows us to introduce structural heterogeneity into the ensemble of predictors. A forward pass of the model is given as:
1 Σ ((()),
1 VN.
5(b)
f+(1))
0 Ar + σ) +
The model's prediction f(+) is an average over k linear predictors. The "measurement matrices" A, E IR N act as linear masks restricting the information about the features available to each member of the ensemble. Subsampling may be implemented by choosing the rows of each A, to coincide with the rows of the identity matrix the row indices corresponding to indices of the sampled features. The feature noise a ~ (000 En) and the readout noises &, ~ A'(0, n/²,, are drawn independently at the execution of each forward pass of the model. Note that while the feature noise is shared across the ensemble, readout noise is drawn independently for each readout: [[G-Er] m,

---

The weight vectors are trained separately in order to minimize an ordinary least-squares loss function with ridge regularization:
w. arg min
Σ
w. Ar + 8"t + 1r VN,
1
(2)
-
Here { represents the readout noise which is present during training, and independently drawn: EH (O, ) E [[ = , As measure of model performance, consider the generaliza- tion enror, given by the neaa-squared-error (MSE) on ensemble-averagen prediction:
[
1
(D) (E.)
w
2/
M
Here, the expectation over the data distribution and noise: &b ~ (O E.), o ~ Eo), & (O n² The generalization error depends on the particular realization of the dataset D through the learned weights (in") }. We may decompose the generalization error as follows:
k
(D)
Σ mx'==
Err
E
1
Virip' A A (1 w²
Enp (D) M
A, w.
w²)
E,
Vpr
(5)
0 A, En A,, + Mö
Computing the generalization error of the model is then matter of calculating E,,, in the cases where r = r' and r * concentrates over randomly drawn datasets D.
# r'. In the asymptotic limit we consider, we expect that the generalization error
2.2 Main Result
calculate the generalization error using the replica trick from statistical physics, and present the calculation in Appendix The result of our calculation is stated in proposition 1.
Proposition 1. Consider the ensembled ridge regression problem described in Section 2.1 Consider
the asymptotic limit where M, P, {Nr}
while the ratios 37 and Uer N.
1,....
remain fixed. Define the following quantities:
Σ
A, E. A
(6)
Verelp'
G, IN, +q, ETr
(7)
'er M(2, + (r) gr))) tr G.
8
Then the terms of the average generalization error (eq. 5 may be written as:
EEpr (D)) Trr 2 + 8p
M w= E
1 7rr
Vrrpt
w
Σ
A,A,G.G A, dr A.G. A,
E,w" (9)
M(1 Yert)
Vrr
Vy'p"
w
Σ, Σ,A, G G ES, G. E,w
M(1 )
Ver
where the pairs of order parameters {gr, 2r) for r I,... K, satisfy the following self-consistent saddle-point equations
1
4r Ap 4r
gr =
M tr G, Emr
(10)
3

---

b
a
10
10
Eg
A
E
1e 06
le 06
100
0.001
100
0.001
1
0.0
0.5
1.0
0.0
0.5
1.0
a = P/M
a = P/M
Figure 1: Comparison of numerical and theoretical learning curves for ensembled linear regression.
Circles represent numerical results averaged over 100 trials; lines indicate theoretical predictions. Error bars represent the standard error of the mean but are often smaller than the markers. (a) Testing of proposition I with M 2000, [El .li-j [Eol,y to (0.3))1--1, <= 0.1, and all 7/, 0.2 and X (see legend). k 3 linear predictors access fixed, randomly selected (with replacement) subsets of the features with fractional sizes V/r = 0.2, 0.4, 0.6. Fixed ground-trute -rrnc weights w are drawn from an isotropic Gaussian distribution. (b) Testing of proposition with M 5000.
= 0.6, 3 0.1. { = 0.1, all 7)- 0.1. and all A (see legend). Ground truth weights sampled as in eq. I with p 0.3. Feature subsets accessed by each readout are mutually exclusive (inset) with fractional sizes Vrr 0.1, 0.5.
Proof. We calculate the terms in the generalization non-rigorous method from the statistical physics of disordered systems. The full derivation may be found in the Appendix (F) When the matrices A, (E, Σg) A, r A ...... k have bounded spectra,
using the replica trick, standard but
this result may be obtained by extending the results of (31] to include readout noise, as shown in Appendix G
make the following remarks:
Remark 1. Implicit in this theorem is the assumption that the relevant matrix contractions and traces which appear in the generalization error (eq. and the surrounding definitions tend a well-defined limit which remains (1)) as M - DC.
Remark 2. This result applies for any (well-behaved) linear masks {A,). We will focus the where each A, implements subsampling of an extensive fraction vrr of the features.
Remark 3. When & 1, our result reduces to the generalization error of a single ridge regression model, as studied in refs. 33 39].
Remark 4. We include "readout noise" which independently corrupts the predictions of each ensemble member. This models sources of variation between ensemble members not otherwise accounted for. For example, ensembles of deep networks will vary due to random initialization of parameters 24
31 36] Readout noise is more directly present in physical neural networks, such an analog biological neural circuits 41 due to their inherent stochasticity.
neural networks
40
In Figure p. we confirm the result of the general calculation by comparing with numerical experi- ments using a synthetic dataset with M = 2000 highly structured features (see caption for details). k 3 readouts see random, fixed subsets of features. Theory curves are calculated by solving the fixed-point equations 10 nomerically for the chosen E., Eg and {A, 1 then evaluating eq.
2.3 Equicorrelated Data
Our general result allows the freedom to tune many important parameters the learning problem: the correlation structure of the dataset, the number of ensemble members, the scales of noise, etc. However, the derived expressions are rather opaque. order to better understand the phenomena captured by these expressions, examine the following special case:

---

Proposition 2. In the setting of section 2.1 and proposition 7
consider the following special case:
w¹ = V1 p w + plM
(11)
w/ NO IM)
(12)
Σ, = s M clMiM]
(13)
3o = I2
(14)
with [0, 1.P E [-1,1]. Label and readout noises S,7h 0 are permitted. Here P1
IM 1 MIMIM M is projection matrix which removes the component of w which is parallel to 1 The matrices {A, have rows consisting of distinct one-hot vectors that each of the k readouts has access to subset of Nr = -rrM features. For r + r', denote by Thps the number of neurons sampled by both A, and A, and let Vyr' TTr /M remain fixed M OC.
M
Define the following quantities:
4.
a Ver S.S.
a s=
S,
(15)
Vyr aq,
The terms of the decomposed generalization error may then be written:
(Ere') D.mg 1 - ((1 - 12) + +1211r)) + Trr (² 8ppi
(16)
1 - 1 'ir''
where we have defined
s(1 ((l-c)(( - s(1 - c)urs s(1 - e) Sp + as(1 - -c)rr' S.Sww)
(17)
<
(18)
ppl
10
c
and where (gr, } may be obtained analytically as the solution (with or > 0) to:
OVvr
Or Ver + agr
4.
(19)
N + 4
In the "ridgeless" limit where all Dy + we may make the following simplifications:
2a
40
Sr + a (a + Vrr + la -vrr))
Trr
(ax + ver + la - Uprl) (a Uerr + lax Myere)
(20)
Proof. Simplifying the fixed-point equations and generalization error formulas in this special case is an exercise in linear algebra. The main tools used are the fact that the data distribution is isotropic in the features so that the form of E,y and Σ,, depend only on the subsampling and overlap fractions Vrr: Mr',' Prr. To aid in computing the necessary matrix contractions we developed custom Mathematica package which handles block matrices of symbolic dimension, with blocks containing matrices of the form M ei + c1l This package and the Mathematica notebook used to derive these results are available online (see Appendix B
formula 42 and the
In this tractable special case, [0, 1] is a parameter which tunes the strength of correlations between features of the data. When c = O, features are independent, and when e = the features are always equivalent. 8 sets the overall scale of the features and the "Data-Tans Task alignment" tunes the alignment of the ground truth weights with the special direction in the covariance matrix (analogous to "task-model alignment [14 2. A table of parameters is provided in Appendix AA] In Figure Dp. we test these results by comparing the theoretical expressions for generalization error with the results of numerical experiments, finding perfect agreement.
With an analytical formula for the generalization error, we can compute the optimal regularization parameters X, which minimize the generalization error. These may, in general, depend on both and the sample size a. Rather than minimizing the error of the ensemble, we may minimize the generalization error of predictions made by the ensemble members independently. We find that this "locally optimal" regularization, denoted A independent of o, generalizing results from 14 to correlated data distributions (see Appendix H.3
5

---

a
b
ii
iii
102
= 0
A .0001
40
2 .01
13 10
Test Eror
100
Teat Error
104
Text Error
x 10
0.6
101
10
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
P/M
P/M
P/M
Figure Subsampling alters the location of the double-descent peak of a linear predictor. (a) Illustrations of subsampled linear predictors with varying subsampling fraction (b) Comparison between experiment and theory for subsampling linear regression on equicorrelated datasets. We choose task parameters as in proposition 2 with
0, 1, and (i) X 0,
learning curves are for a single linear predictor A with subsampling fraction shown in legend. Circles show results of numerical experiment. Lines are analytical prediction.
(ii) = 10-4, (iii) À 10-2
3 Subsampling shifts the double-descent peak of a linear predic- tor
Consider a single linear regressor (k: 1) which connects to subset of vM features in the equicorrelated data setting of proposition Also setting c = 0, and 7) = w = 0 and taking the limit À the generalization error reads:
(1 v) + (a--v)1 v) 1 - 5
S:: if i }
(Eg)
(21)
We thus see that double descent can arise from two possible sources of variance: explicit label noise (if Ç 0) or implicit label noise induced by feature subsampling (v < 1). F.g ~ (a - = generalization eиror diverges when sample size is equal to the number of sampled features. Intuitively, this occurs because subsampling changes the number of parameters of the regression model, and thus its interpolation threshold. To demonstrate this, plot the learning curves for subsampled linear regression on equicorrelated data Figure 2 At small finite ridge the test error no longer diverges when v, but still displays distinctive peak.
4 Heterogeneous connectivity mitigates double-descent
Double-descent over-fitting to noise in the training near a model's interpolation threshold poses a serious risk in practical machine-lcarning applications [22]. Cross-validating the regularization strength against the training set the canonical approach to avoiding double-descent 17
but in practice requires computationally expensive parameter sweep and prior knowledge of the task. In situations where computational resources are limited or hyperparameters are fixed prior to specification of the task, it is natural to seek an alternative solution. Considering again the plots Figure (b). we observe that at any value of a, the double-descent peak can be avoided with an acceptable choice of the subsampling fraction P. This suggests another strategy to mitigate double descent: heterogeneous ensembling. Ensembling over predictors with a heterogeneous distribution of interpolation thresholds, we may expect that when predictor fails due over-fitting, the other members the ensemble compensate with accurate predictions.
In Figure 3 we show that heterogeneous ensembling can guard against double-descent. We define two ensembling strategies: homogeneous ensembling, each of readouts connects a fraction vrr 1 1/k features. In heterogeneous ensembling, the number of features connected by each of the & readouts are drawn from Gamma distribution (R) with mean 1 1/k: and standard deviation a (see Fig. 3p) then re-scaled to sum to (see Appendix Cfro details). All feature subsets are mutually exclusive (Myp () for r of r'). Homogeneous and heterogeneous ensembling are illustrated for k in Figs. a.i and a.i respectively. We test this hypothesis using eq. [6 in 3c. At small regularization

---

-001,, we find that heterogeneity of the distribution of subsampling fractions (a > 0) lowers the double-descent peak of an ensemble of linear predictors, while at larger regularization (A 0.1), there is little difference between homogeneous and heterogeneous learning curves. The asymptotic (a -> 20) error is unaffected by the presence of heterogeneity in the degrees of connectivity, which can be seen as the coincidence of the triangular markers in Fig. Ep. as well from the or - oo limit of eq. [6 (see Appendix H. Fig. 3 also shows the learning curve of a single linear predictor with no feature subsampling and optimal regularization. We see that the feature-subsampling ensemble appreciably outperforms the fully-connected model when c = 0.8 and `j = 0.5, suggesting the important roles of data correlations and readout noise in determining the optimal readout strategy. These roles are further explored in section 5 and fig 4
We also test the effect of heterogeneous ensembling in the realistic classification task. Specifically, we train ensembles of linear classifiers to predict the labels of imagenet [44] images corresponding to 10 different dog breeds (the "Imagewoof" task (45) from their top-hidden-layer representations in a pre-trained ResNext deep network [46] (see Appendix E for details). We characterize the statisties of the resulting M- 2048-dimensional feature set in Fig. ST] This "ResNext-Imaaewooof" classification task has multiple features which make it amenable to learning with a subturl-uubsampling ensemble. First, the ResNext features have a high degree of redundancy 47), allowing classífication to be performed accurately using only a fraction of the available features (see Fig. 3 and (STp). Second, when classifications of multiple predictors are combined by majority vote, there is natural upper bound on the influence of a single erring ensemble member (unlike in regression where predictions can diverge). Calculating learning curves for the imagewoof classification task using homogeneous ensembles, we see sharp double-descent peaks in an ensemble of size k when P = M/k (Fig. i... Using a heterogeneous ensemble mitigates this catastrophic over-fitting, leading to monotonically decreasing error without regularization (Fig. Bii. A single linear predictor with tuned regularization of 0.1 performs only marginally better than the heterogeneous feature- subsampling ensemble with k = 16 or k = 32. This suggests heterogeneous ensembling can be an effective alternative regularization in real-world classification tasks using pre-trained deep learning feature maps.
Note that training Ceature-subsampliin ensemble also benefits from improved computational complexity. Training an estimator of dimension N, involves, in the worst case, inverting an N, x N, matrix, which requires ((N3) operations. Setting N, M/K, see that the number of operations required to train an ensemble of l predictors scales as (kh
5 Correlations, Noise, and Task Structure Dictate the Trade-off
In resource-constrainet settings, one must decide between training a single large predictor or an ensemble of smaller predictors. When the number of weights is constrained, ensembling may benefit generalization by averaging over multiple predictions, but at the expense of each prediction incor- porating fewer features. Intuitively, the presence of correlations between features limits the penalty incurred by subsampling, as measurements from a subset of features will also confer information about the unsampled features. The equicorrelated data model of proposition 2 permits a solvable toy model for these competing effects. We consider the special case of ensembling over k readouts, each connecting the same fraction Vrr = V = 1/k of all features. For simplicity, we set Vrr' = O for r * r'. We asses the learning curves of this toy model in both the ridgeless limit À 0 where double-descent has a large effect on test error, and at 'locally optimal' regularization À = X for which double-descent is eliminated. In these special cases, one can write the generalization error in the following forms (see Appendix HH4 for derivation):
E.lksss.cn....p...... Ey s, s.lk,,,c,,,.u a, = X*) = s(1 - -)- (k, H, W,2) w.w.c.mww........
0) =0)=s11-00 -1- () H,W.Z)
(22)
(23)
where we have defined the effective noise-to-signal ratios:
H
W
Z
(24)
s(1
s(1 c)
s(1 - c)
7

---

Therefore, given fixed parameters s, p. a, the value k" which minimizes error depends on the noise scales, s, and e only through the ratios H. W and Z:
K3 H,,W, p, arg min F.g (k)
arg min (k; p. o, H, H,W.Z)
(25)
KEN
REN
A x- H. W, Z,
p. Z.W.Z.p.W.
arg min Ey (R:
p. arg min Ex pk, H, W, KEN
(26)
KEN
In Fig. 4a, we plot these reduced errors curves E, E" as a function of a for varying ensemble sizes k and reduced readout noise scales H. At zero regularization learning curves diverge at their interpolation threshold. At locally optimal regularization X = X", learning curves decrease monotonically with sample size. Increasing readout noise H raises generalization error more sharply for smaller k. In Fig. 4p we plot the optimal k- in various two-dimensional slices of parameter space in which p is fixed and W = = while a and H vary. The resulting phase diagrams may be divided into three regions. In the signal signal-domintted phase single fully-connedd readout is optimal (k* = 1). In an intermediate phase. 1 < k* < oo minimizes error. And in a noise-dominated phase have determined an analytical expression for the boundary
k' oc. At zero regularization, between the intermediate and noise-dominated phases (see Appendix I..4 and dotted lines in Figs 4 Ab.,dd.. The signal-dominated, intermediate, and noise-dominated phases persist when À = X", removing the effects of double descent. In all panels, an increase in H causes an increase in k This can occur because of decrease in the signal-to-readout noise ratio s/ s/1/² through an increase in the correlation strength c. An increase in also leads to an increase in k", indicating that ensembling is more effective for easier tasks. Figs f..d show analogous phase diagrams where W Z varied. Signal- Signal-dominttdd, intermediate. and noise-dominated regimes are visible in the resulting phase diagrams at zero regularization. However, when optimal regularization always optimal. The presence of regions where k* > can thus be attributed to double-descent at sub-optimal regularization or to the presence of readout noise which is independent across predictors. We chart the parameter-space of the reduced errors and optimal ensemble size k- extensively in Appendix[]
used, ' is
[[] We plot learning curves for the "ResNext- Inmrextommaewoof"" ensembled linear classification task with varying strength of readout noise in Fig. 4e. and phase diagrams of optimal ensemble size k and Appendix E.4.3 for
Fig. 4 finding similar behavior to the toy model. See Figs. S3 S4 S5 further discussion.
6 Conclusion
In this paper, we provided a theory of feature-subsampled linear ridge regression. We identified the special case in which features of the data "equicorrelated" minimal toy model to explore the combined effects of subsampling, ensembling, and different types of noise on generalization error. The resulting learning curves displayed two potentially useful phenomena.
First, we demonstrated that heterogeneous ensembling can mitigate over-fitting, reducing or eliminat- ing the double-descent peak of an under-regularized model. In most machine learning applications, the size of the dataset is known at the outset and suitable regularization may be determined to mitigate double descent, either by selecting highly over wer-parameterized model [22 or by cross-validation techniques (see for example (17).. However, in contexts where a single network architecture is de- signed for an unknown task or variety of tasks with varying structure and noise levels, heterogencous ensembling may be used to smooth out the perils of double-descent.
Next, we described a trade-off between noise reduction due to ensembling and noise amplification due to subsampling in a resource-constrainee setting where the total number of weights fixed. Our analysis suggests that ensembling is particularly useful in neural networks with an inherent noise. Physical neural networks, such as analog neural networks(40) and biological neural circuits 44 present such a resource-constrainel environments where intrinsic noise is significant issue.
Much work remains to achieve a full understanding of the interactions between data correlations. readout noise, and ensembling. In this work, we have given a thorough treatment of the convenient special case where features are equicorrelated. Future work should analyze subsampling and en- sembling for codes with realistic correlation structure, such as the pawer-law spectra (see Fig. ST 27
49 and sparse activation patterns [50]
8

---

a
0
10
0.38
Homogeneous
Fully Connected
10
0-22
0.25
ok 0.5
ak
C
11
iii
iv
== 0.001
= 0.001
2.5
0.001
À 0.001
2.0
C-0
0_0
Test Error
0.4
0.5
1.5
1.5
0.3
1.0
0.15
0.5
0.2
0.0
0.1
0.2
0.3
0.0
0.2
0.0
0.1
0.2
0.3
0.0
0.1
0.2
0.3
Vi
vii
viii
0.20
1.2
0.45
Test Error
1.0
0.40
0.9
0.18
0.35
A.1.1
C=0.8
0.30
0,5
A=0.1
0.16
A=
A=0.1
0.25
0.7
9=0.5
0.20
0.14
0.0
0.1
0.2
0.3
0.0
0.1
0.2
0.3
0.0
0.2
0.3
ii
0.0
0.1
0.2
0.3
P/M
P/M
P/M
P/M
d
P= 1000 M 2048
e
Homogeneous
Heterogeneous
0.0e+00
10a
10a
1.0e-03
Test Error
1.0e-02
Test Error
Test Error
1.0c-01
1.0e+000
1. 0.1
10-1
10-1
10-1
1.0e+02
1.0e
1.0e+04
10
10a
101
10a
10a
101
102
102
N/M
Figure 3: Heterogeneous ensembling mitigates double-descent. (a) We compare (i) homogencous ensembling, in which & readouts connect the same fraction 1/k of features, and (ii) hetero- geneous ensembling (b) In heterogeneous ensembling subsampling fractions are drawn i.i.d. from TR.o(v), shown here for i 10, then re-scaled sum 1. (c) Generalization Error Curves for Homogeneous and Heterogeneous ensembling with & = 10, = 0, p 0.3 and indicated values of X, c, and 7). Blue: homogeneous subsampling. Red, green, and purple show heterogeneous subsampling with = 0.25/k: 0.5/k, 1/k respectively. Dashed lines show learning curves for particular realizations
(V11,... VKR Solid curves show the average over realizations, Gray shows the learning curve for single linear readout with = and optimal regularization (eq. 193). Triangular marks show the asymptotic generalization error (a - oo). with downward-pon poiattdn--ppnngggggo gray triangles indicating asymptotic error of zero. (d,e) Generalization error of linear classifiers applied to the imagewoof dataset with ResNext features averaged over 100 trials. (d) subsampling fraction and regularization À (legend). (e) Generalization error of (i) homogeneous and (ii) heterogeneous (with o = 0.75 7k) ensembles classifiers. Legend indicates is values. = except for gray curves, where À 0.1
100, k varying

---

= p= 0.0.0
A -=0,-=00
p= =
p= -
2.0
100
10
109
90
Eg/s(1
I 1.0
107
70
10
1.0
2.0
3.0
0.0
10
2.0
3.0
= 0. 0, 1.0
p=
0.6, A ==
p
0.6, À
50
2,0
Eg/s(1 c)
101
I 1.0
20
10"
10"1
10
107
0.0
1.0
2.0
3.0 m
a
d
0.6, À
p 0.6, À ==
p
0.6, À
p 0.6, À =
2.0
W
1.0
0.0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3,0
a
i
10a
0
K
A
0
k
Mean Error
1
200
2
0.75
100
10-1
3
0.5
30
4
0.25
5
0.0
25
ii
0.05
20
10
15
X 0.1
15
20
10
Mean Error
25
0.75
5
0.5
4
100
1
À
0
0.4
0.25
3
200
2
0.05
0.41
0.77
1.14
0.0
0.5
1.5
P/M
Figure 4: Task parameters dictate the ensembling-subsampingg trade-off: (a-d) In the setting of proposition 2 in the special case where all u,r I: 5,,, so that feature subsets are mutually exclusive and the total number of weights is conserved. (a) We plot the reduced generalization errors & (for
with O, using eq. 22 and Ex (for À = using eq. 23 of linear ridge ensembles of varying size and H/ 0, (values indicated above plots). Grey lines indicate & dashed black lines + 20, and intermediate values by the colorbar. (b) We plot optimal ensemble size k- (eqs. 25 26 in the parameter space of sample size and reduced readout noise scale H setting W 2 = and white indicates A-- = oc, with intermediate values given by the colorbar. Appended vertical bars show - Dotted black lines show the analytical boundary between the intermediate and noise-dominated phases given by eq. [14 (c) optimal readout phase diagrams in (b) but showing W -ddependeecee with H 2 O. (d) optimal readout A- phase diagrams as in (b) but showing 7-dependence with = W = O. (e) Learning curves for saauure-ubbsampling
Grey indicates
ensembles of linear classifiers combined using majority vote rule the imagewoof classification task (see Appendix E As in (a-d) we set 1py t Error is calculated as the probability of incorrectly classifying a example. À and 1) values are indicated in each panel. (f) Numerical phase diagrams showing the value of & which minimizes test error in the parameter space of sample size P and readout noise scale 7). with regularization (i) (pseudoinverse rule) (ii) X 0.1.
10

---

Acknowledgements
CP is supported by NSF Award RMS-21341577 NSF CAREER Award IIS IIS-2239780, and a Sloan Research Fellowship. This work has been made possible in part by gift from the Chan Zuckerberg Initiative Foundation to establish the Kempner Institute for the Study of Natural and Artificial Intelligence. BSR was also supported by the National Institutes of Health Molecular Biophysics Training Grant NIH/ NIGMS T32 GM008313. We thank Jacob Zavatone- Veth and Blake Bordelon for thoughtful discussion and comments this manuscript.
References
[1] Gautam Kunapuli. Ensemble Methods for Machine Learning. Simon and Schuster, 2023.
[2] Robert Bryll, Ricardo Gutierrez-Osuna and Francis Quek. Attribute bagging: improving accu- racy of classifier ensembles by using random feature subsets. Puttern recognition, 12616 1302, 2003.
[3] Tin Kam Ho. The random subspace method for constructing decision forests. IEEE transactions on pattern analysis and machine intelligence, 20(8):832-844 1998.
[4] Yali Amit and Donald Geman. Shape quantization and recognition with randomized trees. Neural computation, 9(7):1545-1588, 1997.
[5] Gilles Louppe and Pierre Geurts. Ensembles on random patches. In Machine Learning and Knowledge Discovery in Databases: European Conference. ECML PKDD 2012, Bristol, UK, September 24-28. 2012. Proceedings, Part 23, pages 346-361. Springer, 2012.
[6] Marina Skurichina and Robert P. W. Duin. Bagging, boosting and the random subspace method for linear classifiers. Pattern Analysis & Applications, 5(2):121-135, June 2002.
[7] Tin Kam Ho. Random decision forests. In Proceedings of 3rd international conference on document analysis and recognition, volume pages 278-282, IEEE, 1995.
[8] Mare Mézard, Giorgio Parisi, and Miguel Angel Virasoro. Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications. World Scientific Publishing Company, 1987.
[9] M Mezard, G Parisi, and M Virasoro. Spin Glass Theory and Beyond. WORLD SCIENTIFIC, November 1986.
[10] Hyunjume Sebastian Seung, Haim Sompolinsky, and Naftali Tishby. Statistical mechanics of learning from examples. Physical review A, 4568):6056, 1992.
[11] Andreas Engel and Christian Van den Broeck. Statistical mechanics of learning. Cambridge University Press, 2001.
[12] Yasaman Bahri, Jonathan Kadmon, Jeffrey Pennington, Sam S Schoenholz, Jascha Sohl- Dickstein, and Surya Ganguli. Statistical mechanics of deep learning. Annual Review of Condensed Matter Physics, 11:501-528, 2020.
[13] Preetum Nakkiran. More data can hurt linear regression: Sample-wise double descent. arXiv preprint arXiv:1912.07242.. 2019.
[14] Abdulkadir Canatar, Blake Bordelon, and Cengiz Pehlevan. Spectral bias and task-model alignment explain generalization in kernel regression and infinitely wide neural networks. Nature communications, 12(1):2914, 2021.
[15] Fatih Furkan Yilmaz and Reinhard Heckel. Regularization-wiie wise double descent: Why it occurs and how to eliminate In 2022 IEEE International Symposium on Information Theory (ISIT), pages 426-431, IEEE, 2022.

---

[16] Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J. Tibshirani. Surprises in high- dimensional ridgeless least squares interpolation. The Annals of Statistics, 50(2):949 986, 2022.
[17] Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman. The elements of statistical learning: data mining, inference, and prediction, volume 2. Springer, 2009.
[18] Jason W. Rocks and Pankaj Mehta. Bias-variance decomposition of overparameterized regres- sion with random linear features. Physical Review E. 106:025304, Aug 2022.
[19] Jason W. Rocks and Pankaj Mehta. Memorizing without overfitting: Bias, variance, and interpolation in overparameterized models. Physical Review Research, 4:013201, Mar 2022.
[20] Song Mei and Andrea Montanari. The generalization error of random features regression: Precise asymptotics and the double descent curve. Communications on Pure and Applied Mathematics, 75, 2019.
[21] Peter L. Bartlett, Philip M. Long, Gábor Lugosi, and Alexander Tsigler. Benign overfitting in linear regression. Proceedings of the National Academy of Sciences, 107643) 2020.
[22] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal. Reconciling modern machine- learning practice and the classical varianer trade-off. Proceedings of the National Academy of Sciences,
2019.
[23] Hong Hu and Yue M. Lu. Universality laws for high-dimensional learning with random features. IEEE Transactions on Information Theory, 69(3):1932-199 1964, 2023.
[24] Stéphane D' Ascoli, Maria Refinetti, Giulio Biroli, and Florent Krzakala. Double trouble in double descent: Bias and variance(s) in the lazy regime. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 2280-2290. PMLR, 13-18 Jul 2020.
[25] Ben Adlam and Jeffrey Pennington. Understanding double descent requires a fine-grained bias- variance decomposition. In Advances in Neural Information Processing Systems, volume 33, pages 11022-11032 2020.
[26] Jacob A. Zavatone-Vent Veth, William L. Tong, and Cengiz Pehlevan. Contrasting random and learned features in deep bayesian linear regression, 2022.
[27] Blake Bordelon, Abdulkadir Canatar, and Cengiz Pehlevan. Spectrum dependent learning curves in kernel regression and wide neural networks. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 1024-1034. PMLR, 13-18 Jul 2020.
(28) James B. Simon, Madeline Dickens, Dhruva Karkada, and Michael R. Dewees The Eigen- learning framework: A conservation law perspective on kernel regression and wide neural networks. arXiv, 2022.
[29] Francis Bach. High- High-dimensional analysis of double descent for linear regression with random projections, 2023.
[30] Jacob A. Zavatone-Vert and Cengiz Pehlevan. Learning curves for deep structured gaussian feature models, 2023.
[31] Bruno Loureiro, Cedric Gerbelot, Maria Refinetti, Gabriele Sicuro, and Florent Krzakala. Fluctuations, bias, variance & ensemble of learners: Exact asymptotics for convex losses in high-dimension. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 14283-14314. PMLR, 17-23 Jul 2022.
12

---

[32] Daniel LeJeune, Hamid Javadi, and Richard Baraniuk. The implicit regularization of ordinary least squares ensembles. Silvia Chiappa and Roberto Calandra, editors, Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 3525-35 3535. PMIR, 26-28 Aug 2020.
[33] Pratik Patil and Daniel LeJeune. Asymptotically free sketched ridge ensembles: Risks, cross- validation, and tuning, 2023.
[34] Jin-Hong Du, Pratik Patil, and Arun Kumar Kuchibhotla. Subsample ridge ensembles: Equiva- lences and generalized cross- validation, 2023.
[35] Pratik Patil, Jin-Hong Du, and Arun Kumar Kuchibhotla. Bagging overparameterized learning: Risk characterization and risk monotonization, 2022.
[36] Alexander Atanasov, Blake Bordelon, Sabarish Sainathan, and Cengiz Pehlevan. The onset of variance-limited behavior for networks in the lazy and rich regimes. In The Eleventh International Conference on Learning Representations, 2023.
[37] Abdulkadir Canatar, Blake Bordelon, and Cengiz Pehlevan. Out-of-distribure generalization in kernel regression. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages I 12600-12612 Curran Associates, Inc., 2021.
[38] Peter Sollich and Anders Krogh. Learning with ensembles: How overfitting can be useful. Advances in neural information processing systems, 8, 1995.
[39] Bruno Loureiro, Cedric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Mare Mezard, and Lenka Zdeborova. Learning curves of generic features maps for realistic datasets with teacher-student model. In A. Beygelzimer, Y. Dauphin, P. Liang, and Wortman Vaughan, editors, Advances in Neural Information Processing Systems, 2021.
[40] Devon Janke and David V. Anderson. Analyzing the effects of noise and variation the accuracy of analog neural networks. In 2020 IEEE 63rd International Midwest Symposium on Circuits and Systems (MWSCAS). IEEE, August 2020.
[41] A. Aldo Faisal, Lue P. J. Selen, and Daniel M. Wolpert. Noise in the nervous system. Nature Reviews Neuroscience, (4):292-333, April 2008.
[42] Jack Sherman and Winifred J. Morrison. Adjustment of an inverse matrix corresponding to a change in one element of given matrix. The Annals of Mathematical Statistics, 21(1):124-227, March 1950.
[43] Preetum Nakkiran, Prayaag Venkat, Sham Kakade, and Tengyu Ma. Optimal regularization can mitigate double descent. arXiv preprint arXiv:2003.01897, 2020.
[44] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. Imagenet: A large- scale hierarchical image database. In 2009 IEEE conference on computer vision and pattern recognition, pages 248-255. Ieee, 2009.
[45] J Howard. Imagenette. Accessed: 25-10-2025,
[46] Saining Xie, Ross Girshick, Piotr Dollár, Zhuowen Tu, and Kaiming He. Aggregated residual transformations for deep neural networks, 2017.
[47] Diego Doimo, Aldo Glielmo, Sebastian Goldt, and Alessandro Laio. Redundant representations help generalization in wide neural networks. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances Neural Information Processing Systems, 2022.
[48] Alexander Maloney, Daniel A. Roberts, and James Sully. A solvable model of neural scaling laws, 2022.
[49] Blake Bordelon and Cengiz Pehlevan. Population codes enable learning from few examples by shaping inductive bias. eLife, 11, December 2022.
13

---

[50] В OLSHAUSEN and D FIELD. Sparse coding of sensory inputs. Current Opinion in Neurobi- ology, 14(4):481-487, August 2004.
[51] Luc Devroye, Non-Uniform Random Variate Generation. Springer New York, 1986.
[52]
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury. Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems 32, pages 8024-8035. Curran Associates, Inc., 2019.
[53] John Silvester. Determinants of block matrices. The Mathematical Gazette, 84(501)440-4677, November 2000.
[54] Wolfram Research, Inc. Mathematica, Version 13.3. Champaign, II., 2023.
14

---

A Table of Parameters from Proposition 2 and Figures 2,3,4
Data Scale
Correlation Strength
C
Data-Taskk Alignment
P
Readout Noise Scale
7
Feature noise scale
w
Label Noise Scale
C
Ensemble Size
k
Subsampling Fractions
Vrr
r 1
K
Subsampling Overlap Fractions
Vrr'
r 7 r'
Sample Size
P + 00
Data Dimensionality
M + 00 P
Sample Complexity
a M
(finite)
Effective Noise-Tos To-Signal Ratios
H 11-c)
W
w
Z
x1-c)²
x[1-c]
B Code Availability and Compute
All Code used in this paper has been made available online (see Sarniin-Curre--oo---terrgggeeeee
Ensembles.ghn git). This includes code used to perform numerical experiments, calculate theoretical learning curves, and produce plots as well as the custom Mathematica libraries used to simplify the generalization error in the special case of equicorrelated data. The compute time required to do all of the calculations in this paper is approximately 3 GPU days.
C Homogeneous and Heterogeneous Subsampling
In this section, we describe the homoegeneous and heterogeneous subsampling strategies, as used in this work, in detail. This sampling method was used in calculating the theoretical loss curves for heterogeneous subsampling experiments seen in main text Fig. B and in the heterogeneous subsampling experiments applied to the CesNext-featuree-eee
image classification task in Figs. 3
S2
In homogeneous ensembling, the subsampling fractions Vyr N./M are chosen Vrr 1/k: for all r 1,...,k k. In heterogeneous ensembling, the subsampling fractions (V11,..., U2k are generated according the following statistical process:
Each fraction urr is drawn independently as Vrr Tk.o. where Kk represents a Gamma distribution with mean 1 z and variance a²
2. The fractions re-scaled in order to sum I: Ver + Ver/ly) +... + 1*)
Equivalently, the subsampling fractions are drawn from Dirichlet distribution parameterized by the ensemble size k and a chosen variance a as (V11..... Pisk ~ Dir((ak) ( 2.... (ak) -2) [I]]
In main text Fig. 3r, combine this sampling strategy with theoretical learning curves for the equicorrelated data model in a quasi-numerical experiment. At each trial the experiment, we draw particular realization of the subsampling fractions {VI1,.... Vkkk ) then the analytical expression (eq. (6) to calculate the resulting learning curve. Dotted lines show the loss curves for single trials, corresponding to three particular realizations of the subsampling fractions. The solid lines show the average over 100 trials.
Note that have defined our convention for the parameterization of the I' distribution in which the inverse of the mean and the standard deviation are specified. terms of the standard "shape" and "scale" parameters, we have:
I kko (shape = (ko) scale =ko²)
(27)
15

---

D Numerical Linear Regression with Synthetic Datasets
The code used to perform numerical experiments and generate plots has been made publicly available (see section B).
Numerical experiments were performed using the Py PyTorch library [22]
In numerical regression experiments, synthetic datasets with label noise are constructed as described in section 2.1. drawing data randomly from multivariate Gaussian distributions and adding label noise (see "DatasetMaker.p"" in available code). Representing the training set terms of a data matrix T E RMMx Mx P in which column l consist of the training point 2 and the labels with a column vector y such that y y/u the learned weights are calculated as:
v = v
(
Tu
v + AIp)
(28)
In the ridgeless case, pseudoinverse is used:
i = V
(29)
E Ensembled Linear Classification of Imagenet Images
In this section, we provide the details of numerical experiments which demonstrate that qualitative insights gained from our analysis of the linear regression task with Gaussian data carries over to practical machine learning task. In particular, we apply ensembles of linear classifiers to datasets constructed using a pre-trained ResNext 46 specific type of Convolutional Neural Network (CNN).
E.1 Dataset Construction
To construct the dataset, we start with set of n images {
Jn-1 from subset of C' = 10 classes of the imagenet dataset. For each image æt obtain corresponding feature vector zh" RM as the last-hidden-layer activation of the ResNext [46], which has been pre-trained on the imagenet classification task [44] The architecture we produces M 2048 features per image. These features will serve as the data input to the downstream linear classifier. The corresponding labels
8 E RC are one-hot vectors.
construct two datasets using this method, using images from the "Imagenette" and "Imagewoof" datasets [45]. For the "Imagenette" task, construct training set of size ner 9469 and test set of size Teeest = 3925 containing features corresponding to images from 10 unrelated classes (tench, English springer, cassette player, chain saw, church, French horn, garbage truck, gas pump, golf ball, parachute). For the "Imagewoof" task, we a construct a training set of size 74r = 9025 and test of size nestt 3929 containing features corresponding to images of 10 different dog breeds (Australian terrier, Border terrier, Samoyed, Beagle, Shih-Tzu, English foxhound, Rhodesian ridgeback, Dingo, Golden retriever, Old English sheepdog). The imagewoof classification task is naturally more difficult. The statistics of the resulting datasets are described in Fig. [ST] where we plot the data-data covariance matrix, feature covariance matrix, and the eigenvalue spectrum for both the "Imagenette" and "Imagewoof" tasks.
E.2 Model Training
At training time dataset of P examples is constructed: D {zb"" 1 We represent the training set with "design matrix" + E TR. PxM and a label matrix Y E (0, PxC The loss function for each ensemble member is generalized to multi-class regression. With ridge regularization, the objective for each ensemble member r becomes:
W, arg min
Σ
TW, + - --1 +
A. 1 W
W.C
NN
16

---

Imagenette
Imagewoof
a
Data Covariance Matrix
Data Covariance Matrix
0.020
0.015
0.010
0.010
0.005
0.005
0.000
0.000
b
Eigenvalue Spectrum
Eigenvalue Spectrum
100.0
100.0
Eigenvalue (A)
Eigenvalue (A)
10.0
1.0
0.1
Power
Eigenvalues
0.1
Power-law fit
Eigenvalues
0.0
0.0
218. 14 x
-.880
87 74 X 0 0.85
100
Rank (n)
Rank (n)
C
P 1000 M 2048
A
P 1000 M 2048
A
10a
0,0e+00
10a
0,0e+00
1.0e:03
1.0e-03
Test Error
1.0e-02
Test Error
1.0e-02
1.0e-01
1.0e-01
1.0e+00
1.0e+00
10-2
1.0e+01
10-1
1.0e+01
1.0e+02
1.0e+02
1.0e+03
1.0e+03
Trin
1.0e+04
1.0e+04
10
10
10
10
10
10
10
10a
V =
N/M
V = N/M
Figure S1: In numerical experiments, train linear classifiers to predict labels of imagenet images based on their hast-hideen-laye representations a pre-trained RexNext deep learning architecture [46] Here, show the structure of the datasets constructed using the ResNext feature map for the Imagenette task (left), which consists of categorizing images from 10 unrelated categories, and the Imagewoof task (right), which consists of categorizing images from 10 different dog breeds. (a)
c + for data matrix E IR PxM where P is the total size of the dataset. Dataset is sorted by label and tick marks show the boundaries between classes. (b) The covariance eigenspectrum of the ResNext features is well described by power law decay. (c) Generalization error of I Linear classification with single linear predictor with access to fraction v N/M of the ResNext features averaged over 100 trials (see discussion in section E.4.
Gram matrix of the centered ResNext features defined as p. { ক
17

---

where || - || denotes the Frobenius norm and where the readout noise vector &N ~ (O ( is C dimensional readout noise which corrupts the model's prediction. The weights W, can be determined in closed form as follows:
A. A AA + +X21) ( -
Wr VN, N
where we have defined [Erluc = [e When = we instead use a pseudoinverse rule. In all experiments presented here, subsampling of the features, so that the rows and columns of A, consist of one-hot vectors.
measurement matrices A, which implement ordinary
E.3
Model Prediction
Once trained, the learned weights may be used predict the label of new example 1/ as follows. For r = II...,A A we calculate the prediction of each ensemble member by first assigning each class "score". The scores of predictor stored in vector fr C
f+(+)= VN,
AW +
Where &r IC)) drawn randomly at model evaluation. Each ensemble member's "vote" corresponds to the class with the largest score. The prediction of the ensemble is then calculated as majority vote of the ensemble members. Generalization error then calculated as the probability of misclassifying example from the test set.
E.4
Linear Classification Experiments
We apply the described majority-vote linear classifier ensembles the ResNext-lmagenette and ResNext-Imagewof tasks in three different experiments.
E.4.1 Reduncancy of ResNext Features
In the first experiment, we investigate the performance of a single linear classifier (A: = 1) as fraction of features which it has access varies. We set P 1000 and vary logarithmic scale from 10 10-3 to We also vary the regularization strength over 0 and a logarithmic scale from 10 3 to 10a We average over 100 trials. At each trial the particular subset of P = 1000 training examples and the particular subset of v M features is randomly re-sampled. We find that ResNext features
50 values
highly redundant classification accuracy is very robust to subsampling of the features. For example, in the Imagewoof classification task with the best regularization tested (A 0.1), test error increases from about 1% to about 2% as the subsampling fraction decreases from = to / = 0.1 (meaning 90% of the features are ignored) (see Fig. 3d). Similarly, for the imagenette task, test error increases from about 0.1% about 0.2% the subsampling fraction decreases from // = 1 to = 0.1 (see Fig. Slb... Code used to run these experiments may be found in the folder "DeepNettSubsampp" in the GitHub repository.
E.4.2 Heterogeneous Ensembling Mitigates Double-Descent
In the second experiment, we compare learning curves for homogeneous ensembling and hetero- geneous ensembling applied to the ResNext- Imageewof classification task. each trial, we train ensembles of & = {1,4 4, 8, 16, 32} linear predictors whose subsampling fractions {Urr. ... V2k are assigned either with Homogeneous ensembling or with Heterogeneous Ensembling with a = 0.75/K
After subsampling fraction are assigned, the training set is randomly shuffled. We iterate over 50 sample sizes P logarithmically distributed from 400 4000, and then add the values M/k: for each k to the list of P values. We repeat for 100 trials for both X,- (pscudoinverse rule) and À 0.1.
18

---

A = 0, Homogeneous
A 0, Heterogeneous
10
10"
Test Error
Test Error
10-1
10-1
16
32
10
10 1
100
10-2
10
100
P/M
a P/M
A = 0.1, Homogeneous
A = 0.1, Heterogeneous
100
100
Test Error
Test Error
10-1
10-1
16
32
10-2
10 -
10
10
10
100
P/M
P/M
Figure $2: Results of the experiment described in section E.4.2
Learning curves for ResNext-
Imagenewoof classification task linear classifier ensembles. In Homogeneous ensembling (left) all Vr+ 1/k with & indicated by the legend. heterogeneous ensembling Vrr are drawn from Dirichlet distribution as described in section with a = 3/ 3/44). We use regularization À = 0 (top) 0.1. Lines represent an average over 100 trials, shaded regions show
and À 0.1 (bottom). standard deviation. We set Vpr" for r rr.
which found to mitigate double-descent by the parameter sweep in Fig. STc. In main text Fig. Bl we plot the mean learning curves over 100 trials. In Fig. [S2we we show standard deviation over the 100 trials shaded error bars. When À O, heterogeneous ensembling mitigates double-descent, leading to monotonically decreasing learning curve for sufficiently high k. When 0.1, homogeneous and heterogeneous ensembles of size & perform similarly. Code used run these experiments may be found in the folder "DeepNet_Hone Voreeee"""""" in the GitHub repository.
E.4.3 Readout Noise Encourages Ensembling
In the third experiment, we test the effect of a readout noise which is independent across the members of the ensemble on generalization error. We do this by sweeping over the readout noise scale r as defined in sections E.2 E.3 For X 0, 0.1 and = 0, ......... we compute the learning curves 10, 15, 20, 25, 30. 100. 200 and all Vrr 1/k, averaged and
of linear predictors with 1, 2, over 50 trials. These learning curves are shown in Fig. [S3 for both the ResNext-Imagenedt ResNext-Imagewoof task. Figs. [S4p and S55,, we plot the value ke which minimizes error as pase diagrams in the parameter space of a and 7/. analogous to the phase diagrams in Fig. 4p. We see that the qualitative shape of these phase diagrams is similar to equicorrelated model. The differences may be attributed to the nonlinear nature of the classification task. Furthermore, find
19

---

that, in general the optimal h" tends to be higher with the Imagenette dataset than with the Imagewoof dataset, agreement with our finding the equicorrelated regression model that as p increases (making the classification task easier), optimal ensemble size tends to increase (Fig. 4 and [66] SS7
S8p. In Fig S3 we see that there are often a number of values for which test error or near to its lowest. To quantify this, we also plot diagrams of the minimum and maximum values of l: that are within an small margin of the minimum measured error. For the ResNext- Maaaeett-m task, we use 0.001 and for DesNext-Imagewer bring error near-to-minimun
0.01. We see that there is a wide array of k values which practice (Figs S4p.c. c. S5,,c). We also plot the minimum achieved error, and the difference between minimum errors at À 0.1, 0 (Figs. SS4. (Sdd). Code used to run these experiments may be found in the folder "DeepNet PD" the GitHub repository.
20

---

K Values
1
5
15
25
100
2
10
20
30
200
Imagenette
Imagewoof
A 0
0.1
== 0
0.1
10a
100
0.0
1.00
10-1
10
10-2
10-1
10
109
10-1
100
0.2
0.2
10-1
10-1
10
10-2
U
10-1
10"
100
= 0.4
1004 0.44
10-1
10-1
10-1
10 2
10-2
10a
100
00.6
0.6
10-1
10-1
10-1
10-1
10"
0.8
0.8
0.75
0.6
0.50
0.4
0.25
10-1
10
0.2
0.8
0.8
1.0
0.75
0.6
0.8
0.6
0.50
0.4
0.6
l1
0.4
0.4
0.25
0.2
0
0
= P/M
= P/M
P/M
P/M
Figure $3: Learning curves for ensembles linear classifiers with homogeneous subsampling for
= 0, 0.1 and readout noise values indicated in the figure. Results are averaged over 50 trials. Experiments described section
21

---

K imagenette2
K' = imagenette2
1.0
1.0
0.75
0.75
200
0.5
0.5
100
0.25
0.25
0.0
30
b
1.14
1.5
0.05
0.41
25
Kmin Ả imagenette2
Kmin Ả 0.1 imagenette2
1.0
1.0
20
15
10
0.25
0.25
0.0
0.0
5
0.05
0.41
0.77
1.14
1.5
0.05
0.41
0.77
1.14
4
C
Kmax imagenette2
Kmax imagenette2
3
1.0
1.0
0.75
0.75
0.5
0.5
1
0.05
0.41
0.77
1.14
1.5
0.05
0.41
0.77
1.14
a P/M
a P/M
d
Optimal E 0 imagenette2
Optimal E 0.1 imagenette2
1.0
1.0
10
0.5
0.5
10-2
0.05
0.41
0.77
1.14
15
0.05
0.41
0.77
1.14
a = P/M
= P/M
E(A 0.1) E(A 0) imagenette2
10-1
1.0
10-2
0.5
0
-10
-10-1
0.05
0.41
0.77
1.14
1.5
a P/M
Figure S4: Diagrams described in section E.4.3 for the ResNext-Imagann experiment. Using learning curves in Fig. SS (and for additional values of not shown there). we plot (a) Optimal k in the parameter space of and (b) Minimum value of & for which error is within a tolerance <= .001 of its value for k., (b) Maximum value of for which error within a tolerance c = .001 of its value for k*, (d) the value of the minimum error E(k"), and (e) the difference between this optimal error for = 0.1 and À =
22

---

imagewoof2
K' À 0.1 imagewoof2
200
1.0
1.0
0.75
0.75
100
0.5
0.5
30
0.25
0.25
25
0.0
0.41
0.77
1.14
1.5
0.05
0.41
0.77
20
b
Kmin à imagewoof2
Kmin 0.1 imagewoof2
15
1.0
1.0
0.75
0.75
10
0.5
0.5
5
0.25
0.25
4
0.05
0.41
0.77
1.14
1.5
0.05
0.41
0,77
1.14
1.5
3
Kmax À imagewoof2
Kmax à 0.1 imagewoof2
C
2
1.0
1.0
0.75
0.75
1
0.5
0.5
0.05
0.41
0.77
1.14
1.5
0.05
0.41
0.77
1.14
a P/M
P/M
d
Optimal E A = 0 imagewoof2
Optimal E = 0.1 imagewoof2
1.0
1.0
10-1
0.5
0.5
0.0
0.0
10-2
0.05
0.41
0.77
1.14
1.5
0.05
0.41
0.77
111
P/M
P/M
E(A 0.1)
0) imagewoof2
10 10-1
1.0
10-2
0.5
0
-
-10-1
0.05
0.41
07
1.14
1.5
a P/M
Figure S5: Diagrams described in section E.4.3 for the ResNext-Imagewee experiment. Using learning curves in Fig. S3 (and for additional values of not shown there). we plot (a) Optimal k* in the parameter space of and 17. (b) Minimum value of & for which error is within tolerance e == .001 of its value for k*, (b) Maximum value of for which error is within a tolerance e = .001 of its value for k. (d) the value of the minimum error E(k"), ) and (e) the difference between this optimal error for = 0.1 and À =
23

---

F Generalization error of ensembled linear regression from the replica trick
Here we use the Replica Trick from statistical physics to derive analytical expressions for Eypt. We
treat the cases where = r' and + r' separately. Following statistical mechanics approach, calculate the average generalization error over a Gibbs measure with inverse temperature 3
Z
/
II
duo exp
3
MB
Σ
Jrr Ery (wr, (wr;
(30)
2
r.y"
F Σ
w. A, +
yp
+ X A
(31)
Nr
In the limit where 3 + the gibbs measure will concentrate around the weight vector w, which minimizes the regularized loss function. The replica trick relies on the following identity:
(o(()( lim log ((Z))D) n-00 m
(32)
where ( represents an average over all quenched disorder the system. In this case, quenched disorder the disorder which is fixed prior to and throughout training of the weights consists of the selected training examples along with their feature noise and label noise: = { ot, e - The calculation proceeds by first computing the average of the replicated partition function assuming n is a positive integer. Then, in a non-rigorous but standard step, we analytically extend the result to + 0.
F.1 Diagonal Terms
We start by calculating E,, for some fixed choice of r. This derivation partially follows section D.3 from (36), with the addition of readout noise and label noise. Noting that the diagonal terms of the generalization error Eyr only depend on the learned weights w,, and the loss function separates over the readouts, may consider the Gibbs measure over only these weights:
2 = f dw, exp
(-22 21
3 Et
JM3
(-x
Err (w.))
(33)
2
(
II dw ( (ou
8M
Ar (1b, + +00) TT - VM(e" -C4 12
exp
(34)
2X
M
V
n,a
JM3
Σ Eyr (w0
2
Next we must perform the averages over quenched disorder. fa=1.
We first integrate
(ebm, ,. EH, d Noting that the scalars
L A, (bh, a*)) T Tb VM( (t)
1
h
VM
Gaussian random variables (when conditioned on A,) with mean zero and covariance:
24

---

(g - Qab
(35)
1
I(
( A, w²T 1=,
1 (_ A : w"
Qab
M
(vrr
(36)
1 u
T
A, w. +M( + tz))]
V/r
To perform this integral we re- we-write in terms of (H)) where
hrl
hr
H
IR"
(37)
hmm
1 /II II B on exp
(-44
B Σ
H
H
3
Σ
w
JM3
Σ
Err (w0
21
2
2
P
a
a
(38)
Integrating over the H. we get:
(2")D /II
II d exp
log det
B In +
((
9
Σ
JMB
Σ
Ery (w,)
2
2
(39)
Next we integrate over Q" and add constraints. We use the following identity:
1 1= II / dQH
Q
1 M
AA, A, A,A wh +M(s2 v/:i]
I
( A. w E.
(JAA) A, w
w0)
Vvr
(40)
+ 1 w.
Urr
Using the Fourier representation of the delta function, get:
II
/
dQ exp
d
M
Q
1
[ Vrr
we A, w E,
1
A, u
A,
w
4Ti/M
M
Vir
1 Urr
:.
A, A. 1 +M + +n,)]
(41)
Inserting this identity into the replicated partition function and substituting (w) Qi (2 we
find:
25

---

( D
II daidd exp
(
(-2
P
log logder del
B (n+
(
In +
0m))
+
1
Σ
JMB
Σ (2)
2
2
ab
/ II d exp
8
Σ
E 2at ((
(
A.
w²T)
Es
A, A.u w²
Ver
Vur
A, E A : + M(5² + r//) )
+ Vrr
(42)
In order perform the Gaussian integral over the {wf}, we unfold over the replica index a. We first define the following:
w.
(43)
= 3In & IN. + Q"" 8 A, Eo)A7
(44)
Urr
V"
Q" (1n &
A. A.
(45)
Vpr
We then have for the integral over
f fdu, dww, exp
(-,,; Twi, + VrT VTTwi
w.
(46)
exp ((T 2
log det
(47)
We can finally write the replicated partition function as:
) D x
/ II di exp (-5 log det
P
$ In
JM9
I Qm + Σ MQTQ.
Σ
(Qn -22)
2
exp ( (TT
log det((T))
Σ (M + +,² + +w" + "2xw""
2
(48)
We now make the following replica-symmetrie ansatz:
QT = 1- qoan + 90
(49)
Q 3Gdub 38
(50)
which is well-motivated because the loss function is convex. We will verify that the chosen scalings are self-consistent in the zero-temperature limit where oc. We may then rewrite the partition function follows:
26

---

n.M (2")D - exp (-22
1g. 9o. q, dol
(51)
Where the effective action is written:
g 14, 40: 0.900
log(1 + 2)+
Ago
(qq + 3q9o + $9o9) + 3J [(8 9 + 90) -ca]
U4 M 1 u E,A, G- A,E, +
B
1
[og del + 300 ( + 3ā c² + 1/² M == E,w²
M
(52)
Where G = IN, + + and Σ
1 A,, + E,)AT Myr
To determine the values of the order parameters set the derivatives of [g. qo, q, 40). evaluated at zero source (J = 0), to zero:
0g 0 aB
34
a
(53)
Oqu
Xtg 9
x+ q
0g 0 = -Bq +
B
1 M
tr - Σ
tr G
(54)
Og0
M
0g
0
a
By) - q - -340
aq0 3
+ 90
aqu
(55)
0q
x+
(A
+
(X+ + Bq)
0g
0 = 9 - 390 + 362 + Br7 +
1
tr G
340
tr (G-1
27
00
=
M
M
(56)
8 w*T T
2
1 E A G '5G' A,E,
Σg
A G A,E, +
w
M
Vrr
Vrr
- 9u=
1
2 +
2 + M
w" T Σ,
- A G G-
+
G-
w"
1-
Vrr
Ver
(57)
where y= M tr [G-15)2 and = X + q. We retroactively confirm that the chosen 3-scalings are correct by noticing that the saddle-point equations permit a solution where all order parameters remain (1) as B + c.. The error is given as:
1
E,r = lim
8 B J. glg. qo. 4.9 = lim 1 + 90
22
52 =
90
-
(58)
8-500
Where q, 9o: g. go are the solutions to the saddle-point equations 53 54 55 for 40. we obtain:
57
Substituting eq.
57
Epr
1
1
2 Σ, G.A,
G-A,,GGA,,,,,, + P.EAA
SEG G-GG A,E, w0 +
2S2 + nr
1
M
Ver
Ver
1 7
(59)
1 +
1 PE,A, G A,E, G w² +
252 + 7
Σ,
1
M
Vrr
Ver
1
(60)
27

---

F.2 Off-Diagonal Terms
We now calculate Err for # r'. We now must consider the joint Gibbs Measure over wo, and w,
JMB == / dww exp (-2212 2x (F, + F.J) 2
Err (wr, uPrr)
(61)
(62)
D / II drof. dwo E (Wanom),en)
exp
.3M
1
(hp4)2 2 +
8 Σ 2
|2 +
1
JM3
Σ
e (w w
2A
M
2
(63)
Where the hy are defined as before. Next we must perform the averages over quenched disorder. To do so, note that the he Gaussian random variables with covariance structure:
ho = Qab
(64)
1
(
1 w. A, - w=') Eg
Viser ArTwe w³)
M
Ver
(65)
A, A +M
V UypyU
To perform this integral write in terms of (H {H where
H H
H
E R
(66)
/ II dw dw E (cn)
D
(67)
JMB
exp
Σ
H H
Σ
| + /
Err (w, w
2X
2
2
Integrating over H, we get:
(D")D /II dwo
(68)
P log det
8 Σ
JMB
exp
Izn +
/20 2 +
+/23111
Σ
Evr (wo (w, wg))
2
2
Where we have defined the matrix Q so that:
n Q
Qer² Qr'r'
Q
(69)
Next integrate over Q and add constraints. We use the following identity:
28

---

1=
II ab
/ dQnn's (0th
1
" A, w²²) E.
(Vrry ATr - w0)
1
M
(70)
W VeeWe''
A,E
+MG
Using the Fourier representation of the delta function, we get:
(M 0: (0:5 -1
I( Ver
II| ab
4Ti/M di% a'% exp
w A. w E,
Arree
w²
M
Vipp
A, A, +M(?1) :
Mrr
(71)
Inserting this identity and the corresponding statements for Q% and Q into the replicated partition function and substituting E,y (w") Qan 6 we find:
/
(z")
II abriTz
exp
P log det
Izn
+ 2
Q
Σ MQ
JM3
Σ
(Q
-x2)
abrivz
dw exp
/10012 + 120
1 Σ
2ir [(
w. Ar, - w
E.
(1,110 A
w
2
/Vr
/vrz
abrira
Ta
AA,
A,A
wh M122)
Vrr Vrz
12
(72)
Where sums over 71 and 72 over {r,r'}.
In order to perform the Gaussian integral over the {w we unfold in steps. We first define the following:
w
W
=
(73)
(Q"²") lab
Qun'
(74)
1
Σn
A, + A.
(75)
VirePe'p'
T = 35 In @ IN. + QT' 8 EST
(76)
Unfolding the replica indices, we then get:
29

---

(
D x
II
dQab "Q""
abrir2
exp
P log det
I2n
Q
+. E
12 Σ Que -22)
JM3
2
abrive
Σ
Quer (w"
"'.w" + +M
exp
abriTa
J dan, exp
Σ
TTm +
T
Σ
@mrs @ IN., (1n & /
Ar, Ez
Wri
2
1172
(77)
Note that the dimensionality of Trirs varies for different choices of rı and r2. Next, we unfold over the two readouts:
w
(78)
W
Ter Tay
T
TT'
(79)
V (ÔTT +Q IN.
1, -, 200)
(80)
(iQ'r' + Q'')) @ IN.
1, 0 マリマン
A, Eg w
The integral over then becomes:
/ dw exp (-
(-2w w Tw vTw) x exp 1 V log det T
(81)
We are now ready to make replica-symmetrie ansatz. The order parameter that we wish to constrain is Q Overlaps between the weights from different replicas of the system as well as different readouts. The scale of the overlap between two measurements depends on their overlap with each other and with the principal components of the data distribution. An ansatz which is replica-symmetric but makes no assumptions about the overlaps between different measurements is as follows:
Tri 3- 1 + Q
(82)
Q S Sab + Q
(83)
Next step is to plug the RS ansatz into the free energy and simplify. To make calculations more transparent, we re-label the paramters in the RS ansatz as follows:
Q" 8 4 + TQ1T
(84)
Q'' 3 RI1T
(8.5)
QT 3 T
(86)
Q sqI + "Q11
(87)
Q'r AFI + o Rm
(88)
QTT seI + T
(89)
30

---

In order to simplify log det (A + BQ), we note that this is symmetric 2-by-2-block matrix, where each block commutes with all other blocks. We may then use [53] result to simplify.
log det (AInn + +BQ n log ((A + (1+ -) -c) ( (A+ q) (A (2+ -2C
+ On²2)
(A+
(90)
Σ
2a: 2ab = [(100 3
3qQ + )nQQ + (rF+ 3i Br + (cē BCC + (xC++BCC))))
- O(n²)
(91)
Σ (oad - 2) =
1 3 + - -
+ +0(2)
(92)
Q Bn l0 + f 26 On2))
(93)
abrira
G Gr
G Gp
Grr Gr'r
G G
QE., Ö
]]
log det))t())= = n log(3) + log det
+ t
+0
(Ser
RE
(94)
where G.r IN. + E G.y IN., +
G =
=
(95)
(,+ (r+ Eg
G.r Gp
G G
(++
VVT-IV= = =nww"T T
w + O
(,+
(96)
Collecting these terms, we may write the replicated partition function as follows:
= exp (- nM
(-121 8
D=x)p=exp
(97)
Where the free energy is written:
0 a. r, R, R,
(98)
a log ((A + () r) -2) + (A ( + (A +) 20C
(99)
(A (A++)) -
- [ +3Q + 3qQ) + (f B Br = (30+
38C + BcC]]
(100)
+ J(c+ 30 - -3(3)
(101)
+ + f 20]
M w0 Enw" + +2²
(102)
1 M
( Ju,e (f+6)
(
Bwn
G-1
Waaee (f (+ nn
(103)
1
G- Q (
1
+ M
log(3) + let G + Str
(104)
C E'p'r RE
31

---

Grr G
Gyy
where we have defined G
GP
The saddle-point equations for the replica-diagonal order parameters are:
0g 0 8
a(A + r)
(105)
aQ
(A + ()( + r)
-
0g
1
Σ
0
0
3q + 8
tr G-1
(106)
QQ
M
0g
a + 9)
0 = 8
Br
(107)
(A + () A + +r) -
0g 0 - 3r + B
0
Σ ))
tr G-11
(108)
A
M
0g
2ac
0 -
230 + BJ
(109)
OC
(A + g)(A
r) c2
0g 23c +
1
ETTT
G-1
(110)
aC
M
Σ E''r
Note that when 0, the saddle point equations T09 case the remaining saddle-point equations decouple over the readouts (as expected for independently trained ensemble members) giving: For readout
110 are solved by setting e = and in this
q
(111)
(A + q)
0 -4 + M
tr (G=LE
(112)
and for readout r':
(113)
(A+r) 0 rr tr (G-L ( M
(114)
These are equivalent to the gaddle-point equations for single readout given in equation 53 54 as
expected for independently trained readouts. is physically sensible that e = 0 when J = 0. because source, there is term in the replicated system energy function which would distinguish the overlap between two readouts from the same replica from the overlap between two readouts in separate replicas (we expect that the total overlap between readouts is non-zero, as we may still have C>> 0).
The saddle-point equations obtained by setting the derivatives similarly decouple to recover two copies of the diagonal case expressions here as they are not necessary determine the off-diagonal error term Erp:
8g
a will
8g
86
8F
57|55]
We will not re-write the
0g The remaining saddle-point equations are obtained by setting
0g 02
de
0g
20.3C
aC (++ ((++r)
2,3
Ō=
(115)
dc
(A+
cd=J=
32

---

0g
1 "n'TEww E,w² + +32
0 - 23C +
M "0u""
0c
r===J=
23
1
qA, G A. 1 A G.,, Ap
w'
E,w
M
Vrr
V/'r'
(116)
2
1
w E, A, G Ep. E.
M VireWer''
2BB
tr rr Ep
G p'y!
E''r
M
Solving equations 115 and 116 for C, we obtain:
C
2
M 20 TT,""
1-?
w" M.
i G Ar + -A And
E,w²
(117)
M(1 - 7)
Vrr
1
qf
Σ, A G G p'y! A, E,w0
M(1 -
where
7
G G
(118)
(A + +)+)+++r)
We obtain the generalization error as:
E,r- lim 1 0 8-+oe 3 2 R,
1 8+x
lim
(c+ = C (119)
We may then simplify the expression for the Err error as follows:
E,y'
7
2 +
w
Σ,w"
1 7
?
M
1
1
1 G.A.G
w² Eg
QA, G A,
Esw"
(120)
M(1 - 7)
Vrr
Vp'g'
+ M(1 7)
qp
TT E.A G- Epp G A, E,w"
VerMer'
Re-labeling the order parameters: I + Gr. Gre, 7 Yrr and G Gr, we obtain the result given in the main text.
G Derivation of Proposition 1 from [31
In the case where the data and noise covariance matrices E, and Ep have bounded spectra, our main result may be derived using Theorem 4. from Loureiro et. al. (31), with few additional arguments to incorporate readout noise which varies across ensemble members, and to allow for the presence of label noise in the training set but not test time. Rather than reproducing their very lengthy statements here, direct the reader theorem 4.1 and corollary 4.2 in [31].
33

---

G.1 Altered Expectations at Evaluation
In 3T, labels are generated at both training and test time through the same statistical process y Pg (ylz) where / (=/0)
va However, their results may be easily extended to the case where
data are generated through different statistical process for the training set and at evaluation. This will allow the application of their results to the case where label noise is present during training but not at evaluation as studied in this work. We therefore introduce separate distributions of the labels at training and evaluation. During training, we still have ~ Py (ulp). At evaluation, we put y Py (yle). This leads the updated formula:
)}
(WW UU))
(yy
(121)
Vp
when the expectation is over data-label pairs (y, x) at evaluation.
G.2 Rigorous Proof of Proposition
We now restate the the problem setup of main theorem using notation consistent with [31] We study generalization error in an ensemble of estimators {uix). ). k I,..., K. We say th E R for all & ...., K. The weights are trained independently such that each minimizes ridge loss
function:
2
e arg min w
y"
u
Eh
h ==
1, 1....,AA K
VN&
(122)
2
arg min
y
vp u
E
w
Vikk
Visk
So that our results will correspond to the results of [3] after re-scaling the regularizations Ag PRgAp. The training labels are drawn as y" ~ Po 0 x²) where (y]x)) = N(x (²), &k ~
NO n))) and ! represent the ground truth weights. For k = 1,...,k we have a "measurement matrix" A& E R.N. xd and we may set d = p 2 max; N₂ (so that 7 1). The feature vectors us (œ) are then drawn as
Ax + rr) (pr-Na)
(x:)
Ax σr)
where x (O0 E.) and -0, Eo). For convenience, we have defined the auxilary matrices
A
A
E RPxp
0 (p-NNxxxp
By constructing the feature vectors as p-dimensional vectors with only N₁ non-zero components, we may apply the results of (31 while preserving structural heterogeneity (we may have N₂ * Ne. for
k * k') as is present in our main result. Because (uæ): = 0 for all i >> NR. these auxiliary dimensions
will not affect model predictions or generalization error. We then have
J u wū and labels
I
generated evaluation according to
P8
Σ px) where (y(x)) = N(a, &
Σ 2 k
(123)
K
K2
The generalization error may then be decomposed as Fg
R En. Ers where Exs E (.)) ((y (v) [. We will apply eq. 121 separately to calculate Epy in the cases where * r' and r''
34

---

G.2.1 Off-Diagonal Terms
Eq. 121 cannot be used calculate E,r directly in the case where * r' due the presence of noises Sg which vary over elements of the ensemble (indexed by k) in the loss function (eg. 122). We may argue, however, that the presence of readout noise noise has no effect on the expected value of E,r when + r', then apply eq. 121 with form of the minimizer of the loss function (eq, 122
7l = 0. To see this, we examine the analytical
ع = U (URU& + +x (y Exx)
(124)
where we have defined the design matrices [Ux]ip = vp
1 |ene (att) and the vectors [y]n = yn and
[Ex|μ S We then have
EKK'
1
T
(ULUK
+ +
V 1 V (URUR + +x
"UK
vD
y y (URU + +x11) +xx1) -UA)
( El (UN.UT + + UN
vp
+
y
VP
T (UNUT + + -UU)))
Ux
vP
1 ST (URU + +xI Uk
+ El
(URUT + + Un) EK
UU + 1 -U Un
(125)
Taking the expectation value over the readout noise in the training set get:
E(E Exke = V
= -V- y
(URU + AoT) Un)
T (UN.UT + +xx1) Up
(k: # k²)
(126)
This is identical to E when & = 0 for all k. We may therefore calculate the off- off-diagonal error terms by setting & = eq. [22 which gives a problem compatible with the theorem 4 of [31]. To
calculate Eryr, we appeal to theorem 4.1 and corollary 4.2 of 31. The following objects defined in 12T are given the following definitions:
K Σ
r((11
xxk) 2 k=1
20x1
(127)
(v.0(a)) (y )xx))
(128)
E(y.H) 2 k y1l2
(129)
2d zV,,,,,0) = /
Po(y z)dr /
1
exp (- (y (u- 11)2
20
(- 2( + σ)
(130)
2mo
2 + +
ERK (vy,,) 1( pw) ww)) B(u.m) If dydpgeel)) mv mx))
(131)
(v,r)) [1V px)(U px))] +
m
(132)
K2
K
p - [m], - [m]] w + [Q] Jay +
(133)
R2
1=1
Where p. and Q are defined in 31:
(v, ~ N (OL+X)
( m
m
))
(134)
Q
35

---

What remains is to determine the values of the order parameters p. m and Q. We next define the covariance matrices which characterize the feature maps. The feature-feature covariance is:
Sy Ez (2)); (IsisN.JsisNN) 1Ax(2) En) A
(135)
( [Ene]
ji
|Ette] 1ig
(136)
where we have defined Vikk
No. Ens
Venaain A (E. Eo) A. and Ens'
ERe'
Rexp, Note that while in [31, all S2y must have strictly positive cigenvalues, their result can be easily extended to cover the case where some eigenvalues are zero by a continuity argument.
The feature-label covariance is given by:
④ = Ex
(ux((r)a 0), = =1 (A, = [A,
(137)
k
0 lk
$ $
[A&E,0], [Ax.E.0)
(138)
Recalling Q1/ now obtain explicit form for the proximal h:
(u - WV (u w) Ely, w)
h argmin
VV +
(139)
2
The proximal G should not arise in this special case.
Next, we will simplify the saddle-poin equations. Simplifying where possible, we may write:
f V (h w) (I V) (y1 w)
(140)
Ou (I + V)-
(141)
27
P EN
T
0 Σ,
(142)
Wg = m QQ1122
(143)
6o p Q m
0 Σ,0 - m TTQ11 '='1m
(144)
= -a / dy [2 wp, -n) = a VV-
dyz wa, σo)
(145)
V
(146)
We can simplify the prior equation for V to the following set of equations:
VEL
- ir
Dea URRREe Ip + VAL DrM
4: = 1,..,, K
(147)
p
VEL
tr
ΣAS Vir DAcr
k' * k
(148)
p
36

---

Equations 146 and
148 are solved by setting VR V& = 0 for all le o k so that VEN Vedaer
and VI V. The diagonal components then satisty
a 1/2
(149)
14 + V&
1
Vi Eik PRRAEI
(150)
p
Which may be solved separately for each k ,,..., K. Simplifying the remaining channel equations we have:
Q a(I V) 1) [
pp) +) 11 -1m² -m1"
(I+V-
(151)
Q = VAV [c² + my my + Qxk]
(152)
m = alI v
fi
Vi
(153)
(1 + Vv.)
Simplifying the prior equations obtain (through some tedious but straightforward algebra):
QAa Qua.ua + Ver
(154)
m. VeR
(155)
Where have defined:
G2 PREdEl + Era
(156)
JKK
tr Eisk E.
(157)
Agk 0 E. G EE G A,E
(158)
P R4 = 0 E.A G A.E.O
(159)
Solving eq's 152 for Q, we obtain:
Vw
Que
YKk
(52 + mg
mw)
Ana
(k: * k)
(160)
Yak!
Combining these results, we arrive at a formula for ExAC:
ExP - -kok"
(nRRr(2+ + V.Rs Vr Re + V Ve
Σ m
(k: + k')
K2
(161)
1 Vi Jake
where "kkk"
(162)
a
Where the order parameters Vi, = 1,..,, satisfy the fixed-point equations given by eq's
149
130
G.2.2 Diagonal Terms
To calculate Eyr, we cannot ignore the presence of readout noise in the loss function. However, as the loss function is separable over the readouts k. we may calculate Ery using the results of [3]] in
37

---

the special case where K 1, incorporating the readout noise into the label noise. Concretely, we may calculate E,, as the eπror of a single linear predictor under the same setup as the off-diagonal terms except that in the training set yt ~ (ul wp x²) where (y(()) N ç + mg) This
may recovered from the calculation for the off-diagonal terms by setting k: = k' and re-scaling 5² -5² +n giving:
Esk- 1 Tikk
(TAA(S2 + p p-p-p2VRR, + 12
1
Σn
Th
(163)
K2
where the definitions of Yykk. Rg. and Aga can be inherited from the off-diagonal case, as well as the saddle-point equations 149 150
G.2.3 Full Error
The results obtained here equivalent to the results of main theorem, up to reshuffling of additive constants n/. among the error terms, and a trivial re- re-scaling of the order parameters as follows:
VE
VR.
VE
Vv
H Equicorrelated Data Model
To gain an intuition for the joint effects of correlated data, subsampling, ensembling, feature noise, and readout noise, we simplify the formulas for the generalization error in the following special case:
Σ, = s [1-c) + clMl]]]
(164)
Eg = w I
(165)
Here is a parameter which sets the overall scale of the data and e [0, 1] tunes the correlation structure in the data and w² sets the scale of an isotropic feature noise. We consider an ensemble of readouts, each of which sees subset of the features. Due the isotropic nature of the equicorrelated data model and the pairwise decomposition of the generalization error, we expect that
the generalization error will depend on the partition of features among the readout neurons through only:
The number of features sampled by each readout: N, = Myr M, for 1,..., The number of features jointly sampled by each pair of readouts Npy = Ver M for r.r' E kk
Here, we have introduced the subsampling fractions eyr
Nr and the overlap fractions Kyy Ren M
We will average the generalization error over readout weights drawn randomly from the space perpendicular to 1,s with an added spike along the direction of 1
v1 - pTP wn + pl w ~ AO, IM)
(166)
(167)
The projection matrix may be written P- = Ia 1 The two components of the ground truth weights will yield independent contributions to the generalization error in the sense that
Tr.
( (Fpr)) = (1 - (1-p²)F
0) + Fp
(168)
Calculating E,y and E,y, is an exercise linear algebra which straightforward but tedious. To assist with the tedious algebra, we wrote Mathematica package which can handle multiplication, addition, and inversion of matrices of symbolic dimension the specific form encountered this problem. This form consists of block matrices, where the blocks may be written ad IM + bl M v. where a, are scalars and sMN ensures that there only diagonal component for square blocks (when M = N). This package is included as supplemental material to this publication.
38

---

H.1 Diagonal Terms and Saddle-Point Equations
Here, we solve for the dominant values of gr and gr and simplify the expressions for E,, in the case of equicorrelated features described above. In this isotropic setting, Err, 9r, &r will depend on the subsampling only through N, = Mr
may then write, without loss of generality A. (IN. 0) E R.N where denotes matrix of all zeros, of the appropriate dimensionality.
We start by simplifying the saddle-point equations 53 54 Expanding to Lr GG-SS.)) and keeping only leading order terms, the saddle-pon pointe-pptttee equations for ar and dr reduce
9r Vrr (s(1 c) + w²)
(169)
(s(1 c) w2) 24
a
9. À X dr
(170)
Defining = s(1 - c) + w² and solving this system of equations, we find:
4. Va²a² + 20a a)v, + (a + A)2v? - aa + (a A)1,
(171)
2v,
4 aa + 2aa a)v, (a + VY²yZ + aa (a A)V,
(172)
2aA
We have selected the solution with gr > because self-overlaps must at least as large as overlaps between different replicas. This solution to the saddle-poin equatios can be applied to each of the k readouts.
Next, we calculate E,,. Expanding 7r
a Mxx
[G-15)]
to leading order in M, find:
a²av,
7r (A + 2x)²² (agr + vr))
(173)
1
2 Gr E,A, G.'A,E,+ Uhr
1 A, G.',G "',,',,,, P1
(Err) (p=) (p=
ir P1
E,
Tor M
Urr
+
Yer
(²
Ap:
1 7rr
(174)
Σ, 2 A. A,E, + 1 - G
(En) (p= (p=1)
Trr M 1 M
EG; A,E, 1m
Vrr
Vrr
+
Per
2² + n/:
'ir
(175)
With the aid of our custom Mathematica package, we calculate the traces and contractions in these expressions and expand them to leading order in M, finding:
1
s11 ((11c)) (1- (1_ (1- c)shru (g. - c) w²) + 2v,)) 1)
2r + n/
(Err) (p= - Tor
(0 c) + w vr)
'rr
(176)
1
(1 - c) - Upr) + w²) u²)
5 + n
(p = (p))
(177)
Trr
Vrr
1 -rr
In the "ridgeless" limit where x - O, obtain:
39

---

401/r
Per
(178)
(a + Prr + a - vrr
s(1 s(1-c)pp
1 -0(1-o)(a-2u0
ter
a Vrr
D,
Vra-a
(rrla(l-c)tus
(A -)
s(1 el1-eja
sf1 s[1-c)lrr
+ Morl 4a0t. Mpr
if a > Urr
a(1 (11-))+w2
(179)
Vur
afl-n)(1 11 -
+ Ver Mon tame Myr
threth
i < Ver
(Err(p= = w
Vpr
(A-+0) 0)
of -
if > Vrr
a-
Uyr
(180)
H.2 Off-Diagonal Terms
In this section, we calculate the off-diagonal error terms E,r for r * r', again making use of our custom Mathematica package to simplify contractions of block matrices of the prescribed form. By the isotropic nature of the equicorrelated data model, E,y, can only depend on the subsampling scheme through Ver, Merr' and Mert, We can thus, without loss of generality, write:
A,
Ine Onaxm X
On Xnpp On Xn,
On.xna Inxxne
On,xi Onax1
GRN.X.M
(181)
A,
On Xne
Inp Xny On Xxp,
On,eXny Xns I In
Opprxl On, Xl
E RNN
(182)
On, X
where we have defined nt, to be the number of features shared between the readouts, nr N, 11 and ni, = N₂ n, and the count of remaining features = M Thy nl f1..
Then, to leading order in M, we find:
Wor(s(( - c) + w2)² Prr (.) + 9rr) (prr + (s(1 - c) + w²) (Vrg + (s(1 c) + +))g)r)
(183)
Averaging Fry over w ~ NO IM], we get:
(D)) (p =0) = Trr 5²
1 (1
tr (r(PLE,PP))
1 - Tirp"
1 1- Yrr
M
tr PLE,
A.A. G, A,
qr A. A.
E,P1
M -
Vrr
Vydee
9r
1
tr
IP E,A G. EE, A G.
M(1 - ) VeruUsip'
(184)
(E)r D. (0 =1)
Ver 5² 1pp"
1
(1
Tpp
M(1
Tppe)
A G.
Ar 1 qr A
1
G. A.
E,1
M(1 - 7rr))
Upp
Up'r"
Gr
1
+ M(1 -(rr) VurVure'
1 A G.'
(185)
40

---

Calculating these contractions and traces and expanding to leading order in M, we get:
(Err D. (p= s(1 cc)
(1
s(1 c)
s(1 Cr
'irr'
Vrr + (s(1 c) + 19 Vyrr' + (s(1 - c) + w2 Vr.
s(1 c)(( c) + w u²)
Urr + (8(1 - c) + w (UPr + (s(1 -c) + w )r-))n))
Prr' 5²
+ 1 Yrr
(186)
(D)) (p
= (p=
s(1 - c) ((1-c)Vmr'
v.(Vrr) + Ver'
Yer c² (187)
Tow!
Verle's
Taking A we get the ridgeless limit:
40
P/rr
(A
(188)
(a + Uyr + la vrrl) a Up'e + aa vyr
(D)) (p =0=
1
Pola) + 'rr' 7rr
c² (r
(189)
Top"
where
(c- (c-1)s(u, (2a 1)(c 1)s+)>)++=)++-) 2 (c-1))88,p()) Mel(e- 115--11---ppp ررار
if a S Vm S Vyrr"
Fo(a)
(e- (c-1)e(u (c- 1lsur
paclccll+"" -ape c- 1)
if vr S a S 1/'r" if Ver S My'rt S
(190)
((o- 1]ss m
()-) -oo (c- yv, cs +s+w2)
(e- 1ys-w
(D) (p=1) 1- -rr Yer
s111 (1 c)(rrr' - v + wrvrr' )+
Yer
(A=0))
Yer
(191)
H.3 Optimal Regularization
Here, we derive the "locally" optimal regularization which minimizes the prediction error of the ensemble members independentl. Equivalently we derive the optimal regularization for a single readout (k = 1) or the regularization which minimizes the diagonal terms E,, of the generalization еиrоo. By differentiating E,, with respect the regularization X, one can show that for the optimal regularization, we will have
1 dS, d²vrr
Sr ((1 ((1-1110)) p2)10r + p21 + - + p3)s2(1 (ar - -1] = 0
(1 - (rr) dA aal1 - (mr)
(192)
is easy to show that ds, cannot be equal to 0. Setting the term brackets to zero and solving for X, find with the aid of the computer algebra system Mathematica (54)
d\
2 (a (v + + ap² + + c)u - 2p²)) 1).
(0<c51)
(193)
-1- (1 2 / (1 -2) s²
In the limiting of isotropic data (c 0), the optimal regularization reduces to:
41

---

1 _ + 1²) + Vyr
1 (1)))) Vrr
X
(c= 0)
(194)
We note that while generalization error is discontinuous at c 0, this result can be quickly obtained from eq. by noticing that the formula for the generalization emror with 0 S 1 reduces the formula for generalization error with = O when we set c = 0 (when с = 0, the generalization error does not depend on p because there is special direction in the data).
H.4 Homogeneous Ensembling with Resource Constraints
We make further simplifications in the special case where X = X, 7 = 17, Vrr = for all r = -....k k. For simplicity. we also set Byy'' = for all r * r' so that ensemble members sample mutually exclusive subsets of the features. We will consider both the ridgeless limit À - 0 and the case of "locally optimal" regularization these special cases, generalization error has the forms given eqs. 25
X (see eq. 193p. Concretely, we show here that under
26
We start by analyzing the saddle-pe pointt equations point equations will be identical for all r= 1,...,k k. We will therefore suppress the r index. Solving this quadratic system of equations explicitly, we encounter the following radical, which we assign variable name
112 Note that this special case the saddle
r Vaa av + )2)2 4aAv²
(195)
In order to simplify this radical to begin extracting the factor of s(1 - c) that appears in eqs. 25 define reduced regularization:
26
(196)
s(1 - c)
And substitute A, H. W, and Z (recall eq. 24p for X, v), w, and Ç. giving
I = s(1 c) W,A)
(197)
where X'(a, v, W, = v (a(W 1) ((A W 1))² 4A +1)
(198)
making the same substitutions, we can then write
is(1 11(1-))
(199)
20v
where Q = (-a(W (W 1) v(A W + 1) + X
(200)
Continuing to make substitutions in this manner, arrve at:
S= [s(1 -c
(201)
where S= / + Q W)
(202)
Finally, we may express the generalization error in reasonably compact form. In the special case at hand, the generalization error is written:
1 k- 1
F
(Urr Eur
E,
Em K
k Spr:)
(203)
Substituting 201 for S in eq. 16 and making further substitutions necessary, arrive at:
Eyr(Pr 1 ) -10011--@pp,,....AAAHH..
(204)
42

---

where
Err N1 + N22
(205)
D
Ni a(H + 1)k + S(8(W + ()(a +lS WZ + 2a))
(206)
N₂ ap² (1:
(207)
D ak ( +1
(208)
We also obtain
E (Ver' 8pr)) - (.
(209)
k
where
2 20p 1)
E,
2p² +1
(210)
k
Combining these, we obtain the error of the ensemble as:
Eg =11= 1111-c)E(0 p... H.W.
(211)
where:
1 Err
b: - k
& = k
E,p'
(212)
These result are derived in the Mathematica notebook titled "MaqiCoorParameteReectiiiio included with the available code. It follows from eq. 21T that when = error can written in the form of eq. [2 It follows from equations 211 209 that at "locally optimal" regularization X", error can be written the form of eq. 23 To see this, note that the reduced regularization A' which minimizes E,r will only depend on the other arguments of E,,. Full expresions for the generalization error in the case X > 0 and À = X can be found in the mathematica notebooks "EquiCorrPhe Caseerraaaaaaaaaaa ZeroReg.nb" and "EquiCorPPhaseDic" the available code. These equations are long and difficult to interperet nor are they directly used in code- and are omitted here.
LocalReg.nn included with
H.4.1 The Intermediate to Noise-Dominated Transition
The transition between the intermediate regime where 1 < k co and the noise-dominated regime where j-' oo can be studied analytically relatively painlessly in the ridgeless limit -> 0. The strategy we employ determine this phase boundary is to examine the large-k asymptotic expansion of E& to determine whether the error approaches its asymptotic value from below or above. If Es approaches Er from below, then A- must be finite. If E&. approaches E.zo from above, then k- may be infinite however, there still the possibility of bt* < oo F& is non-monotonic in k. In practice, we check the values of Eg for 1,2.... 100 and oc.
Setting = and expanding Eg around oc, we find:
Es s(1 c)
+ + (1+W)P'
H HW +
=
+0
(213)
(1 + W)ack
kz
Setting the coefficient of -1 to zero, we find the phase boundary as:
(1 + W)
a 2(1 - H W)
(214)
This equation explains the shapes of the boundaries between the intermediate and noise-dominated regions in the phase diagrams of fig. 4 with À 0 (see black dotted lines in pancls b, d).
43

---

An analytical formula for the boundary between the intermediate and voise-dominated regimes at locally optimal regularization X = X cannot be easily obtained. To understand why, we can asses the large-k asymptotic expansion of the generalization emor at locally optimal regularization:
H k
Fx x²) ((1-)))))))))))))) p² W
+
O
(215)
12
This shows that at large kt when H > O, error always approaches its asymptotic value from above (recall that when 0 we always have ks so that there is no phase boundary unless H > 0). Thus, determining the phase boundary requires determining the value of k which minimizes Eg. which is not analytically tractable.
H.5 Infinite Data Limit
In this section we consider the behavior of generalization enror in the equicorrelated data model as + oo while keeping the À O(1) For simplicity, we assume Vrr' for r * r', isotropic features (c = 0), no feature noise (w = 0) and uniform readout noise 1/- 1 as in main text Fig. 3. This limit corresponds data-rich learning, where the number of training examples is large relative to the number of model parameters. In this case, the saddle point equations reduce to:
Gr +
(216)
Vrr
gr
(217)
a
In this limit, we find that Yrr - Using this, we can simplify the generalization error as follows:
2 2- 2- (--
n
E
Σ E =
(218)
Interestingly, find that the readout error in this case depends on the subsampling fractions vrr only through their mean. Therefore, with infinite data, there will be distinction between homogeneous and heterogeneous subsampling.
I Theoretical Learning Curves and Optimal Subsampling Phase Diagrams
Here, we provide additional learning curves and phase diagrams of k" such as those in Fig, bb..cc.d exploring more parameter values for the task-model alignment and the Reduced noises H, W, and Z. Generalization errors are calculated for homogeneous ensemble of k linear regressors, as described in sections 5 H.4 We also show diagrams of generalization error Eg..
44

---

a
p 0, 0
p 0.3, λ 0
p = 0.6, A = 0
2.0
2.0
2.0
90
1.0
1.0
1.0
70
60
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0 m
0.0
10
2.0
3.0 m
50
α
α
α
40
0
O, À A
p
0.3, À A
P
0.6, À
A
2.0
2.0
2.0
20
I
1.0
I 1.0
I 1.0
0.0
0.0
0.0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
10
2.0
3.0
α
b
p
O,
0
p 0.3, À 0
p 0.6, À =
2.0
2.0
2.0
1.0
1.0
1.0
1.0
0.8
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0 8
0.0
10
2.0
3.0 m
0.6
α
α
α
p 0, À
A
p
0.3, À A
p
0.6, À =
0.4
I
1.0
I 1.0
1.0
0.0
0.0
0.0
0.0
10
3.0
1.0
3.0
10
3.0
α
α
Figure S6: (a) Optimal ensemble size 1* (eqs.
25
26 in the parameter space of sample size a and
0. Grey indicates k" 1 and white indicates k* - oo, with intermediate values given by the colorbar. Appended vertical bars show - and À indicated in panel titles. X denotes the 'locally optimal" regularization (see section H.3 (b) Optimal generalization error Eg. for the same parameter values in (a).
reduced readout noise scale setting W
45

---

a
p= 0, A =
p= 0.3, À =
p 0.6, À =
2.0
2.0
2.0
100
90
W
10-
1.0
80
60
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
0.0
1.0
20
3.0
α
40
p 0, À
=
0.3, À = -
0.6, À A
2.0
2.0
2.0
20
W
M
10-
1.0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
α
b
p
p 0.3,
p
0.6, À
2
2.0
1.0
10-
1.0
1.0
0.8
00
0.0
0.0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
α
α
p 0, = A
= 0.3, À = A
0.6, À
2.0
2.0
2.0
W
10
1.0
1.0
0.0
0.0
2.0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
α
Figure S7: (a) Optimal ensemble size k* (eqs. 25
26 the parameter space of sample size a and
O. Grey indicates k* 1 and white indicates k' - 3o, with intermediate values given by the colorbar. Appended vertical bars show - co. and À indicated in panel titles. X = X denotes the 'locally optimal" regularization (see section H.3 (b) Optimal generalization error Eg- for the same parameter values in (a).
reduced feature noise scale W setting H Z
46

---

a
p
0
0.3, À
p= 0.6, =
2.0
100
N 1.0
N 1.0
N 1.0
80
0.0
0.0
0.0
60
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
50
α
40
== O, =
p
0.3, À A
0.6, À =
30
2.0
10
N
1.0
N 1.0
N 1.0
0.0
0.0
0.0
0.0
1.0
2.0
3.0
2.0
0.0
1.0
2.0
3.0
α
α
b
==
0
0.3, =
0.6, À =
2.0
2.0
2.0
1.0
N
1.0
N 1.0
N 1.0
0.8
0.0
0.0
0.0
0.6
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
α
p
0, À = A
0.3, À
0.6, À 7
2.0
2.0
2.0
0.2
N
1.0
N 1.0
N 1.0
0.0
0.0
1.0
2.0
3.0
0,0
1.0
2,0
3.0
0,0
1.0
2,0
3,0
α
a
α
Figure S8: (a) Optimal ensemble size k* (eqs. and reduced label noise scale Z setting = W
25
26 in the parameter space of sample size
O. Grey indicates 4 = 1 and white indicates k* + 0o, with intermediate values given by the colorbar. Appended vertical bars show - and À indicated in panel titles. À = X denotes the 'locally optimal" regularization (see section H3 (b) Optimal generalization error Eg- for the same parameter values in (a).
47

---

p= 0, = 0.0
7
0, p = 0.0
0, p 0.6, 0.0
104
10A
104
Eg/s(1 - C)
100
10a
109
100
10-1
11
11
10-1
11
mM
10
10]
90
10
10
109
10
10
10
10
10a
10
10
100
10a
a
a
a
80
0, p
0,
H
1.0
A = 0, p 0.3,
H
1.0
A =
0, p
0.6,
H 1.0
70
10a
10a
10a
60
Eg/5(1
50
40
10-1
10-1
10-1
30
10-2
10-1
10a
10
10
10
10a
10a
10
10-1
10a
10a
20
α
0
p= H 2.0
0, p= 0.3,
2.0
0, p= 0.6
2.0
10
1
Eg/s(1 c)
10a
109
10a
11
100
10
10
10
10
10
10
10
10
10
10a
10
α
a
α
Figure $9: Reduced generalization errors Eg/s(1 c) with X = 0 and W = = 0 (given by eq. for linear ridge ensembles of varying size k. p and H values indicated above plots. Grey lines indicate 1, dashed black lines k + 0o, and intermediate & values by the colorbar.
22
48

---

À
p
0, H
0.0
== 0.3,
0 - p= H- 0.0
10a
10a
10a
Eg/s(1
10
10"1
10
100
10
10-1
10
10
10
10
10
10
10
107
10
10a
90
a
a
a
A ^=
p
H= 1.0
p= 0.3, H= 1.0 A =
p
0.6,
H
1.0
80
10a
10
10
70
Eg/s(1
60
50
40
10
10
10
10
10a
10a
10a
30
1
p
H
2.0
p
0.3, H
2.0 A
p 0.6, H
2.0
20
10
Eg/s(1 a
109
10
1
10-1
10-1
110g
10-1
10-2
10-1
10
10 10-2
10
10"
10a
10
10-1
109
10a
a
a
a
Figure S10: Reduced generalization errors Eg/ c) with = and W Z (given by eq. for linear ridge ensembles of varying size k. and H values indicated above plots. Grey lines indicate dashed black lines k oo, and intermediate & values by the colorbar.
23
49

---

A 0, 0- 0.0.0
À = 0, 0 = 0.0
À = 0, 0.6, =.0.0
10
104
10a
Eg/s(1
10a
10a
10
11
100
100
10-2
10-1
10
10a
10
10
10
101
10
10
109
10a
a
90
1
0, p
O, W
1.0
A
0, p 0.3, W
1.0
0,
0.6, W 1.0
80
Eg/s(1 - c)
101
10
10a
70
60
100
109
10c
50
10-1
11
mm)
40
10
10
10
10a
10
10
10a
10a
10
10
100
10a
30
α
4 = 0, p 0, W 2.0
A 0, p 0.3, W
2.0
A = 0,
0.6, W 2.0
20
101
10
10a
10
1
Eg/s(1
10a
10a
10a
10-1
10
mm
10
m0
mm
10
10
10a
10
10
10
10a
10
10
10
10a
101
α
a
Figure S11: Reduced generalization errors Eg/s(1 - c) with À = 0 and H Z 0 (given by eq. for linear ridge ensembles of varying size k. p and W values indicated above plots. Grey lines indicate = 1, dashed black lines k + co, and intermediate & values by the colorbar.
22
50

---

A = p 0, 0.0.0 Λ= 0.3, W = p= 0.6, W == 0.0
10a
10
10
Eg/s(1
100
10
10
110g
10
10
10
10a
10a
90
P = O, W
1.0
==
P 0.3, W
1.0 1 P 0.6, W
1.0
80
70
Eg/s(1 -c
109
10a
10
60
50
40
10-1
10p
10
10-1
10
10
10-2
10
10
10a
10
10-1
109
10a
30
a
a
a
7 =
= 0, W = 2.0
== p= 0.3, == 2.0 X" P 0.6, W = 2.0
20
10
Eg/s(1 c)
10a
10a
10a
10
10
10
10
10
10
10
10
10
10
10a
10
10
10a
10a
α
Figure S12: Reduced generalization errors Eg (s(s c) with İ = X and H = 2 0 (given by eq.
23 for linear ridge ensembles of varying size k. p and W values indicated above plots. Grey lines indicate k = 1, dashed black lines k + 0o, and intermediate k values by the colorbar.
51

---

0,
p = 0, 0.0
0, p
0.3, Z
0.0
0, p
0.6,
0.0
10
10a
10a
Eg/s(1 -
10
10a
10
10
10-1
10-1
100
10
10
10
10a
10
10
10c
10a
10
10
10
10a
0
90
A 0, p 0,
1.0
A
0, p 0.3, Z
1.0
0,
0.6,
1.0
80
10a
70
Eg/s(1
60
10
10a
100
50
10a1
[[|I]]
10
10"1
11
1110
1110
40
10
10
10
10a
10
10
10
10
10
10
10
10
30
a
α
4 0, p 0,
Z 2.0
À 0, p 0.3,
Z 2.0
À 0, 0.6,
Z
2.0
20
10
101
10a
10
1
Eg/s(1
10a
109
10"
10-1
10-1
mm)
10-1
10g
10-2
10-1
10
10
10-2
10-1
10a
10
10
10-1
109
101
α
a
a
Figure S13: Reduced generalization enrors Eg/s( with and H W (given eq.
22 for linear ridge ensembles of varying size k. and Z values indicated above plots. Grey lines indicate dashed black lines & oo, and intermediate & values by the colorbar
52

---

À
7 p 0, Z 0.0
== p = 0,0 À = P 0.6, Z = 0.0
Eg/s(1 - -c)
10
10a
10a
100
10
11
PTTI
T1/M
10
11
TTIM
10
11
TI/
10
10
10
10
10
10
10
10:
10
10-1
10
10a
90
a
a
7
λ P
0,
Z
1.0
A =
1
p 0.3,
Z 1.0 A
= P
0.6, Z 1.0
80
70
Eg/s(1 -)
10
10
10a
60
50
40
10
10a1
10
10
10
10a
10
10
10
10
10
10
10
10a
30
a
1
P
Z
2.0
1
==
0.3, Z
2.0 A
p
0.6, Z
2.0
20
10
10
10
10a
1
Eg/s(1
10-1
10-1
10
10
2
10
10
1
10
10
10
1
10
101
Figure S14: Reduced generalization errors Eg/s(1 c) with X = A and H W = 0 (given by eq.
23 for linear ridge ensembles of varying size k. p and Z values indicated above plots. Grey lines indicate k = 1, dashed black lines k + 0o, and intermediate k values by the colorbar.
53


Reviews:
[{"review_id": "IrSVjUEbKj", "reviewer": "Reviewer_tCNd", "summary": "This paper provides an asymptotic analysis of ensembles of ridge regressors using varying numbers of subsampled features. The authors consider a Gaussian data model with feature noise and readout noise in addition to label noise. Using the replica method from statistical physics, they obtain precise asymptotics of the limiting ensemble risk for any collection of subsampling operators. They specialize this result in the case of globally correlated isotropic data and numerically plot double descent curves for various non-overlapping ensemble strategies. They find that in poorly regularized regimes, increased heterogeneity of feature subset sizes leads to a significant decrease of worst-case generalization risk, in other words dampening the double descent peak, albeit typically at the expense of increased risk in other data aspect ratios. Departing from heterogeneity to consider instead homogeneous non-overlapping ensembles, the authors further demonstrate that in the presence of readout noise, there are regimes in the phase space of readout noise level and global feature correlation in which it is optimal to use an ensemble rather than a single learner, although without readout noise it is never optimal to ensemble.", "questions": "Instead of the replica method, do the authors think that the same result could be obtained by the more rigorous \"deterministic equivalences\" from random matrix theory [DS21]? Related work has recently successfully applied this technique to random projections, which should be directly applicable to ensembles [LPJBT22], which seems closely related to the future direction raised by the authors at the end of the paper.\n\nIn Figure 4.a.i, I cannot understand the meaning of the white dashed line. I thought it should indicate the boundary between when $k^* = 1$ and $k^* = 2$, but I cannot figure out why there is a purple region above the line, and why the caption says \"there is no region of the phase\nspace where $k^* = 2$.\" Another suggestion is to add markers in the (a) plots to indicate the points in phase space that correspond to the (b) plots.\n\n[DS21] E. Dobriban and Y. Sheng, Distributed linear regression by averaging, The Annals of Statistics, 49\n(2021), pp. 918 – 943, https://doi.org/10.1214/20-AOS1984.\n\n[LPJBT22] D. LeJeune, P. Patil, H. Javadi, R. G. Baraniuk, R. J. Tibshirani, “Asymptotics of the sketched pseudoinverse,” 2022, https://arxiv.org/abs/2211.03751\n\nOther remarks:\n- 1: \"Bagging\" is conventionally short for \"bootstrap aggregating\", but you only subsample, which is often called \"subbagging\" in the literature.\n- 2: Ensembles do reduce variance, but subsampling feature adds bias.\n- 4: \"the the\"\n- 212: \"eta\" should be $\\eta$", "limitations": "While I understand the usefulness of the replica method for obtaining quick and often accurate results, it is unfortunate that it is a non-rigorous approach. The authors have mentioned this in the supplemental material, but I think it needs mention somewhere in the main paper. This is important, as it is not clear whether some of the things done here are \"legal\"---for example, many random matrix theory results which could be applied here to perform similar analysis require covariance matrices with uniformly bounded operator norm, yet the globally correlated model violates this, which would typically require additional care to ensure that the result should hold, rather than simply plugging in the covariance. The reader deserves to know that while the general intuitions are likely true to hold, they should be careful in blindly applying the results.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The biggest strength of this paper in comparison to existing work on linear ensembles of random projections of features is that this analysis applies to arbitrary projection operators that need not be independent. In fact, for the specific cases the authors consider of non-overlapping feature subsets, which would correspond closely with distributed optimization practice (I recommend the authors emphasize this in their introduction), the projection operators are clearly very much not independent. Furthermore, the most general form of their results in Proposition 1 also makes no assumption on the relationship between the projection operators and the covariances, in contrast with most existing work which assumes them to also be independent.\n\nI am also not aware of any other work which explicitly considers heterogeneous linear ensembles, which is a novel contribution.", "weaknesses": "I think that the biggest weakness of this paper is the justification of the data model and of the advantage of using a (heterogenous) ensemble. I think that these could be improved by adding justifying and motivating comments throughout the paper, and if the authors would indicate how they would make such changes for the camera ready version, I would increase my score of the paper.\n\nData model:\n- The features are corrupted with additional noise that is independent of the label, which is an interesting but highly non-standard assumption. Even more non-standard is the assumption that this noise is drawn independently in every forward pass. This assumption really should be motivated---should the reader be imagining this as something like noise due to dropout in a neural network?\n- The model outputs are corrupted with readout noise independent across each member of the ensemble, again drawn independently each forward pass. This seems like an even stranger assumption to me; I could imagine some negligible computational noise, but I struggle to come up with a setting in which we would see the noise levels that are modeled in this paper. It is very important to motivate this, as the experiment in Figure 4 only justifies ensembles in settings where readout noise is high. Without good motivation, it seems like a strawman.\n\nEnsemble advantages:\n- This is a complaint I have for many other papers that deal with double descent, but applies equally to this paper. The \"correct\" solution to double descent is to properly regularize, and other strategies such as the ensembling proposed here can be seen as \"hacks\" that can mitigate some of the worst double descent effects but are typically suboptimal. We see this in this case as well. E.g., consider Figure 3, where the heterogenous ensembles only improve risk in the under-regularized setting. In the setting with more appropriate regularization (c.iv-vi), where the risk is lower than in (c.i-iii) for most values of $\\alpha$, the heterogeneous ensemble is strictly worse than the homogeneous ensemble, which I am sure itself is strictly worse than the optimally tuned full ridge regressor. \n- No comparisons are made to the full ridge regressor. Obviously, the ensemble would have to be worse, as ridge regression is typically the minimum mean squared error estimator in additive noise settings. However, this is a lost opportunity to demonstrate perhaps that bad double descent effects could be mitigated at a fraction of the computational cost (since the ensembles consider small subsets).\n- Essentially, the risks are a linear combination of the individual risks of each ensemble member, plus some cross-terms. Thus in the limit of zero regularization, these ensembles would still have infinte double descent peaks, one for each different subsampling ratio. So to an extent, the authors already rely on regularization to help make the ensemble robust to double descent, although the benefits are realized for fairly small regularization already. Still, the authors should justify the need to regularize even with ensembles and discuss tuning.\n- Related to the previous points, it's not clear why worst case risk is an important metric, since tuning of regularization in some way or another seems necessary.\n- Part of the issues of statistical performance may be due to the non-overlapping partition of features. In the related works, ensembles with optimality guarantees typically require many members with overlapping subsets of features, which the authors do not consider here. \n\nOthers:\nOverall, I found the plots to be very busy and difficult to decipher. I give one suggestion on how to improve Figure 4 in Questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for writing a thorough review of our submission. Please see the global response as well as the below.\n\nWe had two interpretations of the feature noise in mind when creating the problem setup.  The first is as an inherent noise due to stochasticity in a physical neural network, such as an analog neural network or a biological neural circuit.  The second is as a corruption of the input data. In both cases, the noise is present in both training and test data. Whether the noise is introduced through noise in the data itself or introduced through noise in the neural network computation is a matter of interpretation—but leads to an identical mathematical model.  If accepted, we will clarify this in the statement of and discussion following proposition \n\nThank you for your comment that regularization is the “‘correct’ solution to double descent,” which we agree with. We do not argue that heterogeneous ensembling is the best regularizer, but to introduce heterogeneous ensembling as a regularization method and study its properties. We will add to figure 3c comparisons with the learning curves at optimal L2 regularization.  We will also add the following sentences: “in the feature subsampling ensemble, we find that heterogeneity of this type smooths out the double-descent peak of the learning curve. However, when computational resources permit it, a single fully-connected readout layer with optimal l2 regularization outperforms a feature subsampling ensemble, unless additional sources of variance over the ensemble are introduced (see section 2.6).\nThank you for raising the interesting point that error diverges at zero regularization.  Note that the  genuine divergence of the generalization error is an artifact of the infinite-dimensional limit.  We also show in the supplemental material that heterogeneous subsampling can mitigate double-descent without any regularization in a classification setting.  The linear regression model with small positive regularization serves as a fruitful demonstration of the benefits of heterogeneous ensembling.  We will add that “in the infinite-dimensional regression setting, a small regularization  is always necessary to prevent divergence of the generalization error at the interpolation threshold of each ensemble member.”  We will also add to figure 3c comparisons with the learning curves at optimal regularization.. \n\nIn a regression setting, where regularization is necessary to prevent divergence of the generalization error, we propose the following procedure to quickly obtain a robust predictor:\nPick a reasonable but small regularization parameter.\nTrain an ensemble of linear predictors on varying numbers of data features.\nThrough this procedure, one can avoid the sharp edges of double-descent without a computationally expensive parameter sweep.  We will update our discussion to make this clearer.  In the CIFAR10 classification setting, we demonstrate that heterogeneous ensembling smooths out the double-descent curve at zero regularization, as in a classification task percentage error is naturally bounded. \n\nWe have since updated our calculations for proposition 2 to consider overlapping partitions of features.\n\n## On Using the Replica Trick\n\n We do in the remarks made directly after the statement of our main result in the main text state that the replica trick is a “non-rigorous but standard” approach.  We make no claim that our results are fully rigorous, but because the results obtained from the replica trick have been demonstrated to coincide with rigorous results for many similar problems and our analytical results show excellent agreement with numerical simulations, we have good reason to expect that our result is correct. However, In light of the concerns raised here we see it fit to change the statement of our main result to further emphasize that the replica trick is not rigorous, and to establish a rigorous basis for our main result.\n\nIn the case where the data covariance matrix has a bounded spectrum, we believe  that our results may be obtained through a clever special case of the rigorous result of Loureiro [2022].  While their derivation does not explicitly consider ensembles with variation in the number of features viewed by each ensemble member, this type of heterogeneity may be added in post-hoc by choosing data covariance matrices which “zero out” a number of features which varies over the ensemble.  If accepted, we will add a supplemental section discussing this correspondence and including a detailed derivation of our result from the general result of Loureiro [2022].  We will also update the “proof” of the main theorem to the following:\n\n“We calculate the terms in the generalization error using the replica trick,  a standard but non-rigorous method from the statistical physics of disordered systems.  The full derivation may be found in the SI.  In the special case where the covariance matrices $\\Sigma_s, \\Sigma_0$” have bounded spectrum, this result may be obtained as a clever special case of the results of Loureiro [2022] (see SI for derivation).”\n\nWe will also add the following remark after the proof of proposition 2:\n“Note that, as in this case $\\Sigma_s$ does not have a bounded spectrum, this result does not follow from the rigorous results of Loureiro [2022].  However, we find excellent agreement between theory and experiment when data dimension is sufficiently large.”\n\nWe emphasize that, even as our general result may be recovered as a special case of the results of Loureiro [2022], simplifying this general result in the special case of subsampling from globally correlated features is a very tedious calculation which requires significant work.  Our investigation of this general result also differs entirely from Loureiro [2022], which considered only ensembles with the same number of readout weights across ensemble members, and included no study of feature noise or readout noise."}, {"review_id": "WUlyx2OLtO", "reviewer": "Reviewer_JR3L", "summary": "This article characterizes the asymptotic performance curve of a heterogeneous feature ensembling framework for ridge linear regression, in the limit of comparably large numbers of data samples and variables. For having different error peaks in a double-decent performance curve as a function of the data sample size, predictors built of heterogeneous feature sets give rise to an ensemble that is robust to the number of data samples, without a carefully tuned regularization parameter.", "questions": "I wonder if the optimal learning performance is achieved at a uniform subsampling rate (i.e., homogeneous ensembles); and whether the feature ensembling can indeed improve the performance without readout noise, as is the case with readout noise (and under certain conditions on the data model parameters) according to the discussion of Section 2.6.", "limitations": "The limitations are discussed. There does not seem to be any negative societal impact.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "* Clarity of the mathematical setup.\n* Numerical validation of asymptotic results on finite data sets.\n* Discussion of the theoretical implications of the analysis.", "weaknesses": "* The analysis applies only to the square loss while recent related work addressed already a large family of convex losses.\n* Even though the interest of heterogeneous feature ensembling is explained as mitigating the double-descent phenomenon, this problem can be in fact settled by optimizing the regularization parameter on a cross-validation set.\n* There is no experiment on real data sets to see how well the observed consequences such as the mitigated double-descent curve apply to real data.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for writing a helpful review.  Your questions about cross-validated regularization and real-world experiments are addressed in the global rebuttal and attached PDF.\n\n"}, {"review_id": "oU0JHKg4OO", "reviewer": "Reviewer_gGcm", "summary": "This paper introduces a theoretical investigation into ensembling methods applied to linear ridge regression with feature subsampling. It builds upon previous research in this area by extending the analysis to include scenarios with varying readout dimensionality.\n\nBy employing the replica method, the authors derive expressions for the average-case generalization error. These expressions capture the influence of both the data structure and the hyperparameters of the ensembling method employed.\n\nFurthermore, the paper focuses on a simplified scenario to examine how the degree of subsampling and heterogeneity affect the occurrence of the double descent phenomena in learning curves for generalization error. This analysis provides potential insights into the understanding of this well-known and recurring phenomenon in machine learning.", "questions": "- Would the authors consider revising the presentation of their results from Proposition 1 to reflect the heuristic nature of the findings? This would involve explicitly stating that the results are derived using a heuristic method, thereby providing clarity to readers about the level of mathematical rigour employed in the study. Such a revision would enhance transparency and enable readers to appropriately interpret the results.", "limitations": "While the authors do acknowledge certain limitations in their work and propose potential new directions for applying the techniques employed, there are a couple of crucial aspects that remain unaddressed. Firstly, the paper fails to explicitly mention that the replica method utilized lacks mathematical rigour, which is an important consideration for interpreting the validity and reliability of the results. Additionally, the authors do not provide any indications or suggestions on how to potentially establish a more rigorous foundation for their findings. ", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The paper is clearly written, and the experiments presented are of high quality.", "weaknesses": "- The authors use as proof technique for the main theorem the replica method, a heuristic method coming from statistical physics. While this method has been proven to be rigorous in certain contexts, its general correctness has not yet been formally established.  Unfortunately, the paper fails to acknowledge this limitation, which undermines the overall credibility of the research findings.\n- In terms of technical results, the computation leading to preposition 1 appears to be a minor modification of the calculation presented in a previous work by Loureiro [2022]. However, unlike the present study, Loureiro's work accompanies the calculation with a rigorous proof. This discrepancy raises concerns about the rigour and robustness of the current paper's technical contributions.\n\n(Minor) Typos:\n- Line 4: the the\n- Line 53: studied studied\n- Line 212: eta", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for pointing out typos, we will correct them.  We thank the reader for their comments and suggestions.  The weaknesses identified in this review revolve mainly around the fact that we have derived our results using the replica method from statistical physics. We do in the remarks made directly after the statement of our main result in the main text plainly state that the replica trick is a “non-rigorous but standard” approach.  We make no claim that our results are fully rigorous, but because the results obtained from the replica method have been demonstrated to coincide with rigorous results for many similar problems and our results show excellent agreement with numerical experiments,, we have good reason to expect that our result is correct. However, in light of the concerns raised here we see it fit to change the statement of our main result to further emphasize that the replica trick is not rigorous, and to establish a rigorous basis for our main result.\n\nIn the case where the data covariance matrix has a bounded spectrum, we believe  that our results may be obtained through a clever special case of the rigorous result of Loureiro [2022].  While their derivation does not explicitly consider ensembles with variation in the number of features viewed by each ensemble member, this type of heterogeneity may be added in post-hoc by choosing data covariance matrices which “zero out” a number of neurons which varies over the ensemble.  If accepted, we will add a supplemental section discussing this correspondence and including a detailed derivation of our result from the general result of Loureiro [2022].  We will also update the “proof” of the main theorem to the following:\n\n“We calculate the terms in the generalization error using the replica method,  a standard but non-rigorous method from the statistical physics of disordered systems.  The full derivation may be found in the SI.  In the special case where the covariance matrices $\\Sigma_s, \\Sigma_0$ have bounded spectrum, this result may be obtained as a clever special case of the results of Loureiro [2022] (see SI for derivation).”\n\nWe will also add the following remark after the proof of proposition 2:\n“Note that, as in this case $\\Sigma_s$ does not have a bounded spectrum, this result does not follow from the rigorous results of Loureiro [2022].  However, we find excellent agreement between theory and experiment when data dimension is sufficiently large.”\n\nWe emphasize that, even as our general result may be recovered as a special case of the results of Loureiro [2022], simplifying this general result in the special case of subsampling from globally correlated features is a very tedious calculation which requires significant work.  Our investigation of this general result also differs entirely from Loureiro [2022], which considered only ensembles with the same number of readout weights across ensemble members, and included no study of feature noise or readout noise."}, {"review_id": "zH0qGwXYST", "reviewer": "Reviewer_QQ7H", "summary": "The authors provide theoretical results on generalization for ridge regression for the case of an ensemble of regression models trained on feature subsets, with noise, and correlation.\n\nThey relate this to the previously observed double-descent phenomena and characterize how the fraction of subsampled features, the regularization strength, and other properties affect the double descent phenomenon and where it occurs (i.e., how these impact generalization error).  \n\nFinding that double descent occurs at different places for different feature subset amounts, they also propose and analyze using a heterogenous ensemble of different feature subsampling sizes as a way to mitigate double descent without having to tune regularization or model architectures.\n\n\n****Update****\nRevised my score after the responses - see my response comment for rationale.\n", "questions": "\n\"...correlated data by ensembles with heterogeneous readout dimensionality.\"\nThis is not very clear - what is meant by \"readout dimensionality\"?  I'm assuming\nit means heterogeneous input data dimensionality - so why not just say that?  \"Readout dimensionality\" is a term I've almost never heard and sounds like it mean the dimension of the model output, which doesn't make sense in this case.\n\n\nIt would be helpful to explain what \"readout noise\" means in this paper when first introduced - this was not clear to me and I found its use to be confusing in the theoretical setup section.  It's just mentioned without defining what it is or how it related to the model, and it would help to give an example of what it would correspond to.  After several re-reads I concluded it probably corresponds to noisy model outputs during training when using something like dropout.  However, for the most cases in practice for regression, even if noise is used in training, when applying the models in practice deterministic outputs are used (no longer any \"readout noise\").\n\n\n2.6 is hard to follow - again what is mean by each type of \"readout\"?  What is a \"single fully-connected readout\" and what does it mean by \"multiple sparsely connected readouts\"?  I'm assuming the latter is the feature subsampling ensemble, and the former a single model with no feature-subsampling, but this is described in a confusing way so it's not clear and took some though to try to decipher this.\nWhat is mean by correlation code and correlation strength?  How are these defined, and how are they used to generate the models / how are they used in these analyses and results?\n\nAlso how do these conclusions differ from past studies of ensembling approaches?  \n\n\n\n\nWhy in all the examples is P <= M?  It would be nice to explain this as well, and also what happens when P > M.  I.e., after some thought I assume this is to focus on the over-parameterized regime, but this is never really stated and explained, and there are still many practical cases where P >> M.\n\nWhat impact does k have on generalization error?\n\nAlso wondering how the results compare to the non-ensemble case?  I.e., a single model with no feature subset sampling.\nI.e., this shifts the problem to tuning the subset sample fraction sizes\n\n\nAlso it's not clear to me from the results and figures if tuning lambda (regularization) can always result in better performance then using the heterogenous ensemble?\n\n\n\n", "limitations": "Yes - doesn't seem needed in this case.", "rating": 5, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 2, "strengths": "The generalization bounds for the particular setting seem novel.  \n\nThe problem is a highly relevant one and the idea of using feature subset ensembles is interesting and potentially useful- as alluded to in the conclusion - as searching over many architecture and regularization settings for large neural nets may be prohibitive, so this kind of approach could provide an alternative for such cases with a fixed architecture, by creating a feature subsample ensemble instead - however this possibility is just alluded to and not tested, as only theoretical results for simple ridge regression problems are analyzed.\n\nThe authors provide thorough analyses of the results of the theorem for generalization error, to try to elucidate what different aspects of the ridge regression problem and choices for the modeling mean in terms of generalization error.", "weaknesses": "1) The paper is hard to follow and extract any useful meaning or conclusions from.\n\na) Some terms are not defined or given adequate explanation.  As particular example - it would be helpful to explain what \"readout noise\" and readout dimensionality mean in this paper when first introduced - this was not clear to me and I found its use to be confusing in the theoretical setup section and even the related work.  It's just mentioned without defining what it is or how it related to the model, and it would help to give an example of what it would correspond to.\n\nb) A lot of different parameters are introduced and hard to remember when trying to connect to the various results figures and analyses -  it may help to remind the reader and also include simple examples of models describing what the different parameters and values would correspond to.\n\nc) I feel the main results and conclusion need a simpler explanation and deciphering for the reader - as the various plots and analyses are still based on a large number of symbols / parameters that it's hard to remember, and the reader would have to dig through the text each time to recall what they are.  Some simplified example and conclusion, as well as a table may help.  It's hard to come away with any kind of recommendation or firm conclusion of what we could expect to work better in practice as well.\n\nd) \n\n2) It's not clear how useful and impactful the results from the paper and the derived theorems are.  \nIt seems the complicated expression derived is too complex to derive any useful properties, and requires knowing ground truth properties - so the benefit is not clear.\n\nEven in the simplified tractable case - it's not clear what interpretation comes out of it.\n\nIn generally it's not clear if any novel and useful findings have come out of this.\n\nI think a key part of the issue might be - after presenting each set of results / figures, some interpretations and conclusions and take-ways for the reader are needed and currently lacking - results are just presented with little explanation or analysis.\n\n\n3) Related to the previous points, I feel having real data and model experiments and results could really strengthen this work and also more fully explore the idea and benefit proposed of using heterogenous feature sub-sampling ensembles.  \n\n\n4) Novelty is also somewhat unclear.\nObviously, heterogenous feature subsampling ensembles have been studied and applied in the past, and for some specific models there are theoretical results on generalization error as well, such as tree ensembles or even linear model ensembles with different kinds of approaches such as random subspace models.  \nAdditionally, the authors point to some specific related work as well on theoretical study of generalization error as well.  However, clearly spelling out the differences in this work seems lacking.  For example, in related work - they mention other work also studying feature subsampling for ensembles and finding it provides an implicit regularization approach, in particular [20].  However, it's not stated what differentiates this work from this prior work - what novel aspect does this provide?  It's necessary to spell out what limitations of previous work this addresses and overcomes.\n\nI.e., overall it would help to clearly spell out how their work is different and novel with respect to prior work on feature-subsampling ensembles and what new conclusions and approaches result from it.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for writing a thorough review of our submission.  Please see the response to all reviewers for a discussion of the meaning of “readout noise”, a proposed table clarifying the meaning of the parameters referenced in proposition 2 and figures 2, 3, and 4, discussion of a comparison between ensembling and L2 regularization, and a discussion of applications to real data.\n\nThank you also for pointing out that “readout dimensionality” is confusing terminology.  This confusion may be due to the fact that in linear regression the number of input features is equal to the number of model parameters, which we referred to as the “readout dimensionality.”  If accepted, we will remove this term from the paper, and instead refer to the “input data dimensionality” as suggested.  We clarify that the model output is one-dimensional for the setup described in propositions 1 and 2 and all figures of the main text.  Our qualitative finding that heterogeneous input data dimensionality can mitigate double-descent carries over to simulations of the CIFAR10 classification task where the output dimensionality is 10, as shown in the SI.\n\nOn terminology – You have assumed correctly that a “single fully-connected readout”  is a single model with no feature-subsampling, and “multiple sparsely connected readouts” is the feature subsampling ensemble.  We will update this sentence to the following: “We now ask whether ensembling is a fruitful strategy – i.e. whether a feature subsampling ensemble outperforms a single model without feature subsampling.”\n\nThank you for prompting us to clarify the novelty of our work. We introduce and study a new form of heterogeneous subsampling, and demonstrate its benefits in an analytically tractable setting.  Specifically, we consider ensembling over models with heterogeneity in the number of features connected by each ensemble member.  All prior research that we are aware of, including [20], considers only the setting in which all ensemble members have the same size.  We emphasize this distinction in the “Related works” section, writing: “However, [18] and [19] focus their analysis on the case of isotropic data and Gaussian random masks of homogeneous dimensionality.  In contrast, we explicitly consider learning from correlated data by ensembles with heterogeneous [input data dimensionality].” More work is required to establish and study the practical benefits of this form of heterogeneity.  If accepted, we will acknowledge this more clearly in the final version of the paper.  This is also the first work we are aware of which presents a detailed study of the interplay between ensembling and readout noise, or to produce phase diagrams of the type shown in figure 4.\n\nThank you for pointing out that Section 2.6 is hard to follow.  This will be re-written to reflect an updated version of the special case in proposition 2.  The parameters of the phase diagrams in figure 4 will be re-interpreted in terms of an effective signal-to-noise ratio and data-task alignment.  Please see the response to all reviewers for an overview of these developments and the attached updated version of figure 4.\n\nAlso how do these conclusions differ from past studies of ensembling approaches?\nPast studies of ensembling approaches have typically aimed to establish formal equivalences between optimal ensembling and optimal regularization.  The optimal subsampling fraction, though, depends on the size (number of samples) of the training set.  In this work, we present heterogeneous ensembling as a strategy to mitigate double-descent in settings where the size of the training set is unknown.  We are the first study to explicitly analyze the behavior of ensembles with heterogeneous input data dimensionality and to recognize the benefits of this type of heterogeneity to mitigate double-descent.\n\nTo clarify why P<=M in all the figures, we will add the following sentence to the discussion:\n“As ensemble members are trained independently, the learning curve cannot be divided into an under-parameterized and over-parameterized regime.  Rather, at any sample size $P$, each ensemble member with $N_r<P$ is in its over-parameterized regime and each ensemble member with $N_r>P$ is in its under-parameterized regime.”\nBecause of this, the plots shown cover both the under-parameterized and over-parameterized regimes of each member of the ensemble (since $N_r \\leq M$ for each ensemble member)\n\nTo clarify the impact of k on generalization error, we will add the following sentences to the discussion of figure 4 at the end of section 2.6:\n“These plots show the effect of the ensemble size $k$ on the generalization error curves for ‘homogeneous ensembles’, where each ensemble member sees the same number of features, and in the absence of regularization.  We find that increasing $k$ may either help or hurt generalization performance, depending on the parameters of the task, and the size of the training set.  The shape of the learning curve changes with $k$ to shift the double-descent peak to the value $\\alpha = \\frac{1}{k}$.”\nHowever, readers may still wonder how the value of $k$ affects generalization at finite ridge, and in heterogeneous ensembles.  If accepted, we will clarify this by adding panels to figure 3 which plot learning curves for ensembled linear regression at finite ridge $\\lambda$ for ensembles of a handful of $k$ values, as we do for the CIFAR10 classification problem in figure S1.\nTo the question:\n*Also wondering how the results compare to the non-ensemble case? I.e., a single model with no feature subset sampling. I.e., this shifts the problem to tuning the subset sample fraction sizes*\nCalculating the optimal subset sample fraction sizes is precisely what we do in the phase diagrams in figure 4.  Please see the response to all reviewers for an updated version of figure 4 and subsequent analysis.\n"}, {"review_id": "19xxsrHkQ0", "reviewer": "Reviewer_uQUq", "summary": "The authors provide a theoretical analysis of the case of ensemble learning with linear ridge regression for the case where heterogeneous feature subsampling is used. The authors make a number simplifying assumptions about the distribution of the data (Gaussian distribution and noise, linear function), and using the replica trick from statistical physics derive generalization error for the L2-regularized least-squares solution. In simulations, the authors demonstrate that the derived solution coincides with numerically calculated generalization error, and based on analyzing learning curves provide a number of insights about the behavior of ensembled ridge regression, suggesting a novel way for mitigating the double descent phenomenon through the use of heterogenous feature subsampling in ensembling.", "questions": "How computationally demanding is the computation of the generalization error curves?", "limitations": "N/A", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The work provides novel insights about how double-descent behaves based on noise, regularization, and subsampling and provides alternative strategies to regularization for avoiding the peak. These are quite central questions in modern machine learning theory and new ideas here can prove to be significant. The mathematical analysis seems rigorous and the replica trick used for the analysis appears to be a new type of tool for studying this phenomenon, though I do not have suitable mathematical background to follow all the proofs in detail. The basic setting and assumptions are clearly defined, related work seems to be covered well, and novelty of the contribution compared to related similar works is clearly established.", "weaknesses": "The work makes many simplifying assumptions that are unlikely to match most real-world learning problems (e.g. Gaussian data, linear teacher function), and analyzes may not be in practice tractable beyond the simpler cases such as analyzed in Section 2.3. with globally correlated data. It is not so clear whether the results would offer practical tools for an analyst wishing to apply heterogeneous feature-subsampled ridge ensembles in an optimal way, though results in the supplementary materials about CIFAR10 Classification Task suggest the results can have also practical relevance.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*Weaknesses:\nThe work makes many simplifying assumptions that are unlikely to match most real-world learning problems (e.g. Gaussian data, linear teacher function), and analyzes may not be in practice tractable beyond the simpler cases such as analyzed in Section 2.3. with globally correlated data. \nIt is not so clear whether the results would offer practical tools for an analyst wishing to apply heterogeneous feature-subsampled ridge ensembles in an optimal way, though results in the supplementary materials about CIFAR10 Classification Task suggest the results can have also practical relevance.*\n\nWe thank the reader for raising this point.  To address this weakness, we will add the following sentences to the discussion section:\n“A large  line of work has shown that Gaussian data approximations accurately predict the learning curves in many real world problems, particularly when using a least-squares loss function (see, for example, Pesce [2023], Canatar [2021], Gerace [2022], Hu [2022]).” \nReferring to the following papers:\nPesce [2023]: https://doi.org/10.48550/arXiv.2302.08923\nCanatar [2021]: https://www.nature.com/articles/s41467-021-23103-1\nGerace [2022]: https://arxiv.org/abs/2205.13303\n Hu [2022]: https://arxiv.org/abs/2009.07669\n\nWe have also performed our own realistic simulations which are presented in the supplemental material, demonstrating that our main insights carry over  to the CIFAR10 classification task.  We will add the following sentences  to section 2.5 of the main text to highlight these findings:\n“In figure S1, we train linear models to predict the labels of images from the CIFAR10 dataset from subsamples of the features in the top hidden layer of a pre-trained deep neural network.   We find that heterogeneous ensembling (over linear models with varying input dimensionality) prevents catastrophic over-fitting, yielding monotonic learning curves without the need for regularization.  (see SI for details)”\n\n*Questions:\nHow computationally demanding is the computation of the generalization error curves?*\n\nAll of the computations necessary for this paper can be performed in a matter of hours on a single GPU.  Computing the theoretical generalization error curves can be expensive (taking many minutes) in the general case where the result is given as a contraction of matrices, and where solving the fixed point equations for the order parameters requires inverting a large matrix at each step.  However, in the special case of globally correlated (equicorrelated) features, we obtain accurate, fully analytical expressions which can be plotted instantly. If accepted, we will add a note about the computational complexity of the theoretical and empirical calculations of the error curves to the paper.\n"}]

Author rebuttal:
Thank you for your insightful reviews.  Please find below a description of updates to the paper and responses to comments which were raised in multiple reviews.

-We have changed the definition of the “readout noise” so that it is present both during training and evaluation of the model.  This leads to a minor change in the form of the generalization error. Now, the contribution due to the readout noise goes as $\langle E_{rr’} \rangle \sim \frac{\delta_{rr’} \eta_r^2}{1-\gamma_{rr’}}$

-We have changed “globally correlated” to “equicorrelated” to agree with previous literature.

-We have updated proposition 2 to be more general, including an overall scale of the data, isotropic feature noise, and allowing for overlap in the features sampled by different readouts in the ensemble.  In the updated statement, we simplify the general expression for the generalization error from proposition 1 under the following special case:
>$$ \mathbf{w}^* = \sqrt{1-\rho^2} \mathbb{P}_{\perp} \mathbf{w}^*_0 + \rho\mathbf{1}_M $$
$$ \mathbf{w}^*_0 \sim \mathcal{N}(0, \mathbf{I}_M) $$
$$   \mathbf{\Sigma}_s = s \left[(1-c) \mathbf{I}_M + c \mathbf{1}_M \mathbf{1}_M^\top \right] $$
$$ \mathbf{\Sigma}_0 = \omega \mathbf{I}_M $$

-To enhance readability, we will use the extra space to include a table of parameters relevant to proposition 2 and figures 2, 3, and 4.  This table is included in the PDF attachment.

-Figure 4 has been updated (see attached PDF) and the surrounding discussion will be updated as follows:

>In the ridgeless limit, we can then express the error as : $E_g(k) = s(1-c)F(H, k, \rho, \alpha)$, where 
$H \equiv \frac{\eta^2}{s(1-c)}$ is an effective inverse signal-to-noise ratio and 
$F(H, k, \rho, \alpha)$ is a rational function of its arguments (see SI for full expressions).  Thus the value 
$k^*$ which minimizes error depends on 
$\eta$, 
$s$, and 
$c$ only through the ratio 
$H$.  
>Using our analytical theory, we plot the optimal number of readouts 
$k$ in the parameter space of 
$H$ and 
$\rho$ (see Fig. 4a).  The resulting phase diagrams are naturally divided into three regions.  In the signal-dominated phase a single fully-connected readout is optimal ($k^* = 1$).  In an intermediate phase, $1<k^*<\infty$ minimizes error.  And in a noise-dominated phase $k^* = \infty$...  As is evident in these phase diagrams, an increase in H (decrease in SNR or increase in c) or an increase in $\rho$ causes an increase in $k^*$.

-We will clarify our motivation for considering a “readout noise” by adding the following paragraph to the introduction:
>“Subsampling from different sets of features introduces variance between members of an ensemble.  However, other sources of variance may also distinguish members of an ensemble of predictors.  For example, random initialization or weight dropout during training may introduce variance in an ensemble of deep networks trained with gradient descent.  Physical neural networks may also have intrinsically noisy neurons which introduce variance between ensemble members. To capture these effects which are not naturally present in the regression setting, we introduce an  explicit “readout noise” which is drawn independently for each ensemble member.   We provide a detailed analysis of the interplay between feature subsampling ensembles and readout noise in the special case of equicorrelated features.”
We will add the following remark after the statement of our main result:
“Remark 4:  “Readout noise” is a noise which is added to the prediction of each member of the ensemble of readouts before they are averaged.  The readout noise parameter $\eta$ represents all additional sources of variance across predictors not accounted for in our simple model.  For example, in an ensemble of deep neural networks trained with gradient descent, there will be additional variance across predictors due to random initialization of the networks (Atanasov 2023).  Additional variance may also be present in physical neural networks, such as an analog neural network or biological neural circuits due to inherent stochasticity in the physical mechanisms of computation.”
And by adding the following sentence the discussion of Figure 4 in section 2.6:
“The resulting phase diagrams demonstrate that partitioning data features amongst multiple readouts can be beneficial when there are sources of variance between ensemble members beyond that induced by sampling from separate sets of features.”

-On the need for real data and model experiments – we have included in the supplementary material an application of heterogeneous ensembles to the CIFAR10 classification task.  While the classification setting does not correspond to the regression setting of our main theorem, the qualitative behavior carries over to this real-world task – we observe that heterogeneity in the number of features seen by each member of the readout ensemble mitigates double-descent even without regularization.  We plan to significantly expand this portion of the paper to demonstrate the benefits of heterogeneous ensembling for image classification with pre-trained feature maps (for example,  a pre-trained ResNet).

-A single fully connected readout with optimal L2 regularization does perform better than a feature subsampling ensemble. However, training a feature-subsampling ensemble is much faster than training a single fully-connected readout layer. (Assuming cubic scaling of the matrix inverse operation, training $k$ models of size $M/k$ scales as  $1/k^2$). Further, introducing heterogeneity of the type we study into the ensemble gives the added benefit of avoiding catastrophic over-fitting without a task-tuned regularization.  This approach is most likely to find applications in settings where training time or computational resources are limited.  If accepted, we will add these points to our discussion.  We will also add curves to panels c.i-vi which show the learning curves for a single model with optimal regularization.

Meta-review: The paper analyses the asymptotic properties of ensembles of ridge regressors trained from heterogeneous (in terms of content and size) feature subsamples. The paper combines solid theory and some interesting preliminary experiments that support this theory. The reviewers highlighted several limitations of the work (justification of the data model and presentation problems) that were partially addressed by the rebuttals. Despite these limitations, the reviewers and I agree that the work is novel, can certainly generate interesting discussions and lead to further theoretical or more practical developments in the future. I therefore recommend acceptance. The authors are urged to make all promised changes in the final version of their paper.
Decision: Accept (poster)