# Bias-variance decomposition of overparameterized regression with random
linear features

Jason W. Rocks  Department of Physics, Boston University, Boston,
Massachusetts 02215, USA  Pankaj Mehta  Department of Physics, Boston
University, Boston, Massachusetts 02215, USA  Faculty of Computing and Data
Sciences, Boston University, Boston, Massachusetts 02215, USA

###### Abstract

In classical statistics, the bias-variance trade-off describes how varying a
model’s complexity (e.g., number of fit parameters) affects its ability to
make accurate predictions. According to this trade-off, optimal performance is
achieved when a model is expressive enough to capture trends in the data, yet
not so complex that it overfits idiosyncratic features of the training data.
Recently, it has become clear that this classic understanding of the bias-
variance must be fundamentally revisited in light of the incredible predictive
performance of “overparameterized models” – models that avoid overfitting even
when the number of fit parameters is large enough to perfectly fit the
training data. Here, we present results for one of the simplest examples of an
overparameterized model: regression with random linear features (i.e. a two-
layer neural network with a _linear_ activation function). Using the zero-
temperature cavity method, we derive analytic expressions for the training
error, test error, bias, and variance. We show that the linear random features
model exhibits three phase transitions: two different transitions to an
interpolation regime where the training error is zero, along with an
additional transition between regimes with large bias and minimal bias. Using
random matrix theory, we show how each transition arises due to small nonzero
eigenvalues in the Hessian matrix. Finally, we compare and contrast the phase
diagram of the random _linear_ features model to the random _nonlinear_
features model and ordinary regression, highlighting the new phase transitions
that result from the use of linear basis functions.

##  I Introduction

One of the core concepts in modern statistics and supervised learning is the
bias-variance decomposition. It states that the test error, the predictive
performance of the model on new data, can be decomposed into three parts:
bias, variance, and noise Geman  _et al._ (1992); Bishop (2006); Mehta  _et
al._ (2019). The bias captures errors due to underfitting, resulting from the
inability of a statistical model to sufficiently express statistical
relationships present in the data distribution. The variance, on the other
hand, characterizes errors that result from “over-fitting” unrepresentative
aspects of the training data set that do not generalize (e.g., label noise).
Finally, the noise describes irreducible errors in a test data set due to
randomness in the data generating process.

In classical statistics, the bias-variance trade-off suggests that optimal
predictive performance is achieved by utilizing statistical models with
intermediate model complexities, balancing errors due to bias and variance.
While increasing a model’s complexity (e.g., increasing the number of fit
parameters) reduces bias, it comes at the price of increasing variance. One of
the most interesting and surprising empirical results to emerge from deep
learning over the last five years is the realization that this basic intuition
is fundamentally incomplete; it does not apply to “overparameterized” models
where the number of fit parameters is large enough to perfectly fit the
training data (i.e. achieve zero error on the training data set) Zhang  _et
al._ (2017).

While the classic bias-variance trade-off still holds in the
underparameterized regime (i.e., for models that have too few fit parameters
to achieve zero training error), once a model’s complexity is increased passed
the interpolation threshold – the point at which the training error goes to
zero – the test error once again decreases. The resulting combination of a
“U-shaped” test error in the underparameterized regime and the subsequent
decrease in test error in the overparameterized regime is now commonly
referred to as a “double-descent” curve Belkin  _et al._ (2019); Loog  _et
al._ (2020). This double-descent behavior seems to be a generic property of
all overparameterized supervised learning models and for this reason, has
become a major area of research.

An important open question in the field is to understand the double-descent
phenomena in terms of classical ideas of bias and variance. One fruitful
approach has been to analyze analytically tractable models that exhibit the
double-descent phenomena Adlam and Pennington (2020); Advani  _et al._ (2020);
Ba  _et al._ (2020); Barbier  _et al._ (2019); Bartlett  _et al._ (2020);
Belkin  _et al._ (2020); Bibas  _et al._ (2019); Deng  _et al._ (2020);
D’Ascoli  _et al._ (2020a, b); Dereziński  _et al._ (2020); Dhifallah and Lu
(2020); Gerace  _et al._ (2020); Hastie  _et al._ (2019); Jacot  _et al._
(2020); Kini and Thrampoulidis (2020); Lampinen and Ganguli (2019); Li  _et
al._ (2020); Liang and Rakhlin (2020); Liang  _et al._ (2020); Liao  _et al._
(2020); Lin and Dobriban (2021); Mitra (2019); Mei and Montanari (2021);
Muthukumar  _et al._ (2019); Nakkiran (2019); Xu and Hsu (2019); Yang  _et
al._ (2020); Rocks and Mehta (2020). Among, the most popular of these models
are linear regression (ridge regression without basis functions) and the
random nonlinear features model (a two-layer neural network with an arbitrary
nonlinear activation function where the top layer is trained and parameters
for the intermediate layer are chosen to be random but fixed) Rocks and Mehta
(2020). Here, we build upon this previous work by examining a random features
model for the special case of a linear activation function (i.e., the random
linear features model). Using the zero temperature cavity method, we derive
analytic expressions for the bias-variance decomposition and relate these
results to the eigenvalue spectrum of the Hessian matrix.

The random linear features has been treated analytically Advani  _et al._
(2020); Ba  _et al._ (2020); D’Ascoli  _et al._ (2020a); Lin and Dobriban
(2021); Yang  _et al._ (2020), with a subset of these studies attempting to
carrying out bias-variance decompositions Ba  _et al._ (2020); Lin and
Dobriban (2021); Yang  _et al._ (2020). However, these studies use non-
standard definitions of bias and variance that deviate from the traditional
textbook definitions Geman  _et al._ (1992); Bishop (2006). This choice of
definition can lead to qualitatively different and difficult to reconcile
results. For example, the authors of Ref. Ba  _et al._ , 2020 find that the
bias diverges at the interpolation threshold, while the authors of Ref. Lin
and Dobriban, 2021 find no such divergence (see Ref. Rocks and Mehta, 2020 for
an in-depth discussion). For this reason, here we utilize the standard
definitions and carry out the bias-variance decomposition in a manner
consistent with traditional definitions of these quantities in the
underparameterized regime, allowing us to identify which properties stem from
the model architecture versus random sampling of the data. In addition, we use
the zero-temperature cavity method to provide an alternative derivation of the
spectrum of the Hessian matrix of the random linear features model (i.e., the
spectrum of a Wishart product matrix calculated previously in Ref. Dupic and
Castillo, 2014), allowing us to directly relate the eigenvalues of the Hessian
to the double-descent phenomenon.

### Summary of Major Results

We briefly summarize our major results:

  * •

We derive analytic expressions for the test (generalization) error, training
error, bias, and variance for the random linear features model with a
nonlinear data distribution using the zero-temperature cavity method.

  * •

We find that the behavior of this model is characterized by three distinct
regimes: (i) an underparameterized regime with finite training error and large
bias, (ii) a second underparameterized regime with minimal, constant bias, and
(iii) an overparameterized, or interpolation, regime with zero training error.

  * •

We find that the three regimes are separated by three phase transitions with
two transitions to the interpolation regime, each characterized by a
divergence in the test error, and one transition between the large bias and
minimal bias underparameterized regimes. Importantly, we find that the
variance, but not the bias, diverges at the phase transition to the
interpolation regime.

  * •

We explain how each phase transition arises as a result of small nonzero
eigenvalues in the Hessian matrix and demonstrate how this phenomenon is
captured by susceptibilities.

  * •

We explain how the presence of linear features leads to an additional
interpolation phase transition not present in an analogous model with
nonlinear activation functions. We use random matrix theory to argue that the
underlying reason for this difference is that nonlinear basis functions
implicitly regularize small eigenvalues in the design matrix.

##  II Theoretical Setup

In this work, we focus on the supervised learning task of using relationships
learned from a training data set, consisting of labels and associated input
features, to accurately predict the labels of new data points from their input
features. Here, we closely follow the theoretical formalism previously
described in Ref. Rocks and Mehta, 2020.

###  II.1 Data Distribution (Teacher Model)

We consider data points $(y,\vec{\mathbf{x}})$, each consisting of a
continuous label $y$ paired with a set of $N_{f}$ continuous input features
$\vec{\mathbf{x}}$. We assume that the relationship between the input features
and labels (the data distribution or teacher model) can be expressed as

$\displaystyle y(\vec{\mathbf{x}})$
$\displaystyle=y^{*}(\vec{\mathbf{x}};\vec{\boldsymbol{\beta}})+\varepsilon$
(1)

where $\varepsilon$ is the label noise. The unknown function
$y^{*}(\vec{\mathbf{x}};\vec{\boldsymbol{\beta}})$ represents the “true”
labels and depends on a set of $N_{f}$ “ground truth” parameters
$\vec{\boldsymbol{\beta}}$, characterizing the correlations between the
features and labels. Here, we restrict ourselves to a teacher model of the
form

$y^{*}(\vec{\mathbf{x}};\vec{\boldsymbol{\beta}})=\frac{\sigma_{\beta}\sigma_{X}}{\expectationvalue*{f^{\prime}}}f\quantity(\frac{\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}}{\sigma_{X}\sigma_{\beta}}),$
(2)

where the function $f$ is an arbitrary nonlinear function and
${\expectationvalue*{f^{\prime}}=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\differential
he^{-\frac{h^{2}}{2}}f^{\prime}(h)}$ is a normalization constant chosen for
convenience with prime notation used to indicate a derivative. Note that Eq.
(2) reduces to a linear teacher model
$y^{*}(\vec{\mathbf{x}})=\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}$ when
$f(h)=h$.

We draw the input features for each data point independently and identically
from a normal distribution with zero mean and variance $\sigma_{X}^{2}/N_{f}$.
We consider ground truth parameters $\vec{\boldsymbol{\beta}}$ and label noise
$\varepsilon$ that are drawn independently from normal distributions with zero
mean and variances $\sigma_{\beta}^{2}$ and $\sigma_{\varepsilon}^{2}$,
respectively. Furthermore, we assume the labels are centered so that $f$ has
zero mean with respect to its argument.

###  II.2 Model Architectures (Student Models)

We consider a student model of the form

$\displaystyle\hat{y}(\vec{\mathbf{x}})$
$\displaystyle=\vec{\mathbf{z}}(\vec{\mathbf{x}})\cdot\hat{\mathbf{w}},$ (3)

where $\hat{\mathbf{w}}$ is a vector of $N_{p}$ fit parameters. For the random
linear features model, the vector of ‘hidden” features
$\vec{\mathbf{z}}(\vec{\mathbf{x}})$ takes the form

$\displaystyle\vec{\mathbf{z}}(\vec{\mathbf{x}})$
$\displaystyle=W^{T}\vec{\mathbf{x}},$ (4)

where $W$ is a random transformation matrix of size $N_{f}\times N_{p}$, whose
elements are drawn independently from a normal distribution with zero mean and
variance $\sigma_{W}^{2}/N_{p}$.

###  II.3 Fitting Procedure

We train each model on a training data set consisting of $M$ data points,
${\mathcal{D}=\\{(y_{a},\vec{\mathbf{x}}_{a})\\}_{a=1}^{M}}$. For convenience,
we organize the vectors of input features in the training set into an
observation matrix $X$ of size $M\times N_{f}$ and define the length-$M$
vectors of training labels $\vec{\mathbf{y}}$, training label noise
$\vec{\boldsymbol{\varepsilon}}$, and label predictions for the training set
$\hat{\mathbf{y}}$. We also organize the vectors of hidden features evaluated
on the input features of the training set,
$\\{\vec{\mathbf{z}}(\vec{\mathbf{x}}_{a})\\}_{a=1}^{M}$, into the rows of a
hidden feature matrix $Z$ of size $M\times N_{p}$.

Given a set of training data $\mathcal{D}$, we solve for the optimal values of
the fit parameters $\hat{\mathbf{w}}$ by minimizing the standard ridge
regression loss function composed of the mean squared label error with $L_{2}$
regularization,

$\displaystyle L(\hat{\mathbf{w}};\mathcal{D})$
$\displaystyle=\frac{1}{2}\norm*{\Delta\vec{\mathbf{y}}}^{2}+\frac{\lambda}{2}\norm*{\hat{\mathbf{w}}}^{2},$
(5)

where the notation $\norm*{\cdot}$ indicates an $L_{2}$ norm,
${\Delta\vec{\mathbf{y}}=\vec{\mathbf{y}}-\hat{\mathbf{y}}}$ is the vector of
residual label errors for the training set, and $\lambda$ is the
regularization parameter. The exact solution for the fit parameters resulting
from this loss function are

$\displaystyle\hat{\mathbf{w}}$ $\displaystyle=\quantity[\lambda
I_{N_{p}}+Z^{T}Z]^{-1}Z^{T}\vec{\mathbf{y}}.$ (6)

We will often work in the “ridge-less limit” where we take the limit
$\lambda\rightarrow 0$. In this limit, we refer to the matrix $Z^{T}Z$ in the
above expression as the Hessian matrix.

###  II.4 Model Evaluation

To evaluate a model’s prediction accuracy, we measure the training and test
(generalization) errors. We define the training error as the mean squared
residual label error of the training data,

$\displaystyle\mathcal{E}_{\mathrm{train}}$
$\displaystyle=\frac{1}{M}\norm*{\Delta\vec{\mathbf{y}}}^{2}.$ (7)

We define the interpolation threshold as the model complexity at which the
training error becomes exactly zero (in the ridge-less limit). Analogously, we
define the test error as the mean squared error evaluated on a test data set,
${\mathcal{D}^{\prime}=\\{(y_{a}^{\prime},\vec{\mathbf{x}}_{a}^{\prime})\\}_{a=1}^{M^{\prime}}}$,
composed of $M^{\prime}$ new data points drawn independently from the same
data distribution as the training set,

$\displaystyle\mathcal{E}_{\mathrm{test}}$
$\displaystyle=\frac{1}{M^{\prime}}\norm*{\Delta\vec{\mathbf{y}}^{\prime}}^{2},$
(8)

where
$\Delta\vec{\mathbf{y}}^{\prime}=\vec{\mathbf{y}}^{\prime}-\hat{\mathbf{y}}^{\prime}$
is a length-$M^{\prime}$ vector of residual label errors between the vector of
test labels $\vec{\mathbf{y}}^{\prime}$ and their predicted values
$\hat{\mathbf{y}}^{\prime}$. Furthermore, we define the ensemble-averaged
training and test errors, $\expectationvalue{\mathcal{E}_{\mathrm{train}}}$
and $\expectationvalue{\mathcal{E}_{\mathrm{test}}}$, respectively, by taking
averages of the above definitions with respect to all sources of randomness
(e.g., $X$, $\vec{\boldsymbol{\varepsilon}}$, $\vec{\boldsymbol{\beta}}$,
etc.).

###  II.5 Bias-Variance Decomposition

