Exact expressions for double descent and implicit regularization via surrogate random design

Introduction

Classical statistical learning theory asserts that to achieve generalization one must use training sample size that sufficiently exceeds the complexity of the learning model, where the latter is typically represented by the number of parameters (or some related structural parameter; see Friedman et al., 2001). In particular, this seems to suggest the conventional wisdom that one should not use models that fit the training data exactly. However, modern machine learning practice often seems to go against this intuition, using models with so many parameters that the training data can be perfectly interpolated, in which case the training error vanishes. It has been shown that models such as deep neural networks, as well as certain so-called interpolating kernels and decision trees, can generalize well in this regime. In particular, Belkin et al. (2019a) empirically demonstrated a phase transition in generalization performance of learning models which occurs at an interpolation thershold, i.e., a point where training error goes to zero (as one varies the ratio between the model complexity and the sample size). Moving away from this threshold in either direction tends to reduce the generalization error, leading to the so-called double descent curve.

We build on methods from Randomized Numerical Linear Algebra (RandNLA) in order to obtain exact non-asymptotic expressions for the mean squared error (MSE) of the Moore-Penrose estimator (see Theorem 1). This provides a precise characterization of the double descent phenomenon for the linear regression problem. In obtaining these results, we are able to provide precise formulas for the implicit regularization induced by minimum norm solutions of under-determined training samples, relating it to classical ridge regularization (see Theorem 2). To obtain our precise results, we use a somewhat non-standard random design, based on a specially chosen determinantal point process (DPP), which we term surrogate random design. DPPs are a family of non-i.i.d. sampling distributions which are typically used to induce diversity in the produced samples (Kulesza and Taskar, 2012). Our aim in using a DPP as a surrogate design is very different: namely, to make certain quantities (such as the MSE) analytically tractable, while accurately preserving the underlying properties of the original data distribution. This strategy might seem counter-intuitive since DPPs are typically found most useful when they differ from the data distribution. However, we show both theoretically (Theorem 3) and empirically (Section 5), that for many commonly studied data distributions, such as multivariate Gaussians, our DPP-based surrogate design accurately preserves the key properties of the standard i.i.d. design (such as the MSE), and even matches it exactly in the high-dimensional asymptotic limit. In our analysis of the surrogate design, we introduce the concept of determinant preserving random matrices (Section 4), a class of random matrices for which determinant commutes with expectation, which should be of independent interest.

The noise $\xi=y(\mathbf{x})-\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{w}^{*}$ has mean and variance $\sigma^{2}$.

Under Assumptions 1 and 2, we can establish our first main result, stated as the following theorem, where we use $\mathbf{X}^{\dagger}$ to denote the Moore-Penrose inverse of $\mathbf{X}$.

If the response noise is homoscedastic (Assumption 1) and $\mu$ is in general position (Assumption 2), then for ${\bar{\mathbf{X}}}\sim S_{\mu}^{n}$ (Definition 3) and ${\bar{y}}_{i}=y(\bar{\mathbf{x}}_{i})$,

When Assumption 1 holds, then $\mathbf{v}_{\mu,y}=\mathbf{\Sigma}_{\mu}\mathbf{w}^{*}$, however ridge-regularized least squares is well-defined for much more general response models. Our second result makes a direct connection between the (expectation of the) unregularized minimum norm solution on the sample and the global ridge-regularized solution. While the under-determined regime (i.e., $n<d$) is of primary interest to us, for completeness we state this result for arbitrary values of $n$ and $d$. Note that, just like the definition of regularized least squares, this theorem applies more generally than Theorem 1, in that it does not require the responses to follow any linear model as in Assumption 1 (proof in Appendix D).

We illustrate this result in Figure 1b, plotting the norm of the expectation of the Moore-Penrose estimator. As for the MSE, our surrogate theory aligns well with the empirical estimates for i.i.d. Gaussian designs, showing that the shrinkage of the unregularized estimator in the under-determined regime matches the implicit ridge-regularization characterized by Theorem 2. While the shrinkage is a linear function of the sample size $n$ for isotropic features (i.e., $\mathbf{\Sigma}=\mathbf{I}$), it exhibits a non-linear behavior for other spectral decays. Such implicit regularization has been studied previously (see, e.g., Mahoney and Orecchia, 2011; Mahoney, 2012); it has been observed empirically for RandNLA sampling algorithms (Ma et al., 2015); and it has also received attention more generally within the context of neural networks (Neyshabur, 2017). While our implicit regularization result is limited to the Moore-Penrose estimator, this new connection (and others, described below) between the minimum norm solution of an unregularized under-determined system and a ridge-regularized least squares solution offers a simple interpretation for the implicit regularization observed in modern machine learning architectures.

