Stein Risk Estimation for Eigenvalue Shrinkage

Why Estimate Risk From the Observed Covariance?

A covariance estimator can look plausible and still shrink too little or too much. Held-out likelihood is one way to compare choices; Stein risk estimation offers a different route when the sampling model is Gaussian. It uses the observed sample covariance and the way an estimator changes with that covariance to estimate its expected loss.

Here the loss is Stein’s loss, not squared Frobenius loss. It measures relative covariance error through a trace and a log determinant. The distinction matters: a rule minimizing an estimate of Stein risk need not minimize an estimate of Frobenius risk.

Definition: For positive-definite covariance matrix $\Sigma$ and positive-definite estimate $A$, Stein’s loss is $$ L(A,\Sigma)=\operatorname{tr}(\Sigma^{-1}A)-\log|\Sigma^{-1}A|-p. $$ Its risk is $R(\Sigma,\delta)=E_\Sigma[L(\delta(S),\Sigma)]$, where $\delta(S)$ is the covariance estimator computed from the sample covariance $S$.

The determinant expression means $\log|\Sigma^{-1}A|=\log|A|-\log|\Sigma|$. Thus the loss is dimensionless and penalizes both an incorrect overall scale and an incorrect covariance shape. We will estimate its repeated-sampling expectation, not claim that one observed value reveals the realized loss.

Model, Estimator, and Regularity Conditions

Assume $X_1,\ldots,X_n$ are independent $p$-dimensional observations with $X_k\sim N_p(0,\Sigma)$, where the mean is known to be zero and $\Sigma$ is positive definite. Define

S=\frac{1}{n}\sum_{k=1}^n X_kX_k^\mathsf{T}.

This tutorial uses the known-mean model; estimating the mean changes the sampling degrees of freedom and is a separate issue. When $n\ge p$, $S$ is positive definite with probability one. We consider orthogonally equivariant estimators: if $S=U\operatorname{diag}(l_1,\ldots,l_p)U^\mathsf{T}$, then

\delta(S)=U\operatorname{diag}\{d_1(l),\ldots,d_p(l)\}U^\mathsf{T}.

The eigenvalues $l_i$ are positive, and the estimator eigenvalues $d_i(l)$ must also be positive so that Stein loss is defined. The eigenvalues are taken in any fixed order; the functions must transform consistently when their arguments are permuted.

Regularity conditions: The result below is used for $n>p+1$, with $S$ positive definite almost surely. The matrix map $\delta$ must be continuously differentiable on the positive-definite cone, have positive eigenvalues, and satisfy the integrability and boundary conditions for Gaussian integration by parts. In particular, the relevant trace, divergence, and log-determinant terms must have finite expectations, and boundary terms at the edge of the cone must vanish. These conditions hold for the smooth affine shrinkage rules used below. At repeated sample eigenvalues, the spectral map and its divided differences must have continuous extensions.

The Stein Risk Estimate

For distinct sample eigenvalues, define the spectral divergence

D_\delta(l)= \sum_{i=1}^p\frac{\partial d_i}{\partial l_i} +\sum_{i<j}\frac{d_i-d_j}{l_i-l_j}.

The quotient accounts for how the eigenvectors, as well as the eigenvalues, change when the sample covariance is perturbed. At a repeated eigenvalue it is interpreted by continuity. With this convention, the Gaussian Stein–Haff identity for the scaled sample covariance is

E_\Sigma\{\operatorname{tr}(\Sigma^{-1}\delta(S))\} = E_\Sigma\left\{ \frac{n-p-1}{n}\operatorname{tr}(S^{-1}\delta(S)) +\frac{2}{n}D_\delta(l) \right\}.

Theorem: Under the model and regularity conditions above, an unbiased estimate of Stein risk is $$ \widehat R_{\mathrm{Stein}}(S;\delta) = \frac{n-p-1}{n}\operatorname{tr}(S^{-1}\delta(S)) +\frac{2}{n}D_\delta(l) +\log|S|-\log|\delta(S)|-c_{n,p}-p, $$ where $$ c_{n,p}=\sum_{i=1}^p\psi\!\left(\frac{n+1-i}{2}\right)+p\log 2-p\log n. $$ Here $\psi$ is the digamma function, $\psi(x)=\frac{d}{dx}\log\Gamma(x)$. In particular, $E_\Sigma[\widehat R_{\mathrm{Stein}}(S;\delta)]=R(\Sigma,\delta)$.

