Multiple Descent: Design Your Own Generalization Curve

Introduction

The main goal of machine learning methods is to provide an accurate out-of-sample prediction, known as generalization. For a fixed family of models, a common way to select a model from this family is through empirical risk minimization, i.e., algorithmically selecting models that minimize the risk on the training dataset. Given a variably parameterized family of models, the statistical learning theory aims to identify the dependence between model complexity and model performance. The empirical risk usually decreases monotonically as the model complexity increases, and achieves its minimum when the model is rich enough to interpolate the training data, resulting in zero (or near-zero) training error. In contrast, the behaviour of the test error as a function of model complexity is far more complicated. Indeed, in this paper we show how to construct a model family for which the generalization curve can be fully controlled (away from the interpolation threshold) in both under-parameterized and over-parameterized regimes. Classical statistical learning theory supports a U-shaped curve of generalization versus model complexity . Under such a framework, the best model is found at the bottom of the U-shaped curve, which corresponds to appropriately balancing under-fitting and over-fitting the training data. From the view of the bias-variance trade-off, a higher model complexity increases the variance while decreasing the bias. A model with an appropriate level of complexity achieves a relatively low bias while still keeping the variance under control. On the other hand, a model that interpolates the training data is deemed to over-fit and tends to worsen the generalization performance due to the soaring variance.

Although classical statistical theory suggests a pattern of behavior for the generalization curve up to the interpolation threshold, it does not describe what happens beyond the interpolation threshold, commonly referred to as the over-parameterized regime. This is the exact regime where many modern machine learning models, especially deep neural networks, achieved remarkable success. Indeed, neural networks generalize well even when the models are so complex that they have the potential to interpolate all the training data points .

Modern practitioners commonly deploy deep neural networks with hundreds of millions or even billions of parameters. It has become widely accepted that large models achieve performance superior to small models that may be suggested by the classical U-shaped generalization curve . This indicates that the test error decreases again once model complexity grows beyond the interpolation threshold, resulting in the so called double-descent phenomenon described in , which has been broadly supported by empirical evidence and confirmed empirically on modern neural architectures by Nakkiran et al. . On the theoretical side, this phenomenon has been recently addressed by several works on various model settings. In particular, Belkin et al. proved the existence of double-descent phenomenon for linear regression with random feature selection and analyzed the random Fourier feature model . Mei and Montanari also studied the Fourier model and computed the asymptotic test error which captures the double-descent phenomenon. Bartlett et al. , Tsigler and Bartlett analyzed and gave explicit conditions for “benign overfitting” in linear and ridge regression, respectively. Caron and Chretien provided a finite sample analysis of the nonlinear function estimation and showed that the parameter learned through empirical risk minimization converges to the true parameter with high probability as the model complexity tends to infinity, implying the existence of double descent. Liu et al. studied the high dimensional kernel ridge regression in the under- and over-parameterized regimes and showed that the risk curve can be double descent, bell-shaped, and monotonically decreasing.

Among all the aforementioned efforts, one particularly interesting question is whether one can observe more than two descents in the generalization curve. d’Ascoli et al. empirically showed a sample-wise triple-descent phenomenon under the random Fourier feature model. Similar triple-descent was also observed for linear regression . More rigorously, Liang et al. presented an upper bound on the risk of the minimum-norm interpolation versus the data dimension in Reproducing Kernel Hilbert Spaces (RKHS), which exhibits multiple descent. However, a multiple-descent upper bound without a properly matching lower bound does not imply the existence of a multiple-descent generalization curve. In this work, we study the multiple descent phenomenon by addressing the following questions:

Can the existence of a multiple descent generalization curve be rigorously proven?

Can an arbitrary number of descents occur?

Can the generalization curve and the locations of descents be designed?

In this paper, we show that the answer to all three of these questions is yes. Further related work is presented in Section 2.

Our Contribution. We consider the linear regression model and analyze how the risk changes as the dimension of the data grows. In the linear regression setting, the data dimension is equal to the dimension of the parameter space, which reflects the model complexity. We rigorously show that the multiple descent generalization curve exists under this setting. To our best knowledge, this is the first work proving a multiple descent phenomenon.

On the one hand, we show theoretically that the generalization curve is malleable and can be constructed in an arbitrary fashion. On the other hand, we rarely observe complex generalization curves in practice, besides carefully curated constructions. Putting these facts together, we arrive at the conclusion that realistic generalization curves arise from specific interactions between properties of typical data and the inductive biases of algorithms. We should highlight that the nature of these interactions is far from being understood and should be an area of further investigations.

Related Work

Our work is directly related to the recent line of research in the theoretical understanding of the double descent and the multiple descent phenomenon . Here we briefly discuss some other work that is closely related to this paper.

In this paper we focus on the least square linear regression with no regularization. For the regularized least square regression, De Vito et al. proposed a selection procedure for the regularization parameter. Advani and Saxe analyzed the generalization of neural networks with mean squared error under the asymptotic regime where both the sample size and model complexity tend to infinity. Richards et al. proved for least square regression in the asymptotic regime that as the dimension-to-sample-size ratio $d/n$ grows, an additional peak can occur in both the variance and bias due to the covariance structure of the features. As a comparison, in this paper the sample size is fixed and the model complexity increases. Rudi and Rosasco studied kernel ridge regression and gave an upper bound on the number of the random features to reach certain risk level. Our result shows that there exists a natural setting where by manipulating the random features one can control the risk curve.

Over-Parameterization and Interpolation.

