Optimal Shrinkage of Eigenvalues in the Spiked Covariance Model

Introduction

Suppose we observe $p$-dimensional Gaussian vectors $X_{i}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\mathcal{N}(0,\Sigma_{p})$, $i=1,\dots,n$, with $\Sigma=\Sigma_{p}$ the underlying $p$-by-$p$ population covariance matrix. To estimate $\Sigma$, we form the empirical (sample) covariance matrix $S=S_{n,p}=n^{-1}\sum_{i=1}^{n}X_{i}X^{\prime}_{i}$; this is the maximum likelihood estimator. Stein observed that the maximum likelihood estimator $S$ ought to be improvable by eigenvalue shrinkage.

In high dimensional problems, $p$ and $n$ are often of comparable magnitude. There, the maximum likelihood estimator is no longer a reasonable choice for covariance estimation and the need to shrink becomes acute.

In this paper, we consider a popular large $n$, large $p$ setting with $p$ comparable to $n$, and a set of assumptions about $\Sigma$ known as the Spiked Covariance Model . We study a variety of loss functions derived from or inspired by the literature, and show that to each “reasonable” nonlinearity $\eta$ there corresponds a well-defined asymptotic loss.

In the sibling problem of matrix denoising under a similar setting, it has been shown that there exists a unique asymptotically admissible shrinker . The same phenomenon is shown to exist here: for many different loss functions, we show that there exists a unique optimal nonlinearity $\eta^{*}$, which we explicitly provide. Perhaps surprisingly, $\eta^{*}$ is the only asymptotically admissible nonlinearity, namely, it offers equal or better asymptotic loss than that of any other choice of $\eta$, across all possible Spiked Covariance models.

Consider a sequence of covariance estimation problems, satisfying two basic assumptions.

The number of observations $n$ and the number of variables $p_{n}$ in the $n$-th problem follows the proportional-growth limit $p_{n}/n\rightarrow\gamma$, as $n\rightarrow\infty$, for a certain $0<\gamma\leq 1$.

The spiked model exhibits three important phenomena, not seen in classical fixed-$p$ asymptotics, that play an essential role in the construction of optimal estimators. Drawing on results from , we highlight:

Loss functions and optimal estimation. Now consider a class of estimators for the population covariance $\Sigma$, based on individual shrinkage of the sample eigenvalues. Specifically,

assuming such limit exists. If a nonlinearity $\eta^{*}$ satisfies

To give the flavor of results to be developed systematically later, we now look at four error measures in common use. The first three, based on the operator, Frobenius and nuclear norms, use the singular values $\sigma_{j}$ of $\hat{\Sigma}-\Sigma$:

The fourth is Stein’s loss, widely studied in covariance estimation .

Remark. The optimal shrinker also depends on $\gamma$, so we might write $\eta^{*}(\lambda,\gamma)$. In model [Asy($\gamma$)], one can use the same $\gamma$ for each problem size $n$. Alternatively, in the $n$-th problem, one might use $\gamma_{n}=p_{n}/n$. The former choice is simpler, as $\eta^{*}$ can be regarded as a univariate function of $\lambda$, and so we make it in Sections 1–6. The latter choice is preferable technically, and perhaps also in practice, when one has $p$ and $n$, but not $\gamma$. It does, however, require us to treat $\eta(\lambda,c)$ as a bivariate function – see Section 7.

2 Some key observations

The sections to follow construct a framework for evaluating and optimizing the asymptotic loss (1.8). We highlight here some observations that will play an important role. Beforehand, let us introduce a useful modification of (1.7) to a rank-aware shrinkage rule:

where the dimension $r$ of the spiked model is taken as known. While our main results concern estimators $\hat{\Sigma}_{\eta}$ that naturally do not require $r$ to be known in advance, it will be easier conceptually and technically to analyze rank-aware shrinkage rules as a preliminary step.

[Obs. 1] Simultaneous block diagonalization. (Lemmas 1 and 5). There exists a (random) basis $W$ such that

where $A_{i}$ and $B_{i}$ are square blocks of equal size $d_{i}$, and $\sum d_{i}=2r$. (Here and below, $A\oplus B$ denotes a block-diagonal matrix with blocks $A$ and $B$).

[Obs. 2] Decomposable loss functions. The loss functions (1.11) and many others studied below satisfy

or the corresponding equality with sum replaced by max.

[Obs. 3] Asymptotic deterministic loss. (Lemmas 3 and 7). For rank-aware estimators, when $\eta$ and $L$ are suitably continuous, almost surely

[Obs. 4] Asymptotic equivalence of losses. (Proposition 2). Conclusions derived for rank-aware estimators (1.13) carry over to the original estimators (1.7) because, under suitable conditions

3 Organization of the paper

Simultaneous Block-Diagonalization

We first develop [Obs. 1] in the simplest case, $r=1$, assumping a rank-aware shrinker. In general, the estimator $\hat{\Sigma}_{\eta}$ and estimand $\Sigma$ are not simultaneously diagonalizable. However, in the particular case that both are rank-one perturbations of the identity, we will see that simultaneous block diagonalization is possible.

