Multiple Descent: Design Your Own Generalization Curve
Year: 2021

Abstract:
This paper explores the generalization loss of linear regression in variably parameterized families of models, both under-parameterized and over-parameterized. We show that the generalization curve can have an arbitrary number of peaks, and moreover, locations of those peaks can be explicitly controlled. Our results highlight the fact that both classical U-shaped generalization curve and the recently observed double descent curve are not intrinsic properties of the model family. Instead, their emergence is due to the interaction between the properties of the data and the inductive biases of learning algorithms.

heading: 1 Introduction
text: 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 [31, 33]. 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,
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
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 [61, 10, 32, 34].
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 [13, 38, 55, 35, 36]. 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 [9], which has been broadly supported by empirical evidence [49, 48, 29, 30] and confirmed empirically on modern neural architectures by Nakkiran et al. [46]. On the theoretical side, this phenomenon has been recently addressed by several works on various model settings. In particular, Belkin et al. [11] proved the existence of double-descent phenomenon for linear regression with random feature selection and analyzed the random Fourier feature model [50]. Mei and Montanari [44] also studied the Fourier model and computed the asymptotic test error which captures the double-descent phenomenon. Bartlett et al. [8], Tsigler and Bartlett [56] analyzed and gave explicit conditions for “benign overfitting” in linear and ridge regression, respectively. Caron and Chretien [16] 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. [42] 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. [21] empirically showed a samplewise triple-descent phenomenon under the random Fourier feature model. Similar triple-descent was also observed for linear regression [47]. More rigorously, Liang et al. [41] 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.
Our analysis considers both the underparametrized and overparametrized regimes. In the overparametrized regime, we show that one can control where a descent or an ascent occurs in the generalization curve. This is realized through our algorithmic construction of a feature-revealing process. To be more specific, we assume that the data is in RD, where D can be arbitrarily large or even essentially infinite. We view each dimension of the data as a feature. We consider a linear regression problem restricted on the first d features, where d < D. New features are revealed by increasing the dimension of the data. We then show that by specifying the distribution of the newly revealed feature to be either a standard Gaussian or a Gaussian mixture, one can determine where an ascent or a descent occurs. In order to create an ascent when a new feature is revealed, it is sufficient that the feature follows a Gaussian mixture distribution. In order to have a descent, it is sufficient that the new feature follows a standard Gaussian distribution. Therefore, in the overparametrized regime, we can fully control the occurrence of a descent and an ascent. As a comparison, in the underparametrized regime, the generalization loss always increases regardless of the feature distribution. Generally speaking, we show that we are able to design the generalization curve.
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.

heading: 2 Related Work
text: Our work is directly related to the recent line of research in the theoretical understanding of the double descent [11, 34, 60, 44] and the multiple descent phenomenon [41, 39]. Here we briefly discuss some other work that is closely related to this paper.
Least Square Regression. In this paper we focus on the least square linear regression with no regularization. For the regularized least square regression, De Vito et al. [22] proposed a selection procedure for the regularization parameter. Advani and Saxe [1] 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. [52] 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 [53] 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 [15] gave the optimal rate for least square ridge regression via careful selection of the regularization parameter. Belkin et al. [12] 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 [51] showed that the risk is bounded away from zero when the data dimension is fixed with respect to the sample size. Liang and Rakhlin [40] then considered the case when d 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. [59] analyzed AdaBoost and Random Forest from the perspective of interpolation. There has also been a line of work on wide neural networks [4–6, 23, 3, 58, 14, 2, 18, 62, 54].
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. [46] empirically observed the epoch-wise double-descent and sample-wise non-monotonicity for neural networks. Chen et al. [19] and Min et al. [45] identified and proved the sample-wise double descent under the adversarial training setting, and Javanmard et al. [37] discovered double-descent under adversarially robust linear regression. Loog et al. [43] 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. [20] demonstrated the double descent curve of the generalization errors of subspace fitting problems. Fei et al. [28] studied the risk-sample tradeoff in reinforcement learning.