Our exact non-asymptotic expressions in Theorem 1 and our exact implicit regularization results in Theorem 2 are derived for the surrogate design, which is a non-i.i.d. distribution based on a determinantal point process. However, Figure 1 suggests that those expressions accurately describe the MSE (up to lower order terms) also under the standard i.i.d. design $\mathbf{X}\sim\mu^{n}$ when $\mu$ is a multivariate Gaussian. As a third result, we verify that the surrogate expressions for the MSE are asymptotically consistent with the MSE of an i.i.d. design, for a wide class of distributions which include multivariate Gaussians.

with probability one as $d,n\to\infty$ with $n/d\to\bar{c}\in(0,\infty)\setminus\{1\}$.

The above result is particularly remarkable since our surrogate design is a determinantal point process. DPPs are commonly used in ML to ensure that the data points in a sample are well spread-out. However, if the data distribution is sufficiently regular (e.g., a multivariate Gaussian), then the i.i.d. samples are already spread-out reasonably well, so rescaling the distribution by a determinant has a negligible effect that vanishes in the high-dimensional regime. Furthermore, our empirical estimates (Figure 1) suggest that the surrogate expressions are accurate not only in the asymptotic limit, but even for moderately large dimensions. Based on a detailed empirical analysis described in Section 5, we conjecture that the convergence described in Theorem 3 has the rate of $O(1/d)$.

Related work

There is a large body of related work, which for simplicity we cluster into three groups.

Double descent. The double descent phenomenon has been observed empirically in a number of learning models, including neural networks (Belkin et al., 2019a; Geiger et al., 2019), kernel methods (Belkin et al., 2018a, 2019b), nearest neighbor models (Belkin et al., 2018b), and decision trees (Belkin et al., 2019a). The theoretical analysis of double descent, and more broadly the generalization properties of interpolating estimators, have primarily focused on various forms of linear regression (Bartlett et al., 2019; Liang and Rakhlin, 2019; Hastie et al., 2019; Muthukumar et al., 2019). Note that while we analyze the classical mean squared error, many works focus on the squared prediction error. Also, unlike in our work, some of the literature on double descent deals with linear regression in the so-called misspecified setting, where the set of observed features does not match the feature space in which the response model is linear (Belkin et al., 2019c; Hastie et al., 2019; Mitra, 2019; Mei and Montanari, 2019), e.g., when the learner observes a random subset of $d$ features from a larger population.

RandNLA and DPPs. Randomized Numerical Linear Algebra (Drineas and Mahoney, 2016, 2017) has traditionally focused on obtaining purely algorithmic improvements for tasks such as least squares regression, but there has been growing interest in understanding the statistical properties of these randomized methods (Ma et al., 2015; Raskutti and Mahoney, 2016). Determinantal point processes (Kulesza and Taskar, 2012) have been recently shown to combine strong worst-case regression guarantees with elegant statistical properties (Dereziński and Warmuth, 2017). However, these results are limited to the over-determined setting (Dereziński et al., 2018, 2019, 2019) and ridge regression (Dereziński and Warmuth, 2018; Dereziński et al., 2019). Our results are also related to recent work on using DPPs to analyze the expectation of the inverse (Dereziński and Mahoney, 2019) and generalized inverse (Mutný et al., 2019) of a subsampled matrix.

Implicit regularization. The term implicit regularization typically refers to the notion that approximate computation can implicitly lead to statistical regularization. See Mahoney and Orecchia (2011); Perry and Mahoney (2011); Gleich and Mahoney (2014) and references therein for early work on the topic; and see Mahoney (2012) for an overview. More recently, often motivated by neural networks, there has been work on implicit regularization that typically considered SGD-based optimization algorithms. See, e.g., theoretical results (Neyshabur et al., 2014; Neyshabur, 2017; Soudry et al., 2018; Gunasekar et al., 2017; Arora et al., 2019; Kubo et al., 2019) as well as extensive empirical studies (Martin and Mahoney, 2018, 2019). The implicit regularization observed by us is different in that it is not caused by an inexact approximation algorithm (such as SGD) but rather by the selection of one out of many exact solutions (e.g., the minimum norm solution). In this context, most relevant are the asymptotic results of Kobak et al. (2018) and LeJeune et al. (2019).

