Quadratic Shrinkage: Theory & Applications

• Quadratic shrinkage is a family of estimators that improve traditional methods by shrinking estimates toward a target to lower quadratic risk under quadratic loss.
• It encompasses techniques ranging from the James–Stein estimator and its positive-part version to higher-order polynomial shrinkage for improved risk performance.
• The framework extends to predictive densities, matrix estimation, and structured regression penalties, demonstrating versatility across various high-dimensional models.

Quadratic shrinkage denotes a family of estimators that modify a baseline estimator by shrinking it toward zero, a fixed target, a subspace, or a structured prior in order to reduce risk under a quadratic criterion. In the canonical multivariate normal mean model $X\sim N_p(\mu,\sigma^2I_p)$ with loss $L(\mu,\delta)=\|\delta-\mu\|^2$, the usual estimator $X$ is minimax but, for $p\ge3$, inadmissible; this observation motivates James–Stein rules, Bayes formulations based on superharmonic marginals, and a large body of extensions to predictive densities, matrix-valued estimators, quadratic-penalty methods, and high-dimensional models under general quadratic loss (George et al., 2012, Matsuda, 2022).

1. Classical normal-mean formulation

The classical setup observes

$X\sim N_p(\mu,\sigma^2I_p),$

and seeks an estimator $\delta(X)$ of $\mu$ under quadratic loss

$$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$

The usual estimator $\delta_0(X)=X$ has constant risk

$R(\mu,X)=p\sigma^2,$

so its maximum risk over $L(\mu,\delta)=\|\delta-\mu\|^2$0 equals $L(\mu,\delta)=\|\delta-\mu\|^2$1, and it is minimax (George et al., 2012).

For $L(\mu,\delta)=\|\delta-\mu\|^2$2, Stein (1956) showed that $L(\mu,\delta)=\|\delta-\mu\|^2$3 is inadmissible. James and Stein (1961) proposed

$L(\mu,\delta)=\|\delta-\mu\|^2$4

with risk

$L(\mu,\delta)=\|\delta-\mu\|^2$5

Hence $L(\mu,\delta)=\|\delta-\mu\|^2$6 is minimax and strictly better than $L(\mu,\delta)=\|\delta-\mu\|^2$7. If $L(\mu,\delta)=\|\delta-\mu\|^2$8 is unknown, one may replace $L(\mu,\delta)=\|\delta-\mu\|^2$9 by an estimate.

The essential feature is a data-dependent multiplicative factor smaller than $X$0. Risk reduction is obtained not by unbiasedness but by a bias–variance trade-off that lowers quadratic risk uniformly over $X$1.

2. Bayes representation and superharmonicity

The Bayes formulation makes shrinkage explicit. Under the conjugate Gaussian prior

$X$2

the posterior of $X$3 given $X$4 is

$X$5

and the Bayes rule is the posterior mean

$X$6

As $X$7, this tends to $X$8; as $X$9, it tends to $p\ge3$0. One may choose $p\ge3$1 to optimize Bayes risk or replace the shrinkage factor by its unbiased estimate in an empirical-Bayes construction (George et al., 2012).

For a general prior $p\ge3$2, possibly improper, let

$p\ge3$3

Brown (1971) and Stein (1974) showed that the Bayes-posterior mean can be written as

$p\ge3$4

Substituting this representation into the risk and using integration by parts yields

$p\ge3$5

A sufficient condition for minimaxity is either

$p\ge3$6

so that $p\ge3$7 is superharmonic, or the weaker condition

$p\ge3$8

This criterion organizes a broad class of shrinkage priors. The improper harmonic prior

$p\ge3$9

has marginal

$X\sim N_p(\mu,\sigma^2I_p),$0

which is superharmonic for $X\sim N_p(\mu,\sigma^2I_p),$1, so the corresponding Bayes rule $X\sim N_p(\mu,\sigma^2I_p),$2 dominates $X\sim N_p(\mu,\sigma^2I_p),$3. Strawderman’s proper priors $X\sim N_p(\mu,\sigma^2I_p),$4, obtained as scale-mixtures of Gaussians with

$X\sim N_p(\mu,\sigma^2I_p),$5

provide examples for which $X\sim N_p(\mu,\sigma^2I_p),$6 is superharmonic even though $X\sim N_p(\mu,\sigma^2I_p),$7 is not.

The same device yields multiple shrinkage estimators. Given targets or subspaces $X\sim N_p(\mu,\sigma^2I_p),$8, one forms recentered superharmonic marginals $X\sim N_p(\mu,\sigma^2I_p),$9 and then a mixture marginal