Some notation is needed. We denote the eigenvalues and eigenvectors of the spectral decompostion $S_{n,p_{n}}=V\Lambda V^{\prime}$ by

Whenever possible, we supress the index $n$ and write e.g. $S$, $\lambda_{i}$ and $v_{i}$ instead. Similarly, we often write $\Sigma_{p}$ or even $\Sigma$ for $\Sigma_{p_{n}}$.

Let $\Sigma$ and $\hat{\Sigma}$ be (fixed, nonrandom) $p$-by-$p$ symmetric positive definite matrices with

where the fundamental $2\times 2$ matrices $A$ and $B$ are given by

Let $\Delta=\text{diag}(\eta,1,\ldots,1)=I+(\eta-1)e_{1}e_{1}^{\prime}$, where $e_{1}$ denotes the unit vector in the first co-ordinate direction. It is evident that

It is natural, then, to work in the “common” basis of $u_{1}$ and $v_{1}$. We apply one step of Gram-Schmidt if we can, setting

In the second–exceptional–case, $v_{1}=\pm u_{1}$, so we pick a convenient vector orthogonal to $u_{1}$. In either case, the columns of the $p\times 2$ matrix $W_{2}=[u_{1}\ z]$ are orthonormal and their span contains both $u_{1}$ and $v_{1}$. Now fill out $W_{2}$ to an orthogonal matrix $W=[W_{2}\ \,W_{2}^{\perp}]$. Observe now that if $y$ lies in the column span of $W_{2}$ and $\alpha$ is a scalar, then necessarily

The expressions (2.3) – (2.5) now follow from the rank one perturbation forms (2.6) along with

Decomposable Loss Functions

Here and below, by loss function $L_{p}$ we mean a function of two $p$-by-$p$ positive semidefinite matrix arguments obeying $L_{p}\geq 0$, with $L_{p}(A,B)=0$ if and only if $A=B$. A loss family is a sequence $L=\left\{L_{p}\right\}_{p=1}^{\infty}$, one for each matrix size $p$. We often write loss function and refer to the entire family. [Obs. 2] calls out a large class of loss functions which naturally exploit the simultaneously block-diagonalizability property of Lemma 1; we now develop this observation.

Orthogonal Invariance. We say the loss function $L_{p}(A,B)$ is orthogonally invariant if for each orthogonal $p$-by-$p$ matrix $O$,

For given $p$ and a given sequence of block sizes $\{d_{i}\}$ such that $\sum_{i}d_{i}=p$, consider block-diagonal matrix decompositions of $p$ by $p$ matrices $A$ and $B$ into blocks $A^{i}$ and $B^{i}$ of size $d_{i}$:

Sum-Decomposability and Max-Decomposability. We say the loss function $L_{p}(A,B)$ is sum-decomposable if for all decompositions (3.1),

We say that it is max-decomposable if if for all decompositions (3.1),

Clearly, such loss functions can exploit the simultaneous block diagonalization of Lemma 1. Indeed,

Reduction to Two-Dimensional Problem. Consider an orthogonally invariant loss function, $L_{p}$, which is sum- or max-decomposable. Suppose that $\Sigma$ and $\hat{\Sigma}$ satisfy (2.1) and (2.2) respectively. Then

Lemma 1 provides a change of basis $W$ yielding decompositions (2.3) and (2.4). From the invariance and decomposability hypotheses,

Asymptotic Loss in the Spiked Covariance Model

If $\eta$ is continuous, then the convergence results (1.2) and (1.5) imply that the principal block converges as $n\to\infty$. Specifically,

say, with the convergence occurring in all norms on $2\times 2$ matrices.

We say that a scalar function $\eta:[0,\infty)\to[1,\infty)$ is a bulk shrinker if $\eta(\lambda)=1$ when $\lambda\leq\lambda_{+}(\gamma)$, and a neighborhood bulk shrinker if for some $\epsilon>0$, $\eta(\lambda)=1$ whenever $\lambda\leq\lambda_{+}(\gamma)+\epsilon.$

The neighborhood bulk shrinker condition on $\eta$ is rather strong, but does hold for $\eta_{*}^{N}$ in (1.12), for example. (Note that our definitions ignore the lower bulk edge $\lambda_{-}(\gamma)$, which is of less interest in the spiked model.)

Furthermore, if (b) $\eta$ is a neighborhood bulk shrinker, then $L_{p_{n}}(\Sigma_{p_{n}},\hat{\Sigma}_{\eta})$ also has this limit a.s.

Each of the 26 losses considered in this paper satisfies conditions (a).

In the rank-aware case $\hat{\Sigma}_{\eta}=\hat{\Sigma}_{\eta,1}$ satisfies

where the limit on the right hand side follows from convergence (4.1) and the assumed continuity of $L_{2}$.

Examples of Decomposable Loss Functions

Many of the loss functions that appear in the literature are Pivot-Losses. They can be obtained via the following common recipe:

Pivots. A matrix pivot is a matrix-valued function $\Delta(A,B)$ of two real positive definitee matrices $A,B$ such that: (i) $\Delta(A,B)=0$ if and only if $A=B$, (ii) $\Delta$ is orthogonally equivariant and (iii) $\Delta$ respects block structure in the sense that

