Stein Risk Estimates for Series Shrinkage

Why Estimate Risk from the Observed Coefficients?

A series estimator may depend on a tuning choice: how much to shrink the estimated coefficients, or how many coefficients to retain. The ideal choice minimizes risk, but risk depends on the unknown true coefficients. Stein’s unbiased risk estimate (SURE) replaces that unknown risk with a quantity computed from the observed coefficients and the known noise variance.

The distinction between a fixed rule and a rule chosen from the data is essential. SURE is unbiased for each fixed estimator satisfying the conditions below. Minimizing several such estimates is a useful tuning procedure, but the minimum is not automatically an unbiased estimate of the risk of the selected estimator.

Gaussian Series-Coefficient Model and SURE

Consider a finite set of $m$ series coefficients observed in the Gaussian sequence model

Y=\theta+\varepsilon,\qquad \varepsilon\sim N_m(0,\sigma^2 I_m),

where $\theta\in\mathbb{R}^m$ is fixed and $\sigma^2>0$ is known. Thus the coordinates $Y_i$ are independent with $Y_i\sim N(\theta_i,\sigma^2)$. For an estimator $\delta(Y)\in\mathbb{R}^m$, use squared coefficient error

\mathcal R(\theta,\delta)=E_\theta\!\left[\|\delta(Y)-\theta\|^2\right].

This is coefficient-space risk. For an orthonormal series representation, it is the squared $L^2$ error contributed by these $m$ coefficients. Write an estimator as $\delta(y)=y+g(y)$, where $g$ is its adjustment to the observed vector. The divergence is $\operatorname{div}\delta(y)=\sum_{i=1}^m \partial\delta_i(y)/\partial y_i$.

Definition: For a fixed estimator rule $\delta(y)=y+g(y)$, its Stein unbiased risk estimate is $$ \widehat{\mathcal R}_{\mathrm{SURE}}(Y) =\|g(Y)\|^2+2\sigma^2\operatorname{div}\delta(Y)-m\sigma^2. $$ Equivalently, it is $\|\delta(Y)-Y\|^2+2\sigma^2\operatorname{div}\delta(Y)-m\sigma^2$. It is an ordinary unbiased estimate of risk when the regularity and finite-expectation conditions stated below hold.

Stein’s Identity and the Unbiasedness Result

The key identity is a Gaussian integration-by-parts formula. The conditions matter: differentiability alone does not ensure that the risk and the proposed estimate have finite expectations.

Theorem: Suppose $Y\sim N_m(\theta,\sigma^2I_m)$. For each coordinate $i$, assume $g_i$ is absolutely continuous as a function of $y_i$ for almost every fixed value of the other coordinates; $E_\theta[|\partial g_i(Y)/\partial y_i|]<\infty$; and the Gaussian integration-by-parts boundary term vanishes. A sufficient boundary condition is that, for almost every fixed $y_{-i}$, $g_i(y)$ times the $N(\theta_i,\sigma^2)$ density in coordinate $y_i$ tends to zero at both ends of the real line. Also assume $E_\theta[\|g(Y)\|^2]<\infty$. Then the SURE above is integrable and $$ E_\theta[\widehat{\mathcal R}_{\mathrm{SURE}}(Y)] =\mathcal R(\theta,\delta). $$

Proof. Fix $i$ and condition on the other coordinates. Let $\phi_i$ denote the $N(\theta_i,\sigma^2)$ density. Since $\partial\phi_i(y_i)/\partial y_i=-(y_i-\theta_i)\phi_i(y_i)/\sigma^2$, integration by parts gives

\int (y_i-\theta_i)g_i(y)\phi_i(y_i)\,dy_i =\sigma^2\int \frac{\partial g_i(y)}{\partial y_i}\phi_i(y_i)\,dy_i.

The boundary term is zero by assumption. The integrability assumptions allow these conditional integrals to be integrated over $Y_{-i}$. In particular, $E_\theta[|g_i(Y)(Y_i-\theta_i)|]<\infty$ follows from Cauchy–Schwarz, because both $g_i(Y)$ and $Y_i-\theta_i$ have finite second moments. Summing the coordinate identities yields

E_\theta[g(Y)^\mathsf{T}(Y-\theta)] =\sigma^2 E_\theta\!\left[\sum_{i=1}^m \frac{\partial g_i(Y)}{\partial y_i}\right].

Now expand the squared error. Since $\delta(Y)-\theta=(Y-\theta)+g(Y)$,

