Year: 2020

Abstract:
This paper investigates the finite-sample prediction risk of the high-dimensional least squares estimator. We derive the central limit theorem for the prediction risk when both the sample size and the number of features tend to infinity. Furthermore, the finite-sample distribution and the confidence interval of the prediction risk are provided. Our theoretical results demonstrate the sample-wise nonmonotonicity of the prediction risk and confirm “more data hurt” phenomenon.

heading: 1 INTRODUCTION
text: More data hurt refers to the phenomenon that training on more data can hurt the prediction performance of the learned model, especially for some deep learning tasks. Loog et al. (2019) shows that various standard learners can lead to sample-wise non-monotonicity. Nakkiran et al. (2019) experimentally confirms the sample-wise non-monotonicity of the test accuracy on deep neural networks. This challenges the conventional understanding in large sample properties: if an estimator is consistent, more data makes the estimator more stable and improves its finite-sample performance. Nakkiran (2019) considers adding one single data point to a linear regression task and analyzes its marginal effect to the test risk. Dereziński et al. (2019) gives an exact non-asymptotic risk of the high-dimensional least squares estimator and observes the sample-wise non-monotonicity on mean square error. For adversarially robust models, Min et al. (2020) proves that more data may increase the gap between the generalization error of adversarially-trained models and standard models. Chen et al. (2020) shows that more training data causes the generalization error to increase in the strong adversary regime. In this work, we derive the finite-sample distribution of the prediction risk under linear models and prove the “more data hurt” phenomenon from an asymptotic point of view.
Intuitively, the “more data hurt” stems from the “double descent” risk curve: as the model complexity increases, the prediction risk of the learned model first decreases and then increases, and then decreases again. The double descent phenomenon can be precisely quantified for certain simple models (Hastie et al. (2019); Mei & Montanari (2019); Ba et al. (2019); Belkin et al. (2019); Bartlett et al. (2020); Xing et al. (2019)). Among these works, Hastie et al. (2019) and Mei & Montanari (2019) use the tools from random matrix theory and explicitly prove the double descent curve of the asymptotic risk of linear regression and random features regression in high dimensional setup. Ba et al. (2019) gives the asymptotic risk of two-layer neural networks when either the first or the second layer is trained using a gradient flow.
The second decline of the prediction risk in the double descent curve is highly related to the more data hurt phenomenon. In the over-parameterized regime when the model complexity is fixed while the sample size increases, the degree of over-parameterization decreases and becomes close to the interpolation boundary (for example p/n = 1 in Hastie et al. (2019)), in which a high prediction risk is achieved. However, the existing asymptotic results, which focus on the first-order limit of the prediction risk, cannot fully describe the more data hurt phenomenon. In fact, the “double descent” curve is a function of the limiting ratio lim p/n, which may not be able to characterize the empirical prediction risk in finite sample situations. There will be a non-negligible discrepancy between the empirical prediction risk and its limit, especially when the sample size or dimension is small. Finegrained second-order results are thus needed to fully characterize such discrepancy and further, a confidence band for the prediction risk can be constructed to evaluate its finite sample performance. We take Figure 1 as an example to illustrate this. According to the first-order limit, given a fixed dimension p = 100, the prediction risks at sample size n = 90 and n = 98 are about 10.20 and
49.02. More data hurt seems true. However, the 95% confidence interval of the prediction risks with sample size 98 is [4.91, 142.12], which contains the risk for n = 90. Then more data hurt is not statistically significant. Hence, in this work, we characterize the second-order fluctuations of the prediction risk and make attempts to fill this gap. We employ the linear regression task in Hastie et al. (2019) and Nakkiran (2019), and introduce new tools from the random matrix theory, e.g. the central limit theorems for linear spectral statistics in Bai & Silverstein (2004); Bai et al. (2007), to derive the central limit theorem of the prediction risk.
Consider a linear regression task with n data points and p features, the setup of the more data hurt is similar to that in the classical asymptotic analysis in Van der Vaart (2000). According to the classical asymptotic analysis with p fixed and n → ∞, the least square estimator is unbiased and√ n-consistent to the ground truth. This implies that the more data will not hurt and even improve the prediction performance. However, the story is very different in the over-parameterized regime. The prediction risk doesn’t decrease monotonously with n when p > n. More data does hurt in the over-parameterized case. In the following, we will justify this phenomenon by developing the second-order asymptotic results as both n and p tend to infinity. We assume p/n → c, and denote 0 < n1 < n2 < +∞, c1 = p/n1 and c2 = p/n2. Then the direct comparison of the prediction risk between sample sizes n1 and n2 can be decomposed into three parts: (i) the gap between the finitesample risk under n = n1 and the asymptotic risk with c = c1; (ii) the gap between the finite-sample risk under n = n2 and the asymptotic risk with c = c2; (iii) the comparison between two asymptotic risks under c = c1 and c = c2. Theorem 1 and 2 of Hastie et al. (2019) give answers to the task (iii). For (i) and (ii), we develop in this paper the convergence rate and the limiting distribution of the prediction risk as n, p → +∞, p/n → c. Furthermore, the confidence interval of the finite-sample risk can be obtained as well.
The main goal of this paper is to study the second order asymptotic behavior of two different types of conditional prediction risk in the linear regression model. One is RX,β(β̂,β) given both the training data and regression coefficient while the other is RX(β̂,β) given the training data only. We summarize our main results as follows: (1) The regression coefficient is set to be either random or nonrandom to cover more cases. Different convergence rates and limiting distributions of both prediction risk are derived under various scenarios. (2) In particular, the finite-sample distribution of the conditional prediction risk given both the training data and regression coefficient is derived and the sample-wise double descent is characterized in Theorem 4.2 and Theorem 4.5 (see Figure 1). Under certain assumptions, the more data hurt phenomenon can be confirmed by comparing the confidence intervals built via the central limit theorems. (3) Our results incorporate non-Gaussian observations. For Gaussian data, the limiting mean and variance in the central limit theorems have simpler forms, see Section 4.2 and 4.3 for more details.
The rest of this paper is organized as follows. Section 3 introduces the model settings and two different prediction risk. Section 4 presents the main results on CLTs for the two types of risk with
discussion. Section 5 conducts simulation experiments to verify the main results. All the technical proofs and lemmas are relegated to the appendix in the supplementary file.