for any orthogonal matrix $O$ of the appropriate dimension.

Matrix pivots can be symmetric-matrix valued, for example $\Delta(A,B)=A-B$, but need not be, for example $\Delta(A,B)=A^{-1}B-I$.

Pivot-Losses. Let $g$ be a non-negative function of a symmetric matrix variable that is definite: $g(A)=0$ if and only if $A=0$, and orthogonally invariant: $g(O\Delta O^{\prime})=g(\Delta)$ for any orthogonal matrix $O$. A symmetric-matrix valued pivot $\Delta$ induces an orthgonally-invariant pivot loss

More generally, for any matrix pivot $\Delta$, set $|\Delta|=(\Delta^{\prime}\Delta)^{1/2}$ and define

An orthogonally invariant function $g$ depends on its matrix argument $\Delta$ or $|\Delta|$ only through its eigenvalues or singular values $\delta_{1},\ldots,\delta_{p}$. We abuse notation to write $g(\Delta)=g(\delta_{1},\ldots,\delta_{p})$. Observe that if $g$ has either of the forms

for some univariate $g_{1}$, then the pivot loss $L(A,B)=g(\Delta(A,B))$ (symmetric pivot) or $L(A,B)=g(|\Delta|(A,B))$ (general pivot) is respectively sum- or max-decomposable. In case $\Delta$ is symmetric, the two definitions agree so long as $g_{1}$ is an even function of $\delta$.

There are different strategies to derive sum-decomposable pivot-losses. First, we can use statistical discrepancies between the Normal distributions $\mathcal{N}(0,A)$ and $\mathcal{N}(0,B)$:

This is just twice the Kullback distance $D_{KL}(\mathcal{N}(0,B)||\mathcal{N}(0,A))$. Stein’s loss is a pivot-loss with respect to $\Delta(A,B)=A^{-1/2}BA^{-1/2}$ and $g(\Delta)={\rm tr}(\Delta-I)-\log\det(\Delta)=\sum_{i}g_{1}(\delta_{i})$, where $g_{1}(\delta)=\delta-1-\log\delta.$

Entropy/Divergence Losses: Because the Kullback discrepancy is not symmetric in its arguments, we may consider two other losses: reversing the arguments we get Entropy loss $L^{ent}(A,B)=L^{st}(B,A)$ and summing the Stein and Entropy losses gives divergence loss:

see . Each can be shown to be sum-decomposable, following the same argument as above.

This measures the statistical distinguishability of $\mathcal{N}(0,A)$ and $\mathcal{N}(0,B)$ based on independent observations, since $L^{aff}=\frac{1}{2}\log(\int\sqrt{\phi_{A}}\sqrt{\phi_{B}})$ with $\phi_{A}$ and $\phi_{B}$ the densities of $\mathcal{N}(0,A)$ and $\mathcal{N}(0,B)$. Hence convergence of affinity loss to zero is equivalent to convergence of the underlying densities in Hellinger or Variation distance. This is a pivot-loss w.r.t $\Delta(A,B)=A^{-1/2}BA^{-1/2}$ and

as is seen by setting $C=A^{-1/2}(A+B)B^{-1/2}$ and noting that $C^{\prime}C=(2I+\Delta+\Delta^{-1})$. Here, $g_{1}(\delta)=\tfrac{1}{4}\log(2+\delta+\delta^{-1})/4$.

Fréchet Discrepancy : Let $L^{fre}(A,B)={\rm tr}(A+B-2A^{1/2}B^{1/2})$. This measures the minimum possible mean-squared difference between zero-mean random vectors with covariances $A$ and $B$ respectively. This is a pivot-loss w.r.t $\Delta(A,B)=A^{1/2}-B^{1/2}$, and $g(\Delta)={\rm tr}(\Delta^{2})=\sum_{i}g_{1}(\delta_{i})$ with $g_{1}(\delta)=\delta^{2}$.

Second, we may obtain sum-decomposable pivot-losses $L(A,B)=g(\Delta(A,B))$ by simply taking $g$ to be one of the standard matrix norms:

Squared Error Loss : Let $L^{F,1}(A,B)=\|A-B\|_{F}^{2}$. This is a pivot-loss w.r.t $\Delta(A,B)=A-B$ and $g(\Delta)={\rm tr}\Delta^{\prime}\Delta=\sum_{i}g_{1}(\delta_{i})$ with $g_{1}(\delta)=\delta^{2}$.

Squared Error Loss on Precision : Let $L^{F,2}(A,B)=\|A^{-1}-B^{-1}\|_{F}^{2}$. This is a pivot-loss w.r.t $\Delta(A,B)=A^{-1}-B^{-1}$ and $g(\Delta)={\rm tr}\Delta^{\prime}\Delta$.

Nuclear Norm Loss. Let $L^{N,1}(A,B)=\|A-B\|_{*}$ where $\|\Delta\|_{*}$ denotes the nuclear norm of the matrix $\Delta$, i.e. the sum of its singular values. This is a pivot-loss w.r.t $\Delta(A,B)=A-B$ and $g(\Delta)=\sum_{i}|\delta_{i}|$.