Surrogate random designs

In this section, we provide the definition of our surrogate random design $S_{\mu}^{n}$, where $\mu$ is a $d$-variate probability measure and $n$ is the sample size. This distribution is used in place of the standard random design $\mu^{n}$ consisting of $n$ row vectors drawn independently from $\mu$.

The above definition can be interpreted as rescaling the density function of $\mu^{K}$ by the pseudo-determinant, and then renormalizing it. We now construct our surrogate design $S_{\mu}^{n}$ by appropriately selecting the random variable $K$. The obvious choice of $K=n$ does not result in simple closed form expressions for the MSE in the under-determined regime (i.e., $n<d$), which is the regime of primary interest to us. Instead, we derive our random variables $K$ from the Poisson distribution.

if $n>d$, then $K\sim{\operatorname{Poisson}}(\gamma_{n})_{\geq d}$ with $\gamma_{n}=n-d$.

The first non-trivial property of the surrogate design $S_{\mu}^{n}$ is that the expected sample size is in fact always equal to $n$, which we prove in Appendix A.

Our general template for computing expectations under a surrogate design ${\bar{\mathbf{X}}}\sim\mathbf{S}_{\mu}^{n}$ is to use the following expressions based on the i.i.d. random design $\mathbf{X}\sim\mu^{K}$:

These formulas follow from Definitions 2 and 3 because the determinants $\det(\mathbf{X}\mathbf{X}^{\scriptscriptstyle{\top}})$ and $\det(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})$ are non-zero precisely in the regimes $n\leq d$ and $n\geq d$, respectively, which is why we can drop the restrictions on the range of the Poisson distribution. We compute the normalization constants by introducing the concept of determinant preserving random matrices, discussed in Section 4.

We focus here on the under-determined regime (i.e., $n<d$), highlighting the key new expectation formulas we develop to derive the MSE expressions for surrogate designs. A standard decomposition of the MSE yields:

Thus, our task is to find closed form expressions for the two expectations above. The latter, which is the expected projection onto the complement of the row-span of ${\bar{\mathbf{X}}}$, is proven in Appendix D.

No such expectation formula is known for i.i.d. designs, except when $\mu$ is an isotropic Gaussian. In Appendix D, we also prove a generalization of Lemma 2 which is then used to establish our implicit regularization result (Theorem 2). We next give an expectation formula for the trace of the Moore-Penrose inverse of the covariance matrix for a surrogate design (proof in Appendix C).

Determinant preserving random matrices

In this section, we introduce the key tool for computing expectation formulas of matrix determinants. It is used in our analysis of the surrogate design, and it should be of independent interest.

The key question motivating the following definition is: When does taking expectation commute with computing a determinant for a square random matrix?

A random $d\times d$ matrix $\mathbf{A}$ is called determinant preserving (d.p.), if

We next give a few simple examples to provide some intuition. First, note that every $1\times 1$ random matrix is determinant preserving simply because taking a determinant is an identity transfomation in one dimension. Similarly, every fixed matrix is determinant preserving because in this case taking the expectation is an identity transformation. In all other cases, however, Definition 4 has to be verified more carefully. Further examples (positive and negative) follow.

In fact, it can be shown that all random matrices with independent entries are determinant preserving. However, this is not a necessary condition.

To construct more complex examples, we show that determinant preserving random matrices are closed under addition and multiplication. The proof of this result is an extension of an existing argument, given by Dereziński and Mahoney (2019) in the proof of Lemma 7, for computing the expected determinant of the sum of rank-1 random matrices (proof in Appendix B).

If $\mathbf{A}$ and $\mathbf{B}$ are independent and determinant preserving, then:

$\mathbf{A}+\mathbf{B}$ is determinant preserving,

$\mathbf{A}\mathbf{B}$ is determinant preserving.

Next, we introduce another important class of d.p. matrices: a sum of i.i.d. rank-1 random matrices with the number of i.i.d. samples being a Poisson random variable. Our use of the Poisson distribution is crucial for the below result to hold. It is an extension of an expectation formula given by Dereziński (2019) for sampling from discrete distributions (proof in Appendix B).

