Prediction Risk and Estimation Risk of the Ridgeless Least Squares Estimator under General Assumptions on Regression Errors

## Abstract

Abstract. In recent years, there has been a significant growth in research focusing on minimum ℓ 2 subscript ℓ 2 \ell_{2} roman_ℓ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT norm (ridgeless) interpolation least squares estimators. However, the majority of these analyses have been limited to an unrealistic regression error structure, assuming independent and identically distributed errors with zero mean and common variance. In this paper, we explore prediction risk as well as estimation risk under more general regression error assumptions, highlighting the benefits of overparameterization in a more realistic setting that allows for clustered or serial dependence. Notably, we establish that the estimation difficulties associated with the variance components of both risks can be summarized through the trace of the variance-covariance matrix of the regression errors. Our findings suggest that the benefits of overparameterization can extend to time series, panel and grouped data.

## 1. Introduction

Recent years have witnessed a fast growing body of work that analyzes minimum $\ell_{2}$ norm (ridgeless) interpolation least squares estimators (see, e.g., Bartlett
et al., ; Hastie et al., ; Tsigler and
Bartlett, , and references therein).
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”
(e.g., Belkin
et al., , , ; Zou
et al., ; Mei and
Montanari, , among many others).
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

$$ $\displaystyle y_{i}=x_{i}^{\top}\beta+\varepsilon_{i},\;i=1,\ldots,n,$ $$

$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

$$ $\displaystyle\mathbb{E}\Big{[}(x_{0}^{\top}\hat{\beta}-x_{0}^{\top}\beta)^{2} \mid x_{1},\ldots,x_{n}\Big{]},$ $$

where $x_{0}$ is a test observation that is identically distributed as $x_{i}$ but independent of the training data.
For example, Dobriban and
Wager (); Wu and Xu (); Richards
et al. (); Hastie et al. () 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.
Bartlett
et al. (); Kobak
et al. (); Tsigler and
Bartlett () 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,

$$ (1) $(y_{i},x_{i}^{\top})^{\top}\sim$ i.i.d. with $\mathbb{E}[\varepsilon_{i}]=0$ , $\mathbb{E}[\varepsilon_{i}^{2}]=\sigma^{2}<\infty$ and $\varepsilon_{i}$ is independent of $x_{i}$ . $$

This assumption, although convenient, is likely to be unrealistic in various real-world examples.
For instance, Liao
et al. () 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 (1) holds in these applications. As another example, Spiess
et al. () 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 Abadie
et al. () dataset. Considering the geographical aspects of the U.S. states, it is unlikely that the regression errors underlying the synthetic control estimators adhere to (1). In short, it is desirable to go beyond the simple but unrealistic regression error assumption given in (1).

Figure: Figure 1. Comparison of in-sample and out-of-sample mean squared error (MSE) across various degrees of clustered noise. The vertical line indicates $p=n\;(=1,415)$.
Refer to caption: x1.png

To further motivate, we start with our own real-data example from American Community Survey (ACS) 2018, extracted from IPUMS USA (Ruggles et al., ). 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 A 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 1 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:

- •
How to analyze the out-of-sample prediction risk of the ridgeless least squares estimator under *general* assumptions on the regression errors?
- •
Why does not the prediction risk seem to be affected by the degrees of dependence across observations?

To delve into the prediction risk, suppose that $\Sigma:=\mathbb{E}[x_{0}x_{0}^{\top}]$ is finite and positive definite. Then,

$$ $\displaystyle\mathbb{E}\Big{[}(x_{0}^{\top}\hat{\beta}-x_{0}^{\top}\beta)^{2} \mid x_{1},\ldots,x_{n}\Big{]}=\mathbb{E}\Big{[}(\hat{\beta}-\beta)^{\top} \Sigma(\hat{\beta}-\beta)\mid x_{1},\ldots,x_{n}\Big{]}.$ $$

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 and
Stein (), 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.,
Hansen () 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 and
Wager () called the former “predictive risk” and the latter “estimation risk” in high-dimensional linear models;
Berthier
et al. () 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 (1), our work is not the first paper to consider more general regression errors in the overparameterized regime.
Chinot
et al. (); Chinot and
Lerasle () analyzed minimum norm interpolation estimators as well as regularized empirical risk minimizers in linear models without any conditions on the regression errors. Specifically,
Chinot and
Lerasle () 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.
Chinot and
Lerasle () 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 Chinot
et al. () and Chinot and
Lerasle () 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., Hansen () 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 1, 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 and
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 and
Wager ().
Finally, we connect our findings to the existing results on the prediction risk (e.g., Hastie et al., ) 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.

## 2. The Framework under General Assumptions on Regression Errors

We first describe the minimum $\ell_{2}$ norm (ridgeless) interpolation least squares estimator in the overparameterized case ($p>n$).
Define

$$ $\displaystyle y$ $\displaystyle:=[y_{1},y_{2},\cdots,y_{n}]^{\top}\in\mathbb{R}^{n},\;\$ $\displaystyle\varepsilon$ $\displaystyle:=[\varepsilon_{1},\varepsilon_{2},\cdots,\varepsilon_{n}]^{\top} \in\mathbb{R}^{n},\;\$ $\displaystyle X^{\top}$ $\displaystyle:=\begin{bmatrix}x_{1},x_{2},\cdots,x_{n}\end{bmatrix}\in\mathbb{ R}^{p\times n},$ $$

so that $y=X\beta+\varepsilon$.
The estimator we consider is

$$ $\displaystyle\hat{\beta}$ $\displaystyle:=\operatorname*{arg\,min}_{b\in\mathbb{R}^{p}}\{\|b\|:Xb=y\}=(X^ {\top}X)^{\dagger}X^{\top}y=X^{\dagger}y,$ $$

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.

###### Assumption 2.1 .

(i) $y=X\beta+\varepsilon$, where
$\varepsilon$ is independent of $X$,
and
$\mathbb{E}[\varepsilon]=0$.
(ii) $\Omega:=\mathbb{E}[\varepsilon\varepsilon^{\top}]$ is finite and positive definite (but not necessarily spherical).

We emphasize that Assumption 2.1 is more general than the standard assumption in the literature on benign overfitting that typically assumes that $\Omega\equiv\sigma^{2}I$.
Assumption 2.1 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.

###### Example 2.1 (AR(1) Errors) .

Suppose that the regressor error follows an autoregressive process:

$$ (2) $\displaystyle\varepsilon_{i}=\rho\varepsilon_{i-1}+\eta_{i},$ $$

where $\rho\in(-1,1)$ is an autoregressive parameter,
$\eta_{i}$ is independent and identically distributed with mean zero and variance $\sigma^{2}(0<\sigma^{2}<\infty)$ and is independent of $X$. Then, the $(i,j)$ element of $\Omega$ is

$$ $\displaystyle\Omega_{ij}=\frac{\sigma^{2}}{1-\rho^{2}}\rho^{|i-j|}.$ $$

Note that $\Omega_{ij}\neq 0$ as long as $\rho\neq 0$.

###### Example 2.2 (Clustered Errors) .