\mathcal R(\theta,\delta) =m\sigma^2+E_\theta[\|g(Y)\|^2] +2E_\theta[g(Y)^\mathsf{T}(Y-\theta)].

Substitute the identity just proved. Because $\operatorname{div}\delta=m+\operatorname{div}g$, the result becomes

\mathcal R(\theta,\delta) =E_\theta\!\left[\|g(Y)\|^2+ 2\sigma^2\operatorname{div}\delta(Y)-m\sigma^2\right].

Finally, $E_\theta[\|g(Y)\|^2]<\infty$ and the assumed integrability of the derivatives imply that the SURE has finite absolute expectation. The risk is finite because $Y-\theta$ and $g(Y)$ have finite second moments. Thus the equality is ordinary finite-mean unbiasedness, not merely an identity involving infinite expectations. $\square$

Fixed Linear Shrinkage

For fixed weights $a_1,\ldots,a_m$, consider $\delta_i(y)=a_i y_i$. The weights are fixed constants here, not functions of the observed vector. Since $\partial\delta_i/\partial y_i=a_i$, substitution into the definition gives the following result.

Result: For fixed $a_i$, the SURE for coordinatewise linear shrinkage is $$ \widehat{\mathcal R}_{\mathrm{SURE}}(Y) =\sum_{i=1}^m\left\{(1-a_i)^2Y_i^2+2\sigma^2a_i-\sigma^2\right\}. $$ Its risk is $$ \mathcal R(\theta,\delta) =\sum_{i=1}^m\left\{(1-a_i)^2\theta_i^2+\sigma^2a_i^2\right\}. $$

Derivation of the risk formula. For each coordinate, $\delta_i(Y)-\theta_i=a_i(Y_i-\theta_i)+(a_i-1)\theta_i$. The two terms have zero covariance because $E_\theta[Y_i-\theta_i]=0$, so the expected squared error is $a_i^2\sigma^2+(1-a_i)^2\theta_i^2$. Summing gives the displayed risk. For SURE, $g_i(y)=(a_i-1)y_i$, so $\|g(Y)\|^2=\sum_i(1-a_i)^2Y_i^2$ and $\operatorname{div}\delta=\sum_i a_i$. Substitution gives the stated formula. Fixed linear rules satisfy the theorem’s conditions: their adjustments and derivatives are integrable, and the Gaussian boundary terms vanish.

Worked Example: SURE for a Fixed Shrinkage Rule

Take $m=3$, $\sigma^2=1$, $Y=(2,-1,3)$, and use the fixed common weight $a_i=1/2$. First calculate $\sum_iY_i^2=4+1+9=14$. The residual term is $\sum_i(1-a_i)^2Y_i^2=(1/4)(14)=3.5$, while $2\sigma^2\operatorname{div}\delta-m\sigma^2=2(1)(3/2)-3=0$. Hence the observed SURE is $3.5$.

For comparison, if the true vector were $\theta=(1,0,2)$, the rule’s risk would be $(1/4)(1^2+0^2+2^2)+3(1)(1/2)^2=1.25+0.75=2$. The observed SURE need not equal the risk for this particular observation. Unbiasedness concerns its expectation over repeated observations at the fixed $\theta$, not equality for each data set.

Using SURE to Choose a Shrinkage Level

Suppose all $m$ coefficients use a common fixed shrinkage factor $a$. For $Q=\sum_iY_i^2$, the SURE as a function of a candidate $a$ is

\widehat{\mathcal R}_{\mathrm{SURE}}(a) =(1-a)^2Q+2m\sigma^2a-m\sigma^2.

For a fixed candidate $a$, this is an unbiased risk estimate under the theorem’s conditions. Minimizing the displayed quadratic over $0\leq a\leq1$ gives a data-based choice. If $Q>0$, its derivative is $-2(1-a)Q+2m\sigma^2$, so the unconstrained minimizer is $1-m\sigma^2/Q$. Restricting to the interval gives

\hat a=\max\left\{0,\,1-\frac{m\sigma^2}{Q}\right\},

with $\hat a=0$ when $Q=0$. The minimization uses the observed coefficients, so the resulting procedure is data-dependent. The formula is a tuning rule; the fixed-$a$ unbiasedness statement alone does not prove that the minimized SURE is an unbiased estimate of the selected rule’s risk.

Worked Example: Choosing a Common Shrinkage Factor