$\delta(X)$0

which remains superharmonic. The resulting Bayes rule

$\delta(X)$1

is minimax and adapts to whichever target is most supported by the data.

An empirical-Bayes analogue of the conjugate-normal rule uses the unbiased estimate $\delta(X)$2 of the shrinkage factor and truncates it at zero, producing the positive-part estimator

$\delta(X)$3

described as minimax and admissible.

3. Second-order and higher-order polynomial shrinkage

A different use of the adjective “quadratic” appears in polynomial shrinkage estimators of a multivariate normal mean. With

$\delta(X)$4

and balanced squared-error loss

$\delta(X)$5

one considers estimators of the form

$\delta(X)$6

The second-degree, or quadratic, shrinkage function is

$\delta(X)$7

so that

$\delta(X)$8

In the balanced-loss setting the optimal first- and second-order coefficients are

$\delta(X)$9

The estimator is minimax relative to the MLE $\mu$0 if and only if

$\mu$1

When $\mu$2, this reduces to

$\mu$3

(Benkhaled et al., 2021).

The explicit risk formula for the quadratic estimator is

$\mu$4

Choosing

$\mu$5

gives

$\mu$6

For $\mu$7, the quadratic rule therefore strictly dominates both the usual James–Stein estimator and the MLE.

The paper’s figures plot the risk ratios

$\mu$8

against $\mu$9 for $$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$0 and $$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$1; all curves lie below $$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$2, and each higher-order polynomial pushes the ratio down further, especially when $$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$3 and $$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$4 are small. The same Stein-lemma argument extends to third- and higher-order rules with

$$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$5

but one needs $$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$6 to control $$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$7 and guarantee domination. The practical gain diminishes because $$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$8 for large $$L(\mu,\delta)=\|\delta-\mu\|^2, \qquad R(\mu,\delta)=E_\mu\|\delta(X)-\mu\|^2.$$9 and high $\delta_0(X)=X$0, and the complexity and variance of estimating higher-order inverse moments become unattractive.

4. Predictive density, matrix, and singular-value generalizations

Quadratic shrinkage extends beyond point estimation of a single normal mean. In the normal linear model,

$\delta_0(X)=X$1

independent, the Bayes predictive density under $\delta_0(X)=X$2 can be written in the shrinkage form

$\delta_0(X)=X$3

where $\delta_0(X)=X$4 is the uniform-prior predictor. Diagonalizing the covariances reduces the problem to the canonical mean setting, and the same superharmonic-marginal conditions ensure minimaxity. In nonparametric regression, the Gaussian sequence model

$\delta_0(X)=X$5

with $\delta_0(X)=X$6 constrained in an ellipsoid, has minimax KL-risk asymptotically attained by Bayes predictive densities under Gaussian priors, in direct analogy with Pinsker’s theorem for $\delta_0(X)=X$7-risk (George et al., 2012).

A matrix-valued version arises when estimating an $\delta_0(X)=X$8 mean matrix $\delta_0(X)=X$9 from an $R(\mu,X)=p\sigma^2,$0 normal data matrix $R(\mu,X)=p\sigma^2,$1 under the matrix quadratic loss

$R(\mu,X)=p\sigma^2,$2

The usual MLE is $R(\mu,X)=p\sigma^2,$3, with risk

$R(\mu,X)=p\sigma^2,$4

Applying vector James–Stein shrinkage separately to each column gives

$R(\mu,X)=p\sigma^2,$5

or equivalently

$R(\mu,X)=p\sigma^2,$6

Abu-Shanab, Kent and Strawderman showed that if the columnwise shrinkers satisfy the vector cross-product inequality, then $R(\mu,X)=p\sigma^2,$7 strictly dominates $R(\mu,X)=p\sigma^2,$8 under the matrix loss for

$R(\mu,X)=p\sigma^2,$9

The same argument extends to positive-part James–Stein, empirical Bayes, subspace shrinkage, and superharmonic-prior rules, and also to unknown $L(\mu,\delta)=\|\delta-\mu\|^2$00 and correlated columns after whitening (Abu-Shanab et al., 2011).

General quadratic loss can also be encoded by a positive-definite matrix $L(\mu,\delta)=\|\delta-\mu\|^2$01. In the multivariate Gaussian sequence and matrix-denoising settings,

$L(\mu,\delta)=\|\delta-\mu\|^2$02

