{'fraction_non_alphanumeric': 0.05895646060362577, 'fraction_numerical': 0.01983628785215688, 'mean_word_length': 3.9085531180675375, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
James-Stein estimation of the first principal component

The Stein paradox has played an influential role in the field of high dimensional statistics. This result warns that the sample mean, classically regarded as the"usual estimator", may be suboptimal in high dimensions. The development of the James-Stein estimator, that addresses this paradox, has by now inspired a large literature on the theme of"shrinkage"in statistics. In this direction, we develop a James-Stein type estimator for the first principal component of a high dimension and low sample size data set. This estimator shrinks the usual estimator, an eigenvector of a sample covariance matrix under a spiked covariance model, and yields superior asymptotic guarantees. Our derivation draws a close connection to the original James-Stein formula so that the motivation and recipe for shrinkage is intuited in a natural way.

1. Introduction. The Stein paradox has played an influential role in the field of high dimensional statistics. This result warns that the sample mean, classically regarded as the "usual estimator", may be suboptimal in high dimensions. In particular, Stein (1956) showed that the usual estimator of a location parameter θ ∈ R p from uncorrelated Gaussian observations becomes inadmissible when p > 2 under a mean-squared error criterion. That is, an estimator with a uniformly lower risk must exist. That estimator was established by James & Stein (1961) and eponymously named.
Among the numerous perspectives that motivate the James-Stein estimator, 1 the empirical Bayes perspective (see Efron & Morris (1975)) is particularly elegant. Letting η denote the sample mean computed with n measurements of an unknown θ ∈ R and assuming an additive, normally distributed error w that has a zero mean and a variance ν 2 (e.g., ν = δ/ √ n where each measurement has standard error δ), we write η = θ + w. (1) Taking a Gaussian prior on the unknown θ, that is independent of w, implies that the bivariate-normal conditional expectation formula. While, by definition of conditional expectation, E(θ | η) is the best estimator of θ in the sense of mean-squared error, it cannot be implemented directly as the first two moments of η are unknown. 2 Stein's paradox now amounts to the fact that "good" substitutes for E(η) and Var(η) are available only in higher dimensions; precisely, when θ ∈ R p with p > 2.
Formula (2) extends easily to the multivariate 3 case and motivates the estimator where m is an estimate (or guess) of the expected value of η ∈ R p and c ∈ (0, 1) is a shrinkage parameter. In words, (3) attempts to center the 1 A few examples include the Galtonian regression perspective promoted by Stigler (1990), the purely frequentist development of the estimator in Gupta & Peña (1991) and the geometrical explanation in Brown & Zhao (2012) that builds on Stein's original heuristic argument (Stein 1956, Section 1). 2 Often, ν is assumed to be known but estimatesν can also be used as done in James & Stein (1961).
3 See formulas in Anderson (2003, Section 2.5) which can be used to design James-Stein estimators of an unkown θ ∈ R p when w has a general covariance Σ as in Bock (1975). entries of η, shrinks the resulting entries and recenters at m. Assuming ν is known and p > 2, setting where s 2 (η) = p i=1 (η i -m i ) 2 /p yields the James-Stein estimator. Remarkably, any fixed m ∈ R p (e.g., Stein (1956) considers the origin) 4 results in an estimator (3) with a strictly smaller mean-squared error than η (Efron & Morris 1975, Section 1). While three is provably the critical dimension 5 , Stein (1956, Section 1) heuristically argued for higher performance in higher dimensions and for relaxing the normality.
The James-Stein estimator has inspired a rich literature on the theme of "shrinkage" in statistics. Just a small sampling of examples includes ridge regression (Hoerl & Kennard 1970), the LASSO (Tibshirani 1996), the Ledoit-Wolf covariance estimator (Ledoit & Wolf 2004) and the Elastic Net (Zou & Hastie 2005). Excellent textbook treatments of the ideas behind the Stein paradox and James-Stein shrinkage include Gruber (2017) and Fourdrinier, Strawderman & Wells (2018). In this paper, we leverage these ideas to develop and analyze a James-Stein estimator for the first principal component of a sample covariance matrix. The results again prove the efficacy of James-Stein estimation, and do so for one the cornerstone methods in high-dimensional statistics, principal component analysis (Jolliffe & Cadima 2016).
Consider a p × p sample covariance matrix S that is based on n observations of some random vector y ∈ R p . Without loss of generality, we write S = s 2 p hh + G (5) for G = Ss 2 p hh and h, the sample eigenvector with the largest eigenvalue, i.e., Sh = s 2 p h and s 2 p = max |z|=1 z, Sz .
By convention h has unit length and corresponds to a direction along which the variance of S (i.e., s 2 p ) is maximum (i.e., the first principal component). A substantial and rapidly growing literature exists to study the (p and/or n) asymptotic behavior of the eigenpair (s 2 p , h) and the remaining eigenstructure in order to quantify either the estimation or the 4 A natural choice for m takes the sample mean of η in each entry, referred to as shrinkage toward the "grand mean". However, this would require p > 3 (Efron & Morris 1975)). 5 The theme of a critical dimension is encountered frequently in statistics and probability. Brown (1971), for example, derives a close mathematical relationship between the admissibility of the James-Stein estimator and the transience of the Brownian motion in R p , which also requires p > 2. empirical error. See Wang & Fan (2017) for recent results and a systematic discussion of this literature. The topic of shrinkage estimators arises naturally in this context and was raised by Stein (1986), who suggested improving the usual estimate S via eigenvalue shrinkage. Indeed, there is by now a large literature on estimators that adjust the eigenvalues of sample covariance matrices to improve their performance with respect to some loss function (Donoho, Gavish & Johnstone 2018).
The vector h JS is the James-Stein estimator of the first principal component of the data. The numerator contains the shrinkage formula (3) while the divisor normalizes the shrunk vector to a unit length (by convention). The relationship to the shrinkage parameter c in (4) is evident by treating p as large. The estimateν 2 corresponds to the bulk of the eigenvalue spectrum, and may be viewed the "noise" in the context of a signal-to-noise ratio that plays a prominent role of the results in Sections 2 & 3.
It is reasonable to suspect that a James-Stein type shrinkage of the principal component h, a high dimensional vector, could improve either the convergence rate or accuracy of the limit in some appropriate asymptotic regime. However, the standard orthonormal transformation and the eigengap partition of the sample eigenvectors, that is typically leveraged 6 The HDLS framework, as introduced in Hall, Marron & Neeman (2005) and Ahn, Marron, Muller & Chi (2007), is increasingly relevant for data science (Aoshima, Shen, Shen, Yata, Zhou & Marron 2018). 7 Roughly speaking q is the number of factors (or spikes) in the data, after which a sufficiently large eigengap (between the qth and the next eigenvalue) is observed (see Fan, Guo & Zheng (2020)). 8 With a slight abuse of the notation, u -x = (u 1 -x, . . . , u p -x) for u ∈ R p and x ∈ R.
by their asymptotic analyses (e.g., Paul (2007), Shen, Shen, Zhu & Marron (2016) and Wang & Fan (2017)), can obscure the systematic nature of the sample bias. As sensibly pointed out by Wang & Fan (2017) in reference to the partition of the sample eigenvector, the "two parts intertwine in such a way that correction for the biases of estimating eigenvectors is almost impossible." However, in the original (untransformed) coordinate system and the HDLS asymptotic regime, the bias can in fact be identified, characterized and (partially) corrected. This program was carried out in Goldberg, Papanicalaou & Shkolnik (2021), who adopt a factor model in an HDLS regime and utilize a portfolio theory application to motivate their analysis.
The main results of this paper rederive the adjustment of Goldberg et al. (2021) but within a James-Stein type framework. In particular, we establish identity (1) in which η = s p h and θ is related to the associated population eigenvector. From here, the James-Stein shrinkage acts on the perturbation w so that the estimator h JS outperforms h on the mean-squared error and angle metrics as p ↑ ∞. The theoretical guarantees provided here are new and their proofs rely on a different set of mathematical tools than Goldberg et al. (2021). In particular, the new approach leverages Weyl's inequality and the Davis-Kahan theorem from matrix perturbation theory to give simpler proofs and potentially expand the scope of applicability of the resulting estimator. The HDLS regime, in which the number of variables p grows to infinity, and the number of observations n to be fixed, plays a crucial role in the analysis.
The paper is organized as follows. Section 2 defines the spiked covariance model underlying our results. Section 3 develops the James-Stein estimator h JS and Section 4 proves the theoretical guarantees for this estimator. Appendix A contains proofs of the auxiliary results. The following notation is used throughout. Let u, v denote the standard inner product of u, v ∈ R d so that |u| = u, u and m(u) = u, e /d where e = (1, . . . , 1) are the length and mean. Set s 2 (u) = |u -m(u)| 2 /d and cov(u, v) = u -m(u), v -m(v) /d (see the notation of footnote 8). We use a subscript 1 ≤ p ≤ ∞ to highlight the dependence on p of various quantities, e.g., m(η) = m p (η) for η ∈ R p and m ∞ (η) is the limit lim p↑∞ m p (η) when it exists.