Let $m=2$, $\sigma^2=1$, and $Y=(2,1)$, so $Q=2^2+1^2=5$. The candidate SURE is $\widehat{\mathcal R}_{\mathrm{SURE}}(a)=5(1-a)^2+4a-2=5a^2-6a+3$. Its derivative is $10a-6$, which is zero at $a=0.6$; the second derivative is $10>0$, so this is the minimum on $[0,1]$. The minimum SURE is $5(0.6)^2-6(0.6)+3=1.2$.

The same answer follows from $1-m\sigma^2/Q=1-2/5=0.6$. This calculation selects a shrinkage factor from the observations. It does not turn the selected value $1.2$ into an unbiased estimate of the selected estimator’s risk without further analysis of the selection step.

Truncating a Series at a Candidate Level

For a fixed truncation level $K\in\{0,\ldots,m\}$, retain the first $K$ coefficients and set the rest to zero: $\delta_i^{(K)}(Y)=Y_i$ for $i\leq K$, and $0$ for $i>K$. This is a fixed linear rule for each specified $K$. Its divergence is $K$, and its residual vector has entries zero for $i\leq K$ and $-Y_i$ for $i>K$. Therefore,

Result: The SURE for a fixed truncation level is $$ \widehat{\mathcal R}_{\mathrm{SURE}}(K) =\sum_{i>K}Y_i^2+2K\sigma^2-m\sigma^2. $$

Worked Example: Comparing Fixed Truncation Levels

Let $m=4$, $\sigma^2=1$, and $Y=(2,1.5,1,0.5)$. The squared observations are $4,2.25,1,0.25$, with total $7.5$. Calculate the SURE at each level:

For $K=0$, it is $7.5+0-4=3.5$. For $K=1$, it is $2.25+1+0.25+2-4=1.5$. For $K=2$, it is $1+0.25+4-4=1.25$. For $K=3$, it is $0.25+6-4=2.25$. For $K=4$, it is $0+8-4=4$.

Among these candidates, the smallest observed SURE is $1.25$, at $K=2$. Thus SURE selects the rule that keeps the first two coefficients. For every fixed $K$, the displayed criterion is unbiased for that rule’s risk. The act of choosing $K$ by minimizing the criteria is an additional data-dependent step.

What Unbiasedness Does—and Does Not—Say

For a finite collection of fixed candidate rules, SURE provides a comparable estimate of each candidate’s risk using the same observations. Choosing the candidate with the smallest value is a principled tuning method. But the minimum of unbiased estimates generally has a different expectation from the minimum of their target risks: taking a minimum favors candidates whose estimated risk happened to be unusually low. Consequently, candidate-by-candidate unbiasedness does not imply unbiasedness after selection.

If the tuning parameter is itself a function of $Y$, the selected estimator is a single data-dependent function. Its risk can sometimes be estimated by applying Stein’s formula directly to that full function, when it meets the required regularity and integrability conditions. In that calculation, the divergence must account for the data dependence of the tuning parameter. Simply inserting the selected parameter into a formula derived for a fixed parameter omits that issue. Discontinuous selection rules require particular care with the assumptions and boundary contributions.

Check finiteness: a valid integration-by-parts identity is not enough for ordinary unbiasedness. Require finite second moment of $g(Y)$, integrable derivatives, and the stated boundary condition.
Keep the rule fixed when using the basic formula: fixed weights and fixed truncation levels have the divergences used above.
Interpret the output correctly: a SURE value is an estimate of risk, not the realized squared error and not a guarantee about the selected rule’s risk.

Key takeaway: In a finite Gaussian series model, SURE follows from Stein’s identity and estimates the risk of a fixed sufficiently regular estimator when the relevant expectations are finite. It supplies a practical criterion for tuning shrinkage or truncation, but the data-driven selection step must not be mistaken for an automatic extension of fixed-rule unbiasedness.

Check Your Understanding

Use the Gaussian sequence model and conditions stated above.

For a fixed diagonal rule $\delta_i(Y)=a_iY_i$, identify its divergence and write its SURE.
Why does finite expected absolute derivative alone not guarantee ordinary finite-mean unbiasedness?
For $m=3$, what is the SURE for a fixed truncation level $K=1$?
Explain why minimizing SURE over candidate fixed rules does not by itself establish unbiasedness for the risk of the selected rule.
For common shrinkage on $[0,1]$, derive the minimizing value when $Q>m\sigma^2$, and state the value when $Q\leq m\sigma^2$.