and

$L(\mu,\delta)=\|\delta-\mu\|^2$03

The Efron–Morris singular-value shrinker

$L(\mu,\delta)=\|\delta-\mu\|^2$04

shrinks each singular value and dominates the MLE under Frobenius loss. Under arbitrary $L(\mu,\delta)=\|\delta-\mu\|^2$05-loss it satisfies the oracle inequality

$L(\mu,\delta)=\|\delta-\mu\|^2$06

This leads to a multivariate Pinsker linear shrinker on Sobolev ellipsoids and to a blockwise Efron–Morris construction that is stated to be exactly adaptive minimax over multivariate Sobolev ellipsoids under the corresponding quadratic loss (Matsuda, 2022).

5. Quadratic penalties and structured regression shrinkage

Quadratic shrinkage may be implemented as a penalty term rather than as a multiplicative factor. In sparse Laplacian shrinkage for high-dimensional regression, one minimizes

$L(\mu,\delta)=\|\delta-\mu\|^2$07

where $L(\mu,\delta)=\|\delta-\mu\|^2$08 is the minimax concave penalty and $L(\mu,\delta)=\|\delta-\mu\|^2$09 is the graph Laplacian derived from a signed adjacency matrix. The quadratic form is

$L(\mu,\delta)=\|\delta-\mu\|^2$10

This term promotes local smoothness among coefficients associated with connected predictors; unlike a global ridge penalty, the shrinkage is local on the graph. Under the stated sub-Gaussian, restricted-eigenvalue, penalty-size, and minimum-signal conditions, the estimator is selection-consistent and equal to the oracle Laplacian shrinkage estimator with probability at least $L(\mu,\delta)=\|\delta-\mu\|^2$11. Proposition 5.2 further gives a generalized grouping property within graph cliques (Huang et al., 2011).

A related quadratic penalization appears in shrinkage for categorical regressors and group means. If $L(\mu,\delta)=\|\delta-\mu\|^2$12 denotes the sample mean in group $L(\mu,\delta)=\|\delta-\mu\|^2$13, shrinkage toward first-stage targets $L(\mu,\delta)=\|\delta-\mu\|^2$14 is obtained by minimizing

$L(\mu,\delta)=\|\delta-\mu\|^2$15

Writing $L(\mu,\delta)=\|\delta-\mu\|^2$16 and $L(\mu,\delta)=\|\delta-\mu\|^2$17, the unique minimizer is

$L(\mu,\delta)=\|\delta-\mu\|^2$18

Under the local-to-zero framework

$L(\mu,\delta)=\|\delta-\mu\|^2$19

shrinkage introduces asymptotic bias

$L(\mu,\delta)=\|\delta-\mu\|^2$20

but reduces variance. For every bounded $L(\mu,\delta)=\|\delta-\mu\|^2$21, provided $L(\mu,\delta)=\|\delta-\mu\|^2$22, the plug-in shrinkage estimator has asymptotic risk strictly below OLS uniformly over $L(\mu,\delta)=\|\delta-\mu\|^2$23. In Monte Carlo experiments with $L(\mu,\delta)=\|\delta-\mu\|^2$24 and $L(\mu,\delta)=\|\delta-\mu\|^2$25, the plug-in quadratic-shrinkage estimator never exceeds OLS risk and often improves by up to $L(\mu,\delta)=\|\delta-\mu\|^2$26 (Heiler et al., 2019).

6. High-dimensional, spectral, and heteroscedastic extensions

In high-dimensional mean estimation with unknown covariance, Wang et al. consider i.i.d. observations

$L(\mu,\delta)=\|\delta-\mu\|^2$27

with $L(\mu,\delta)=\|\delta-\mu\|^2$28, $L(\mu,\delta)=\|\delta-\mu\|^2$29, finite fourth moments, and no normality assumption. Accuracy is measured by

$L(\mu,\delta)=\|\delta-\mu\|^2$30

where $L(\mu,\delta)=\|\delta-\mu\|^2$31 is known positive definite. Their estimator

$L(\mu,\delta)=\|\delta-\mu\|^2$32

shrinks the sample mean toward $L(\mu,\delta)=\|\delta-\mu\|^2$33, with data-driven coefficients $L(\mu,\delta)=\|\delta-\mu\|^2$34 formed from four $L(\mu,\delta)=\|\delta-\mu\|^2$35-statistics $L(\mu,\delta)=\|\delta-\mu\|^2$36. Under Assumption 3.1 and the asymptotic condition $L(\mu,\delta)=\|\delta-\mu\|^2$37, one has