The bias-variance decomposition separates test error into components stemming
from three distinct sources: bias, variance, and noise. Here, we utilize the
standard definitions of bias and variance Geman  _et al._ (1992); Bishop
(2006),

$\displaystyle\mathrm{Bias}\quantity[\hat{y}(\vec{\mathbf{x}})]$
$\displaystyle=\mathrm{E}_{\mathcal{D}}\quantity[\hat{y}(\vec{\mathbf{x}})]-y^{*}(\vec{\mathbf{x}})$
(9) $\displaystyle\mathrm{Var}\quantity[\hat{y}(\vec{\mathbf{x}})]$
$\displaystyle=\mathrm{E}_{\mathcal{D}}\quantity[\hat{y}^{2}(\vec{\mathbf{x}})]-\mathrm{E}_{\mathcal{D}}\quantity[\hat{y}(\vec{\mathbf{x}})]^{2},$
(10)

where the subscript $\mathcal{D}$ denotes the sampling average with respect to
the training set (i.e., with respect to the input features $X$ and label noise
$\vec{\boldsymbol{\varepsilon}}$, but not the ground truth parameters
$\vec{\boldsymbol{\beta}}$).

In order to incorporate other sources of randomness (e.g.,
$\vec{\boldsymbol{\beta}}$ and $W$), we define the more general ensemble-
averaged squared bias and variance, respectively, as

$\displaystyle\expectationvalue*{\mathrm{Bias}^{2}[\hat{y}]}$
$\displaystyle=\mathrm{E}_{\vec{\boldsymbol{\beta}},W,\vec{\mathbf{x}}}\quantity[\mathrm{Bias}[\hat{y}(\vec{\mathbf{x}})]^{2}]$
(11) $\displaystyle\expectationvalue*{\mathrm{Var}[\hat{y}]}$
$\displaystyle=\mathrm{E}_{\vec{\boldsymbol{\beta}},W,\vec{\mathbf{x}}}\quantity[\mathrm{Var}[\hat{y}(\vec{\mathbf{x}})]].$
(12)

Using these definitions, we define the ensemble-averaged bias-variance
decomposition of the test error,

$\displaystyle\expectationvalue{\mathcal{E}_{\mathrm{test}}}$
$\displaystyle=\expectationvalue*{\mathrm{Bias}^{2}[\hat{y}]}+\expectationvalue*{\mathrm{Var}[\hat{y}]}+\sigma_{\varepsilon}^{2},$
(13)

which we utilize throughout this work.

###  II.6 Derivation of Closed-Form Solutions

Following the derivations in Ref. Rocks and Mehta, 2020, we utilize the zero-
temperature cavity method to derive closed-form expressions for the training
error, test error, bias, and variance. In this derivation, we work in the
thermodynamic limit, where $N_{f},M,N_{p}\rightarrow\infty$, but their ratios,
$\alpha_{f}=N_{f}/M$ and $\alpha_{p}=N_{p}/M$, remain finite. Our results are
exact in this limit. Furthermore, we utilize the procedure described in Ref.
Cui  _et al._ , 2020 to reproduce the closed-form solution for the eigenvalues
spectrum of the Hessian matrix for this model. We refer the reader to the
Appendix for further details on these calculations.

Figure 1: Random Linear Features Model (Two-layer Linear Neural Network).
Analytic solutions plotted as a function of $\alpha_{p}=N_{p}/M$ with fixed
$\alpha_{f}=N_{f}/M$ for (a) less input features than training data points
($\alpha_{f}=1/2$) and (b) more input features than training data points
($\alpha_{f}=4$). Shown are the ensemble-averaged training error (blue
squares), test error (black circles), squared bias (green triangles), and
variance (red diamonds). Analytic solutions are indicated as dashed lines with
numerical results shown as points. In (a), a black dashed vertical line marks
the boundary between the large bias and minimal bias underparameterized
regimes at $\alpha_{p}=\alpha_{f}$, while in (b), a similar line marks the
boundary between the under- and overparameterized regimes at $\alpha_{p}=1$.
Analytic solutions as a function of $\alpha_{p}$ and $\alpha_{f}$ are also
shown for the the ensemble-averaged (c) training error, (d) test error, (e)
squared bias, and (f) variance. Results are shown for a linear teacher model
$y(\vec{\mathbf{x}})=\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}+\varepsilon$,
a signal-to-noise ratio of
$\sigma_{\beta}^{2}\sigma_{X}^{2}/\sigma_{\varepsilon}^{2}=10$, and have been
scaled by the variance of the training set labels
$\sigma_{y}^{2}=\sigma_{\beta}^{2}\sigma_{X}^{2}+\varepsilon^{2}$. In each
panel, black dashed lines indicate boundaries between different regimes of the
solutions depending on which is the smallest of the quantities $M$, $N_{f}$,
or $N_{p}$. See Appendix for additional numerical details.

##  III Analytic Expressions

We find that the closed-form solutions for the random linear features model
are characterized by three distinct regimes, each defined by which of the
following three quantities is the smallest: the number of input features
$N_{f}$, the number of fit parameters (hidden features) $N_{p}$, or the size
of the training set $M$. In terms of $\alpha_{f}=N_{f}/M$ and
$\alpha_{p}=N_{p}/M$, the expressions for ensemble-averaged the training
error, test error, bias, and variance are

