Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
41,206 characters · 12 sections · 27 citation commands
Prediction Risk and Estimation Risk of the Ridgeless Least Squares Estimator under General Assumptions on Regression Errors
\onehalfspacing
Recent years have witnessed a fast growing body of work that analyzes minimum $\ell_2$ norm (ridgeless) interpolation least squares estimators bartlett2020benign,hastie2022surprises,tsigler2023benign. Researchers in this field were inspired by the ability of deep neural networks to accurately predict noisy training data with perfect fits, a phenomenon known as “double descent” or “benign overfitting” belkin2018understand,belkin2019reconciling,belkin2020two,zou2021benign,mei2022generalization. They discovered that to achieve this phenomenon, overparameterization is critical.
In the setting of linear regression, we have the training data $\{(x_i,y_i) \in\mathbb{R}^p\times \mathbb{R}: i=1,\cdots,n \}$, where the outcome variable $y_i$ is generated from
$x_i$ is a vector of features (or regressors), $\beta$ is a vector of unknown parameters, and $\varepsilon_i$ is a regression error. Here, $n$ is the sample size of the training data and $p$ is the dimension of the parameter vector $\beta$.
In the literature, the main object for the theoretical analyses has been mainly on the out-of-sample prediction risk. That is, for the ridge or interpolation estimator $\hat{\beta}$, the literature has focused on
where $x_0$ is a test observation that is identically distributed as $x_i$ but independent of the training data. For example, Dobriban:Wager,wu2020optimal,richards2021asymptotics,hastie2022surprises analyzed the predictive risk of ridge(less) regression and obtained exact asymptotic expressions under the assumption that $p/n$ converges to some constant as both $p$ and $n$ go to infinity. Overall, they found the double descent behavior of the ridgeless least squares estimator in terms of the prediction risk. bartlett2020benign,kobak2020optimal,tsigler2023benign characterized the phenomenon of benign overfitting in a different setting.
To the best of our knowledge, a vast majority of the theoretical analyses have been {confined to a simple data generating process, namely, the observations are independent and identically distributed (i.i.d.), and the regression errors have mean zero, have the common variance, and are independent of the feature vectors.} That is,
{This assumption, although convenient, is likely to be unrealistic in various real-world examples. For instance, liao2023economic adopted high-dimensional linear models to examine the double descent phenomenon in economic forecasts. In their applications, the outcome variables include S&P firms' earnings, U.S. equity premium, U.S. unemployment rate, and countries’ GDP growth rate. As in their applications, economic forecasts are associated with time series or panel data. As a result, it is improbable that ((ref)) holds in these applications. As another example, spiess2023double examined the performance of high-dimensional synthetic control estimators with many control units. The outcome variable in their application is the state-level smoking rates in the ADH:2010 dataset. Considering the geographical aspects of the U.S. states, it is unlikely that the regression errors underlying the synthetic control estimators adhere to ((ref)). In short, it is desirable to go beyond the simple but unrealistic regression error assumption given in ((ref)).}
To further motivate, we start with our own real-data example from American Community Survey (ACS) 2018, extracted from IPUMS USA IPUMS. The ACS is an ongoing annual survey by the US Census Bureau that provides key information about the US population. To have a relatively homogeneous population, the sample extract is restricted to white males residing in California with at least a bachelor’s degree. We consider a demographic group defined by their age, the type of degree, and the field of degree. Then, we compute the average of log hourly wages for each age-degree-field group, treat each group average as the outcome variable, and predict group wages by various group-level regression models where the regressors are constructed using the indicator variables of age, degree, and field as well as their interactions. We consider 7 specifications ranging from 209 to 2,182 regressors. To understand the role of non-i.i.d. regressor errors, we add the artificial noise to the training sample. See Appendix (ref) for details regarding how to generate the artificial noise. In the experiment, the constant $c$ varies such that $c=0$ corresponds to no clustered dependence across observations but as a positive $c$ gets larger, the noise has a larger share of clustered errors but the variance of the overall regression errors remains the same regardless of the value of $c$. Figure (ref) shows the in-sample (train) vs. out-of-sample (test) mean squared error (MSE) for various values of $c \in \{0, 0.25, 0.5, 0.75\}$. It can be seen that the experimental results are almost identical across different values of $c$ especially when $p > n$, suggesting that the double descent phenomenon might be universal for various degrees of clustered dependence, provided that the overall variance of the regression errors remains the same. It is our main goal to provide a firm foundation for this empirical phenomenon. To do so, we articulate the following research questions:
To delve into the prediction risk, suppose that $\Sigma :=\mathbb{E}[x_0 x_0^\top]$ is finite and positive definite. Then,
If $\Sigma = I$ (i.e., the case of isotropic features), where $I$ is the identity matrix, the mean squared error of the estimator defined by $\mathbb{E} [ \| \hat\beta-\beta \|^2 ]$, where $\| \cdot \|$ is the usual Euclidean norm, is the same as the expectation of the prediction risk defined above. However, if $\Sigma \neq I$, the link between the two quantities is less intimate. One may regard the prediction risk as the $\Sigma$-weighted mean squared error of the estimator; whereas $\mathbb{E} [ \| \hat\beta-\beta \|^2 ]$ can be viewed as an “unweighted” version, even if $\Sigma \neq I$. In other words, regardless of the variance-covariance structure of the feature vector, $\mathbb{E} [ \| \hat\beta-\beta \|^2 ]$ treats each component of $\beta$ “equally.” The mean squared error of the estimator is arguably one of the most standard criteria to evaluate the quality of the estimator in statistics. For instance, in the celebrated work by james:stein:61, the mean squared error criterion is used to show that the sample mean vector is not necessarily optimal even for standard normal vectors (so-called “Stein's paradox”). Many follow-up papers used the same criterion; e.g., hansen2016risk compared the mean-squared error of ordinary least squares, James–Stein, and Lasso estimators in an underparameterized regime. Both $\Sigma$-weighted and unweighted versions of the mean squared error are interesting objects to study. For example, Dobriban:Wager called the former “predictive risk” and the latter “estimation risk” in high-dimensional linear models; berthier2020tight called the former “generalization error” and the latter “reconstruction error” in the context of stochastic gradient descent for the least squares problem using the noiseless linear model. In this paper, we analyze both weighted and unweighted mean squared errors of the ridgeless estimator under general assumptions on the data-generating processes, not to mention anisotropic features. Furthermore, our focus is on the finite-sample analysis, that is, both $p$ and $n$ are fixed but $p > n$.
Although most of the existing papers consider the simple setting as in ((ref)), our work is not the first paper to consider more general regression errors in the overparameterized regime. Chinot2022AoS,chinot2020robustness analyzed minimum norm interpolation estimators as well as regularized empirical risk minimizers in linear models without any conditions on the regression errors. Specifically, chinot2020robustness showed that, with high probability, without assumption on the regression errors, for the minimum norm interpolation estimator, $(\hat{\beta}- \beta)^\top \Sigma (\hat{\beta}- \beta)$ is bounded from above by $\left( \|\beta\|^2 \sum_{i \geq c \cdot n} \lambda_i(\Sigma) \vee \sum_{i=1}^n \varepsilon_i^2 \right) / n$, where $c$ is an absolute constant and $\lambda_i (\Sigma)$ is the eigenvalues of $\Sigma$ in descending order. chinot2020robustness also obtained the bounds on the estimation error $(\hat{\beta}- \beta)^\top (\hat{\beta}- \beta)$. Our work is distinct and complements these papers in the sense that we allow for a general variance-covariance matrix of the regression errors. The main motivation of not making any assumptions on $\varepsilon_i$ in Chinot2022AoS and chinot2020robustness is to allow for potentially adversarial errors. We aim to allow for a general variance-covariance matrix of the regression errors to accommodate time series and clustered data, which are common in applications. See, e.g., hansen2022 for a textbook treatment (see Chapter 14 for time series and Section 4.21 for clustered data).
The main contribution of this paper is that we provide exact finite-sample characterization of the variance component of the prediction and estimation risks under the assumption that $X=[x_1,x_2,\cdots,x_n]^\top$ is left-spherical (e.g., $x_i$'s can be i.i.d. normal with mean zero but more general); $\varepsilon_i$'s can be correlated and have non-identical variances; and $\varepsilon_i$'s are independent of $x_i$'s. Specifically, the variance term can be factorized into a product between two terms: one term depends only on the trace of the variance-covariance matrix, say $\Omega$, of $\varepsilon_i$'s; the other term is solely determined by the distribution of $x_i$'s. Interestingly, we find that although $\Omega$ may contain non-zero off-diagonal elements, only the trace of $\Omega$ matters, as hinted by Figure (ref), and further demonstrate our finding via numerical experiments. In addition, we obtain exact finite-sample expression for the bias terms when the regression coefficients follow the random-effects hypothesis Dobriban:Wager. Our finite-sample findings offer a distinct viewpoint on the prediction and estimation risks, contrasting with the asymptotic inverse relationship (for optimally chosen ridge estimators) between the predictive and estimation risks uncovered by Dobriban:Wager. Finally, we connect our findings to the existing results on the prediction risk hastie2022surprises by considering the asymptotic behavior of estimation risk.
One of the limitations of our theoretical analysis is that the design matrix $X$ is assumed to be left-spherical, although it is more general than i.i.d. normal with mean zero. We not only view this as a convenient assumption but also expect that our findings will hold at least approximately even if $X$ does not follow the left-spherical distribution. It is a topic for future research to formally investigate this conjecture.
We first describe the minimum $\ell_2$ norm (ridgeless) interpolation least squares estimator in the overparameterized case ($p>n$). Define
so that $y = X \beta + \varepsilon$. The estimator we consider is
where $A^\dagger$ denotes the Moore–Penrose inverse of a matrix $A$.
The main object of interest in this paper is the prediction and estimation risks of $\hat \beta$ under the data scenario such that the regression error $\varepsilon_i$ may not be i.i.d. Formally, we make the following assumptions.
We emphasize that Assumption (ref) is more general than the standard assumption in the literature on benign overfitting that typically assumes that $\Omega \equiv \sigma^2 I$. Assumption (ref) allows for non-identical variances across the elements of $\varepsilon$ because the diagonal elements of $\Omega$ can be different among each other. Furthermore, it allows for non-zero off-diagonal elements in $\Omega$. It is difficult to assume that the regression errors are independent among each other with time series or clustered data; thus, in these settings, it is important to allow for general $\Omega \neq \sigma^2 I$. Below we present a couple of such examples.
For vector $a$ and square matrix $A$, let $\| a \|_A^2 := a^\top A a$. Conditional on $X$ and given $A$, we define
and we write $\mathrm{Var}=\mathrm{Var}_I$ and $\mathrm{Bias}=\mathrm{Bias}_I$ for the sake of brevity in notation.
The mean squared prediction error for an unseen test observation $x_0$ with the positive definite covariance matrix $\Sigma := \mathbb{E}[x_0 x_0^\top]$ (assuming that $x_0$ is independent of the training data $X$) and the mean squared estimation error of $\hat \beta$ conditional on $X$ can be written as:
In what follows, we obtain exact finite-sample expressions for prediction and estimation risks:
We first analyze the variance terms for both risks and then study the bias terms.
We rewrite the variance component of prediction risk as follows:
where positive definite symmetric matrices $S:=\Sigma^{1/2}$ and $T:=\Omega^{1/2}$ are the square root matrices of the positive definite matrices $\Sigma$ and $\Omega$, respectively. To compute the above Frobenius norm of the matrix $SX^{\dagger}T$, we need to compute the alignment of the right-singular vectors of $B:=SX^{\dagger}\in\mathbb{R}^{p\times n}$ with the left-eigenvectors of $T\in\mathbb{R}^{n\times n}$. Here, $B$ is a random matrix while $T$ is fixed. Therefore, we need the distribution of the right-singular vectors of the random matrix $B$.
Perhaps surprisingly, to compute the expected variance $\mathbb{E}_X[\mathrm{Var}_\Sigma(\hat\beta\mid X)]$, it turns out that we do not need the distribution of the singular vectors if we make a minimal assumption (the left-spherical symmetry of $X$) which is weaker than the assumption that $\{x_i\}_{i=1}^n$ is i.i.d. normal with $\mathbb{E}[x_1]=0$.
For the isotropic error case ($\Omega=I$), we have $\mathbb{E}_X[\mathrm{Var}_\Sigma(\hat\beta\mid X)]=\mathbb{E}_X[\operatorname{Tr}((X^\top X)^\dagger \Sigma)]$ directly from ((ref)) since $X^\dagger X^{\dagger\top}=(X^\top X)^\dagger$. Moreover, for the arbitrary error, the left-spherical symmetry of $X$ plays a critical role to factor out the same $\mathbb{E}_X[\operatorname{Tr}((X^\top X)^\dagger \Sigma)]$ and the trace of the variance-covariance matrix of the regression errors, $\operatorname{Tr}(\Omega)$, from the variance after the expectation over $X$.
The proofs of Lemma (ref) and Theorem (ref) are in the supplementary appendix.
For the expected variance $\mathbb{E}_X[\mathrm{Var}(\hat\beta\mid X)]$ of the estimation risk, a similar argument still holds if plugging-in $B=X^\dagger$ instead of $B=\Sigma^{1/2}X^\dagger$.
In this section, we validate our theory with some numerical experiments of Examples (ref) and (ref), especially how the expected variance is related to the general covariance $\Omega$ of the regressor error $\varepsilon$. In the both examples, we sample $\{x_i\}_{i=1}^n$ from $\mathcal{N}(0,\Sigma)$ with a general feature covariance $\Sigma=U_{\Sigma} D_{\Sigma} U_{\Sigma}^\top$ for an orthogonal matrix $U_{\Sigma} \in O(p)$ and a diagonal matrix $D_{\Sigma}\succ 0$. In this setting, we have $\mathrm{rank}(XX^\top)=n$ and $\Lambda^\dagger=\Lambda^{-1}$ almost everywhere.
As shown in Example (ref), when the regressor error follows an autoregressive process in ((ref)), we have $\Omega_{ij}=\sigma^2\rho^{|i-j|}/(1-\rho^2)$ and $\operatorname{Tr}(\Omega)/n=\sigma^2/(1-\rho^2)$. Therefore, for pairs of $(\sigma^2,\rho^2)$ with the same $\operatorname{Tr}(\Omega)/n$, they are expected to yield the same variances of the prediction and estimation risk from Theorem (ref) and (ref) even though they have different off-diagonal elements in $\Omega$. To be specific, the pairs $(\sigma^2,\rho^2)$ on a line $\{(\sigma^2,\rho^2): \sigma^2/\kappa^2+\rho^2 = 1\}$ have the same $\operatorname{Tr}(\Omega)/n$ and the same expected variance which gets larger for the line with respect to a larger $\kappa^2$.
Figure (ref) (left) shows the contour plots of $\mathbb{E}_X[\mathrm{Var}_\Sigma(\hat\beta\mid X)]$ and $\frac{1}{n}\operatorname{Tr}(\Omega)\mathbb{E}_X[\operatorname{Tr}((X^\top X)^\dagger \Sigma)]$ for different pairs of $(\sigma^2, \rho^2)$ in Example (ref). They have different slopes $-\kappa^{-2}$ according to the value of $\kappa^2=\operatorname{Tr}(\Omega)/n$. The right panel shows equivalent contour plots for estimation risk.
Now consider the block diagonal covariance matrix $\Omega=\operatorname{diag}(\Omega_1,\Omega_2,\cdots,\Omega_G)$ in Example (ref), where $\Omega_g$ is an $n_g\times n_g$ matrix with $(\Omega_g)_{ii} = \sigma_g^2$ and $(\Omega_g)_{ij} = \rho_g$ ($i\neq j$) for each $i, j=1,2,\cdots,n_g$ and $g=1,2,\cdots,G$. Let $n = \sum_{g=1}^G n_g$. We then have $\operatorname{Tr}(\Omega)/n=\sum_{g=1}^G \operatorname{Tr}(\Omega_g)/n=\sum_{g=1}^G(n_g/n) \sigma_g^2$. Therefore, given a partition $\{n_g\}_{g=1}^G$ of the $n$ observations, the covariance matrices $\Omega$ with different $\{\sigma^2_g\}_{g=1}^G$ have the same $\operatorname{Tr}(\Omega)/n$ if $(\sigma_1^2,\sigma_2^2,\cdots,\sigma_G^2)\in\mathbb{R}^G$ are on the same hyperplane $\frac{n_1}{n}\sigma_1^2+\frac{n_2}{n}\sigma_2^2+\cdots+\frac{n_G}{n}\sigma_G^2=\kappa^2$ for some $\kappa^2>0$.
Figure (ref) (left) shows the contour plots of $\mathbb{E}_X[\mathrm{Var}_\Sigma(\hat\beta\mid X)]$ and $\frac{1}{n}\operatorname{Tr}(\Omega)\mathbb{E}_X[\operatorname{Tr}((X^\top X)^\dagger \Sigma)]$ for different pairs of $(\sigma_1^2,\sigma_2^2)$ for a simple two-clusters example ($G=2$) of Example (ref) with $(n_1,n_2)=(5,15)$. Here, we use a fixed value of $\rho_1=\rho_2=0.05$, but the results are the same regardless of their values, as shown in the appendix. Unlike Example (ref), the hyperplanes are orthogonal to $v=[n_1,n_2]^\top$ regardless of the value of $\kappa^2=\operatorname{Tr}(\Omega)/n$. Again, the right panel shows equivalent contour plots for estimation risk.
Our main contribution is to allow for general assumptions on the regression errors, and thus the bias parts remain the same as they do not change with respect to the regression errors. For completeness, in this section, we briefly summarize the results on the bias components. First, we make the following assumption for a constant rank deficiency of $X^\top X$ which holds, for example, each $x_i$ has a positive definite covariance matrix and is independent of each other.
The bias term of prediction risk can be expressed as follows:
where $\hat\Sigma := X^\top X/n$. Now, in order to obtain an exact closed form solution, we make the following assumption:
A similar assumption (see Assumption (ref)) has been shown to be useful to obtain closed-form expressions in the literature Dobriban:Wager,richards2021asymptotics,li2021asymptotic,chen2023sketched.
Under this assumption, since $[\mathrm{Bias}_\Sigma(\hat\beta\mid X)]^2 = \operatorname{Tr}[S\beta(S\beta)^\top \lim_{\lambda\searrow 0} \lambda^2 (S^{-1}\hat\Sigma S+\lambda I)^{-2}]$ from ((ref)), we have the expected bias (conditional on $X$) as follows:
where $\tilde s_i$ are the eigenvalues of $S^{-1}\hat\Sigma S\in\mathbb{R}^{p\times p}$ and $\mathrm{rank}(S^{-1}\hat\Sigma S) =\mathrm{rank}(X)=n$ almost everywhere under Assumption (ref). This bias is independent of the distribution of $X$ or the spectral density of $S^{-1}\hat \Sigma S$, but only depending on the rank deficiency of the realization of $X$.
Finally, the prediction risk $R_P(\hat\beta)$ can be summarized as follows:
For the bias component of estimation risk, we can obtain a similar result with (ref) as follows:
Under Assumption (ref), we have the expected bias (conditional on $X$) as follows:
where $s_i$ are the eigenvalues of $\hat\Sigma\in\mathbb{R}^{p\times p}$ and $\mathrm{rank}(\hat\Sigma)=\mathrm{rank}(X)=n$ under Assumption (ref).
Thanks to Theorem (ref) and ((ref)), we obtain the following corollary for estimation risk.
The proof of Corollary (ref) is in the appendix.
To study the asymptotic behavior of estimation risk, we follow the previous approaches Dobriban:Wager,hastie2022surprises. First, we define the Stieltjes transform as follows:
We are now ready to investigate the asymptotic behavior of the mean squared estimation error with the following theorem:
This theorem is a direct consequence of Theorem 1.1 in silverstein1995empirical. Then, from Corollary (ref), we can write the limit of estimation risk as follows:
Here, the limit $s^\ast$ of the Stieltjes transform $s_F$ is highly connected with the shape of the spectral distribution of $\Sigma$. For example, in the case of isotropic features ($\Sigma=I$), i.e., $dH(\tau)=\delta(\tau-1) d\tau$, we have $s^\ast_\text{iso} =(\gamma-1)^{-1}$ from $1-\frac{1}{\gamma}=\frac{1}{1+s^\ast_\text{iso}}$. In addition, if $\Omega = \sigma^2 I$, then the limit of the mean squared error is exactly the same as the expression for $\gamma > 1$ in equation (10) of hastie2022surprises. This is because prediction risk is the same as estimation risk when $\Sigma = I$.
We conclude this paper by plotting the “descent curve” in the overparameterization regime in Figure (ref). On one hand, the expected variance ($\times$) perfectly matches its theoretical counterpart ($\medbullet$) and goes to zero as $\gamma$ gets large. On the other hand, the bias term is bounded even if $\gamma \rightarrow \infty$. The appendix contains the experimental details for all the figures.