2.1. Oracle Linear Estimators

Let A be an arbitrary q × q matrix. This section considers the candidate linear estimator XA for M. The quadratic risk (3) of XA is

where B = A′. This risk is a strictly convex function of B. Expanded algebraically,

We find the matrix

that minimizes (15). Let

The matrix derivatives

(see Rao and Toutenberg, 1995, Sect. A.15) entail that

Setting this risk derivative equal to zero and simplifying yields the minimizing matrix

Applying (9) to (16) and (19) yields the expressions

with π_j as in (10). The matrix

is evidently a symmetric q × q matrix whose eigenvalues all lie in [0,1].

Thus, the linear estimator XA with smallest quadratic risk is

. It is an oracle estimator because its definition involves the unknown parametric function W and σ². By substitution of

into (14), the risk of the oracle linear estimator is seen to be (13). Section 3 constructs an adaptive estimator whose risk converges asymptotically to this desirable value.

2.2. Oracle Symmetric Affine Shrinkage Estimators

Let

denote the set of all q × q symmetric matrices whose eigenvalues lie in [0,1]. By the analysis of Section 2.1, it is reasonable to limit the search for low risk candidate linear estimators to the symmetric affine shrinkage class

. Let

denote those symmetric matrices whose eigenvalues are either 0 or 1. This is the class of orthogonal projections into R^q. For 0 ≤ k ≤ q, let

denote the orthogonal projections that have eigenvalue 1 with multiplicity k and eigenvalue 0 with multiplicity q − k. Evidently

For every

, the risk (14) of the linear estimator XA simplifies to

where

Of interest for subsequent developments are the following oracle estimators, obtained by minimizing risk over various closed subsets of

:

[bull ] The oracle TS estimator of M is

, where

[bull ] The oracle order k TLS estimator of M is

, where

[bull ] The oracle TLS estimator of M is

, where

The labels given to the oracle estimators are justified by their linkage to adaptive TLS and TS estimators in Sections 3 and 4. The next theorem provides explicit formulas for these oracle estimators and their risks. The following properties of the {π_j : 1 ≤ j ≤ q}, defined by (10), assist understanding of this theorem:

1. The {π_j} are nonincreasing in j and all lie in [0,1].
2. pσ²(1 − π_j)⁻¹ = pσ² + l_j².
3. p⁻¹|l_ju_jv_j′|²/σ² = (p⁻¹l_j²)/σ² = π_j/(1 − π_j) is the mean square signal-to-noise ratio of the jth summand in the singular value decomposition of M.
4. π_j > ½ if and only if the foregoing signal-to-noise ratio exceeds 1.
5. π_j = 0 if and only if the foregoing signal-to-noise ratio is 0.

THEOREM 2.1. The following expressions hold:

with

defined in (19), and

and

Proof. The first expression in (27) follows from (23) and (20) or from Section 2.1.

Consider symmetric matrices A that take the form

where each f_j ∈ [0,1]. For such A, (20) and the second equality of (23) imply

The second expression in (27) follows by applying (31) to the first expression in (27).

For every

, (23) gives

This and (25) entail

In this equation, rank(AW^1/2) ≤ k under the constraint rank(A) = k. By the Eckart–Young matrix approximation theorem and (20), the best approximation in norm of rank ≤ k to W^1/2 is

This quantity equals

when

is given by the first expression in (29), thereby justifying that formula. The second expression follows from (31) because

has the form (30).

Because

minimizes risk over the matrices

, finding it amounts to minimizing the risk (31) over {f_j} values that are either 0 or 1. The minimum risk is achieved by setting f_j = 1 when π_j > ½ and f_j = 0 otherwise. The expressions in (28) follow. In other words, the optimal choice of k in constructing the oracle TLS estimator is

. Note that the first

summands in the singular value decomposition of M are precisely those summands whose signal-to-noise ratio exceeds 1. This makes

an intuitively plausible thresholding estimator. █

3.1. Estimated Risk and Adaptation

Let

be an L₁-consistent estimator of σ². Two constructions of

are described in Section 4.1. It will be seen in proving Theorem 4.3 that lim_p→∞sup_{p⁻¹|M|²≤c} E|p⁻¹X′X − W − σ²I_q|² = 0 for every finite c > 0. This suggests estimating W by

Let

Applying singular value decomposition (5) yields

for

as in (7). The estimated risk of the candidate symmetric affine shrinkage estimator XA is defined to be

Corresponding to the oracle estimators discussed in Section 2.2 are the following adaptive estimators, obtained by minimizing estimated risk over various closed subsets of

:

[bull ] The adaptive TS estimator of M is

, where

[bull ] The adaptive order k TLS estimator of M is

, where

[bull ] The adaptive TLS estimator of M is

, where