2.
A scarcely sampled spiked model. We use a spiked covariance model, borrowed from the HDLS literature. We also restrict ourselves to a single unbounded spike in the "boundary case" (see Jung, Sen & Marron (2012)) wherein the largest eigenvalue of the covariance matrix grows linearly in p. In particular, consider a mean-zero (w.l.o.g.) random vector y ∈ R p with a p × p covariance matrix Σ = Var(y) and let, for a symmetric, positive-semidefinite p × p matrix Γ and vector β ∈ R p . The following affirms that b = β/|β| is an eigenvector of Σ with eigenvalue β, β .
To state our additional assumptions on the model, we project the data vector y ∈ R p onto the eigenvectors of Σ. More precisely, define It is immediate that E(ψ p ) = 0 and Var(ψ p ) = 1. For every eigenvalue α i of Γ, let Note, E(φ i ) = 0 and Var(φ i ) = α i . As the dimension p grows we obtain a sequence {φ i } i≥1 . As a technical remark, φ i = φ i,p and α i = α i,p depend on p. We consider a sequence of models (7) constructed from sequences {β i } i≥1 and {Γ p }.
Condition (i) imposes regularity on the sequence {β i } i≥1 and implies that the largest eigenvalue of Σ (i.e., β, β ) grows linearly with p. The random variable ψ ∞ in (ii) is closely related to a principal component score (in the limit p ↑ ∞), and it captures the randomness along the first (population) principal component. Conditions (iii) and (iv) are related to certain requirements on a measure of sphericity of the model (e.g., tr(Σ) 2 /(tr(Σ 2 )p)) and summarize the conclusions of Jung & Marron (2009, Theorem 1). In particular, (iii) may be viewed as laws of large numbers for {φ i } i≥1 and {ϕ i } i≥1 and suggests pm(γ i ) 1 whereas (i) implies the first principal component has √ pm(b) 1. This highlights that the spike eigenvector b differs from those of Γ both in terms of the magnitude of its eigenvalue as well as the structure of the vector itself. According to (iv), the average eigenvalue of Γ (unlike the spike eigenvalue β, β ) is bounded in p, and per (iii), the eigenvectors γ i do not have entries biased towards a nonzero mean √ pm(γ i ), unlike b (an exception is μ = 0, which is the case when the James-Stein estimator will turn out to be ineffective; see Remark 2.8).
The following assumption is a standard one in statistical data analysis but may be (and often is) relaxed in the HDLS setup (e.g., Jung et al. (2012)). Assumption 2.3. There are a fixed n ≥ 2 i.i.d. observations of y ∈ R p .
Our forthcoming results hold even when only 2 observations are available, hence, a scarcely sampled model. Let Y be the p × n data matrix with the kth column containing the kth observation of y, and define the sample covariance matrix S by S = YY /n.
We let h ∈ R p denote the eigenvector of S with the largest eigenvalue s 2 p (see (6)). It is unique only up to sign (and |h| = 1), motivating the following (c.f. Assumption 2.1).
We write S = G + ηη in analogy to (7) (taking G = Sηη ) and set η = s p h.
Next, define the following measure of finite-sample distortion, χ 2 n = X , X /n and X = lim The latter limit exists by Assumption 2.2 while Assumption 2.3 implies that X ∈ R n has i.i.d. entries (distributed as ψ ∞ ). Consequently, we have χ 2 n → 1 as n ↑ ∞. We can measure the error in any estimator η of θ by the mean-squared error, as would be consistent with the James-Stein framework.
For SNR, we regard σχ n as a distorted signal, and δ/ √ n as noise that vanishes as n grows (also, σχ n → σ). The signal-incoherence r ∞ is the limit of r p = r p (β) = s p (β)/|β| (per Assumption 2.2) determined by the signal-to-noise ratio μ/σ of the vector β. A large value of r ∞ corresponds to more variation, or "incoherence" in the entries {β i } i≥1 A more standard way to evaluate the goodness of a sample eigenvector h is via its angle away from its population counterpart b = β/|β|. To this end, let Proposition 2.7. Let θ = χ n β and suppose Assumptions 2.1-2.4 hold. Then, Remark 2.8. Clearly, (17) also holds with θ = β.