Let $L^{F,3}(A,B)=\|A^{-1}B-I\|_{F}^{2}$. This is a pivot-loss w.r.t $\Delta(A,B)=A^{-1}B-I$. It was studied in and later work.

Let $L^{F,7}(A,B)=\|\log(A^{-1/2}BA^{-1/2})\|_{F}^{2}$, where $\log()$ denotes the matrix logarithmThe matrix logarithm transfers the matrices from the Riemannian manifold of symmetric positive semidefinite matrices to its tangent space at $A$. It can be shown that $L^{F,7}$ is the squared geodesic distance in this manifold. This metric between covariances has attracted attention, for example, in diffusion tensor MRI . . This is a pivot-loss w.r.t

2 Examples of Max-Decomposable Losses

Max-decomposable losses arise by applying the operator norm (the maximal singular value or eigenvalue of a matrix) to a suitable pivot. Here are a few examples:

Operator Norm Loss : Let $L^{O,1}(A,B)=\|A-B\|_{op}$. This is a pivot-loss w.r.t $\Delta(A,B)=A-B$ and $g(\Delta)=\|\Delta\|_{op}=\max_{i}\delta_{i}$.

Operator Norm Loss on Precision: Let $L^{O,2}(A,B)=\|A^{-1}-B^{-1}\|_{op}$. This is a pivot-loss w.r.t. $\Delta(A,B)=A^{-1}-B^{-1}$.

Condition Number Loss: Let $L^{O,7}(A,B)=\|\log(A^{-1/2}BA^{-1/2})\|_{op}$. This is a pivot-loss w.r.t $\Delta(A,B)=\log(A^{-1/2}BA^{-1/2})$, related to . In the spiked model discussed below, $L^{O,7}$ effectively measures the condition number of $A^{-1/2}BA^{-1/2}$.

We adopt the systematic naming scheme $L^{\mbox{norm},\mbox{pivot}}$ where $\mbox{norm}\in\{F,O,N\}$, and $\mbox{pivot}\in\{1,\dots,7\}$. This set of 21 combinations covers the previous matrix norm examples and adds some more. Together with Stein’s loss and the others based on statistical discrepancy mentioned above, we arrive at a set of 26 loss functions, Table 1, to be studied in this paper.

Optimal Shrinkage for Decomposable Losses

Below, we call formally optimal shrinkers simply “optimal”. By definition, the optimal shrinkage rule $\eta^{*}(\lambda\,;\,\gamma,L)$ is the unique admissible rule, in the asymptotic sense, among rules of the form $\hat{\Sigma}_{\eta}(S_{n,p})=V\eta(\Lambda)V^{\prime}$ in the single-spike model. In the single spiked model (and as we show later, generally in the spiked model) one never regrets using the optimal shrinker over any other (reasonably regular) univariate shrinker. In light of our results so far, an obvious characterization of an optimal shrinker is as follows.

Characterization of Optimal Shrinker. Let $L=\left\{L_{p}\right\}_{p=1}^{\infty}$ be a loss family. Define

Many of the 26 loss families discussed in Section 3 admit a closed form expression for the optimal shrinker; see Table 2. For others, we computed the optimal shrinker numerically, by implementing in software a solver for the simple scalar optimization problem (6.3). Figure 3 portrays the optimal shrinkers for our 26 loss functions. We refer readers interested in computing specific individual shrinkers to our reproducibility advisory at the bottom of this paper, and invite the reader to explore the code supplement , consisting of online resources and code we offer.

2 Optimal Shrinkers Collapse the Bulk

We first observe that, for any of the 26 losses considered, the optimal shrinker collapses the bulk to $1$. The following lemma is proved in the supplemental article :

3 Optimal Shrinkers by Computer

The scalar optimization problem (6.3) is easy to solve numerically, so that one can always compute the optimal shrinker at any desired value $\lambda$. In the code supplement we provide Matlab code to compute the optimal nonlinearity for each of the 26 loss families discussed. In the sibling problem of singular value shrinkage for matrix denoising, demonstrates numerical evaluation of optimal shrinkers for the Schatten-$p$ norm, where analytical derivation of optimal shrinkers appears to be impossible.

4 Optimal Shrinkers in Closed Form

If $\lambda\leq\lambda_{+}(\gamma)$ set $\eta^{*}(\lambda)=1$. Otherwise:

This asymptotic relationship reflects the classical fact that in finite samples, the top empirical eigenvalue is always biased upwards of the underlying population eigenvalue . Formally defining the (asymptotic) bias as

On the other hand, within the bottom branch, the effect is to shrink the bulk to 1. In terms of Definition 3 we see that $\eta^{*}$ is a bulk shrinker, but not a neighborhood bulk shrinker.

One might expect asymptotic debiasing from every loss function, but, perhaps surprisingly, precise asymptotic debiasing is exceptional. In fact, none of the other optimal nonlinearities in Table 2 is precisely debiasing.

In the supplemental article we also provide a detailed investigation of the large-$\lambda$ asymptotics of the optimal shrinkers, including their asymptotic slopes, asymptotic shifts and asymptotic percent improvement.