The next theorem shows that these adaptive estimators coincide with TLS and TS estimators discussed in the Introduction and provides formulas for their estimated risks. The following properties of the

, defined by (7), assist understanding of this theorem:

1. The
   
   are nonincreasing in j with values that do not exceed 1 but can be negative, unlike the values of the {π_j}.
2. 
   .
3. 
   is an asymptotically unbiased estimator for p⁻¹l_j² as p tends to infinity (cf. the reasoning for Theorem 4.3 and, more generally, Van Huffel and Vandewalle (1991)).
4. Hence,
   
   estimates the mean square signal-to-noise ratio (p⁻¹l_j²)/σ² = p⁻¹|l_ju_jv_j′|² = π_j/(1 − π_j) of the jth summand in the singular value decomposition of M.
5. 
   if and only if the foregoing estimated signal-to-noise exceeds 1.
6. 
   if and only if the foregoing estimated signal-to-noise ratio is nonnegative.

Unlike its oracle counterpart

, the matrix

defined in (37) has eigenvalues less than or equal to 1 rather than in [0,1]. Define

Note that

,

, and

. The adaptive TS estimator involves

.

THEOREM 3.1. The following expressions hold:

and

and

Consequently,

Proof. Let

From the second equality in (38),

Observe that

For every

,

because A is positive semidefinite. Similarly,

. Moreover,

. Hence, (49) implies

On the other hand,

From (51) and (52),

, and the minimum is achieved when

. This proves the first expression in (43).

Consider symmetric matrices A that take the form

where each f_j ∈ [0,1]. For such A, (37) and the second equality in (38) imply

The second expression in (43) follows by applying (54).

For every

, the first equality in (38) and (47) imply

This and (40) entail

In this equation,

under the constraint rank(A) = k. By the Eckart–Young matrix approximation theorem and (47), the best approximation in norm of rank ≤ k to

is

This quantity equals

when

is given by the first expression in (45), thereby justifying that formula. The second expression in (45) follows from (54) because

has the form (53).

Because

minimizes estimated risk over the matrices

, finding it amounts to minimizing the estimated risk (54) over {f_j} values that are either 0 or 1. The minimum risk is achieved by setting f_j = 1 when

and f_j = 0 otherwise. The expressions in (44) follow. In other words, the optimal choice of k in constructing the adaptive TLS estimator is

. Note that the first

summands in the singular value decomposition of M are precisely those summands whose estimated signal-to-noise ratio exceeds 1. This makes

an intuitively plausible thresholding estimator.

The expressions (46) for

and

follow from their definitions and the already proved parts of Theorem 3.1. █

3.2. Total Shrinkage and the Efron–Morris Estimator

In the case of Gaussian errors, empirical Bayes estimators of M may be constructed as follows. Suppose that the conditional distribution of x_i given m_i is N(m_i,σ²I_q) and that the {x_i} are conditionally independent. Suppose that the {m_i} are independent and that the distribution of m_i is N(0,W). Under quadratic loss, the Bayes estimator of M is then the oracle estimator

.

Because

is a consistent estimator of

in the Bayes model, the adaptive TS estimator

derived in Theorem 3.1, is an empirical Bayes estimator of M. Another consistent estimator of

in the Bayes model is

, which need not belong to

for finite p. This generates the empirical Bayes estimator

Related estimators, albeit more complex, have been used by Green, Berman, Switzer, and Craig (1985) to analyze multiband satellite image data. The Efron and Morris (1976) estimator of M, assuming σ² known, is

The estimator

generalizes the positive-part James and Stein (1961) estimator to q-dimensional means, and the estimator

in (59) similarly generalizes the James–Stein estimator. The Efron–Morris estimator

is a refinement of

whose efficacy in reducing risk depends on the Gaussian error assumption. Under Assumptions A1 and A2 of Section 4, which do not assume Gaussian errors, the risks of

,

, and

converge as p tends to infinity, as do the estimators themselves in terms of the normalized loss metric. Each of these symmetric affine shrinkage estimators has a computationally stable expression in terms of the singular value decomposition of X. Moreover, the singular value representations reveal the relationship between these shrinkage estimators and the adaptive TLS estimator

.

4.1. Convergence of Loss, Risk, and Estimated Risk

Let

denote the set of all symmetric matrices A = {a_ij} such that max_i,j|a_ij| ≤ 1. Evidently

. The following assumptions support the asymptotic results:

A1. The error vectors {e_i : 1 ≤ i ≤ p} are independent, identically distributed, such that E(e_i) = 0, Cov(e_i) = σ²I_q, where σ² > 0 is unknown.

A2. Under the model of Assumption A1, the variance estimator

satisfies

for every finite c > 0.

THEOREM 4.1. Suppose that Assumptions A1 and A2 hold. Let T(A) denote either the loss L(XA,M) or the estimated risk