$\displaystyle\expectationvalue{\mathcal{E}_{\mathrm{train}}}$
$\displaystyle=\left\\{\begin{array}[]{c}(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})(1-\alpha_{f})\\\
\sigma_{\beta}^{2}\sigma_{X}^{2}\frac{(1-\alpha_{p})(\alpha_{f}-\alpha_{p})}{\alpha_{f}}+(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})(1-\alpha_{p})\\\ 0\end{array}\right.$
$\displaystyle\begin{array}[]{l}\mbox{\quad if\quad}N_{f}<N_{p},M\\\
\mbox{\quad if\quad}N_{p}<N_{f},M\\\ \mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}$ (20)
$\displaystyle\expectationvalue{\mathcal{E}_{\mathrm{test}}}$
$\displaystyle=\left\\{\begin{array}[]{c}(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{1}{(1-\alpha_{f})}\\\
\sigma_{\beta}^{2}\sigma_{X}^{2}\frac{(\alpha_{f}-\alpha_{p})}{\alpha_{f}(1-\alpha_{p})}+(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{1}{(1-\alpha_{p})}\\\
\sigma_{\beta}^{2}\sigma_{X}^{2}\frac{\alpha_{p}(\alpha_{f}-1)}{\alpha_{f}(\alpha_{p}-1)}+(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{(\alpha_{f}\alpha_{p}-1)}{(\alpha_{f}-1)(\alpha_{p}-1)}\end{array}\right.$
$\displaystyle\begin{array}[]{l}\mbox{\quad if\quad}N_{f}<N_{p},M\\\
\mbox{\quad if\quad}N_{p}<N_{f},M\\\ \mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}$ (27)
$\displaystyle\expectationvalue*{\mathrm{Bias}^{2}[\hat{y}]}$
$\displaystyle=\left\\{\begin{array}[]{c}\sigma_{\delta y^{*}}^{2}\\\
\sigma_{\beta}^{2}\sigma_{X}^{2}\frac{(\alpha_{f}-\alpha_{p})}{\alpha_{f}}+\sigma_{\delta
y^{*}}^{2}\\\
\sigma_{\beta}^{2}\sigma_{X}^{2}\frac{\alpha_{p}(\alpha_{f}-1)^{2}}{\alpha_{f}(\alpha_{f}\alpha_{p}-1)}+\sigma_{\delta
y^{*}}^{2}\end{array}\right.$ $\displaystyle\begin{array}[]{l}\mbox{\quad
if\quad}N_{f}<N_{p},M\\\ \mbox{\quad if\quad}N_{p}<N_{f},M\\\ \mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}$ (34)
$\displaystyle\expectationvalue*{\mathrm{Var}[\hat{y}]}$
$\displaystyle=\left\\{\begin{array}[]{c}(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{\alpha_{f}}{(1-\alpha_{f})}\\\
\sigma_{\beta}^{2}\sigma_{X}^{2}\frac{\alpha_{p}(\alpha_{f}-\alpha_{p})}{\alpha_{f}(1-\alpha_{p})}+(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{\alpha_{p}}{(1-\alpha_{p})}\\\
\sigma_{\beta}^{2}\sigma_{X}^{2}\frac{\alpha_{p}(\alpha_{f}-1)(\alpha_{f}-1+\alpha_{p}-1)}{\alpha_{f}(\alpha_{p}-1)(\alpha_{f}\alpha_{p}-1)}+(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{(\alpha_{f}-1+\alpha_{p}-1)}{(\alpha_{f}-1)(\alpha_{p}-1)}\end{array}\right.$
$\displaystyle\begin{array}[]{l}\mbox{\quad if\quad}N_{f}<N_{p},M\\\
\mbox{\quad if\quad}N_{p}<N_{f},M\\\ \mbox{\quad
if\quad}M<N_{f},N_{p},\end{array}$ (41)

where we have taken the limit $\lambda\rightarrow 0$ (with leading order terms
of order $\lambda^{2}$ reported in the Appendix for quantities reported here
as zero). The quantity $\sigma_{\delta y^{*}}^{2}$ is the statistical variance
of the nonlinear components of the true labels, measured via their deviation
from a linear teacher model, $\sigma_{\delta
y^{*}}^{2}=\mathrm{E}_{\vec{\mathbf{x}}}[(y^{*}(\vec{\mathbf{x}})-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}]$.
In the thermodynamic limit, we find that this quantity evaluates to

$\begin{gathered}\sigma_{\delta
y^{*}}^{2}=\sigma_{\beta}^{2}\sigma_{X}^{2}\Delta f,\qquad\Delta
f=\frac{\expectationvalue*{f^{2}}-\expectationvalue*{f^{\prime}}^{2}}{\expectationvalue*{f^{\prime}}^{2}}\\\
\expectationvalue*{f^{2}}=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}\differential
he^{-\frac{1}{2}h^{2}}f^{2}(h)\\\
\expectationvalue*{f^{\prime}}=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}\differential
he^{-\frac{1}{2}h^{2}}f^{\prime}(h).\end{gathered}$ (42)

In Figs. 1(a) and (b), we plot the training error, test error, bias, and
variance as a function of $\alpha_{p}=N_{p}/M$ for fixed $\alpha_{f}=N_{f}/M$
for the two cases $\alpha_{f}<1$ and $\alpha_{f}>1$, respectively. To gain a
better grasp of the full set of solutions, we also plot all quantities in Eqs.
(20)-(41) as a function of both $\alpha_{p}$ and $\alpha_{f}$ in Figs.
1(c)-(f). All solutions are shown for a linear teacher model
${y^{*}(\vec{\mathbf{x}})=\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}}$
(${\sigma_{\delta y^{*}}^{2}=0}$).

The three regimes we observe in the closed-form solutions are separated by
three distinct phase transitions. Examining the training error in Fig. 1(c),
we find that it goes to zero at two of the transitions, $\alpha_{p}=1$ with
$\alpha_{f}\geq 1$ and $\alpha_{f}=1$ with $\alpha_{p}\geq 1$, giving rise to
an interpolation boundary. On one side of this boundary, where $\alpha_{p}<1$
or $\alpha_{f}<1$ (there are less data points $M$ than fit parameters $N_{p}$
or input features $N_{f}$), the model is underparameterized, while beyond this
boundary the model is overparameterized and in the “interpolation” regime. We
note that the interpolation transition for this model is markedly different
from that of the random nonlinear features model (nonlinear activation
function) where the interpolation threshold occurs at $\alpha_{p}=1$
independently of $\alpha_{f}$ (see Sec. V and Ref. Rocks and Mehta, 2020).

We find that the test error and variance diverge along the entire
interpolation boundary, while the bias remains finite. This is in contrast
with previous studies which employed non-standard definitions of bias and
variance and found that both variance and bias diverge at the interpolation
threshold Ba  _et al._ (2020). Examining Figs. 1(a) and (b), we find that the
test error (and similarly, the variance) exhibits very different behavior as a
function of $\alpha_{p}$, depending on whether $\alpha_{f}<1$ (less input
features than training data points $N_{f}<M$) or $\alpha_{f}>1$ (more input
features than training data points $N_{f}>M$). When $\alpha_{f}>1$, the test
error diverges at $\alpha_{p}=1$ and decreases monotonically in the
overparameterized regime. In contrast, when $\alpha_{f}<1$, the test error
monotonically decreases to a small, constant value at
$\alpha_{p}\geq\alpha_{f}$. Although this model does not display the full
canonical double-descent behavior in either case, the test error in the
overparameterized regime is always at least as small as – if not smaller than
– that of the underparameterized regime for fixed $\alpha_{f}$.

Examining the bias in Fig 1(c) reveals that there is an additional phase
transition in the underparameterized regime located at the boundary
$\alpha_{f}=\alpha_{p}$ for $\alpha_{p}\leq 1$ and $\alpha_{f}\leq 1$ (i.e.,
when the number of input features $N_{f}$ equals the number of hidden features
$N_{p}$, with both $N_{f}$ and $N_{p}$ less than the number of data points
$M$). This transition divides the non-interpolation solutions into two pieces.
When $\alpha_{p}<\alpha_{f}$ (less fit parameters than input features
$N_{p}<N_{f}$), the model exhibits a large bias because there are not enough
fit parameters (or hidden features) to fully express the input features in the
data [see Eq. (34)]. In contrast, when $\alpha_{p}>\alpha_{f}$ (more fit
parameters than input features $N_{p}>N_{f}$), the only contribution to the
bias is a small constant $\sigma_{\delta y^{*}}^{2}$ stemming from the
nonlinear components of the labels. For the special case of a linear teacher
model shown in the figures, the bias is identically zero in this regime.

Interestingly, we also observe that $\sigma_{\delta y^{*}}^{2}$ appears as an
additive component to the label noise $\sigma_{\varepsilon}^{2}$ in the
training error, test error, and variance, indicating that the model interprets
the nonlinear components of the labels as effective noise Rocks and Mehta
(2020).

Figure 2:  Susceptibilities and Eigenvalue Spectra (a)-(b) Analytic solutions
for the susceptibility $\nu$ which measures the sensitivity of the fit
parameters to small perturbations in the gradient. In the small $\lambda$
limit, ${\nu\approx\lambda^{-1}\nu_{-1}+\nu_{0}}$. (a) The coefficient
$\nu_{-1}$ counts the fraction of unutilized fit parameters, or the fraction
of parameters beyond that needed to attain minimal training error. (b) The
coefficient $\nu_{0}$ diverges at each phase transition when $Z^{T}Z$ has a
small eigenvalue. (c) Analytic solution for the minimum nonzero eigenvalue
$\sigma_{\min}^{2}$ of the Hessian matrix $Z^{T}Z$. (i)-(ix) Examples of the
full eigenvalue spectrum are shown for each of the corresponding points in
(c). See Appendix for additional simulation details.

##  IV Phase transitions, susceptiblities, and eigenvalue spectra

In the previous section, we found that the analytic solutions for the random
linear features model are characterized by three distinct regimes separated by
three different phase transitions. As a natural byproduct of our cavity
derivations, we find that each of these phase transitions is marked by a
diverging susceptibility. In particular, setting the gradient of the loss
function in Eq. (5) equal to a small nonzero field $\vec{\boldsymbol{\eta}}$,
such that $\partialderivative*{L}{\hat{\mathbf{w}}}=\vec{\boldsymbol{\eta}}$,
we define the susceptibility matrix
$\partialderivative*{\hat{\mathbf{w}}}{\vec{\boldsymbol{\eta}}}$. This
quantity measures the sensitivity of the fit parameters to small perturbation
in the gradient and can be shown to be equivalent to the inverse Hessian of
the loss function. Taking the trace of this matrix, we define the scalar
susceptibility

$\nu=\frac{1}{N_{p}}\Tr\partialderivative{\hat{\mathbf{w}}}{\vec{\boldsymbol{\eta}}}=\frac{1}{N_{p}}\Tr\quantity[\lambda
I_{N_{p}}+Z^{T}Z]^{-1}.$ (43)

In the small $\lambda$ limit, we make the approximation
${\nu\approx\lambda^{-1}\nu_{-1}+\nu_{0}}$. In exact matrix form, we find that
the two coefficients are

$\displaystyle\nu_{-1}=1-\frac{1}{N_{p}}\rank(Z^{T}Z),\quad\nu_{0}=\frac{1}{N_{p}}\Tr[Z^{T}Z]^{+},$
(44)

where + denotes a Moore-Penrose pseudoinverse.

In Figs. 2(a) and (b), we plot the analytic closed-form expressions for these
two quantities in the thermodynamic limit as a function of $\alpha_{f}$ and
$\alpha_{p}$ (see Appendix for expressions). We find that the first
coefficient $\nu_{-1}$ counts the fraction of fit parameters that go beyond
the minimum needed to attain minimal training error. In contrast, the second
coefficient $\nu_{0}$ diverges along each phase boundary. Based on the exact
matrix form of $\nu_{0}$ in Eq. (44), we infer that these divergences can be
attributed to small eigenvalues in the Hessian matrix $Z^{T}Z$.

To illustrate this connection between the eigenvalues of the Hessian and the
susceptibility $\nu$, we note that for this problem, the inverse Hessian is
equivalent to the Green’s function and can be used to extract the eigenvalue
spectrum Cui  _et al._ (2020) (see Appendix for derivation). In Fig. 2(c), we
show the analytic solution for the minimum nonzero eigenvalue
$\sigma_{\min}^{2}$ of $Z^{T}Z$. Consistent with $\nu_{0}$, we find that
$\sigma_{\min}^{2}$ goes to zero along each phase transition. In Figs.
2(i)-(ix), we also plot the eigenvalue distributions for the points indicated
in Fig. 2(c). While $\nu_{0}$ captures the distribution of nonzero
eigenvalues, $\nu_{-1}$ captures the weight of the delta function at zero in
the overparameterized and minimal bias regimes. In each regime and along each
phase boundary, these distributions are qualitatively similar to the
Marchenko-Pastur distribution Marčenko and Pastur (1967). Along each phase
transition, the gap in the distribution goes to zero, while the gap is finite
away from each boundary. The presence of this eigenvalue gap was previously
shown to be the root cause of the decrease in variance in the
overparameterized regime Rocks and Mehta (2020).

To understand the source of these small eigenvalue gaps, we note that the
types of random matrices we consider in this work typically exhibit
infinitesimally small eigenvalues in the thermodynamic limit if they contain
an equal number of rows and columns. This fact allows us to identify which
matrix is the root cause of each transition. Since $Z$ is a product of $X$ and
$W$, this phenomenon arises in two forms. First, $Z$ can exhibit a small
eigenvalue if either $X$ or $W$ is square and the expression of its input
feature space is not limited by its product with the other matrix (e.g., if
$X$ is square and $N_{p}\geq N_{f}=M$ or $W$ is square and $M\geq
N_{f}=N_{p}$). This behavior explains the interpolation transition at
$\alpha_{f}=1$ which arises due to small eigenvalues in $X$, but does not
extend below $\alpha_{p}=1$ when the rank of $W$ becomes too low to preserve
every direction in the space of input features encoded in $X$. Similarly, the
minimal bias transition at $\alpha_{f}=\alpha_{p}$ occurs due to small
eigenvalues in $W$, disappearing above $\alpha_{f}=1$ when the rank of $X$ is
too low to fully express the input feature space of $W$. Second, $Z$ can
exhibit a small eigenvalue if it is square and full rank, giving rise to the
transition at $\alpha_{p}=1$, but only when $\alpha_{f}\geq 1$.

Finally, we observe that the test error only diverges at the two phase
transitions to the interpolation regime, but not at the minimal bias
transition, despite $\sigma_{\min}^{2}$ going to zero in all three cases. This
lack of divergence is explained by the fact that the minimal bias transition
arises due to small eigenvalues in $W$, which is used to transform both the
training data and the test data. Since both data sets are transformed in the
same way, predictions by the model for the test set based on the training set
will not be limited by small eigenvalues in $W$, and the test error will not
diverge when $W$ is square ($N_{f}=N_{p}$).

Figure 3:  Comparison of the Random Linear Features Model to Linear Regression
(No Basis Functions) and the Random Nonlinear Features Model. (a) The training
error, test error, bias, and variance as a function of $\alpha_{p}$ (or
equivalently, $\alpha_{f}$) for linear regression (no basis functions) , and
(b) the same quantities for the random linear features model for the special
case $\alpha_{f}=\alpha_{p}$. The training error, test error, bias, and
variance as a function of $\alpha_{p}$ and $\alpha_{f}$ for (c)-(f) the random
nonlinear features model with ReLU activation, $\varphi(h)=\max(0,h)$, and
(g)-(j) the random linear features model. In all panels, black dashed lines
indicate phase transitions. Results for linear regression and the random
nonlinear features model are reproduced from Ref. Rocks and Mehta, 2020.

##  V Comparison to Linear Regression and the Random Nonlinear Features Model

One of the more surprising results of our analysis is that the phase diagram
for the random _linear_ features model is qualitatively different from the
random _nonlinear_ features model. On other hand, we find that the random
linear feature model is qualitatively similar to ordinary ridge-less
regression when the number of hidden features matches the input features. To
better understand the similarities and differences, we have reproduced the
phase diagrams for all three models in Fig. 3 (see Ref. Rocks and Mehta (2020)
for a detailed analysis of ridge regression and the random nonlinear features
model).

First, we compare the random linear features model to linear regression in
which the number of hidden features matches the input features,

$\vec{\mathbf{z}}(\vec{\mathbf{x}})=\vec{\mathbf{x}}.$ (45)

Fig. 3(a) shows the training error, test error, bias, and variance for linear
regression. Since linear regression lacks basis functions, in Fig. 3(b), we
plot the same quantities for the random linear features model for the special
case where the number of input features equals the number of hidden features
($\alpha_{f}=\alpha_{p}$).

We observe that along this cut of the phase diagram, the random linear
features model behaves qualitatively similar to linear regression, with
variance first increasing as one approaches the interpolation threshold
($\alpha_{p}=1$) and then decreasing monotonically beyond the threshold.
Meanwhile, the bias is zero below the interpolation threshold and then
increases once one crosses the interpolation threshold. As discussed in detail
in Ref. Rocks and Mehta (2020), the underlying reason for the increase in bias
for $\alpha_{p}>1$ is that in this regime, the model does not have enough
training data points to sample the entire input feature space. Therefore, any
predictions made about these unsampled directions represent implicit
assumptions of the model.

Next, we compare to the random nonlinear features model in which the hidden
features take the form

$\vec{\mathbf{z}}(\vec{\mathbf{x}})=\frac{1}{\expectationvalue*{\varphi^{\prime}}}\frac{\sigma_{W}\sigma_{X}}{\sqrt{N_{p}}}\varphi\quantity(\frac{\sqrt{N_{p}}}{\sigma_{W}\sigma_{X}}W^{T}\vec{\mathbf{x}})$
(46)

where $\varphi$ is a nonlinear activation function that acts separately on
each element of its input and
${\expectationvalue*{\varphi^{\prime}}=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}dhe^{-\frac{h^{2}}{2}}\varphi^{\prime}(h)}$
is a normalization constant. Figs. 2(c)-(f) show the training error, test
error, bias, and variance for this model for the case of ReLU activation,
$\varphi(h)=\max(0,h)$, as a function of both $\alpha_{p}$ and $\alpha_{f}$,
while Figs. 2(g)-(j) depict the same for the random linear features model.

We observe that while the interpolation transition boundary at $\alpha_{p}=1$
for $\alpha_{f}\geq 1$ remains the same, the addition of a nonlinear
activation function suppresses the interpolation transition at $\alpha_{f}=1$
for $\alpha_{p}>1$, along with the transition to a minimal bias regime at
$\alpha_{f}=\alpha_{p}$ for $\alpha_{p}<1$. At the same time, the
interpolation transition at $\alpha_{p}=1$ is extended to all values of
$\alpha_{f}$.

The two changes to the shape of the interpolation boundary can be attributed
to the behavior of small eigenvalues in the Hessian. Upon the introduction of
a nonlinear activation, the interpolation transition at $\alpha_{p}=1$ for
$\alpha_{f}<1$ arises due to the creation of new small eigenvalues in $Z$. The
nonlinear transformation promotes $Z$ to full rank when it is square
($N_{p}=M$), even if the product $XW$ is not full rank. As a result, $Z$
exhibits small eigenvalues at this transition whether or not $XW$ has a small
eigenvalue, translating to small eigenvalues in the Hessian and a divergence
in the test error. In contrast, our random matrix theory analysis suggests
that the presence of a nonlinear activation function suppresses the transition
at $\alpha_{f}=1$ with $\alpha_{p}>1$ by serving as an implicit regularizer of
small eigenvalues in the Hessian (this behavior was previously observed in
Ref. D’Ascoli  _et al._ , 2020a). In particular, the use of nonlinear basis
functions masks divergences arising from small eigenvalues that arise when the
design matrix $X$ is square ($N_{f}=M$).

Finally, to account for the removal of the minimal bias transition when
$\alpha_{f}=\alpha_{p}$ with $\alpha_{p}<1$, we note that the training data is
generated using a teacher model that depends directly on the input features,
while the nonlinear model first applies a nonlinear transformation. This
nonlinear basis masks properties of the underlying input feature space like
its dimension, introducing additional bias. Therefore, the bias does not
approach a minimal value at $\alpha_{p}=\alpha_{f}$, even if there are in
principle enough hidden features to fully encode the full space of input
features. Instead, we observe that the bias only reaches a minimum in the
limit $\alpha_{p}\rightarrow\infty$ for fixed $\alpha_{f}$.

##  VI Conclusions

A central question in machine learning is understanding why complicated models
with many more parameters than training data points can still make accurate
predictions. Here, we have tackled this problem by analyzing one of the
simplest examples of non-trivial supervised learning: regression with random
linear features. Despite the simplicity of the model, it exhibits remarkably
rich behavior with multiple phase transitions.

We found that the phase diagram of the model has three distinct phases: (i) an
underparameterized regime with finite training error and large bias, (ii) a
second underparameterized regime with minimal bias, and (iii) an
overparameterized, or interpolation, regime with zero training error. We also
showed that at the transition to the interpolation regime, the variance but
not the bias diverges. For this reason, while the classical bias-variance
trade-off captures much of the behavior of the model in the underparameterized
regime, it fails to describe the interpolation regime where the variance
decreases with increasing model complexity.

We showed that the divergence of the variance is due to the presence of small
eigenvalues in the Hessian matrix. This is consistent with the general picture
advocated in Refs. Advani  _et al._ , 2020; Rocks and Mehta, 2020 that large
test errors are associated with the closing of a spectral gap near the
interpolation transition. On both sides of the transition, when the spectral
gap is large, it is easy to distinguish noise from poorly sampled directions
in feature space. However, when the gap closes this is no longer possible,
accounting for the large variance.

Our work suggests that many of the fundamental features of double-descent can
be understood even when considering simple convex models. An important
question is how to generalize the intuitions developed here to more complex
settings. In contrast with the random linear features model, modern deep
learning methods are often non-convex and attempt to learn meaningful features
directly from data. In the future, it will be interesting to see how this
changes the understanding of the bias-variance trade-off developed here
Chaudhari  _et al._ (2019); Baldassi  _et al._ (2020); Pittorino  _et al._
(2021).

## Acknowledgments

This work was supported by NIH NIGMS grant 1R35GM119461 and a Simons
Investigator in the Mathematical Modeling of Living Systems (MMLS) award to
PM. The authors also acknowledge support from the Shared Computing Cluster
administered by Boston University Research Computing Services.

## References

  * Geman  _et al._ (1992) Stuart Geman, Elie Bienenstock, and René Doursat, “Neural Networks and the Bias/Variance Dilemma,” Neural Computation 4, 1–58 (1992).
  * Bishop (2006) Christopher M. Bishop,  _Pattern Recognition and Machine Learning_ (Springer, 2006).
  * Mehta  _et al._ (2019) Pankaj Mehta, Marin Bukov, Ching Hao Wang, Alexandre G.R. Day, Clint Richardson, Charles K. Fisher, and David J. Schwab, “A high-bias, low-variance introduction to Machine Learning for physicists,” Physics Reports 810, 1–124 (2019).
  * Zhang  _et al._ (2017) Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals, “Understanding Deep Learning Requires Re-thinking Generalization,” International Conference on Learning Representations (ICLR) (2017).
  * Belkin  _et al._ (2019) Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal, “Reconciling modern machine-learning practice and the classical bias–variance trade-off,” Proceedings of the National Academy of Sciences 116, 15849–15854 (2019).
  * Loog  _et al._ (2020) Marco Loog, Tom Viering, Alexander Mey, Jesse H. Krijthe, and David M. J. Tax, “A brief prehistory of double descent,” Proceedings of the National Academy of Sciences 117, 10625–10626 (2020).
  * Adlam and Pennington (2020) Ben Adlam and Jeffrey Pennington, “Understanding double descent requires a fine-grained bias-variance decomposition,” Advances in Neural Information Processing Systems (NeurIPS) 33, 11022–11032 (2020).
  * Advani  _et al._ (2020) Madhu S. Advani, Andrew M. Saxe, and Haim Sompolinsky, “High-dimensional dynamics of generalization error in neural networks,” Neural Networks 132, 428–446 (2020).
  * Ba  _et al._ (2020) Jimmy Ba, Murat Erdogdu, Taiji Suzuki, Denny Wu, and Tianzong Zhang, “Generalization of Two-layer Neural Networks: An Asymptotic Viewpoint,” International Conference on Learning Representations (ICLR) (2020).
  * Barbier  _et al._ (2019) Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, and Lenka Zdeborová, “Optimal errors and phase transitions in high-dimensional generalized linear models,” Proceedings of the National Academy of Sciences 116, 5451–5460 (2019).
  * Bartlett  _et al._ (2020) Peter L. Bartlett, Philip M. Long, Gábor Lugosi, and Alexander Tsigler, “Benign overfitting in linear regression,” Proceedings of the National Academy of Sciences 117, 30063–30070 (2020).
  * Belkin  _et al._ (2020) Mikhail Belkin, Daniel Hsu, and Ji Xu, “Two Models of Double Descent for Weak Features,” SIAM Journal on Mathematics of Data Science 2, 1167–1180 (2020).
  * Bibas  _et al._ (2019) Koby Bibas, Yaniv Fogel, and Meir Feder, “A New Look at an Old Problem: A Universal Learning Approach to Linear Regression,” IEEE International Symposium on Information Theory (ISIT) , 2304–2308 (2019).
  * Deng  _et al._ (2020) Zeyu Deng, Abla Kammoun, and Christos Thrampoulidis, “A Model of Double Descent for High-Dimensional Logistic Regression,” IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) , 4267–4271 (2020).
  * D’Ascoli  _et al._ (2020a) Stéphane D’Ascoli, Levent Sagun, and Giulio Biroli, “Triple descent and the two kinds of overfitting: Where & why do they appear?” Advances in Neural Information Processing Systems (NeurIPS) 33 (2020a).
  * D’Ascoli  _et al._ (2020b) Stéphane D’Ascoli, Maria Refinetti, Giulio Biroli, and Florent Krzakala, “Double Trouble in Double Descent: Bias and Variance(s) in the Lazy Regime,” Proceedings of the 37th International Conference on Machine Learning (ICML) PMLR, 119, 2280–2290 (2020b).
  * Dereziński  _et al._ (2020) Michał Dereziński, Feynman Liang, and Michael W. Mahoney, “Exact expressions for double descent and implicit regularization via surrogate random design,” Advances in Neural Information Processing Systems (NeurIPS) 33, 5152–5164 (2020).
  * Dhifallah and Lu (2020) Oussama Dhifallah and Yue M. Lu, “A Precise Performance Analysis of Learning with Random Features,”  (2020), arXiv:2008.11904 .
  * Gerace  _et al._ (2020) Federica Gerace, Bruno Loureiro, Florent Krzakala, Marc Mézard, and Lenka Zdeborová, “Generalisation error in learning with random features and the hidden manifold model,” Proceedings of the 37th International Conference on Machine Learning (ICML) PMLR, 119, 3452–3462 (2020).
  * Hastie  _et al._ (2019) Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J. Tibshirani, “Surprises in High-Dimensional Ridgeless Least Squares Interpolation,”  (2019), arXiv:1903.08560 .
  * Jacot  _et al._ (2020) Arthur Jacot, Berfin Şimşek, Francesco Spadaro, Clément Hongler, and Franck Gabriel, “Implicit regularization of random feature models,” Proceedings of the 37th International Conference on Machine Learning (ICML) PMLR, 119, 4631–4640 (2020).
  * Kini and Thrampoulidis (2020) Ganesh Ramachandra Kini and Christos Thrampoulidis, “Analytic Study of Double Descent in Binary Classification: The Impact of Loss,” IEEE International Symposium on Information Theory (ISIT) , 2527–2532 (2020).
  * Lampinen and Ganguli (2019) Andrew K. Lampinen and Surya Ganguli, “An analytic theory of generalization dynamics and transfer learning in deep linear networks,” International Conference on Learning Representations (ICLR) (2019).
  * Li  _et al._ (2020) Zhu Li, Weijie J. Su, and Dino Sejdinovic, “Benign overfitting and noisy features,” (2020), arXiv:2008.02901 .
  * Liang and Rakhlin (2020) Tengyuan Liang and Alexander Rakhlin, “Just interpolate: Kernel “Ridgeless” regression can generalize,” Annals of Statistics 48, 1329–1347 (2020).
  * Liang  _et al._ (2020) Tengyuan Liang, Alexander Rakhlin, and Xiyu Zhai, “On the Multiple Descent of Minimum-Norm Interpolants and Restricted Lower Isometry of Kernels,” Proceedings of Thirty Third Conference on Learning Theory PMLR, 125, 2683–2711 (2020).
  * Liao  _et al._ (2020) Zhenyu Liao, Romain Couillet, and Michael W Mahoney, “A random matrix analysis of random Fourier features: beyond the Gaussian kernel, a precise phase transition, and the corresponding double descent,” Advances in Neural Information Processing Systems (NeurIPS) 33 (2020).
  * Lin and Dobriban (2021) Licong Lin and Edgar Dobriban, “What causes the test error? Going beyond bias-variance via ANOVA,” Journal of Machine Learning Research 22, 1–82 (2021).
  * Mitra (2019) Partha P Mitra, “Understanding overfitting peaks in generalization error: Analytical risk curves for $l_{2}$ and $l_{1}$ penalized interpolation,”  (2019), arXiv:1906.03667 .
  * Mei and Montanari (2021) 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 (2021), 10.1002/cpa.22008.
  * Muthukumar  _et al._ (2019) Vidya Muthukumar, Kailas Vodrahalli, and Anant Sahai, “Harmless interpolation of noisy data in regression,” IEEE International Symposium on Information Theory (ISIT) , 2299–2303 (2019).
  * Nakkiran (2019) Preetum Nakkiran, “More Data Can Hurt for Linear Regression: Sample-wise Double Descent,”  (2019), arXiv:1912.07242 .
  * Xu and Hsu (2019) Ji Xu and Daniel Hsu, “On the number of variables to use in principal component regression,” Advances in Neural Information Processing Systems (NeurIPS) 32 (2019).
  * Yang  _et al._ (2020) Zitong Yang, Yaodong Yu, Chong You, Jacob Steinhardt, and Yi Ma, “Rethinking Bias-Variance Trade-off for Generalization of Neural Networks,” Proceedings of the 37th International Conference on Machine Learning (ICML) PMLR, 119, 10767–10777 (2020).
  * Rocks and Mehta (2020) Jason W. Rocks and Pankaj Mehta, “Memorizing without overfitting: Bias, variance, and interpolation in over-parameterized models,” arXiv:2010.13933 (2020), arXiv:2010.13933 .
  * Dupic and Castillo (2014) Thomas Dupic and Isaac Pérez Castillo, “Spectral density of products of Wishart dilute random matrices. Part I: the dense case,”  (2014), arXiv:1401.7802 .
  * Cui  _et al._ (2020) Wenping Cui, Jason W. Rocks, and Pankaj Mehta, “The perturbative resolvent method: Spectral densities of random matrix ensembles via perturbation theory,” (2020), arXiv:2012.00663 .
  * Marčenko and Pastur (1967) Vladimir Alexandrovich Marčenko and Leonid Andreevich Pastur, “Distribution of eigenvalues for some sets of random matrices,” Mathematics of the USSR-Sbornik 1, 457–483 (1967).
  * Chaudhari  _et al._ (2019) Pratik Chaudhari, Anna Choromanska, Stefano Soatto, Yann LeCun, Carlo Baldassi, Christian Borgs, Jennifer Chayes, Levent Sagun, and Riccardo Zecchina, “Entropy-SGD: biasing gradient descent into wide valleys,” Journal of Statistical Mechanics: Theory and Experiment 2019, 124018 (2019).
  * Baldassi  _et al._ (2020) Carlo Baldassi, Enrico M. Malatesta, Matteo Negri, and Riccardo Zecchina, “Wide flat minima and optimal generalization in classifying high-dimensional Gaussian mixtures,” Journal of Statistical Mechanics: Theory and Experiment 2020, 124012 (2020).
  * Pittorino  _et al._ (2021) Fabrizio Pittorino, Carlo Lucibello, Christoph Feinauer, Enrico M. Malatesta, Gabriele Perugini, Carlo Baldassi, Matteo Negri, Elizaveta Demyanenko, and Riccardo Zecchina, “Entropic gradient descent algorithms and wide flat minima,” International Conference on Learning Representations (ICLR) (2021).

##  Appendix A Cavity Derivations

In this section, we provide detailed derivations of all closed-form solutions
for the random linear features model. These calculations follow the general
procedure laid out in Ref. Rocks and Mehta, 2020.

###  A.1 Notational Conventions

  * •

We define $M$ as the number of points in the training data set, $N_{f}$ as the
number of input features, and $N_{p}$ as the number of fit parameters/hidden
features. We define the ratios $\alpha_{f}=N_{f}/M$ and $\alpha_{p}=N_{p}/M$.

  * •

Unless otherwise specified, the type of symbol used for an index label (e.g.,
$\Delta y_{a}$) or as a summation index (e.g., $\sum_{a}$) implies its range.
The symbols $a$, $b$, or $c$ imply ranges over the training data points from
$1$ to $M$, the symbols $j$, $k$, or $l$ imply ranges over the input features
from $1$ to $N_{f}$, and the symbols $J$, $K$, or $L$ imply ranges over the
fit parameters/hidden features from $1$ to $N_{p}$.

  * •

The notation $\mathrm{E}_{x}\quantity[\cdot]$,
$\mathrm{Var}_{x}\quantity[\cdot]$ and
$\mathrm{Cov}_{x}\quantity[\cdot,\cdot]$ represent the mean, variance, and
covariance, respectively, with respect to one or more random variables $x$. A
lack of subscript implies averages taken with respect to the total ensemble
distribution, i.e., taken over all possible sources of randomness. A subscript
$0$ implies averages taken with respect to random variables containing one or
more $0$-valued indices (e.g., $X_{a0}$, $X_{0j}$, $W_{0J}$, or $W_{j0}$).

###  A.2 Nonlinear Label Decomposition

In order to calculate the statistical properties of the nonlinear teacher
model, we first decompose the labels into their linear and nonlinear
components as follows:

$y(\vec{\mathbf{x}})=\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})+\varepsilon,\qquad\vec{\boldsymbol{\beta}}\equiv\Sigma_{\vec{\mathbf{x}}}^{-1}\mathrm{Cov}_{\vec{\mathbf{x}}}[\vec{\mathbf{x}},y^{*}(\vec{\mathbf{x}})].$
(47)