The double descent occurs when the model complexity reaches and increases beyond the interpolation threshold. Most previous works focused on proving an upper bound or optimal rate for the risk. Caponnetto and De Vito gave the optimal rate for least square ridge regression via careful selection of the regularization parameter. Belkin et al. showed that the optimal rate for risk can be achieved by a model that interpolates the training data. In a series of work on kernel regression with regularization parameter tending to zero (a.k.a. kernel ridgeless regression), Rakhlin and Zhai showed that the risk is bounded away from zero when the data dimension is fixed with respect to the sample size. Liang and Rakhlin then considered the case when $d\asymp n$, showed empirically the multiple descent phenomenon and proved a risk upper bound that can be small given favorable data and kernel assumptions. Instead of giving a bound, our paper presents an exact computation of risk in the cases of underparametrized and overparametrized linear regression, and proves the existence of the multiple descent phenomenon. Wyner et al. analyzed AdaBoost and Random Forest from the perspective of interpolation. There has also been a line of work on wide neural networks .

Sample-wise Double Descent and Non-monotonicity.

There has also been recent development beyond the model-complexity double-descent phenomenon. For example, regarding sample-wise non-monotonicity, Nakkiran et al. empirically observed the epoch-wise double-descent and sample-wise non-monotonicity for neural networks. Chen et al. and Min et al. identified and proved the sample-wise double descent under the adversarial training setting, and Javanmard et al. discovered double-descent under adversarially robust linear regression. Loog et al. showed that empirical risk minimization can lead to sample-wise non-monotonicity in the standard linear model setting under various loss functions including the absolute loss and the squared loss, which covers the range from classification to regression. We also refer the reader to their discussion of the earlier work on non-monotonicity of generalization curves. Dar et al. demonstrated the double descent curve of the generalization errors of subspace fitting problems. Fei et al. studied the risk-sample tradeoff in reinforcement learning.

Preliminaries and Problem Formulation

Distributions. Let $\cN(\mu,\sigma^{2})$ ($\mu,\sigma\in\bR$) and $\cN(\mu,\Sigma)$ ($\mu\in\bR^{n}$, $\Sigma\in\bR^{n\times n}$) denote the univariate and multivariate Gaussian distributions, respectively, where $\mu\in\bR^{n}$ and $\Sigma\in\bR^{n\times n}$ is a positive semi-definite matrix. We define a family of trimodal Gaussian mixture distributions as follows

Let $\chi^{2}(k,\lambda)$ denote the noncentral chi-squared distribution with $k$ degrees of freedom and the non-centrality parameter $\lambda$. For example, if $X_{i}\sim\cN(\mu_{i},1)$ (for $i=1,2,\dots,k$) are independent Gaussian random variables, we have $\sum_{i=1}^{k}X_{i}^{2}\sim\chi^{2}(k,\lambda)$, where $\lambda=\sum_{i=1}^{k}\mu_{i}^{2}$. We also denote by $\chi^{2}(k)$ the (central) chi-squared distribution with $k$ degrees and the $F$-distribution by $F(d_{1},d_{2})$ where $d_{1}$ and $d_{2}$ are the degrees of freedom.

Problem Setup. Let $x_{1},\dots,x_{n}\in\bR^{D}$ be column vectors that represent the training data of size $n$ and let $x_{\textnormal{test}}\in\bR^{D}$ be a column vector that represents the test data. We assume that they are all independently drawn from a distribution

where the noise $\varepsilon_{i}\sim\cN(0,\eta^{2})$. We use the same setup as in (see Equations (1) and (2) in ). Moreover, in another closely related work , if the kernel is set to the linear kernel, it is equivalent to our setup.

where $y=x^{\top}\beta+\varepsilon_{\textnormal{test}}$ and $\varepsilon_{\textnormal{test}}\sim\cN(0,\eta^{2})$. We call the term $\bE\left[(x^{\top}(A^{+}A-I)\beta)^{2}\right]$ the bias and call the term $\eta^{2}\bE\left\|(A^{\top})^{+}x\right\|^{2}$ the variance.

The next remark shows that in the underparametrized regime, the bias vanishes. The vanishing bias in the underparametrized regime is also observed by Hastie et al. and shown in their Proposition 2.

In the underparametrized regime, if $\cD$ is a continous distribution (our construction presented later satisfies this condition), the matrix $A$ has independent column almost surely. In this case, we have $A^{+}A=I$ and therefore the bias $\bE\left[(x^{\top}(A^{+}A-I)\beta)^{2}\right]$ vanishes irrespective of $\beta$. In other words, in the underparametrized regime, $L_{d}$ equals $\eta^{2}\bE\|(A^{\top})^{+}x\|^{2}$.

According to Remark 1, we have $L_{d}=\eta^{2}\bE\|(A^{\top})^{+}x\|^{2}$ in the underparametrized regime. It also holds in the overparametrized regime when $\beta=0$. Without loss of generality, we assume $\eta=1$ in the underparametrized regime (for all $\beta$). In the overparametrized regime, we also assume $\eta=1$ for the $\beta=0$ case. In this case, we have

We assume a general $\eta$ (i.e., not necessarily being $1$) in the overparametrized regime when $\beta$ is non-zero.

We would like to study the change in the loss caused by the growth in the number of features revealed. Recall $L_{d}=\bE\|(A^{\top})^{+}x\|^{2}$. Once we reveal a new feature, which adds a new row $b^{\top}$ to $A^{\top}$ and a new component $a_{1}$ to $x$, we have $L_{d+1}=\bE\left\|\begin{bmatrix}A^{\top}\\ b^{\top}\end{bmatrix}^{+}\begin{bmatrix}x\\ a_{1}\end{bmatrix}\right\|^{2}$.

Local Maximum and Multiple Descent. Throughout the paper, we say that a local maximum occurs at a dimension $d\geq 1$ if $L_{d-1}<L_{d}$ and $L_{d}>L_{d+1}$. Intuitively, a local maximum occurs if there is an increasing stage of the generalization loss, followed by a decreasing stage, as the dimension $d$ grows. Additionally, we define $L_{0}\triangleq-\infty$. If the generalization loss exhibits a single descent, based on our definition, a unique local maximum occurs at $d=1$. For a double-descent generalization curve, a local maximum occurs at two different dimensions. In general, if we observe local maxima at multiple dimensions, we say there is a multiple descent.

Underparametrized Regime