If $K$ is a Poisson random variable and $\mathbf{A},\mathbf{B}$ are random $K\times d$ matrices whose rows are sampled as an i.i.d. sequence of joint pairs of random vectors, then $\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{B}$ is d.p., and so:

Finally, we show the expectation formula needed for obtaining the normalization constant of the under-determined surrogate design, given in (1). The below result is more general than the normalization constant requires, because it allows the matrices $\mathbf{A}$ and $\mathbf{B}$ to be different (the constant is obtained by setting $\mathbf{A}=\mathbf{B}=\mathbf{X}\sim\mu^{K}$). In fact, we use this more general statement to show Theorems 1 and 2. The proof uses Lemmas 4 and 5 (see Appendix B).

If $K$ is a Poisson random variable and $\mathbf{A}$, $\mathbf{B}$ are random $K\times d$ matrices whose rows are sampled as an i.i.d. sequence of joint pairs of random vectors, then

Empirical evaluation of asymptotic consistency

In this section, we empirically quantify the convergence rates for the asymptotic result of Theorem 3. We focus on the under-determined regime (i.e., $n<d$) and separate the evaluation into the bias and variance terms, following the MSE decomposition given in (2). Consider $\mathbf{X}=\mathbf{Z}\mathbf{\Sigma}^{1/2}$, where the entries of $\mathbf{Z}$ are i.i.d. standard Gaussian, and define:

When $\mu$ is a centered multivariate Gaussian and its covariance has a constant condition number, then, for $n/d$ fixed, the surrogate MSE satisfies: \big{|}\frac{\textnormal{MSE}[\mathbf{X}^{\dagger}\mathbf{y}]}{\mathcal{M}}-1\big{|}=O(1/d).

Conclusions

We derived exact non-asymptotic expressions for the MSE of the Moore-Penrose estimator in the linear regression task, reproducing the double descent phenomenon as the sample size crosses between the under- and over-determined regime. To achieve this, we modified the standard i.i.d. random design distribution using a determinantal point process to obtain a surrogate design which admits exact MSE expressions, while capturing the key properties of the i.i.d. design. We also provided a result that relates the expected value of the Moore-Penrose estimator of a training sample in the under-determined regime (i.e., the minimum norm solution) to the ridge-regularized least squares solution for the population distribution, thereby providing an interpretation for the implicit regularization resulting from over-parameterization.

We would like to acknowledge ARO, DARPA, NSF, ONR, and GFSD for providing partial support of this work. We also thank Zhenyu Liao for pointing out fruitful connections between our results and the asymptotic analysis of random matrix resolvents.

References

Appendix A Proof of Lemma 1

We first record an important property of the design $S_{\mu}^{d}$ which can be used to construct an over-determined design for any $n>d$. A similar version of this result was also previously shown by Dereziński et al. (2019) for a different determinantal design.

Let ${\bar{\mathbf{X}}}\sim S_{\mu}^{d}$ and $\mathbf{X}\sim\mu^{K}$, where $K\sim{\operatorname{Poisson}}(\gamma)$. Then the matrix composed of a random permutation of the rows from ${\bar{\mathbf{X}}}$ and $\mathbf{X}$ is distributed according to $S_{\mu}^{d+\gamma}$.

We now proceed with the proof of Lemma 1, where we establish that the expected sample size of $S_{\mu}^{n}$ is indeed $n$.

Proof of Lemma 1 The result is obvious when $n=d$, whereas for $n>d$ it is an immediate consequence of Lemma 7. Finally, for $n<d$ the expected sample size follows as a corollary of Lemma 2, which states that

Appendix B Proofs for Section 4

We use $\operatorname{\textnormal{adj}}(\mathbf{A})$ to denote the adjugate of $\mathbf{A}$, defined as follows: the $(i,j)$th entry of $\operatorname{\textnormal{adj}}(\mathbf{A})$ is $(-1)^{i+j}\det(\mathbf{A}_{[n]\backslash\{j\},[n]\backslash\{i\}})$. We will use two useful identities related to the adjugate: (1) $\operatorname{\textnormal{adj}}(\mathbf{A})=\det(\mathbf{A})\mathbf{A}^{-1}$ for invertible $\mathbf{A}$, and (2) $\det(\mathbf{A}+\mathbf{u}\mathbf{v}^{\scriptscriptstyle{\top}})=\det(\mathbf{A})+\mathbf{v}^{\scriptscriptstyle{\top}}\!\operatorname{\textnormal{adj}}(\mathbf{A})\mathbf{u}$ (see Fact 2.14.2 in Bernstein, 2011).