The first term in this decomposition captures the linear correlations between
the labels $y^{*}$ and the input features $\vec{\mathbf{x}}$ via the ground
truth parameters $\vec{\boldsymbol{\beta}}$. The second term, defined as
${\delta y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})\equiv
y^{*}(\vec{\mathbf{x}})-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}}$,
represents the remaining nonlinear component of the labels. By defining the
ground truth parameters as shown above, where
$\Sigma_{\vec{\mathbf{x}}}=\mathrm{Cov}_{\vec{\mathbf{x}}}[\vec{\mathbf{x}},\vec{\mathbf{x}}^{T}]$
is the covariance matrix of the input features (assumed to be invertible), it
can be proven that the linear and nonlinear components are statistically
independent. Using the definition of the nonlinear component of the labels
$\delta y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})$, it is straightforward to show
that their mean and covariance in the thermodynamic limit evaluate to

$\mathrm{E}\quantity[\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{a})]=0,\qquad\mathrm{Cov}\quantity[\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{a}),\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{b})]=\sigma_{\delta
y^{*}}^{2}\delta_{ab},$ (48)

where $\vec{\mathbf{x}}_{a}$ and $\vec{\mathbf{x}}_{b}$ are to independent
data points and we have defined the variance $\sigma_{\delta y^{*}}^{2}$ of
the nonlinear components as

$\begin{gathered}\sigma_{\delta
y^{*}}^{2}=\sigma_{\beta}^{2}\sigma_{X}^{2}\Delta f,\qquad\Delta
f=\frac{\expectationvalue*{f^{2}}-\expectationvalue*{f^{\prime}}^{2}}{\expectationvalue*{f^{\prime}}^{2}}\\\
\expectationvalue*{f^{2}}=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}dhe^{-\frac{h^{2}}{2}}f^{2}(h),\qquad\expectationvalue*{f^{\prime}}=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}dhe^{-\frac{h^{2}}{2}}f^{\prime}(h).\end{gathered}$
(49)

The decomposition in Eq. (47) and its statistical properties are derived in
detail in Ref. Rocks and Mehta, 2020.

###  A.3 General Solutions

Next, we derive some useful formulas for the ensemble-averaged quantities we
wish to calculate. First, we express the ensemble-averaged training error as

$\expectationvalue{\mathcal{E}_{\mathrm{train}}}=\expectationvalue*{\Delta
y^{2}},\qquad\expectationvalue*{\Delta
y^{2}}=\mathrm{E}\quantity[\frac{1}{M}\sum_{b}\Delta y_{b}^{2}]$ (50)

where we have defined $\expectationvalue*{\Delta y^{2}}$ as the mean squared
label error for the training data.

Next, we evaluate the average over the test data set in the ensemble average
of the test error,

$\expectationvalue{\mathcal{E}_{\mathrm{test}}}=\sigma_{X}^{2}\expectationvalue*{\Delta\beta^{2}}+\sigma_{\delta
y^{*}}^{2}+\sigma_{\varepsilon}^{2}.$ (51)

To obtain this expression, we have defined the set of ground truth parameters
estimated by the model as $\hat{\boldsymbol{\beta}}\equiv W\hat{\mathbf{w}}$
and the corresponding residual parameter error
$\Delta\vec{\boldsymbol{\beta}}\equiv\vec{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}$.
The quantity $\expectationvalue*{\Delta\beta^{2}}$ is then the mean squared
residual parameter error,

$\expectationvalue*{\Delta\beta^{2}}=\mathrm{E}\quantity[\frac{1}{N_{f}}\sum_{k}\Delta\beta_{k}^{2}].$
(52)

To correctly calculate the ensemble average of the squared bias, we make use
of the following trick: we reinterpret the square of the average over
$\mathcal{D}$ as two separate averages over uncorrelated training data sets,

$\displaystyle\mathrm{Bias}^{2}[\hat{y}(\vec{\mathbf{x}})]$
$\displaystyle=(\mathrm{E}_{\mathcal{D}}\quantity[\hat{y}(\vec{\mathbf{x}})]-y^{*}(\vec{\mathbf{x}}))^{2}$
(53)
$\displaystyle=\mathrm{E}_{\mathcal{D}_{1},\mathcal{D}_{2}}\quantity[\quantity(y(\vec{\mathbf{x}})-\hat{y}_{1}(\vec{\mathbf{x}}))\quantity(y(\vec{\mathbf{x}})-\hat{y}_{2}(\vec{\mathbf{x}}))].$