First, we present our main theorem for the underparametrized regime below, whose proof is deferred to the end of Section 4. It states that the generalization loss $L_{d}$ is always non-decreasing as $d$ grows. Moreover, it is possible to have an arbitrarily large ascent, i.e., $L_{d+1}-L_{d}>C$ for any $C>0$.

If $d<n$, we have $L_{d+1}\geq L_{d}$ irrespective of the data distribution. Moreover, for any $C>0$, there exists a distribution $\cD$ such that $L_{d+1}-L_{d}>C$.

The first part of Theorem 1 holds irrespective of the data distribution. For the second part of the theorem ( i.e., for any $C>0$ there exists a distribution such that $L_{d+1}-L_{d}>C$) to hold, one extremely simple and elegant choice of the distribution $\cD$ is a product distribution $\cD=\cD_{1}\times\dots\times\cD_{D}$ such that $x_{i,j}\stackrel{{\scriptstyle iid}}{{\sim}}\cD_{j}$ for all $1\leq i\leq n$, where $\cD_{j}$ is a Gaussian mixture $\cN^{\textnormal{mix}}_{\sigma_{j},1}$ for some $\sigma_{j}>0$. Since the second part of Theorem 1 is of independent interest, the result is summarized by Theorem 4.

In light of Remark 2, $\cD$ can be chosen to be a product distribution that consists $\cN^{\textnormal{mix}}_{\sigma_{j}}$. Note that one can simulate $\cN^{\textnormal{mix}}_{\sigma,1}$ with $\cN(0,1)$ through the inverse transform sampling. To see this, let $F_{\cN(0,1)}$ and $F_{\cN^{\textnormal{mix}}_{\sigma,1}}$ be the cdf of $\cN(0,1)$ and $\cN^{\textnormal{mix}}_{\sigma,1}$, respectively. If $X\sim\cN(0,1)$, we have $F_{\cN(0,1)}(X)\sim\operatorname{Unif}((0,1))$ and therefore $\varphi_{\sigma}(X)\triangleq F^{-1}_{\cN^{\textnormal{mix}}_{\sigma,1}}(F_{\cN(0,1)}(X))\sim\cN^{\textnormal{mix}}_{\sigma,1}$. In fact, we can use a multivariate Gaussian $\cD^{\prime}=\cN(0,I_{D\times D})$ and a sequence of non-linear kernels $k^{[1:d]}(x,x^{\prime})\triangleq\langle\phi^{[1:d]}(x),\phi^{[1:d]}(x^{\prime})\rangle$, where the feature map is $\phi^{[1:d]}(x)\triangleq[\phi_{1}(x_{1}),\phi_{2}(x_{2}),\dots,\phi_{d}(x_{d})]^{\top}\in\bR^{d}$. Here is a simple rule for defining $\phi_{j}$: if $\cD_{j}=\cN^{\textnormal{mix}}_{\sigma_{j}}$, we set $\phi_{j}$ to $\varphi_{\sigma_{j}}$. Thus, the problem becomes a kernel regression problem on the standard Gaussian data.

The first part of Theorem 1, which says that $L_{d}$ is increasing (or more precisely, non-decreasing), agrees with Figure 1 of and Proposition 2 of . In , they proved that the risk increases with $\gamma=d/n$. Note that, at first glance, Theorem 1 may look counterintuitive since it does not obey the classical U-shaped generalization curve. However, we would like to emphasize that the U-shaped curve does not always occur. In Figure 1 and Proposition 2 of these two papers respectively, there is no U-shaped curve. The intuition behind Theorem 1 is that in the underparametrized setting, the bias is always zero and as $d$ approaches $n$, the variance keeps increasing.