Beyond Formal Optimality

Assume that $\Sigma$ and $\hat{\Sigma}$ are fixed matrices with

Let $U_{r}$ and $V_{r}$ denote the $p$-by-$r$ matrices consisting of the top $r$ eigenvectors of $\Sigma$ and $\hat{\Sigma}$ respectively. Suppose that $[U_{r}\ V_{r}]$ has full rank $2r$, and consider the $QR$ decomposition

where $Q$ has $2r$ orthonormal columns and the $2r\times 2r$ matrix $R$ is upper triangular. Let $R_{2}$ denote the $2r\times r$ submatrix formed by the last $r$ columns of $R$. Fill out $Q$ to an orthogonal matrix $W=[Q\ Q^{\perp}]$. Then in the transformed basis we have the simultaneous block decompositions

We start with observations about the structure of $Q$ and $R$. Since the first $r$ columns of $Q$ are identically those of $U_{r}$, we let $Z_{r}$ be the $n$-by-$r$ matrix such that $Q=[U_{r}\,Z_{r}]$. For the same reason, $R$ has the block structure

where the matrices $R_{12}$ and $R_{22}$ satisfy $V_{r}=U_{r}R_{12}+Z_{r}R_{22}\,,$ so that

Since $V_{r}$ has orthogonal columns, we have

Let $H$ be a $p\times r$ matrix whose columns lie in the column span of $Q$ and let $\Delta$ be an $r\times r$ diagonal matrix. Observe that

say, since the columns of $Q^{\perp}$ are orthogonal to those of $H$.

and so both of the form $I+H\Delta H^{\prime}$, with $H=U_{r}$ and $V_{r}$ respectively. We find that

We can then compute the value of $C_{2r}$ in the two cases to be given by $\Sigma_{2r}^{\circ}$ and $\hat{\Sigma}_{2r}^{\circ}$ respectively, which establishes (7.1) and (7.2), and hence the lemma. ∎

We intend to apply Lemma 5 to $\Sigma$ and $\hat{\Sigma}=\hat{\Sigma}_{\eta,r}$, the “rank-aware” modification (1.13) of the estimator $\hat{\Sigma}_{\eta}$ in (1.7). Assume now that $\hat{\Sigma}$ and the $p\times r$ matrix $V_{r,n}$ formed by the top eigenvectors of $V$ are random.

The rank of $[U_{r}\ V_{r,n}]$ equals $2r$ almost surely.

Let $\Pi_{r}(V)$ be the projection that picks out the first $r$ columns of an orthogonal matrix $V$. For a fixed $r$-frame $U_{r}$, we consider the event

where the $2r\times 2r$ matrices $\Sigma_{2r},\hat{\Sigma}_{2r}$ satisfy

Suppose also that the family $L=\left\{L_{p}\right\}$ of loss functions is orthogonally invariant and sum- or max- decomposable, and that $B\to L_{2r}(A,B)$ is continuous. Then

If $\eta$ is a neighborhood bulk shrinker, then $L_{p}(\Sigma,\hat{\Sigma}_{\eta})$ also has this limit a.s.

This is the rank $r$ analog of Lemma 3. The optimal nonlinearity $\eta^{*}$ is continuous on $[0,\infty)$ for all losses except the operator norm ones, for which $\eta^{*}$ is continuous except at $\lambda=\lambda_{+}(\gamma)$. Our result (7.7) requires only continuity on $(\lambda_{+}(\gamma),\infty)$ and so is valid for all 26 loss functions, as is the deterministic limit (7.8) for the rank-aware $\hat{\Sigma}_{\eta,r}$. However, as we saw earlier, only the nuclear norm based loss functions yield optimal functions that are neighborhood bulk shrinkers. To show that (7.8) holds for $L_{p}(\Sigma,\hat{\Sigma}_{\eta})$ for most other important shrinkage functions, some further work is needed – see Section 7.1 below.

We apply Lemma 5 to $\Sigma$ and $\hat{\Sigma}_{\eta,r}$ on the set of probability 1 provided by Lemma 6. First, we rewrite (7.2) to show the subblocks of $R$:

To rewrite the limit in block diagonal form, let $\Pi_{2r}$ be the permutation matrix corresponding to the permutation defined by

Permuting rows and columns in (7.1) and (7.10) using $\Pi_{2r}$ to obtain

we obtain (7.7). Using (7.6), the orthogonal invariance and sum/max decomposability, along with the continuity of $L_{2r}(A,\cdot)$, we have

In this section we prove Proposition 2 below, whereby the asymtotic losses coincide for a given estimator sequence $\hat{\Sigma}_{\eta}$ and the rank-aware versions $\hat{\Sigma}_{\eta,r}$. This result is plausible because of two observations:

Null eigenvalues stick to the bulk, i.e. for $i\geq r+1$, most eigenvalues $\lambda_{in}\leq\lambda_{+}(\gamma)$ and the few exceptions are not much larger. Hence, if $\eta$ is a continuous bulk shrinker, we expect $\hat{\Sigma}_{\eta}$ to be close to $\hat{\Sigma}_{\eta,r}$,