Now, instead of a single regression problem trained on a single data set
$\mathcal{D}$, we consider two separate regression problems each trained
independently on different training sets, $\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$, drawn from the same distribution with the same ground truth
parameters $\vec{\boldsymbol{\beta}}$. These regression problems will also
share all other random variables including the test data point
$(y,\vec{\mathbf{x}})$, $W$, etc.

Next, we apply the ensemble average and explicitly average over the test data
point $\vec{\mathbf{x}}$ to obtain

$\expectationvalue*{\mathrm{Bias}^{2}[\hat{y}(\vec{\mathbf{x}})]}=\sigma_{X}^{2}\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}},\qquad\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}=\mathrm{E}\quantity[\frac{1}{N_{f}}\sum_{k}\Delta\beta_{1,k}\Delta_{2,k}]$
(54)

where we have defined $\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}$ as
the covariance of the residual label errors between the two models trained on
data sets $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$.

Finally, we find an expression for the variance by subtracting the bias and
noise ($\sigma_{\varepsilon}^{2}$) from the test error,

$\expectationvalue*{\mathrm{Var}[\hat{y}(\vec{\mathbf{x}})]}=\sigma_{X}^{2}\quantity(\expectationvalue*{\Delta\beta^{2}}-\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}).$
(55)

Based on these expressions, we find that the training error, test error, bias,
and variance depend on three key ensemble-averaged quantities:
$\expectationvalue*{\Delta y^{2}}$, $\expectationvalue*{\Delta\beta^{2}}$, and
$\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}$. We aim to calculate
these quantities in the remainder of this derivation.

###  A.4 Linear System of Equations

In this section, we derive a linear system of equations to which we will apply
the cavity method. To do this, we first evaluate the gradient of the loss
function in Eq. (5) with respect to the fit parameters,

$0=\partialderivative{L(\hat{\mathbf{w}})}{\hat{w}_{J}}=-\sum_{b}\Delta
y_{b}Z_{bJ}+\lambda\hat{w}_{J}.$ (56)

In addition to this gradient equation, we will also need the equations for the
residual label errors for the training set,

$\Delta
y_{a}=y^{*}(\vec{\mathbf{x}}_{a})+\varepsilon_{a}-\sum_{K}\hat{w}_{K}Z_{aJ}.$
(57)

Next, we decompose these two sets of equations such that they are linear in
the random matrices $W$ and $X$, resulting in four different sets of
equations,

$\displaystyle\lambda\hat{w}_{J}$
$\displaystyle=\sum_{k}\hat{u}_{k}W_{kJ}+\eta_{J}$ (58)
$\displaystyle\hat{u}_{j}$ $\displaystyle=\sum_{b}\Delta y_{b}X_{bj}+\psi_{j}$
$\displaystyle\Delta y_{a}$
$\displaystyle=\sum_{k}\Delta\beta_{k}X_{ak}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{a})+\varepsilon_{a}+\xi_{a}$
$\displaystyle\Delta\beta_{j}$
$\displaystyle=\beta_{j}-\sum_{K}\hat{w}_{K}W_{jK}+\zeta_{j},$

where we have also utilized Eq. (47) to decompose the training labels into
their linear and nonlinear components. We have also added a small auxiliary
field, $\eta_{J}$, $\psi_{j}$, $\xi_{a}$, or $\zeta_{j}$, to each equation. We
will use these extra fields to define perturbations about the solutions to
these equations with the intent of setting the fields to zero at the end of
the derivation. The quantities $\hat{u}_{j}$ can be interpreted as
representations of the fit parameters in the space of input features.

###  A.5 Cavity Expansion

Next, we add an additional variable of each type, resulting in a total of
$M+1$ data points, $N_{f}+1$ input features and $N_{p}+1$ fit parameters. Each
additional variable is represented using an index value of $0$, written as
$\hat{w}_{0}$, $\hat{u}_{0}$, $\Delta y_{0}$, and $\Delta\beta_{0}$. After
including these new unknown quantities, the four equations become

$\displaystyle\lambda\hat{w}_{J}$
$\displaystyle=\sum_{k}\hat{u}_{k}W_{kJ}+\eta_{J}+\hat{u}_{0}W_{0J}$ (59)
$\displaystyle\hat{u}_{j}$ $\displaystyle=\sum_{b}\Delta
y_{b}X_{bj}+\psi_{j}+\Delta y_{0}X_{0j}$ $\displaystyle\Delta y_{a}$
$\displaystyle=\sum_{k}\Delta\beta_{k}X_{ak}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{a})+\varepsilon_{a}+\xi_{a}+X_{a0}\Delta\beta_{0}$
$\displaystyle\Delta\beta_{j}$
$\displaystyle=\beta_{j}-\sum_{K}\hat{w}_{K}W_{jK}+\zeta_{j}-\hat{w}_{0}W_{j0},$

with each new variable described by a new equation,

$\displaystyle\lambda\hat{w}_{0}$
$\displaystyle=\sum_{k}\hat{u}_{k}W_{k0}+\eta_{0}+\hat{u}_{0}W_{00}$ (60)
$\displaystyle\hat{u}_{0}$ $\displaystyle=\sum_{b}\Delta
y_{b}X_{b0}+\psi_{0}+\Delta y_{0}X_{00}$ $\displaystyle\Delta y_{0}$
$\displaystyle=\sum_{k}\Delta\beta_{k}X_{0k}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{0})+\varepsilon_{0}+\xi_{0}+\Delta\beta_{0}X_{00}$
$\displaystyle\Delta\beta_{0}$
$\displaystyle=\beta_{0}-\sum_{K}\hat{w}_{K}W_{0K}+\zeta_{0}-\hat{w}_{0}W_{00}.$

Now we take the thermodynamic limit in which $M$, $N_{f}$, and $N_{p}$ tend
towards infinity, but their ratios, $\alpha_{f}=N_{f}/M$ and
$\alpha_{p}=N_{p}/M$, remain fixed. We interpret the extra terms in Eq. (59)
as small perturbations to the auxiliary fields,

$\delta\eta_{J}=\hat{u}_{0}W_{0J},\qquad\delta\psi_{j}=\Delta
y_{0}X_{0j},\qquad\delta\xi_{a}=\Delta\beta_{0}X_{a0},\qquad\delta\zeta_{j}=-\hat{w}_{0}W_{j0},$
(61)

allowing us to expand each unknown quantity about its solution in the absence
of the $0$-indexed variables (i.e., the solutions for $M$ data points, $N_{f}$
input features, and $N_{p}$ fit parameters),

$\displaystyle\hat{w}_{J}$ $\displaystyle\approx\hat{w}_{J\setminus
0}+\sum_{K}\nu^{\hat{w}}_{JK}\delta\eta_{K}+\sum_{k}\phi^{\hat{w}}_{Jk}\delta\psi_{k}+\sum_{b}\chi^{\hat{w}}_{Jb}\delta\xi_{b}+\sum_{k}\omega^{\hat{w}}_{Jk}\delta\zeta_{k}$
(62) $\displaystyle\hat{u}_{j}$ $\displaystyle\approx\hat{u}_{j\setminus
0}+\sum_{K}\nu^{\hat{u}}_{jK}\delta\eta_{K}+\sum_{k}\phi^{\hat{u}}_{jk}\delta\psi_{k}+\sum_{b}\chi^{\hat{u}}_{jb}\delta\xi_{b}+\sum_{k}\omega^{\hat{u}}_{jk}\delta\zeta_{k}$
$\displaystyle\Delta y_{a}$ $\displaystyle\approx\Delta y_{a\setminus
0}+\sum_{K}\nu^{\Delta y}_{aK}\delta\eta_{K}+\sum_{k}\phi^{\Delta
y}_{ak}\delta\psi_{k}+\sum_{b}\chi^{\Delta
y}_{ab}\delta\xi_{b}+\sum_{k}\omega^{\Delta y}_{ak}\delta\zeta_{k}$
$\displaystyle\Delta\beta_{j}$ $\displaystyle\approx\Delta\beta_{j\setminus
0}+\sum_{K}\nu^{\Delta\beta}_{jK}\delta\eta_{K}+\sum_{k}\phi^{\Delta\beta}_{jk}\delta\psi_{k}+\sum_{b}\chi^{\Delta\beta}_{jb}\delta\xi_{b}+\sum_{k}\omega^{\Delta\beta}_{jk}\delta\zeta_{k}.$

We define each of the susceptibility matrices as a derivative of a variable
with respect to an auxiliary field,

$\displaystyle\nu^{\hat{w}}_{JK}$
$\displaystyle=\partialderivative{\hat{w}_{J}}{\eta_{K}},\qquad$
$\displaystyle\phi^{\hat{w}}_{Jk}$
$\displaystyle=\partialderivative{\hat{w}_{J}}{\psi_{k}},\qquad$
$\displaystyle\chi^{\hat{w}}_{Jb}$
$\displaystyle=\partialderivative{\hat{w}_{J}}{\xi_{b}},\qquad$
$\displaystyle\omega^{\hat{w}}_{Jk}$
$\displaystyle=\partialderivative{\hat{w}_{J}}{\zeta_{k}},$ (63)
$\displaystyle\nu^{\hat{u}}_{jK}$
$\displaystyle=\partialderivative{\hat{u}_{j}}{\eta_{K}},\qquad$
$\displaystyle\phi^{\hat{u}}_{jk}$
$\displaystyle=\partialderivative{\hat{u}_{j}}{\psi_{k}},\qquad$
$\displaystyle\chi^{\hat{u}}_{jb}$
$\displaystyle=\partialderivative{\hat{u}_{j}}{\xi_{b}},\qquad$
$\displaystyle\omega^{\hat{u}}_{jk}$
$\displaystyle=\partialderivative{\hat{u}_{j}}{\zeta_{k}},$
$\displaystyle\nu^{\Delta y}_{aK}$ $\displaystyle=\partialderivative{\Delta
y_{a}}{\eta_{K}},\qquad$ $\displaystyle\phi^{\Delta y}_{ak}$
$\displaystyle=\partialderivative{\Delta y_{a}}{\psi_{k}},\qquad$
$\displaystyle\chi^{\Delta y}_{ab}$ $\displaystyle=\partialderivative{\Delta
y_{a}}{\xi_{b}},\qquad$ $\displaystyle\omega^{\Delta y}_{ak}$
$\displaystyle=\partialderivative{\Delta y_{a}}{\zeta_{k}},$
$\displaystyle\nu^{\Delta\beta}_{jK}$
$\displaystyle=\partialderivative{\Delta\beta_{j}}{\eta_{K}},\qquad$
$\displaystyle\phi^{\Delta\beta}_{jk}$
$\displaystyle=\partialderivative{\Delta\beta_{j}}{\psi_{k}},\qquad$
$\displaystyle\chi^{\Delta\beta}_{jb}$
$\displaystyle=\partialderivative{\Delta\beta_{j}}{\xi_{b}},\qquad$
$\displaystyle\omega^{\Delta\beta}_{jk}$
$\displaystyle=\partialderivative{\Delta\beta_{j}}{\zeta_{k}}.$

####  A.5.1 Central Limit Theorem

Substituting the expansions in Eq. (62) into the $0$-indexed equations in Eq.
(60), we find that each of the resulting sums contains a thermodynamically
large number of statistically uncorrelated terms. This means that each sum
satisfies the conditions necessary to apply the central limit theorem,
allowing us to express each in terms of a single normally-distributed random
variable described by just its mean and its variance.

First, we approximate each of the sums containing one of the unperturbed
quantities, $\hat{w}_{J\setminus 0}$, $\hat{u}_{j\setminus 0}$, $\Delta
y_{a\setminus 0}$, or $\Delta\beta_{j\setminus 0}$. The unperturbed quantities
in each of these sums are statistically independent of all elements of both
$X$ and $W$ with a $0$-valued index. Using this fact, we find

$\displaystyle\sum_{k}\hat{u}_{k\setminus 0}W_{k0}$
$\displaystyle\approx\sigma_{\hat{w}}z_{\hat{w}},\qquad$
$\displaystyle\sigma_{\hat{w}}^{2}$
$\displaystyle=\sigma_{W}^{2}\frac{\alpha_{f}}{\alpha_{p}}\expectationvalue*{\hat{u}^{2}},\qquad$
$\displaystyle\expectationvalue*{\hat{u}^{2}}$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\hat{u}_{k\setminus 0}^{2}$ (64)
$\displaystyle\sum_{b}\Delta y_{b\setminus 0}X_{b0}$
$\displaystyle\approx\sigma_{\hat{u}}z_{\hat{u}},\qquad$
$\displaystyle\sigma_{\hat{u}}^{2}$
$\displaystyle=\sigma_{X}^{2}\alpha_{f}^{-1}\expectationvalue*{\Delta
y^{2}},\qquad$ $\displaystyle\expectationvalue*{\Delta y^{2}}$
$\displaystyle=\frac{1}{M}\sum_{b}\Delta y_{b\setminus 0}^{2}$
$\displaystyle\sum_{k}\Delta\beta_{k\setminus 0}X_{0k}$
$\displaystyle\approx\sigma_{\Delta y}z_{\Delta y},\qquad$
$\displaystyle\sigma_{\Delta y}^{2}$
$\displaystyle=\sigma_{X}^{2}\expectationvalue*{\Delta\beta^{2}},\qquad$
$\displaystyle\expectationvalue*{\Delta\beta^{2}}$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\Delta\beta_{k\setminus 0}^{2}$
$\displaystyle\sum_{K}\hat{w}_{K\setminus 0}W_{0K}$
$\displaystyle\approx\sigma_{\Delta\beta}z_{\Delta\beta},\qquad$
$\displaystyle\sigma_{\Delta\beta}^{2}$
$\displaystyle=\sigma_{W}^{2}\expectationvalue*{\hat{w}^{2}},\qquad$
$\displaystyle\expectationvalue*{\hat{w}^{2}}$
$\displaystyle=\frac{1}{N_{p}}\sum_{K}\hat{w}_{K\setminus 0}^{2},$

where $z_{\hat{w}}$, $z_{\hat{u}}$, $z_{\Delta y}$, and $z_{\Delta\beta}$ are
all random variables with zero mean and unit variance and can easily be shown
to be statistically independent from one another.

Note that we have used the same notation, $\expectationvalue*{\Delta y^{2}}$
and $\expectationvalue*{\Delta\beta^{2}}$, for the two averages defined
previously in Sec. A.3 even though they each lack an ensemble average. In
doing so, we have employed the ansatz that these sums will converge to their
ensemble averages in the thermodynamic limit. This assumption is typical of
the cavity method.

Next, we approximate each of the sums containing one of the square
susceptibility matrices. Using the fact that all of the susceptibility
matrices are statistically independent of all elements of both $X$ and $W$
with a $0$-valued index, we find that each of these sums is dominated by its
mean with its variance going to zero in the thermodynamic limit,