James-Stein estimation of sample eigenvectors.
Having established that the sample eigenvector corresponding to the spike (i.e., the largest eigenvalue) carries finite sample error, it is natural to ask whether James-Stein shrinkage can improve this "usual" estimator. The key to this question is the (to be established) identity η = θ + w and θ = χ n β (18) for a random vector w ∈ R p specified in (27) of Section 4.1. As in definition (11), we have η = s p h where h is the sample eigenvector with the largest eigenvalue, s 2 p . The perturbation w of θ turns out to be such that the shrinkage of η is effective. We remark that the recipe of Section 1 extends our derivation of the James-Stein estimator below to the case of multiple (there, q) spikes in a natural way. This extension is effective because the eigenvectors corresponding to the spikes are mutually orthogonal, but it is suboptimal. An optimal estimator in the multi-spiked setup is left for future work.
based on (4) withν, an estimate of the "noise". It is reasonable to assign the latter to be the average of the non-spiked, non-zero eigenvalues of S. That is,ν where the scaling by p turns out to be necessary due to the counterintuitive behavior of the HDLS asymptotics. When p < n the divisor n -1 must be replaced by p -1.
This paves the way for the James-Stein sample eigenvector estimate, of the unnormalized eigenvector and by convention, we take unit length version, as the James-Stein estimator of the population eigenvector b = β/|β|.
The following James-Stein type theorems characterize the improvement due to shrinkage in the original mean-squared sense as well as in the angle metrics.
where c ∞ ∈ (0, 1) is the limit of c p = c in (19) (22) and with SNR and r ∞ in (15), we have Related results may be found in Goldberg, Papanicolaou, Shkolnik & Ulucam (2020) and Goldberg et al. (2021) but the metrics there are motivated by solutions of certain quadratic programs that are useful in finance and portfolio theory.