Suppose that regression errors are mutually independent across clusters but they can be arbitrarily correlated within the same cluster. For instance, students in the same school may affect each other and also have the same teachers; thus it would be difficult to assume independence across student test scores within the same school. However, it might be reasonable that student test scores are independent across different schools.
For example, assume that (i) if the regression error
$\varepsilon_{i}$ belongs to cluster $g$,
where $g=1,\ldots,G$ and $G$ is the number of clusters,
$\mathbb{E}[\varepsilon_{i}^{2}]=\sigma_{g}^{2}$ for some constant $\sigma_{g}^{2}>0$ that can vary over $g$;
(ii) if the regression errors
$\varepsilon_{i}$ and $\varepsilon_{j}$ $(i\neq j)$ belong to the same cluster $g$,
$\mathbb{E}[\varepsilon_{i}\varepsilon_{j}]=\rho_{g}$ for some constant $\rho_{g}\neq 0$ that can be different across $g$;
and
(iii) if the regression errors
$\varepsilon_{i}$ and $\varepsilon_{j}$ $(i\neq j)$ do not belong to the same cluster,
$\mathbb{E}[\varepsilon_{i}\varepsilon_{j}]=0$.
Then, $\Omega$ is block diagonal with possibly non-identical blocks.

For vector $a$ and square matrix $A$, let $\|a\|_{A}^{2}:=a^{\top}Aa$. Conditional on $X$ and given $A$,
we define

$$ $\displaystyle\mathrm{Bias}_{A}(\hat{\beta}\mid X):=\|\mathbb{E}[\hat{\beta} \mid X]-\beta\|_{A}\;\;\text{ and }\;\;\mathrm{Var}_{A}(\hat{\beta}\mid X):= \operatorname{Tr}(\mathrm{Cov}(\hat{\beta}\mid X)A),$ $$

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:

$$ $\displaystyle R_{P}(\hat{\beta}\mid X)$ $\displaystyle:=\mathbb{E}\big{[}(x_{0}^{\top}\hat{\beta}-x_{0}^{\top}\beta)^{2 }\mid X\big{]}=[\mathrm{Bias}_{\Sigma}(\hat{\beta}\mid X)]^{2}+\mathrm{Var}_{ \Sigma}(\hat{\beta}\mid X),$ $\displaystyle R_{E}(\hat{\beta}\mid X)$ $\displaystyle:=\mathbb{E}\big{[}\|\hat{\beta}-\beta\|^{2}\mid X\big{]}=[ \mathrm{Bias}(\hat{\beta}\mid X)]^{2}+\mathrm{Var}(\hat{\beta}\mid X).$ $$

In what follows, we obtain exact finite-sample expressions for prediction and estimation risks:

$$ $\displaystyle R_{P}(\hat{\beta}):=\mathbb{E}_{X}[R_{P}(\hat{\beta}\mid X)]\;\; \text{ and }\;\;R_{E}(\hat{\beta}):=\mathbb{E}_{X}[R_{E}(\hat{\beta}\mid X)].$ $$

We first analyze the variance terms for both risks and then study the bias terms.

###### Assumption 2.1 .

###### Example 2.1 (AR(1) Errors) .

###### Example 2.2 (Clustered Errors) .

## 3. The Variance Components of Prediction and Estimation Risks

### 3.1. The variance component of prediction risk

We rewrite the variance component of prediction risk as follows:

$$ (3) $\displaystyle\mathrm{Var}_{\Sigma}(\hat{\beta}\mid X)=\operatorname{Tr}( \mathrm{Cov}(\hat{\beta}\mid X)\Sigma)$ $\displaystyle=\operatorname{Tr}(X^{\dagger}\Omega X^{\dagger\top}\Sigma)=\|SX^ {\dagger}T\|_{F}^{2},$ $$

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$.

###### Definition 3.1 (Left-Spherical Symmetry (Dawid, 1977 , 1978 , 1981 ; Gupta and Nagar, 1999 ) ) .

A random matrix $Z$ or its distribution is called to be left-spherical if $OZ$ and $Z$ have the same distribution ($OZ\overset{d}{=}Z$) for any fixed orthogonal matrix $O\in O(n):=\{A\in\mathbb{R}^{n\times n}:AA^{\top}=A^{\top}A=I\}$.

###### Assumption 3.1 .

The design matrix $X$ is left-spherical.

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 (3) 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$.

###### Lemma 3.1 .

For a subset ${\mathcal{S}}\subset\mathbb{R}^{m\times m}$ satisfying $C^{-1}\in{\mathcal{S}}$ for all $C\in{\mathcal{S}}$, if matrix-valued random variables $Z$ and $AZ$ have the same distribution measure $\mu_{Z}$ for any $A\in{\mathcal{S}}$,
then we have

$$ $\mathbb{E}_{Z}[f(Z)]=\mathbb{E}_{Z}[f(AZ)]=\mathbb{E}_{Z}[\mathbb{E}_{A^{ \prime}\sim\nu}[f(A^{\prime}Z)]]$ $$

for any function $f\in L^{1}(\mu_{Z})$ and any probability density function $\nu$
on ${\mathcal{S}}$.

###### Theorem 3.2 .

Let Assumptions
2.1, and  3.1 hold.
Then,
we have

$$ $\displaystyle\mathbb{E}_{X}[\mathrm{Var}_{\Sigma}(\hat{\beta}\mid X)]$ $\displaystyle=\frac{1}{n}\operatorname{Tr}(\Omega)\mathbb{E}_{X}[\operatorname {Tr}((X^{\top}X)^{\dagger}\Sigma)].$ $$

###### Sketch of Proof.

With $B=\Sigma^{1/2}X^{\dagger}$ and $T=\Omega^{1/2}$, we can rewrite the variance as follows:

$$ $\displaystyle\mathrm{Var}_{\Sigma}(\hat{\beta}\mid X)=\|BT\|_{F}^{2}=\|UDV^{ \top}U_{T}D_{T}V_{T}^{\top}\|_{F}^{2}{=\|DV^{\top}U_{T}D_{T}\|_{F}^{2}}$ $$

from the singular value decompositions $B=UDV^{\top}$ and $T=U_{T}D_{T}V_{T}^{\top}$ with orthogonal matrices $U,V,U_{T},V_{T}$, and diagonal matrices $D,D_{T}$.
Then, we need to compute the alignment $V^{\top}U_{T}$ of the right-singular vectors of $B$ with the left-eigenvectors of $T$ because

$$ $\displaystyle{\|DV^{\top}U_{T}D_{T}\|_{F}^{2}}$ $\displaystyle=\lambda\left((X^{\top}X)^{\dagger}\Sigma\right)^{\top}\Gamma(X) \lambda(\Omega)=a(X)^{\top}\Gamma(X)b,$ $$

where $v^{(i)}:=V_{:i}$, $u^{(j)}:=(U_{T})_{:j}$, $\gamma_{ij}:=\langle v^{(i)},u^{(j)}\rangle^{2}\geq 0$, $\Gamma(X):=(\gamma_{ij})_{i,j}\in\mathbb{R}^{n\times n}$
and $\lambda(A)\in\mathbb{R}^{n}$ is a vector where its elements are the eigenvalues of $A$.

Now, we want to compute the expected variance. To do so, from Lemma 3.1 with ${\mathcal{S}}=O(n)$ and the left-spherical symmetry of $X$, we can obtain

$$ $\displaystyle\mathbb{E}_{X}[a(X)^{\top}\Gamma(X)b]$ $\displaystyle=\mathbb{E}_{X}\left[\mathbb{E}_{O\sim\nu}[a(OX)^{\top}\Gamma(OX) b]\right]=\mathbb{E}_{X}\left[a(X)^{\top}\mathbb{E}_{O\sim\nu}[\Gamma(OX)]b \right],$ $$