Coming to the second part of Theorem 1, we now discuss how we will construct such a distribution $\cD$ inductively to satisfy $L_{d+1}-L_{d}>C$. We fix $d$. Again, denote the first $d$ features of $x_{\textnormal{test}}$ by $x\triangleq x_{\textnormal{test}}[1:d]$. Let us add an additional component to the training data $x_{1}[1:d],\dots,x_{n}[1:d]$ and test data $x$ so that the dimension $d$ is incremented by 1. Let $b_{i}\in\bR$ denote the additional component that we add to the vector $x_{i}$ (so that the new vector is given as $[x_{i}[1:d]^{\top},b_{i}]^{\top}$. Similarly, let $a_{1}\in\bR$ denote the additional component that we add to the test vector $x$. We form the column vector $b=[b_{1},\dots,b_{n}]^{\top}\in\bR^{n}$ that collects all additional components that we add to the training data.

We consider the change in the generalization loss as follows

Note that the components $b_{1},\dots,b_{n},a_{1}$ are i.i.d. The proof of Theorem 1 starts with Lemma 2 which relates the pseudo-inverse of $[A,b]^{\top}$ to that of $A^{\top}$. In this way, we can decompose $\left\|\begin{bmatrix}A^{\top}\\ b^{\top}\end{bmatrix}^{+}\begin{bmatrix}x\\ a_{1}\end{bmatrix}\right\|^{2}$ into multiple terms for further careful analysis in the proofs hereinafter.

Let $A\in\bR^{n\times d}$ and $0\neq b\in\bR^{n\times 1}$, where $n\geq d+1$. Additionally, let $P=AA^{+}$ and $Q=bb^{+}=\frac{bb^{\top}}{\|b\|^{2}}$, and define $z\triangleq\frac{b^{\top}(I-P)b}{\|b\|^{2}}$. If $z\neq 0$ and the columnwise partitioned matrix $[A,b]$ has linearly independent columns, we have

In our construction of $\cD$, the components $\cD_{j}$ are all continuous distributions. The matrix $I-P$ is an orthogonal projection matrix and therefore $\operatorname{rank}(I-P)=n-d$. As a result, it holds almost surely that $b\neq 0$, $z\neq 0$, and $[A,b]$ has linearly independent columns. Thus the assumptions of Lemma 2 are satisfied almost surely. In the sequel, we assume that these assumptions are always fulfilled.

Theorem 3 guarantees that if $L_{d}=\bE\left\|(A^{+})^{\top}x\right\|^{2}$ is finite and the $(d+1)$-th features $b_{1},\dots,b_{n},a_{1}$ are i.i.d. sampled from $\cN(0,1)$ or $\cN^{\textnormal{mix}}_{\sigma,1}$, $L_{d+1}=\bE\left\|\begin{bmatrix}A^{\top}\\ b^{\top}\end{bmatrix}^{+}\begin{bmatrix}x\\ a_{1}\end{bmatrix}\right\|^{2}$ is also finite.

Let $z$ be as defined in Lemma 2. If $b_{1},\dots,b_{n},a_{1}$ are i.i.d. and follow a distribution with mean zero, conditioned on $A$ and $x$, we have

In particular, if $d+2<n$ and $b_{1},\dots,b_{n},a_{1}\stackrel{{\scriptstyle iid}}{{\sim}}\cN(0,1)$, conditioned on $A$ and $x$, we have

If $d+2<n$ and $b_{1},\dots,b_{n},a_{1}\stackrel{{\scriptstyle iid}}{{\sim}}\cN^{\textnormal{mix}}_{\sigma,1}$, conditioned on $A$ and $x$, we have

Using Theorem 3, we can show inductively (on $d$) that $L_{d}$ is finite for every $d$. Provided that we are able to guarantee finite $L_{1}$, Theorem 3 implies that $L_{d}$ is finite for every $d$ if the components are always sampled from $\cN(0,1)$ or $\cN^{\textnormal{mix}}_{\sigma,1}$.

Making a large $L_{d}$ can be achieved by adding an entry sampled from $\cN^{\textnormal{mix}}_{\sigma,1}$ when the data dimension increases from $d-1$ to $d$ in the previous step. Theorem 4 shows that adding a $\cN^{\textnormal{mix}}_{\sigma,1}$ feature can increase the loss by arbitrary amount, which in turn implies the second part of Theorem 1.

For any $C>0$ and $\bE\left\|(A^{+})^{\top}x\right\|^{2}<+\infty$, there exists a $\sigma>0$ such that if $b_{1},\dots,b_{n},a_{1}\stackrel{{\scriptstyle iid}}{{\sim}}\cN^{\textnormal{mix}}_{\sigma,1}$, we have

We follow the notation convention in (3):

Recall $d<n$ and the matrix $B^{\prime}\triangleq\begin{bmatrix}A^{\top}\\ b^{\top}\end{bmatrix}$ is of size $(d+1)\times n$. Both matrices $B^{\prime}$ and $B\triangleq A^{\top}$ are fat matrices. As a result, if $x^{\prime}\triangleq\begin{bmatrix}x\\ a_{1}\end{bmatrix}$, we have

Since $\{z\mid B^{\prime}z=x^{\prime}\}\subseteq\{z\mid Bz=x\}$, we get $\|B^{\prime+}x^{\prime}\|^{2}\geq\|B^{+}x\|^{2}$. Therefore, we obtain $L_{d+1}\geq L_{d}$. The second part follows from Theorem 4. ∎

Overparametrized Regime

In this section, we study the multiple decent phenomenon in the overparametrized regime. Note that as stated in Section 3, we consider the minimum-norm solution here. We first consider the case where the model $\beta=0$ and $L_{d}$ is as defined in (2). Then we discuss the setting $\beta\neq 0$.

As stated in the following theorem, we require $d\geq n+8$. This is merely a technical requirement and we can still say that $d$ starts at roughly the same order as $n$. In other words, the result covers almost the entire spectrum of the overparametrized regime.

Let $n<D-9$. Given any sequence $\Delta_{n+8},\Delta_{n+9},\dots$, $\Delta_{D-1}$ where $\Delta_{d}\in\{{\uparrow},{\downarrow}\}$, there exists a distribution $\cD$ such that for every $n+8\leq d\leq D-1$, we have

In Theorem 5, the sequence $\Delta_{n+8}$, $\Delta_{n+9}$, $\cdots$, $\Delta_{D-1}$ is just used to specify the increasing/decreasing behavior of the $L_{d}$ sequence for $d>n+8$. Compared to Theorem 1 for the underparametrized regime, where $L_{d}$ always increases, Theorem 5 indicates that one is able to fully control both ascents and descents in the overparametrized regime. Fig. 2 is an illustration.

We now present tools for proving Theorem 5. Lemma 6 gives the pseudo-inverse of $A$ when $d>n$.

Let $A\in\bR^{n\times d}$ and $b\in\bR^{n\times 1}$, where $n\leq d$. Assume that matrix $A$ and the columnwise partitioned matrix $B\triangleq[A,b]$ have linearly independent rows. Let $G\triangleq(AA^{\top})^{-1}\in\bR^{n\times n}$ and $u\triangleq\frac{b^{\top}G}{1+b^{\top}Gb}\in\bR^{1\times n}$. We have

Lemma 7 establishes finite expectation for several random variables. These finite expectation results are necessary for Theorem 8 and Theorem 9 to hold. Technically, they are the dominating random variables needed in Lebesgue’s dominated convergence theorem. Lemma 7 indicates that to guarantee these finite expectations, it suffices to set the first $n+8$ distributions to the standard normal distribution and then set $\cD_{n+8},\dots,\cD_{D}$ to either a Gaussian or a Gaussian mixture distribution. In fact, in Theorem 8 and Theorem 9, we always add a Gaussian distribution or a Gaussian mixture.

Let $\cD=\cD_{1}\times\cdots\times\cD_{D}$ be a product distribution where

$\cD_{d}=\cN(0,1)$ if $d=1,\dots,n+8$; and

$\cD_{d}$ is either $\cN(0,\sigma_{d}^{2})$ or $\cN^{\textnormal{mix}}_{\sigma_{d},\mu_{d}}$ for $d>n+8$.

Let $\cD_{[1:d]}$ denote $\cD_{1}\times\cdots\times\cD_{d}$. Assume that every row of $A\in\bR^{n\times d}$ and $x\in\bR^{d\times 1}$ are i.i.d. and follow $\cD_{[1:d]}$. For any $d$ such that $n+8\leq d\leq D$, all of the followings hold:

Theorems 8 and 9 are the key technical results for constructing multiple descent in the overparametrized regime. One can create a descent ($L_{d+1}<L_{d}$) by adding a Gaussian feature (Theorem 8) and create an ascent ($L_{d+1}>L_{d}$) by adding a Gaussian mixture feature (Theorem 9).

If $\bE[\|(A^{\top}A)^{+}x\|^{2}]>0$ and all equations in (4) hold, there exists $\sigma>0$ such that if $a_{1},b_{1},\dots,b_{n}\stackrel{{\scriptstyle iid}}{{\sim}}\cN(0,\sigma^{2})$, we have

Theorem 9 shows that adding a Gaussian mixture feature can make $L_{d+1}>L_{d}$.

Assume $\bE\|(A^{+})^{\top}x\|^{2}<+\infty$. For any $C>0$, there exist $\mu$, $\sigma>0$ such that if $a_{1},b_{1},\dots,b_{n}\stackrel{{\scriptstyle iid}}{{\sim}}\cN^{\textnormal{mix}}_{\sigma,\mu}$, we have

The proof of Theorem 5 immediately follows from Theorem 8 and Theorem 9.

We construct the product distribution $\cD=\prod_{d=1}^{D}\cD_{d}$. We set $\cD_{d}=\cN(0,1)$ for $d=1,\dots,n+8$. For $n+8<d\leq D$, $\cD_{d}$ is either $\cN(0,\sigma_{d}^{2})$ or $\cN^{\textnormal{mix}}_{\sigma_{d},\mu_{d}}$ depending on $\Delta_{d}$ being either $\downarrow$ or $\uparrow$.

First we show that for each step $d$, the assumption $\bE[\|(A^{\top}A)^{+}x\|^{2}]>0$ of Theorem 8 is satisfied. If $\bE[\|(A^{\top}A)^{+}x\|^{2}]=0$, we know that $(A^{\top}A)^{+}x=0$ almost surely. Since $\cD$ is a continuous distribution, the matrix $A$ has full row rank almost surely. Therefore, $\operatorname{rank}((A^{\top}A)^{+})=\operatorname{rank}(A^{\top}A)=n$ almost surely. Thus $\dim\ker(A^{\top}A)^{+}=d-n\leq d-1$ almost surely, which implies $x\notin\ker(A^{\top}A)^{+}$. In other words, $(A^{\top}A)^{+}x\neq 0$ almost surely. We reach a contradiction. Moreover, by Lemma 7, the assumption $\bE\|(A^{+})^{\top}x\|^{2}<+\infty$ of Theorem 9 is also satisfied.

If $\Delta_{d-1}={\downarrow}$, by Theorem 8, there exists $\sigma_{d}>0$ such that if $\cD_{d}=\cN(0,\sigma^{2}_{d})$, then $L_{d}<L_{d-1}$. Similarly if $\Delta_{d-1}={\uparrow}$, by Theorem 9, there exists $\sigma_{d}$ and $\mu_{d}$ such that $\cD_{d}=\cN^{\textnormal{mix}}_{\sigma_{d},\mu_{d}}$ guarantees $L_{d}>L_{d-1}$.

Gaussian $\beta$ setting. In what follows, we study the case where the model $\beta$ is non-zero. In particular, we consider a setting where each entry of $\beta$ is i.i.d. $\cN(0,\rho^{2})$. Recalling (1), define the biases

where $\beta\sim\cN(0,\rho^{2}I_{d})$ and $\beta_{1}\sim\cN(0,\rho^{2})$. The second term in $L^{\textnormal{exp}}_{d}$ and $L^{\textnormal{exp}}_{d+1}$ is the variance term. Note that $L^{\textnormal{exp}}_{d}$ is the expected value of $L_{d}$ in (1) and averages over $\beta$. Theorem 10 shows that one can add a Gaussian mixture feature in order to make $L^{\textnormal{exp}}_{d+1}>L^{\textnormal{exp}}_{d}$, and add a Gaussian feature in order to make $L^{\textnormal{exp}}_{d+1}<L^{\textnormal{exp}}_{d}$.

Let $a_{1},\beta_{1}\in\bR$, $x\in\bR^{d\times 1}$, $\beta\in\bR^{d\times 1}$, $A\in\bR^{n\times d}$ and $b\in\bR^{n\times 1}$, where $n\leq d$. Assume that $x,a_{1},\beta_{1},\beta,A,b$ are jointly independent, $[\beta^{\top},\beta_{1}]^{\top}\sim\cN(0,\rho^{2}I_{d+1})$. Moreover, assume that the matrix $[A,b]$ has linearly independent rows almost surely. The following statements hold:

If $a_{1},b_{1},\dots,b_{n}\stackrel{{\scriptstyle iid}}{{\sim}}\cN^{\textnormal{mix}}_{\sigma,\mu}$, for any $C>0$, there exist $\mu,\sigma$ such that $L^{\textnormal{exp}}_{d+1}-L^{\textnormal{exp}}_{d}>C$.

If $a_{1},b_{1},\dots,b_{n}\stackrel{{\scriptstyle iid}}{{\sim}}\cN(0,\sigma^{2})$, there exists $\sigma>0$ such that for all

we have $L^{\textnormal{exp}}_{d+1}<L^{\textnormal{exp}}_{d}$.

Theorem 10 indicates that for $\beta$ obeying a normal distribution, one can still construct a generalization curve as desired by adding a Gaussian or Gaussian mixture feature properly. We make this construction explicit for any desired generalization curve in (the proof of) Theorem 11. Similar to the construction in the underparametrized regime (for all $\beta$) and overparametrization regime (for $\beta=0$), the distribution $\cD$ can be made a product distribution.

Let $n<D-9$. Given any sequence $\Delta_{n+8},\Delta_{n+9},\dots$, $\Delta_{D-1}$ where $\Delta_{d}\in\{{\uparrow},{\downarrow}\}$, there exists $\rho>0$ and a distribution $\cD$ such that for $\beta\sim\cN(0,\rho^{2})$ and every $n+8\leq d\leq D-1$, we have

Define the design matrix $A_{d}\triangleq[x_{1}[1:d],\dots,x_{n}[1:d]]^{\top}\in\bR^{n\times d}$. Similar to the proof of Theorem 5, we construct the product distribution $\cD=\prod_{d=1}^{D}\cD_{d}$. We set $\cD_{d}=\cN(0,1)$ for $d=1,\dots,n+8$. For $n+8<d\leq D$, $\cD_{d}$ is either $\cN(0,\sigma_{d}^{2})$ or $\cN^{\textnormal{mix}}_{\sigma_{d},\mu_{d}}$ depending on $\Delta_{d}$ being either $\downarrow$ or $\uparrow$.

If $\Delta_{d-1}={\uparrow}$, by Theorem 10, there exists $\sigma_{d}$ and $\mu_{d}$ such that $\cD_{d}=\cN^{\textnormal{mix}}_{\sigma_{d},\mu_{d}}$ guarantees $L^{\textnormal{exp}}_{d}>L^{\textnormal{exp}}_{d-1}$. If $\Delta_{d-1}={\downarrow}$, define

By Theorem 10, there exists $\sigma_{d}>0$ such that if $\rho\leq\rho_{d}$ and $\cD_{d}=\cN(0,\sigma^{2}_{d})$, then $L^{\textnormal{exp}}_{d}<L^{\textnormal{exp}}_{d-1}$. We take

Conclusion

Our work proves that the expected risk of linear regression can manifest multiple descents when the number of features increases and sample size is fixed. This is carried out through an algorithmic construction of a feature-revealing process where the newly revealed feature follows either a Gaussian distribution or a Gaussian mixture distribution. Notably, the construction also enables us to control local maxima in the underparametrized regime and control ascents/descents freely in the overparametrized regime. Overall, this allows us to design the generalization curve away from the interpolation threshold.

We believe that our analysis of linear regression in this paper is a good starting point for explaining non-monotonic generalization curves observed in machine learning studies. Extending these results to more complex problem setups would be a meaningful future direction.

Funding Transparency Statement

LC: Funding in direct support of this work: postdoctoral research fellowship by the Simons Institute for the Theory of Computing, University of California, Berkeley, and Google PhD Fellowship by Google. Additional revenues related to this work: internships at Google.

MB acknowledges support from NSF IIS-1815697, and the support of the NSF and the Simons Foundation for the Collaboration on the Theoretical Foundations of Deep Learning through awards DMS-2031883 and #814639.

AK: Funding in direct support of this work: NSF (IIS-1845032) and ONR (N00014-19-1-2406).

References

Appendix A Almost Sure Convergence of Sequence of Normal Random Variables

In this paper, we need a sequence of random variables $\{X_{n}\}_{n\geq 1}$ such that $X_{n}\sim\cN(0,\sigma_{n}^{2})$, $\lim_{n\to+\infty}\sigma_{n}=0$, and $X_{n}\to 0$ almost surely. The following lemma shows the existence of such a sequence.

There exist a sequence of random variables $\{X_{n}\}_{n\geq 1}$ such that $X_{n}\sim\cN(0,\sigma_{n}^{2})$, $\lim_{n\to+\infty}\sigma_{n}=0$, and $X_{n}\to 0$ almost surely.

Let $\sigma_{n}=1/n^{2}$ and $X_{n}\sim\cN(0,\sigma_{n}^{2})$. Define the event $E_{n}\triangleq\{|X_{n}|>\varepsilon\}$. We have

By the Borel–Cantelli lemma, we have $\bP(\limsup_{n\to+\infty}E_{n})=0$, which implies that $X_{n}\to 0$ almost surely. ∎

Appendix B Proofs for Underparametrized Regime

Define $r\triangleq A^{\top}b\in\bR^{d}$. Since $A$ has linearly independent columns, the Gram matrix $G=A^{\top}A$ is non-singular. The Sherman-Morrison formula gives

where we use the facts $r=A^{\top}b$ and $AG^{-1}=(A^{+})^{\top}$ in the last equality. Therefore, we deduce

Therefore, we obtain the desired expression.

B.2 Proof of Theorem 3

First, we rewrite the expression as follows

where $P,Q,z$ are defined in Lemma 2. Since $a_{1}$ has mean 0 and is independent of other random variables, so that the cross term vanishes under expectation over $b$ and $a_{1}$:

where $\langle\cdot,\cdot\rangle$ denotes the inner product. Therefore taking the expectation of (6) over $b$ and $a_{1}$ yields

We simplify the third term. Recall that $I-P=I-AA^{+}$ is an orthogonal projection matrix and thus idempotent

We consider the first and second terms. We write $v=(A^{+})^{\top}x$ and define $z=\frac{b^{\top}(I-P)b}{\|b\|^{2}}$. The sum of the first and second terms equals

The rank of $M$ is at most $2$. To see this, we re-write $M$ in the following way

Notice that $\operatorname{rank}(M_{1})\leq\operatorname{rank}(Q)$, $\operatorname{rank}(M_{2})\leq\operatorname{rank}(Q)$, and $\operatorname{rank}(Q)=1$.

It follows that $\operatorname{rank}(M)\leq\operatorname{rank}(M_{1})+\operatorname{rank}(M_{2})=2$. The matrix $M$ has at least $n-2$ zero eigenvalues. We claim that $M$ has two non-zero eigenvalues and they are $1-1/z<0$ and $1$.

thus $PQ$ has a unique non-zero eigenvalue $1-z$. Let $u\neq 0$ denote the corresponding eigenvector such that $PQu=(1-z)u$. Since $u\in\operatorname{im}P$ and $P$ is a projection, we have $Pu=u$. Therefore we can verify that

To show that the other non-zero eigenvalue of $M$ is $1$, we compute the trace of $M$

where we use the fact that $\operatorname{tr}(Q)=1$, $\operatorname{tr}(PQ)=1-z$,

We have shown that $M$ has eigenvalue $1-1/z$ and $M$ has at most two non-zero eigenvalues. Therefore, the other non-zero eigenvalue is $\operatorname{tr}(M)-(1-1/z)=1$.

We are now in a position to upper bound (13) as follows:

Putting all three terms of the change in the dimension-normalized generalization loss yields

For $b_{1},\dots,b_{n},a_{1}\stackrel{{\scriptstyle iid}}{{\sim}}\cN(0,1)$, we have $\bE[a_{1}^{2}]=1$. Moreover, $b^{\top}(I-P)b$ follows $\chi^{2}(n-d)$ a distribution. Thus $\frac{1}{b^{\top}(I-P)b}$ follows an inverse-chi-squared distribution with mean $\frac{1}{n-d-2}$. Therefore the expectation $\bE[\frac{a_{1}^{2}}{b^{\top}(I-P)b}]=\frac{1}{n-d-2}$.

Notice that $1/z$ follows a $1+\frac{d}{n-d}F(d,n-d)$ distribution and thus $\bE[1/z]=1+\frac{d}{n-d-2}$.

For $b_{1},\dots,b_{n},a_{1}\stackrel{{\scriptstyle iid}}{{\sim}}\cN^{\textnormal{mix}}_{\sigma,1}$, we need the following lemma.

Assume $d$, $n>d+2$ and $P$ are fixed, where $P\in\bR^{n\times n}$ is an orthogonal projection matrix whose rank is $d$. Define $z\triangleq\frac{b^{\top}(I-P)b}{\|b\|^{2}}$, where $b=[b_{1},\dots,b_{n}]^{\top}\in\bR^{n}$. If $a_{1},\ b_{1},\cdots,\ b_{n}\stackrel{{\scriptstyle iid}}{{\sim}}\cN^{\textnormal{mix}}_{\sigma,1}$, we have $\bE[1/z]\leq\frac{n-2+\sqrt{d}}{n-d-2}$ and $\bE[a_{1}^{2}/b^{\top}(I-P)b]\leq\frac{2/(3\sigma^{2})+1}{n-d-2}$.

B.3 Proof of Lemma 13

Lemma 14 shows that a noncentral $\chi^{2}$ distribution first-order stochastically dominates a central $\chi^{2}$ distribution of the same degree of freedom. It will be needed in the proof of Lemma 13.

Assume that random variables $X\sim\chi^{2}(k,\lambda)$ and $Y\sim\chi^{2}(k)$, where $\lambda>0$. For any $c>0$, we have

In other words, the random variable $X$ (first-order) stochastically dominates $Y$.

Let $Y_{1},X_{2},\dots,X_{k}\stackrel{{\scriptstyle iid}}{{\sim}}\cN(0,1)$ and $X_{1}\sim\cN(\sqrt{\lambda},1)$ and all these random variables are jointly independent. Then $X^{\prime}\triangleq\sum_{i=1}^{k}X_{i}^{2}\sim\chi^{2}(k,\lambda)$ and $Y^{\prime}\triangleq Y_{1}^{2}+\sum_{i=2}^{k}X_{i}^{2}\sim\chi^{2}(k)$.

Since $b_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\cN^{\textnormal{mix}}_{\sigma,1}$, we can rewrite $b=u+w$ where $w\sim\cN(0,\sigma^{2}I_{n})$ and the entries of $u$ satisfy $u_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\operatorname{Unif}(\{-1,0,1\})$. Furthermore, $u$ and $w$ are independent. Similarly, we can write $a_{1}=\hat{u}+\hat{w}$, where $\hat{u}\sim\operatorname{Unif}(\{-1,0,1\})$ and $\hat{w}\sim\cN(0,\sigma^{2})$ are independent. To bound $\bE[a_{1}^{2}]$, we have

Since $P$ is an orthogonal projection, there exists an orthogonal transformation $O$ depending only on $P$ such that

and that these two quantities are independent. It follows that

Putting the numerator and denominator together yields

B.4 Proof of Theorem 4

We start from (12). Taking expectation over all random variables gives

Our strategy is to choose $\sigma$ so that $\bE\left[\frac{a_{1}^{2}}{\sum_{i=1}^{n}b_{i}^{2}}\right]$ is sufficiently large. This is indeed possible as we immediately show. Define independent random variables $u\sim\operatorname{Unif}(\{-1,0,1\})$ and $w\sim\cN(0,\sigma^{2})$. Since $a_{1}$ has the same distribution as $u+w$, we have

Appendix C Proofs for Overparametrized Regime

Since $A$ and $B$ have full row rank, $(AA^{\top})^{-1}$ and $(BB^{\top})^{-1}$ exist. Therefore we have

Transposing the above equation yields to the promised equation.

C.2 Proof of Lemma 7

First note that by Cauchy-Schwarz inequality, it suffices to show there exists $\cD$ such that $\bE[\lambda^{4}_{\textnormal{max}}(G)]<+\infty$ and $\bE\|v\|^{4}<+\infty$.

We define $A_{d}\in\bR^{n\times d}$ to be the submatrix of $A$ that consists of all $n$ rows and first $d$ columns. Denote

We will prove $\bE[\lambda_{\textnormal{max}}^{4}(G)]<+\infty$ by induction.

The base step is $d=n+8$. Recall $\cD_{[1:n+8]}=\cN(0,I_{n+8})$. We first show $\bE[\lambda_{\textnormal{max}}(G_{n+8})]^{4}<+\infty$. Note that since $G_{n+8}$ is almost surely positive definite,

By our choice of $\cD_{[1:n+8]}$, the matrix $(A_{n+8}A_{n+8}^{\top})^{-1}$ is an inverse Wishart matrix of size $n\times n$ with $(n+8)$ degrees of freedom, and thus has finite fourth moment (see, for example, Theorem 4.1 in ). It then follows that

For the inductive step, assume $\bE[\lambda_{\textnormal{max}}(G_{d})]^{4}<+\infty$ for some $d\geq n+8$. We claim that

under the Loewner order, where $b\in\bR^{n\times 1}$ is the $(d+1)$-th column of $A$. Therefore, we have

and by induction, we conclude that $\bE[\lambda_{\textnormal{max}}^{4}(G)]<+\infty$ for all $d\geq n+8$.

Now we proceed to show $\bE\|v\|^{4}<+\infty$. We have

where the last equality uses the fact that $A^{\top}(AA^{\top})^{-2}A$ is positive semidefinite. Moreover, we deduce

Using the fact that $A_{d}A_{d}^{\top}\preccurlyeq A_{d+1}A_{d+1}^{\top}$ established above, induction gives

where again we use that fact that inverse Wishart matrix $\left(A_{n+8}A_{n+8}^{\top}\right)^{-1}$ has finite second moment.

Next, we demonstrate $\bE\|x\|^{4}<+\infty$. Recall that every $\cD_{i}$ is either a Gaussian or a Gaussian mixture distribution. Therefore, every entry of $x$ has a subgaussian tail, and thus $\bE\|x\|^{4}<+\infty$. Together with (14) and the fact that $x$ and $A$ are independent, we conclude that

C.3 Proof of Theorem 8

The randomness comes from $A,x,a_{1}$ and $b$. We first condition on $A$ and $x$ being fixed.

Let $G\triangleq(AA^{\top})^{-1}\in\bR^{n\times n}$ and $u\triangleq\frac{b^{\top}G}{1+b^{\top}Gb}\in\bR^{1\times n}$. Define

We compute the left-hand side but take the expectation over only $a_{1}$ for the moment

Let us first consider the first and third terms of the above equation:

Write $G=V\Lambda V^{\top}$, where $\Lambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{n})\in\bR^{n\times n}$ is a diagonal matrix ($\lambda_{i}>0$) and $V\in\bR^{n\times n}$ is an orthogonal matrix. Recall $b\sim\cN(0,\sigma^{2}I_{n})$. Therefore $w\triangleq V^{\top}b\sim\cN(0,\sigma^{2}I_{n})$. Taking the expectation over $b$, we have