heading: 2 RELATED WORK
text: Double Descent The double descent curve describes how generalization ability changes as model capacity increases. It subsumes the classical bias-variance trade-off, a U-shape curve, and further show that the test error exhibits a second drop when the model capacity exceeds the interpolation threshold (Belkin et al. (2018); Geiger et al. (2019); Spigler et al. (2019); Advani & Saxe (2017)). The double descent phenomenon has been quantified for certain models, including two-layer neural networks via non-asymptotic bounds or asymptotic risk (Belkin et al. (2019); Muthukumar et al. (2020); Hastie et al. (2019); Mei & Montanari (2019); Ba et al. (2019)). As our results are based on the linear regression, we mainly focus on the literature of linear models. Muthukumar et al. (2020) and Bartlett et al. (2020) derive the generalization bounds for over-parametrized linear models and show the benefits of the interpolation. Hastie et al. (2019) gives the first-order limit of the generalization error for linear regressions as n, p → +∞. Dereziński et al. (2019) provides an exact non-asymptotic expression for double descent of the high-dimensional least square estimator. Wu & Xu (2020) extends the first-order limit of the prediction error of the generalized weighted ridge estimator to more general case with anisotropic features and signals. Montanari et al. (2019), Deng et al. (2019) and Kini & Thrampoulidis (2020) investigate the sharp asymptotic of binary classification tasks with the max-margin solution and the maximum likelihood solution. Emami et al. (2020) and Gerbelot et al. (2020a) consider the double descent in generalized linear models. Furthermore, the double descent phenomenon is also observed on linear tasks with various problems and assumptions, e.g. LeJeune et al. (2020); Gerbelot et al. (2020b); Javanmard et al. (2020); Dar & Baraniuk (2020); Xu & Hsu (2019); Dar et al. (2020). Xing et al. (2019) sharply quantifies the benefit of interpolation in the nearest neighbors algorithm. Mei & Montanari (2019) derives the limiting risk on the random features model and shows that minimum generalization error is achieved by highly overparametrized interpolators. Ba et al. (2019) gives the limiting risk of the regression problem under two-layer neural networks. However, the existing asymptotic results focus on the first-order limit of the prediction risk and do not indicate the convergence rate. There are very few second-order results in the literature, Shen & Bellec (2020) establishes the asymptotic normality for the derivatives of 2-layers neural network, but not the exact limiting distribution of the risk. In this work, we are the first to develop results on second-order fluctuations of the prediction risk in linear regressions and provide its corresponding confidence intervals. The more data hurt phenomenon is further justified from the asymptotic point of view.
Random Matrix Theory The primary tool for analyzing the second-order fluctuations of prediction risk comes from random matrix theory. In particular, Bai & Silverstein (2004) refines the central limit theorem for linear spectral statistics of large dimensional sample covariance matrix with general population and the population is not necessary to be Gaussian. Similar central limit theorems are also developed for other random matrix ensembles, see Sinai & Soshnikov (1998); Bai & Yao (2005); Zheng (2012). Other than the central limit theorem for linear spectral statistics, Bai et al. (2007) and Pan & Zhou (2008) study the asymptotic fluctuation of eigenvectors of sample covariance matrices. Bai & Yao (2008) considers the fluctuation of quadratic forms. All these technical tools and results are adopted and fully utilized in this paper, especially those related to Stieltjes transform, which are closely connected to the prediction risk studied in this paper.

heading: 3 PRELIMINARIES