under a suitable continuity assumption on the loss functions $L_{p}$, $L(\Sigma,\hat{\Sigma}_{\eta})$ should then be close to $L(\Sigma,\hat{\Sigma}_{\eta,r})$.

where the $(\mu_{in})$ are the eigenvalues of a white Wishart matrix $W_{p_{n}-r}(n,I)$.

The second step is a bound on eigenvalues of a white Wishart that exit the bulk. Before stating it, we return to an important detail introduced in the Remark concluding Section 1.1.

Definition 3 of a bulk shrinker depends on the parameter $\gamma=\lim p/n$ through $\lambda_{+}(\gamma)$. Making that dependence explicit, we obtain a bivariate function $\eta(\lambda,c)$. In model [Asy($\gamma$)]and in the $n$-th problem, we might use $\eta(\lambda,c_{n})$ either with $c_{n}=\gamma$ or $c_{n}=p/n$. For Proposition 1 below, it will be more natural to use the latter choice. We also modify Definition 3 as follows.

We call $\eta:[0,\infty)\times(0,1]\to[1,\infty)$ a jointly continuous bulk shrinker if $\eta(\lambda,c)$ is jointly continuous in $\lambda$ and $c$, satisfies $\eta(\lambda,c)=1$ for $\lambda\leq\lambda_{+}(c)$ and is dominated: $\eta(\lambda,c)\leq M\lambda$ for some $M$ and all $\lambda$.

The following result is proved in [58, Theorem 2(a)].

Let $(\mu_{in})_{i=1}^{N}$ denote the sample eigenvalues of a matrix distributed as $W_{N}(n,I)$, with $N/n\to\gamma>0$. Suppose that $\eta(\lambda,c)$ is a jointly continuous bulk shrinker and that $c_{n}-N/n=O(n^{-2/3})$. Then for $q>0$,

In the next proposition we adopt the convention that estimators $\hat{\Sigma}_{\eta}$ of (1.7) and $\hat{\Sigma}_{\eta,r}$ of (1.13) are constructed with a jointly continuous bulk shrinker, which we denote $\eta(\lambda,c_{n})$.

and so $L_{p}(\Sigma,\hat{\Sigma}_{\eta})$ converges in probability to the deterministic asymptotic loss (7.8).

In the left side of (7.13), substitute $A=\Sigma,B_{1}=\hat{\Sigma}_{\eta}$ and $B_{2}=\hat{\Sigma}_{\eta,r}$. By definition, $\hat{\Sigma}_{\eta}$ and $\hat{\Sigma}_{\eta,r}$ share the same eigenvectors. The components of $\eta_{1}-\eta_{2}$ then satisfy

We now use (7.11) to compare the eigenvalues $\lambda_{in}$ of the spiked model to those of a suitable white Wishart matrix to which Proposition 1 applies. The function $\eta^{\uparrow}(\mu,c)=\max\{\eta(\lambda,c),1\leq\lambda\leq\mu\}$ and is non-decreasing and jointly continuous. Hence $\eta(\lambda_{in},c_{n})\leq\eta^{\uparrow}(\lambda_{in},c_{n})\leq\eta^{\uparrow}(\mu_{i-r,n},c_{n})$, and so

with a corresponding bound for $q=\infty$. From continuity condition (7.13),

2 Asymptotic loss for discontinuous optimal shrinkers

where $W$ has a two point distribution in which

For the proof, write $\|\cdot\|$ for $\|\cdot\|_{\infty}$. Let $W=[W_{1}\ W_{2}]$ be the orthogonal change of basis matrix constructed in Lemma 7, with $W_{1}$ containing the first $2r$ columns. We treat the two losses $L^{O,1}$ and $L^{O,2}$ at once using an exponent $a=\pm 1$, and write $\eta^{a}(\lambda)$ for $\eta^{a}(\lambda,\gamma_{n})$. Thus, let

lies in the column span of $W_{1}$. We have $\hat{\Sigma}^{a}_{\eta}-\Sigma^{a}=\Psi_{n}+\Delta_{n}$, and the main task will be to show that for $a=\pm 1$,

Assuming the truth of this for now, let us derive the proposition. The quantities of interest in (7.14), (7.15) become

The rescaled noise eigenvalue $p^{2/3}(\lambda_{r+1,n}-\lambda_{+}(\gamma_{n}))\stackrel{{\scriptstyle\mathcal{D}}}{{\to}}\sigma(\gamma)W$ has a limiting real Tracy-Widom distribution with scale factor $\sigma(\gamma)>0$ [60, Prop. 5.8]. Hence, using the discontinuity of the optimal shrinker $\eta^{*}$, and the square root singularity from above

which leads to (7.15) and hence the main result.

It remains to prove (7.16). For a symmetric block matrix,

Apply this to $W^{\prime}(\Psi+\Delta)W$ with

We now show that $\|\Delta W_{1}\|\to_{P}0$. Using notation from Lemma 5,

Since $\Delta v_{k}=0$ for $k=1,\ldots,r$,

From (7.18) we have $\|\Delta_{n}\|=O_{P}(1)$. Since each $v_{i},i>r$ is uniformly distributed on $S^{p-1}$, a simple union bound based on (7.23) below yields