First, note that from the definition of an adjugate matrix it immediately follows that if $\mathbf{A}$ is determinant preserving then adjugate commutes with expectation for this matrix:

where $(*)$ used (4), i.e., the fact that for d.p. matrices, adjugate commutes with expectation. Crucially, through the definition of an adjugate this step implicitly relies on the assumption that all the square submatrices of $\mathbf{A}_{\mathcal{I},\mathcal{J}}$ are also determinant preserving. Iterating this, we get that $\mathbf{A}+\mathbf{Z}$ is d.p. for any fixed $\mathbf{Z}$. We now show the same for $\mathbf{A}+\mathbf{B}$:

where $(*)$ uses the fact that after conditioning on $\mathbf{B}$ we can treat it as a fixed matrix. Next, we show that $\mathbf{A}\mathbf{B}$ is determinant preserving via the Cauchy-Binet formula:

where recall that $\mathbf{A}_{\mathcal{I},*}$ denotes the submatrix of $\mathbf{A}$ consisting of its (entire) rows indexed by $\mathcal{I}$.

To prove Lemma 5, we will use the following lemma, many variants of which appeared in the literature (e.g., van der Vaart, 1965). We use the one given by Dereziński et al. (2019).

If the rows of random $k\times d$ matrices $\mathbf{A},\mathbf{B}$ are sampled as an i.i.d. sequence of $k\geq d$ pairs of joint random vectors, then

Here, we use the following standard shorthand: ${k^{\underline{d}}}=\frac{k!}{(k-d)!}=k\,(k-1)\dotsm(k-d+1)$. Note that the above result almost looks like we are claiming that the matrix $\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{B}$ is d.p., but in fact it is not because $k^{d}\neq{k^{\underline{d}}}$. The difference in those factors is precisely what we are going to correct with the Poisson random variable. We now present the proof of Lemma 5.

To prove Lemma 6, we use the following standard determinantal formula which is used to derive the normalization constant of a discrete determinantal point process.

For any $k\times d$ matrices $\mathbf{A},\mathbf{B}$ we have

Proof of Lemma 6 By Lemma 5, the matrix $\mathbf{B}^{\scriptscriptstyle{\top}}\mathbf{A}$ is determinant preserving. Applying Lemma 4 we conclude that $\mathbf{I}+\mathbf{B}^{\scriptscriptstyle{\top}}\mathbf{A}$ is also d.p., so

where $(*)$ follows from the exchangeability of the rows of $\mathbf{A}$ and $\mathbf{B}$, which implies that the distribution of $\mathbf{A}_{S,*}\mathbf{B}_{S,*}^{\scriptscriptstyle{\top}}$ is the same for all subsets $S$ of a fixed size $k$.

Appendix C Proof of Theorem 1

In this section we use $Z_{\mu}^{n}$ to denote the normalization constant that appears in (1) when computing an expectation for surrogate design $S_{\mu}^{n}$. We first prove Lemma 3.

If ${\bar{\mathbf{X}}}\sim S_{\mu}^{n}$ for $n<d$, then we have

Proof Let $\mathbf{X}\sim\mu^{K}$ for $K\sim{\operatorname{Poisson}}({\gamma_{n}})$. Note that if $\det(\mathbf{X}\mathbf{X}^{\scriptscriptstyle{\top}})>0$ then using the fact that $\det(\mathbf{A})\mathbf{A}^{-1}=\operatorname{\textnormal{adj}}(\mathbf{A})$ for any invertible matrix $\mathbf{A}$, we can write:

In the over-determined regime, a more general matrix expectation formula can be shown (omitting the trace). The following result is related to an expectation formula derived by Dereziński et al. (2019), however they use a slightly different determinantal design so the results are incomparable.

If ${\bar{\mathbf{X}}}\sim S_{\mu}^{n}$ and $n>d$, then we have

Proof Let $\mathbf{X}\sim\mu^{K}$ for $K\sim{\operatorname{Poisson}}(\gamma_{n})$. Assumption 2 implies that for $K\neq d-1$ we have

however when $k=d-1$ then (6) does not hold because $\det(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})=0$ while $\operatorname{\textnormal{adj}}(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})$ may be non-zero. It follows that:

where the first term in $(*)$ follows from Lemma 6 and (4), whereas the second term comes from Lemma 2.3 of Dereziński et al. (2019). Dividing both sides by $Z_{\mu}^{n}=\det(\gamma_{n}\mathbf{\Sigma}_{\mu})$ completes the proof.

Applying the closed form expressions from Lemmas 2, 3 and 11, we derive the formula for the MSE and prove Theorem 1 (we defer the proof of Lemma 2 to Appendix D).

Proof of Theorem 1 First, assume that $n<d$, in which case we have $\gamma_{n}=\frac{1}{\lambda_{n}}$ and moreover

While the expression given after $(a)$ is simpler than the one after $(b)$, the latter better illustrates how the MSE depends on the sample size $n$ and the dimension $d$. Now, assume that $n>d$. In this case, we have $\gamma_{n}=n-d$ and apply Lemma 11:

The case of $n=d$ was shown in Theorem 2.12 of Dereziński et al. (2019). This concludes the proof.

Appendix D Proof of Theorem 2

As in the previous section, we use $Z_{\mu}^{n}$ to denote the normalization constant that appears in (1) when computing an expectation for surrogate design $S_{\mu}^{n}$. Recall that our goal is to compute the expected value of ${\bar{\mathbf{X}}}^{\dagger}\bar{\mathbf{y}}$ under the surrogate design $S_{\mu}^{n}$. Similarly as for Theorem 1, the case of $n=d$ was shown in Theorem 2.10 of Dereziński et al. (2019). We break the rest down into the under-determined case $(n<d)$ and the over-determined case ($n>d$), starting with the former. Recall that we do not require any modeling assumptions on the responses.

Proof Let $\mathbf{X}\sim\mu^{K}$ for $K\sim{\operatorname{Poisson}}(\gamma_{n})$ and denote $y(\mathbf{x}_{i})$ as $y_{i}$. Note that when $\det(\mathbf{X}\mathbf{X}^{\scriptscriptstyle{\top}})>0$, then the $j$th entry of $\mathbf{X}^{\dagger}\mathbf{y}$ equals $\mathbf{f}_{j}^{\scriptscriptstyle{\top}}(\mathbf{X}\mathbf{X}^{\scriptscriptstyle{\top}})^{-1}\mathbf{y}$, where $\mathbf{f}_{j}$ is the $j$th column of $\mathbf{X}$, so:

If $\det(\mathbf{X}\mathbf{X}^{\scriptscriptstyle{\top}})=0$, then also $\det(\mathbf{X}\mathbf{X}^{\scriptscriptstyle{\top}}+\mathbf{y}\mathbf{f}_{j}^{\scriptscriptstyle{\top}})=0$, so we can write:

Proof of Lemma 2 We let $y(\mathbf{x})=x_{j}$ where $j\in[d]$ and apply Lemma 12 for each $j$, obtaining:

from which the result follows by simple algebraic manipulation.

We move on to the over-determined case, where the ridge regularization of adding the identity to $\mathbf{\Sigma}_{\mu}$ vanishes. Recall that we assume throughout the paper that $\mathbf{\Sigma}_{\mu}$ is invertible.

Proof Let $\mathbf{X}\sim\mu^{K}$ for $K\sim{\operatorname{Poisson}}(\gamma_{n})$ and denote $y_{i}=y(\mathbf{x}_{i})$. Similarly as in the proof of Lemma 12, we note that when $\det(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})>0$, then the $j$th entry of $\mathbf{X}^{\dagger}\mathbf{y}$ equals $\mathbf{e}_{j}^{\scriptscriptstyle{\top}}(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-1}\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{y}$, where $\mathbf{e}_{j}$ is the $j$th standard basis vector, so:

If $\det(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})=0$, then also $\det(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X}+\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{y}\mathbf{e}_{j}^{\scriptscriptstyle{\top}})=0$. We proceed to compute the expectation:

where $(*)$ uses Lemma 5 twice (the first time, with $\mathbf{A}=\mathbf{X}$ and $\mathbf{B}=\mathbf{X}+\mathbf{y}\mathbf{e}_{j}^{\scriptscriptstyle{\top}}$). Dividing both sides by $Z_{\mu}^{n}=\det(\gamma_{n}\mathbf{\Sigma}_{\mu})$ concludes the proof.