heading: 3.1 PROBLEM, DATA AND ESTIMATOR
text: Suppose that the training data {(xi, yi) ∈ Rp×R, i = 1, 2, . . . , n} is generated independently from the model (ground truth or teacher model):
yi = β Txi + i, and (xi, i) ∼ (Px, P ), i = 1, 2, . . . , n. (1)
Here, Px is a distribution on Rp such that E(xi) = 0, Cov(xi) = Σ, and P is a distribution on R such that E( i) = 0, Var( i) = σ2. In particular, the coordinates of xi are not necessarily
independent, that is, Σ is not restricted to be diagonal. To proceed further, we denote
Xn×p = (x1,x2, . . . ,xn) T, y = (y1, y2, . . . , yn) T.
The minimum `2 norm (min-norm) least squares estimator, of y on X, is defined by
β̂ = arg min β
‖y −Xβ‖2 = (XTX)+XTy, (2)
where (XTX)+ denotes the Moore-Penrose pseudoinverse of XTX.

heading: 3.2 BIAS, VARIANCE AND RISK
text: Similar to Hastie et al. (2019), we define two different types of out-of-sample prediction risk. The first one is given by
RX(β̂,β) = E [ (xT0 β̂ − xT0β)2 ∣∣X] = E[‖β̂ − β‖2Σ∣∣X], (3) where x0 ∼ Px is a test point and is independent of the training data, and the notation ‖β‖2Σ stands for βTΣβ. Here β is assumed to be a random vector independent of x0. In this definition, the expectation E stands for the conditional expectation for x0, β̂ and β when X is given. According to the bias-variance decomposition, we have RX(β̂,β) := BX(β̂,β) + VX(β̂,β), where
BX(β̂,β) = E { ‖E(β̂|X)− β‖2Σ ∣∣X} and VX(β̂,β) = Tr{Cov(β̂|X)Σ}. (4) Plugging the model (1) into the min-norm estimator (2), the bias and variance terms can be rewritten as
BX(β̂,β) = E { βTΠΣΠβ ∣∣X} and VX(β̂,β) = σ2 n Tr(Σ̂+Σ),
where Σ̂ = XTX/n is the (uncentered) sample covariance matrix of X, and Π = Ip − Σ̂+Σ̂ is the projection onto the null space of X.
The second type of out-of-sample prediction risk is defined as RX,β(β̂,β) = E [ (xT0 β̂ − xT0β)2 ∣∣X,β] = E[‖β̂ − β‖2Σ∣∣X,β], (5) where
BX,β(β̂,β) = β TΠΣΠβ and VX,β(β̂,β) = VX(β̂,β) =
σ2
n Tr(Σ̂+Σ).
In this definition, the parameter β is assumed to be given. The expectation E is the conditional expectation for x0 and β̂ when X and β are given. This is consistent with the commonly-used testing procedure, in which a trained model is evaluated by the average loss on unseen testing data.

heading: 4 MAIN RESULTS
text: Before stating our main results, we briefly highlight the challenges we faced in proving the more data hurt phenomenon. First, the finite-sample behaviors of the prediction risk is required. Hastie et al. (2019) gives the first-order limits of both RX,β(β̂,β) and RX(β̂,β) as n, p → +∞ and p/n → c ∈ (0,+∞). However, to prove the more data hurt phenomenon, we should fix p and investigate the finite-sample risk with sample size n varies. This implies that only knowing the firstorder limit is not enough, the convergence rate is also needed. To solve this problem, we have derived the central limit theorems for both RX,β(β̂,β) and RX(β̂,β), respectively, which characterize the second-order fluctuations of the risk. Then we can figure out the finite-sample behavior of the risk by computing the gap between the risk and its limit. The confidence intervals of the risk can be further obtained. Second, the parameter β also contributes randomness to the finite-sample risk, which further influences the convergence rate. To analyze the contribution of β, we need to make use of the technical tools and asymptotic results for eigenvectors and quadratic forms developed in Bai et al. (2007) and Bai & Yao (2008). Another interesting finding is that, in the over-parameterized regime such that p > n, the two types of out-of-sample prediction risk RX,β(β̂,β) and RX(β̂,β) enjoy different convergence rates.

heading: 4.1 ASSUMPTIONS AND MORE NOTATIONS
text: As follows are some notations used in this paper. The p× p identity matrix is denoted by Ip. For a symmetric matrixA ∈ Rp×p, we define its empirical spectral distribution as
FA(x) = 1
p p∑ i=1 1{λi(A) ≤ x}
where 1{·} is the indicator function and λi(A), i = 1, 2, . . . p are the eigenvalues ofA. The notation d−→ stands for the convergence in distribution. Zα/2 is the α/2 upper quantile of the standard normal distribution, λmax(A) and λmin(A) denote the largest and smallest eigenvalues ofA, respectively.
Here we list all the assumptions for X and β needed under different scenarios:
(A) xj ∼ Px is of the form xj = Σ1/2zj , where zj is a p-length random vector with i.i.d. entries that have zero mean, unit variance, and a finite 4-th order moment E(z4ij) = ν4, i = 1, · · · , p, j = 1, · · · , n.
(B1) Σ is a deterministic positive definite matrix, such that 0 < c0 ≤ λmin(Σ) ≤ λmax(Σ) ≤ c1, for all n, p and some constants c0, c1. As p→∞, we assume that the empirical spectral distribution FΣ converges weakly to a probability measure H.
(B2) Σ is an identity matrix, Σ = Ip.
(C1) β is a nonrandom constant vector, and ‖β‖22 = βTβ = r2. (C2) β ∼ Pβ is independent of X and follows multivariate Gaussian distribution Np(0, r 2
p Ip).
Throughout this paper, we consider the limiting distributions and the convergence rates of the outof- sample prediction risk when n, p→∞ such that p/n = cn → c ∈ (0,∞). If c > 1, the sample size n is smaller than the number of parameters p, we call this case “over-parametrized”. Otherwise when c < 1, we call it “under-parameterized”.

heading: 4.2 UNDER-PARAMETRIZED ASYMPTOTICS
text: In this section, we focus on the risk of the min-norm estimator (2) in the under-parametrized regime. According to Theorem 1 of Hastie et al. (2019), both BX,β(β̂,β) and BX(β̂,β) converge to σ2c/(1 − c) almost surely. The following Theorem 4.1 and 4.2 show that both BX(β̂,β) and BX,β(β̂,β) converge to σ2c/(1− c) at the rate 1/p. Furthermore, the limiting distributions are derived by making use of the CLT for linear spectral statistics of large-dimensional sample covariance matrices. Theorem 4.1. Suppose that the training data is generated from the model (1), and the assumptions (A) and (B1) hold. Then the first type of out-of-sample prediction risk RX(β̂,β) of the min-norm estimator (2) satisfies that, as n, p→∞ such that p/n = cn → c < 1,
p ( RX(β̂,β)− cnσ 2
1− cn
) d−→ N(µc, σ2c ), (6)
where
µc = c2σ2 (c− 1)2 + σ2c2(ν4 − 3) 1− c and σ2c = 2c3σ4 (c− 1)4 + c3σ4(ν4 − 3) (1− c)2 .
Conclusively, P (Lα,c ≤ RX(β̂,β) ≤ Uα,c)→ 1− α, (7) where 1− α is the confidence level and
Lα,c = cnσ
2
1− cn +
1 p (µc − Zα/2σc), Uα,c =
cnσ 2
1− cn +
1 p (µc + Zα/2σc).
Under the assumptions of Theorem 4.1, we know that Π = Ip − Σ̂+Σ̂ = 0 and
BX(β̂,β) = BX,β(β̂,β) = 0, VX(β̂,β) = VX,β(β̂,β) = σ2
n Tr(Σ̂+Σ).
Thus RX(β̂,β) equals to RX,β(β̂,β) and the two risk share the same asymptotic limit.
Theorem 4.2. Under the assumptions of Theorem 4.1, the second type of out-of-sample prediction risk RX,β(β̂,β) of the min-norm estimator (2) satisfies that, as n, p → ∞ such that p/n = cn → c < 1,
p ( RX,β(β̂,β)− cnσ 2 1− cn ) d−→ N(µc, σ2c ),
and P (Lα,c ≤ RX,β(β̂,β) ≤ Uα,c)→ 1− α,
where µc, σ2c , Lα,c and Uα,c are the same as those in Theorem 4.1.

heading: 4.3 OVER-PARAMETRIZED ASYMPTOTICS
text: In this section, we consider the min-norm estimator (2) in the over-parametrized case c > 1. The bias term, either BX(β̂,β) or BX,β(β̂,β), is generally nonzero. According to Lemma 2 in Hastie et al. (2019), both BX(β̂,β) and BX,β(β̂,β) converge to r2(1− 1/c) as n, p→ +∞ and p/n→ c > 1. This implies that the bias term can influence the asymptotic behavior of the prediction risk, including the convergence rate. Hence to derive the CLT of the out-of-sample prediction risk, we need to consider both the bias and variance terms in (4).
In the following, we investigate the asymptotic properties of the two prediction risk RX(β̂,β) and RX,β(β̂,β) under various combinations of the assumptions (A1), (B2) for X and scenarios (C1), (C2) for both random and nonrandom β. We start with the case when β is a constant vector.
Theorem 4.3. Suppose that the training data is generated from the model (1), and the assumptions (A), (B2) and (C1) hold. Then the first type of out-of-sample prediction risk RX(β̂,β) of the minnorm estimator (2) satisfies that, as n, p→∞ such that p/n = cn → c > 1,
√ p { RX(β̂,β)− (1− 1
cn )r2 − σ
2
cn − 1
} d−→ N(µc,1, σ2c,1), (8)
where µc,1 = 0 and σ2c,1 = 2(c−1) c2 r 4. A more practical version is to replace µc,1 and σ2c,1 with
µ̃c,1 = 1 √ p { cσ2 (1− c)2 + σ2(ν4 − 3) c− 1 } and σ̃2c,1 = 2(c− 1) c2 r4 + 1 p { 2c3σ4 (1− c)4 + cσ4(ν4 − 3) (c− 1)2 } .
Conclusively, P (Lα,c ≤ RX(β̂,β) ≤ Uα,c)→ 1− α, (9)
where 1− α is the confidence level and
Lα,c = (1− 1
cn )r2 +
σ2
cn − 1 + 1 √ p (µ̃c,1 − Zα/2σ̃c,1),
Uα,c = (1− 1
cn )r2 +
σ2
cn − 1 + 1 √ p (µ̃c,1 + Zα/2σ̃c,1).
Remark 4.1. Under assumption (C1), BX(β̂,β) = BX,β(β̂,β) and RX(β̂,β) = RX,β(β̂,β). Thus Theorem 4.3 still holds if we replace RX(β̂,β) with RX,β(β̂,β).
Remark 4.2. Under Assumption (B2), the eigenvector of Σ̂ is asymptotically Haar distributed. Therefore, the bias term BX(β̂,β) is only related to the length of β. However, in the anisotropic settings with general Σ, the eigenvector of the Σ̂ is no longer asymptotically Haar distributed. The limiting behavior of BX(β̂,β) heavily relies on the interaction between β and the eigenvectors of Σ̂. Therefore, we conjecture that there is no universal convergence rate for the bias term BX(β̂,β) that can cover arbitrary non-random β and anisotropic Σ in the over-parametrized case, not to mention the prediction risk RX(β̂,β). A small simulation experiment is conducted in Appendix G to confirm our conjecture on this point.
Next we consider the case when β is a random vector that follows assumption (C2).
Theorem 4.4. Suppose that the training data is generated from the model (1), and the assumptions (A), (B2) and (C2) hold. Then as n, p → ∞ such that p/n = cn → c > 1, the first type of out-of-sample prediction risk RX(β̂,β) of the min-norm estimator (2) satisfies,
p { RX(β̂,β)− (1− 1
cn )r2 − σ
2
cn − 1
} d−→ N(µc,2, σ2c,2),
where
µc,2 = cσ2 (1− c)2 + σ2(ν4 − 3) c− 1 and σ2c,2 = 2c3σ4 (1− c)4 + cσ4(ν4 − 3) (c− 1)2 .
Hence we have P (Lα,c ≤ RX(β̂,β) ≤ Uα,c)→ 1− α, where
Lα,c = σ2 cn − 1 + (1− 1 cn )r2 + 1 p (µc,2 − Zα/2σc,2), Uα,c = σ2
cn − 1 + (1− 1 cn )r2 + 1 p (µc,2 + Zα/2σc,2).
As for RX,β(β̂,β), we have the following theorem. Theorem 4.5. Suppose that the training data is generated from the model (1), and the assumptions (A), (B2) and (C2) hold. Then, as n, p → ∞ such that p/n = cn → c > 1, the second type of out-of-sample prediction risk RX,β(β̂,β) of the min-norm estimator (2) satisfies,
√ p { RX,β(β̂,β)− (1− 1
cn )r2 − σ
2
cn − 1
} d−→ N(µc,3, σ2c,3), (10)
where µc,3 = 0 and σ2c,3 = 2(1− 1c )r 4. A more practical version is to replace µc,3 and σ2c,3 with
µ̃c,3 = 1 √ p
{ cσ2
(1− c)2 + σ2(ν4 − 3) c− 1
} and σ̃2c,3 = 2(1− 1
c )r4 +
1
p
{ 2c3σ4
(1− c)4 + cσ4(ν4 − 3) (c− 1)2
} .
The corresponding (1− α)-confidence interval is given by P (Lα,c ≤ RX,β(β̂,β) ≤ Uα,c)→ 1− α, (11)
with
Lα,c = σ2 cn − 1 + (1− 1 cn )r2 + 1 √ p (µ̃c,3 − Zα/2σ̃c,3),
Uα,c = σ2 cn − 1 + (1− 1 cn )r2 + 1 √ p (µ̃c,3 + Zα/2σ̃c,3).
Remark 4.3. Note that besides the leading constants in (µc,3, σc,3), the version (µ̃c,3, σ̃c,3) also contains smaller order terms, including terms of order O(1/ √ p) in µ̃c,3 and terms of order O(1/p) in σ̃c,3. These smaller order terms will vanish when p and n grow very large, but for finite sample situations, these smaller order terms will provide a finer approximation for the finite sample distribution of RX,β(β̂,β). As shown in the following experiments, these terms have indeed made non-negligible contributions to fitting the empirical distribution of RX,β(β̂,β), which sheds new lights for practitioners.
Remark 4.4. If we compare the results in Theorem 4.3 and 4.5, we will find out that RX(β̂,β) with constant β and RX,β(β̂,β) with random β share the same first-order limit and second-order error rate O(p−1/2). This is quite intuitive because both risk treat β as a constant. Their differences are reflected in their limiting variances. Nevertheless, it’s very interesting to observe from Theorem 4.4 that, RX(β̂,β) with random β under the over-parametrized case has a smaller second-order error rate O(p−1). It enjoys the same rate as the under-parametrized case in Theorem 4.1. A possible explanation would be that averaging over the randomness in β can partially offset the curse of dimensionality so that RX(β̂,β) achieves the same error rate for all p, n combinations. Remark 4.5. It’s worth mentioning that the only assumption regarding data distribution is assumption (A), where only finite fourth order moment is required. Non-Gaussianity allows our theoretical results more widely applied.

heading: 4.4 DISCUSSION
text: In this section, we first make a short conclusion of what we have done theoretically in this paper and further discuss some possible directions of extension.
We have systematically investigated the second-order fluctuations of two types of prediction risk, RX(β̂,β) and RX,β(β̂,β), for the high dimensional least square estimators β̂. Theorem 4.1 and 4.4 are for RX(β̂,β) while Theorem 4.2 and 4.5 are for RX,β(β̂,β). Both fixed effect and random effect of the regression coefficients β are discussed following the settings in Hastie et al. (2019). RX(β̂,β) and RX,β(β̂,β) are the same when β is nonrandom as established in Theorem 4.3. Asymptotic results are categorized into the under-parametrized case (p < n) and the over-parametrized case (p > n).
Although the first-order limits of the prediction risk in high dimensional linear models have already been well studied in recent years, including general extensions to anisotropic features and signals in Wu & Xu (2020), the “double descent” risk curve is just a function of the limiting ratio lim p/n. There is still a non-negligible discrepancy between the finite-sample prediction risk and its first-order limit on the “double descent” curve. How large is this discrepancy? How fast does the risk converge to its limit? Our CLTs provide answers to such questions and give a fine-grained characterization of the second-order fluctuations of the prediction risks. Not only explicit forms of the leading constants in the limiting means and variances are shown in the main theorems, smaller order terms are also derived to improve the empirical performance for practitioners.
It is also important to recognize the limitations of our results. First, the present paper only concerns linear regression tasks since the linear regression task is simple but important as well. Some recent works linearize neural networks at the initialization and employ Neural Tangent Kernels (Jacot et al. (2018)) to approximate the training procedure of a strongly over-parametrized neural network by solving a linear regression task, e.g. Du et al. (2018); Arora et al. (2019); Lee et al. (2019). Though the setting considered in this paper is simple and limited, the problem has not however been fully understood so far in the literature. Therefore, we are among the first to take the task and develop the second-order fluctuation results for the prediction risk. Second, we assume general covariance Σ and non-Gaussianity for the under-parametrized case, which fits the most updated and realistic settings in the literature, however we only investigate the isotropic settings while still allow for non-Gaussianity under over-parametrization. We haven’t extended it to the more general anisotropic settings yet. The reasons are two-fold. On the one hand, according to Wu & Xu (2020), the first-order limits depend on the Stieltjes transforms of the unknown spectral distribution of Σ. Since Σ is unknown, we cannot obtain any explicit characterization of the first-order limits, not to mention the second-order fluctuations. The CLTs would only be written as certain complicated implicit functions of Σ and would be too abstract to evaluate practically. More restrictions would be imposed on Σ to guarantee the second-order convergence. On the other hand, from the technical perspective, the techniques required for anisotropic over-parametrized cases are very different from the isotropic cases due to difference in the bias-variance decomposition in (4). The tools in random matrix theory have not been fully developed yet for anisotropic cases. Since we have considered various scenarios in this paper, including random and nonrandom signals β for both conditional and unconditional risks, it will take great efforts and continuous work to extend all of them to the most general settings, which would lead to many subsequent works in the field of machine learning and random matrix theory literature.

heading: 5 EXPERIMENTS
text: In this section, we carry out simulation experiments to examine the central limit theorems and the corresponding confidence intervals in Theorem 4.2 and Theorem 4.5. We generate data points from the linear model (1) and directly compute the prediction risk via the bias-variance decomposition in (4). The generative distribution Px is taken to be the standard normal distribution. The noise distribution P is taken to beN(0, 1). In the following, we present the gap between the finite-sample distribution of the prediction risk and the corresponding limiting distribution to check the central limit theorems and use the cover rate to measure the effectiveness of the confidence intervals. More simulation results, including cases with non-Gaussian distributions for Px and P are relegated to the Appendix due to space limitations.
Example 1. This example examines the results in Theorem 4.2. We define a statistic
Tn = p
σc
( RX,β(β̂,β)− σ2
cn 1− cn ) − µc σc .
According to Theorem 4.2, Tn weakly converges to the standard normal distribution as n, p → ∞. In this example, c = 1/2 and p = 50, 100, 200. The finite-sample distribution of Tn is presented by the histogram of Tn in Figure 2 with 1000 repetitions, where the solid blue curve stands for standard normal density function. It can be seen that the finite-sample distribution of Tn is very consistent with the standard normal distribution, especially when n, p become larger. When α = 0.05, the empirical cover rates of the 95%-confidence interval are 94.2%, 93.5% and 95.3% for p = 50, 100 and 200, respectively. All these experiments verify the correctness of our theoretical results.
Example 2. This example verifies the results in Theorem 4.5. Here we define two statistics:
Tn,0 =
√ p
σc,3
{ RX,β(β̂,β)− (1− 1
cn )r2 − σ
2
cn − 1 } − µc,3 σc,3 ,
Tn,1 =
√ p
σ̃c,3
{ RX,β(β̂,β)− (1− 1
cn )r2 − σ
2
cn − 1 } − µ̃c,3 σ̃c,3 .
According to Theorem 4.5, both Tn,0 and Tn,1 weakly converge to the standard normal distribution as n, p→ +∞. Compared to Tn,0, Tn,1 provides a better approximation for the finite sample distribution of RX,β(β̂,β) because it contains smaller order terms in the asymptotic mean and variance. We take c = 3/2 and p = 150, 300, 450. Similarly the finite-sample distributions of Tn,0 and Tn,1 are presented by the histogram of Tn,0 and Tn,1 with 1000 repetitions. The comparison between these two statistics is shown in Figure 3. It can also be seen that the finite sample distributions of Tn,0 and Tn,1 both match the standard normal distribution quite well, especially Tn,1 with more precise characterization. The empirical cover rates of the 95%-confidence interval (11) are 93.8%, 94.7% and 94.4% for p = 150, 300 and 600 respectively, which further shows the validity of our theoretical results.

heading: A PROOF OF THEOREM 4.1 AND THEOREM 4.2
text: Let X = ZΣ1/2. According to the Bai-Yin theorem (Bai & Yin (2008)), the smallest eigenvalue of ZTZ/n is almost surely larger than (1− √ c)2/2 for sufficiently large n. Thus
λmin( 1
n XTX) ≥ c0λmin(
1 n ZTZ) ≥ c0 2 (1− √ c)2,
which implies that the matrix XTX/n is almost surely invertible for large n. By Section 3.2, Π = 0, BX(β̂,β) = BX,β(β̂,β) = 0 and VX(β̂,β) = VX,β(β̂,β). Thus the CLT of RX(β̂,β) is same to that of RX,β(β̂,β). For simplicity, we focus on RX(β̂,β) in the following. Notice that
VX(β̂,β) = σ2
n Tr(Σ̂−1Σ)
= σ2 n Tr ( Σ−1/2 (ZTZ n )−1 Σ−1/2Σ ) = σ2
n p∑ i=1 1 si = σ2p n
1 s dFZ(s)
where FZ is the spectral measure of ZTZ/n. According to Theorem 1 of Hastie et al. (2019), as n, p → ∞ such that p/n = cn → c ∈ (0,∞), FZ(x) weakly converges to the standard MarcenkoPastur law Fc(x) and
VX(β̂,β)→ σ2c
1 s dFc(s) = σ 2 c 1− c .
Here the standard Marcenko-Pastur law Fc(x) has a density function
pc(x) =
{ 1
2πcx
√ (b− x)(x− a), if a ≤ x ≤ b,
0, o.w.,
where a = (1− √ c)2, b = (1 + √ c)2 and pc(x) has a point mass 1− 1c at the origin if c > 1. Hence
RX(β̂,β)− σ2 cn
1− cn =
σ2p
n
1 s dFZ(s)− σ2cn 1 s dFcn(s)
= σ2cn
1
s
( dFZ(s)− dFcn(s) ) .
According to Theorem 1.1 of Bai & Silverstein (2004), p ( RX(β̂,β)− σ2
cn 1− cn
) d−→ N(µc, σ2c ), (12)
where
µc = − σ2c
2πi
γ
1
z
cm(z)3(1 +m(z))−3
{1− cm(z)2(1 +m(z))−2}2 dz (13)
−σ 2c(ν4 − 3)
2πi
γ
1
z
cm(z)3(1 +m(z))−3
1− cm(z)2(1 +m(z))−2 dz,
σ2c = − σ4c2
2π2
C1 C2
1
z1z2
1 (m(z1)−m(z2))2 d dz1 m(z1) d dz2 m(z2)dz1dz2 (14)
−σ 4c3(ν4 − 3)
4π2
C1 C2
1
z1z2
1
(1 +m(z1))2(1 +m(z2))2 dm(z1)dm(z2).
Here the contours in (13) and (14) are closed and taken in the positive direction in the complex plane, enclosing the support of FZ, i.e. [(1− √ c)2, (1 + √ c)2]. The Stieltjes transform m(z) satisfies the equation
z = − 1 m + c 1 +m .
To further simplify the integrations in µc and σc, let z = 1 + √ c(rξ + 1rξ ) + c and perform change of variables, then we have
m(z) = − 1 1 + √ crξ , dz = √ c(r − 1 rξ2 )dξ, dm =
√ cr
(1 + √ crξ)2 dξ
and when ξ moves along the unit circle |ξ| = 1 on the complex plane, z will orbit around the center point 1 + c along an ellipse which enclosing the support of FZ. Thus
µc = − σ2c
2πi
γ
1
z
cm(z)3(1 +m(z))−3
(1− cm(z)2(1 +m(z))−2)2 dz
−σ 2c(ν4 − 3)
2πi
γ
1
z
cm(z)3(1 +m(z))−3
1− cm(z)2(1 +m(z))−2 dz
= σ2c
2πi
|ξ|=1
1
r( √ c+ rξ)(1 + √ crξ)(ξ − 1r )(ξ + 1 r ) dξ
+ σ2c(ν4 − 3)
2πi
|ξ|=1
1
rξ2( √ c+ rξ)(1 + √ crξ) dξ
= σ2c2 (c− 1)2 + σ2c2(ν4 − 3) 1− c .
As for σ2c , note that
1
2πi
γ1
1
z1(m1 −m2)2 dm1
= 1
2πi
|ξ1|=1
1
1 + √ c(r1ξ1 +
1 r1ξ1
) + c ·
√ c r1
(m2 + 1 1+ √ cr1ξ1
)2(1 + √ cr1ξ1)2 dξ1
= 1
2πi
|ξ1|=1
√ c r1ξ1
(ξ1 + √ c
r1 )(r1ξ1
√ c+ 1) {(r1ξ1 √ c+ 1)m2 + 1} 2 dξ1
= c
(c− 1) {(c− 1)m2 − 1} 2 ,
therefore
−σ 4c2
2π2
1
z1z2(m1 −m2)2 dm1dm2
= 2σ4c2
2πi
|ξ2|=1
c
z2(c− 1) {(c− 1)m2 − 1} 2 dm2
= 2σ4c2
2πi
|ξ2|=1
√ c r22ξ2
(c− 1)(1 + √ c r2ξ2)( √ c+ r2ξ2)3
dξ2 = 2c3σ4
(c− 1)4 .
Meanwhile, 1
2πi
γ1
1
z1
1
(1 +m(z1))2 dm(z1)
= 1
2πi
|ξ|=1
1√ cξ(1 + √ crξ)( √ c+ rξ) dξ = 1 c− 1 ,
hence
−σ 4c3(ν4 − 3)
4π2
C1 C2
1
z1z2
1
(1 +m(z1))2(1 +m(z2))2 dm(z1)dm(z2) = σ4c3(ν4 − 3) (1− c)2 ,
and
σ2c = 2c3σ4 (c− 1)4 + σ4c3(ν4 − 3) (1− c)2 .
Let Tn = p
σc
( RX(β̂,β)− σ2
cn 1− cn − µc p
) .
According to (12), we have
P (Lα,c ≤ RX,β(β̂,β) ≤ Uα,c) = P (−Zα/2 ≤ Tn ≤ Zα/2)→ 1− α, where
Lα,c = σ 2 cn
1− cn +
1 p (µc − Zα/2σc),
Uα,c = σ 2 cn
1− cn +
1 p (µc + Zα/2σc).

heading: B PROOF OF THEOREM 4.3
text: Notice that
BX(β̂,β) = β T(Ip − Σ̂+Σ̂)β
= lim z→0+
βT ( Ip − (Σ̂ + zIp)−1Σ̂ ) β
= lim z→0+
zβT(Σ̂ + zIp) −1β.
Since β is a constant vector, we can make use of the results in Theorem 3 in Bai et al. (2007) and Theorem 1.3 in Pan & Zhou (2008) regarding eigenvectors. Their works investigate the sample covariance matrix Ap = T 1/2 p XTpXpT 1/2 p /n, where Tp is an p × p nonnegative definite Hermitian matrix with a square root T 1/2p and Xp is an n× p matrix with i.i.d. entries (xij)n×p. LetUpΛpUTp denote the spectral decomposition of Ap where Λp = diag(λ1, · · · , λp) and Up is a unitary matrix consisting of the orthonormal eigenvectors of Ap. Assume that xp is an arbitrary nonrandom unit vector and y = (y1, y2, · · · , yp)T = UTp xp, two empirical distribution functions based on eigenvectors and eigenvalues are defined as
F Ap 1 (x) = p∑ i=1 |yi|21(λi ≤ x), FAp(x) = 1 p p∑ i=1 1(λi ≤ x).
Then for a bounded continuous function g(x), we have p∑ j=1 |yj |2g(λj)− 1 p p∑ j=1 g(λj) = g(x)dF Ap 1 (x)− g(x)dFAp(x).
The results in Bai et al. (2007) and Pan & Zhou (2008) show that Lemma B.1. (Theorem 3 Bai et al. (2007) and Theorem 1.3 Pan & Zhou (2008)) Suppose that
(1) xij’s are i.i.d. satisfying E(xij) = 0, E(|xij |2) = 1 and E(|xij |4) <∞;
(2) xp ∈ Cp, ‖xp‖ = 1, limn,p→∞ p/n = c ∈ (0,∞);
(3) Tp is nonrandom Hermitian non-negative definite with with its spectral norm bounded in p,
with Hp = FTp d−→ H a proper distribution function and xTp (Tp − zIp)−1xp → mFH (z), where mFH (z) denotes the Stieltjes transform of H(t);
(4) g1, · · · , gk are analytic functions on an open region of the complex plain which contains the real interval[
lim inf p
λmin(Tp)1(0,1)(c)(1− √ c)2, lim sup
p λmax(Tp)1(0,1)(c)(1 +
√ c)2 ] ;
(5) as n, p→∞,
sup z √ n ∥∥∥xTp(mF cn,Hp (z)Tp − Ip)−1xp − 11 + tmF cn,Hp (z)dHn(t) ∥∥∥→ 0. Define Gp(x) = √ n(F
Ap 1 (x)− FAp(x)), then the random vectors(
g1(x)dGp(x), · · · , gk(x)dGp(x) ) forms a tight sequence and converges weakly to a Gaussian vector xg1 , · · · , xgk with mean zero and covariance function
Cov(xg1 , xg2) = − 1
2π2
C1 C2 g1(z1)g2(z2)
(z2m2 − z1m1)2
c2z1z2(z2 − z1)(m2 −m1) dz1dz2.
The contours C1, C2 are disjoint, both contained in the analytic region for the functions (g1, · · · , gk) and enclose the support of F cn,Hp for all large p.
(6) If H(x) satisfies dH(t)
(1 + tm(z1))(1 + tm(z2)) =
dH(t)
1 + tm(z1)
dH(t)
1 + tm(z2) ,
then the covariance function can be further simplified to
Cov(xg1 , xg2) = 2
c ( g1(x)g2(x)dF
c,H(x)− g1(x)dF c,H(x) g2(x)dF c,H(x)).
Recall that BX(β̂,β) = limz→0+ zβT(Σ̂ + zIp)−1β. Let g(x) = 1/(x+ z) and xp = β/r. Then we have
g(x)dGn(x) = √ n ( 1 r2 βT(Σ̂ + zIp) −1β − g(x)dFcn(x) ) ,
where Fcn(x) is the standard Marcenko-Pastur law. It is not difficult to check that under Assumptions (A1), (B1) and (C1), all the conditions (1)-(6) in Lemma B.1 are satisfied. To proceed further, denote a = (1− √ c)2, b = (1+ √ c)2. If c is replaced by cn, a and b are denoted by an and bn respectively. By some algebraic calculations, we have g(x)dFcn(x) = (1− 1
cn ) · 1 z + bn an 1 x+ z · 1 2πcnx
√ (bn − x)(x− an)dx
= (1− 1 cn ) · 1 z − −1 + cn + z −
√ c2n + 2cn(z − 1) + (1 + z)2
2cnz ,
and
Var(xg) = 2
c
( {g(x)}2dFc(x)− { g(x)dFc(x) }2) = 2
c
( (1− 1
c ) · 1 z2 + b a
1 (x+ z)2 · 1 2πcx
√ (b− x)(x− a)dx ) −2 c ( (1− 1 c ) · 1 z + b a 1 x+ z · 1 2πcx √ (b− x)(x− a)dx )2 .
Therefore,
lim z→0+ z g(x)dFcn(x) = 1−
1
cn and lim z→0+ z2Var(xg) = 2(c− 1) c3 .
Furthermore, as n, p→∞, p/n = cn → c > 1, √ n ( BX(β̂,β)− (1− 1 cn )r2 ) d−→ N ( 0, 2(c− 1) c3 r4 ) . This can be rewritten as √ p ( BX(β̂,β)− (1− 1 cn )r2 ) d−→ N ( 0, 2(c− 1) c2 r4 ) .
Next we deal with the variance term VX(β̂,β). According to the Assumption (B1), the variance term is
VX(β̂,β) = σ2 n Tr{Σ̂+} = σ 2 n n∑ i=1 1 si ,
where si, i = 1, . . . , n are the nonzero eigenvalues of XTX/n. Let {ti, i = 1, . . . n} denote the non-zero eigenvalues of XXT/p, then we have
VX(β̂,β) = σ2
p n∑ i=1 1 ti = σ2n p
1 t dFXXT/p(t)→
σ2
c− 1 .
By interchanging the role of p and n, from the result in Theorem 4.1, as n, p → ∞, p/n = cn → c > 1, we have,
n∑ i=1 1 ti − n 1− c′n d−→ N ( c′ (c′ − 1)2 + c′(ν4 − 3) 1− c′ ,
2c′ (c′ − 1)4 + c′(ν4 − 3) (1− c′)2
) .
where c′n = n/p = 1/cn, c ′ = 1/c. This result can be rewritten as
n∑ i=1 1 ti − p cn − 1 d−→ N ( c (1− c)2 + (ν4 − 3) c− 1 , 2c3 (1− c)4 + c(ν4 − 3) (c− 1)2 ) .
Hence the CLT of VX(β̂,β) is given by p ( VX(β̂,β)− σ2
cn − 1
) d−→ N ( cσ2 (1− c)2 + σ2(ν4 − 3) c− 1 , 2c3σ4 (1− c)4 + cσ4(ν4 − 3) (c− 1)2 ) .
Notice that Cov ( BX(β̂,β), VX(β̂,β) ) = 0. According to the consistency rate and the limiting
distribution of BX(β̂,β) and VX(β̂,β), we know that the bias BX(β̂,β) is the leading term of RX(β̂,β). This implies that
√ p { RX(β̂,β)− (1− 1
cn )‖β‖22 −
σ2
cn − 1
} d−→ N ( 0, σ2c,1 ) ,
where σ2c,1 = 2(c− 1)r4/c2. A practical version of this CLT is given by √ p { RX(β̂,β)− (1− 1
cn )‖β‖22 −
σ2
cn − 1
} d−→ N ( µ̃c,1, σ̃ 2 c,1 ) ,
where
µ̃c,1 = 1 √ p { cσ2 (1− c)2 + σ2(ν4 − 3) c− 1 } ,
σ̃2c,1 = 2(c− 1) c2 r4 + 1 p { 2c3σ4 (1− c)4 + cσ4(ν4 − 3) (c− 1)2 } .

heading: C PROOF OF THEOREM 4.4
text: First we consider the bias term BX(β̂,β). By Assumption (A1), (B1), and (C2),
BX(β̂,β) = E[βTΠΣΠβ|X] = E[βTΠβ|X] = Tr { (Ip − Σ̂+Σ̂)E(ββT|X) }
= r2
p Tr{Ip − Σ̂+Σ̂} = r2(1− n/p).
Alternatively, we can rewrite the bias as
BX(β̂,β) = lim z→0+
E[βT(Ip − (Σ̂ + zIp)−1Σ̂ ) β|X]
= lim z→0+
E[zβT(Σ̂ + zIp)−1β|X]
= lim z→0+
z r2
p Tr(Σ̂ + zIp)
−1.
Define that fn(z) = z r 2 p Tr(Σ̂+zIp) −1. Notice that |fn(z)| and |f ′n(z)| are bounded above. By the Arzela-Ascoli theorem, we deduce that fn(z) converges uniformly to its limit. Under Assumption (C2), by the Moore-Osgood theorem, almost surely,
lim n,p→∞ BX(β̂,β) = lim z→0+ lim n,p→∞
z r2
p Tr(Σ̂ + zIp)
−1
= lim z→0+ lim n,p→∞
z r2
p Tr
( 1
n XXT + zIn
)−1 ,
In fact, lim
n,p→∞ BX(β̂,β) = r 2 lim z→0+ lim n,p→∞ zmn(−z),
where mn(z) is the Stieltjes transform of empirical spectral distribution of Σ̂ = XTX/n. According to Theorem 2.1 in Zheng et al. (2015) and Lemma 1.1 in Bai & Silverstein (2004), the truncated version of p(mn(z)−m(z)) converges weakly to a two-dimensional Gaussian process M(·) satisfying
E[M(z)] = cm3(1 +m)
{(1 +m)2 − cm2}2 +
c(ν4 − 3)m3
(1 +m) {(1 +m)2 − cm2} ,
and
Cov ( M(z1),M(z2) ) = 2 { m′(z1)m′(z2) (m(z1)−m(z2))2 − 1 (z1 − z2)2 } +
c(ν4 − 3)m′(z1)m′(z2) (1 +m(z1))2(1 +m(z2))2 ,
where m = m(z) represents the Stieltjes transform of limiting spectral distribution of companion matrix XXT/n satisfying the equation
z = − 1 m + c 1 +m , m(z) = −1− c z + cm(z).
When p > n, we can actually solve m(z) equation and obtain that
m(z) = −1 + c− z + √ −4z + (1− c+ z)2 2z ,
m(z) = 1− c− z + √ −4z + (1− c+ z)2 2cz .
Therefore, by some algebraic calculations, we have
lim n,p→∞ BX(β̂,β) = lim n,p→∞ r2 lim z→0+ zmn(−z) = r2 lim z→0+
{ zm(−z) + z(1− 1
c ) 1 z } = lim
n,p→∞ r2 lim z→0+ z n p mn(z) = r 2 1 c lim z→0+ zm(−z)
= r2(1− 1 c ).
Moreover, Var ( M(z) ) = lim
z1→z2=z Cov
( M(z1),M(z2) ) = 2m′(z)m′′′(z)− 3(m′′(z))2
6(m′(z))2 + c(ν4 − 3)(m′(z))2 (1 +m(z))4 .
By substituting of the explicit form of m(z), we can easily derive that
lim z→0+ zE[M(−z)] = 0, lim z→0+ z2Var(M(−z)) = 0,
which means that the second-order limit of BX(β̂,β) is still r2(1 − 1/c). All in all, BX(β̂,β) is identical with a constant r2(1− 1/c) in distribution. On the other hand, by Assumption (B1),
VX(β̂,β) = σ2 n Tr{Σ̂+} = σ 2 n n∑ i=1 1 si ,
where si, i = 1, . . . , n are the nonzero eigenvalues of XTX/n. Similar to the proof of Theorem 4.3, the CLT of VX(β̂,β) is given by
p ( VX(β̂,β)− σ2
cn − 1
) d−→ N ( cσ2 (1− c)2 + σ2(ν4 − 3) c− 1 , 2c3σ4 (1− c)4 + cσ4(ν4 − 3) (c− 1)2 ) .
Combining the results of BX(β̂,β) and VX(β̂,β), we have p { RX(β̂,β)− r2(1− 1
cn )− σ
2
cn − 1
} d−→ N(µc,2, σ2c,2),
where
µc,2 = cσ2 (1− c)2 + σ2(ν4 − 3) c− 1 , σ2c,2 = 2c3σ4 (1− c)4 + cσ4(ν4 − 3) (c− 1)2 .

heading: D PROOF OF THEOREM 4.5
text: Note that under Assumption (B1) and (C2), BX,β(β̂,β) = βTΠβ = βT(Ip − Σ̂+Σ̂)β. If we directly consider βT(Ip − Σ̂+Σ̂)β, we can make use of the asymptotic results for quadratic forms Theorem 7.2 in Bai & Yao (2008) stated as follows. Lemma D.1. (Theorem 7.2 in Bai & Yao (2008)) Let {An = [aij(n)]} be a sequence of n × n real symmetric matrices, {xi}i∈N be a sequence of i.i.d. K dimensional real random vectors, with E(xi) = 0, E(xixTi ) = (γij)K×K and E[‖xi‖4] <∞. Denote
xi = (x`i)K×1, X(`) = (x`1, · · · , x`n)T, ` = 1, · · · ,K, i = 1, · · · , n, assume the following limits exist
ω = lim n→∞
1
n n∑ i=1 a2ii(n), θ = lim n→∞ 1 n TrA2n.
Then the K-dimensional random vectors
zn = (zn,`)K×1, zn,` = 1√ n
( X(`)TAnX(`)− γ`` Tr{An} ) , 1 ≤ ` ≤ K,
converge weakly to a zero-mean Gaussian vector with covariance matrixD = D1 +D2 where [D1]``′ = ω { E(x2`1x2`′1)− γ``γ`′`′ } , [D2]``′ = (θ − ω)(γ``′γ`′` + γ2``′), 1 ≤ `, `′ ≤ K.
According to the results in Lemma D.1, letAn = Π = Ip − Σ̂+Σ̂, then we have, as p→∞, √ p { βTΠβ − r 2
p Tr(Π)
} d−→ N(0, d2 = d21 + d22),
where
ω = lim p→∞
1
p p∑ i=1 Π2ii, θ = lim p→∞ 1 p Tr(Π2) = 1− 1 c ,
and
d21 = ω { E(x2`1x2`1)− γ2`` } = ω (p2 r4 E(β4i )− 1 ) r4,
d22 = (θ − ω)(γ2`` + γ2``) = 2(θ − ω)r4. Since in the proof of Theorem 4.4, we have already shown that
r2
p Tr(Π) = r2(1− n p ).
In particular, if β follows multivariate Gaussian distribution, i.e. β ∼ Np(0, r 2 p Ip), then as p→∞,
√ p { BX,β(β̂,β)− r2(1− n p ) } d−→ N ( 0, 2(1− 1 c )r4 ) .
Moreover, VX,β(β̂,β) = VX(β̂,β), we have already proved in Theorem 4.4 that
p(VX,β(β̂,β)− σ2
cn − 1 ) d−→ N ( cσ2 (1− c)2 + σ2(ν4 − 3) c− 1 , 2c3σ4 (1− c)4 + cσ4(ν4 − 3) (c− 1)2 ) .
Note that Cov(BX,β(β̂,β), VX,β(β̂,β)) = 0. According to the consistency rate of BX,β(β̂,β) and VX,β(β̂,β), we know that the bias BX(β̂,β) is the leading term of RX,β(β̂,β). This implies that
√ p { RX,β(β̂,β)− r2(1− 1
cn )− σ
2
cn − 1
} d−→ N(0, σ2c,3),
where σ2c,3 = 2r 4(1− 1/c). A practical version of this CLT is given by
√ p { RX,β(β̂,β)− r2(1− 1
cn )− σ
2
cn − 1
} d−→ N(µ̃c,3, σ̃2c,3),
where
µ̃c,3 = 1 √ p { cσ2 (1− c)2 + σ2(ν4 − 3) c− 1 } ,
σ̃2c,3 = 2(1− 1
c )r4 +
1
p { 2c3σ4 (1− c)4 + cσ4(ν4 − 3) (c− 1)2 } .

heading: E MORE EXPERIMENTS
text: E.1 MORE RESULTS OF EXAMPLE 1
This example checks Theorem 4.2. We define a statistic
Tn = p
σc
( RX(β̂,β)− σ2
cn 1− cn ) − µc σc .
According to Theorem 4.2, Tn weakly converges to the standard normal distribution as n, p→∞. In this example, c = 1/2 and p = 50, 100, 200. To make sure the assumption (A) holds, the generative distribution Px is taken to be the standard normal distribution, the centered gamma with shape 4.0 and scale 0.5, and the normalized Student-t distribution with 6.0 degree of freedom. The finitesample distribution of Tn is estimated by the histogram of Tn under 1000 repetitions. The results are presented in Figure 4. One can find that the finite-sample distribution of Tn tends to the standard normal distribution as n, p→ +∞.Whenα = 0.05, the empirical cover rates of the 95%-confidence interval are reported in Figure 5.
E.2 MORE RESULTS OF EXAMPLE 2
The Example 2 checks Theorem 4.5. Here we consider the standardized statistics:
Tn,0 =
√ p
σc,3
{ RX(β̂,β)− (1− 1
cn )r2 − σ
2
cn − 1 } − µc,3 σc,3 ,
Tn,1 =
√ p
σ̃c,3
{ RX(β̂,β)− (1− 1
cn )r2 − σ
2
cn − 1 } − µ̃c,3 σ̃c,3 .
According to the central limit theorem (10) and its practical version, both Tn,0 and Tn,1 weakly converge to the standard normal distribution as n, p→ +∞. We take c = 2 and p = 100, 200, 400. The finite-sample distributions of Tn,0 and Tn,1 are estimated by the histogram of Tn,0 and Tn,1 under 1000 repetitions. The results are presented in Figure 6 and Figure 7. When α = 0.05, the empirical cover rates of the 95%-confidence interval (11) are reported in Figure 8.
E.3 EXAMPLE 3
This example checks Theorem 4.3. To proceed further, we denote two statistics:
Tn,2 =
√ p
σc,1
{ RX(β̂,β)− (1− 1
cn )r2 − σ
2
cn − 1 } − µc,1 σc,1 ,
Tn,3 =
√ p
σ̃c,1
{ RX(β̂,β)− (1− 1
cn )r2 − σ
2
cn − 1 } − µ̃c,1 σ̃c,1 .
According to the central limit theorem (8) and its practical version, both Tn,2 and Tn,3 weakly converge to the standard normal distribution as n, p→ +∞. We take c = 2 and p = 100, 200, 400. The finite-sample distributions of Tn,2 and Tn,3 are estimated by the histogram of Tn,2 and Tn,3 under 1000 repetitions. The results are presented at Figure 9 and Figure 10. One can see that the finite-sample distributions of Tn,2 and Tn,3 are close to the standard normal distribution, and the finite-sample performance of Tn,3 is better than that of Tn,2. When α = 0.05, the empirical cover rates of the 95%-confidence interval (9) are reported in Figure 11.

heading: F AN EXAMPLE FROM FIGURE 1

heading: G AN ANISOTROPIC EXAMPLE FOR REMARK 4.2
text: In the over-parameterized case, the bias termBX(β̂,β) = βTΠΣΠβ is non-zero while the variance term VX(β̂,β) remains the same as under-parameterized case. Therefore in this section, we conduct a small simulation to examine the fluctuation of the bias BX for both isotropic and anisotropic Σ in the over-parameterized case with non-random β satisfying Assumption (C1). In particular, in the following we set r = 1.
We consider both localized and delocalized β such that
1. Localized case: β1 = (1, 0, · · · , 0); 2. Delocalized case: β2 = 1√p (1, · · · , 1);
and both the isotropic and anisotropic Σ
3. Identity case: Σ1 = Ip; 4. Compound symmetric case: Σ2 = 0.5Ip + 0.51p1Tp .
Then we fix p/n = 2 and let p vary from 10 to 300, we present in Figure 13 the empirical variance of √ p ∗BX and p ∗BX under various combinations of Σ and β with 1000 replications. From the plot on the top left panel in Figure 13, we can see that the variance of √ p ∗ BX for both β1 and β2 remain constant as p grows, which indicates that the convergence rate of BX is 1/ √ p under the isotropic case regardless of localized or delocalized β. As for the anisotropic case on the top right corner, the variance of √ p∗BX stabilizes for β1, while decays for β2, which indicates that convergence rate of BX under (Σ2,β2) and (Σ2,β1) are different.
This simulation result further confirms our conjecture that in the over-parameterized case, there is no universal CLT for the prediction risk RX(β̂,β) under the anisotropic setting for non-random β.

Review summaries:
In this article, the authors characterized the second-order fluctuation of the prediction risk of the (min-norm) least square estimator, by assuming an underlying noisy teacher model $y_i = \beta^T x_i + \epsilon_i$, in the regime where the data dimension $p$ and the number of training samples $n$ grow large at the same pace. Results in both the under-parameterized (Theorem 4.1 and 4.2) and over-parameterized regimes (Theorem 4.3-4.5) were provided, under the statistical model where the data $x_i$ are zero-mean random vectors with **generic** i.i.d. entries and then "rotated" to have some possible covariance structures. Numerical experiments for relatively small values of $n,p$ were conducted to support the theoretical assessment.

This paper investigates the phenomenon of double descent, also referred to as "more data hurts", in high dimensional linear regression using the least square estimator.  In the same setup, previous sharp results were already established in the asymptotic regime. Non-asymptotic results are also known but are less precise. The authors of this paper try to provide a new type of results that fill the gap between the two regimes (asymptotic vs non-asymptotic). To do so they have managed to derive second order (CLT type) asymptotic results for different risks based on more refined random matrix theory results.