It remains to bound $N_{n}$. From the interlacing inequality (7.11),

From (7.21) and the preceding two paragraphs, we conclude that $\|\Delta U_{r}\|=O_{P}(p^{-1/2}\sqrt{\log p})$ and so $\|\Delta W_{1}\|\to_{P}0$.

Returning to (7.20), we deduce now that $\|B_{n}\|\leq\|\Delta W_{1}\|\to_{P}0$. From the definition of $W_{1}$ we have $\|W_{1}^{\prime}\Psi W_{1}\|=\|\Psi\|$ and hence the inequalities

Now observe that $\|C_{n}\|\leq\|\Delta_{n}\|$. Apply (7.19) to $W^{\prime}\Delta W$ to get

and hence that $\|C_{n}\|\geq\|\Delta_{n}\|-o_{P}(1)$. Thus $\|C_{n}\|=\|\Delta_{n}\|+o_{P}(1)$. Inserting these results into (7.20), we obtain

which completes the proof of (7.16), and hence of Proposition 3. ∎

Finally, we record a concentration bound for the uniform distribution on spheres. While more sophisticated results are known , an elementary bound suffices for us.

If $U$ is uniformly distributed on $S^{n-1}$ and $u\in S^{n-1}$ is fixed, then for $M>0$ and $n\geq 4$,

Since $U_{1}^{2}:=\langle U,u\rangle^{2}$ has the $\text{Beta}(\frac{1}{2},\frac{n-1}{2})$ distribution,

where by Gautschi’s inequality [62, 63, (5.6.4)]

Since $(1-x/m)^{m}<e^{-x}$ for $x,m>0$, and $4/n\geq 2/(n-2)$ for $n\geq 4$,

Optimality Among Equivariant Procedures

The notion of optimality in asymptotic loss, with which we have been concerned so far, is relatively weak. Also, the class of covariance estimators we have considered, namely procedures that apply the same univariate shrinker to all empirical eigenvalues, is fairly restricted.

Consider the much broader class of orthogonally-equivariant procedures for covariance estimation , in which estimates take the form $\hat{\Sigma}=V\,\Delta\,V^{\prime}$. Here, $\Delta=\Delta(\Lambda)$ is any diagonal matrix that depends on the empirical eigenvalues $\Lambda$ in possibly a more complex way than the simple scalar element-wise shrinkage $\eta(\Lambda)$ we have considered so far. One might imagine that the extra freedom available with more general shrinkage rules would lead to improvements in loss, relative to our optimal scalar nonlinearity; certainly the proposals of are of this more general type.

The smallest achievable loss by any orthogonally equivariant procedure is obtained with the “oracle” procedure $\hat{\Sigma}^{oracle}=V\,\Delta^{oracle}\,V^{\prime}$, where

the minimum being taken over diagonal matrices with diagonal entries $\geq 1$. Clearly, this optimal performance is not attainable, since the minimization problem explicitly demands perfect knowledge of $\Sigma$, precisely the object that we aim to recover. This knowledge is never available to us in practice – hence the label oracleThe oracle procedure does not attain zero loss since it is “doomed” to use the eigenbasis of the empirical covariance, which is a random basis corrupted by noise, to estimate the population covariance.. Nevertheless, this optimal performance is a legitimate benchmark.

where $\eta^{*}$ is the optimal shrinker for the losses $L^{F,1}$ or $L^{st}$ in Table 2.

In short, the shrinker $\eta^{*}()$, which has been designed to minimize the limiting loss, asymptotically delivers the same performance as the oracle procedure, which has the lowest possible loss, in finite-$n$, over the entire class of covariance estimators by arbitrary high-dimensional shrinkage rules. On the other hand, by definition, the oracle procedure outperforms every orthogonally-equivariant statistical estimator. We conclude that $\eta^{*}$ – as one such orthogonally-invariant estimator – is indeed optimal (in the sense of having the lowest limiting loss) among all orthogonally invariant procedures. While we only discuss the cases $L^{F,1}$ and $L^{st}$, we suspect that this theorem holds true for many of the 26 loss functions considered.

For both $L=L^{F,1}$ and $L^{st}$, we establish a decomposition

Here, $a$ is a constant depending only on the loss function,

Together (8.6) and (8.7) establish the Theorem.

Turning to the details, we begin by showing (8.3). For Frobenius loss, we have from our definitions and (8.2) that

It remains to verify (8.6) and (8.7). Theorem 1 says that for $1\leq i\leq r$,

which yields (8.6). From (8.5), we observe that in our two cases

From the previous two displays, we conclude

which is (8.7), and so completes the full proof. ∎

Observe that for each of the loss families we consider, $L_{p}(\sigma^{2}A,\sigma^{2}B)=\sigma^{2\kappa}L_{p}(A,B)$, where $\kappa\in\{-2,-1,0,1,2\}$ depends on the family $\{L_{p}\}$ alone. Hence

Again for each of the loss families we consider, almost surely,