To justify the theorem, first consider the Wishart density of $S$, proportional on the positive-definite cone to $\lvert S\rvert^{(n-p-1)/2}\exp\{-\tfrac n2\operatorname{tr}(\Sigma^{-1}S)\}$. Integration by parts in the symmetric-matrix coordinates, with the stated boundary conditions, gives the Stein–Haff identity. In spectral coordinates the derivative contribution is $D_\delta$: a diagonal perturbation changes $d_i$ at rate $\partial d_i/\partial l_i$, while an off-diagonal perturbation mixing eigenvectors $i$ and $j$ contributes the divided difference $(d_i-d_j)/(l_i-l_j)$. The scaling $S=W/n$ gives the factors $1/n$ and $2/n$ in the identity.

Next, $nS$ has a Wishart distribution with $n$ degrees of freedom and scale $\Sigma$. Its Bartlett decomposition gives $\log|nS|=\log|\Sigma|+\sum_{i=1}^p\log Q_i$, where the independent $Q_i$ have chi-square distributions with $n+1-i$ degrees of freedom. Since $E[\log Q_i]=\psi((n+1-i)/2)+\log 2$, it follows that

E_\Sigma[\log|S|]=\log|\Sigma|+c_{n,p}.

Finally, expand the loss: $L(\delta(S),\Sigma)=\operatorname{tr}(\Sigma^{-1}\delta(S))-\log|\delta(S)|+\log|\Sigma|-p$. Replace its trace expectation using Stein–Haff and its $\log|\Sigma|$ expectation using $E[\log|S|]-c_{n,p}$. The result is exactly the displayed risk estimate, so its expectation equals the risk. $\square$

Worked Examples

Worked Example: Choosing Among Three Shrinkage Levels

Suppose $p=2$, $n=8$, and the observed sample eigenvalues are $l_1=4$ and $l_2=1$. Consider shrinking linearly toward the average sample eigenvalue: $d_i=(1-\lambda)l_i+\lambda\bar l$, where $\bar l=(l_1+l_2)/2=2.5$. Compare $\lambda=0,\frac12,1$.

For this rule, $\partial d_i/\partial l_i=1-\lambda+\lambda/2$, and $(d_1-d_2)/(l_1-l_2)=1-\lambda$. Therefore $D_\delta=3-2\lambda$. Also, $\operatorname{tr}(S^{-1}\delta)=(1-\lambda)2+\lambda(2.5)(1/4+1)=2+1.125\lambda$. For these $n,p$, $c_{8,2}\approx-0.413315$, so the constant $-c_{8,2}-p$ is approximately $-1.586685$.

At $\lambda=0$, $\delta=S$, the trace term is $2$, the divergence is $3$, and the log-determinant difference is zero. Hence $\widehat R=(5/8)(2)+(2/8)(3)-1.586685=2-1.586685\approx0.413315$.

At $\lambda=\frac12$, the estimate eigenvalues are $3.25$ and $1.75$, with determinant $5.6875$. The trace term is $2+1.125/2=2.5625$, and $D_\delta=2$. The log-determinant difference is $\log(4/5.6875)=\log(64/91)\approx-0.351976$. Thus \[ \widehat R=(5/8)(2.5625)+(2/8)(2)-0.351976-1.586685 \approx0.162901. \]

At $\lambda=1$, both estimate eigenvalues are $2.5$, so $\operatorname{tr}(S^{-1}\delta)=2.5(1/4+1)=3.125$, $D_\delta=1$, and $\log(|S|/|\delta|)=\log(4/6.25)=\log(0.64)\approx-0.446287$. Therefore \[ \widehat R=(5/8)(3.125)+(2/8)(1)-0.446287-1.586685 \approx0.170153. \]