of candidate estimator XA. Then, for every finite c > 0 and every finite σ² > 0,

This theorem shows that the loss, risk, and estimated risk of candidate estimator XA converge asymptotically. The uniformity of this convergence over all

makes estimated risk a trustworthy surrogate for the true loss or risk.

THEOREM 4.2. Suppose that Assumptions A1 and A2 hold. Then, for every finite c > 0 and for every finite σ² > 0,

Moreover, for S equal to either the loss

or the risk

,

Replacing the subscript TS by TLS in these statements yields valid assertions about

.

By (63), the risk of the adaptive TS estimator

converges to the risk of the oracle TS estimator

, which achieves minimum risk over the class of symmetric affine shrinkage estimators

. By (64), the plug-in estimator

of the risk of

converges to the actual risk or loss of

. In this manner, we can gauge directly from the data how well adaptation has controlled risk. Theorem 3.1 gives a convenient, computationally stable expression for

. Analogous statements hold for the parts of Theorems 4.2 and 3.1 that concern the adaptive TLS estimator

.

Variance Estimation.

It remains to construct suitably consistent estimators of σ². When enough repeated observations are available on rows of M, the least squares estimator of σ² satisfies Assumption A2. In the absence of replication, we may consider the smallest singular value estimator

The next theorem shows that

is asymptotically biased upward in general and satisfies Assumption A2 under a restriction on the magnitude of the smallest singular value of M.

THEOREM 4.3. Suppose that Assumption A1 holds and that the q components of each e_i have finite fourth moments. Then, for every finite c > 0 and for every finite σ² > 0,

Let {ε_p ≥ 0 : p ≥ 1} be any sequence such that lim_p→∞ ε_p = 0. Then

Thus,

is a consistent estimator of σ² provided that l_q²/p tends to zero as p tends to infinity. If

is used to construct the adaptive TS and TLS estimators, then Theorems 4.1 and 4.2 continue to hold with one modification: the supremum over all M such that p⁻¹|M|² ≤ c must be replaced by the supremum over M such that p⁻¹|M|² ≤ c, p⁻¹l_q² ≤ ε_p. This recognizes that (67) holds for

rather than Assumption A2.

When σ² is estimated by

, then the

in (7) that defines the adaptive TLS estimator simplifies elegantly to

. The interpretation of

remains as described in the proof of Theorem 3.1: the first

summands in the singular value decomposition of M are precisely those summands whose estimated signal-to-noise ratio exceeds 1.

4.2. Numerical Experiment

The pertinence of the asymptotic theory to finite sample sizes can be checked on artificial data. We report one experiment with

and

. Let {u_i : 1 ≤ i ≤ 4} be independent 100 × 1 random vectors, the elements of each being independent uniform [0,1] random variables. Set

Let M be the result of rounding each entry in the matrix [c₁,c₂,c₃,c₄,c₅] to the nearest integer. Let E be a 100 × 5 matrix of independent normal(0,σ²) random variables, with σ² = 25, each rounded to the nearest integer. Set X = M + E. This construction satisfies Assumption A1 with p = 100 and q = 5.

For an instance of such artificial data constructed in S-Plus (code available from the author), the variance estimate

approximates σ² = 25, the vectors

are close, and

The estimated singular values exhibit an upward bias that is caused by the errors E. Theorem 4.3 explains this for the smallest singular value, and similar analysis is available for the others.

Because the data have been constructed, it is possible to compute oracle estimators and to compare their loss and risk with that of the unbiased estimator X. In particular, using risk expressions in Theorem 2.1,

The rank of the oracle TLS estimator is 3. As expected from limits in Theorem 4.2, the loss and risk are close to one another for each oracle estimator and are smaller for the oracle TS estimator than for the oracle TLS estimator. Both the oracle TLS estimator and the oracle TS estimator dominate the unbiased estimator X substantially.

Using estimated risk expressions in Theorem 3.1, the loss and estimated risk of the adaptive estimators are

The rank of the adaptive TLS estimator is 3, like that of its oracle counterpart. The loss and estimated risk approximate one another for each adaptive estimator and are smaller for the adaptive TS estimator than for the adaptive TLS estimator. Moreover, the loss of the adaptive TS or TLS estimator is only slightly larger than the loss of its oracle counterpart. These findings are again consistent with limits in Theorem 4.2. Both the adaptive TLS estimator and the adaptive TS estimator dominate the unbiased estimator X substantially.

4.3. Proofs

We first state three propositions that will be used in proving Theorem 4.1. Let the {λ_i(B)} denote the eigenvalues of any symmetric matrix B and let λ_max(B) denote the largest eigenvalue. The spectral norm of any matrix B, symmetric or not, is defined by