Let $R\triangleq\bE_{w}\left[\frac{ww^{\top}\Lambda+\Lambda ww^{\top}}{1+w^{\top}\Lambda w}\right]$. We have

Notice that for any $w$ and $j$, it has the same distribution if we replace $w_{j}$ by $-w_{j}$. As a result,

Thus the matrix $R$ is a diagonal matrix and

Moreover, by the monotone convergence theorem, we deduce

Again, by the monotone convergence theorem, we have

We apply a similar method to the term $\frac{\|Gb\|^{2}}{r^{2}}$. We deduce

where $\bE[\operatorname{tr}(G^{2})]\leq n\bE[\lambda^{2}_{\textnormal{max}}((AA^{\top})^{-1})]<{}+\infty$.

Putting all three terms together, we have as $\sigma\to 0^{+}$

Therefore, there exists $\sigma>0$ such that $L_{d+1}-L_{d}<0$.

C.4 Proof of Theorem 9

Again we first condition on $A$ and $x$ being fixed. Let $G\triangleq(AA^{\top})^{-1}\in\bR^{n\times n}$ and $u\triangleq\frac{b^{\top}G}{1+b^{\top}Gb}\in\bR^{1\times n}$ as defined in Lemma 6. We also define the following variables:

We compute $L_{d+1}-L_{d}$ but take the expectation over only $a_{1}$ for the moment