$\displaystyle\sum_{jk}\omega^{\hat{u}}_{jk}W_{j0}W_{k0}$
$\displaystyle\approx\sigma_{W}^{2}\frac{\alpha_{f}}{\alpha_{p}}\omega,\qquad$
$\displaystyle\omega$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\omega^{\hat{u}}_{kk}$ (65)
$\displaystyle\sum_{ab}\chi^{\Delta y}_{ab}X_{a0}X_{b0}$
$\displaystyle\approx\sigma_{X}^{2}\alpha_{f}^{-1}\chi,\qquad$
$\displaystyle\chi$ $\displaystyle=\frac{1}{M}\sum_{b}\chi^{\Delta y}_{bb}$
$\displaystyle\sum_{jk}\phi^{\Delta\beta}_{jk}X_{0j}X_{0k}$
$\displaystyle\approx\sigma_{X}^{2}\phi,\qquad$ $\displaystyle\phi$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\phi^{\Delta\beta}_{kk}$
$\displaystyle\sum_{JK}\nu^{\hat{w}}_{JK}W_{0J}W_{0K}$
$\displaystyle\approx\sigma_{W}^{2}\nu,\qquad$ $\displaystyle\nu$
$\displaystyle=\frac{1}{N_{p}}\sum_{K}\nu^{\hat{w}}_{KK}$

where $\omega$, $\chi$, $\phi$, and $\nu$ can be interpreted as a set of
scalar susceptibilities.

Finally, it is straightforward to show that both the mean and variance each of
the sums containing a rectangular susceptibility matrix goes to zero in the
thermodynamic limit and can therefore be neglected.

####  A.5.2 Self-consistency Equations

Applying the approximations from the previous section, we find a set of self-
consistent equations for $\hat{w}_{0}$, $\hat{u}_{0}$, $\Delta y_{0}$, and
$\Delta\beta_{0}$,

$\displaystyle\lambda\hat{w}_{0}$
$\displaystyle\approx\sigma_{\hat{w}}z_{\hat{w}}-\hat{w}_{0}\sigma_{W}^{2}\frac{\alpha_{f}}{\alpha_{p}}\omega+\eta_{0}$
(66) $\displaystyle\hat{u}_{0}$
$\displaystyle\approx\sigma_{\hat{u}}z_{\hat{u}}+\Delta\beta_{0}\sigma_{X}^{2}\alpha_{f}^{-1}\chi+\psi_{0}$
$\displaystyle\Delta y_{0}$ $\displaystyle\approx\sigma_{\Delta y}z_{\Delta
y}+\Delta y_{0}\sigma_{X}^{2}\phi+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{0})+\varepsilon_{0}+\xi_{0}$
$\displaystyle\Delta\beta_{0}$
$\displaystyle\approx\beta_{0}-\sigma_{\Delta\beta}z_{\Delta\beta}-\hat{u}_{0}\sigma_{W}^{2}\nu+\zeta_{0},$

where have also made use of the fact that the terms including $X_{00}$ or
$W_{00}$ are infinitesimally small in the thermodynamic limit with zero mean
and variances of $\order{1/N_{f}}$ and $\order{1/N_{p}}$, respectively.
Solving these equations for the $0$-indexed variables, we find

$\displaystyle\hat{w}_{0}$
$\displaystyle=\frac{\sigma_{\hat{w}}z_{\hat{w}}+\eta_{0}}{\lambda+\sigma_{W}^{2}\frac{\alpha_{f}}{\alpha_{p}}\omega}$
(67) $\displaystyle\hat{u}_{0}$
$\displaystyle=\frac{\sigma_{\hat{u}}z_{\hat{u}}+\psi_{0}+\sigma_{X}^{2}\alpha_{f}^{-1}\chi\quantity(\beta_{0}-\sigma_{\Delta\beta}z_{\Delta\beta}+\zeta_{0})}{1+\sigma_{W}^{2}\sigma_{X}^{2}\alpha_{f}^{-1}\chi\nu}$
$\displaystyle\Delta y_{0}$ $\displaystyle=\frac{\sigma_{\Delta y}z_{\Delta
y}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{0})+\varepsilon_{0}+\xi_{0}}{1-\sigma_{X}^{2}\phi}$
$\displaystyle\Delta\beta_{0}$
$\displaystyle=\frac{\beta_{0}-\sigma_{\Delta\beta}z_{\Delta\beta}+\zeta_{0}-\sigma_{W}^{2}\nu^{2}\quantity(\sigma_{\hat{u}}z_{\hat{u}}+\psi_{0})}{1+\sigma_{W}^{2}\sigma_{X}^{2}\alpha_{f}^{-1}\chi\nu}.$

Next, we derive a set of self-consistent equations for the scalar
susceptibilities by taking appropriate derivatives of these variables with
respect to the auxiliary fields,

$\displaystyle\nu$
$\displaystyle=\frac{1}{N_{p}}\sum_{K}\nu^{\hat{w}}_{KK}\approx\mathrm{E}\quantity[\nu^{\hat{w}}_{00}]=\mathrm{E}\quantity[\partialderivative{\hat{w}_{0}}{\eta_{0}}]=\frac{1}{\lambda+\sigma_{W}^{2}\frac{\alpha_{f}}{\alpha_{p}}\omega}$
(68) $\displaystyle\omega$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\omega^{\hat{u}}_{kk}\approx\mathrm{E}\quantity[\omega^{\hat{u}}_{00}]=\mathrm{E}\quantity[\partialderivative{\hat{u}_{0}}{\zeta_{0}}]=\frac{\sigma_{X}^{2}\alpha_{f}^{-1}\chi}{1+\sigma_{W}^{2}\sigma_{X}^{2}\alpha_{f}^{-1}\chi\nu}$
$\displaystyle\chi$ $\displaystyle=\frac{1}{M}\sum_{b}\chi^{\Delta
y}_{bb}\approx\mathrm{E}\quantity[\chi^{\Delta
y}_{00}]=\mathrm{E}\quantity[\partialderivative{\Delta
y_{0}}{\xi_{0}}]=\frac{1}{1-\sigma_{X}^{2}\phi}$ $\displaystyle\phi$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\phi^{\Delta\beta}_{kk}\approx\mathrm{E}\quantity[\phi^{\Delta\beta}_{00}]=\mathrm{E}\quantity[\partialderivative{\Delta\beta_{0}}{\psi_{0}}]=-\frac{\sigma_{W}^{2}\nu}{1+\sigma_{W}^{2}\sigma_{X}^{2}\alpha_{f}^{-1}\chi\nu}.$

Furthermore, we note that there are two additional derivatives that have not
yet appeared in the calculation up to this point,
$\partial\hat{u}_{j}/\partial\psi_{j}$ and
$\partial\Delta\beta_{j}/\partial\zeta_{j}$. It is clear to see from the
equations for $\hat{u}_{0}$ and $\Delta\beta_{0}$ that these two additional
derivatives are equivalent. Evaluating these derivatives, we define a fifth
scalar susceptibility,

$\kappa=\frac{1}{N_{f}}\sum_{k}\phi^{\hat{u}}_{kk}=\frac{1}{N_{f}}\sum_{k}\omega^{\Delta\beta}_{kk}\approx\mathrm{E}\quantity[\partialderivative{\hat{u}_{0}}{\psi_{0}}]=\mathrm{E}\quantity[\partialderivative{\Delta\beta_{0}}{\zeta_{0}}]=\frac{1}{1+\sigma_{W}^{2}\sigma_{X}^{2}\alpha_{f}^{-1}\chi\nu}.$
(69)

Using this formula for $\kappa$, we re-express the four other susceptibilities
as

$\omega=\sigma_{X}^{2}\alpha_{f}^{-1}\chi\kappa,\qquad\phi=-\sigma_{W}^{2}\nu\kappa,\qquad\chi=\frac{1}{1+\sigma_{W}^{2}\sigma_{X}^{2}\nu\kappa},\qquad\nu=\frac{1}{\lambda+\sigma_{W}^{2}\sigma_{X}^{2}\alpha_{p}^{-1}\chi\kappa}.$
(70)

Finally, we square and average each of the expressions in Eq. (67) to find
self-consistent equations for the four mean squared averages (setting the
auxiliary fields to zero),

$\displaystyle\expectationvalue*{\hat{w}^{2}}$
$\displaystyle=\frac{1}{N_{p}}\sum_{K}\hat{w}_{K\setminus
0}^{2}\approx\mathrm{E}\quantity[\hat{w}_{0}^{2}]=\nu^{2}\sigma_{W}^{2}\frac{\alpha_{f}}{\alpha_{p}}\expectationvalue*{\hat{u}^{2}}$
(71) $\displaystyle\expectationvalue*{\hat{u}^{2}}$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\hat{u}_{k\setminus
0}^{2}\approx\mathrm{E}\quantity[\hat{u}_{0}^{2}]=\kappa^{2}\sigma_{X}^{2}\alpha_{f}^{-1}\expectationvalue*{\Delta
y^{2}}+\omega^{2}\quantity(\sigma_{\beta}^{2}+\sigma_{W}^{2}\expectationvalue*{\hat{w}^{2}})$
$\displaystyle\expectationvalue*{\Delta y^{2}}$
$\displaystyle=\frac{1}{M}\sum_{b}\Delta y_{b\setminus
0}^{2}\approx\mathrm{E}\quantity[\Delta
y_{0}^{2}]=\chi^{2}\quantity(\sigma_{X}^{2}\expectationvalue*{\Delta\beta^{2}}+\sigma_{\delta
y^{*}}^{2}+\sigma_{\varepsilon}^{2})$
$\displaystyle\expectationvalue*{\Delta\beta^{2}}$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\Delta\beta_{k\setminus
0}^{2}\approx\mathrm{E}\quantity[\Delta\beta_{0}^{2}]=\kappa^{2}\quantity(\sigma_{\beta}^{2}+\sigma_{W}^{2}\expectationvalue*{\hat{w}^{2}})+\phi^{2}\sigma_{X}^{2}\alpha_{f}^{-1}\expectationvalue*{\Delta
y^{2}}.$

####  A.5.3 Solution with Finite Regularization ($\lambda\sim 1$)

Next, we derive the solutions when the regularization parameter $\lambda$ is
finite. To do this, we combine the self-consistency equations for the
susceptibilities in Eqs. (69) and (70) to derive a cubic equation for $\chi$,

$0=\chi^{3}+(\alpha_{f}+\alpha_{p}-2)\chi^{2}+\quantity[\quantity(\alpha_{f}-1)\quantity(\alpha_{p}-1)+\alpha_{f}\alpha_{p}\bar{\lambda}]\chi-\alpha_{f}\alpha_{p}\bar{\lambda},$
(72)

where we have defined the dimensionless regularization parameter

$\bar{\lambda}=\frac{\lambda}{\sigma_{W}^{2}\sigma_{X}^{2}}.$ (73)

This cubic equation indicates that we should expect three different solutions
for $\chi$. Using these solutions, we can derive the associated solutions for
the rest of the susceptibilities. Furthermore, we solve Eq. (71) to find

