Year: 2021

Abstract:
In this paper, we study the generalization risk of ridge and ridgeless linear regression. We assume that the data features follow a multivariate normal distribution and that the spectrum of the covariance matrix consists of a given set of eigenvalues of proportionally growing multiplicity. We characterize the limiting bias and variance when the dimension and the number of training samples tend to infinity proportionally. Exact formulae for the bias and variance are derived using the random matrix theory and convex Gaussian min-max theorem. Based on these formulae, we study the sample-wise multiple descent phenomenon of the generalization risk curve, i.e., with more data, the generalization risk can be non-monotone, and specifically, can increase and then decrease multiple times with more training data samples. We prove that sample-wise multiple descent occurs when the spectrum of the covariance matrix is highly ill-conditioned. We also present numerical results to confirm the values of the bias and variance predicted by our theory and illustrate the multiple descent of the generalization risk curve. Moreover, we theoretically show that the ridge estimator with optimal regularization can result in a monotone generalization risk curve and thereby eliminate multiple descent under some assumptions.

1 INTRODUCTION
The double/multiple descent phenomenon attracted recent research attention due to (Belkin et al., 2019). This line of work focuses on the parameter-wise double/multiple descent phenomenon of the risk curve (Bartlett et al., 2020; Tsigler & Bartlett, 2020; Belkin et al., 2019; 2020; Chen et al., 2020a; Liang et al., 2020; Advani et al., 2020; Bös & Opper, 1998; Krogh & Hertz, 1992; Le Cun et al., 1991; Mei & Montanari, 2019; Opper et al., 1990; Vallet et al., 1989; Watkin et al., 1993). The classical learning theory shows that when the number of parameters (which reflects the model complexity) increases, the test error (generalization risk) first decreases due to more fitting power, and then increases due to overfitting. The generalization risk attains a peak at the interpolation threshold (the number of parameters equals the number of data points so that the model interpolates the data). This results in a U-shaped risk curve if we plot the test error versus the number of parameters. The double descent risk curve posits that the risk will decrease (again) if one further increases the model complexity beyond the interpolation threshold (Belkin et al., 2019). Thus there is a second descent in addition to the first one in the U-shaped stage of the curve. Belkin et al. (2019) presented empirical results and showed the existence of such double descent behavior in the random Fourier features model, the fully connected neural network, and the random forest model. Prior to (Belkin et al., 2019), earlier studies of the shape and features of the risk curve in a number of contexts include (Vallet et al., 1989; Opper et al., 1990; Le Cun et al., 1991; Krogh & Hertz, 1992; Bös & Opper, 1998; Watkin et al., 1993; Advani et al., 2020). Loog et al. (2020) presented a prehistory of the double descent phenomenon. Belkin et al. (2020) proved the double descent curve in the Gaussian model and the Fourier series model. Mei & Montanari (2019) theoretically established the double descent curve of the random features regression. Bartlett et al. (2020); Tsigler & Bartlett (2020) characterized the conditions for ridgeless and ridge linear regression problems, respectively, under which the minimum-norm interpolants achieve near-optimal generalization risk. Liang et al. (2020) showed that the test error of the minimum-norm interpolator of data in reproducing kernel Hilbert space is upper bounded by a multiple descent curve as the model complexity increases. They also presented a numerical result supporting that the test error itself exhibits a multiple descent curve.
Chen et al. (2020a) proved that the multiple descent curve does exist for the minimum-norm interpolator in linear regression and that the curve can be even designed.
Following the parameter-wise double descent, research interest extended to epoch-wise and samplewise double descent (Nakkiran et al., 2020; Chen et al., 2020b; Min et al., 2021; Nakkiran et al., 2021). Nakkiran et al. (2020) observed from their numerical result that the generalization risk experiences a double descent as one keeps the model size fixed and increases the training time. They called this observation epoch-wise double descent. Nakkiran et al. (2020) also noted sample-wise non-monotonicity, which means that more data can hurt generalization. Nakkiran et al. (2021) proved that for isotropic features, optimally regularized ridge regression yields a monotonic generalization risk curve with more samples. Nakkiran et al. (2021) also showed that if the features are formed by projecting high-dimensional isotropic data to a random low-dimensional space (say, ddimensional), the optimally regularized ridge regression has a monotonic generalization risk curve with increasing d (the model size). Sample-wise non-monotonicity and double descent was also observed in (Chen et al., 2020b; Min et al., 2021) in adversarially trained models. C ompared to (Wu & Xu, 2020; ichi Amari et al., 2021; Dobriban & Wager, 2018; Richards et al., 2021), in what follows, we highlight our contributions and the differences from them. First, our major contribution is providing a rigorous proof for the existence of sample-wise (test error vs. the number of training samples) double and multiple descent in linear regression. However, (Richards et al., 2021) only mentioned parameter-wise double descent (test error vs. model capacity) in their related work. (ichi Amari et al., 2021) only mentioned epoch-wise (test error vs. training time) double descent in Appendix A.2. Neither (Richards et al., 2021) nor (ichi Amari et al., 2021) mentioned multiple descent. Second, we made and theoretically proved the observation that an ill-conditioned covariance matrix is a sufficient condition for the existence of sample-wise multiple descent. To the best of our knowledge, our work is the first paper that pointed this out. Third, we solved the Stieltjes transform explicitly and derived explicit formulae for the risk and variance in our setup. In addition, we also provided rigorous treatment to the ridgeless setting and also obtained explicit formulae for it. Fourth, there is another difference between our paper and the papers that the reviewer mentioned. (Wu & Xu, 2020; ichi Amari et al., 2021; Dobriban & Wager, 2018; Richards et al., 2021) assumed a prior on the true linear model and takes expectation over the prior. In our paper, we do not assume a prior on the true linear model and our risk does not take the expectation over a random true linear model.
In the setting of generally anisotropic features, this paper gives an asymptotic characterization of the generalization risk curve with more samples. The asymptotic regime is an approximation for large n, d and can also shed light on practical machine learning problems. We first introduce our problem setup.
1.1 PROBLEM SETUP
Data Distribution Let Σ ∈ Rd×d be a positive semi-definite matrix which is termed the covariance matrix, and let θ∗ ∈ Rd. The eigenvalues of Σ are λ1, . . . , λm with multiplicity d1, . . . , dm, respectively. We have d = ∑m i=1 di. Assume that λ1, . . . , λm are fixed, distinct, all positive, and do not depend on d (i.e., for all d, the eigenvalue of Σ are always λ1, . . . , λm). We assume the following data distribution D for (x, y) ∈ Rd × R:
x ∼ N (0,Σ) , y = x>θ∗ + ,
where x and are independent and ∼ N (0, σ2). In practice, there are natural random variables x that satisfy our assumption. For example, assume that we want to use machine A to measure the length of several objects and use machine B to measure their temperature. The measured lengths and temperatures follow an i.i.d. Gaussian distribution. However, the variance of measurement of machine A is different from that of machine B. Then we consider the random vector formed by the measurements x = (l1, . . . , ln, t1, . . . , tn), where li and ti are the length and temperature of object i, respectively. This results in a block-structured covariance matrix. When we measure more objects, the size of the covariance matrix tends to infinity. Second, the motivation came from (Nakkiran et al., 2021). (Nakkiran et al., 2021) observed empirically in their Figure 2 that when the covariance matrix has a block structure (specifically, there are only two fixed different eigenvalues 10 and 1), the expected excess risk exhibits multiple descent. We quantitatively studied this observation and
obtained the related formulae. The excess risk of an estimator θ ∈ Rd is given by R(θ) = Ex,y∼D [( y − x>θ )2 − (y − x>θ∗)2] . Assume that the training data {(xi, yi)}ni=1 ⊆ Rd × R is drawn i.i.d. from D. Write
X =  x > 1 ... x>n  ∈ Rn×d , y =  y1... yn  ∈ Rn . (1) We have y = Xθ∗ + , where ∼ N (0, σ2In).
Ridge Estimator and Minimum-Norm Estimator
Definition 1 (Ridge estimator). The ridge estimator θ̂λ,n,d ∈ Rd (λ > 0) solves the following minimization problem
min θ∈Rd
1 n ‖Xθ − y‖22 + λ ‖θ‖ 2 2 .
Definition 2 (Minimum-norm estimator). The minimum-norm estimator (also known as the ridgeless estimator) θ̂0,n,d ∈ Rd solves the following minimization problem
min θ∈Rd ‖θ‖2 such that ‖Xθ − y‖2 = min θ∈Rd ‖Xθ − y‖2 .
We are interested in the expected excess risk of θ̂λ,n,d, which is given by Rλ,n,d = E [ R ( θ̂λ,n,d )] .
The expectation is taken over the randomness of the training data {(xi, yi)}ni=1.
Asymptotic Regime Let Πi ∈ Rd×d be the orthogonal projection to the eigenspace of λi. This paper focuses on the asymptotic behavior of the expected excess risk of θ̂λ,n,d where n, di → +∞, di/n→ zi (zi is a fixed positive constant), and ‖Πiθ∗‖2 → ηi. In other words, we are interested in
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi Rλ,n,d .
1.2 OUR CONTRIBUTIONS
Our contributions are summarized as follows.
1. We obtain the formulae for the limiting bias and variance, and thereby the limiting risk. We use two methods to obtain these formulae. Specifically, we obtain the limiting bias and variance by solving the Stieltjes transform and computing its derivatives and antiderivatives. We also use convex Gaussian min-max theorem (CGMT) (Thrampoulidis et al., 2015) to compute the limiting variance. The advantage of the CGMT method is that it is more mathematically tractable for the ridgeless estimator. Through the CGMT approach, we obtain a closed-form formula for the variance in the underparameterized regime and simplify the formula for the variance in the overparameterized regime. Moreover, based on the simplified formula, we deduce a closed-form expression for the variance if the covariance matrix of the data distribution has two different eigenvalues.
2. We find and theoretically prove that sample-wise multiple descent happens when the covariance matrix has eigenvalues of very different orders of magnitude (thus the covariance matrix is highly ill-conditioned).
3. We show that if the true linear model θ∗ satisfies ‖Πiθ∗‖2 = √ di d , optimal regulariza-
tion (i.e., pick λ that minimizes the generalization risk of θ̂λ,n,d) results in a monotone generalization risk curve—in other words, with optimal regularization, more data samples
always improve generalization. Thus there is no sample-wise double or multiple descent. This provides a theoretical proof of a phenomenon observed in (Nakkiran et al., 2021) that optimal regularization can mitigate double descent for anisotropic data. Note that without regularization, there will be a blow-up in expected excess risk when n = d (the linear model exactly interpolates the data) and therefore, there is no samplewise descent across the under- and over-parameterized regimes.
2 PRELIMINARIES
Notation Write [m] for {1, 2, . . . ,m}. Let i denote the imaginary unit. If x ∈ Rn and Σ ∈ Rn×n is a positive semidefinite matrix, write ‖x‖Σ , √ x>Σx. For a vector x, let ‖·‖1 and ‖·‖2 denote the `1 and `2 norm, respectively. Let denote the Hadamard (entry-wise) product between vectors. Write ‖ · ‖2 and ‖ · ‖F for the spectral matrix norm and Frobenius matrix norm, respectively. Let 4 denotes the Loewner order. For two square matrices A and B of the same size, write A 4 B if B − A is positive semidefinite. Define spec (A) as the set of all eigenvalues of A. Let O(d) ={ A ∈ Rd×d | AA> = A>A = Id } denote the set of d × d orthogonal matrices. Define Sd−1(r) , {x ∈ Rd | ‖x‖2 = r}. Denote almost sure convergence by a.s.→, and convergence in probability plim and P→.
Ridge Estimator and Minimum-Norm Estimator We begin with the equivalent characterizations of the ridge and minimum-norm estimator. An equivalent characterization of the ridge estimator θ̂λ,n,d is
θ̂λ,n,d = ( X>X + λnId )−1 X>y = X> ( λnIn +XX >)−1 y . (2) The second equality in Equation (2) is because of the Sherman–Morrison–Woodbury formula. A proof of Equation (2) can be found in (Tsigler & Bartlett, 2020).
An equivalent definition of the minimum-norm estimator θ̂0,n,d is that θ̂0,n,d solves the following minimization problem
min θ∈Rd
‖θ‖2 such that X >Xθ = X>y .
Thus we have θ̂0,n,d = ( X>X )+ X>y = X> ( XX> )+ y = X+y ,
where A+ denotes the pseudo-inverse of A. The second and third equalities are because of the identity X+ = ( X>X )+ X> = X> ( XX> )+ . The minimum-norm estimator is the limit of the
ridge estimator θ̂λ,n,d as λ→ 0+: θ̂0,n,d = lim
λ→0+ θ̂λ,n,d . This is because of the identity limλ→0+ ( X>X + λnId )−1 X> =
limλ→0+ X > (λnIn +XX>)−1 = X+.
Bias-Variance Decomposition of Expected Excess Risk We first show that the excess risk of an estimator θ equals the norm of θ − θ∗:
R(θ) = E(x,y)∼D [( y − x>θ )2 − (y − x>θ∗)] = Ex [(x> (θ∗ − θ))2] = E [ (θ∗ − θ)> Σ (θ∗ − θ) ] = E [ ‖θ∗ − θ‖2Σ ] .
For the ridge estimator, the expected excess risk is Rλ,d,n =E [ ‖θ∗ −X>(nλIn +XX>)−1(Xθ∗ + )‖2Σ ] =E [ ‖(Id −X>(nλIn +XX>)−1X)θ∗ −X>(nλIn +XX>)−1 ‖2Σ
] =E [ ‖(Id −X>(nλIn +XX>)−1X)θ∗‖2Σ ] + E
[∥∥X>(nλIn +XX>)−1 ∥∥2Σ] =E [ ‖(Id −X>(nλIn +XX>)−1X)θ∗‖2Σ ] + σ2E tr [ XΣX>(nλIn +XX >)−2 ]
,Bλ,d,n + Vλ,d,n . (3)
For the minimum-norm estimator, the expected excess risk is R0,d,n = E [ ‖θ∗ −X+(Xθ∗ + )‖2Σ ] = E
[∥∥(Id −X+X) θ∗ −X+ ∥∥2Σ] = E
[∥∥(Id −X+X) θ∗∥∥2Σ]+ E [∥∥X+ ∥∥2Σ] = E
[∥∥(Id −X+X) θ∗∥∥2Σ]+ σ2E tr [(X+)>ΣX+] , B0,d,n + V0,d,n . (4)
We call Bλ,d,n and B0,d,n the bias term, and call Vλ,d,n and V0,d,n the variance term. The bias and variance for the minimum-norm estimator are the limit of their counterpart for the ridge estimator as λ → 0+, i.e., limλ→0+ Bλ,d,n = B0,d,n and limλ→0+ Vλ,d,n = V0,d,n (this can be shown by Lebesgue’s dominated convergence theorem, see our proof in Lemma 5 and Lemma 6, respectively).
3 MAIN RESULTS
3.1 LIMITING RISK AND SAMPLE-WISE MULTIPLE DESCENT
We study the limiting bias and variance for a linear regression problem in which the data distribution follows a multivariate normal distribution, the spectrum of the covariance matrix exhibits a block structure and tends to a discrete distribution. Thanks to the random matrix theory, we obtain the formulae (presented in Theorem 1) for the limiting bias and variance, and thereby the total risk.
We use two methods to obtain these formulae. The first method is through the Stieltjes transform of the matrix 1nXX
>. The central quantity for computing the limiting bias and variance through the first method is the solution ρ∗ to the optimization problem Equation (5) in Item 1 of Theorem 1. Item 1 guarantees the existence of a solution and determines its optimality condition Equation (6). Item 2 computes the Jacobian matrix of ρ∗ with respect to λi and provides a closed-form formula to compute the Jacobian matrix. Equation (9) and Equation (10) in Item 4 give the formulae for the limiting bias obtained by the first method. Equation (11) and Equation (12) give the limiting variance.
The second method is through the convex Gaussian min-max theorem (CGMT) (Thrampoulidis et al., 2015). The central quantity is the solution r∗ to the minimax optimization problem Equation (8) in Item 3. We use CGMT to obtain the formulae for the variance term. They are presented in Equation (13) and Equation (14) in Item 4. Theorem 1. The following statements hold:
1. There exists a minimizer ρ ∈ Rm+ that solves
inf ρ∈Rm+
log λ+ m∑
j=1
λjρj + m∑ j=1 ( ρj − zj(log ρj zj + 1) ) . (5) The minimizer ρ∗ satisfies
λi λ+ ∑m j=1 λjρ ∗ j + 1− zi ρi = 0 , ∀i ∈ [m] . (6)
2. Let ρ∗ ∈ Rm be a minimizer of Equation (5) and J = ∂ρ ∗
∂λ ∈ R m×m be the Jacobian
matrix Jij = ∂ρ∗i ∂λj . Then J is given by
J = ( diag (λ) + ( λ+ λ>ρ∗ ) Im − (z− ρ∗)λ> )−1 ( (z− ρ∗) ρ∗> − diag (ρ∗) ) and the matrix ( diag (λ) + ( λ+ λ>ρ∗ ) Im − (z− ρ∗)λ> ) is always invertible.
3. Define r = (r1, . . . , rm), λ = (λ1, . . . , λm), and
ϑ(rt, r,λ) = 2rt
√ 1 + ∑ i∈[m] r2i − 2rt ∑ i∈[m] √ ziri + ∑ i∈[m] 1 λi r2i − λr2t . (7)
For any Kt ≥ 2λ and Ku ≥ 2λ+(2+
√ γ)
λ , we have
max 0≤rt≤Kt min 0≤ri≤Ku ϑ(rt, r,λ) = min 0≤ri≤Ku max 0≤rt≤Kt ϑ(rt, r,λ) = max rt≥0 min ri≥0 ϑ(rt, r,λ) = min ri≥0 max rt≥0 ϑ(rt, r,λ) (8) and the above optimization problem has a solution.
4. Let r∗ = (r∗1 , . . . , r ∗ m) solve Equation (8). Define q = ( η21/z1, . . . , η 2 m/zm )> and view
λ = (λ1, . . . , λm) > as a column vector. The limiting bias is given by
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi
Bλ,d,n = q >(λ ρ∗ + Jλ 2) , (9)
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi B0,d,n = lim λ→0+
q>(λ ρ∗ + Jλ 2) . (10)
The limiting variance is given by
lim n,di→+∞ di/n→zi
Vλ,d,n = σ 2λ 2>(ρ∗ + J>λ)
(λ+ λ>ρ∗)2 , (11)
lim n,di→+∞ di/n→zi V0,d,n = σ 2 lim λ→0+
λ 2>(ρ∗ + J>λ)
(λ+ λ>ρ∗)2 , (12)
lim n,di→+∞ di/n→zi Vλ,d,n = σ 2 m∑ i=1 r∗2i , (13)
lim n,di→+∞ di/n→zi V0,d,n = σ 2 lim λ→0+ m∑ i=1 r∗2i . (14)
Figure 1 illustrates the theoretical and numerical values of the bias, variance, and total risk. We observe a triple descent in Figure 1a where the covariance matrix has three blocks, and a quadruple descent in Figure 1b where the covariance has four blocks. In the three-block example, we set λ3 λ2 λ1 (λ1 = 1, λ2 = 100, λ3 = 1000). In the four-block example, we set λ4 λ3 λ2 λ1 (λ1 = 1, λ2 = 100, λ3 = 104, λ4 = 107). For the values of other parameters, please refer to the caption of Figure 1 Our findings provide an explanation for the occurrence of sample-wise multiple descent: it occurs when the covariance matrix is highly ill-conditioned. Moreover, we find that the generalization risk curve is continuous in ridge regression (λ > 0) while it blows up at n = d in ridgeless regression (λ = 0). We can see the singularity (at n = d = 200) of the ridgeless generalization risk curve in Figure 2a.
Following Theorem 1, we focus on the variance in the ridgeless case (λ = 0) and further study the expressions in Equation (13) and Equation (14). We find that the variance exhibits sharply different behaviors in the underparameterized and overparameterized regimes. Recall that we will let n, di → +∞ and keep di/n → zi. Then d/n → ∑ i∈[m] zi. If lim d/n = ∑ i∈[m] zi > 1, we are in the underparameterized regime. In this regime, the bias vanishes and therefore the risk equals the variance. If lim d/n < 1, we are in the overparameterized regime.
Theorem 2. If d/n→ ∑ i∈[m] zi > 1 and r ∗ = (r∗1 , . . . , r ∗ m) solves
min ri≥0 ∑ i∈[m] 1 λi r2i subject to √∑ i∈[m] r2i + 1 = ∑ i∈[m] √ ziri ,
then we have an optimality condition for r∗:
r∗i r∗j = λi λj · √ ziA ∗ − r∗i√ zjA∗ − r∗j , i, j ∈ [m] , (15)
where A∗ = √∑
i∈[m] r ∗2 i + 1. Moreover, we have limn,di→+∞
di/n→zi V0,d,n = σ
2 limλ→0+ ∑m i=1 r ∗2 i .
If d/n→ ∑ i∈[m] zi < 1, then we have
lim n,di→+∞ di/n→zi
V0,d,n = σ 2
∑ i∈[m] zi
1− ∑ i∈[m] zi .
.
Corollary 1. If m = 1 and d/n→ z1 > 1, we have limn,di→+∞ di/n→zi V0,d,n = σ 2 1 z1−1 .
Proof. In the case m = 1, we have r∗1 solves minr1≥0 1 λ1 r21 subject to √ r21 + 1 = √ z1r1. The equality constraint gives r∗21 = 1 z1−1 . Then by Theorem 2, the limiting variance is σ 2r∗21 = σ 2 1 z1−1 .
In Theorem 2, we find that in the underparameterized regime, r∗ solves an equality-constrained minimization problem. In the proof of Theorem 2, we see that the equality constraint is feasible in the underparameterized regime but infeasible in the overparameterized regime. Moreover, we present an optimality condition for r∗, which will be used in Theorem 3 to study the two-block (m = 2) case. If the data distribution is isotropic (which means that the covariance matrix is a scalar matrix), Collorary 1 shows that the limiting variance is σ2 1z1−1 , which agrees with (Hastie et al., 2019, Theorem 1).
In the overparameterized regime, however, we find that the limiting variance does not depend on the spectrum {λ1, . . . , λm}of the covariance matrix and only depends on the noise intensity σ and the ratios zi = lim di/n. This agrees with (Hastie et al., 2019, Proposition 2).
In Theorem 3, we study the case m = 2 and present a concrete closed-form formula for the limiting variance in the overparameterized regime. Recall that the limiting variance in the underparameterized regime has a closed-form σ2 ∑ i∈[m] zi
1− ∑ i∈[m] zi for general m, as shown in Theorem 2.
Theorem 3. If m = 2 and d/n→ z1 + z2 > 1, we have
lim n,di→+∞ di/n→zi
V0,d,n = σ 2 q
2 + 1
q2(z1 − 1) + 2q √ z1z2 + z2 − 1 .
where
q = λ1 (z1 − 1) + λ2 (1− z2) + √ (λ1 (z1 − 1) + λ2 (1− z2)) 2 + 4λ1λ2z1z2 2λ2 √ z1z2 . (16)
We illustrate the theoretical values predicted by Theorem 3 (overparameterized regime) and Theorem 2 (underparameterized regime) in Figure 2a and compare it to the numerical values. Corollary 2 (Triple descent in the two-block case). Assume m = 2, z1 = z2, d/n→ ζ = 2z1, and λ2/λ1 = %. Define f%(ζ) = limn,di→+∞
di/n→zi V0,d,n. We have
f(ζ) , lim %→+∞ fρ(ζ) =  σ2 ζ1−ζ ζ < 1 , σ2 ( 1 ζ−1 + 2 2−ζ − 1 ) 1 < ζ < 2
σ2 2ζ−2 ζ > 2
. (17)
There exists ζ1, ζ2, ζ3, ζ4 and %0 such that for all % > %0, we have f ′%(ζ1) < 0, f ′ %(ζ2) > 0, f ′%(ζ3) < 0, and f ′ %(ζ4) < 0.
Proof. The case ζ < 1 is already given in Theorem 2. In the sequel, assume ζ > 1. Define q as in Equation (16). We have
q = ζ +
√ ζ2(%+ 1)2 − 4ζ(%− 1)2 + 4(%− 1)2 − (ζ − 2)%− 2
2ζ% .
Recall Theorem 3, we get
f%(ζ) = lim n,di→+∞ di/n→zi
V0,d,n = 2 ( q2 + 1 ) ζ (q + 1) 2 − 2 (q2 + 1) =
2σ2
ζ (q+1) 2 q2+1 − 2 .
Direct calculation yields
lim %→+∞ f%(ζ) =  σ2 ζ1−ζ ζ < 1 , σ2 ( 1 ζ−1 + 2 2−ζ − 1 ) 1 < ζ < 2
σ2 2ζ−2 ζ > 2 .
,
g(ζ) , lim %→+∞
f ′%(ζ) =  σ2 1(ζ−1)2 ζ < 1 , σ2 ζ 2−2 (ζ2−3ζ+2)2 1 < ζ < 2
σ2 −2(ζ−2)2 ζ > 2 .
,
The function g(ζ) > 0 if ζ ∈ (√ 2, 2 ) and we have g(ζ) < 0 if ζ < √
2 or ζ > 2. Pick ζ1 > 2 > ζ2 > √ 2 > ζ3 > 1 > ζ4. Then we have g(ζ1) < 0, g(ζ2) > 0, g(ζ3) < 0, and g(ζ4) > 0. There exists %0 such that for all % > %0, we have f ′%(ζ1) < 0, f ′ %(ζ2) > 0, f ′ %(ζ3) < 0, and f ′ %(ζ4) < 0.
Collorary 2 theoretically proves that there exists triple descent when m = 2 and λ2 λ1. Note that a larger ζ = lim d/n reflects a relatively smaller n. If f ′%(ζ) < 0, then f%(ζ) decreases on a neighborhood of ζ and therefore the limiting variance increases with a relatively larger n. As n becomes relatively larger, we see an increasing stage, a decreasing stage, and finally an increasing stage in order in the overparameterized regime (n < d). When we further increase n and enter the underparameterized regime, we observe a decreasing stage. We illustrate f(ζ) in Figure 2b. In Figure 2b, we observe two singularities at ζ = 1 and ζ = 2.
3.2 OPTIMAL REGULARIZATION MONOTONIZES GENERALIZATION RISK CURVE
Recall the definition of the ridge estimator in Definition 1. Since this subsection concerns samplewise monotonicity, we add a subscript n toX and y (they are defined by Equation (1) in Section 1.1) to emphasize that they consist of n data items. Therefore we write
θ̂λ,n,d , arg min θ
1 n ‖yn −Xnθ‖22 + λ ‖θ‖ 2 2 .
In this subsection, under an assumption, we show that optimal regularization (i.e., pick λ that minimizes the generalization risk of θ̂λ,n,d) results in a monotone generalization risk curve—in other words, with optimal regularization, more data always reduces the generalization risk. The assump-
tion is that ‖Πiθ∗‖2 = √ di d , i.e., the squared norm of the projection of θ
∗ onto each eigenspace of the covariance matrix is proportional to the dimension of that eigenspace. (Nakkiran et al., 2021) showed by numerical results that optimal regularization can mitigate double descent for anisotropic data distribution. We give a partial theoretical proof of their observed phenomenon.
To ease the notation, we use γi , limn/di rather than zi , lim di/n in Theorem 4 because a larger γ reflects a relatively larger n (in the limit). Theorem 4 shows that with the optimal regularization, the limiting risk is an increasing function of γ1, . . . , γm.
Theorem 4 (Optimal regularization). If ‖Πiθ∗‖2 = √ di d , (18)
then there exists a function g(γ1, . . . , γm) such that g(γ1, . . . , γm) is increasing in every γi and
lim n,di→∞ n/di→γi inf λ>0
EXn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2
Σ = g(γ1, . . . , γm) .
4 CONCLUSION
We studied the generalization risk (test error) versus the number of training samples in ridgeless regression. Under the assumption that the data distribution is Gaussian and the spectrum distribution of its covariance matrix converges to a discrete distribution, we obtained the exact formulae for the limiting bias and variance terms using the random matrix theory when the dimension and the number of training samples go to infinity in a proportional manner. Using these formulae, we proved the sample-wise multiple descent phenomenon of the generalization risk curve.Moreover, we theoretically showed that the ridge estimator with optimal regularization can result in a monotone generalization risk curve and thereby eliminate multiple descent under some assumptions.
REFERENCES Madhu S Advani, Andrew M Saxe, and Haim Sompolinsky. High-dimensional dynamics of gener-
alization error in neural networks. Neural Networks, 132:428–446, 2020.
Zhi-Dong Bai and Yong-Qua Yin. Limit of the smallest eigenvalue of a large dimensional sample covariance matrix. In Advances In Statistics, pp. 108–127. World Scientific, 2008.
Zhidong Bai and Jack W Silverstein. Spectral analysis of large dimensional random matrices, volume 20. Springer, 2010.
Peter L Bartlett, Philip M Long, Gábor Lugosi, and Alexander Tsigler. Benign overfitting in linear regression. Proceedings of the National Academy of Sciences, 117(48):30063–30070, 2020.
Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal. Reconciling modern machinelearning practice and the classical bias–variance trade-off. Proceedings of the National Academy of Sciences, 116(32):15849–15854, 2019.
Mikhail Belkin, Daniel Hsu, and Ji Xu. Two models of double descent for weak features. SIAM Journal on Mathematics of Data Science, 2(4):1167–1180, 2020.
Siegfried Bös and Manfred Opper. Dynamics of batch training in a perceptron. Journal of Physics A: Mathematical and General, 31(21):4835, 1998.
Lin Chen, Yifei Min, Mikhail Belkin, and Amin Karbasi. Multiple descent: Design your own generalization curve. arXiv preprint arXiv:2008.01036, 2020a.
Lin Chen, Yifei Min, Mingrui Zhang, and Amin Karbasi. More data can expand the generalization gap between adversarially robust and standard models. In International Conference on Machine Learning, pp. 1670–1680. PMLR, 2020b.
Donald L Cohn. Measure theory. Springer, 2013.
Edgar Dobriban and Stefan Wager. High-dimensional asymptotics of prediction: Ridge regression and classification. The Annals of Statistics, 46(1):247–279, 2018.
Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J Tibshirani. Surprises in highdimensional ridgeless least squares interpolation. arXiv preprint arXiv:1903.08560, 2019.
Shun ichi Amari, Jimmy Ba, Roger Baker Grosse, Xuechen Li, Atsushi Nitanda, Taiji Suzuki, Denny Wu, and Ji Xu. When does preconditioning help or hurt generalization? In International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id= S724o4_WB3.
Anders Krogh and John A Hertz. Generalization in a linear perceptron in the presence of noise. Journal of Physics A: Mathematical and General, 25(5):1135, 1992.
Yann Le Cun, Ido Kanter, and Sara A Solla. Eigenvalues of covariance matrices: Application to neural-network learning. Physical Review Letters, 66(18):2396, 1991.
Tengyuan Liang, Alexander Rakhlin, and Xiyu Zhai. On the multiple descent of minimum-norm interpolants and restricted lower isometry of kernels. In Conference on Learning Theory, pp. 2683–2711. PMLR, 2020.
F. Liese and K.J. Miescke. Statistical Decision Theory: Estimation, Testing, and Selection. Springer Series in Statistics. Springer New York, 2008. ISBN 9780387731940.
Marco Loog, Tom Viering, Alexander Mey, Jesse H Krijthe, and David MJ Tax. A brief prehistory of double descent. Proceedings of the National Academy of Sciences, 117(20):10625–10626, 2020.
Song Mei and Andrea Montanari. The generalization error of random features regression: Precise asymptotics and the double descent curve. Communications on Pure and Applied Mathematics, 2019.
Yifei Min, Lin Chen, and Amin Karbasi. The curious case of adversarially robust models: More data can help, double descend, or hurt generalization. In Conference on Uncertainty in Artificial Intelligence, 2021.
Preetum Nakkiran, Gal Kaplun, Yamini Bansal, Tristan Yang, Boaz Barak, and Ilya Sutskever. Deep double descent: Where bigger models and more data hurt. In International Conference on Learning Representations, 2020. URL https://openreview.net/forum?id=B1g5sA4twr.
Preetum Nakkiran, Prayaag Venkat, Sham M. Kakade, and Tengyu Ma. Optimal regularization can mitigate double descent. In International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id=7R7fAoUygoa.
M Opper, W Kinzel, J Kleinz, and R Nehl. On the ability of the optimal perceptron to generalise. Journal of Physics A: Mathematical and General, 23(11):L581, 1990.
Dominic Richards, Jaouad Mourtada, and Lorenzo Rosasco. Asymptotics of ridge (less) regression under general source condition. In International Conference on Artificial Intelligence and Statistics, pp. 3889–3897. PMLR, 2021.
Walter Rudin. Principles of mathematical analysis. New York, NY: McGraw-Hill, Inc., 3 edition, 1976.
Christos Thrampoulidis, Samet Oymak, and Babak Hassibi. Regularized linear regression: A precise analysis of the estimation error. Proceedings of Machine Learning Research, 40:1683–1709, 2015.
Alexander Tsigler and Peter L Bartlett. Benign overfitting in ridge regression. arXiv preprint arXiv:2009.14286, 2020.
F Vallet, J-G Cailton, and Ph Refregier. Linear and nonlinear extension of the pseudo-inverse solution for learning boolean functions. EPL (Europhysics Letters), 9(4):315, 1989.
Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv preprint arXiv:1011.3027, 2010.
Roman Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018.
Timothy LH Watkin, Albrecht Rau, and Michael Biehl. The statistical mechanics of learning a rule. Reviews of Modern Physics, 65(2):499, 1993.
Denny Wu and Ji Xu. On the optimal weighted ` 2 regularization in overparameterized linear regression. arXiv preprint arXiv:2006.05800, 2020.
A EIGENDECOMPOSITION AND MORE NOTATION
Write Σ = PΛP>, where P is an orthogonal matrix and Λ = diag (λ1Id1 , . . . , λmIdm) ∈ Rd×d is a diagonal matrix. Write λ− = mini∈[m] λi and λ+ = maxi∈[m] λi. We can generate x1, . . . , xn from standard normal random vector zi ∼ N (0, Id) by setting xi = PΛ1/2zi. Therefore, if Z = ( z1 . . . zn ) ∈ Rd×n,we get
X> = ( x1 . . . xn ) = PΛ 1/2 ( z1 . . . zn ) = PΛ 1/2Z .
Take the transpose gives X = Z>Λ1/2P>. Note that every entry of Z ∈ Rd×n follows i.i.d. N (0, 1). Write Z in a row-partitioned form
Z =  Z > 1 ... Z>m  , where Zi ∈ Rn×di . Write P in a column-partitioned form
P = ( P1 . . . Pm ) ,
where Pi ∈ Rd×di . Recall that Πi ∈ Rd×d denotes the orthogonal projection to the eigenspace of λi. We have Πi = PiP>i . Define θ ′ , P>θ∗ and write it in a row-partitioned form
θ′ =  P > 1 θ ∗
... P>mθ ∗
 =  θ
′ 1 ... θ′m  , (19) where θ′i ∈ Rdi . Then ‖θ′i‖2 =
∥∥P>i θ∗∥∥2 = ∥∥PiP>i θ∗∥∥2 = ‖Πiθ∗‖2. We summarize part of the notation above in Table 1.
B BIAS AND VARIANCE UNDER EIGENDECOMPOSITION
Lemma 1 characterizes the smallest and largest eigenvalue of Z >Z d (if n/d → γ < 1) and ZZ>
n (if n/d→ γ > 1). Recall that we study the asymptotic regime di/n→ zi. Therefore γ = 1∑
j∈[m] zj .
Lemma 1 ((Bai & Yin, 2008, Theorem 2)). Let Z ∈ Rd×n be a random matrix whose entries are i.i.d. N (0, 1) random variables. As n, d→ +∞, n/d→ γ ∈ (0, 1), we have
limλmin
( Z>Z
d
) = (1−√γ)2 , limλmax ( Z>Z
d
) = (1 + √ γ) 2
almost surely. If γ ∈ (1,∞), as n, d→ +∞, n/d→ γ, we have
limλmin
( ZZ>
n
) = ( 1− √ 1/γ )2
, limλmax
( ZZ>
n
) = ( 1 + √ 1/γ )2
almost surely.
Lemma 2 (Corollary 5.35 (Vershynin, 2010)). Let A be an N × n matrix whose entries are independent standard normal random variables. Then for every t ≥ 0, with probability at least 1− 2 exp ( −t2/2 ) one has √ N − √ n− t ≤ smin(A) ≤ smax(A) ≤ √ N + √ n+ t ,
where smin(A) and smax(A) are the smallest and largest singular value of A. Lemma 3. Let Z ∈ Rd×n be a random matrix whose entries are i.i.d. N (0, 1) random variables, where d = d(n) satisfies limn→∞ nd(n) = γ. There exists universal positive constants C1, C2, N such that for all n > N , we have
0 < C1 < 1
n s2min(Z) ≤
1 n s2max(Z) < C2 .
Proof. Since d(n) n, with loss of generality, we assume n/d → γ ∈ (0, 1). Take t = c1 √ n in Lemma 2, where c1 = 12 ( 1√ γ − 1 ) > 0. With probability at least 1− 2e−c21n/2, we have
√ d− √ n− c1 √ n ≤ smin(Z) ≤ smax(Z) ≤ √ d+ √ n+ c1 √ n .
Therefore, we deduce(√ d
n − 1− c1 )2 ≤ 1 n s2min(Z) ≤ 1 n s2max(Z) ≤ (√ d n + 1 + c1 )2 .
Define C1 = 18 ( 1√ γ − 1 )2 > 0 and C2 = ( 3√ γ + 1 )2 . Then there exists a universal constant N1 such that for all n > N1, with probability at least 1− 2e−c 2 1n/2, we have
0 < C1 < 1
n s2min(Z) ≤
1 n s2max(Z) < C2 .
Define event En = { C1 < 1 ns 2 min(Z) ≤ 1ns 2 max(Z) < C2 }c . Then we have Pr {En} ≤ 2e−c 2 1n/2.
Since ∑ n≥1 Pr {En} ≤ ∑ n≥1 2e
−c21n/2 < ∞, then the probability that infinitely many of En occur is 0, i.e.,
Pr { lim sup
n En
} = 0 .
Therefore, there exists a universal constant N2 such that for all n > N2, En does not happen, in other words,
0 < C1 < 1
n s2min(Z) ≤
1 n s2max(Z) < C2
holds.
Lemma 4. Let Z ∈ Rd×n be a random matrix whose entries are i.i.d. N (0, 1) random variables, and let p be a fixed positive integer which is viewed as a constant and hidden in .. If n d, we have E tr ( ZZ> ) n2, E tr ( ZZ> )2 n3, and E ‖Z‖p2 . np/2. Proof. We have
E tr ( ZZ> ) = E ‖Z‖2F = ∑ i∈[d],j∈[n] Ez2ij = nd n2 .
Write Z =  z > 1 ... z>d , where zi ∈ Rn and zi ∼ N (0, In). We have E (z>i zi)2 = E ‖zi‖42 = n(n + 2). For i 6= j, we deduce E ( z>i zj )2 = E ( ‖zi‖2 ‖zj‖2 u >v )2 where u, v ∼ Unif ( Sn−1
) and ‖zi‖2 , ‖zj‖2 , u, v are independent. Then we get
E ( ‖zi‖ ‖zj‖ s>i sj )2 = E ‖zi‖22 ‖zj‖ 2 2 ( s>i sj )2 = n2Eu21 = n2 · 1
n = n .
As a result, we have E tr ( ZZ> )2 = E ∥∥ZZ>∥∥2 F = ∑ i,j∈[d] E ( z>i zj )2 = dn(n+ 2) + ( d2 − d ) n n3 .
By (Vershynin, 2018), there exists a universal constant C > 0 such that for any t > 0, P { ‖Z‖2 > C (√ n+ √ d+ t )} < 2e−t 2 . Define K = C (√ n+ √ d ) . Then we have
P {‖Z‖2 > K + t} < 2e −t2/C2 . (20)
Recall Γ(z) = ∫∞
0 xz−1e−xdx. Setting t = C √ u in the equation below yields∫ ∞
0
e−t 2/C2tp−1dt . ∫ ∞ 0 e−uu p 2−1du = Γ (p 2 ) 1 .
Then we can bound the following integral∫ ∞ K P {‖Z‖2 ≥ t} pt p−1dt
= ∫ ∞ 0 P {‖Z‖2 ≥ K + t} p (t+K) p−1 dt
. ∫ ∞
0
e−t 2/C2 (t+K) p−1 dt
. ∫ ∞
0
e−t 2/C2 ( tp−1 +Kp−1 ) dt
= ∫ ∞ 0 e−t 2/C2tp−1dt+Kp−1 ∫ ∞ 0 e−t 2/C2dt .n p−1 2 ,
where the first inequality is because of Equation (20). We are in a position to bound E ‖Z‖p2: E ‖Z‖p2
= ∫ ∞ 0 P {‖Z‖p2 ≥ u} du
= ∫ ∞ 0 P {‖Z‖2 ≥ t} pt p−1dt
= ∫ K 0 P {‖Z‖2 ≥ t} pt p−1dt+ ∫ ∞ K P {‖Z‖2 ≥ t} pt p−1dt .np/2 + n(p−1)/2
.np/2 ,
where the first inequality is because∫ K 0 P {‖Z‖2 ≥ t} pt p−1dt ≤ ∫ K 0 ptp−1dt = Kp . np/2 .
Lemma 5. The following equation for the bias term Bλ,d,n (defined in Equation (3)) holds Bλ,d,n =E [ ‖Λ1/2 ( Id − Λ1/2Z ( nλIn + Z >ΛZ )−1 Z>Λ1/2 ) θ′‖22 ] (21)
=E [ ‖Λ1/2 ( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 θ′‖22 ] . (22)
Moreover, we have |Bλ,d,n| . ‖θ∗‖22 and limλ→0+ Bλ,d,n = B0,d,n. For all sufficiently large n and d such that n/d → γ ∈ (0, 1), we have 0 ≤ ddλBλ,d,n . ‖θ
∗‖22. Therefore, {Bλ,d,n} is uniformly bounded and uniformly equicontinuous with respect to λ ∈ (0,∞).
Proof. Introduce the shorthand notation M = Λ1/2ZZ>Λ1/2 ∈ Rd×d, A = Id + 1nλM ∈ R d×d, N = nλIn + Z >ΛZ ∈ Rn×n, and Q = Id − Λ1/2ZN−1Z>Λ1/2 ∈ Rd×d. Because
X>(nλIn +XX >)−1X = PΛ1/2Z ( nλIn + Z >ΛZ )−1 Z>Λ1/2P> ,
we have
Bλ,d,n =E [ ‖(Id − PΛ1/2Z ( nλIn + Z >ΛZ )−1 Z>Λ1/2P>)θ∗‖2PΛP> ]
=E [ ‖Λ1/2 ( Id − Λ1/2Z ( nλIn + Z >ΛZ )−1 Z>Λ1/2 ) θ′‖22 ] =E [ ‖Λ1/2Qθ′‖22 ] .
Using the Sherman–Morrison-Woodbury formula yields
N−1 = 1
nλ In −
1
(nλ) 2Z >Λ1/2
( I + 1
nλ Λ1/2ZZ>Λ1/2
)−1 Λ1/2Z
= 1
nλ
( In − Z>Λ1/2 ( nλId + Λ 1/2ZZ>Λ1/2 )−1 Λ1/2Z ) = 1
nλ
( In − Z>Λ1/2 (nλId +M)−1 Λ1/2Z ) . (23)
It follows that
Q =Id − Λ1/2ZN−1Z>Λ1/2
=Id − 1
nλ Λ1/2Z
( In − Z>Λ1/2 (nλId +M)−1 Λ1/2Z ) Z>Λ1/2
=Id − M
nλ
( Id − (nλId +M)−1M ) =Id − M
nλ
( Id − (nλId +M)−1 (nλId +M − nλId) ) =Id −M (nλId +M)−1
= ( Id + 1
nλ M )−1 =A−1 .
Therefore, we deduce
Bλ,d,n = E [ ‖Λ1/2 ( Id + 1
nλ M
)−1 θ′‖22 ] = E [ ‖Λ1/2A−1θ′‖22 ] .
Because ∥∥Λ1/2∥∥
2 . 1 and ∥∥∥(Id + 1nλΛ1/2ZZ>Λ1/2)−1∥∥∥ 2 ≤ 1, we have
‖Λ1/2A−1θ′‖22 . ‖θ′‖ 2 2 = ‖θ ∗‖22 .
Therefore |Bλ,d,n| . ‖θ∗‖22. Moreover, by the dominated convergence theorem,
lim λ→0+ Bλ,d,n = B0,d,n .
We compute the derivative of A−1:
dA−1
dλ = −A−1 dA dλ A−1 =
MA−2
nλ2 .
The matrix Mn = Λ1/2ZZ>Λ1/2 n ∈ R d×d is positive semidefinite and its d − n smallest eigenvalues are zeros. Its non-zero eigenvalues are the same as the non-zero eigenvalues of Z >ΛZ n . Because all eigenvalues of Z >ΛZ n are positive almost surely, the spectrum of M n consists of d− n zeros and the
spectrum of Z >ΛZ n . We study the range of the spectrum of Z>ΛZ n . Because λ− Z>Z n 4 Z>ΛZ n 4 λ+ Z>Z n , we deduce
λmin
( Z>ΛZ
n
) ≥ λ−λmin ( Z>Z
n
) → λ− ( 1− √ 1/γ )2
(24)
λmax
( Z>ΛZ
n
) ≤ λ+λmax ( Z>Z
n
) → λ+ ( 1 + √ 1/γ )2 . (25)
Define L1 = λ−λmin ( Z>Z n ) and L2 = λ+λmax ( Z>Z n ) . We get limn,d→+∞
n/d→γ<1 L1 =
λ−
( 1− √ 1/γ )2
, limn,d→+∞ n/d→γ<1 L2 = λ+
( 1 + √ 1/γ )2 and
spec
( Z>ΛZ
n
) ⊆ [L1, L2] .
We bound ∥∥MA−3∥∥
2 ∥∥MA−3∥∥ 2 =n ∥∥∥∥∥Mn ( Id + M nλ )−3∥∥∥∥∥ 2
=n max s∈spec(Mn )
s
(1 + s/λ) 3
=n max s∈{0}∪spec ( Z>ΛZ n . ) s (1 + s/λ) 3
=n max s∈spec ( Z>ΛZ n . ) s (1 + s/λ) 3
≤n max s∈[L1,L2]
s
(1 + s/λ) 3 .
We compute ddλ‖Λ 1/2A−1θ′‖22:
d
dλ ‖Λ1/2A−1θ′‖22
= 1 nλ2 θ′> ( A−1ΛMA−2 +MA−2ΛA−1 ) θ′ = 1
nλ2 ( A−1θ′ )> ( ΛMA−1 +MA−1Λ ) ( A−1θ′ ) Next, we bound ∣∣ d dλ‖Λ
1/2A−1θ′‖22 ∣∣:∣∣∣∣ ddλBλ,d,n
∣∣∣∣ ≤ 1 nλ2 ∥∥MA−2ΛA−1 +A−1ΛMA−2∥∥ 2 ‖θ′‖22 ≤ 2 nλ2 ∥∥MA−2ΛA−1∥∥ 2 ‖θ′‖22 = 2
nλ2 ∥∥MA−3AΛA−1∥∥ 2 ‖θ′‖22
≤ 2 nλ2 ∥∥MA−3∥∥ 2 ∥∥AΛA−1∥∥ 2 ‖θ′‖22 . 1
λ2 max s∈[L1,L2]
s
(1 + s/λ) 3 ‖θ ′‖22 ,
where the last inequality is because ∥∥AΛA−1∥∥
2 = ‖Λ‖2 ≤ λ+ . 1. Define f(s) = s (1+s/λ)3 .
Because f ′(s) = λ 3(λ−2s) (λ+s)4 , the function f is increasing on [0, λ/2] and decreasing on [λ/2,+∞).
If λ ≤ 2L1, we have max
s∈[L1,L2]
s
(1 + s/λ) 3 =
L1
(1 + L1/λ) 3 .
It follows that
1 λ2 · L1
(1 + L1/λ) 3 =
L1λ
(λ+ L1) 3 ≤ max λ∈[0,2L1]
L1λ
(λ+ L1) 3 .
1
L1 .
If λ ≥ 2L2, we get
1
λ2 max s∈[L1,L2]
s
(1 + s/λ) 3 =
1 λ2 · L2 (1 + L2/λ) 3 ≤ max λ∈[2L2,∞)
L2λ
(λ+ L2) 3 .
1 L2 ≤ 1 L1 .
If 2L1 < λ < 2L2, we obtain
1
λ2 max s∈[L1,L2]
s
(1 + s/λ) 3 .
1 λ . 1 L1 .
In all three cases, we show that 1λ2 maxs∈[L1,L2] s (1+s/λ)3 . 1L1 . It follows that∣∣∣∣ ddλ‖Λ1/2A−1θ′‖22 ∣∣∣∣ . 1L1 ‖θ′‖22 = ‖θ ′‖22
λ−λmin ( Z>Z n ) ‖θ′‖22 λmin ( Z>Z n ) . By Lemma 3, there exists a universal constant n0 such that for all n > n0, one has 1
λmin ( Z>Z n ) . 1. Thus we conclude that ∣∣∣∣ ddλ‖Λ1/2A−1θ′‖22
∣∣∣∣ . ‖θ′‖22 . We can exchange differentiation and expectation and get
d
dλ Bλ,d,n = E
[ d
dλ ‖Λ1/2A−1θ′‖22 ] and ∣∣∣∣ ddλBλ,d,n ∣∣∣∣ = E [∣∣∣∣ ddλ‖Λ1/2A−1θ′‖22 ∣∣∣∣] . ‖θ′‖22 .
Lemma 6. The following equation for the variance term holds
Vλ,d,n =σ 2E‖Λ1/2 ( λnId + Λ 1/2ZZ>Λ1/2 )−1 Λ1/2Z‖22
=σ2E‖ΛZ ( λnIn + Z >ΛZ )−1 ‖22 .
Moreover, for all sufficiently large n and d such that n/d → γ 6= 1, we have limλ→0+ Vλ,d,n = V0,d,n, |Vλ,d,n| . 1 and ∣∣ d dλVλ,d,n ∣∣ . 1. Therefore, {Vλ,d,n} is uniformly bounded and uniformly equicontinuous with respect to λ ∈ (0,∞).
Proof. As in the proof of Lemma 5, defineM = Λ1/2ZZ>Λ1/2 ∈ Rd×d andN = nλIn+Z>ΛZ ∈ Rn×n. Recalling Σ = PΛP> and X = Z>Λ1/2P>, we have
Vλ,d,n =σ 2E tr [ XΣX>(nλIn +XX >)−2 ]
=σ2E tr [ Z>Λ2ZN−2 ] =σ2E tr [ N−1Z>Λ2ZN−1
] =σ2E ∥∥ΛZN−1∥∥2 F .
Recalling Equation (23) yields
ΛZN−1 = 1
nλ ΛZ
( In − Z>Λ1/2 (nλId +M)−1 Λ1/2Z ) = 1
nλ Λ1/2
( In −M (nλId +M)−1 ) Λ1/2Z
=Λ1/2 (nλId +M) −1 Λ1/2Z .
Define R = Λ1/2Z ∈ Rd×n. We get
(nλId +M) −1 Λ1/2Z = ( nλId +RR >)−1R . Notice that if 0 < a < b, then aId +RR> 4 bId +RR>. We deduce(
bId +RR >)2 − (aId +RR>)2 = (b2 − a2) Id + 2 (b− a)RR> < 0 .
Thus ( bId +RR >)2 < (aId +RR>)2, which implies (bId +RR>)−2 4 (aId +RR>)−2. We get
R> ( bId +RR >)−2R 4 R> (aId +RR>)−2R , tr ( R> ( bId +RR
>)−2R) ≤ tr(R> (aId +RR>)−2R) Let λ0 (·) denote the smallest non-zero eigenvalue of a positive semidefinite matrix. We bound the Frobenius norm ∥∥∥(nλId +M)−1 Λ1/2Z∥∥∥2
F = tr ( R> ( nλId +RR >)−1R) ≤ tr ( lim λ→0+ R> ( nλId +RR
>)−2R) = tr ( R>R
)+ = tr ( Z>ΛZ
)+ . tr ( Z>Z )+ .
It follows that∥∥ΛZN−1∥∥2 F = ∥∥∥Λ1/2 (nλId +M)−1 Λ1/2Z∥∥∥2 F . ∥∥∥(nλId +M)−1 Λ1/2Z∥∥∥2 F . tr ( Z>Z )+ = tr ( ZZ> )+ . If n/d → γ < 1, the matrix Z>Z is full-rank almost surely. Then, using the formula for the mean of inverse Wishart distribution, we have E tr ( Z>Z )+ = trE ( Z>Z )−1 = tr ( In
d−n−1
) 1. If
n/d → γ > 1, the matrix ZZ> is full-rank almost surely. Similarly, we have E tr ( ZZ> )+ =
trE ( ZZ> )−1 1. By the dominated convergence theorem, we have limλ→0+ Vλ,d,n = V0,d,n. Moreover, Vλ,d,n . E ∥∥ΛZN−1∥∥2 F . 1.
Next we bound ddλV (θ̂). Because dN−1 dλ = −N −1 dN dλN −1 = −nN−2, we deduce
d
dλ
∥∥ΛZN−1∥∥2 2 = d dλ tr ( N−1Z>Λ2ZN−1 ) = −2n tr ( Z>Λ2ZN−3 ) ≤ 0 .
On the other hand, we have tr ( Z>Λ2ZN−3 ) = tr ( N−3/2Z>Λ2ZN−3/2
) . tr ( N−3/2Z>ΛZN−3/2
) = tr ( Z>ΛZN−3
) =
∑ s∈spec(Z>ΛZ)
s
(λn+ s) 3 .
Because the number of non-zero eigenvalues of Z>ΛZ equals rank ( Z>ΛZ ) = n ∧ d n, we get∣∣∣∣ ddλ ∥∥ΛZN−1∥∥22 ∣∣∣∣ n tr (Z>Λ2ZN−3) . n2 max s∈spec(Z>ΛZ) s (λn+ s) 3 = max s∈spec ( Z>ΛZ n ) \{0} s (λ+ s) 3 .
If γ < 1, the matrix Z >ΛZ n is full-rank almost surely. By Equation (24) and Equation (25) in the proof of Lemma 5, there exists universal positive constants C1 and C2 such that spec ( Z>ΛZ n ) ⊆ [C1, C2] for all sufficiently large n and d such that n/d→ γ < 1. If γ > 1, the non-zero eigenvalues of Z
>ΛZ n and M n are the same. The matrix M n is full-rank almost surely. Thus spec ( Z>ΛZ n ) \{0} =
spec ( M n ) . Because y>Λ−1y . y>y,
λmin
( M
n
) = min
x6=0
x> Λ 1/2ZZ>Λ1/2
n x
x>x = min y 6=0
y> ZZ >
n y
y>Λ−1y & min y 6=0
y> ZZ >
n y
y>y = λmin
( ZZ>
n
) .
Similarly, we get
λmax
( M
n
) . λmax ( ZZ>
n
) .
Therefore, there exists universal positive constants C1 and C2 such that
spec
( Z>ΛZ
n
) \ {0} ⊆ [ C1λmin ( Z>Z
n
) , C2λmax ( Z>Z
n
)] .
Thus in both cases, we have shown that there exists universal positive constants C1 and C2 such that
spec
( Z>ΛZ
n
) \ {0} ⊆ [ C1λmin ( Z>Z
n
) , C2λmax ( Z>Z
n
)] .
Define L1 = C1λmin ( Z>Z n ) and L2 = C2λmax ( Z>Z n ) . As a result, we get∣∣∣∣ ddλ ∥∥ΛZN−1∥∥22
∣∣∣∣ . maxs∈[L1,L2] s(λ+ s)3 . Define f(s) = s
(λ+s)3 . Because f ′(s) = λ−2s(λ+s)4 , the function f is increasing on [0, λ/2] and
decreasing on [λ/2,+∞). If λ ≥ 2L2 or λ ≤ 2L1, we get
max s∈[L1,L2]
s
(λ+ s) 3 ≤
L1
(λ+ L1) 3 ∨
L2
(λ+ L2) 3 ≤
1
L21 .
If 2L1 < λ < 2L2, we get
max s∈[L1,L2]
s
(λ+ s) 3 .
1
λ2 .
1
L21 .
As a result, for all sufficiently large n, we have∣∣∣∣ ddλ ∥∥ΛZN−1∥∥22 ∣∣∣∣ = maxs∈[L1,L2] s(λ+ s)3 . 1L21 . 1λ2min (Z>Zn ) . 1 ,
where the final inequality is because of Lemma 3. We can exchange the expectation and differentiation and obtain
d
dλ Vλ,d,n = σ
2E d
dλ
∥∥ΛZN−1∥∥2 2
and ∣∣∣∣ ddλVλ,d,n ∣∣∣∣ ≤ σ2E ∣∣∣∣ ddλ ∥∥ΛZN−1∥∥22 ∣∣∣∣ . 1 .
C LEMMAS ON STIELTJES TRANSFORM
Definition 3 (Stieltjes transform). The Stieltjes transform of a distribution with cumulative distribution function F is defined by
sF (z) =
∫ 1
λ− z dF (λ) (z ∈ H , {z ∈ C | =z > 0}) .
Lemma 7 (Theorem 4.3 (Bai & Silverstein, 2010)). Suppose that the entries of Xn ∈ Cn×p are complex random variables that are independent for each n and identically distributed for all n and satisfy E [ |x11 − Ex11|2 ] = 1. Also, assume that Tn = diag(τ1, . . . , τp), τi is real, and the
empirical distribution function of {τ1, . . . , τp} converges almost surely to a probability distribution function H as n → ∞. The entries of both Xn and Tn may depend on n, which is suppressed for brevity. Set Bn = An + 1nXnTnX ∗ n, where X ∗ n is the conjugate transpose of Xn, An is Hermitian, n× n satisfying FAn → FA almost surely, where FA is a distribution function (possibly defective) on the real line. Assume also that Xn, Tn, and An are independent. When p = p(n) with p/n → y > 0 as n → ∞, then, almost surely, FBn , the empirical spectral distribution of the eigenvalues of Bn, converges vaguely, as n → ∞, to a (nonrandom) distribution function F , where for any z ∈ C+ = {z ∈ C | =z > 0}, its Stieltjes transform s = s(z) is the unique solution in C+ to the equation
s = sA
( z − y ∫ τdH(τ)
1 + τs
) ,
where sA is the Stieltjes transform of FA.
Lemma 8. If the functions fα, gα : I → R satisfy fα(x)− gα(x)→ 0 uniformly as α→ +∞, then limα→+∞ (infx∈I f(x)− infx∈I g(x)) = 0.
Proof. Because fα(x)− gα(x)→ 0 uniformly as α→ +∞, we have for ∀ > 0, there exists N( ) such that for ∀α > N( ) and ∀x ∈ I , it holds that |fα(x)− gα(x)| < . Therefore, we get
gα(x)− < fα(x) < gα(x) + .
Thus we obtain
inf x∈I
fα(x) ≤ fα(x) < gα(x) +
inf x∈I
gα(x)− ≤ gα(x)− < fα(x) ,
which in turn implies
inf x∈I fα(x) ≤ inf x∈I gα(x) +
inf x∈I gα(x)− ≤ inf x∈I fα(x) .
It follows that |infx∈I fα(x)− infx∈I gα(x)| ≤ . In other words, we proved
lim α→+∞ ( inf x∈I f(x)− inf x∈I g(x) ) = 0 .
Lemma 9. Define N = λnIn + Z>ΛZ. Then we have
lim n,di→+∞ di/n→zi
tr ( N−1 ) = d
dλ inf ρ∈Rm+
log λ+ ∑
i∈[m]
λiρi + ∑ i∈[m] ( ρi − zi ( log ρi zi + 1 )) , (26)
lim n,di→+∞ di/n→zi
1 n log det N n = inf ρ∈Rm+
log λ+ ∑
i∈[m]
λiρi + ∑ i∈[m] ( ρi − zi ( log ρi zi + 1 )) . (27)
Proof. Proof of Equation (26). We apply Lemma 7 with An = 0n×n, Xn = Z> ∈ Rn×d, Tn = Λ, and Bn = 1nZ
>ΛZ. The distribution function of 0n×n converges to 1t≤0 and its Stieltjes transform is sA(z) =
∫ 1
λ−zd1λ≤0 = − 1 z . The empirical distribution function
of {λ1, . . . , λ1︸ ︷︷ ︸ d1 , . . . , λm . . . , λm︸ ︷︷ ︸ dm
} is Hn,di(t) = ∑ i∈[m] di d 1t≤λi . Recall di/n → zi. Thus
di/d → zi/K, where d/n → y = ∑ j∈[m] zj . The empirical distribution function converges to
H(t) = ∑ i∈[m] zi y 1t≤λi . Then the empirical spectral distribution of the eigenvalues of 1 nZ >ΛZ converges vaguely to a nonrandom distribution function F and its Stieltjes transform is
s = s(z) = lim n,di→+∞ di/n=zi
1 n tr
( 1
n Z>ΛZ − zIn )−1 = lim n,di→+∞ di/n=zi tr ( Z>ΛZ − znIn )−1 (this is because of (Bai & Silverstein, 2010, Theorem B.9)). By Lemma 7, s(z) is the unique solution in C+ to the equation
s(z) = sA
( z − y ∫ τdH(τ)
1 + τs
) = − 1
z − ∑ i∈[m] λizi 1+λis(z) ,
which gives
s(z) z − ∑ i∈[m] λizi 1 + λis(z)  = −1 . We want to prove Equation (26) first. The lefthand side of Equation (26) equals
lim n,di→+∞ di/n=zi
tr ( λnIn + Z >ΛZ )−1 = s(−λ) .
Because the matrix 1nZ >ΛZ is positive semidefinite and thereby all of its eigenvalues are nonnegative, its limiting spectral distribution is supported on [0,∞). The Stieltjes transform s(z) of the limiting spectral distribution can be continuously extended to (−∞, 0). Therefore, for ∀λ > 0, s(−λ) is the unique solution to the following equation
s(−λ) ( λ+
m∑ i=1 λizi 1 + λis(−λ)
) = 1 . (28)
We will verify that
d
dλ inf ρ∈Rm+
log λ+ m∑
j=1
λjρj + m∑ j=1 ( ρj − zj(log ρj zj + 1) ) satisfies Equation (28). Take a minimizer ρ∗ of Equation (5). Using the envelope theorem yields
d
dλ inf ρ∈Rm+
log λ+ m∑
j=1
λjρj + m∑ j=1 ( ρj − zj(log ρj zj + 1) ) = 1 λ+ ∑m j=1 λjρ ∗ j . (29)
Plugging the righthand side of Equation (29) into Equation (28), we get
1 λ+ ∑m j=1 λjρ ∗ j
( λ+
m∑ i=1
λizi 1 + λi · 1λ+∑mj=1 λjρ∗j ) = 1 .
Rewriting the above equation yields m∑ i=1
λizi
1 + λi · 1λ+∑mj=1 λjρ∗j = m∑ i=1 λiρ ∗ i .
It suffices to show that each summand on the lefthand side equals its counterpart on the righthand side
λizi
1 + λi · 1λ+∑mj=1 λjρ∗j = λiρ ∗ i .
We need to show zi ρ∗i = 1 + λi · 1 λ+ ∑m j=1 λjρ ∗ j ,
which is equivalent to Equation (6) and therefore holds. Hence we have proved Equation (26).
Proof of Equation (27). We use α to denote the indices n, di. Define
h(λ) = inf ρ∈Rm+
log λ+ m∑
j=1
λjρj + m∑ j=1 ( ρj − zj(log ρj zj + 1) ) . First, we want to show that limλ0→+∞ (h(λ0)− log λ0) = 0. Define
lλ0(ρ) = log 1 + 1 λ0 m∑ j=1 λjρj + m∑ j=1 ( ρj − zj ( log ρj zj + 1 )) ,
q(ρ) = m∑ j=1 ( ρj − zj ( log ρj zj + 1 )) .
The Hessian matrix of q(p) is diag ( z1 ρ21 , . . . , zmρ2m ) , which is positive definite since zi, ρi > 0. Therefore, q(p) is convex and the minimum of q(ρ) on Rm+ is attained at ρ = z, where z = (z1, . . . , zm) >. The minimum is infρ∈Rm+ q(ρ) = q(z) = 0. Because lim‖ρ‖2→+∞ lλ0(ρ) = +∞, there exists a universal constant K1 > ‖z‖2 > 0 such that lλ0(ρ) > lλ0(z) for all ‖ρ‖2 > K1. Define E = { ρ ∈ Rm+ | ‖ρ‖2 ≤ K1 } . We have z ∈ E , infρ∈E lλ0(ρ) = infρ∈Rm+ lλ0(ρ), and infρ∈E q(ρ) = infρ∈Rm+ q(ρ) = 0. Therefore, we get
h(λ0)− log λ0 = inf ρ∈Rm+ lλ0(ρ) = inf ρ∈E lλ0(ρ)− inf ρ∈E q(ρ) . (30)
On E, there exists a universal constant K2 > 0 such that ∑ j∈[m] λjρj < K2. Thus on E, we deduce
0 < lλ0(ρ)− q(ρ) = log 1 + 1 λ0 m∑ j=1 λjρj  < log(1 + K2 λ0 ) .
The right-hand side log (
1 + K2λ0 ) → 0 as λ0 → +∞. Thus limλ0→+∞ (lλ0(ρ)− q(ρ)) = 0
uniformly for ρ ∈ E. By Lemma 8, we get
lim λ0→+∞ ( inf ρ∈E lλ0(ρ)− inf ρ∈E q(ρ) ) = 0 .
Recalling Equation (30) yields
lim λ0→+∞
(h(λ0)− log λ0) = 0 . (31)
Define fα(λ) = 1n log det N n . Second, we want to show limα fα(λ) = h(λ), where limα means
limn,di→+∞ di/n=zi
. We have fα(λ)− fα(λ0) = ∫ λ λ0 f ′α(x)dx for ∀λ, λ0 > 0. It follows that
|fα(λ)− h(λ)| ≤ |fα(λ)− h(λ) + h(λ0)− fα(λ0) + fα(λ0)− log λ0 + log λ0 − h(λ0)| ≤ |fα(λ)− h(λ) + h(λ0)− fα(λ0)|+ |fα(λ0)− log λ0|+ |log λ0 − h(λ0)|
= ∣∣∣∣∣ ∫ λ λ0 f ′α(x)dx− (h(λ)− h(λ0)) ∣∣∣∣∣+ |fα(λ0)− log λ0|+ |log λ0 − h(λ0)| . Taking lim supα on both sides gives
lim sup α |fα(λ)− h(λ)| ≤ lim sup α ∣∣∣∣∣ ∫ λ λ0 f ′α(x)dx− (h(λ)− h(λ0)) ∣∣∣∣∣+lim supα |fα(λ0)− log λ0|+|log λ0 − h(λ0)| . (32)
Recall f ′α(λ) = trN −1 and limα f ′α(λ) = h ′(λ) (this is exactly Equation (26)). Because∣∣trN−1∣∣ = trN−1 ≤ 1λ and ∫ λλ0 1xdx < +∞, by the dominated convergence theorem, we have lim α ∫ λ λ0 f ′α(x)dx = ∫ λ λ0 h′(x)dx = h(λ)− h(λ0) .
It follows that
lim sup α ∣∣∣∣∣ ∫ λ λ0 f ′α(x)dx− (h(λ)− h(λ0)) ∣∣∣∣∣ = limα ∣∣∣∣∣ ∫ λ λ0 f ′α(x)dx− (h(λ)− h(λ0)) ∣∣∣∣∣ = 0 . (33) Since
fα(λ0)−log λ0 = 1
n log det
( λ0In + 1
n Z>ΛZ ) − 1 n log det (λ0In) = 1 n log det ( In + 1 nλ0 Z>ΛZ ) and the matrix 1nλ0Z >ΛZ is positive semidefinite, we have
fα(λ0)− log λ0 ≥ 0 .
We have
fα(λ0)− log λ0
= 1
n log det
( In + 1
nλ0 Z>ΛZ ) ≤ 1 n log det ( In + λ+ nλ0 Z>Z
) ≤ log ( 1 +
λ+ λ0 λmax
( Z>Z
n )) ≤λ+ λ0 λmax ( Z>Z n ) .
Then taking lim supα, we get
lim sup α |fα(λ0)− log λ0| = lim sup α (fα(λ0)− log λ0) ≤ λ+ λ0 lim sup α λmax
( Z>Z
n
) . 1
λ0 ,
(34) where the last inequality is because lim supα λmax ( Z>Z n ) = ( 1 + √ γ ∨ 1γ )2 1 by Lemma 1. Using Equation (32), Equation (33) and Equation (34) gives
lim sup α |fα(λ)− h(λ)| .
1
λ0 + |log λ0 − h(λ0)| .
Then taking limλ0→+∞ and recalling Equation (31) yields
lim α |fα(λ)− h(λ)| = lim sup α |fα(λ)− h(λ)| = 0 .
Therefore, we conclude limα fα(λ) = h(λ).
Lemma 10. Define N = λnIn+Z>ΛZ = λnIn+ ∑ i∈[m] λiZiZ > i . The following equation holds
lim n,di→+∞ di/n→zi
E [ ∂
∂λi
1 n log det N n
] = ∂
∂λi inf ρ∈Rm+
log λ+ ∑
i∈[m]
λiρi + ∑ i∈[m] ( ρi − zi ( log ρi zi + 1 )) , (35)
lim n,di→+∞ di/n→zi
E [ ∂2
∂λj∂λi
1 n log det N n
] = ∂2
∂λj∂λi inf ρ∈Rm+
log λ+ ∑
i∈[m]
λiρi + ∑ i∈[m] ( ρi − zi ( log ρi zi + 1 )) . (36)
Proof. Proof of Equation (35). We use α to denote the indices n, di and use limα to denote limn,di→+∞
di/n=zi
. Define fα(λi) = E [ 1 n log det N n ] , f ′α(λi) = ∂ ∂λi E [ 1 n log det N n ] , and
h(λi) = inf ρ∈Rm+
log λ+ ∑
i∈[m]
λiρi + ∑ i∈[m] ( ρi − zi ( log ρi zi + 1 )) . We have∣∣∣∣ 1n log det Nn ∣∣∣∣ ≤ 1n log det ( λIn + λ+ Z>Z n ) = log λ+ 1 n log det ( In + λ+ nλ Z>Z ) .
By Lemma 3, there exists a universal constant C > 0 such that for all sufficiently large n,
1 n log det
( In +
λ+ nλ Z>Z
) ≤ log ( 1 + C
λ
) .
Therefore, we get ∣∣∣∣ 1n log det Nn ∣∣∣∣ ≤ log (λ+ C) .
By the dominated convergence theorem and Lemma 9 (specifically, Equation (27)), we obtain
lim α fα(λi) = h(λi) . (37)
Because ∣∣∣∣ ∂∂λi 1n log det Nn ∣∣∣∣ = 1n tr (Z>i N−1Zi) ≤ 1λn2 tr (Z>i Zi)
and E [
1 λn2 tr ( Z>i Zi )] < +∞, we can interchange the differentiation and the expectation and get
f ′α(λi) = ∂ ∂λi E [ 1 n log det N n ] = E [ ∂ ∂λi 1 n log det N n ] . (38)
Thus we deduce∣∣∣∣ ∂∂λiE [ 1 n log det N n ]∣∣∣∣ = ∣∣∣∣E [ ∂∂λi 1n log det Nn ]∣∣∣∣ ≤ E ∣∣∣∣ ∂∂λi 1n log det Nn ∣∣∣∣ ≤ E [ 1λn2 tr (Z>i Zi) ] .
By Lemma 4, E tr ( Z>i Zi ) n2 and therefore E [ 1 λn2 tr ( Z>i Zi )] . 1λ . The function sequence {f ′α} is uniformly bounded. Then we want to show that {f ′α} is uniformly equicontinuous by showing that {f ′′α} is uniformly bounded. Because∣∣∣∣ ∂2∂λ2i 1n log det Nn ∣∣∣∣ = 1n tr (Z>i N−1Zi)2 ≤ 1nλ2 tr ( Z>i Zi n
)2 and E [ 1 nλ2 tr ( Z>i Zi n )2] < +∞, we can interchange the differentiation and the expectation and
get ∂2 ∂λ2i E [ 1 n log det N n ] = ∂ ∂λi E [ ∂ ∂λi 1 n log det N n ] = E [ ∂2 ∂λ2i 1 n log det N n ] . Therefore, we deduce∣∣∣∣ ∂2∂λ2i E [ 1 n log det N n ]∣∣∣∣ ≤ E ∣∣∣∣ ∂2∂λ2i 1n log det Nn ∣∣∣∣ ≤ 1nλ2E tr ( Z>i Zi n )2 .
Again, by Lemma 4, tr ( Z>i Zi n )2 n. It follows that 1nλ2E tr ( Z>i Zi n )2 . 1λ2 . Therefore {f ′ α} is uniformly equicontinuous.
We want to show limα f ′α(λi) = h ′(λi) by contradiction. If it is not true, there exists > 0 and a subsequence {f ′αk} such that ∣∣f ′αk(λi)− h′(λi)∣∣ ≥ . Let E = [a, b] 3 λi (b > a > 0) be
a closed interval that contains λi. The subsequence {f ′αk} is uniformly bounded and uniformly equicontinuous. By the Arzela-Ascoli theorem, there exists a subsequence { f ′αkj } that converges uniformly on λi ∈ E. Recall limα fα(λi) = h(λi) (Equation (37)). Thus limj fαkj (λi) = h(λi). By (Rudin, 1976, Theorem 7.17), for λi ∈ E, we have
lim j f ′αkj (λi) = h ′(λi) .
This is a contradiction. Hence, we have shown that limα f ′α(λi) = h ′(λi), which is exactly Equation (35) (recall f ′α(λi) = ∂ ∂λi E [ 1 n log det N n ] = E [ ∂ ∂λi 1 n log det N n ] in Equation (38)).
Proof of Equation (36). Define gα(λj) = ∂∂λiE [ 1 n log det N n ] = E [ ∂ ∂λi 1 n log det N n ] . Then g′α(λj) = ∂2 ∂λj∂λi E [ 1 n log det N n ] = ∂∂λjE [ ∂ ∂λi 1 n log det N n ] . We have
∣∣∣∣ ∂2∂λj∂λi 1n log det Nn ∣∣∣∣
= 1 n tr ( ZiZ > i N −1ZjZ > j N −1) = 1
n tr ( Z>i N −1ZjZ > j N −1Zi )
= 1
n ∥∥Z>j N−1Zi∥∥2F ≤ 1 n ‖Zj‖22 ‖Zi‖ 2 2 ∥∥N−1∥∥2 F ≤ 1 λ2n2 ‖Zj‖22 ‖Zi‖ 2 2 .
where the last inequality is because ∥∥N−1∥∥2 F ≤ ∥∥ 1 λnIn ∥∥2 F = 1λ2n . If i 6= j, by Lemma 4, we have
1
λ2n2 E ‖Zj‖22 ‖Zi‖ 2 2 =
1
λ2n2 E ‖Zj‖22 · E ‖Zi‖ 2 2 .
1
λ2 .
If i = j, by Lemma 4, we have
1
λ2n2 E ‖Zi‖42 .
1 λ2n2 · n2 = 1 λ2 .
As a result, we get
∣∣∣∣ ∂2∂λj∂λi 1n log det Nn ∣∣∣∣ . 1n · n2 · 1λ2n = 1λ2 .
Thus we can interexchange ∂∂λj and expectation, and get g ′ α(λj) = E
[ ∂2
∂λj∂λi 1 n log det N n
] . Be-
cause |g′α(λj)| ≤ E ∣∣∣ ∂2∂λj∂λi 1n log det Nn ∣∣∣ . 1λ2 , the function sequence {g′α} is uniformly bounded for λj .
Define L = Z>j N −1Zi and W = Z>j N −1Zj . We have∣∣∣∣∣ ∂3∂λ2j∂λi 1n log det Nn ∣∣∣∣∣
= 2 n tr ( L>WL ) . 1
λn2 tr ( L>Z>j ZjL ) = 1 λn2 tr ( Z>i N
−1 (ZjZ>j )2N−1Zi) . = 1
λn2 ∥∥ZjZ>j N−1Zi∥∥2F
≤ 1 λn2 ∥∥N−1∥∥2 F ∥∥ZjZ>j ∥∥22 ‖Zi‖22 ≤ 1 λ3n3
∥∥ZjZ>j ∥∥22 ‖Zi‖22 = 1
λ3n3 ‖Zj‖42 ‖Zi‖ 2 2 ,
where the first inequality is becauseW 4 1λnZ > j Zj and the third inequality is becauseN −1 4 1λnIn and then ∥∥N−1∥∥2 F ≤ ∥∥ 1 λnIn ∥∥2 F ≤ 1λ2n . By Lemma 4, we have E ‖Zj‖ 4 2 . n
2 and E ‖Zi‖22 . n. If i 6= j, then Zj and Zi are independent, and we deduce
1
λ3n3 E ‖Zj‖42 ‖Zi‖ 2 2 .
1
λ3 .
If i = j, we have 1
λ3n3 E ‖Zi‖42 ‖Zi‖ 2 2 =
1
λ3n3 E ‖Zi‖62 .
1
λ3 .
As a result, we deduce E [ ∂3
∂λ2j∂λi 1 n log det N n ] = ∂∂λjE [ ∂2 ∂λj∂λi 1 n log det N n ] = g′′α(λj). Moreover,
we have
|g′′α(λj)| ≤ E ∣∣∣∣∣ ∂3∂λ2j∂λi 1n log det Nn ∣∣∣∣∣ . 1λ3 .
Therefore {g′α} is uniformly equicontinuous. Define
w(λj) = ∂
∂λi inf ρ∈Rm+
log λ+ ∑
i∈[m]
λiρi + ∑ i∈[m] ( ρi − zi ( log ρi zi + 1 )) . We want to show by contradiction that limα g′α(λj) = w
′(λj). Assume that it is not true. Then there exists > 0 and a subsequence { g′αk } such that ∣∣g′αk(λj)− w′(λj)∣∣ > . Since {g′αk} is uniformly bounded and uniformly equicontinuous, by the Arzela-Ascoli theorem, there is a subsequence { g′αkr } that converges uniformly on a closed intervalE containing λj . Equation (35) shows that limα gα(λj) = w(λj). It follows that limr gαkr (λj) = w(λj). By (Rudin, 1976, Theorem 7.17), for λi ∈ E, we have
lim r g′αkr (λj) = w ′(λj) ,
which is a contradiction. Therefore, we have shown that limα g′α(λj) = w ′(λj), which is exactly Equation (36).
D PROOF OF THEOREM 1
D.1 PROOF OF ITEM 1 Define g(ρ) = log ( λ+ ∑m j=1 λjρj ) + ∑m j=1 ( ρj − zj(log ρjzj + 1) ) . The function g(ρ) is continuously differentiable on Rm+ . The boundary of R+m is ∂Rm+ = {ρ ∈ Rm | (∀i ∈ [m], ρi ≥ 0) ∧ (∃i ∈ [m], ρi = 0)}. Because limRm+3ρ→ρ0∈∂Rm+ g(ρ) = limRm+3ρ→∞ g(ρ) = +∞, there exists a minimizer ρ ∗ ∈ Rm+ of g(ρ).
Taking the derivative with respect to ρi gives
∂
∂ρi
log λ+ m∑
j=1
λjρj + m∑ j=1 ( ρj − zj(log ρj zj + 1) ) = λi λ+ ∑m j=1 λjρj + 1− zi ρi .
Setting it to zero gives Equation (6).
D.2 PROOF OF ITEM 2
Recall Equation (6) λi
λ+ ∑m j=1 λjρ ∗ j + 1− zi ρ∗i = 0 , ∀i ∈ [m] .
Rewriting the above equation gives
(zi − ρ∗i )
( λ+
m∑ k=1 λkρ ∗ k
) = λiρ ∗ i , ∀i ∈ [m] .
Rewriting it in the linear algebraic form yields (z− ρ∗) ( λ+ λ>ρ∗ ) = λ ρ∗ .
Applying ∂∂λ to both sides and using the implicit function theorem, we get (z− ρ∗) ( ρ∗> + λ>J ) − J ( λ+ λ>ρ∗ ) = diag (λ) J + diag (ρ∗) .
Arranging the above equation yields( diag (λ) + ( λ+ λ>ρ∗ ) Im − (z− ρ∗)λ> ) J = (z− ρ∗) ρ∗> − diag (ρ∗) .
Define a = λ + λ>ρ∗, A = diag (λ) + ( λ+ λ>ρ∗ ) Im = diag (λ) + aIm and B = diag (λ) +(
λ+ λ>ρ∗ ) Im − (z− ρ∗)λ> = A− (z− ρ∗)λ>. The matrix determinant lemma gives
det (B) = ( 1− λ>A−1 (z− ρ∗) ) det (A) .
Recall Equation (6) again and we have
λi + a = zia
ρ∗i .
We have
a− ∑ i∈[m] λiρ ∗ i (1− ρ∗i /zi)
= λ+ ∑ i∈[m] λiρ ∗ i − ∑ i∈[m] λiρ ∗ i (1− ρ∗i /zi)
=λ+ ∑ i∈[m] λi (ρ∗i ) 2 zi > 0 .
It follows that ∑ i∈[m] λiρ
∗ i (1− ρ∗i /zi) a < 1 .
Then we compute λ>A−1 (z− ρ∗):
λ>A−1 (z− ρ∗)
= ∑ i∈[m] λi (zi − ρ∗i ) λi + a
= ∑ i∈[m] λi (zi − ρ∗i ) zia ρ∗i
=
∑ i∈[m] λiρ
∗ i (1− ρ∗i /zi) a < 1 .
Thus we get 1− λ>A−1 (z− ρ∗) > 0. Therefore, detB 6= 0 and the matrix B is invertible.
D.3 PROOF OF ITEM 3 Lemma 11. Define N = λnIn + Z>ΛZ, γ = ∑ i∈[m] zi, r = (r1, . . . , rm), λ = (λ1, . . . , λm), and
ϑ(rt, r,λ) = 2rt
√ 1 + ∑ i∈[m] r2i − 2rt ∑ i∈[m] √ ziri + ∑ i∈[m] 1 λi r2i − λr2t .
For any Kt ≥ 2λ and Ku ≥ 2λ+(2+
√ γ)
λ , we have
lim n,di→+∞ di/n→zi trN−1 = lim n,di→+∞ di/n→zi E trN−1
= max 0≤rt≤Kt min 0≤ri≤Ku ϑ(rt, r,λ) = min 0≤ri≤Ku max 0≤rt≤Kt ϑ(rt, r,λ) = max rt≥0 min ri≥0 ϑ(rt, r,λ) = min ri≥0 max rt≥0 ϑ(rt, r,λ) .
(39)
If r∗ is a solution to the optimization problem in Equation (39), then
1 + m∑ j=1 r∗2j =  m∑ j=1 r∗j √ zj + λr ∗ t 2 , (40) r∗t
r∗i√ 1 + ∑m j=1 r ∗2 j = r∗t √ zi − r∗i λi . (41)
Moreover, we have ∂
∂λi max rt≥0 min ri≥0 ϑ(rt, r,λ) = − r∗2i λ2i .
Proof. Let g ∼ N (0, In) be a multivariate standard normal random vector. We have
trN−1
=Egg>N−1g =Eg sup
t∈Rn
( 2g>t− t>Nt ) =Eg sup
t∈Rn
( 2g>t− t>Z>ΛZt− nλ ‖t‖22 ) =Eg sup
t∈Rn inf u∈Rd
( 2g>t− 2u>ΛZt+ u>Λu− nλ ‖t‖22 ) =− 2Eg inf
t∈Rn sup u∈Rd
( u>ΛZt− g>t− 1
2 u>Λu+
1 2 nλ ‖t‖22 ) =− 2Eg inf
t∈Rn sup u∈Rd
( u>Zt− g>t− 1
2 u>Λ−1u+
1 2 nλ ‖t‖22
)
We view
inf t∈Rn sup u∈Rd
( u>Zt− g>t− 1
2 u>Λ−1u+
1 2 nλ ‖t‖22
) (42)
as the primal optimization (PO) problem in the convex Gaussian min-max theorem (CGMT) (Thrampoulidis et al., 2015).
The KKT conditions for Equation (42) give
Z>u− g + nλt = 0 , Zt− Λ−1u = 0 .
Solving the above equations gives
t = N−1g , u = ΛZN−1g .
With probability at least 1− 4 exp(−cn) (c > 0 is a universal constant), we have ‖g‖2 ≤ 2 √ n and
‖Z‖ ≤ √ d+ 2 √ n ≤ ( 2 + √ γ )√ n. Therefore, we get
‖t‖2 ≤ ∥∥N−1∥∥ ‖g‖2 ≤ 1λn · 2√n = 2λ√n ,
‖u‖2 ≤ λ+ ‖Z‖ ‖t‖2 ≤ λ+ (2 + √ γ) √ n · 2
λ √ n
= 2λ+
( 2 + √ γ )
λ .
Write u =  u1... um , where ui ∈ Rdi . For all Kt ≥ 2λ , Ku ≥ 2λ+(2+√γ)λ , the optimal solutions t∗ and u∗ to Equation (42) satisfy √ n ‖t∗‖2 ≤ Kt and ‖ui‖2 ≤ Ku for all
i ∈ [m] with probability at least 1 − 4 exp(−cn). Define St = {t ∈ Rn | √ n ‖t‖2 ≤ Kt} and Su = { u ∈ Rd | ‖ui‖ ≤ Ku,∀i ∈ [m] } . We use α to denote the indices n, di and use limα to denote limn,di→+∞ di/n=zi . Define event
Eα = { inf t∈Rn sup u∈Rd ( u>Zt− g>t− 1 2 u>Λ−1u+ 1 2 nλ ‖t‖22 ) = inf t∈St sup u∈Su ( u>Zt− g>t− 1 2 u>Λ−1u+ 1 2 nλ ‖t‖22 )} .
Then with probability at least 1 − 4 exp(−cn), we have t∗ ∈ St and u∗ ∈ Su. Therefore the event Eα occurs with probability at least 1− 4 exp(−cn), which yields
P {Ecα} ≤ 4 exp(−cn) . Since ∑ n≥1 4 exp(−cn) < +∞, by Borel-Cantelli lemma, we have
P {
lim sup α
Ecα
} = P {( lim inf
α Eα
)c} = 0 .
Then with probability 1, all but finitely many Eα occur. Then almost surely there exists n0 such that for all n > n0, Eα occurs.
The auxiliary optimization (AO) problem is
inf t∈St sup u∈Su
( ‖t‖2 g > 1 u+ ‖u‖2 g > 2 t− g>t− 1
2 u>Λ−1u+
1 2 nλ ‖t‖22
)
= inf 0≤rt≤Kt sup 0≤ri≤Ku
− ∥∥∥g −√∑i∈[m] r2i g2∥∥∥
2√ n rt + rt ∑ i∈[m] ‖g1,i‖2√ n ri − 1 2 ∑ i∈[m] 1 λi r2i + 1 2 λr2t  = inf
0≤rt≤Kt sup
0≤ri≤Ku −√1 + ∑ i∈[m] r2i ‖g3‖2√ n rt + rt ∑ i∈[m] ‖g1,i‖2√ n ri − 1 2 ∑ i∈[m] 1 λi r2i + 1 2 λr2t  , where g1 ∼ N (0, Id), g2 ∼ N (0, In), and g3 ∼ N (0, In).
Taking n, di → +∞ with di/n→ zi constant, the strong law of large numbers gives√ 1 + ∑ i∈[m] r2i ‖g3‖2√ n a.s.→ √ 1 + ∑ i∈[m] r2i ,
‖g1,j‖2√ n = √ dj n ‖g1,j‖2√ dj a.s.→ √zj .
Define Xα (rt, r) = − √ 1 + ∑ i∈[m] r2i ‖g3‖2√ n rt + rt ∑ i∈[m] ‖g1,i‖2√ n ri − 1 2 ∑ i∈[m] 1 λi r2i + 1 2 λr2t .
It is a stochastic process on (rt, r) ∈ [0,Kt]× [0,Ku]m. We have
lim α Xα (rt, r) = X (rt, r) := −rt
√ 1 + ∑ i∈[m] r2i + rt ∑ i∈[m] √ ziri − 1 2 ∑ i∈[m] 1 λi r2i + 1 2 λr2t
almost surely. Since √
1 + x2 is convex and increasing and the function ‖r‖2 is convex, thus √ 1 + ‖r‖22 is convex in r and then − √ 1 + ∑ i∈[m] r 2 i ‖g3‖2√ n rt = − √ 1 + ‖r‖22 ‖g3‖2√ n rt
is concave in r. Because − 12 ∑ i∈[m] 1 λi r2i is concave in r and rt ∑ i∈[m] ‖g1,i‖2√ n ri is linear in r, we deduce that Xα (rt, r) is concave in r. By (Liese & Miescke, 2008, Lemma 7.75), supr∈[0,Ku]m |Xα (rt, r)−X (rt, r)| → 0 almost surely. Then for ∀ > 0, there exists n0( ), d0,i( ), δ0,i( ) such that for all n > n0( ), di > d0,i( ), |di/n− zi| < δ0,i( ) and for all r ∈ [0,Ku]m, we have
X (rt, r)− < Xα(rt, r) < X (rt, r) + .
Thus we obtain
X (rt, r)− < Xα(rt, r) ≤ sup r∈[0,Ku]m Xα(rt, r)
Xα(rt, r) < X (rt, r) + ≤ sup r∈[0,Ku]m X(rt, r) + ,
which in turn implies
sup r∈[0,Ku]m X (rt, r)− ≤ sup r∈[0,Ku]m Xα(rt, r)
sup r∈[0,Ku]m Xα(rt, r) ≤ sup r∈[0,Ku]m X(rt, r) + .
It follows that ∣∣∣supr∈[0,Ku]m Xα(rt, r)− supr∈[0,Ku]m X(rt, r)∣∣∣ ≤ . In other words, we showed
|Yα (rt)− Y (rt)| → 0
almost surely, where Y (rt) := supr∈[0,Ku]m Xα(rt, r) and Y (rt) := supr∈[0,Ku]m X(rt, r).
Because Xα (rt, r) is convex in rt, then Y (rt) = supr∈[0,Ku]m Xα(rt, r) is convex in rt. By (Liese & Miescke, 2008, Lemma 7.75) again, suprt∈[0,Kt] |Yα (rt)− Y (rt)| → 0 almost surely. A similar argument shows that∣∣∣∣ infrt∈[0,Kt]Yα (rt)− infrt∈[0,Kt]Y (rt) ∣∣∣∣ = ∣∣∣∣∣ infrt∈[0,Kt] supr∈[0,Ku]mXα(rt, r)− infrt∈[0,Kt] supr∈[0,Ku]mX(rt, r)
∣∣∣∣∣→ 0 almost surely.
Therefore, we obtain
inf t∈St sup u∈Su
( ‖t‖2 g > 1 u+ ‖u‖2 g > 2 t− g>t− 1
2 u>Λ−1u+
1 2 nλ ‖t‖22
)
= inf 0≤rt≤Kt sup 0≤ri≤Ku −√1 + ∑ i∈[m] r2i ‖g3‖2√ n rt + rt ∑ i∈[m] ‖g1,i‖2√ n ri − 1 2 ∑ i∈[m] 1 λi r2i + 1 2 λr2t  a.s.→ inf
0≤rt≤Kt sup
0≤ri≤Ku −rt√1 + ∑ i∈[m] r2i + rt ∑ i∈[m] √ ziri − 1 2 ∑ i∈[m] 1 λi r2i + 1 2 λr2t  (43) =:µ .
Define event
Aα = {∣∣∣∣ inft∈Rn supu∈Rd ( u>Zt− g>t− 1 2 u>Λ−1u+ 1 2 nλ ‖t‖22 ) − µ ∣∣∣∣ > τ} , Bα = {∣∣∣∣ inft∈St supu∈Su ( u>Zt− g>t− 1 2 u>Λ−1u+ 1 2 nλ ‖t‖22 ) − µ
∣∣∣∣ > τ} , Cα = {∣∣∣∣ inft∈St supu∈Su ( ‖t‖2 g > 1 u+ ‖u‖2 g > 2 t− g>t− 1 2 u>Λ−1u+ 1 2 nλ ‖t‖22 ) − µ
∣∣∣∣ > τ} . Recall
Eα = { inf t∈Rn sup u∈Rd ( u>Zt− g>t− 1 2 u>Λ−1u+ 1 2 nλ ‖t‖22 ) = inf t∈St sup u∈Su ( u>Zt− g>t− 1 2 u>Λ−1u+ 1 2 nλ ‖t‖22 )} .
We have Aα ∩ Eα ⊆ Bα. Equation (43) gives limα P {Cα} = 0 for any τ > 0 because almost sure convergence implies convergence in probability. By the convex Gaussian min-max theorem (Thrampoulidis et al., 2015), we have
P {Bα} ≤ 2P {Cα} . It follows that
P {Aα} ≤ P {Aα ∩ Eα}+ P {Ecα} ≤ P {Bα}+ P {Ecα} ≤ 2P {Cα}+ P {Ecα} . Taking lim supα on both sides, because lim supα P {Bα} ≤ 2 lim supα P {Cα} = 0, we get
lim sup α P {Aα} ≤ lim sup α
P {Ecα} ≤ P {
lim sup α
Ecα
} = 0 ,
where the second inequality is because of the reverse Fatou’s lemma. Thus
inf t∈Rn sup u∈Rd
( u>Zt− g>t− 1
2 u>Λ−1u+
1 2 nλ ‖t‖22
) P→ µ .
Therefore, we deduce
g>N−1g P→ −2 inf
0≤rt≤Kt sup
0≤ri≤Ku −rt√1 + ∑ i∈[m] r2i + rt ∑ i∈[m] √ ziri − 1 2 ∑ i∈[m] 1 λi r2i + 1 2 λr2t  = sup
0≤rt≤Kt inf 0≤ri≤Ku 2rt√1 + ∑ i∈[m] r2i − 2rt ∑ i∈[m] √ ziri + ∑ i∈[m] 1 λi r2i − λr2t  = sup
0≤rt≤Kt inf 0≤ri≤Ku ϑ(rt, r,λ) . (44)
Because ∣∣g>N−1g∣∣ ≤ 1λn ‖g‖22 and E 1λn ‖g‖22 = 1λ < ∞, by the dominated convergence theorem for convergence in probability (Cohn, 2013, Proposition 3.1.6), we get
lim α trN−1 = lim α
Eg [ g>N−1g ] = max
0≤rt≤Kt min 0≤ri≤Ku ϑ(rt, r,λ) . (45)
Note that 2rt √ 1 + ∑ i∈[m] r 2 i is convex in r, −2rt ∑ i∈[m] √ ziri is linear in r, and ∑ i∈[m] 1 λi r2i
is strongly convex in r. Thus ϑ is strongly convex in r. Note that 2rt √ 1 + ∑ i∈[m] r 2 i −
2rt ∑ i∈[m] √ ziri is linear in rt and that −λr2t is strongly concave in rt. Thus ϑ is strongly concave in rt. Then ϑ has a unique saddle point (r∗t , r ∗) on [0,Kt]× [0,Ku]m that satisfies
max rt∈[0,Kt] min r∈[0,Ku]m ϑ (rt, r) = min r∈[0,Ku]m max rt∈[0,Kt]
ϑ (rt, r) = ϑ (r ∗ t , r ∗) , (46)
where the first equality is due to Sion’s minimax theorem. Since ∣∣trN−1∣∣ ≤ 1λ , using the dominated convergence theorem and combining Equation (45) and Equation (46) yields
lim α E trN−1 = max 0≤rt≤Kt min 0≤ri≤Ku ϑ(rt, r,λ) = min 0≤ri≤Ku max 0≤rt≤Kt ϑ(rt, r,λ) .
By the uniqueness of the limit, the right-hand side max0≤rt≤Kt min0≤ri≤Ku ϑ(rt, r,λ) and min0≤ri≤Ku max0≤rt≤Kt ϑ(rt, r,λ) do not depend on Kt and Ku as long as Kt ≥ 2λandKu ≥ 2λ+(2+ √ γ)
λ . Thus we have
lim α E trN−1 = max rt≥0 min ri≥0 ϑ(rt, r,λ) = min ri≥0 max rt≥0 ϑ(rt, r,λ) .
If r∗t = 0, then ϑ (0, r ∗) = minr∈[0,Ku]m ∑ i∈[m] 1 λi r2i = 0. Thus r
∗ must be zero. However, ϑ (
1 2λ ,0 ) = 34λ > ϑ (0, r ∗). Therefore r∗t > 0. We compute the partial derivative
∂ϑ ∂ri = 2rt ri√ 1 + ∑ i∈[m] r 2 i − 2rt √ zi + 2 ri λi .
If r∗i = 0, we have ∂ϑ
∂ri ∣∣∣ ri=0,rt=r∗t = −2r∗t √ zi < 0 .
Therefore, one can increase r∗i and make maxrt∈[0,Kt] minr∈[0,Ku]m ϑ (rt, r) smaller, which results in a contradiction. Thus r∗i > 0. Thus the minimax value is attained when rt, ri > 0 for all i ∈ [m]. To obtain the optimality condition, we compute the partial derivatives
∂ϑ ∂rt = 2
√ 1 + ∑ i∈[m] r2i − 2 ∑ i∈[m] ri √ zi − 2λrt , ∂ϑ ∂ri = 2rt ri√ 1 + ∑ i∈[m] r 2 i − 2rt √ zi + 2 ri λi .
Setting them to zero gives the optimality condition for r∗t , r ∗ 1 , . . . , r ∗ m and yields Equation (40) and Equation (41).
Using the envelope theorem, we get ∂
∂λi max rt∈[0,Kt] min r∈[0,Ku]m ϑ(rt, r,λ)
= ∂ϑ(r∗t , r ∗, λ1, . . . , λm)
∂λi
=− r ∗2 i
λ2i .
Lemma 12. Define N = λnIn + Z>ΛZ. The following equation holds
lim n,di→+∞ di/n→zi
E [ ∂ ∂λi tr ( N−1 )] = ∂2 ∂λi∂λ inf ρ∈Rm+ log λ+ m∑
j=1
λjρj + m∑ j=1 ( ρj − zj(log ρj zj + 1) ) = −r∗2i λ2i ,
where r∗ is a solution to suprt>0 infr1,...,rm>0 ϑ(rt, r1, . . . , rm, λ1, . . . , λm) and
ϑ(rt, r1, . . . , rm, λ1, . . . , λm) = 2rt
√ 1 + ∑ i∈[m] r2i − 2rt ∑ i∈[m] √ ziri + ∑ i∈[m] 1 λi r2i − λr2t .
Proof. Since ∣∣∣∣ ∂∂λi tr (N−1) ∣∣∣∣ = tr (Z>i N−2Zi) ≤ 1(λn)2 tr (Z>i Zi) (47)
and E 1 (λn)2 tr ( Z>i Zi ) 1λ2 (by Lemma 4), using the dominated convergence theorem gives
E [ ∂ ∂λi tr ( N−1 )] = ∂ ∂λi E tr ( N−1 ) .
We use α to denote the indices n, di and use limα to denote limn,di→+∞ di/n=zi . Define fα(λi) = E tr ( N−1 ) , g(λi) = ∂∂λ infρ∈Rm+ [ log ( λ+ ∑ i∈[m] λiρi ) + ∑ i∈[m] ( ρi − zi ( log ρizi + 1 ))] ,
and h(λi) = suprt>0 infr1,...,rm>0 ϑ(rt, r1, . . . , rm, λ1, . . . , λm). Because ∣∣tr (N−1)∣∣ ≤
tr (
1 λnIn ) ≤ 1λ and limα tr ( N−1 ) = g(λi) (by Lemma 9), we have Lemma 9 limα fα(λi) =
limα E tr ( N−1 ) = g(λi). Lemma 11 shows limα fα(λi) = h(λi). Therefore limα fα(λi) = g(λi) = h(λi).
Because of Equation (47), we have f ′α(λi) = ∂ ∂λi
E tr ( N−1 ) and
|f ′α(λi)| = ∣∣∣∣ ∂∂λiE tr (N−1) ∣∣∣∣ = ∣∣∣∣E [ ∂∂λi tr (N−1) ]∣∣∣∣ ≤ E ∣∣∣∣ ∂∂λi tr (N−1) ∣∣∣∣ . 1λ2 and therefore {f ′α} is uniformly bounded for λi. Because∣∣∣∣ ∂2∂λ2i tr (N−1)
∣∣∣∣ =2 tr ( N−1ZiZ > i N −1ZiZ > i N
−1) . 1 λn tr ( N−1 ( ZiZ > i )2 N−1
) = 1
λn tr ( ZiZ > i N −2ZiZ > i ) ≤ 1
(λn) 3 tr
( ZiZ > i )2 ,
and E 1 (λn)3
tr ( ZiZ > i )2 1λ3 (by Lemma 4), using the dominated convergence theorem yields E [ ∂2
∂λ2i tr ( N−1 )] = ∂ ∂λi E [ ∂ ∂λi tr ( N−1 )] = ∂2 ∂λ2i E [ tr ( N−1 )] = f ′′α(λi) .
Moreover, we have
|f ′′α(λi)| ≤ E ∣∣∣∣ ∂2∂λ2i tr (N−1) ∣∣∣∣ . 1λ3 . Thus {f ′α} is uniformly equicontinuous for λi. We want to show that limα f ′α(λi) = g′(λi) by contradiction. Assume that it is not true. Then there exists > 0 and a subsequence { f ′αk }
such that ∣∣f ′αk(λi)− g′(λi)∣∣ > . Since {f ′αk} is uniformly bounded and uniformly equicontinuous for λi ∈ E (E is any closed finite interval containing λi), by the Arzela-Ascoli theorem, there exists a subsequence { f ′αkr } that converges uniformly on E. Since limα fαkr (λi) = g(λi), by (Rudin, 1976, Thoerem 7.17), we have lim r f ′αkr (λi) = g ′(λi) , which yields a contradiction. Therefore, we have limα f ′α(λi) = g ′(λi). Recall g(λi) = h(λi) for any λi > 0. Then by the final part of Lemma 11, we have limα f ′α(λi) = g ′(λi) = h ′(λi) = − r ∗2 i
λ2i .
D.4 BIAS Lemma 13. Suppose that U ∼ ⊕
i∈[m] Unif (O (di)) and V are two independent d × d random matrices such that V d= UV U>, where d = ∑m i=1 di. Let θ ∈ Rd be a fixed vector. Write
θ =  θ1... θm , where θi ∈ Rdi . Let φ ∼
⊕ i∈[m] Unif ( Sdi−1 (‖θi‖2) ) be a random vector independent of V and let Λ = diag (λ1Id1 , . . . , λmIdm) ∈ Rd×d. Then we have
E [∥∥∥Λ1/2V θ∥∥∥2
2
] = E [∥∥∥Λ1/2V φ∥∥∥2 2 ] .
Proof. Recall UΛU> = Λ and noticing U>θ d= φ, we get E [∥∥∥Λ1/2V θ∥∥∥2
2 ] =E [∥∥∥Λ1/2UV U>θ∥∥∥2 2
] =E [ θ>UV >U>ΛUV U>θ
] =E [ θ>UV >ΛV U>θ
] =E [∥∥∥Λ1/2V U>θ∥∥∥2 2
] =E [∥∥∥Λ1/2V φ∥∥∥2 2 ] .
Lemma 14. Define Θ̃ = diag ( ‖θ′1‖22/d1Id1 , . . . , ‖θ′m‖22/dmIdm ) and S =
Λ1/2Z ( nλIn + Z >ΛZ )−1 Z>Λ1/2. Then we have
Bλ,d,n = ‖θ∗‖2Σ − 2E tr ( ΛSΘ̃ ) + E tr ( SΛSΘ̃ ) .
Proof. Recall Equation (22) in Lemma 5
Bλ,d,n = E [ ‖Λ1/2 ( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 θ′‖22 ] .
Let U ∼ ⊕
i∈[m] Unif (O(di)) be a random matrix independent of Z. Because UZ d = Z, we have
Id + 1
nλ Λ1/2ZZ>Λ1/2 d = Id +
1
nλ Λ1/2UZZ>U>Λ1/2 = U
( Id + 1
nλ Λ1/2ZZ>Λ1/2
) U> .
Define θ̃ ∼ ⊕
i∈[m] Unif(Sdi−1(‖θ′i‖2)). Lemma 13 gives
Bλ,d,n =E [ ‖Λ1/2 ( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 θ̃‖22 ] =E [ ‖Λ1/2 (Id − S) θ̃‖22
] =E [ ‖ (Id − S) θ̃‖2Λ
] =E ∥∥∥θ̃∥∥∥2 Λ − E [ θ̃>ΛSθ̃ ] − E [ θ̃>SΛθ̃ ] + E [ θ̃>SΛSθ̃ ] .
Notice that ∥∥∥θ̃∥∥∥2
Λ = ‖θ′‖2Λ and Θ̃ = E
[ θ̃θ̃> ] . Because Θ̃ commutes with Λ, we have tr ( SΛΘ̃ ) =
tr ( SΘ̃Λ ) = tr ( ΛSΘ̃ ) . In light of these, we deduce
Bλ,d,n = ‖θ′‖ 2 Λ − E tr
( ΛSΘ̃ ) − E tr ( SΛΘ̃ ) + E tr ( SΛSΘ̃ ) = ‖θ′‖2Λ − 2E tr ( ΛSΘ̃ ) + E tr ( SΛSΘ̃ ) .
Lemma 14 expresses the bias Bλ,d,n as the sum of three terms.
Computing ‖θ′‖2Λ Note that ‖θ′‖ 2 Λ = θ ′>Λθ′ = ∑ i∈[m] λi ‖θ′i‖ 2 2. Therefore,
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi
‖θ′‖2Λ = q > (λ z) .
Computing E tr ( ΛSΘ̃ ) Define N = λnIn + Z>ΛZ = λnIn + ∑ i∈[m] λiZiZ > i . We have
E tr ( ΛSΘ̃ )
=E tr ( Z>Λ1/2Θ̃Λ3/2Z ( nλIn + Z >ΛZ )−1)
=E tr ( Z>Λ2Θ̃ZN−1 ) = ∑ i∈[m] λ2i ‖θ′i‖22 di E tr ( ZiZ > i N
−1) = ∑ i∈[m] λ2i ‖θ′i‖22 n di E [ ∂ ∂λi 1 n log det N n ]
= ∑ i∈[m] λ2i ‖θ′i‖22 n di ∂ ∂λi E [ 1 n log det N n ] ,
where the second inequality is because Θ̃ commutes with Λ3/2 and the final equality is because of Equation (38). Taking lim n,di→+∞
di/n→zi ‖Πiθ∗‖2→ηi
and using Lemma 10 gives
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi
E tr ( ΛSΘ̃ ) = ∑ i∈[m] λ2i η 2 i zi ∂ ∂λi inf ρ∈Rm+
log λ+ ∑
i∈[m]
λiρi + ∑ i∈[m] ( ρi − zi ( log ρi zi + 1 )) . Using the envelope theorem yields
∂
∂λi inf ρ∈Rm+
log λ+ ∑
i∈[m]
λiρi + ∑ i∈[m] ( ρi − zi ( log ρi zi + 1 )) = ρ∗i λ+ ∑ i∈[m] λiρ ∗ i = zi − ρ∗i λi ,
where the final equality is because of Equation (6) in Item 1. Therefore, we deduce
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi
E tr ( ΛSΘ̃ ) = ∑ i∈[m] λ2i η 2 i zi zi − ρ∗i λi = ∑ i∈[m] λiη 2 i ( 1− ρ ∗ i zi ) = q> (λ (z− ρ∗)) .
Computing E tr ( SΛSΘ̃ ) We have
E tr ( SΛSΘ̃ ) =E tr [ Λ1/2ZN−1Z>Λ2ZN−1Z>Λ1/2Θ̃
] =E tr [ Z>Λ1/2Θ̃Λ1/2ZN−1Z>Λ2ZN−1
] =E tr [ Z>ΛΘ̃ZN−1Z>Λ2ZN−1
] = ∑ i∈[m] λi‖θ′i‖22 di ∑ j∈[m] λ2jE tr [ ZiZ > i N −1ZjZ > j N
−1] =−
∑ i∈[m] λi‖θ′i‖22n di ∑ j∈[m] λ2jE [ ∂2 ∂λj∂λi 1 n log det N n ] ,
where the third equality is because Θ̃ commutes with Λ1/2. Taking lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi and using Lemma 10 gives
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi
E tr ( SΛSΘ̃ ) = − ∑ i∈[m] λiη 2 i zi ∑ j∈[m] λ2j ∂2 ∂λj∂λi inf ρ∈Rm+
log λ+ ∑
l∈[m]
λlρl + ∑ l∈[m] ( ρl − zl ( log ρl zl + 1 )) . Write λ = (λ1, . . . , λm)> and z = (z1, . . . , zm)>. Let ρ∗ ∈ Rm be a minimizer of Equation (5) and J = ∂ρ ∗
∂λ ∈ R m×m be the Jacobian matrix Jij = ∂ρ∗i ∂λj . Recall Item 2( diag (λ) + ( λ+ λ>ρ∗ ) Im − (z− ρ∗)λ> ) J = (z− ρ∗) ρ∗> − diag (ρ∗) .
Using the envelope theorem, we have
∂
∂λi inf ρ∈Rm+
log λ+ ∑
l∈[m]
λlρl + ∑ l∈[m] ( ρl − zl ( log ρl zl + 1 )) = ρ∗i λ+ ∑ l∈[m] λlρ ∗ l = ρ∗i λ+ λ>ρ∗ .
Recall Equation (6) yields ρ∗i
λ+ λ>ρ∗ = zi − ρ∗i λi .
Differentiating the above equation with respect to λj gives
∂2
∂λj∂λi inf ρ∈Rm+
log λ+ ∑
l∈[m]
λlρl + ∑ l∈[m] ( ρl − zl ( log ρl zl + 1 )) = ∂
∂λj zi − ρ∗i λi
= −λiJij − (zi − ρ∗i ) δij
λ2i .
It follows that
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi
E tr ( SΛSΘ̃ )
= ∑ i∈[m] λiη 2 i zi ∑ j∈[m] λ2j λiJij + (zi − ρ∗i ) δij λ2i
= ∑
i,j∈[m]
qiλ 2 j ( Jij +
(zi − ρ∗i ) δij λi ) =q> ( λ (z− ρ∗) + Jλ 2 )
Putting all three terms together, we have
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi
Bλ,d,n = q > (λ z)−2q> (λ (z− ρ∗))+q> ( λ (z− ρ∗) + Jλ 2 ) = q> ( λ ρ∗ + Jλ 2 ) .
Since {Bλ,d,n} is uniformly bounded and uniformly equicontinuous for λ ∈ (0, 1] by Lemma 5, {Bλ,d,n} can be extended continuously to [0, 1] and the family of extended functions is still uniformly bounded and uniformly equicontinuous for λ ∈ [0, 1]. By the Arzela-Ascoli theorem, {Bλ,d,n} converges uniformly to the limit. By the Moore-Osgood theorem, we can exchange the two limits lim n,di→+∞
di/n→zi ‖Πiθ∗‖2→ηi
and limλ→0+ and get
lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi B0,d,n = lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi lim λ→0+ Bλ,d,n = lim λ→0+ lim n,di→+∞ di/n→zi ‖Πiθ∗‖2→ηi
Bλ,d,n = q > (λ ρ∗ + Jλ 2) |λ=0 .
D.5 VARIANCE
Define N = nλIn + Z>ΛZ. Recalling Lemma 6 gives Vλ,d,n
=σ2E‖ΛZN−1‖22
=σ2 m∑ i=1 λ2iE tr ( ZiZ > i N −2) =− σ2
m∑ i=1 λ2iE [ ∂ ∂λi tr ( N−1 )] .
Using Lemma 12, we get lim
n,di→∞ di/n→zi
Vλ,d,n
=− σ2 m∑ i=1 λ2i lim n,di→∞ di/n→zi E [ ∂ ∂λi tr ( N−1 )]
=− σ2 m∑ i=1 λ2i lim n,di→∞ di/n→zi E [ ∂ ∂λi tr ( N−1 )] =− σ2 m∑ i=1 λ2i ∂2 ∂λi∂λ inf ρ∈Rm+ log λ+ m∑ j=1 λjρj + m∑ j=1 ( ρj − zj(log ρj zj + 1)
) . Using the envelope theorem, we deduce
∂
∂λ inf ρ∈Rm+
log λ+ m∑
j=1
λjρj + m∑ j=1 ( ρj − zj ( log ρj zj + 1 )) = 1 λ+ ∑m j=1 λjρ ∗ j .
Then we take ∂∂λi and obtain
∂2
∂λi∂λ inf ρ∈Rm+
log λ+ m∑
j=1
λjρj + m∑ j=1 ( ρj − zj(log ρj zj + 1) ) = ∂
∂λi
1 λ+ ∑m j=1 λjρ ∗ j
=− ρ∗i + ∑ j∈[m] λjJji(
λ+ ∑ j∈[m] λjρ ∗ j )2 .
As a result,
lim n,di→∞ di/n→zi Vλ,d,n = σ 2 m∑ i=1 λ2i ρ∗i + ∑ j∈[m] λjJji( λ+ ∑ j∈[m] λjρ ∗ j )2 = σ2 ( λ 2 )> ( ρ∗ + J>λ ) (λ+ λ>ρ∗) 2 .
By Lemma 12, the variance is given by
lim n,di→∞ di/n→zi Vλ,d,n = lim n,di→∞ di/n→zi −σ2 m∑ i=1 λ2iE [ ∂ ∂λi tr ( N−1 )] = σ2 m∑ i=1 r∗2i ,
where r∗ solves
sup rt>0 inf r1,...,rm>0 2rt√1 + ∑ i∈[m] r2i − 2rt ∑ i∈[m] √ ziri + ∑ i∈[m] 1 λi r2i − λr2t  . Since {Vλ,d,n} is uniformly bounded and uniformly equicontinuous with respect to λ ∈ (0, 1] by Lemma 5, {Vλ,d,n} can be extended continuously to [0, 1] and the family of extended functions is still uniformly bounded and uniformly equicontinuous. By the Arzela-Ascoli theorem, {Vλ,d,n} converges uniformly to the limit. By the Moore-Osgood theorem, we can exchange the two limits limn,di→∞
di/n→zi and limλ→0+ and get
lim n,di→∞ di/n→zi lim λ→0+ Vλ,d,n = lim λ→0+ lim n,di→∞ di/n→zi Vλ,d,n = σ 2 m∑ i=1 r∗2i |λ=0 .
E PROOF OF THEOREM 2
We use Theorem 1 to prove Theorem 2. As in Theorem 1, let r∗ solve minri≥0 maxrt≥0 ϑ(rt, r,λ), where ϑ is defined in Equation (7). Note that ϑ is a quadratic function of rt. Define A =√∑
i∈[m] r 2 i + 1, B = ∑ i∈[m] √ ziri, A∗ = √∑ i∈[m] r ∗2 i + 1, and B ∗ = ∑ i∈[m] √ zir ∗ i . Then
r∗t = A−B λ and we get
min ri≥0 max rt≥0 ϑ(rt, r,λ) = min ri≥0  (A−B)2 λ + ∑ i∈[m] 1 λi r2i  = min ri≥0  (A−B)2 λ + ∑ i∈[m] 1 λi r2i  . Taking the partial derivative with respect to ri gives
∂
∂ri
 (A−B)2 λ + ∑ i∈[m] 1 λi r2i  = 2 · A−B λ (ri A − √ zi ) + 2 · ri λi .
Setting it to zero gives the optimality condition for r∗i :
A∗ −B∗
λ
( r∗i A∗ − √ zi ) = −r ∗ i λi , i ∈ [m] . (48)
It follows that r∗i A∗ − √ zi
r∗j A∗ − √ zj
= r∗i /λi r∗j /λj , i, j ∈ [m] .
Some algebraic manipulation in the above equation yields
r∗i r∗j = λi λj · √ ziA ∗ − r∗i√ zjA∗ − r∗j , i, j ∈ [m] .
Define z = (z1, . . . , zm). Then ‖z‖1 = ∑ i∈[m] zi. By Cauchy–Schwarz inequality, if d/n →∑
i∈[m] zi < 1
B ≤ √∑ i∈[m] zi ‖r‖2 < ‖r‖2 < √ ‖r‖22 + 1 = A .
Thus there does not exist r such that A = B. If d/n → ∑ i∈[m] zi > 1, then A = B is feasible for r. For example, set
r = 1√
(‖z‖1 − 1) ‖z‖1
√ z .
We have
B = 〈 r, √ z 〉 = √ ‖z‖1 ‖z‖1 − 1
A = √ 1 + ‖r‖22 = √ 1 +
1
‖z‖1 − 1 = B .
If ‖z‖1 > 1, since A = B is feasible, then
lim λ→0+ min ri≥0  (A−B)2 λ + ∑ i∈[m] 1 λi r2i  = min ri≥0 A=B ∑ i∈[m] 1 λi r2i .
If ‖z‖1 < 1, then A−B always holds. To be precise, we have A−B ≥ (√ ‖r‖22 + 1− ‖r‖2 ) ∨ (( 1− √ ‖z‖1 ) ‖r‖2 ) . If ‖r‖2 > 1, then ( 1− √ ‖z‖1 ) ‖r‖2 > 1− √ ‖z‖1. If ‖r‖2 ≤ 1, then √ ‖r‖22 + 1−‖r‖2 ≥ √ 2−1.
Thus there exists a universal constant C0 = ( 1− √ ‖z‖1 ) ∨ (√ 2− 1 ) > 0 such that
A−B ≥ C0 .
Recall Equation (48). We have (A∗ −B∗) ( r∗i A∗ − √ zi ) = −λr ∗ i λi , i ∈ [m] .
Taking limλ→0+ , since A∗ −B∗ ≥ C0 does not go to zero, we have
r∗i A∗ − √ zi = 0 , i ∈ [m] .
Then we get
r∗2i 1 + ∑ j∈[m] r ∗2 j = zi , i ∈ [m] .
Summing all i ∈ [m] yields
‖z‖1 = ∑ i∈[m] r ∗2 i
1 + ∑ i∈[m] r ∗2 i .
Therefore, we have
lim n,di→+∞ di/n→zi V0,d,n = σ 2 lim λ→0+ m∑ i=1 r∗2i = σ 2
∑ i∈[m] zi
1− ∑ i∈[m] zi .
F PROOF OF THEOREM 3 Define A∗ = √∑
i∈[m] r ∗2 i + 1 and B ∗ = ∑ i∈[m] √ zir ∗ i . Equation (15) in Theorem 2 yields
r∗1 r∗2 = λ1 λ2 · √ z1A ∗ − r∗1√ z2A∗ − r∗2 .
Using the constraint A∗ = B∗, we get
r∗1 r∗2
= λ1 (√ z1B ∗ − r∗1 )
λ2 (√ z2B∗ − r∗2 ) . Define q = r ∗ 1
r∗2 . We have the following equation q = λ1 ( q(z1 − 1) + √ z1z2 ) λ2 ( q √ z1z2 + z2 − 1
) . Solving the above equation yields
q = λ1 (z1 − 1) + λ2 (1− z2) + √ (λ1 (z1 − 1) + λ2 (1− z2)) 2 + 4λ1λ2z1z2 2λ2 √ z1z2 . (49)
Here we discard the negative root. Let x = r∗21 + r ∗2 2 = r ∗2 2
( 1 + q2 ) . ?? yields
1 + x = r∗22 (q √ z1 + √ z2) 2 = x 1 + q2 (q √ z1 + √ z2) 2 .
Solving x from the above equation gives
x = q2 + 1
q2(z1 − 1) + 2q √ z1z2 + z2 − 1 .
Therefore,
lim n,di→+∞ di/n→zi
V0,d,n = q2 + 1
q2(z1 − 1) + 2q √ z1z2 + z2 − 1 .
G PROOF OF THEOREM 4
Instead of considering the θ∗ specified in Equation (18), we first consider a Bayesian setting where θ∗ ∼ N ( 0, 1dId ) . Later, we will show that the setup in Equation (18) is asymptotically (as di →∞) equivalent to this Bayesian setting. The precise meaning of equivalence will also be presented later. Our strategy can be divided into two steps. The first step is to show that the Bayes risk of the Bayes estimator is monotonically decreasing in the sample size n. The second step is to translate the sample-wise monotonicity of the Bayes estimator to the excess risk of the optimally regularized estimator θ̂λ,d,n in the setup of Equation (18).
Recall that since we are interested in sample-wise monotonicity, we add a subscript n to X and y (they are defined by Equation (1) in Section 1.1) to emphasize that they consist of n data items. In this Bayesian setting, the likelihood function of θ∗ is L (θ∗ | Xn,yn) = ∏ i∈[n] L (θ∗ | xi, yi) ∝ exp ( − ∑ i∈[n] (yi − 〈θ∗, xi〉) 2 2σ2 ) = exp ( − ‖Xnθ∗ − yn‖22 2σ2 ) .
The density of the prior of θ∗ is proportional to exp ( −d2 ‖θ ∗‖22 )
. Therefore, the posterior density of θ∗ is given by
p (θ∗ | Xn,yn) ∝ exp ( − d ‖θ∗‖22
2 − ‖Xnθ∗ − yn‖22 2σ2
) .
As a result, the posterior distribution of θ∗ is Gaussian. The Bayes estimator is
θ̂Bayes (Xn,yn) = arg min θ
Eθ∗∼p(θ∗|Xn,yn) ‖θ − θ ∗‖2Σ .
Taking the derivative with respect to θ gives
∂
∂θ Eθ∗∼p(θ∗|Xn,yn) ‖θ − θ ∗‖2Σ = 2Σ (θ − θ ∗) .
Setting the above equation to zero yields Σ ( θ̂Bayes (Xn,yn)− Eθ∗∼p(θ∗|Xn,yn)θ∗ ) = 0 and there-
fore
θ̂Bayes (Xn,yn) = Eθ∗∼p(θ∗|Xn,yn)θ ∗ = E [θ∗ | Xn,yn] = arg min
θ
( d ‖θ∗‖22 + ‖Xnθ∗ − yn‖22 σ2 ) .
The final equality is because the posterior mean of a Gaussian distribution equals its mode.
Define the Bayes risk
Rn , Eθ∗∼N (0, 1d Id),Xn,yn [∥∥∥θ̂Bayes (Xn,yn)− θ∗∥∥∥2 Σ ] .
Write X = Rd and Y = R. Define
R′n , inf θ̂:Xn×Yn→R
Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥θ̂ (Xn,yn)− θ∗∥∥∥2 Σ .
We have
R′n = inf θ̂:Xn×Yn→Rd
Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥Σ1/2θ̂ (Xn,yn)− Σ1/2θ∗∥∥∥2
2
= inf θ̂:Xn×Yn→Rd
Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥θ̂ (Xn,yn)− Σ1/2θ∗∥∥∥2
2 = Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥E [Σ1/2θ∗ | Xn,yn]− Σ1/2θ∗∥∥∥2
2 = Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥Σ1/2θ̂Bayes (Xn,yn)− Σ1/2θ∗∥∥∥2
2
= Rn .
where the third equality is because the conditional expectation minimizes the `2 loss. Next, we want to show that Rn+1 ≤ Rn, i.e., the Bayes risk of the Bayes estimator is monotonically decreasing in the sample size n.
Rn+1 = inf θ̂:Xn+1×Yn+1→Rd Eθ∗∼N (0, 1d Id),Xn+1,yn+1 [∥∥∥θ̂ (Xn+1,yn+1)− θ∗∥∥∥2 Σ ] ≤ inf θ̂:Xn×Yn→R Eθ∗∼N (0, 1d Id),Xn+1,yn+1 [∥∥∥θ̂ (Xn,yn)− θ∗∥∥∥2 Σ
] = Rn .
Then we want to show that Rn equals the Bayes risk of the optimally regularized estimator θ̂λ,n,d:
Rn = inf λ≥0 Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ .
Since Rn = inf θ̂:Xn×Yn→Rd Eθ∗,Xn,yn ∥∥∥θ̂ (Xn,yn)− θ∗∥∥∥2
Σ , we get
Rn ≤ inf λ≥0 Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ .
On the other hand, recalling θ̂Bayes (Xn,yn) = arg minθ ( d ‖θ‖22 + ‖Xnθ−yn‖22 σ2 ) = θ̂σ2d
n ,n,d and
Rn = Eθ∗∼N (0, 1d Id),Xn,yn [∥∥∥θ̂Bayes (Xn,yn)− θ∗∥∥∥2 Σ ] , we deduce
Rn ≥ inf λ≥0 Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ .
Therefore we deduce Rn = infλ≥0 Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ . As a result, we establish sample-wise monotonicity of the Bayes risk of optimal regularized θ̂λ,n,d:
Rn+1 = inf λ≥0
Eθ∗,Xn+1,yn+1 ∥∥∥θ̂λ,n+1,d − θ∗∥∥∥2
Σ ≤ inf λ≥0
Eθ∗,Xn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2
Σ = Rn . (50)
In what follows, we show that if θ∗ is given by Equation (18), the excess risk of θ̂λ,n,d is asymptotically equal to its Bayes risk when θ∗ ∼ N (0, 1dId):
lim di→∞ ∣∣∣∣EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2Σ − Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2Σ ∣∣∣∣ = 0 .
We abuse the notation in the above equation. The θ∗ in EXn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2
Σ satisfies Equation (18), while the θ∗ in Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ follows a normal distribution N (0, 1dId). By Lemma 5 and Lemma 6, if Σ = PΛP> and θ′ = P>θ∗ are as defined in Table 1 (where P is an orthogonal matrix and Λ = diag(λ1Id1 , . . . , λmIdm) ∈ Rd×d is a diagonal matrix), for fixed θ∗ we have
EXn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2
Σ =EXn,yn
[ ‖Λ1/2 ( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 θ′‖22 ] + σ2EXn,yn [ ‖ΛZ ( λnIn + Z >ΛZ )−1 ‖22] ,
where every entry ofZ ∈ Rd×n follows i.i.d. N (0, 1). If θ∗ ∼ N (0, 1dId), we have θ ′ ∼ N (0, 1dId). Since the variance term σ2EXn,yn [ ‖ΛZ ( λnIn + Z >ΛZ )−1 ‖22] does not depend on θ∗, the two
variance terms cancel out and we get EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2
Σ − Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ
=EXn,yn
[ ‖Λ1/2 ( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 θ′‖22 ]
− EXn,yn,θ′∼N (0, 1d Id)
[ ‖Λ1/2 ( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 θ′‖22 ]
For U ∼ ⊕
i∈[m] Unif (O (di)), we have( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 d = ( Id + 1
nλ Λ1/2UZZ>U>Λ1/2
)−1 = U ( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 U> .
By Lemma 13, for θ∗ (and thereby θ′) specified in Equation (18), we get
EXn,yn
[ ‖Λ1/2 ( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 θ′‖22 ] = EXn,yn,φ [ ‖Λ1/2 ( Id + 1
nλ Λ1/2ZZ>Λ1/2
)−1 φ‖22 ] ,
where φ ∼ ⊕ i∈[m] Unif ( Sdi−1 (‖θ′i‖2) ) = ⊕ i∈[m] Unif ( Sdi−1 (√ di/d )) .
In the Bayesian setting, if θ′ ∼ N (0, 1dId), then U >θ′ ∼ N (0, 1dId). We have Eθ′∼N (0, 1d Id),Xn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2 Σ = Eθ′∼N (0, 1d Id),Xn,yn ∥∥∥∥∥Λ1/2U ( Id + 1 nλ Λ1/2ZZ>Λ1/2 )−1 U>θ′ ∥∥∥∥∥ 2
2  = EXn,yn,ψ ∥∥∥∥∥Λ1/2 ( Id + 1 nλ Λ1/2ZZ>Λ1/2 )−1 ψ ∥∥∥∥∥ 2
2  , where ψ = U>θ′ ∼ N (0, 1dId).
Next, we want to couple φ and ψ. Let si i.i.d.∼ Unif ( Sdi−1(1) ) , hi i.i.d.∼ χ2(di), and define
φ =  √ d1/ds1
...√ dm/dsm
 , ψ =  √ h1/ds1
...√ hm/dsm  . We have ‖φ‖2 = 1 and
‖ψ‖2 = √√√√ m∑ i=1 hi d = √√√√ m∑ i=1 di d · hi di ,
‖φ− ψ‖2 = √√√√ m∑ i=1 di d ( 1− √ hi di )2 .
By the strong law of large numbers, limdi→+∞ hi/di = 1 almost surely. Thus we get limdi→+∞,di/d→νi ‖ψ‖2 = √∑m i=1 νi and limdi→+∞,di/d→νi ‖φ− ψ‖2 = 0 almost surely (recall that we will let di → +∞ and di/d → νi for some constant νi > 0. ). Because∥∥Λ1/2∥∥ 2 . 1 and ∥∥∥(Id + 1nλΛ1/2ZZ>Λ1/2)−1∥∥∥
2 ≤ ‖Id‖2 = 1, we bound the norm of Q , Λ1/2 ( Id + 1 nλΛ 1/2ZZ>Λ1/2 )−1 as follows
‖Q‖2 ≤ ∥∥∥Λ1/2∥∥∥
2 ∥∥∥∥∥ ( Id + 1 nλ Λ1/2ZZ>Λ1/2 )−1∥∥∥∥∥ 2 . 1 .
It follows that ∣∣∣EXn,yn,φ [‖Qφ‖22]− EXn,yn,ψ [‖Qψ‖22]∣∣∣ ≤EXn,yn,φ,ψ
∣∣∣‖Qφ‖22 − ‖Qψ‖22∣∣∣ =EXn,ynφ,ψ (‖Qφ‖2 + ‖Qψ‖2) |‖Qφ‖2 − ‖Qψ‖2| .EXn,ynφ,ψ [(‖φ‖2 + ‖ψ‖2) ‖Q(φ− ψ)‖2] .EXn,ynφ,ψ ‖φ− ψ‖2 ,
where the last inequality is because ‖φ‖2 + ‖ψ‖2 . 1 for all sufficiently large di. We know that limdi→+∞,di/d→νi ‖φ− ψ‖2 = 0 almost surely. To apply Lebesgue’s dominated convergence theorem, we need to find a dominating integrable random variable. In fact, 1 + ‖ψ‖2 dominates ‖φ− ψ‖2:
‖φ− ψ‖2 ≤ ‖φ‖2 + ‖ψ‖2 = 1 + ‖ψ‖2 . It is integrable because E ‖ψ‖2 = E [√ χ2(d) d ] ≤ √ E[χ2(d)] d = 1. Application of Lebesgue’s
dominated convergence theorem yields
lim di→+∞,di/d→νi
∣∣∣EXn,yn,φ [‖Qφ‖22]− EXn,yn,ψ [‖Qψ‖22]∣∣∣ = 0 .
Therefore, we conclude that
lim di→+∞,di/d→νi ∣∣∣∣EXn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2Σ − Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2Σ ∣∣∣∣ = 0
and this convergence is uniform in n and λ ∈ (0,∞). It follows that
lim n,di→∞ n/di→γi
∣∣∣∣EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ − Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ ∣∣∣∣ = 0 and this convergence is uniform in λ ∈ (0,∞). By Lemma 8 (the proof is similar when we replace α→ +∞ by n, di →∞, n/di → γi), we have
lim n,di→∞ n/di→γi
∣∣∣∣ infλ>0EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2Σ − infλ>0Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2Σ ∣∣∣∣ = 0 . (51)
Define fα(λ) = Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ . We use α to denote the indices n, di. By Lemma 5 and Lemma 6, we have
Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ . Eθ∗∼N (0, 1d Id) ‖θ ∗‖22 = 1 .
Therefore {fα(λ)} is uniformly bounded for λ > 0. Since ∣∣∣∣ ddλEXn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2
Σ ∣∣∣∣ . ‖θ∗‖22 and Eθ∗∼N (0, 1d Id) ‖θ
∗‖22 = 1, we have∣∣∣∣ ddλEθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2Σ ∣∣∣∣ = ∣∣∣∣Eθ∗∼N (0, 1d Id) ddλEXn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2Σ ∣∣∣∣ . 1 . As a result, {fα(λ)} is uniformly equicontinuous for λ > 0, and in particular λ ∈ (0,M ] for any M > 0. Therefore {fα(λ)} can be extended continuously to [0,M ] and the family of extended functions is still uniformly bounded and uniformly equicontinuous. Recall that if θ∗ ∼ N (0, 1dId), we have θ′ ∼ N (0, 1dId). As in Equation (19), write θ ′ in a row-partitioned form
θ′ =  θ ′ 1 ... θ′m  , where θ′i ∈ Rdi . Then ‖Πiθ∗‖2 = ‖θ′i‖2 ∼ √ χ2(di) d = √ χ2(di) di · did → √ νi as n, di → +∞
and n/di → γi, where νi = ( γi ∑ j∈[m] 1 γj )−1 . By Theorem 1, {fα(λ)} converges pointwise, say, to h(λ, γ1, . . . , γm). By the Arzela-Ascoli theorem, limα fα(λ) = h(λ, γ1, . . . , γm) uniformly on λ ∈ [0,M ]. Therefore, as n, di →∞ and n/di → γi, by Lemma 8, we have
inf λ∈[0,M ] Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2 Σ → inf λ∈[0,M ] h(λ, γ1, . . . , γm) .
Recalling θ̂Bayes (Xn,yn) = arg minθ ( d ‖θ‖22 + ‖Xnθ−yn‖22 σ2 ) = θ̂σ2d
n ,n and
Rn = Eθ∗∼N (0, 1d Id),Xn,yn [∥∥∥θ̂Bayes (Xn,yn)− θ∗∥∥∥2 Σ ] = inf λ≥0 Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ .
For all M > CM := 2σ2 ∑ i∈[m] 1 γi ≥ σ 2d n (recall d n → ∑ i∈[m] 1 γi ), we have
inf λ∈[0,M ] Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ = inf λ≥0 Eθ∗∼N (0, 1d Id)EXn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ
→ inf λ∈[0,M ] h(λ, γ1, . . . , γm) .
The uniqueness of limits implies that infλ∈[0,M ] h(λ, γ1, . . . , γm) is independent of M as long as M > σ2. As a result, if M > CM , we have infλ∈[0,M ] h(λ, γ1, . . . , γm) = infλ≥0 h(λ, γ1, . . . , γm), which yields
inf λ≥0 Eθ∗∼N (0, 1d Id),Xn,yn ∥∥∥θ̂λ,n,d − θ∗∥∥∥2 Σ → inf λ≥0 h(λ, γ1, . . . , γm) . (52)
Equation (50) implies infλ≥0 h(λ, γ1, . . . , γm) is decreasing in every γi. Combining Equation (51) and Equation (52) gives
inf λ≥0
EXn,yn ∥∥∥θ̂λ,n − θ∗∥∥∥2
Σ → inf λ≥0 h(λ, γ1, . . . , γm) .

Review summaries:
This paper studies the generalization risk of ridge and ridgeless linear regression under Gaussian data and Gaussian noise settings. The limiting bias and variance are described in closed-form. The analysis leverages random matrix theory and CGMT. The authors show that sample-wise multiple descent phenomenon can take place or be avoided if optimal regularization is applied.

Authors study the generalization risk of ridge and ridgeless linear regression. It is assumed that the data follows a multivariate normal distribution and there is also an assumption on the spectrum of the covariance matrix. Among the contributions, the authors show a formula for the limiting bias and variance, show that sample-wise multiple descent happens when the covariance matrix is highly ill-conditioned, and study optimal regularization.

The paper builds on a very recent line of work studying the “multiple descent” phenomenon.  Multiple descent is known to occur along many axes: (i) number of parameters, (ii) train time, (iii) amount of train data. The current paper focuses on studying the theoretical rates for (iii) in the asymptotic regime where the number of samples and the dimensionality of data tend to infinity at similar rates.  This covers both under and over-parameterized settings. In this limiting regime, the main results of the paper are: 1.	Computing formulae for the excess risk of the least squares estimator and the ridge regression estimator in a specific setting. 2.	Using 1, they show that if the covariance matrix of data is ill-conditioned, then multiple descent can happen. 3.	They also show that under an assumption on the scaling of the projections of the true parameter to the eigen-spaces of the eigenvalues of the covariance matrix, an optimal amount of regularization ridge regression estimator avoids multiple descent even in anisotropic settings which is novel.