Our strategy is to make $\bE[a_{1}^{2}\frac{\|Gb\|^{2}}{r^{2}}]$ arbitrarily large. To this end, by the independence of $a_{1}$ and $b$ we have

By definition of $\cN^{\textnormal{mix}}_{\sigma,\mu}$, with probability $2/3$, $a_{1}$ is sampled from either $\cN(\mu,\sigma^{2})$ or $\cN(-\mu,\sigma^{2})$, which implies $\bE[a_{1}^{2}]\geq\frac{1}{3}\mu^{2}$. For each $b_{i}$, we have

Also note that $G$ is positive definite. It follows that

where we switch the order of expectation and limit using the monotone convergence theorem. Taking full expectation over $A,x,b$ and $a_{1}$ of (15) and using the assumption that $\bE\|v\|^{2}<+\infty$ we have

C.5 Proof of Theorem 10

If we define $G\triangleq(AA^{\top})^{-1}\in\bR^{n\times n}$ and $u\triangleq\frac{b^{\top}G}{1+b^{\top}Gb}\in\bR^{1\times n}$, Lemma 6 implies

We obtain the expression for $\mathcal{E}_{d+1}$:

If $a_{1},b_{1},\dots,b_{n}\stackrel{{\scriptstyle iid}}{{\sim}}\cN^{\textnormal{mix}}_{\sigma,\mu}$ or $a_{1},b_{1},\dots,b_{n}\stackrel{{\scriptstyle iid}}{{\sim}}\cN(0,\sigma^{2})$, it holds that $\bE[a_{1}]=0\in\bR$, $\bE[x]=0\in\bR^{d}$, and $\bE[b]=0\in\bR^{n\times 1}$. Therefore we have

where the second equality is because $\beta$ is independent from the remaining random variables and the third step is because of $\beta\sim\cN(0,\rho^{2}I)$. Recalling that $w=A^{+}b$ and $A^{+}AA^{+}=A^{+}$, we have

If $a_{1},b_{1},\dots,b_{n}\stackrel{{\scriptstyle iid}}{{\sim}}\cN^{\textnormal{mix}}_{\sigma,\mu}$, Theorem 9 implies that for any $C>0$, there exist $\mu,\sigma$ such that