Among these three candidates, $\lambda=\frac12$ has the smallest estimated Stein risk. This selects a candidate from the stated grid; it does not assert that the continuous risk is minimized exactly at one half.

Worked Example: Checking the Divergence for the Identity Rule

For the unshrunk estimator $\delta(S)=S$, $d_i=l_i$. Each derivative $\partial d_i/\partial l_i$ is $1$, and each divided difference is also $1$. Thus $D_\delta=p+\binom{p}{2}=p(p+1)/2$. This check catches a common factor-of-two error: in dimension two the divergence is $2+1=3$, not $2$.

For $p=2,n=8$ and eigenvalues $4,1$, the identity estimate has $\operatorname{tr}(S^{-1}\delta)=2$, and its log-determinant difference is zero. Substitution gives $\widehat R=(5/8)(2)+(2/8)(3)-1.586685=0.413315$. The value is an unbiased risk estimate for this observed sample; it is not required to be nonnegative, even though the underlying loss is nonnegative.

Worked Example: Evaluating a Shrinkage Rule at a Second Sample Size

Let $p=2,n=10$, let the observed eigenvalues be $9$ and $1$, and use the same average-target rule with $\lambda=1/4$. Here $\bar l=5$, so the estimate eigenvalues are $d_1=(3/4)9+(1/4)5=8$ and $d_2=(3/4)1+(1/4)5=2$.

The trace is $\operatorname{tr}(S^{-1}\delta)=8/9+2=26/9$. The divergence formula gives $D_\delta=3-2(1/4)=2.5$, and the determinant ratio is $|S|/|\delta|=9/(8\cdot2)=9/16$, so its logarithm is approximately $-0.575364$.

The constant is $c_{10,2}=\psi(5)+\psi(4.5)+2\log2-2\log10\approx-0.323887$. Consequently, \[ \widehat R=(7/10)(26/9)+(2/10)(2.5)-0.575364-(-0.323887)-2 \approx0.270745. \] Equivalently the constant contribution $-c_{10,2}-p$ is approximately $-1.676113$. Both forms give the same estimate.

Using the Estimate and Avoiding Misinterpretation

To tune a smooth eigenvalue-shrinkage rule, compute its estimated risk at each candidate tuning value and choose the value with the smallest estimate. This procedure compares the rules under the specified Gaussian model and Stein loss. It does not make the estimate a pointwise loss, nor does it guarantee low risk if the model or regularity conditions fail. As with any data-based tuning criterion, candidates should be chosen and interpreted in light of the estimation problem.

Do not confuse loss and risk. The loss depends on the unknown $\Sigma$; the displayed estimate is random and unbiased for its expectation over repeated samples.
Do not omit the divergence. The estimator's sensitivity to $S$ contributes to risk. Omitting it generally destroys unbiasedness.
Do not change the loss silently. This formula estimates Stein loss, not Frobenius loss or held-out likelihood.
Check positivity and regularity. A rule producing a zero or negative estimate eigenvalue has undefined log determinant and is outside this result.
Check the sampling convention. The coefficients above use $S=n^{-1}\sum X_kX_k^\mathsf{T}$ and a known zero mean. Changing the divisor or estimating the mean requires corresponding changes.

Key takeaway: Under Gaussian sampling, Stein–Haff integration by parts estimates the trace component of Stein loss, while the Wishart log-determinant calibration estimates its unknown normalization. Together they provide a computable, unbiased risk criterion for regular eigenvalue-shrinkage rules.

Check Your Understanding

Use the definitions and conditions above to answer the following.

Why must the estimator eigenvalues be positive for the displayed Stein risk estimate?
For the identity rule in dimension $p$, show how the two parts of $D_\delta$ add to $p(p+1)/2$.
In the first worked example, which of the three candidates has the smallest estimated risk, and what limitation applies to that conclusion?
What quantity does $c_{n,p}$ calibrate, and why does it enter with a minus sign in the risk estimate?
Why is it incorrect to describe $\widehat R_{\mathrm{Stein}}(S;\delta)$ as the realized loss for the observed sample?