Bias-variance decomposition of overparameterized regression with random linear features

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.

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.

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

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

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)

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

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,

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

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,

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,

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),

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

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

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.

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

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).

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

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

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}$).

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,

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

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

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.

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}$.

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:

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

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

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,

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,

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,

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

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,

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,

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

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,

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 , 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

with each new variable described by a new equation,

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,

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

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

Substituting the expansions in Eq. (62) into the -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 -valued index. Using this fact, we find

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 -valued index, we find that each of these sums is dominated by its mean with its variance going to zero in the thermodynamic limit,

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}$,

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 -indexed variables, we find

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,

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,

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

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),

A.5.3 Solution with Finite Regularization (λ∼1similar-to𝜆1\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$,

where we have defined the dimensionless regularization parameter

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

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 (λ→0→𝜆0\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}$,

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.

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

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

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

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

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

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

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

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

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}}$,

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

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

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

The fraction of eigenvalues at zero is then

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

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

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}$,

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

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

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.

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

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,

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.