where $\nu$ is the unique uniform distribution (the Haar measure) over the orthogonal matrices $O(n)$.

Here, we can show that $\mathbb{E}_{O\sim\nu}[\Gamma(OX)]=\frac{1}{n}J$,
where $J$ is the all-ones matrix with
$J_{ij}=1(i,j=1,2,\cdots,n)$.
Therefore, we have the expected variance as follows:

$$ $\displaystyle\mathbb{E}_{X}[\mathrm{Var}_{\Sigma}(\hat{\beta}\mid X)]$ $\displaystyle=\mathbb{E}_{X}\left[a(X)^{\top}\frac{1}{n}Jb\right]=\frac{1}{n} \sum_{i,j=1}^{n}\mathbb{E}_{X}[a_{i}(X)]b_{j}=\frac{1}{n}\mathbb{E}_{X}[ \operatorname{Tr}((X^{\top}X)^{\dagger}\Sigma)]\operatorname{Tr}(\Omega).$ $$

∎

The proofs of Lemma 3.1 and Theorem 3.2 are in the supplementary appendix.

###### Definition 3.1 (Left-Spherical Symmetry (Dawid, 1977 , 1978 , 1981 ; Gupta and
Nagar, 1999 ) ) .

###### Assumption 3.1 .

###### Lemma 3.1 .

###### Theorem 3.2 .

###### Sketch of Proof.

### 3.2. The variance component of estimation risk

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}$.

###### Theorem 3.3 .

Let Assumptions
2.1, and 3.1 hold.
Then,
we have

$$ $\displaystyle\mathbb{E}_{X}[\mathrm{Var}(\hat{\beta}\mid X)]$ $\displaystyle=\frac{1}{np}\operatorname{Tr}(\Omega)\mathbb{E}_{X}[ \operatorname{Tr}(\Lambda^{\dagger})],$ $$

where
$XX^{\top}/p=U\Lambda U^{\top}$ for some orthogonal matrix $U\in O(n)$.

Figure: Figure 2. Our theory (dashed lines) matches the expected variances (solid lines) of the prediction (left) and estimation risks (right) in Example 2.1 (AR(1) Errors). Each point $(\sigma^{2},\rho^{2})$ represents a different noise covariance matrix $\Omega$, but with the same $\operatorname{Tr}(\Omega)$ along each line $\{(\sigma^{2},\rho^{2}):\sigma^{2}/\kappa^{2}+\rho^{2}=1\}$ for some $\kappa^{2}>0$, they have the same expected variance. We set $n=50,p=100$, and evaluate on 100 samples of $X$ and 100 samples of $\varepsilon$ (for each realization of $X$) to approximate the expectations.
Refer to caption: x2.png

###### Theorem 3.3 .

### 3.3. Numerical experiments

In this section, we validate our theory with some numerical experiments of Examples 2.1 and 2.2, 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.

#### 3.3.1. AR(1) Errors