$\matrixquantity(\expectationvalue*{\hat{w}^{2}}\\\
\expectationvalue*{\hat{u}^{2}}\\\ \expectationvalue*{\Delta y^{2}}\\\
\expectationvalue*{\Delta\beta^{2}})=\matrixquantity(1&-\sigma_{W}^{2}\frac{\alpha_{f}}{\alpha_{p}}\nu^{2}&0&0\\\
-\sigma_{W}^{2}\omega^{2}&1&-\sigma_{X}^{2}\alpha_{f}^{-1}\kappa^{2}&0\\\
0&0&1&-\sigma_{X}^{2}\chi^{2}\\\
-\sigma_{W}^{2}\kappa^{2}&0&-\sigma_{X}^{2}\alpha_{f}^{-1}\phi^{2}&1)^{-1}\matrixquantity(0\\\
\sigma_{\beta}^{2}\omega^{2}\\\ (\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\chi^{2}\\\ \sigma_{\beta}^{2}\kappa^{2}).$ (74)

In combination with the solutions for the five scalar susceptibilities, these
solutions are exact in the thermodynamic limit.

####  A.5.4 Solutions in Ridge-less Limit ($\lambda\rightarrow 0$)

Next, we take the ridge-less limit in which $\lambda\rightarrow 0$. Based on
the cubic equation for $\chi$ in Eq. (72), we make the ansatz that the lowest
order contribution to $\chi$ is $\order{1}$ in small $\bar{\lambda}$,

$\chi\approx\chi_{0}+\bar{\lambda}\chi_{1}.$ (75)

We then expand Eq. (72) in orders of $\lambda$ to find solutions for
$\chi_{0}$ and $\chi_{1}$. Using these solutions, we solve for the following
coefficients for the remaining susceptibilities.

$\displaystyle\nu$
$\displaystyle\approx\frac{1}{\bar{\lambda}}\nu_{-1}+\nu_{0}$ (76)
$\displaystyle\kappa$ $\displaystyle\approx\kappa_{0}+\bar{\lambda}\kappa_{1}$
$\displaystyle\phi$
$\displaystyle\approx\frac{1}{\bar{\lambda}}\phi_{-1}+\phi_{0}$
$\displaystyle\omega$
$\displaystyle\approx\omega_{0}+\bar{\lambda}\omega_{1},$

Finally, we expand the mean squared averages in small $\lambda$ as

$\displaystyle\expectationvalue*{\hat{w}^{2}}$
$\displaystyle\approx\expectationvalue*{\hat{w}^{2}}_{0}+\bar{\lambda}^{2}\expectationvalue*{\hat{w}^{2}}_{2}$
(77) $\displaystyle\expectationvalue*{\hat{u}^{2}}$
$\displaystyle\approx\expectationvalue*{\hat{u}^{2}}_{0}+\bar{\lambda}^{2}\expectationvalue*{\hat{u}^{2}}_{2}$
$\displaystyle\expectationvalue*{\Delta y^{2}}$
$\displaystyle\approx\expectationvalue*{\Delta
y^{2}}_{0}+\bar{\lambda}^{2}\expectationvalue*{\Delta y^{2}}_{2}$
$\displaystyle\expectationvalue*{\Delta\beta^{2}}$
$\displaystyle\approx\expectationvalue*{\Delta\beta^{2}}_{0}+\bar{\lambda}^{2}\expectationvalue*{\Delta\beta^{2}}_{2},$

and then use the solutions for the susceptibilities to solve Eq. (71) for
these coefficients.

We find three sets of solutions for all quantities, corresponding to the three
regimes of the random linear features model. To determine when each solution
applies, we use the fact that each of the ensemble-averaged quantities
$\expectationvalue*{\hat{w}^{2}}$, $\expectationvalue*{\hat{u}^{2}}$,
$\expectationvalue*{\Delta y^{2}}$, and $\expectationvalue*{\Delta\beta^{2}}$
must be positive. All together, we find the solutions for the ensemble-
averaged squared quantities in the $\lambda\rightarrow 0$ limit to be

$\displaystyle\expectationvalue*{\hat{w}^{2}}$
$\displaystyle=\left\\{\begin{array}[]{cl}\frac{\sigma_{\beta}^{2}}{\sigma_{W}^{2}}\frac{\alpha_{f}}{(\alpha_{p}-\alpha_{f})}+\frac{(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})}{\sigma_{W}^{2}\sigma_{X}^{2}}\frac{\alpha_{f}^{2}}{(1-\alpha_{f})(\alpha_{p}-\alpha_{f})}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
\frac{\sigma_{\beta}^{2}}{\sigma_{W}^{2}}\frac{(1-\alpha_{p}+\alpha_{f}-\alpha_{p})}{(1-\alpha_{p})(\alpha_{f}-\alpha_{p})}+\frac{(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})}{\sigma_{W}^{2}\sigma_{X}^{2}}\frac{\alpha_{f}\alpha_{p}}{(1-\alpha_{p})(\alpha_{f}-\alpha_{p})}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\
\frac{\sigma_{\beta}^{2}}{\sigma_{W}^{2}}\frac{1}{(\alpha_{p}-1)}+\frac{(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})}{\sigma_{W}^{2}\sigma_{X}^{2}}\frac{\alpha_{f}}{(\alpha_{f}-1)(\alpha_{p}-1)}&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.$ (81)
$\displaystyle\expectationvalue*{\hat{u}^{2}}$
$\displaystyle=\left\\{\begin{array}[]{cl}\frac{\lambda^{2}}{\sigma_{X}^{4}\sigma_{W}^{4}}\quantity[\sigma_{\beta}^{2}\sigma_{X}^{4}\frac{\alpha_{p}^{3}}{(\alpha_{p}-\alpha_{f})^{3}}+\sigma_{X}^{2}(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{\alpha_{f}\alpha_{p}^{3}}{(1-\alpha_{f})(\alpha_{p}-\alpha_{f})^{3}}]&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
\sigma_{\beta}^{2}\sigma_{X}^{4}\frac{(1-\alpha_{p})(\alpha_{f}-\alpha_{p})(1-\alpha_{p}+\alpha_{f}-\alpha_{p})}{\alpha_{f}^{3}}+\sigma_{X}^{2}(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{(1-\alpha_{p})(\alpha_{f}-\alpha_{p})}{\alpha_{f}^{2}}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\
\frac{\lambda^{2}}{\sigma_{X}^{4}\sigma_{W}^{4}}\quantity[\sigma_{\beta}^{2}\sigma_{X}^{4}\frac{\alpha_{p}^{3}}{\alpha_{f}(\alpha_{p}-1)^{3}}+\sigma_{X}^{2}(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{\alpha_{p}^{3}}{(\alpha_{f}-1)(\alpha_{p}-1)^{3}}]&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.$ (85)
$\displaystyle\expectationvalue*{\Delta y^{2}}$
$\displaystyle=\left\\{\begin{array}[]{cl}(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})(1-\alpha_{f})&\mbox{\quad if\quad}N_{f}<N_{p},M\\\
\sigma_{\beta}^{2}\sigma_{X}^{2}\frac{(1-\alpha_{p})(\alpha_{f}-\alpha_{p})}{\alpha_{f}}+(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})(1-\alpha_{p})&\mbox{\quad if\quad}N_{p}<N_{f},M\\\
\frac{\lambda^{2}}{\sigma_{X}^{4}\sigma_{W}^{4}}\quantity[\sigma_{\beta}^{2}\sigma_{X}^{2}\frac{\alpha_{f}\alpha_{p}^{3}}{(\alpha_{f}-1)(\alpha_{p}-1)^{3}}+(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})\frac{\alpha_{f}^{2}\alpha_{p}^{2}(\alpha_{f}-1+\alpha_{p}-1)}{(\alpha_{f}-1)^{3}(\alpha_{p}-1)^{3}}]&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.$ (89)
$\displaystyle\expectationvalue*{\Delta\beta^{2}}$
$\displaystyle=\left\\{\begin{array}[]{cl}\frac{(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})}{\sigma_{X}^{2}}\frac{\alpha_{f}}{(1-\alpha_{f})}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
\sigma_{\beta}^{2}\frac{(\alpha_{f}-\alpha_{p})}{\alpha_{f}(1-\alpha_{p})}+\frac{(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})}{\sigma_{X}^{2}}\frac{\alpha_{p}}{(1-\alpha_{p})}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\
\sigma_{\beta}^{2}\frac{\alpha_{p}(\alpha_{f}-1)}{\alpha_{f}(\alpha_{p}-1)}+\frac{(\sigma_{\varepsilon}^{2}+\sigma_{\delta
y^{*}}^{2})}{\sigma_{X}^{2}}\frac{(\alpha_{f}-1+\alpha_{p}-1)}{(\alpha_{f}-1)(\alpha_{p}-1)}&\mbox{\quad
if\quad}M<N_{f},N_{p}.\end{array}\right.$ (93)

In addition, to lowest order in small $\lambda$, the five scalar
susceptibilities are

$\displaystyle\chi$
$\displaystyle=\left\\{\begin{array}[]{cl}1-\alpha_{f}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\ 1-\alpha_{p}&\mbox{\quad if\quad}N_{p}<N_{f},M\\\
\frac{\lambda}{\sigma_{W}^{2}\sigma_{X}^{2}}\frac{\alpha_{f}\alpha_{p}}{(1-\alpha_{f})(1-\alpha_{p})}&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.$ (97) $\displaystyle\nu$
$\displaystyle=\left\\{\begin{array}[]{cl}\frac{1}{\lambda}\frac{(\alpha_{p}-\alpha_{f})}{\alpha_{p}}+\frac{1}{\sigma_{W}^{2}\sigma_{X}^{2}}\frac{\alpha_{f}^{2}}{\quantity(1-\alpha_{f})\quantity(\alpha_{p}-\alpha_{f})}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
\frac{1}{\sigma_{W}^{2}\sigma_{X}^{2}}\frac{\alpha_{f}\alpha_{p}}{(1-\alpha_{p})(\alpha_{f}-\alpha_{p})}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\
\frac{1}{\lambda}\frac{(\alpha_{p}-1)}{\alpha_{p}}+\frac{1}{\sigma_{W}^{2}\sigma_{X}^{2}}\frac{\alpha_{f}}{\quantity(\alpha_{f}-1)\quantity(\alpha_{p}-1)}&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.$ (101) $\displaystyle\kappa$
$\displaystyle=\left\\{\begin{array}[]{cl}\frac{\lambda}{\sigma_{W}^{2}\sigma_{X}^{2}}\frac{\alpha_{f}\alpha_{p}}{(1-\alpha_{f})(\alpha_{p}-\alpha_{f})}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
\frac{(\alpha_{f}-\alpha_{p})}{\alpha_{f}}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\ \frac{(\alpha_{f}-1)}{\alpha_{f}}&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.$ (105) $\displaystyle\omega$
$\displaystyle=\left\\{\begin{array}[]{cl}\frac{\lambda}{\sigma_{W}^{2}}\frac{\alpha_{p}}{(\alpha_{p}-\alpha_{f})}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
\sigma_{X}^{2}\frac{(1-\alpha_{p})(\alpha_{f}-\alpha_{p})}{\alpha_{f}^{2}}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\
\frac{\lambda}{\sigma_{W}^{2}}\frac{\alpha_{p}}{\alpha_{f}(\alpha_{p}-1)}&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.$ (109) $\displaystyle\phi$
$\displaystyle=\left\\{\begin{array}[]{cl}-\frac{1}{\sigma_{X}^{2}}\frac{\alpha_{f}}{(1-\alpha_{f})}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
-\frac{1}{\sigma_{X}^{2}}\frac{\alpha_{p}}{(1-\alpha_{p})}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\
-\frac{\sigma_{W}^{2}}{\lambda}\frac{(\alpha_{f}-1)(\alpha_{p}-1)}{\alpha_{f}\alpha_{p}}&\mbox{\quad
if\quad}M<N_{f},N_{p}.\end{array}\right.$ (113)

We use the quantities $\expectationvalue*{\Delta y^{2}}$ and
$\expectationvalue*{\Delta\beta^{2}}$ above in combination with formulas for
the training and test error in Sec. A.3 to obtain the expressions in Eqs. (20)
and (27).

####  A.5.5 Bias-Variance Decomposition

Next, we derive the bias and variance. According to the general solutions in
Eqs. (54) and (55), we require the quantity
$\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}$. To calculate
$\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}$, we apply the self-
consistent equations, Eq. (67), to two models each trained separately on one
of two independent training sets, with all other random variables held in
common. We specify which quantities depend on each of the training sets using
a subscript $1$ or $2$ for data sets $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$,
respectively. For training set $\mathcal{D}_{1}$, these equations are

$\displaystyle\hat{w}_{1,0}$
$\displaystyle=\nu\sigma_{\hat{w}}z_{\hat{w}_{1}}$ (114)
$\displaystyle\hat{u}_{1,0}$
$\displaystyle=\kappa\sigma_{\hat{u}}z_{\hat{u}_{1}}+\omega\quantity(\beta_{0}-\sigma_{\Delta\beta}z_{\Delta\beta_{1}})$
$\displaystyle\Delta y_{1,0}$ $\displaystyle=\chi\quantity(\sigma_{\Delta
y}z_{\Delta y_{1}}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{1,0})+\varepsilon_{1,0})$
$\displaystyle\Delta\beta_{1,0}$
$\displaystyle=\kappa\quantity(\beta_{0}-\sigma_{\Delta\beta}z_{\Delta\beta_{1}})+\phi\sigma_{\hat{u}}z_{\hat{u}_{1}},$

while for training set $\mathcal{D}_{2}$, they are

$\displaystyle\hat{w}_{2,0}$
$\displaystyle=\nu\sigma_{\hat{w}}z_{\hat{w}_{2}}$ (115)
$\displaystyle\hat{u}_{2,0}$
$\displaystyle=\kappa\sigma_{\hat{u}}z_{\hat{u}_{2}}+\omega\quantity(\beta_{0}-\sigma_{\Delta\beta}z_{\Delta\beta_{2}})$
$\displaystyle\Delta y_{2,0}$ $\displaystyle=\chi\quantity(\sigma_{\Delta
y}z_{\Delta y_{2}}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}_{2,0})+\varepsilon_{2,0})$
$\displaystyle\Delta\beta_{2,0}$
$\displaystyle=\kappa\quantity(\beta_{0}-\sigma_{\Delta\beta}z_{\Delta\beta_{2}})+\phi\sigma_{\hat{u}}z_{\hat{u}_{2}}.$

Multiplying these equations and making the self-averaging approximation, we
find equations for the covariance of each of the unknown variables,

$\displaystyle\expectationvalue*{\hat{w}_{1}\hat{w}_{2}}$
$\displaystyle=\frac{1}{N_{p}}\sum_{K}\hat{w}_{1,K}\hat{w}_{2,K}\approx\mathrm{E}\quantity[\hat{w}_{1,0}\hat{w}_{2,0}]=\nu^{2}\mathrm{E}\quantity[\sigma_{\hat{w}}^{2}z_{\hat{w}_{1}}z_{\hat{w}_{2}}]$
(116) $\displaystyle\expectationvalue*{\hat{u}_{1}\hat{u}_{2}}$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\hat{u}_{1,k}\hat{u}_{2,k}\approx\mathrm{E}\quantity[\hat{u}_{1,0}\hat{u}_{2,0}]=\kappa^{2}\mathrm{E}\quantity[\sigma_{\hat{u}}^{2}z_{\hat{u}_{1}}z_{\hat{u}_{2}}]+\omega^{2}\quantity(\sigma_{\beta}^{2}+\mathrm{E}\quantity[\sigma_{\Delta\beta}^{2}z_{\Delta\beta_{1}}z_{\Delta\beta_{2}}])$
$\displaystyle\expectationvalue*{\Delta y_{1}\Delta y_{2}}$
$\displaystyle=\frac{1}{M}\sum_{b}\Delta y_{1,b}\Delta
y_{2,b}\approx\mathrm{E}\quantity[\Delta y_{1,0}\Delta
y_{2,0}]=\chi^{2}\mathrm{E}\quantity[\sigma_{\Delta y}^{2}z_{\Delta
y_{1}}z_{\Delta y_{2}}]$
$\displaystyle\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}$
$\displaystyle=\frac{1}{N_{f}}\sum_{k}\Delta\beta_{1,k}\Delta\beta_{2,k}\approx\mathrm{E}\quantity[\Delta\beta_{1,0}\Delta\beta_{2,0}]=\kappa^{2}\quantity(\sigma_{\beta}^{2}+\mathrm{E}\quantity[\sigma_{\Delta\beta}^{2}z_{\Delta\beta_{1}}z_{\Delta\beta_{2}}])+\phi^{2}\mathrm{E}\quantity[\sigma_{\hat{u}}^{2}z_{\hat{u}_{1}}z_{\hat{u}_{2}}].$

Next, we calculate each of the four resulting expectation values of products
of random variables. Converting each of the random variables
$z_{\hat{w}_{1}}$, $z_{\Delta\beta_{1}}$, etc., back into their forms as sums,
we use the independence of elements of the random matrices and the other
variables to find

$\displaystyle\mathrm{E}\quantity[\sigma_{\hat{w}}^{2}z_{\hat{w}_{1}}z_{\hat{w}_{2}}]$
$\displaystyle\approx\mathrm{E}\quantity[\sum_{jk}\hat{u}_{1,j\setminus
0}\hat{u}_{2,k\setminus
0}W_{j0}W_{k0}]=\sigma_{W}^{2}\frac{\alpha_{f}}{\alpha_{p}}\expectationvalue*{\hat{u}_{1}\hat{u}_{2}}$
(117)
$\displaystyle\mathrm{E}\quantity[\sigma_{\Delta\beta}^{2}z_{\Delta\beta_{1}}z_{\Delta\beta_{2}}]$
$\displaystyle\approx\mathrm{E}\quantity[\sum_{JK}\hat{w}_{1,J\setminus
0}\hat{w}_{2,K\setminus
0}W_{0J}W_{0K}]=\sigma_{W}^{2}\expectationvalue*{\hat{w}_{1}\hat{w}_{2}}$
$\displaystyle\mathrm{E}\quantity[\sigma_{\hat{u}}^{2}z_{\hat{u}_{1}}z_{\hat{u}_{2}}]$
$\displaystyle\approx\mathrm{E}\quantity[\sum_{ab}\Delta y_{1,a\setminus
0}\Delta y_{2,b\setminus 0}X_{1,a0}X_{2,b0}]=0$
$\displaystyle\mathrm{E}\quantity[\sigma_{\Delta y}^{2}z_{\Delta
y_{1}}z_{\Delta y_{2}}]$
$\displaystyle\approx\mathrm{E}\quantity[\sum_{jk}\Delta\beta_{1,j\setminus
0}\Delta\beta_{2,k\setminus 0}X_{1,0j}X_{2,0k}]=0.$