We conclude that, using (9.2), any consistent sequence of estimators $\hat{\sigma}_{n}$ yields a sequence of shrinkers with the same asymptotic loss as the optimal shrinker for known $\sigma^{2}$. In other words, at least inasmuch as the asymptotic loss is concerned, under the spiked model, there is no penalty for not knowing $\sigma^{2}$.

Define, for a symmetric $p$-by-$p$ positive definite matrix $S$ with eigenvalues $\lambda_{1},\ldots,\lambda_{p}$ the quantity

where $\lambda_{med}$ is a median of $\lambda_{1},\ldots,\lambda_{p}$ and $\mu_{\gamma}$ is the median of the Marčenko-Pastur distribution, namely, the unique solution in $\lambda_{-}(\gamma)\leq x\leq\lambda_{+}(\gamma)$ to the equation

where as before $\lambda_{\pm}(\gamma)=(1\pm\sqrt{\gamma})^{2}$. Note that the median $\mu_{\gamma}$ is not available analytically but can easily be obtained numerically, for example using remarks on the Marčenko-Pastur cumulative distribution function included in SI. Now for a sequence $\left\{S_{n,p_{n}}\right\}$ of sample covariance matrices, define the sequence of estimators

In summary, using (9.1) (for $\sigma^{2}$ known) or (9.2) with (9.4) (for $\sigma^{2}$ unknown) one can use the optimal shrinkers for each of the loss families discussed above, designed for the case $\sigma=1$, to construct a shrinker that is optimal, for the same loss family, under the spiked model with common variance $\sigma^{2}\neq 1$.

Discussion

In this paper, we considered covariance estimation in high dimensions, where the dimension $p$ is comparable to the number of observations $n$. We chose a fixed-rank principal subspace, and let the dimension of the problem grow large. A different asymptotic framework for covariance estimation would choose a principal subspace whose rank is a fixed fraction of the problem dimension; i.e. the rank of the principal subspace is growing rather than fixed. (In the sibling problem of matrix denoising, compare the “spiked” setup with the “fixed fraction” setup of .)

In the fixed fraction framework, some of underlying phenomena remain qualitatively similar to those governing the spiked model, while new effects appear. Importantly, the relationships used in this paper, predicting the location of the top empirical eigenvalues, as well as the displacement of empirical eigenvectors, in terms of the top theoretical eigenvalues, no longer hold. Instead, a complex nonlinear relation exists between the limiting distribution of the empirical eigenvalues and the limiting distribution of the theoretical eigenvalues, as expressed by the Marčenko-Pastur (MP) relation between their Stieltjes transforms .

Covariance shrinkage in the proportional rank model should then, naturally, make use of the so-called MP Equation. Noureddine El Karoui proposed a method for debiasing the empirical eigenvalues, namely, for estimating (in a certain specific sense) their corresponding population eigenvalues; Olivier Ledoit and Sandrine Peché developed analytic tools to also account for the inaccuracy of empirical eigenvectors, and Ledoit and Michael Wolf have implemented such tools and applied them in this setting.

The proportional rank case is indeed subtle and beautiful. Yet, the fixed-rank case deserves to be worked out carefully. In particular, the shrinkers we have obtained here in the fixed-rank case are extremely simple to implement, requiring just a few code lines in any scientific computing language. In comparison, the covariance estimation ideas of , based on powerful and deep insights from MP theory, require a delicate, nontrivial effort to implement in software, and call for expertise in numerical analysis and optimization. As a result, the simple shrinkage rules we propose here may be more likely to be applied correctly in practice, and to work as expected, even in relatively small sample sizes.

An analogy can be made to shrinkage in the normal means problem, for example . In that problem, often a full Bayesian model applies, and in principle a Bayesian shrinkage would provide an optimal result . Yet, in applications one often wants a simple method which is easy to implement correctly, and which is able to deliver much of the benefit of the full Bayesian approach. In literally thousands of cases, simple methods of shrinkage - such as thresholding - have been chosen over the full Bayesian method for precisely that reason.

Reproducible Research

In the code supplement we offer a Matlab software library that includes:

A function to compute the value of each of the 26 optimal shrinkers discussed to high precision.

A function to test the correctness of each of the 18 analytic shrinker fomulas provided.

Scripts that generate each of the figures in this paper, or subsets of them for specified loss functions.

Acknowledgements

We thank Amit Singer, Andrea Montanari, Sourav Chatterjee and Boaz Nadler for helpful discussions. We also thank the anonymous referees for significantly improving the manuscript through their helpful comments. This work was partially supported by NSF DMS-0906812 (ARRA). MG was partially supported by a William R. and Sara Hart Kimball Stanford Graduate Fellowship.

Proofs and Additional Results

In the supplementary material we provide proofs omitted from the main text for space considerations and auxiliary lemmas used in various proofs. Notably, we prove Lemma 4, and provide detailed derivations of the 17 explicit formulas for optimal shrinkers, as summarized in Table 2. In addition, in the supplementary material we offer a detailed study of the large-$\lambda$ asymptotics (asymptotic slope and asymptotic shift) of the optimal shrinkers discovered in this paper, and tabulate the asymptotic behavior of each optimal shrinker. We also study the asymptotic percent improvement of the optimal shrinkers over naive hard thresholding of the sample covariance eigenvalues.

References