heading: 3 Preliminaries and Problem Formulation
text: Notation. For x ∈ RD and d ≤ D, we let x[1 : d] ∈ Rd denote a d-dimensional vector with x[1 : d]i = xi for all 1 ≤ i ≤ d. For a matrix A ∈ Rn×d, we denote its Moore-Penrose
pseudoinverse by A+ ∈ Rd×n and denote its spectral norm by ‖A‖ , supx6=0 ‖Ax‖2 ‖x‖2 , where ‖ · ‖2 is the Euclidean norm for vectors. If v is a vector, its spectral norm ‖v‖ agrees with the Euclidean norm ‖v‖2. Therefore, we write ‖v‖ for ‖v‖2 to simplify the notation. We use the big O notation O and write variables in the subscript of O if the implicit constant depends on them. For example, On,d,σ(1) is a constant that only depends on n, d, and σ. If f(σ) and g(σ) are functions of σ, write f(σ) ∼ g(σ) if lim f(σ)g(σ) = 1. It will be given in the context how we take the limit.
Distributions. Let N (µ, σ2) (µ, σ ∈ R) and N (µ,Σ) (µ ∈ Rn, Σ ∈ Rn×n) denote the univariate and multivariate Gaussian distributions, respectively, where µ ∈ Rn and Σ ∈ Rn×n is a positive semi-definite matrix. We define a family of trimodal Gaussian mixture distributions as follows
Nmixσ,µ , 1 3 N (0, σ2) + 1 3 N (−µ, σ2) + 1 3 N (µ, σ2) .
For an illustration, please see Fig. 1.
Let χ2(k, λ) denote the noncentral chi-squared distribution with k degrees of freedom and the non-centrality parameter λ. For example, if Xi ∼ N (µi, 1) (for i = 1, 2, . . . , k) are independent Gaussian random variables, we have ∑k i=1X 2 i ∼ χ2(k, λ), where λ = ∑k i=1 µ 2 i . We also denote by χ2(k) the (central) chi-squared distribution with k degrees and the F -distribution by F (d1, d2) where d1 and d2 are the degrees of freedom.
Problem Setup. Let x1, . . . , xn ∈ RD be column vectors that represent the training data of size n and let xtest ∈ RD be a column vector that represents the test data. We assume that they are all independently drawn from a distribution
x1, . . . , xn, xtest iid∼ D .
Let us consider a linear regression problem on the first d features, where d ≤ D for some arbitrary large D. Here, d can be viewed as the number of features revealed. Then the feature vectors are x̃1, . . . , x̃n, where x̃i = xi[1 : d] ∈ Rd denotes the first d entries of xi. The corresponding response variable yi satisfies
yi = x̃ > i β + εi , i = 1, . . . , n ,
where the noise εi ∼ N (0, η2). We use the same setup as in [34] (see Equations (1) and (2) in [34]). Moreover, in another closely related work [41], if the kernel is set to the linear kernel, it is equivalent to our setup.
Next, we introduce the estimate β̂ of β and its excess generalization loss. Let ε = (ε1, . . . , εn)> ∈ Rn denote the noise vector. The design matrix A equals [x̃1, . . . , x̃n]> ∈ Rn×d. Let x = xtest[1 : d] denote the first d features of the test data. For the underparametrized regime where d < n, the least square solution on the training data is A+(Aβ + ε). For the overparametrized regime where d > n, A+(Aβ + ε) is the minimum-norm solution. In both regimes we consider the solution
β̂ , A+(Aβ + ε). The excess generalization loss on the test data is then given by Ld , E [( y − x>β̂ )2 − ( y − x>β )2] = E [( x>(β̂ − β)
)2] = E [( x> ( (A+A− I)β +A+ε
))2] = E [ (x>(A+A− I)β)2 ] + E [ (x>A+ε)2
] = E [ (x>(A+A− I)β)2 ] + η2E
∥∥(A>)+x∥∥2 , (1) where y = x>β + εtest and εtest ∼ N (0, η2). We call the term E [ (x>(A+A− I)β)2 ] the bias and
call the term η2E ∥∥(A>)+x∥∥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. [34] and shown in their Proposition 2. Remark 1. In the underparametrized regime, if D 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 E [ (x>(A+A− I)β)2 ] vanishes irrespective of β. In other words, in the underparametrized regime, Ld equals η2E‖(A>)+x‖2.
According to Remark 1, we have Ld = η2E‖(A>)+x‖2 in the underparametrized regime. It also holds in the overparametrized regime when β = 0. Without loss of generality, we assume η = 1 in the underparametrized regime (for all β). In the overparametrized regime, we also assume η = 1 for the β = 0 case. In this case, we have
Ld = E‖(A>)+x‖2 . (2) We assume a general η (i.e., not necessarily being 1) in the overparametrized regime when β is non-zero.
We would like to study the change in the loss caused by the growth in the number of features revealed. Recall Ld = E‖(A>)+x‖2. Once we reveal a new feature, which adds a new row b> to A> and a
new component a1 to x, we have Ld+1 = E ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2 .
Local Maximum and Multiple Descent. Throughout the paper, we say that a local maximum occurs at a dimension d ≥ 1 if Ld−1 < Ld and Ld > Ld+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 L0 , −∞. 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.