The geometry of Stein's paradox.
We shed insight into the James-Stein estimator in (21) by deriving general conditions under which Theorems 3.1 and 3.2 hold. Our analysis adopts the pure frequentist perspective in Gupta & Peña (1991) and supplements it by illustrating the Euclidean and the spherical geometry of the estimator. The two geometries reflect the definitions of the error (MSE and SPH) in the two theorems.
We define the signal-to-noise ratio SNR and the signal-incoherence r ∞ as SNR = ξ ν and r ∞ = 1 which are compatible with (15) upon taking θ = χ n β and ν = δ/ √ n. The following result establishes the conclusions of Theorems 3.1 and 3.2 in our abstract setting.

Proofs the main results.
We proceed by the following three main steps.
(2) We derive the convergence properties in the HDLS regime of the eigenvalues and eigenvectors of S under our spiked covariance model setting in Section 4.2.
(3) We verify that θ, w in (1) andν in (20) satisfy the conditions of Assumption 3.5, which leads to the guarantees for James-Stein shrinkage in Proposition 3.7.
Theorems 3.1 and 3.2 and then corollaries of Proposition 3.7, which is proved in Appendix A. We will make use of two classic results in matrix perturbation theory.

Establishing the key identity.
A key tool for random matrix theory in the HDLS regime is the dual sample covariance matrix. This is n × n matrix (n ≥ 2 is fixed), The next result is well known and relates the spectra of S = YY /n and L.
Proof. Multiplying the identity Lu = 2 u by Y from both sides, we obtain YLu = 2 Yu ⇒ SYu = 2 p n Yu.
The spike model Σ = Γ + ββ has a full basis of eigenvectors given by b = β/|β| and {γ i } p-1 i=1 , the latter corresponding to the nonzero eigenvalues {α} p-1 i=1 of Γ. Thus, and ψ = ψ p = β, y / β, β as in (8). Consequently, letting Y denote the p × n matrix of i.i.d. observations of y ∈ R p , we have where X = Y β/ β, β ∈ R n consists of i.i.d. observations of ψ and E is a p × n matrix with i.i.d. columns consisting of the observations of ε as defined above.
By orthogonality, we obtain that L = Y Y/p = ( β, β /p) XX + EE /p. Let x p be the eigenvector of L with the largest eigenvalue, 2 p , and x ∞ = X /χ n ∈ R n (the unit length normalization of X ), where X and χ n are defined in (12). By Lemma 4.1, We deduce that η = s p h = χ n β + w as required by (18) where 4.2. Convergence of the eigen-values/vectors. It is not difficult to establish that limit of L = L (p) as p ↑ ∞ (in any norm on R n×n ) takes the following form L (∞) = (σ 2 + μ 2 )(nχ 2 n ) x ∞ x ∞ + δ 2 I. The first term is the limit of ( β, β /p)XX under Assumption 2.2 (note that β, β = s 2 (β) + m 2 (β)) and the definitions of X and x ∞ above. By Assumption 2.3, the columns of E are i.i.d. copies of ε with E(ε) = 0 p by definition, the strong law of large numbers, confirms the that the offdiagonal entries of the second term are zero. That δ 2 determines all the diagonal entries is again a consequence of Assumption 2.2. This entails proving that ε, ε /p converges almost surely to δ 2 as is done in Section 4.3 item (ii). The matrix L (∞) has an easily described spectrum. Its largest eigenvalue is given by 2 ∞ = (σ 2 + μ 2 )(nχ 2 n ) + δ 2 and has the eigenvector x ∞ . All remaining eigenvalues equal δ 2 . Since L (p) converges (in any norm) to L (∞) , the largest eigenvalue 2 p converges to 2 ∞ almost surely. All the remaining eigenvalues converge to δ 2 . By Weyl's inequality, setting A = Ω and B = (L -Ω) so that A + B = L, we have | 2 p-i+1 -((σ 2 + μ 2 )(nχ 2 n )1 {i=1} + δ 2 )| ≤ |L (p) -L (∞) | 1 ≤ i ≤ n.
We immediately deduce (by Lemma 4.1) the following result, variants of which appear in Jung et al. (2012), Shen et al. (2016) and Goldberg et al. (2021). Let s 2 p-i+1 denote the ith largest eigenvalue of S (for i > min(p, n) all eigenvalues are zero).
Note the condition x p , x ∞ ≥ 0 is without loss of generality as the orientation of the eigenvectors is always arbitrary. The following result follows immediately.