Substituting these results back into Eq. (116), we find the self-consistent
equations

$\displaystyle\expectationvalue*{\hat{w}_{1}\hat{w}_{2}}$
$\displaystyle=\nu^{2}\sigma_{W}^{2}\frac{\alpha_{f}}{\alpha_{p}}\expectationvalue*{\hat{u}_{1}\hat{u}_{2}}$
(118) $\displaystyle\expectationvalue*{\hat{u}_{1}\hat{u}_{2}}$
$\displaystyle=\omega^{2}\quantity(\sigma_{\beta}^{2}+\sigma_{W}^{2}\expectationvalue*{\hat{w}_{1}\hat{w}_{2}})$
$\displaystyle\expectationvalue*{\Delta y_{1}\Delta y_{2}}$ $\displaystyle=0$
$\displaystyle\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}$
$\displaystyle=\kappa^{2}\quantity(\sigma_{\beta}^{2}+\sigma_{W}^{2}\expectationvalue*{\hat{w}_{1}\hat{w}_{2}}).$

Next, we make the ansatz that the ensemble-averaged covariances are
$\order{1}$ in small $\bar{\lambda}$ with the next order terms at
$\order{\bar{\lambda}^{2}}$,

$\displaystyle\expectationvalue*{\hat{w}_{1}\hat{w}_{2}}$
$\displaystyle\approx\expectationvalue*{\hat{w}_{1}\hat{w}_{2}}_{0}+\bar{\lambda}^{2}\expectationvalue*{\hat{w}_{1}\hat{w}_{2}}_{2}$
(119) $\displaystyle\expectationvalue*{\hat{u}_{1}\hat{u}_{2}}$
$\displaystyle\approx\expectationvalue*{\hat{u}_{1}\hat{u}_{2}}_{0}+\bar{\lambda}^{2}\expectationvalue*{\hat{u}_{1}\hat{u}_{2}}_{2}$
$\displaystyle\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}$
$\displaystyle\approx\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}_{0}+\bar{\lambda}^{2}\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}_{2}.$

All together, the covariances in the limit $\lambda\rightarrow 0$ are

$\displaystyle\expectationvalue*{\hat{w}_{1}\hat{w}_{2}}$
$\displaystyle=\left\\{\begin{array}[]{cl}\frac{\sigma_{\beta}^{2}}{\sigma_{W}^{2}}\frac{\alpha_{f}}{(\alpha_{p}-\alpha_{f})}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
\frac{\sigma_{\beta}^{2}}{\sigma_{W}^{2}}\frac{\alpha_{p}}{(\alpha_{f}-\alpha_{p})}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\
\frac{\sigma_{\beta}^{2}}{\sigma_{W}^{2}}\frac{1}{(\alpha_{f}\alpha_{p}-1)}&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.$ (123)
$\displaystyle\expectationvalue*{\hat{u}_{1}\hat{u}_{2}}$
$\displaystyle=\left\\{\begin{array}[]{cl}\frac{\lambda^{2}}{\sigma_{X}^{4}\sigma_{W}^{4}}\sigma_{\beta}^{2}\sigma_{X}^{4}\frac{\alpha_{p}^{3}}{(\alpha_{p}-\alpha_{f})^{3}}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
\sigma_{\beta}^{2}\sigma_{X}^{4}\frac{(1-\alpha_{p})^{2}(\alpha_{f}-\alpha_{p})}{\alpha_{f}^{3}}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\
\frac{\lambda^{2}}{\sigma_{X}^{4}\sigma_{W}^{4}}\sigma_{\beta}^{2}\sigma_{X}^{4}\frac{\alpha_{p}^{3}}{\alpha_{f}(\alpha_{p}-1)^{2}(\alpha_{f}\alpha_{p}-1)}&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.$ (127)
$\displaystyle\expectationvalue*{\Delta y_{1}\Delta y_{2}}$ $\displaystyle=0$
(128) $\displaystyle\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}$
$\displaystyle=\left\\{\begin{array}[]{cl}\frac{\lambda^{2}}{\sigma_{X}^{4}\sigma_{W}^{4}}\sigma_{\beta}^{2}\frac{\alpha_{f}^{2}\alpha_{p}^{3}}{(1-\alpha_{f})^{2}(\alpha_{p}-\alpha_{f})^{3}}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\
\sigma_{\beta}^{2}\frac{(\alpha_{f}-\alpha_{p})}{\alpha_{f}}&\mbox{\quad
if\quad}N_{p}<N_{f},M\\\
\sigma_{\beta}^{2}\frac{\alpha_{p}(\alpha_{f}-1)^{2}}{\alpha_{f}(\alpha_{f}\alpha_{p}-1)}&\mbox{\quad
if\quad}M<N_{f},N_{p}.\end{array}\right.$ (132)

Finally, we use the solution for
$\expectationvalue*{\Delta\beta_{1}\Delta\beta_{2}}$ to find the bias and
variance according to Eqs. (54) and (55), resulting in Eqs. (34) and (41).

##  Appendix B Spectral Densities of Kernel Matrices

Here, we derive the spectral densities for the kernel matrix $Z^{T}Z$ using
the technique laid out in Ref. Cui  _et al._ , 2020. According to this
formalism, the spectral density of the kernel can be written in terms of the
scalar susceptibility $\nu$, defined in the previous section, using the
formula

$\displaystyle\rho(x)$
$\displaystyle=-\frac{1}{\pi}\lim_{\varepsilon\rightarrow
0^{+}}\imaginary\nu(-x+i\varepsilon).$ (133)

In addition, we expect there to be delta function of eigenvalues located at
zero. Although the above formula can in principle be used to obtain the
fraction of eigenvalues at zero, for convenience, we instead use the scalar
susceptibility $\chi$, which can be shown to be exactly

$\chi=\frac{1}{M}\Tr\partialderivative{\Delta\vec{\mathbf{y}}}{\vec{\boldsymbol{\xi}}}=\frac{1}{M}\Tr\quantity[I_{M}-ZZ^{+}]=1-\frac{1}{M}\rank(Z^{T}Z).$
(134)

The fraction of eigenvalues at zero is then

$\displaystyle
f_{\mathrm{zero}}=1-\frac{1}{N_{p}}\rank(Z^{T}Z)=\frac{\chi+\alpha_{p}-1}{\alpha_{p}}.$
(135)

Next, we define dimensional versions of $\nu$ and $\lambda$,

$\bar{\nu}=\sigma_{W}^{2}\sigma_{X}^{2}\nu,\qquad\bar{\lambda}=\frac{\lambda}{\sigma_{W}^{2}\sigma_{X}^{2}}.$
(136)

Using the self-consistent equations for the scalar susceptibilities in Eqs.
(69) and (70), we find a cubic equation for $\bar{\nu}$,

$0=(\alpha_{p}\bar{\lambda}\bar{\nu})^{3}+\quantity[1-\alpha_{p}+\alpha_{f}-\alpha_{p}](\alpha_{p}\bar{\lambda}\bar{\nu})^{2}+\quantity[(1-\alpha_{p})(\alpha_{f}-\alpha_{p})+\alpha_{f}\alpha_{p}\bar{\lambda}](\alpha_{p}\bar{\lambda}\bar{\nu})-\alpha_{f}\alpha_{p}^{2}\bar{\lambda}.$
(137)

Solving this cubic equation analytically is very involved, so we refer to the
solution in Ref. Dupic and Castillo, 2014. Instead, we solve this equation
numerically for the negative imaginary roots of $\nu(\lambda)$ with
$\lambda=-x$, according to Eq. (133). However, we also need to find the
interval over which the eigenvalue spectrum is positive. To do this, we
rewrite the equation in general form for $\alpha_{p}\bar{\lambda}\bar{\nu}$,

$(\alpha_{p}\bar{\lambda}\bar{\nu})^{3}+a_{2}(\alpha_{p}\bar{\lambda}\bar{\nu})^{2}+a_{1}(\alpha_{p}\bar{\lambda}\bar{\nu})+a_{0}=0,$
(138)

where the coefficients are

$\displaystyle a_{0}$ $\displaystyle=-\alpha_{f}\alpha_{p}^{2}\bar{\lambda}$
(139) $\displaystyle a_{1}$
$\displaystyle=(1-\alpha_{p})(\alpha_{f}-\alpha_{p})+\alpha_{f}\alpha_{p}\bar{\lambda}$
$\displaystyle a_{2}$ $\displaystyle=1-\alpha_{p}+\alpha_{f}-\alpha_{p}.$

The discriminant for a cubic equation is expressed in terms of these
coefficients as

$D(\lambda)=R^{2}-Q^{3}$ (140)

with

$\displaystyle Q$ $\displaystyle=\frac{1}{9}\quantity(a_{2}^{2}-3a_{1})$ (141)
$\displaystyle R$
$\displaystyle=\frac{1}{54}\quantity(9a_{2}a_{1}-27a_{0}-2a_{2}^{3}).$

To find the limiting eigenvalues, we then solve the equation $D(\lambda)=0$
(with $\lambda=-x$) numerically for the largest and smallest non-negative real
roots.

To find the weight of the delta function component at zero, we use the
solution for $\chi$ that we found previously, giving us

$f_{\mathrm{zero}}=\left\\{\begin{array}[]{cl}1-\frac{\alpha_{f}}{\alpha_{p}}&\mbox{\quad
if\quad}N_{f}<N_{p},M\\\ 0&\mbox{\quad if\quad}N_{p}<N_{f},M\\\
1-\alpha_{p}^{-1}&\mbox{\quad
if\quad}M<N_{f},N_{p}\end{array}\right.=\max\quantity(0,1-\frac{\alpha_{f}}{\alpha_{p}},1-\alpha_{p}^{-1}).$
(142)

##  Appendix C Numerical Simulation Details

In this section, we explain our procedures for generating numerical results.
Fig. 4 provides comparisons to numerical results for the training error, test
error, bias, and variance.

###  C.1 General Details

In all plots of training error, test error, bias, and variance, each point (or
pixel for 2$d$ plots) is averaged over $1000$ independent simulations, unless
located exactly at a phase transition, in which case, each point is averaged
over $150000$ simulations. Small error bars are shown each plot, representing
the error on the mean. We also scale the error in each plot by the variance of
the labels $\sigma_{y}^{2}=\sigma_{\beta}^{2}\sigma_{X}^{2}+\sigma_{\delta
y^{*}}^{2}+\sigma_{\varepsilon}^{2}$. In all simulations, we use training and
test sets of size $M=M^{\prime}=512$, a signal-to-noise ratio of
$(\sigma_{\beta}^{2}\sigma_{X}^{2}+\sigma_{\delta
y^{*}}^{2})/\sigma_{\varepsilon}^{2}=10$, and a regularization parameter of
$\lambda=10^{-6}$. We use a linear teacher model
$y^{*}(\vec{\mathbf{x}})=\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}$
($\sigma_{\delta y^{*}}^{2}=0$) for all plots.

To find the solution for a particular regression problem, we solve a different
(but equivalent) system of equations depending on whether $N_{p}<M$ or
$N_{p}>M$, allowing us to reduce the size of the linear system we need to
solve. If $N_{p}<M$, we solve the system of $N_{p}$ equations

$\displaystyle\quantity[\lambda I_{N_{p}}+Z^{T}Z]\hat{\mathbf{w}}$
$\displaystyle=Z^{T}\vec{\mathbf{y}}$ (143)

for the $N_{p}$ unknown fit parameters $\hat{\mathbf{w}}$ where $I_{N_{p}}$ is
the $N_{p}\times N_{p}$ identity matrix. This equation is identical to that in
Eq. (6) in the main text.

Alternatively, if $N_{p}>M$ we solve a system of $M$ equations,

$\displaystyle\quantity[\lambda I_{M}+ZZ^{T}]\hat{\mathbf{a}}$
$\displaystyle=\vec{\mathbf{y}},$ (144)

for the $M$ unknowns $\hat{\mathbf{a}}$ where $I_{M}$ is the $M\times M$
identity matrix. We then convert to fit parameters via the formula
$\hat{\mathbf{w}}=Z^{T}\hat{\mathbf{a}}$.

###  C.2 Bias-Variance Decompositions

To efficiently calculate the ensemble-averaged bias and variance, we take
inspiration from Eq. (53). During each simulation, we independently generate
two training data sets $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$. Using the
results from the first training set, we calculate the training and test error.
To calculate the bias, we also calculate the label predictions for both
training sets for an identical test set, $\hat{\mathbf{y}}_{1}$ and
$\hat{\mathbf{y}}_{2}$, and record the residual label errors between these
predictions and the true labels of the test set $\vec{\mathbf{y}}^{*\prime}$
We then record the dot product
$(\hat{\mathbf{y}}_{1}-\vec{\mathbf{y}}^{*\prime})\cdot(\hat{\mathbf{y}}_{2}-\vec{\mathbf{y}}^{*\prime})$.
When averaged over many simulations, this quantity approximates the bias. We
can then subtract this quantity from the average test error to find the
variance.

###  C.3 Eigenvalue Decompositions of Kernel Matrices

For each of the numerical eigenvalue distributions for the kernel matrices
presented in the main text, we choose $M=4096$. We then average over the
distributions for 10 independently sampled matrices when $\alpha_{p}=1$ or
$\alpha_{p}=8$ and over 80 matrices when $\alpha_{p}=1/8$. In this way, we
ensure that the same number of non-zero eigenvalues is present in the part of
the histograms corresponding to the bulk of the distributions (the
distribution excluding the delta function at zero). For $M<N_{p}$ we calculate
the eigenvalues of $Z^{T}Z$, while for $M>N_{p}$ we instead calculate the
eigenvalues of $ZZ^{T}$ since this matrix is smaller and contains the same
non-zero eigenvalues. In the later case, we then manally append an additional
$N_{p}-M$ zero-valued eigenvalues to the distribution.

Figure 4: Comparison of analytic and numerical results for the random linear
features model: Training error and bias-variance decomposition. (Top Row)
Analytic solutions and (Bottom Row) numerical results are shown as a function
of $\alpha_{p}=N_{p}/M$ and $\alpha_{f}=N_{f}/M$. Plotted are the ensemble-
averaged (a) training error, (b) test error, (c) squared bias, and (d)
variance. In each panel, black dashed lines show boundaries between different
regimes of solutions depending on which is smallest of the quantities $M$,
$N_{f}$, or $N_{p}$. The vertical and horizontal lines bound the
interpolation, or overparameterized, regime, located at $\alpha_{p}>1$ and
$\alpha_{f}>1$, while the diagonal line marks the boundary between the large
bias and minimal bias underparameterized regimes for a linear teacher model.
All solutions have been scaled by the variance of the training set labels
$\sigma_{y}^{2}=\sigma_{\beta}^{2}\sigma_{X}^{2}+\sigma_{\varepsilon}^{2}$.