heading: 4 Underparametrized Regime
text: 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 Ld is always non-decreasing as d grows. Moreover, it is possible to have an arbitrarily large ascent, i.e., Ld+1 − Ld > C for any C > 0. Theorem 1 (Proof in Section 4.1). If d < n, we have Ld+1 ≥ Ld irrespective of the data distribution. Moreover, for any C > 0, there exists a distribution D such that Ld+1 − Ld > C. Remark 2 (D can be a product distribution). 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 Ld+1 − Ld > C) to hold, one extremely simple and elegant choice of the distribution D is a product distribution D = D1 × · · · × DD such that xi,j
iid∼ Dj for all 1 ≤ i ≤ n, where Dj is a Gaussian mixture Nmixσj ,1 for some σj > 0. Since the second part of Theorem 1 is of independent interest, the result is summarized by Theorem 4.
Remark 3 (Kernel regression on Gaussian data). In light of Remark 2, D can be chosen to be a product distribution that consists Nmixσj . Note that one can simulate N mix σ,1 with N (0, 1) through the inverse transform sampling. To see this, let FN (0,1) and FNmixσ,1 be the cdf of N (0, 1) and Nmixσ,1 , respectively. If X ∼ N (0, 1), we have FN (0,1)(X) ∼ Unif((0, 1)) and therefore ϕσ(X) , F−1Nmixσ,1 (FN (0,1)(X)) ∼ Nmixσ,1 . In fact, we can use a multivariate Gaussian D′ = N (0, ID×D) and a sequence of non-linear kernels k[1:d](x, x′) , 〈φ[1:d](x), φ[1:d](x′)〉, where the feature map is φ[1:d](x) , [φ1(x1), φ2(x2), . . . , φd(xd)]> ∈ Rd. Here is a simple rule for defining φj : if Dj = Nmixσj , we set φj to ϕσj . Thus, the problem becomes a kernel regression problem on the standard Gaussian data.
The first part of Theorem 1, which says that Ld is increasing (or more precisely, non-decreasing), agrees with Figure 1 of [11] and Proposition 2 of [34]. In [34], they proved that the risk increases with γ = 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 D inductively to satisfy Ld+1 − Ld > C. We fix d. Again, denote the first d features of xtest by x , xtest[1 : d]. Let us add an additional component to the training data x1[1 : d], . . . , xn[1 : d] and test data x so that the dimension d is incremented by 1. Let bi ∈ R denote the additional component that we add to the vector xi (so that the new vector is given as [xi[1 : d]>, bi]>. Similarly, let a1 ∈ R denote the additional component that we add to the test vector x. We form the column vector b = [b1, . . . , bn]
> ∈ Rn that collects all additional components that we add to the training data. We consider the change in the generalization loss as follows
Ld+1 − Ld = E ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2 − ∥∥(A+)>x∥∥2  . (3) Note that the components b1, . . . , bn, a1 are i.i.d. The proof of Theorem 1 starts with Lemma 2 which
relates the pseudo-inverse of [A, b]> to that of A>. In this way, we can decompose ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2
into multiple terms for further careful analysis in the proofs hereinafter.
Lemma 2 (Proof in Appendix B.1). Let A ∈ Rn×d and 0 6= b ∈ Rn×1, where n ≥ d + 1. Additionally, let P = AA+ and Q = bb+ = bb >
‖b‖2 , and define z , b>(I−P )b ‖b‖2 . If z 6= 0 and the
columnwise partitioned matrix [A, b] has linearly independent columns, we have[ A>
b>
]+ = [( I − bb >
‖b‖2
)( I + AA +bb>
‖b‖2−b>AA+b
) (A+)>, (I−AA
+)b ‖b‖2−b>AA+b ] = [ (I −Q)(I + PQ1−tr(PQ) )(A +)>, (I−P )b b>(I−P )b
] = [ (I −Q)(I + PQz )(A +)>, (I−P )b b>(I−P )b ] .
In our construction of D, the components Dj are all continuous distributions. The matrix I − P is an orthogonal projection matrix and therefore rank(I − P ) = n− d. As a result, it holds almost surely that b 6= 0, z 6= 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 Ld = E ∥∥(A+)>x∥∥2 is finite and the (d+ 1)-th features b1, . . . , bn, a1
are i.i.d. sampled from N (0, 1) or Nmixσ,1 , Ld+1 = E ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2 is also finite.
Theorem 3 (Proof in Appendix B.2). Let z be as defined in Lemma 2. If b1, . . . , bn, a1 are i.i.d. and follow a distribution with mean zero, conditioned on A and x, we have
Eb,a1 ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2  ≤ Eb,a1 [1z ∥∥(A+)>x∥∥2 + a21b>(I − P )b ] .
In particular, if d+ 2 < n and b1, . . . , bn, a1 iid∼ N (0, 1), conditioned on A and x, we have
Eb,a1 ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2  ≤ (n− 2)∥∥(A+)>x∥∥2 + 1 n− d− 2 .
If d+ 2 < n and b1, . . . , bn, a1 iid∼ Nmixσ,1 , conditioned on A and x, we have
Eb,a1 ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2 ≤ (n− 2 + √ d) ∥∥(A+)>x∥∥2 + 2/(3σ2) + 1 n− d− 2 .
Using Theorem 3, we can show inductively (on d) that Ld is finite for every d. Provided that we are able to guarantee finite L1, Theorem 3 implies that Ld is finite for every d if the components are always sampled from N (0, 1) or Nmixσ,1 .
Making a large Ld can be achieved by adding an entry sampled from Nmixσ,1 when the data dimension increases from d − 1 to d in the previous step. Theorem 4 shows that adding a Nmixσ,1 feature can increase the loss by arbitrary amount, which in turn implies the second part of Theorem 1.
Theorem 4 (Proof in Appendix B.4). For any C > 0 and E ∥∥(A+)>x∥∥2 < +∞, there exists a
σ > 0 such that if b1, . . . , bn, a1 iid∼ Nmixσ,1 , we have
E ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2 − ∥∥(A+)>x∥∥2  > C . We are now ready to prove Theorem 1.

heading: 4.1 Proof of Theorem 1
text: Proof. We follow the notation convention in (3):
Ld+1 − Ld = E ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2 − ∥∥(A>)+x∥∥2  . Recall d < n and the matrix B′ , [ A>
b>
] is of size (d+ 1)× n. Both matrices B′ and B , A> are
fat matrices. As a result, if x′ , [ x a1 ] , we have
‖B′+x′‖2 = min z:B′z=x′ ‖z‖2 , ‖B+x‖2 = min z:Bz=x ‖z‖2 .
Since {z | B′z = x′} ⊆ {z | Bz = x}, we get ‖B′+x′‖2 ≥ ‖B+x‖2. Therefore, we obtain Ld+1 ≥ Ld. The second part follows from Theorem 4.
Remark 4. Remark 2 and the proof of Theorem 4 indicate that D = D1 × · · · × DD is a product distribution. The construction in the proof also shows that the generalization curve is determined by the specific choice of the Di’s. Note that permuting the order of Di’s is equivalent to changing the order by which the features are being revealed (i.e., permuting the entries of the data xi’s). Therefore, given the same data points x1, · · · , xn ∈ RD, one can create different generalization curves simply by changing the order of the feature-revealing process.

heading: 5 Overparametrized Regime
text: 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 β = 0 and Ld is as defined in (2). Then we discuss the setting β 6= 0. As stated in the following theorem, we require d ≥ 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.
Theorem 5 (Overparametrized regime, β = 0). Let n < D − 9. Given any sequence ∆n+8,∆n+9, . . . , ∆D−1 where ∆d ∈ {↑, ↓}, there exists a distribution D such that for every n+ 8 ≤ d ≤ D − 1, we have
Ld+1 { > Ld, if ∆d = ↑ < Ld, if ∆d = ↓ .
In Theorem 5, the sequence ∆n+8, ∆n+9, · · · , ∆D−1 is just used to specify the increasing/decreasing behavior of the Ld sequence for d > n + 8. Compared to Theorem 1 for the underparametrized regime, where Ld 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.
Lemma 6 (Proof in Appendix C.1). Let A ∈ Rn×d and b ∈ Rn×1, where n ≤ d. Assume that matrix A and the columnwise partitioned matrix B , [A, b] have linearly independent rows. Let G , (AA>)−1 ∈ Rn×n and u , b
>G 1+b>Gb ∈ R1×n. We have[ A>
b>
]+ = [ (I − bu)>(A+)>, u> ] .
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 Dn+8, . . . ,DD 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.
Lemma 7 (Proof in Appendix C.2). Let D = D1 × · · · × DD be a product distribution where
(a) Dd = N (0, 1) if d = 1, . . . , n+ 8; and
(b) Dd is either N (0, σ2d) or Nmixσd,µd for d > n+ 8.
Let D[1:d] denote D1 × · · · × Dd. Assume that every row of A ∈ Rn×d and x ∈ Rd×1 are i.i.d. and follow D[1:d]. For any d such that n+ 8 ≤ d ≤ D, all of the followings hold:
E[‖(A+)>x‖2] < +∞ , E[λ2max((AA>)−1)] < +∞ , E[λmax((AA>)−1)‖(A+)>x‖2] < +∞ , E[λ2max((AA>)−1)‖(A+)>x‖2] < +∞ . (4)
Theorems 8 and 9 are the key technical results for constructing multiple descent in the overparametrized regime. One can create a descent (Ld+1 < Ld) by adding a Gaussian feature (Theorem 8) and create an ascent (Ld+1 > Ld) by adding a Gaussian mixture feature (Theorem 9). Theorem 8 (Proof in Appendix C.3). If E[‖(A>A)+x‖2] > 0 and all equations in (4) hold, there exists σ > 0 such that if a1, b1, . . . , bn iid∼ N (0, σ2), we have
Ld+1 − Ld = E ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2 − E ∥∥(A+)>x∥∥2 < 0 .
Theorem 9 shows that adding a Gaussian mixture feature can make Ld+1 > Ld. Theorem 9 (Proof in Appendix C.4). Assume E‖(A+)>x‖2 < +∞. For any C > 0, there exist µ, σ > 0 such that if a1, b1, . . . , bn iid∼ Nmixσ,µ , we have
Ld+1 − Ld = E ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2 − E ∥∥(A+)>x∥∥2 > C .
The proof of Theorem 5 immediately follows from Theorem 8 and Theorem 9. Proof of Theorem 5. We construct the product distribution D = ∏D d=1Dd. We set Dd = N (0, 1) for d = 1, . . . , n + 8. For n + 8 < d ≤ D, Dd is either N (0, σ2d) or Nmixσd,µd depending on ∆d being either ↓ or ↑. First we show that for each step d, the assumption E[‖(A>A)+x‖2] > 0 of Theorem 8 is satisfied. If E[‖(A>A)+x‖2] = 0, we know that (A>A)+x = 0 almost surely. Since D is a continuous distribution, the matrix A has full row rank almost surely. Therefore, rank((A>A)+) = rank(A>A) = n almost surely. Thus dim ker(A>A)+ = d−n ≤ d−1 almost surely, which implies x /∈ ker(A>A)+. In other words, (A>A)+x 6= 0 almost surely. We reach a contradiction. Moreover, by Lemma 7, the assumption E‖(A+)>x‖2 < +∞ of Theorem 9 is also satisfied. If ∆d−1 = ↓, by Theorem 8, there exists σd > 0 such that if Dd = N (0, σ2d), then Ld < Ld−1. Similarly if ∆d−1 = ↑, by Theorem 9, there exists σd and µd such that Dd = Nmixσd,µd guarantees Ld > Ld−1.
Gaussian β setting. In what follows, we study the case where the model β is non-zero. In particular, we consider a setting where each entry of β is i.i.d. N (0, ρ2). Recalling (1), define the biases
Ed , (x>(A+A− I)β)2, Ed+1 , ( [x>, a1]([A, b] +[A, b]− I) [ β β1 ])2 ,
and the expected risks
Lexpd , E[Ed] + η 2E ∥∥(A>)+x∥∥2 , Lexpd+1 , E[Ed+1] + η2E ∥∥∥∥∥ [ A> b> ]+ [ x a1 ]∥∥∥∥∥ 2 , (5)
where β ∼ N (0, ρ2Id) and β1 ∼ N (0, ρ2). The second term in Lexpd and L exp d+1 is the variance term. Note that Lexpd is the expected value of Ld in (1) and averages over β. Theorem 10 shows that one can add a Gaussian mixture feature in order to make Lexpd+1 > L exp d , and add a Gaussian feature in order to make Lexpd+1 < L exp d .
Theorem 10 (Proof in Appendix C.5). Let a1, β1 ∈ R, x ∈ Rd×1, β ∈ Rd×1, A ∈ Rn×d and b ∈ Rn×1, where n ≤ d. Assume that x, a1, β1, β, A, b are jointly independent, [β>, β1]> ∼ N (0, ρ2Id+1). Moreover, assume that the matrix [A, b] has linearly independent rows almost surely. The following statements hold:
(a) If a1, b1, . . . , bn iid∼ Nmixσ,µ , for any C > 0, there exist µ, σ such that L exp d+1 − L exp d > C.
(b) If a1, b1, . . . , bn iid∼ N (0, σ2), there exists σ > 0 such that for all
ρ ≤ η √ E[‖(A>A)+x‖2] E‖A+>x‖2 + 1 ,
we have Lexpd+1 < L exp d .
Theorem 10 indicates that for β 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 β) and overparametrization regime (for β = 0), the distribution D can be made a product distribution. Theorem 11 (Overparametrized regime, β being Gaussian). Let n < D − 9. Given any sequence ∆n+8,∆n+9, . . . , ∆D−1 where ∆d ∈ {↑, ↓}, there exists ρ > 0 and a distribution D such that for β ∼ N (0, ρ2) and every n+ 8 ≤ d ≤ D − 1, we have
Lexpd+1 { > Lexpd , if ∆d = ↑ < Lexpd , if ∆d = ↓ .
Proof of Theorem 11. Define the design matrixAd , [x1[1 : d], . . . , xn[1 : d]]> ∈ Rn×d. Similar to the proof of Theorem 5, we construct the product distribution D = ∏D d=1Dd. We set Dd = N (0, 1) for d = 1, . . . , n+ 8. For n+ 8 < d ≤ D, Dd is either N (0, σ2d) or Nmixσd,µd depending on ∆d being either ↓ or ↑. If ∆d−1 = ↑, by Theorem 10, there exists σd and µd such thatDd = Nmixσd,µd guarantees L exp d > L exp d−1. If ∆d−1 = ↓, define
ρd , η
√ E[‖(A>d−1Ad−1)+xtest[1 : d− 1]‖2]
E‖A+>d−1xtest[1 : d− 1]‖2 + 1 .
By Theorem 10, there exists σd > 0 such that if ρ ≤ ρd and Dd = N (0, σ2d), then L exp d < L exp d−1. We take ρ = min
d:∆d−1=↓ ρd .

heading: 6 Conclusion
text: 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).

Review summaries:
================================= FINAL UPDATE AFTER DISCUSSION PHASE  As the authors addressed the points that came up during the discussion phase (which made me initially lower the score), I revert back to my initial scoring of accept.   =================================== UPDATE AFTER INITIAL DISCUSSION PHASE  The authors addressed my own points, but following the discussion with the other reviewers I decided to lower my score. It was pointed out that the data generation process depends on the model-dimension, and thus the resulting investigation is not about 'generalization curves' as claimed by the authors, but rather about the generalization behavior of different interpolation regimes. While I believe the contribution is still valuable, some explanations and need be changed (which may be possible in a camera ready version).  ( For example I don't believe anymore that the paper indeed answers the question: 'Can the existence of a multiple descent generalization curve be rigorously proven?'.)  ==================================  The paper analyses the behavior of generalization curves of linear regression wrt an increasing feature dimensionality. The paper constructs specific examples that can generate any ascending/ descending behavior of the curve, for both the under and overparametrized case. Their construction highlights that ascending/ descending behavior of the generalization curve is not inherent to the model family, but rather a result of the data distribution in combination with the model.

I have read other reviewers' comments and the author response. I lower my score due to the following reasons: 1) changing the feature dimension $d$ doesn't change the groundtruth function, which makes the problem setting questionable. 2) the response to the issue "\sigma should be taken far smaller than O(1/5^n)" is unsatisfactory, which makes the result unattainable.  ============================================= This paper studies the multiple descent of the excess risk for linear regression in under/over-parameterized regimes. In the under-parameterized regime, Theorem 3.1 includes two parts: one is showing that the loss $L_d$ non-decreases with $d$ independent of the data distribution; the other is $L_{d+1} - L_d > C$. In the over-parameterized regime, the authors provide the multiple descents in Theorem 5 under $\beta = 0$ and Theorem 11 with $\beta$ being Gaussian.

This paper studies the generalization error of the min-l2 norm solutions of linear regression. Authors claim that their analytical results can explain the multiple peaks in the generalization curve with respect to the number of features. However, the system model and some fundamental steps of the analysis are problematic, which makes me doubt the correctness of the paper.