We combine Lemmas 12 and 13 to obtain the proof of Theorem 2.

Proof of Theorem 2 The case of $n=d$ follows directly from Theorem 2.10 of Dereziński et al. (2019). Assume that $n<d$. Then we have $\gamma_{n}=\frac{1}{\lambda_{n}}$, so the result follows from Lemma 12. If $n>d$, then the result follows from Lemma 13.

Appendix E Proof of Theorem 3

The proof of Theorem 3 follows the standard decomposition of MSE in Equation 2, and in the process, establishes consistency of the variance and bias terms independently. To this end, we introduce the following two useful lemmas that capture the limiting behavior of the variance and bias terms, respectively.

Under the setting of Theorem 3, we have, as $n,d\to\infty$ with $n/d\to\bar{c}\in(0,\infty)\setminus\{1\}$ that

while $\mathbf{I}-\mathbf{X}^{\dagger}\mathbf{X}=0$ for $\bar{c}>1$.

For $n<d$, we first establish (1) $\liminf_{n}\lambda_{n}>0$ and (2) $\alpha_{n}\to 0$. To prove (1), by hypothesis $\mathbf{\Sigma}\succeq c\mathbf{I}$ for all $d$. Since $\frac{n}{d}<1$, we have (by definition of $\lambda_{n}$) for some $\delta>0$

Rearranging, we have $\lambda_{n}>\frac{\delta c}{1-\delta}>0$. For (2), let $(\tau_{i})_{i\in[d]}$ denote the eigenvalues of $\mathbf{\Sigma}$. Since $1-x\leq e^{-x}$ and $C\mathbf{I}\succeq\mathbf{\Sigma}\succeq c\mathbf{I}$ for all $d$,

and since $\lambda_{n}>0$ eventually as $d\to\infty$ we have $\alpha_{n}\to 0$ so that $(1-\alpha_{n})\lambda_{n}^{-1}-\lambda_{n}^{-1}\to 0$.

almost surely for some $m(z)$ to be defined; (ii) show that both $\delta(z)$ and its derivate $\delta^{\prime}(z)$ are uniformly bounded (by some quantity independent of $z>0$) so that by Arzela-Ascoli theorem, $\delta(z)$ converges uniformly to its limit and we are allowed to take $z\to 0^{+}$ in (10) and state

almost surely, given that the limit $\lim_{z\to 0^{+}}m(z)\equiv m(0)$ exists and eventually (iii) exchange the two limits in (11) with Moore-Osgood theorem, to reach

Step (i) follows from Silverstein and Bai (1995) that, we have, for $z>0$ that

almost surely as $n,d\to\infty$, for $m(z)$ the unique positive solution to

For the above step (ii), we use the assumption $\mathbf{\Sigma}\succeq c\mathbf{I}\succ 0$ for all $d$ large, so that with $\mathbf{X}=\mathbf{Z}\mathbf{\Sigma}^{\frac{1}{2}}$, we have for large enough $n,d$ that

almost surely, where we used Bai-Yin theorem Bai et al. (1993), which states that the minimum eigenvalue of $\mathbf{Z}\mathbf{Z}^{\scriptscriptstyle{\top}}/n$ is almost surely larger than $(\sqrt{\bar{c}}-1)^{2}/2$ for $n<d$ sufficiently large. Note that here the case $\bar{c}=1$ is excluded.

and similarly for its derivative, so that we are allowed to take the $z\to 0^{+}$ limit. Note that the existence of the $\lim_{z\to 0^{+}}m(z)$ for $m(z)$ defined in (12) is well known, see for example Ledoit and Péché (2011). Then, by Moore-Osgood theorem we finish step (iii) and by concluding that

E.1.2 The c¯∈(1,∞)¯𝑐1\bar{c}\in(1,\infty) case

First note that as $n,d\to\infty$ with $n>d$, we have $\beta_{n}=e^{d-n}\to 0$ and it it suffices to show

In the $\bar{c}\in(1,\infty)$ case, it is more convenient to work on the following co-resolvent

almost surely as $n,d\to\infty$, for $m(z)$ the unique solution to

so that for $d<n$ by taking $z=0$ we have

The steps (ii) and (iii) follow exactly the same line of arguments as the $\bar{c}<1$ case and are thus omitted.

E.2 Proof of lemma 15