As shown in Example 2.1,
when the regressor error follows an autoregressive process in (2), 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 3.2 and 3.3 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 2 (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 2.1.
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.

#### 3.3.2. Clustered Errors

Now consider the block diagonal covariance matrix $\Omega=\operatorname{diag}(\Omega_{1},\Omega_{2},\cdots,\Omega_{G})$ in Example 2.2,
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: Figure 3. Our theory (dashed lines) matches the expected variances (solid lines) of the prediction (left) and estimation risks (right) in Example 2.2 (Clustered Errors). Each point $(\sigma^{2},\rho^{2})$ represents a different noise covariance matrix $\Omega$, but with the same $\operatorname{Tr}(\Omega)$ along each line $\{(\sigma_{1}^{2},\sigma_{2}^{2}):\frac{n_{1}}{n}\sigma_{1}^{2}+\frac{n_{2}}{n }\sigma_{2}^{2}=\kappa^{2}\}$ for some $\kappa^{2}>0$, they have the same expected variance. We set $G=2,(n_{1}=5,n_{2}=15),n=20,p=40,\rho_{1}=\rho_{2}=0.05$, and evaluate on 100 samples of $X$ and 100 samples of $\varepsilon$ (for each realization of $X$) to approximate the expectations.
Refer to caption: x4.png

Figure 3 (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 2.2 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 2.1, 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.

## 4. The Bias Components of Prediction and Estimation Risks

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.

###### Assumption 4.1 .

$\mathrm{rank}(X)=n$ almost everywhere.

###### Assumption 4.1 .

### 4.1. The bias component of prediction risk

The bias term of prediction risk can be expressed as follows:

$$ (4) $\displaystyle[\mathrm{Bias}_{\Sigma}(\hat{\beta}\mid X)]^{2}$ $\displaystyle=(S\beta)^{\top}\lim_{\lambda\searrow 0}\lambda^{2}(S^{-1}\hat{ \Sigma}S+\lambda I)^{-2}S\beta,$ $$

where $\hat{\Sigma}:=X^{\top}X/n$.
Now, in order to obtain an exact closed form solution, we make the following assumption:

###### Assumption 4.2 .

$\mathbb{E}_{\beta}[S\beta(S\beta)^{\top}]=r_{\Sigma}^{2}I/p$,
where $r_{\Sigma}^{2}:=\mathbb{E}_{\beta}[\|\beta\|_{\Sigma}^{2}]<\infty$ and $\beta$ is independent of $X$.

A similar assumption (see Assumption 4.3) has been shown to be useful to obtain closed-form expressions in the literature (e.g., Dobriban and
Wager, ; Richards
et al., ; Li
et al., ; Chen
et al., ).

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 (4), we have the expected bias (conditional on $X$) as follows:

$$ $\displaystyle\mathbb{E}_{\beta}[\mathrm{Bias}_{\Sigma}(\hat{\beta}\mid X)^{2} \mid X]$ $\displaystyle=\frac{r_{\Sigma}^{2}}{p}\lim_{\lambda\searrow 0}\sum_{i=1}^{p} \frac{\lambda^{2}}{(\tilde{s}_{i}+\lambda)^{2}}=\frac{r_{\Sigma}^{2}}{p}\left| \{i\in[p]:\tilde{s}_{i}=0\}\right|=r_{\Sigma}^{2}\frac{p-n}{p},$ $$

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 4.1.
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:

###### Corollary 4.1 .

Let Assumptions
2.1, 3.1,
4.1,
and  4.2
hold. Then, we have

$$ $\displaystyle R_{P}(\hat{\beta})$ $\displaystyle=r_{\Sigma}^{2}\left(1-\frac{n}{p}\right)+\frac{\operatorname{Tr} (\Omega)}{n}\mathbb{E}_{X}\left[\operatorname{Tr}((X^{\top}X)^{\dagger}\Sigma) \right].$ $$

###### Assumption 4.2 .

###### Corollary 4.1 .

### 4.2. The bias component of estimation risk

For the bias component of estimation risk, we can obtain a similar result with 4.1 as follows:

$$ $\displaystyle[\mathrm{Bias}(\hat{\beta}\mid X)]^{2}=\beta^{\top}(I-\hat{\Sigma }^{\dagger}\hat{\Sigma})\beta$ $\displaystyle=\lim_{\lambda\searrow 0}\beta^{\top}\lambda(\hat{\Sigma}+\lambda I )^{-1}\beta.$ $$

###### Assumption 4.3 .

$\mathbb{E}_{\beta}[\beta\beta^{\top}]=r^{2}I/p$,
where $r^{2}:=\mathbb{E}_{\beta}[\|\beta\|^{2}]<\infty$ and $\beta$ is independent of $X$.

Under Assumption 4.3,
we have the expected bias (conditional on $X$) as follows:

$$ (5) $\displaystyle\mathbb{E}_{\beta}[\mathrm{Bias}(\hat{\beta}\mid X)^{2}\mid X]$ $\displaystyle=\frac{r^{2}}{p}\lim_{\lambda\searrow 0}\sum_{i=1}^{p}\frac{ \lambda}{s_{i}+\lambda}=\frac{r^{2}}{p}\left|\{i\in[p]:s_{i}=0\}\right|=r^{2} \frac{p-n}{p},$ $$

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 4.1.

Thanks to Theorem 3.3
and (5),
we obtain the following corollary for estimation risk.

###### Corollary 4.2 .

Let Assumptions
2.1,
3.1,
4.1,
and  4.3
hold. Then, we have

$$ $\displaystyle R_{E}(\hat{\beta})$ $\displaystyle=r^{2}\left(1-\frac{n}{p}\right)+\frac{\operatorname{Tr}(\Omega)} {n}\mathbb{E}_{X}\left[\int\frac{1}{s}dF^{XX^{\top}/n}(s)\right],$ $$

where $F^{A}(s):=\frac{1}{n}\sum_{i=1}^{n}1\{\lambda_{i}(A)\leq s\}$ is the empirical spectral distribution of a matrix $A$ and $\lambda_{1}(A),\lambda_{2}(A),\cdots,\lambda_{n}(A)$
are the eigenvalues of $A$.

The proof of Corollary 4.2 is in the appendix.

###### Assumption 4.3 .

###### Corollary 4.2 .

#### 4.2.1. Asymptotic analysis of estimation risk

To study the asymptotic behavior of estimation risk, we follow the previous approaches (Dobriban and
Wager, ; Hastie et al., ).
First, we define the Stieltjes transform as follows:

###### Definition 4.1 .

The Stieltjes transform $s_{F}(z)$ of a df $F$ is defined as:

$$ $s_{F}(z):=\int\frac{1}{x-z}dF(x),\ \text{for}\ z\in\mathbb{C}\setminus\mathrm{ supp}(F).$ $$

Figure: Figure 4. The “descent curve” in the overparameterization regime for prediction risk (left) and estimation risk (right). We test $\Omega$’s with $\operatorname{Tr}(\Omega)/n=1,2,4$ in black, blue, red, respectively. For the anisotropic feature, the expected variance ($\times$) and its theoretical expression ($\medbullet$) are $\Theta\left(\frac{\operatorname{Tr}(\Omega)/n}{\gamma-1}\right)$ and larger than that in the high-dimensional asymptotics for the isotropic $\Sigma=I$. For the isotropic $\Sigma=I$, the variance terms (dotted) and the bias terms (dashed) in the high-dimensional asymptotics are $\frac{1}{\gamma-1}\lim_{n\rightarrow\infty}\frac{\operatorname{Tr}(\Omega)}{n}$ and $r^{2}\left(1-\frac{1}{\gamma}\right)$, respectively.
Refer to caption: x6.png

We are now ready to investigate the asymptotic behavior of the mean squared estimation error
with the following theorem:

###### Theorem 4.3 .

(Silverstein and
Bai, , Theorem 1.1)
Suppose that the rows $\{x_{i}\}_{i=1}^{n}$ in $X$ are i.i.d. centered random vectors with $\mathbb{E}[x_{1}x_{1}^{\top}]=\Sigma$ and that the empirical spectral distribution
$F^{\Sigma}(s)=\frac{1}{p}\sum_{i=1}^{p}1\{\tau_{i}\leq s\}$ of $\Sigma$
converges almost surely to a probability distribution function $H$ as $p\rightarrow\infty$.
When $p/n\rightarrow\gamma>0$ as $n,p\rightarrow\infty$, then a.s., $F^{XX^{\top}/n}$ converges vaguely to a df $F$ and
the limit
$s^{\ast}:=\lim_{z\searrow 0}s_{F}(z)$
of its Stieltjes transform $s_{F}$
is the unique solution to the equation:

$$ (6) $\displaystyle 1-\frac{1}{\gamma}=\int\frac{1}{1+\tau s^{\ast}}dH(\tau).$ $$

This theorem is a direct consequence of Theorem 1.1 in Silverstein and
Bai ().
Then, from Corollary 4.2, we can write the limit of estimation risk as follows:

###### Corollary 4.4 .

Let Assumptions 2.1,
3.1,
4.1,
and  4.3
hold. Then, under the same assumption as Theorem 4.3, as $n,p\rightarrow\infty$ and $p/n\rightarrow\gamma$, where $1<\gamma<\infty$ is a constant,
we have

$$ $\displaystyle R_{E}(\hat{\beta})$ $\displaystyle=\mathbb{E}\big{[}\|\hat{\beta}-\beta\|^{2}\big{]}\rightarrow r^{ 2}\left(1-\frac{1}{\gamma}\right)+s^{\ast}\lim_{n\rightarrow\infty}\frac{ \operatorname{Tr}(\Omega)}{n}.$ $$

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 Hastie et al. (, Theorem 1). This is because prediction risk is the same as estimation risk when $\Sigma=I$.

###### Remark 4.1 .

Generally, if the support of $H$ is bounded within $[c_{H},C_{H}]\subset\mathbb{R}$ for some positive constants $0<c_{H}\leq C_{H}<\infty$, then we can observe the double descent phenomenon in the overparameterization regime with
$\lim_{\gamma\searrow 1}s^{\ast}=\infty$ and $\lim_{\gamma\rightarrow\infty}s^{\ast}=0$ with $s^{\ast}=\Theta\left(\frac{1}{\gamma-1}\right)$ from the following inequalities:

$$ (7) $\displaystyle C_{H}^{-1}\frac{1}{\gamma-1}\leq s^{\ast}\leq c_{H}^{-1}\frac{1} {\gamma-1}.$ $$

In fact, a tighter lower bound is available:

$$ (8) $\displaystyle s^{\ast}\geq\mu_{H}^{-1}(\gamma-1)^{-1},$ $$

where $\mu_{H}:=\mathbb{E}_{\tau\sim H}[\tau]$, i.e., the mean of distribution $H$.
The proofs of (7) and (8) are given in the supplementary appendix.

We conclude this paper by plotting the “descent curve” in the overparameterization regime in Figure 4.
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.

###### Definition 4.1 .

###### Theorem 4.3 .

###### Corollary 4.4 .

###### Remark 4.1 .

## Appendix A Details for drawing Figure 1

To draw Figure 1, we use a sample extract from American Community Survey (ACS) 2018. 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 in years (between 25 and 70), the type of degree (bachelor’s, master’s, professional, and doctoral), and the field of degree (172 unique values).
Then, we compute the average of log hourly wages for each age-degree-field group (all together 7,073 unique groups in the sample). We treat each group average as the outcome variable (say, $y_{a,d,f}$) 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: that is,

$$ $\displaystyle y_{a,d,f}=x_{a,d,f}^{\top}\,\beta+\varepsilon_{a,d,f}.$ $$

For the regressors $x_{a,d,f}$, we consider 7 specifications ranging from 209 to 2,183 regressors:

- •
Spec. 1 $(p=209)$: dummy variables for age (say, $x_{a}$) $+$ dummy variables for the type of degree (say, $x_{d}$) $+$ dummy variables for the field of degree (say, $x_{f}$),
- •
Spec. 2 $(p=391)$: Spec. 1 $+$ all interactions between $x_{d}$ and $x_{a}$,
- •
Spec. 3 $(p=598)$: Spec. 1 $+$ all interactions between $x_{d}$ and $x_{f}$,
- •
Spec. 4 $(p=778)$: Spec. 1 $+$ all interactions between $x_{d}$ and $x_{a}$ $+$ all interactions between $x_{d}$ and $x_{f}$,
- •
Spec. 5 $(p=1640)$: Spec. 1 $+$ all interactions between $x_{d}$ and $x_{a}$ $+$ all interactions between $x_{a}$ and $x_{f}$,
- •
Spec. 6 $(p=1754)$: Spec. 1 $+$ all interactions between $x_{d}$ and $x_{f}$ $+$ all interactions between $x_{a}$ and $x_{f}$,
- •
Spec. 7 $(p=2182)$: Spec. 1 $+$ all three-way interactions among $x_{a}$, $x_{d}$ and $x_{f}$.

Here, the dummy variable are constructed using one-hot encoding.
We randomly split the sample into the train and test samples with a ratio of $1:4$. The resulting sample sizes are 1,415 and 5,658, respectively. To understand the role of non-i.i.d. regressor errors, we add the artificial noise to the training sample: that is, we compute the ridgeless least squares estimator using the training sample of $(\tilde{y}_{a,d,f},x_{a,d,f}^{\top})^{\top}$, where $\tilde{y}_{a,d,f}=y_{a,d,f}+u_{a,d,f}$. Here,
the artificial noise $u_{a,d,f}$ has the form

$$ $\displaystyle u_{a,d,f}\equiv\frac{(1-c)e_{a,d,f}+c\cdot e_{f}}{\sqrt{(1-c)^{2 }+c^{2}}},$ $$

where $e_{a,d,f}\sim N(0,\sigma^{2})$, independently across age ($a$), degree ($d$) and field ($f$); $e_{f}$ is the average of another independent $N(0,\sigma^{2})$ variable within $f$ (hence, $e_{f}$ is identical for each value of $f$) and thus the source of clustered errors; and $c\in\{0,0.25,0.5,0.75\}$ is a constant that will be varied across the experiment. As $c$ gets larger, the noise has a larger share of clustered errors but the variance of the overall regression errors ($u_{a,d,f}$) remains the same: in other words, $\text{var}(u_{a,d,f})=\sigma^{2}$ for each value of $c$.
Figure 1 was generated with $\sigma=0.5$ by generating the artificial noise only once.

## Appendix B Details for drawing Figures 2 , 3 , and 4

To draw Figure 2, 3, and 4, we sample $\{x_{i}\}_{i=1}^{n}$ from $\mathcal{N}(0,\Sigma)$ with $\Sigma=U_{\Sigma}D_{\Sigma}U_{\Sigma}^{\top}$
where $U_{\Sigma}$ is an orthogonal matrix random variable, drawn from the uniform (Haar) distribution on $O(p)$, and $D_{\Sigma}$ is a diagonal matrix with its elements $d_{i}=|z_{i}|/\sum_{i=1}^{p}|z_{i}|$ being sampled with $z_{i}\sim\mathcal{N}(0,1)$ for each $i=1,2,\cdots,p$.
With this general anisotropic $\Sigma$,
the term $\mathbb{E}_{X}[\operatorname{Tr}(\Lambda^{-1})]/p$
is somewhat larger than $\mu_{H}^{-1}s^{\ast}_{\text{iso}}=(\gamma-1)^{-1}$ which is $1$ in Figure 2 and 3 since $\mu_{H}=1$ and $\gamma=2$.
For example, in Figure 2, when $\sigma^{2}=1,\rho^{2}=0$, we have $\operatorname{Tr}(\Omega)/n=1$ but $\operatorname{Tr}(\Omega)\mathbb{E}_{X}[\operatorname{Tr}(\Lambda^{-1})]/(np)>1$.

In Figure 4,
we fix $n=50$ and use $p=n\gamma$ for $\gamma\in[1,100]$.

To compute the expectations of $\mathbb{E}_{X}[\mathrm{Var}(\hat{\beta}|X)]$ and $\mathbb{E}_{X}[\operatorname{Tr}(\Lambda^{-1})]$ over $X$, we sample $N_{X}$ samples of $X$’s, $X_{1},X_{2},\cdots,X_{N_{X}}$.
Moreover, to compute the expectation over $\varepsilon$ in

$$ $\mathrm{Var}(\hat{\beta}|X_{i})\equiv\operatorname{Tr}\left(\mathbb{E}_{ \varepsilon}[\hat{\beta}\hat{\beta}^{\top}]-\mathbb{E}_{\varepsilon}[\hat{ \beta}]\mathbb{E}_{\varepsilon}[\hat{\beta}]^{\top}\right),$ $$

we sample $N_{\varepsilon}$ samples of $\varepsilon$’s, $\varepsilon_{1},\varepsilon_{2},\cdots,\varepsilon_{N_{\varepsilon}}$ for each realization $X_{i}$.
To be specific,

$$ $\displaystyle\mathbb{E}_{X}[\mathrm{Var}(\hat{\beta}|X)]$ $\displaystyle\approx\frac{1}{N_{X}}\sum_{i=1}^{N_{X}}\mathrm{Var}(\hat{\beta}| X_{i})\approx\frac{1}{N_{X}}\sum_{i=1}^{N_{X}}\operatorname{Tr}\left(\frac{1}{ N_{\varepsilon}}\sum_{j=1}^{N_{\varepsilon}}\hat{\beta}_{i,j}\hat{\beta}_{i,j} ^{\top}-\frac{1}{N_{\varepsilon}}\sum_{j=1}^{N_{\varepsilon}}\hat{\beta}_{i,j} \frac{1}{N_{\varepsilon}}\sum_{j=1}^{N_{\varepsilon}}\hat{\beta}_{i,j}^{\top}\right)$ $\displaystyle\frac{1}{p}\mathbb{E}_{X}[\operatorname{Tr}(\Lambda^{-1})]$ $\displaystyle\approx\frac{1}{N_{X}}\sum_{i=1}^{N_{X}}\operatorname{Tr}((X_{i}X _{i}^{\top})^{-1})=\frac{1}{N_{X}}\sum_{i=1}^{N_{X}}\sum_{k=1}^{n}\frac{1}{ \lambda_{k}(X_{i}X_{i}^{\top})},$ $$

where $\hat{\beta}_{i,j}=\operatorname*{arg\,min}_{\beta}\{\|b\|:X_{i}b-y_{i,j}=0\}$, $y_{i,j}=X_{i}\beta+\varepsilon_{j}$,
and $\lambda_{k}(X_{i}X_{i}^{\top})$ is the $k$-th eigenvalue of $X_{i}X_{i}^{\top}$.
We can do similarly for the variance part of the prediction risk.

Figure 5 shows an additional experimental result.

Figure: Figure 5. We use the same setting as Figure 3, except uniformly sample each $\rho_{i}$ from $[0,0.05]$ for each experiment with the pairs $(\sigma_{1}^{2},\sigma_{1}^{2})$. As expected, the off-diagonal elements $\rho_{i}$ of $\Omega$ do not affect the expected variances.
Refer to caption: x8.png

## Appendix C Proofs omitted in the main text

###### Proof of Lemma 3.1 .

For a given $A\in{\mathcal{S}}$, since $A^{-1}\in{\mathcal{S}}$, we have $Z\overset{d}{=}A^{-1}Z:=\tilde{Z}$ and

$$ $\displaystyle\mathbb{E}_{Z}[f(Z)]=\mathbb{E}_{A^{-1}Z}[f(Z)]=\mathbb{E}_{ \tilde{Z}}[f(A\tilde{Z})]=\mathbb{E}_{Z}[f(AZ)].$ $$

This naturally leads to

$$ $\mathbb{E}_{Z}[\mathbb{E}_{A^{\prime}\sim\nu}[f(A^{\prime}Z)]]=\mathbb{E}_{A^{ \prime}\sim\nu}[\mathbb{E}_{Z}[f(A^{\prime}Z)]]=\mathbb{E}_{A^{\prime}\sim\nu} [\mathbb{E}_{Z}[f(Z)]]=\mathbb{E}_{Z}[f(Z)]$ $$

where the first equality comes from Fubini’s theorem and the integrability of $f$.
∎

###### Proof of Theorem 3.2 .

Since $\hat{\beta}=X^{\dagger}y$, we have
$\mathrm{Cov}(\hat{\beta}\mid X)=X^{\dagger}\mathrm{Cov}(y\mid X)X^{\dagger\top
}=X^{\dagger}\Omega X^{\dagger\top}$, which leads to the following expression for the variance component of prediction risk:

$$ $\displaystyle\mathrm{Var}_{\Sigma}(\hat{\beta}\mid X)$ $\displaystyle=\operatorname{Tr}(\mathrm{Cov}(\hat{\beta}\mid X)\Sigma)= \operatorname{Tr}(X^{\dagger}\Omega X^{\dagger\top}\Sigma)=\|SX^{\dagger}T\|_{ F}^{2}=\|BT\|_{F}^{2},$ $$

where $S=\Sigma^{1/2},T=\Omega^{1/2}$, and $B=SX^{\dagger}$.
Using the singular value decomposition (SVD) of $B$ and $T$, respectively, we can rewrite this as follows:

$$ $\displaystyle\|BT\|_{F}^{2}=\|UDV^{\top}U_{T}D_{T}V_{T}^{\top}\|_{F}^{2}{=\|DV ^{\top}U_{T}D_{T}\|_{F}^{2}},$ $$

where $B=UDV^{\top}$ and $T=U_{T}D_{T}V_{T}^{\top}$ with orthogonal matrices $U,V,U_{T},V_{T}$, and diagonal matrices $D,D_{T}$.
Now we need to compute the alignment $V^{\top}U_{T}$ of the right-singular vectors of $B$ with the left-eigenvectors of $T$.

$$ $\displaystyle{\|DV^{\top}U_{T}D_{T}\|_{F}^{2}}$ $\displaystyle=\sum_{i,j=1}^{n}\left(D_{ii}\sum\nolimits_{k=1}^{n}V^{\top}_{ik} (U_{T})_{kj}(D_{T})_{jj}\right)^{2}$ $\displaystyle{=\sum\nolimits_{i,j=1}^{n}\lambda_{i}(B)^{2}\lambda_{j}(T)^{2} \gamma_{ij}}$ $\displaystyle{=\sum\nolimits_{i,j=1}^{n}\lambda_{i}\left((X^{\top}X)^{\dagger} \Sigma\right)\lambda_{j}(\Omega)\gamma_{ij}}$ $\displaystyle=\underbracket{\lambda\left((X^{\top}X)^{\dagger}\Sigma\right)^{ \top}}_{1\times n}\underbracket{\Gamma(X)}_{n\times n}\underbracket{\lambda( \Omega)}_{n\times 1},$ $$

where $\gamma_{ij}:=\langle V_{:i},(U_{T})_{:j}\rangle^{2}\geq 0$, $\Gamma(X):=(\gamma_{ij})_{i,j}\in\mathbb{R}^{n\times n}$
and $\lambda(A)\in\mathbb{R}^{n}$ is a vector with its element $\lambda_{i}(A)$ as the $i$-th largest eigenvalue of $A$.

Therefore, we can rewrite the variance as $\mathrm{Var}_{\Sigma}(\hat{\beta}\mid X)=a(X)^{\top}\Gamma(X)b$ with

$$ $\displaystyle a(X)$ $\displaystyle:=\lambda\left((X^{\top}X)^{\dagger}\Sigma\right)\in\mathbb{R}^{n},$ $\displaystyle b$ $\displaystyle:=\lambda(\Omega)\in\mathbb{R}^{n},$ $\displaystyle\Gamma(X)_{ij}$ $\displaystyle=\gamma_{ij}=\langle v^{(i)},u^{(j)}\rangle^{2},$ $$

where $v^{(i)}:=V_{:i}$ and $u^{(j)}:=(U_{T})_{:j}$.
Note that the alignment matrix $\Gamma(X)$ is a doubly stochastic matrix since $\sum_{j}\gamma_{ij}=\sum_{i}\gamma_{ij}=1$ and $0\leq\gamma_{ij}\leq 1$.

Now, we want to compute the expected variance. To do so, from Lemma 3.1 with ${\mathcal{S}}=O(n)$, we can obtain

$$ $\displaystyle\mathbb{E}_{X}[a(X)^{\top}\Gamma(X)b]$ $\displaystyle=\mathbb{E}_{X}\left[\mathbb{E}_{O\sim\nu}[a(OX)^{\top}\Gamma(OX) b]\right]=\mathbb{E}_{X}\left[a(X)^{\top}\mathbb{E}_{O\sim\nu}[\Gamma(OX)]b \right],$ $$

where $\nu$ is the unique uniform distribution (the Haar measure) over the orthogonal matrices $O(n)$.
For an orthogonal matrix $O\in O(n)$,
we have

$$ $\displaystyle\Gamma(OX)_{ij}$ $\displaystyle=\langle Ov^{(i)},u^{(j)}\rangle^{2}=(v^{(i)\top}O^{\top}u^{(j)}) ^{2},$ $$

since $S(OX)^{\dagger}=SX^{\dagger}O^{\top}=BO^{\top}=UD(OV)^{\top}$.
Here, $(OX)^{\dagger}=X^{\dagger}O^{\top}$ follows from the orthogonality of $O\in O(n)$.
Since the Haar measure is invariant under the matrix multiplication in $O(n)$, if we take the expectation over
the Haar measure,
then we have

$$ $\displaystyle\bar{\Gamma}(X)_{ij}$ $\displaystyle:=\mathbb{E}_{O\sim\nu}[\Gamma(OX)_{ij}]=\mathbb{E}_{O\sim\nu}[(v ^{(i)\top}O^{\top}u^{(j)})^{2}]=\mathbb{E}_{O\sim\nu}[(v^{(i)\top}O^{\top}O^{( j)\top}u^{(j)})^{2}].$ $$

Here, for a given $j$, we can choose a matrix $O^{(j)}\in O(n)$ such that its first column is $u^{(j)}$ and $O^{(j)\top}u^{(j)}=e_{1}$, then $\bar{\Gamma}(X)_{ij}$ is independent of $j$ (say $\bar{\Gamma}(X)_{ij}=\alpha_{i}$). Since $\Gamma(X)$ is doubly stochastic, so is $\bar{\Gamma}(X)$ and we have
$\sum\nolimits_{j=1}^{n}\bar{\Gamma}(X)_{ij}=n\alpha_{i}=1$ which yields $\bar{\Gamma}(X)_{ij}=\alpha_{i}=1/n$, regardless of the distribution of $V$;
thus, $\bar{\Gamma}(X)=\frac{1}{n}J$, where $J_{ij}=1(i,j=1,2,\cdots,n)$.

Therefore, we have the expected variance as follows:

$$ $\displaystyle\mathbb{E}_{X}[\mathrm{Var}_{\Sigma}(\hat{\beta}\mid X)]$ $\displaystyle=\mathbb{E}_{X}[a(X)^{\top}\frac{1}{n}Jb]=\frac{1}{n}\sum_{i,j=1} ^{n}\mathbb{E}_{X}[a_{i}(X)]b_{j}=\frac{1}{n}\mathbb{E}_{X}[\operatorname{Tr}( (X^{\top}X)^{\dagger}\Sigma)]\operatorname{Tr}(\Omega).$ $$

∎

###### Proof of Corollary 4.2 .

Note that

$$ $\displaystyle\mathbb{E}_{X}[\mathrm{Var}(\hat{\beta}|X)]$ $\displaystyle=\frac{\operatorname{Tr}(\Omega)}{p}\mathbb{E}_{X}\left[\frac{1}{ n}\sum_{i}\frac{1}{\lambda_{i}}\right]$ $\displaystyle=\frac{\operatorname{Tr}(\Omega)}{p}\mathbb{E}_{X}\left[\int\frac {1}{s}dF^{XX^{\top}/p}(s)\right]$ $\displaystyle=\frac{\operatorname{Tr}(\Omega)}{n}\mathbb{E}_{X}\left[\int\frac {1}{s}dF^{XX^{\top}/n}(s)\right].$ $$

Then, the desired result follows directly from (5).
∎

###### Proof of ( 4 ).

The bias term of the prediction risk can be expressed as follows:

$$ $\displaystyle[\mathrm{Bias}_{\Sigma}(\hat{\beta}\mid X)]^{2}$ $\displaystyle=\|\mathbb{E}[\hat{\beta}\mid X]-\beta\|_{\Sigma}^{2}$ $\displaystyle=\|(\hat{\Sigma}^{\dagger}\hat{\Sigma}-I)\beta\|_{\Sigma}^{2}$ $\displaystyle=\beta^{\top}(I-\hat{\Sigma}^{\dagger}\hat{\Sigma})\Sigma(I-\hat{ \Sigma}^{\dagger}\hat{\Sigma})\beta$ $\displaystyle=\beta^{\top}\lim_{\lambda\searrow 0}\lambda(\hat{\Sigma}+\lambda I )^{-1}\Sigma\lim_{\lambda\searrow 0}\lambda(\hat{\Sigma}+\lambda I)^{-1}\beta$ $\displaystyle=(S\beta)^{\top}\lim_{\lambda\searrow 0}\lambda^{2}(S^{-1}\hat{ \Sigma}S+\lambda I)^{-2}S\beta,$ $$

where $\hat{\Sigma}=X^{\top}X/n$.
Here, the fourth equality comes from the equation

$$ $\displaystyle I-\hat{\Sigma}^{\dagger}\hat{\Sigma}$ $\displaystyle=\lim_{\lambda\searrow 0}I-(\hat{\Sigma}+\lambda I)^{-1}\hat{\Sigma}$ $\displaystyle=\lim_{\lambda\searrow 0}I-(\hat{\Sigma}+\lambda I)^{-1}(\hat{ \Sigma}+\lambda I-\lambda I)$ $\displaystyle=\lim_{\lambda\searrow 0}\lambda(\hat{\Sigma}+\lambda I)^{-1}.$ $$

∎

###### Proof of ( 7 ).

The RHS of (6) is bounded above by $\int\frac{1}{1+c_{H}s^{\ast}}dH(\tau)=\frac{1}{1+c_{H}s^{\ast}}$, and thus $1-\frac{1}{\gamma}\leq\frac{1}{1+c_{H}s^{\ast}}$, which yields $s^{\ast}\leq c_{H}^{-1}\frac{1}{\gamma-1}$.
We can similarly prove the other inequality in (7) with a lower bound $\frac{1}{1+C_{H}s^{\ast}}$ on the RHS of (6).
∎

###### Proof of ( 8 ).

To further explore the inequalities (7),
we rewrite (6) from Theorem 4.3 as follows:

$$ $\displaystyle 1-\frac{1}{\gamma}=\mathbb{E}_{\tau\sim H}\left[g(\tau;s^{\ast}) \right],\quad\text{where}\ g(t;s):=\frac{1}{1+ts}\ \text{for}\ t,s>0.$ $$

Here, since $g(t;s)$ is convex with respect to $t>0$ for a given $s>0$, by Jensen’s inequality, we then have

$$ $\mathbb{E}_{\tau\sim H}[g(\tau;\mu_{H}^{-1}s^{\ast}_{\text{iso}})]\geq g\left( \mu_{H};\mu_{H}^{-1}s^{\ast}_{\text{iso}}\right)=g(1;s^{\ast}_{\text{iso}})=1- \gamma^{-1}$ $$

where $\mu_{H}=\mathbb{E}_{\tau\sim H}[\tau]$.
Therefore, the limit Stieltjes transform $s^{\ast}$ in the anisotropic case should be larger than $\mu_{H}^{-1}s^{\ast}_{\text{iso}}$ of the isotropic case to satisfy $\mathbb{E}_{\tau\sim H}[g(\tau;s^{\ast})]=1-{\gamma}^{-1}$ since $g(t;s)$ is a decreasing function with respect to $s\geq 0$ when $t>0$.
This leads to a tighter lower bound $s^{\ast}\geq\mu_{H}^{-1}s^{\ast}_{\text{iso}}=\mu_{H}^{-1}(\gamma-1)^{-1}$ than (7) because $\mu_{H}\leq C_{H}$.
∎

###### Proof of Lemma 3.1 .

###### Proof of Theorem 3.2 .

###### Proof of Corollary 4.2 .

###### Proof of ( 4 ).

###### Proof of ( 7 ).

###### Proof of ( 8 ).

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: I have several comments and suggestions for the authors to address.


1. The paper is not complete. Please add the conclusion section in the revision. 

2. It is worthwhile to extend the proposed approach from the case  p>n   to the ultra-high dimensional case with p=exp(n). 

3. It is of interest for the authors to  compare   proposed method with existing ones  (including  Chinot et al. (2022) and Chinot & Lerasle (2023), etc.)  in the experiments.

4. The new idea about benefits of over-parameterization could be extended to time series, panel and grouped data.  It is helpful  for the authors to  elaborate it by proving more details.

5. There are some  typos, grammatical errors,  etc. in the paper. Please check it in the revision carefully. 

in page 1, line 022, "can extend" -> "can be extended". 

in page 5, line 312, "appendix" -> "Appendix".

in page 9, line 432, "appendix" -> "Appendix".

in page 10, line 500, "appendix" -> "Appendix".

in page 14, line 704, "Figure" -> "Figures".

in page 17, line 905, "6" -> "(6)".
soundness: 3
strengths: The paper is written clearly.  I appreciate the significant contribution to high-dimensial data analysis.  The new approach is  promising for 
 more  broad framework. This paper attacked the very challenging problem involved in least-squares estimators beyond the  assumiption with  i.i.d.  errors. The new idea about  benefits of  over-parameterization could be extended  to time series, panel and grouped data, etc.  with the broad impact.
summary: In this paper, the authors investigate prediction risk and  estimation risk under more general regression error assumptions beyond i.i.d. errors. In particlar, they explore the benefits of overparameterization in a more realistic setting, which allows for clustered or serial dependence. This paper  demonstrated   that the estimation difficulties associated with the variance components  can be summarized through the trace of the variance-covariance matrix of the regression errors.
weaknesses: The paper proved several  good theoretical results and properties. The experiemts of data set may not be comprehensive, Since there is no comparison of  the  proposed method with existing methods.   In addition,  the proposed approach only works for  the case  p>n. In practice,    the ultra-high dimensional case with p=exp(n) is very  common.  This is a more realistic setting in the big data era.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 4
questions: 1. Include discussions on relevant regularized least squares estimators, and
   provide scenarios where the ridgeless least squares estimator excels or
   underperforms.
2. Provide comparative experiments with other regularization techniques to
   demonstrate where ridgeless least squares excels or underperforms.
3. Add discussions on the limitations of the current investigation, particularly
   in terms of the assumptions.
4. Provide discussions connecting the results in the high-dimensional setting to
   those in the low-dimensional setting.
soundness: 4
strengths: Originality:

This paper extends the analysis of ridgeless least squares estimators by
relaxing the i.i.d. assumption on the regression errors. Correlated and
heteroscedastic errors are frequently encountered in practice but are not
investigated in prior works in the high-dimensional setting. This work fills a
gap in understanding the performance of least squares estimators, extending the
observation of double descent to more general settings.

Quality:

The paper is technically sound and comprehensive in its treatment of both
prediction and estimation risks. The finite-sample and asymptotic behaviors
finite-sample results are rigorously presented. The numerical result in Section
1 demonstrates the motivation of the generalized assumptions and the results in
Section 3.3 help validate the theoretical findings.

Clarity:

The paper is overall well written with clear presentation.

Significance

The paper extends the theory of least squares estimators to include the
situation with non i.i.d. errors in the context of overparameterization. It fill
a gap in the literature, so the results deserve to be documented.
summary: The paper extends the analysis of ridgeless least squares estimators to more
realistic error structures beyond the traditional i.i.d. assumptions. It
addresses both prediction and estimation risks in settings where regression
errors may be correlated or exhibit heteroscedasticity.

The authors provide exact finite-sample expressions for both types of risk,
which are decomposed into bias and variance components. The variance in
prediction risk is shown to be the product of two distinct terms: one related to
the error covariance matrix and the other dependent on the feature distribution.

The paper also conducts a systematic asymptotic analysis of prediction and
estimation risks. Numerical experiments with autoregressive and clustered
regression errors illustrate the theoretical findings.
weaknesses: For the high-dimensional setting considered in the paper, ridgeless least
squares estimators are seldom applied in practice. Regularization is almost
always helpful, and there is a large body of literature demonstrating the
advantages of regularized least squares. The paper should at least discuss this
line of research. In particular, it has been shown that optimized ridge
regression avoids bias inflation. Why, then, would practitioners care about
ridgeless least squares estimators in this scenario?

While the assumptions are more general than the i.i.d. case, the scenario
considered is still much simpler than what is encountered with real, practical
data. For example, the assumption of left-spherical symmetry for the design
matrix is limiting in practical scenarios. Most real-world datasets feature
asymmetrical or skewed distributions, and many variables in real data are
actually categorical, rather than numerical. It would help to add discussions on
the limitations of the current investigation in these contexts.

Most of the theoretical results, especially the finite-sample results, appear
similar to the low-dimensional or fixed-$p$ case. It would be helpful to include
discussions on this point. If the results are not the same, do they reduce to
the low-dimensional results? Even if they have the same expression, it would be
insightful to explain why it is nontrivial to derive these results in the
high-dimensional setting.

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 3
questions: I believe this paper could improve its significance if it can answer the following questions:

1. What new insights could one gain from Theorem 3.4, 3.5?

2. Related to Q1, the true noise covariance $\Omega$ is not observable in reality. Could the authors provide any examples or algorithms to approximate such noise covariance in real-world datasets? If it is impractical to approximate it, how could we still analyse the risk with your risk expression in $\Omega$? 

3. Could the authors extend their results to more general settings like kernel ridge regression or kernel gradient flow? Or at least explain what are the technical difficulties that might hinder the extension?

4. One Central idea of regularization is to balance out the effect of the noise in the labels. With the ridgeless linear regression setting, this paper unfortunately misses the opportunity to discuss the important interplay between the non-i.i.d. noise and regularization, which I find quite disappointing.
soundness: 3
strengths: This paper provides rigorous proof and experiments to validate their claim. The assumption on the noise and the input is more general than previous works.
summary: This paper studies the prediction and estimation risks of the ridgeless least square regression under a more general assumption on the noise. Specifically, consider the data
$$
y_i = x_i^\top\beta + \epsilon_i,\ i=1,...,n
$$
with target $\beta\in\mathbb{R}^p$ and $ \epsilon_i$ is the noise independent of $x$ with $\mathbb{E}[\epsilon]=0$ and $\mathbb{E}[\epsilon\epsilon^\top]=\Omega$ finite and positive definite. This includes the case of i.i.d. Gaussian noise, autoregressive noise and cluster noise.

This paper applies the classical bias-variance decomposition onto the risks and finds a closed form expression for the variance term for both risks.
weaknesses: However, my biggest concern is the significance of the contribution provided by this paper. 

As mentioned in line 132-133, this paper is not the first to consider non-i.i.d. Gaussian noise. Although this paper requires less assumptions on the noise than [1], it does not contain enough illustrations or interpretations on their main results (Theorem 3.4, 3.5) for "allowing potentially adversarial errors" as promised in line 142-143. Indeed, this paper does not explain how their main results could recover previous result with i.i.d. Gaussian noise or gain new insights with more general noise.

Also, it seems that the techniques used in the main results are rather standard and can be extended easily to more general settings like kernel ridge regression or kernel gradient flow, which could potentially increase the significance of this work.


Reference:
[1] Geoffrey Chinot and Matthieu Lerasle. On the robustness of the minimum ℓ2 interpolator.
Bernoulli, 2023.

Meta-review: This paper studies both the prediction risk as well as the estimation risk in non-trivial minimum l2 norm or 'ridgeless' least squares problems. In particular, they consider a general setting that leads to a better understanding of non-iid settings, including a relationship between estimating variance components and a trace representation. This treatment applies to several types of grouped data, as well as for error structures arising from longitudinal correlation. Reviewers are largely positive and I agree with their assessment of the strengths in this article.

Decision: Accept (Poster)