Let l_max(B) denote the largest singular value of B.

PROPOSITION 4.4.

These properties follow readily from the definitions given.

PROPOSITION 4.5. Let {Y_p : p ≥ 1} be random elements of C([−1,1]^d), the set of all continuous real-valued functions on [−1,1]^d. Let plim_p→∞ denote the limit in probability as p → ∞. Suppose that

and that

Then

The assumptions of Proposition 4.5 imply the weak convergence in C([−1,1]^d) of {Y_p : p ≥ 1} to the zero element. From this, (76) follows. See Wichura (1971) for a short proof of the weak convergence.

PROPOSITION 4.6. Suppose that {y_p : p ≥ 1}, y, {z_p : p ≥ 1}, and z are nonnegative random variables such that plim_p→∞y_p = y, y_p ≤ z_p a.s. for every p, Ez < ∞, and lim_p→∞ E|z_p − z| = 0. Then lim_p→∞ E|y_p − y| = 0.

Indeed, the conditions on the {z_p} and z imply that the {z_p} are uniformly integrable. Hence, so are the {y_p}. The result follows (cf. Neveu, 1965, p. 52).

Proof of Theorem 4.1. For any symmetric q × q matrix A = {a_ij}, let a = {{a_ij : 1 ≤ i ≤ j},1 ≤ j ≤ q} denote the diagonal and subdiagonal elements of A, organized as a column vector of dimension d = q(q + 1)/2. The function A(a) that maps a ∈ [−1,1]^d to

is one-to-one and onto. In the proof, we identify A with A(a) and treat the loss, risk, and estimated risk of a candidate estimator XA as functions of a. For convenient notation, let r(a) = r(A,W,σ²) as defined in (23), let

as defined in (38), and let T(a) = T(A) as in the theorem statement.

The strategy is to show that T(a) − r(a) converges in probability to zero for every a ∈ [−1,1]^d, and then use Proposition 4.5 to show that sup_{a∈[−1,1]^k}|T(a) − r(a)| converges in probability to zero, and finally invoke Proposition 4.6 to establish (62). Repeatedly used are the constraint p⁻¹|M|² ≤ c and the following properties of B = B(a) = A² and B = B(a) = (I_q − A)². For i ≤ j, the operator ∇_ij denotes partial differentiation with respect to the element a_ij of a.

These bounds use Proposition 4.4 and the identity ∇_ijB = ∇_ijA·A + A·∇_ijA.

We first prove the case

of (62). Define

by replacing

with σ² in the definition (38) of

. Because of the inequality

and Assumption A2, we may replace

by r(a) in the subsequent argument.

Pointwise Consistency. Let

By the law of large numbers,

. Moreover,

and

Thus, for every a ∈ [−1,1]^d,

Uniform Consistency. For any a, b in [−1,1]^d,

where

for some

on the line segment that joins a and b. Thus

Moreover, using Proposition 4.4,

Applying Markov's inequality and (78) to the right-hand side of (85) establishes existence of a finite constant C, not depending on p, such that

Hence,

Limits (82) and (87) plus Proposition 4.5 yield

L₁ Uniform Consistency. Let y_p = sup_{a∈[−1,1]^d}|Y_p(a)|. Using (80),

Let

. It follows by the Cauchy–Schwarz inequality and (77) that

By the law of large numbers, w_p converges in probability to its expectation qσ². By Vitali's theorem,

Let z be the constant obtained from z_p, the right-hand side of (90), by replacing each instance of w_p with qσ². It follows from (91) that lim_p→∞ E|z_p − z| = 0. This convergence, inequality (90), and Proposition 4.6 imply that (88) can be strengthened to

This completes the proof of (62) when

.

The argument for the case T(A) = L(XA,M) of (62) is similar. The loss of XA is

where D = B − A. Let Y_p(a) = L(XA,M) − r(a). Then

Because the right-hand side of (94) has the same structure as the last expression in (80), an argument parallel to the one that follows (80) completes the proof. Indeed, using (77) and (78),

█

Proof of Theorem 4.2. We show that (62) implies

where Z can be

or

or

. The three limits to be proved in (63) and (64) are immediate consequences of (96).

First, (62) with

entails

Hence, (96) holds for

and

Second, (62) with T(A) = L(XA,M) gives

These limits together with (98) establish the remaining two cases of (96). The argument continues to hold when the subscript TS is replaced by TLS. █

Proof of Theorem 4.3. First we show that

Write M = {m_ij} and E = {e_ij}. Let δ_jk denote the Kronecker delta. Then,

Under Assumption A1 and the finite fourth moment condition on the components of each e_i,

In addition,

Limit (100) now follows.

Next,

the last inequality using parts (iv) and (vi) of Proposition 4.4. This inequality plus (100) implies limit (66) and hence also (67). █