Since $\mathbf{X}^{\dagger}\mathbf{X}=\mathbf{X}^{\scriptscriptstyle{\top}}(\mathbf{X}\mathbf{X}^{\scriptscriptstyle{\top}})^{\dagger}\mathbf{X}$, to prove lemma 15, we are interested in the limiting behavior of the following quadratic form

For the proof of Lemma 15 we follow the same protocol as that of Lemma 14, namely: (i) we consider, for fixed $z>0$, the limiting behavior of $\frac{1}{n}\mathbf{w}^{\scriptscriptstyle{\top}}\mathbf{X}^{\scriptscriptstyle{\top}}\left(\frac{1}{n}\mathbf{X}\mathbf{X}^{\scriptscriptstyle{\top}}+z\mathbf{I}_{n}\right)^{-1}\mathbf{X}\mathbf{w}$. Note that

and it remains to work on the second $z\mathbf{w}^{\scriptscriptstyle{\top}}\left(\frac{1}{n}\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X}+z\mathbf{I}_{d}\right)^{-1}\mathbf{w}$ term. It follows from Hachem et al. (2013) that

almost surely as $n,d\to\infty$, where we recall $m(z)$ is the unique solution to (12).

We move on to step (ii), under the assumption that $c\leq\lambda_{\min}(\mathbf{\Sigma})\leq\lambda_{\max}(\mathbf{\Sigma})\leq C$ and $\|\mathbf{w}\|\leq C^{\prime}$, we have

so that $\delta(z)$ remains bounded and similarly for its derivative $\delta^{\prime}(z)$, which, by Arzela-Ascoli theorem, yields uniform convergence and we are allowed to take the $z\to 0^{+}$ limit. Ultimately, in step (iii) we exchange the two limits with Moore-Osgood theorem, concluding the proof.

E.3 Finishing the proof of Theorem 3

To finish the proof of Theorem 3, it remains to write

Appendix F Additional details for empirical evaluation

For analyzing the rate of decay of variance and bias discrepancies (as defined in Section 5), it suffices to only consider diagonal covariance matrices $\mathbf{\Sigma}$. This is because if $\mathbf{\Sigma}=\mathbf{Q}\mathbf{D}\mathbf{Q}^{\scriptscriptstyle{\top}}$ is its eigendecomposition and $\mathbf{X}\sim\mathcal{N}_{n,d}(\mathbf{0},\mathbf{I}_{n}\otimes\mathbf{Q}\mathbf{D}\mathbf{Q}^{\scriptscriptstyle{\top}})$, then we have for $\mathbf{W}\sim\mathcal{P}\mathcal{W}(\mathbf{\Sigma},n)$ that $\mathbf{W}\overset{d}{=}\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X}$ and hence, defining $\widetilde{\mathbf{X}}\sim\mathcal{N}_{n,d}(\mathbf{0},\mathbf{I}_{n}\otimes\mathbf{D})$, by linearity and unitary invariance of trace,

and apply existing methods for constructing bootstrapped operator norm confidence intervals described in Lopes et al. (2019). To ensure that estimation noise is sufficiently small, we continually increase the number of Monte Carlo samples until the bootstrap confidence intervals are within $\pm 12.5\%$ of the measured discrepancies. We found that while variance discrepancy required a relatively small number of trials (up to one thousand), estimation noise was much larger for the bias discrepancy, and it necessitated over two million trials to obtain good estimates near $d=100$.

Letting $\lambda_{i}(\mathbf{\Sigma})$ be the $i$th largest eigenvalue of $\mathbf{\Sigma}$, we consider the following eigenvalue profiles (visualized in Figure 3):

diag_linear: linear decay, $\lambda_{i}(\mathbf{\Sigma})=b-ai$;

diag_exp: exponential decay, $\lambda_{i}(\mathbf{\Sigma})=b\,10^{-ai}$;

diag_poly: fixed-degree polynomial decay, $\lambda_{i}(\mathbf{\Sigma})=(b-ai)^{2}$;

diag_poly_2: variable-degree polynomial decay, $\lambda_{i}(\mathbf{\Sigma})=bi^{-a}$.

The constants $a$ and $b$ are chosen to ensure $\lambda_{\text{max}}(\mathbf{\Sigma})=1$ and $\lambda_{\text{min}}(\mathbf{\Sigma})=10^{-4}$ (i.e., the condition number $\kappa(\mathbf{\Sigma})=10^{4}$ remains constant).