Worked Comparison of Shrinkage and Thresholding

Why Compare These Estimators?

When the dimension is large relative to the sample size, a sample covariance matrix can be unstable even when the observations are Gaussian. Different regularization methods impose different structures: identity shrinkage pulls eigenvalues toward a common value, nonlinear eigenvalue shrinkage adjusts them separately, and entrywise thresholding removes selected off-diagonal entries. Their errors need not agree, because a loss function determines which features of the error matter most.

We compare these methods using two losses. The squared Frobenius loss adds squared errors across all entries, so it reflects aggregate matrix error. The operator-norm loss measures the largest absolute eigenvalue of the error for symmetric matrices, so it reflects the error in the most affected direction. Neither loss is universally preferable; the relevant choice depends on the purpose of the covariance estimate.

Three Estimators

Suppose \(X_1,\ldots,X_n\) are independent Gaussian vectors in \(\mathbb{R}^p\) with mean zero and covariance \(\Sigma\). With the mean known to be zero, take the sample covariance to be \(S=n^{-1}\sum_{k=1}^n X_kX_k^\top\). If the mean is estimated instead, the usual centered sample covariance can be used; the comparison below concerns the estimators applied to the resulting \(S\).

Identity shrinkage combines \(S\) with a scalar multiple of the identity matrix. For a chosen target scale \(\mu\) and weight \(\alpha\in[0,1]\), it is \(\widehat{\Sigma}_{\mathrm{Id}}(\alpha)=\alpha S+(1-\alpha)\mu I_p\). A common scale is \(\mu=\operatorname{tr}(S)/p\), which preserves the average sample eigenvalue. The single weight shrinks every eigenvalue toward the same value.

For nonlinear eigenvalue shrinkage, write the spectral decomposition \(S=Q\operatorname{diag}(\lambda_1,\ldots,\lambda_p)Q^\top\), with orthonormal eigenvectors in \(Q\). Replace each sample eigenvalue by a nonnegative adjusted value \(d_i\), giving \(\widehat{\Sigma}_{\mathrm{NL}}=Q\operatorname{diag}(d_1,\ldots,d_p)Q^\top\). The adjustment can depend nonlinearly on the sample eigenvalues and on the sampling regime. This definition describes the estimator's form; it does not prescribe one particular eigenvalue-adjustment rule.

Entrywise thresholding acts on covariance entries rather than on eigenvalues. With threshold \(\tau\geq0\) and the diagonal retained, define the hard-threshold estimator by \(\widehat{\Sigma}_{\mathrm{TH},ij}=S_{ij}\) if \(i=j\), and \(\widehat{\Sigma}_{\mathrm{TH},ij}=S_{ij}\mathbf{1}\{|S_{ij}|\geq\tau\}\) if \(i\ne j\). The separate diagonal rule matters: the threshold is applied only to off-diagonal entries.

Two Useful Loss Calculations

The first result identifies the best weight along a fixed identity-shrinkage path when the target and truth are treated as known. It explains precisely what “best weight” means under squared Frobenius loss; in practice, \(\Sigma\) is unknown, so the result is an oracle benchmark rather than a directly available tuning rule.

To prove the result when \(S\ne T\), put \(D=S-T\) and \(R=T-\Sigma\). The error is \(R+\alpha D\). Expanding the squared norm gives \(\|R+\alpha D\|_F^2=\|R\|_F^2+2\alpha\langle D,R\rangle_F+\alpha^2\|D\|_F^2\). Because \(\|D\|_F^2>0\), this is a strictly convex quadratic in \(\alpha\). Its derivative is \(2\langle D,R\rangle_F+2\alpha\|D\|_F^2\), which vanishes at \(\alpha=-\langle D,R\rangle_F/\|D\|_F^2=\langle D,\Sigma-T\rangle_F/\|D\|_F^2\). The constrained minimum on \([0,1]\) is that unconstrained value if it lies in the interval, and otherwise the nearer endpoint. This is exactly projection onto \([0,1]\). \(\square\)

The next result clarifies what eigenvalue adjustment can achieve if the sample eigenvectors are held fixed. It is an oracle statement: the diagonal quantities involving \(\Sigma\) are unknown in applications, and choosing \(d_i\) from the data requires an estimation method.

Orthogonal invariance of the Frobenius norm gives \(\|QDQ^\top-\Sigma\|_F^2=\|D-Q^\top\Sigma Q\|_F^2\). Write \(A=Q^\top\Sigma Q\). Expanding entry by entry, \(\|D-A\|_F^2=\sum_i(d_i-A_{ii})^2+\sum_{i\ne j}A_{ij}^2\). The second sum does not depend on \(D\), and each term in the first sum is minimized by \(d_i=A_{ii}=q_i^\top\Sigma q_i\). If \(d_i\geq0\) is required, the closest nonnegative number to \(A_{ii}\) is \(\max(A_{ii},0)\), which minimizes the corresponding squared term. \(\square\)

Reading the Comparison Carefully

The two losses weight error geometry differently. Frobenius loss accumulates discrepancies over all entries; operator norm is controlled by the largest error direction. Consequently, the ordering of two estimates can change with the loss, even when both estimates are positive definite. For example, take \(\Sigma=I_3\), \(A=\operatorname{diag}(1.4,1,1)\), and \(B=\operatorname{diag}(1.3,1.3,1)\). Both estimates are positive definite. Their squared Frobenius losses are \(0.4^2=0.16\) and \(0.3^2+0.3^2=0.18\), so \(A\) is better under that loss. Their operator-norm losses are \(0.4\) and \(0.3\), so \(B\) is better under operator norm.

Thresholding exploits approximate sparsity: if many true off-diagonal entries are zero or small, removing noisy sample entries can reduce aggregate error. Its risk reflects both the variance removed and the bias introduced by discarded nonzero entries. Shrinkage instead changes the spectrum or pulls entries toward a structured target, which can be preferable when the covariance is not sparse but sample eigenvalues are noisy. These are structural assumptions, not properties guaranteed by the Gaussian sampling model alone.