$L(\mu,\delta)=\|\delta-\mu\|^2$38

For the percentage relative improvement

$L(\mu,\delta)=\|\delta-\mu\|^2$39

Corollary 3.1 states that if $L(\mu,\delta)=\|\delta-\mu\|^2$40, then $L(\mu,\delta)=\|\delta-\mu\|^2$41; if $L(\mu,\delta)=\|\delta-\mu\|^2$42, then $L(\mu,\delta)=\|\delta-\mu\|^2$43 some $L(\mu,\delta)=\|\delta-\mu\|^2$44; and if $L(\mu,\delta)=\|\delta-\mu\|^2$45, then $L(\mu,\delta)=\|\delta-\mu\|^2$46. Simulations and leukemia microarray analysis are reported to show practical improvement over the sample mean and competing shrinkage rules, especially when $L(\mu,\delta)=\|\delta-\mu\|^2$47 is non-diagonal and training sizes are very small (Wang et al., 2012).

In spectral matrix estimation for partial coherencies, one observes

$L(\mu,\delta)=\|\delta-\mu\|^2$48

and a standard shrinkage estimator is

$L(\mu,\delta)=\|\delta-\mu\|^2$49

Schneider-Luftman and Walden study the quadratic-loss criterion

$L(\mu,\delta)=\|\delta-\mu\|^2$50

derive oracle rules QLa and QLb for spectral-matrix shrinkage, and also formulate precision-matrix shrinkage

$L(\mu,\delta)=\|\delta-\mu\|^2$51

A central methodological point is that Hilbert–Schmidt shrinkage can “over-shrink” and wipe out true large partial coherencies. In EEG-derived simulations, HS-based shrinkage can increase error when $L(\mu,\delta)=\|\delta-\mu\|^2$52 only slightly exceeds $L(\mu,\delta)=\|\delta-\mu\|^2$53, whereas oracle QLb gives the best PRISE, $L(\mu,\delta)=\|\delta-\mu\|^2$54 for $L(\mu,\delta)=\|\delta-\mu\|^2$55. Among the feasible methods, QLP-est retains almost all of the oracle QLP gain, $L(\mu,\delta)=\|\delta-\mu\|^2$56 at $L(\mu,\delta)=\|\delta-\mu\|^2$57 and $L(\mu,\delta)=\|\delta-\mu\|^2$58 at $L(\mu,\delta)=\|\delta-\mu\|^2$59, with much less variance across frequency; it is recommended as the full-estimation method (Schneider-Luftman et al., 2015).

A further generalization treats families of distributions with quadratic variance function

$L(\mu,\delta)=\|\delta-\mu\|^2$60

Estimators

$L(\mu,\delta)=\|\delta-\mu\|^2$61

then have oracle weight

$L(\mu,\delta)=\|\delta-\mu\|^2$62

or, when a sample-size factor $L(\mu,\delta)=\|\delta-\mu\|^2$63 appears,

$L(\mu,\delta)=\|\delta-\mu\|^2$64

Xie, Kou and Brown construct semiparametric URE-based shrinkage estimators

$L(\mu,\delta)=\|\delta-\mu\|^2$65

with monotonicity $L(\mu,\delta)=\|\delta-\mu\|^2$66 whenever $L(\mu,\delta)=\|\delta-\mu\|^2$67, and show that the URE rule is asymptotically first-order risk-optimal in its class. Their simulations report $L(\mu,\delta)=\|\delta-\mu\|^2$68–$L(\mu,\delta)=\|\delta-\mu\|^2$69 lower MSE than classical James–Stein and empirical-Bayes estimators when variances and means are correlated or heteroscedastic (Xie et al., 2016).

Across these formulations, quadratic shrinkage is unified less by a single estimator than by a common risk logic: a controlled introduction of bias to secure lower quadratic risk. The resulting theory is dimension-sensitive and model-sensitive. James–Stein dominance requires $L(\mu,\delta)=\|\delta-\mu\|^2$70; second-order polynomial improvement requires $L(\mu,\delta)=\|\delta-\mu\|^2$71; higher-order polynomial rules require $L(\mu,\delta)=\|\delta-\mu\|^2$72; and, in spectral estimation, methods based on inverse-power traces may be unstable even when their oracle versions are attractive. The main analytical tools are minimaxity, Bayes posterior means, superharmonicity, oracle inequalities, and explicit bias–variance decompositions.