{'fraction_non_alphanumeric': 0.045751945370811495, 'fraction_numerical': 0.01842147054152771, 'mean_word_length': 4.206796758723334, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spectrum Estimation: A Unified Framework for Covariance Matrix Estimation and PCA in Large Dimensions

Covariance matrix estimation and principal component analysis (PCA) are two cornerstones of multivariate analysis. Classic textbook solutions perform poorly when the dimension of the data is of a magnitude similar to the sample size, or even larger. In such settings, there is a common remedy for both statistical problems: nonlinear shrinkage of the eigenvalues of the sample covariance matrix. The optimal nonlinear shrinkage formula depends on unknown population quantities and is thus not available. It is, however, possible to consistently estimate an oracle nonlinear shrinkage, which is motivated on asymptotic grounds. A key tool to this end is consistent estimation of the set of eigenvalues of the population covariance matrix (also known as the spectrum), an interesting and challenging problem in its own right. Extensive Monte Carlo simulations demonstrate that our methods have desirable finite-sample properties and outperform previous proposals.


Introduction
This paper tackles three important problems in multivariate statistics: 1) the estimation of the eigenvalues of the covariance matrix; 2) the estimation of the covariance matrix itself; and 3) principal component analysis (PCA). In many modern applications, the matrix dimension is not negligible with respect to the sample size, so textbook solutions based on classic (fixed-dimension) asymptotics are no longer appropriate. A better-suited framework is large-dimensional asymptotics, where the matrix dimension and the sample size go to infinity together, while their ratio -called the concentration -converges to a finite, nonzero limit. Under large-dimensional asymptotics, the sample covariance matrix is no longer consistent, and neither are its eigenvalues nor its eigenvectors.
One of the interesting features of large-dimensional asymptotics is that principal component analysis can no longer be conducted using covariance matrix eigenvalues. The variation explained by a principal component is not equal to the corresponding sample eigenvalue andperhaps more surprisingly -it is not equal to the corresponding population eigenvalue either. To the best of our knowledge, this fact has not been noticed before. The variation explained by a principal component is obtained instead by applying a nonlinear shrinkage formula to the corresponding sample eigenvalue. This nonlinear shrinkage formula depends on the unobservable population covariance matrix, but thankfully it can approximated by an oracle shrinkage formula which depends 'only' on the unobservable eigenvalues of the population covariance matrix. This is the connection with the first of the three problems mentioned above. Once we have a consistent estimator of the population eigenvalues, we can use it to derive a consistent estimator of the oracle shrinkage.
The connection with the second problem, the estimation of the whole covariance matrix, is that the nonlinear shrinkage formula that gives the variation explained by a principal component also yields the optimal rotation-equivariant estimator of the covariance matrix according to the Frobenius norm. Thus, if we can consistently estimate the population eigenvalues, and if we plug them into the oracle shrinkage formula, we can address the problems of PCA and covariance matrix estimation in a unified framework.
It needs to pointed out there that in a rotation-equivariant framework, consistent (or even improved) estimators of the population eigenvalues are not available; instead, one needs to retain the sample eigenvectors. As a consequence, consistent estimation of the population covariance matrix itself is not possible. Nevertheless, a rotation-equivariant estimator can still be very useful for practical applications, as evidenced by the popularity of the previous proposal of Ledoit and Wolf (2004). An alternative approach that would allow for consistent estimation of the population covariance matrix is to impose additional structure on the estimation problem, such as sparseness, a graph model, or an (approximate) factor model. But whether such structure does indeed exist is something that cannot be verified from the data. Therefore, at least in some applications, a structure-free approach will be preferred by applied researchers. This is the problem that we address, aiming to further improve upon Ledoit and Wolf (2004).
Of course estimating population eigenvalues consistently under large-dimensional asymp-totics is no trivial matter. Until recently, most researchers in the field even feared it might be impossible because deducing population eigenvalues from sample eigenvalues showed some symptoms of ill-posedness. This means that small estimation errors in the sample eigenvalues would be amplified by the specific mathematical structure of the asymptotic relationship between sample and population eigenvalues. But two recent articles by Mestre (2008) and El Karoui (2008) challenged this widely-held belief and gave some hope that it might be possible after all to estimate the population eigenvalues consistently. Still, a general satisfactory solution is not available to date. The work of Mestre (2008) only applies when the number of distinct population eigenvalues remains finite as matrix dimension goes to infinity. In practice, this means that the number of distinct eigenvalues must be negligible with respect to the total number of eigenvalues. As a further restriction, the number of distinct eigenvalues and their multiplicities must be known. The only unknown quantities to be estimated are the locations of the eigenvalues; of course, this is still a difficult task. Yao et al. (2012) propose a more general estimation procedure that does not require knowledge of the multiplicities, though it still requires knowledge of the number of distinct population eigenvalues. This setting is too restrictive for many applications.
The method developed by El Karoui (2008) allows for an arbitrary set of population eigenvalues, but does not appear to have good finite-sample properties. In fact, our simulations seem to indicate that this estimator is not even consistent; see Section 5.1.1.
The first contribution of the present paper is, therefore, to develop an estimator of the population eigenvalues that is consistent under large-dimensional asymptotics regardless of whether or not they are clustered, and that also performs well in finite sample. This is achieved through a more precise characterization of the asymptotic behavior of sample eigenvalues. Whereas existing results only specify how the eigenvalues behave on average, namely, how many fall in any given interval, we determine individual limits.
Our second contribution is to show how this consistent estimator of population eigenvalues can be used for improved estimation of the covariance matrix when the dimension is large compared to the sample size. This was already considered in Ledoit and Wolf (2012), but only in the limited setup where the dimension is smaller than the sample size. Thanks to the advances introduced in the present paper, we can also handle the more difficult case where the dimension exceeds the sample size and the sample covariance matrix is singular.
Our third and final contribution is to show how the same nonlinear shrinkage formula can be used to estimate the fraction of variation explained by a given collection of principal components in PCA, which is key in deciding how many principal components to retain.
The remainder of the paper is organized as follows. Section 2 presents our estimator of the eigenvalues of the population covariance matrix under large-dimensional asymptotics. Section 3 discusses covariance matrix estimation, and Section 4 principal component analysis. Section 5 studies finite-sample performance via Monte Carlo simulations. Section 6 provides a brief empirical application of PCA to stock return data. Section 7 concludes. The proofs of all mathematical results are collected in the appendix.

Large-Dimensional Asymptotics and Basic Framework
Let n denote the sample size and p . . = p(n) the number of variables. It is assumed that the ratio p/n converges as n → ∞ to a limit c ∈ (0, 1) ∪ (1, ∞) called the concentration. The case c = 1 is ruled out for technical reasons. We make the following assumptions.
(A1) The population covariance matrix Σ n is a nonrandom p-dimensional positive definite matrix.
(A2) X n is an n × p matrix of real independent and identically distributed (i.i.d.) random variables with zero mean, unit variance, and finite fourth moment. One only observes Y n . . = X n Σ 1/2 n , so neither X n nor Σ n are observed on their own.
(A3) τ n . . = (τ n,1 , . . . , τ n,p ) ′ denotes a system of eigenvalues of Σ n , sorted in increasing order, and (v n,1 , . . . , v n,p ) denotes an associated system of eigenvectors. The empirical distribution function (e.d.f.) of the population eigenvalues is defined as: where ½ denotes the indicator function of a set. H n is called the spectral distribution (function). It is assumed that H n converges weakly to a limit law H, called the limiting spectral distribution (function).
(A4) Supp(H), the support of H, is the union of a finite number of closed intervals, bounded away from zero and infinity. Furthermore, there exists a compact interval in (0, ∞) that contains Supp(H n ) for all n large enough.
Let λ n . . = (λ n,1 , . . . , λ n,p ) ′ denote a system of eigenvalues of the sample covariance matrix S n . . = n −1 Y ′ n Y n = n −1 Σ 1/2 n X ′ n X n Σ 1/2 n , sorted in increasing order, and let (u n,1 , . . . , u n,p ) denote an associated system of eigenvectors. The first subscript, n, may be omitted when no confusion is possible. The e.d.f. of the sample eigenvalues is defined as: ∀t ∈ R, F n (t) .
The literature on the eigenvalues of sample covariance matrices under large-dimensional asymptotics -also known as random matrix theory (RMT) literature -is based on a foundational result due to Marčenko and Pastur (1967). It has been strengthened and broadened by subsequent authors including Silverstein (1995), Silverstein and Bai (1995), Silverstein and Choi (1995), and Silverstein (1998, 1999), among others. These articles imply that there exists a limiting sample spectral distribution F such that In other words, the average number of sample eigenvalues falling in any given interval is known asymptotically.
In addition, the existing literature has unearthed important information about the limiting distribution F . Silverstein and Choi (1995) show that F is everywhere continuous except (potentially) at zero, and that the mass that F places at zero is given by Furthermore, there is a seminal equation relating F to H and c. Some additional notation is required to present this equation.
For any nondecreasing function G on the real line, m G denotes the Stieltjes transform of G: where C + denotes the half-plane of complex numbers with strictly positive imaginary part. The Stieltjes transform admits a well-known inversion formula: if G is continuous at a and b. Here, and in the remainder of the paper, we shall use the notations Re(z) and Im(z) for the real and imaginary parts, respectively, of a complex number z, so that The most elegant version of the equation relating F to H and c, due to Silverstein (1995), states that m . . = m F (z) is the unique solution in the set to the equation As explained, the Stieltjes transform of F , m F , is a function whose domain is the upper half of the complex plane. It can be extended to the real line, since Silverstein and Choi (1995) show that: ∀λ ∈ R \ {0}, lim z∈C + →λ m F (z) = . .m F (λ) exists. When c < 1,m F (0) also exists and F has a continuous derivative F ′ = π −1 Im [m F ] on all of R with F ′ ≡ 0 on (−∞, 0]. (One should remember that although the argument ofm F is real-valued now, the output of the function is still a complex number.) For purposes that will become apparent later, it is useful to reformulate equation (2.5). The limiting e.d.f. of the eigenvalues of n −1 Y ′ n Y n = n −1 Σ 1/2 n X ′ n X n Σ 1/2 n was defined as F . In addition, define the limiting e.d.f. of the eigenvalues of n −1 Y n Y ′ n = n −1 X n Σ n X ′ n as F ; note that the eigenvalues of n −1 Y ′ n Y n and n −1 Y n Y ′ n only differ by |n − p| zero eigenvalues. It then holds: (2.9) (Recall here that F has mass (c−1)/c at zero when c > 1, so that both F and F are nonnegative functions indeed for any value c > 0.) With this notation, equation (1.13) of Marčenko and Pastur (1967) reframes equation (2.5) as: for each z ∈ C + , m . . = m F (z) is the unique solution in C + to the equation (2.10) While in the case c < 1,m F (0) exists and F is continuously differentiable on all of R, as mentioned above, in the case c > 1,m F (0) exists and F is continuously differentiable on all of R.

Individual Behavior of Sample Eigenvalues: the QuEST Function
We introduce a nonrandom multivariate function called the Quantized Eigenvalues Sampling Transform, or QuEST for short, which discretizes, or quantizes, the relationship between F , H, and c defined in equations (2.1)-(2.3). For any positive integers n and p, the QuEST function, denoted by Q n,p , is defined as where ∀z ∈ C + m . . = m t n,p (z) is the unique solution in the set to the equation however, that the 'smoothed' version q i n,p (t) given in equation (2.17) leads to improved accuracy, higher stability, and faster computations of our numerical algorithm, to be detailed below, in practice.
Since F n is an empirical distribution (function), its quantiles are not uniquely defined. For example, the statistical software R offers nine different versions of sample quantiles in its function quantile; version 5 corresponds to our convention of considering λ n,i as the (i−0.5)/p quantile of F n .
Consequently, a set of (i − 0.5)/p quantiles (i = 1, . . . , p) is given by Q n,p (t) for F t n,p and is given by λ n for F n . The relationship between Q n,p (t) and λ n is further elucidated by the following theorem.
Theorem 2.1 states that the sample eigenvalues converge individually to their nonrandom QuEST function counterparts. This individual notion of convergence is defined as the Euclidian distance between the vectors λ n and Q n,p (τ n ), normalized by the matrix dimension p. It is the appropriate normalization because, as p goes to infinity, the left-hand side of equation (2.18) approximates the L 2 distance between the functions F −1 n and F τ n n,p −1 . This metric can be thought of as a 'cross-sectional' mean squared error, in the same way that F n is a cross-sectional distribution function.
Theorem 2.1 improves over the well-known results from the random matrix theory literature reviewed in Section 2.1 in two significant ways.
1) It is based on the p population eigenvalues τ n , not the limiting spectral distribution H.
Dealing with τ n (or, equivalently, H n ) is straightforward because it is integral to the actual data-generating process; whereas dealing with H is more delicate because we do not know how H n converges to H. Also there are potentially different H's that H n could converge to, depending on what we assume will happen as the dimension increases.
2) Theorem 2.1 characterizes the individual behavior of the sample eigenvalues, whereas equation (2.1) only characterizes their average behavior, namely, what proportion falls in any given interval. Individual results are more precise than average results. Thus, Theorem 2.1 shows that the sample eigenvalues are better behaved under large-dimensional asymptotics than previously thought.
Both of these improvements are made possible thanks to the introduction of the QuEST function. In spite of the apparent complexity of the mathematical definition of the QuEST function, it can be computed quickly and efficiently along with its analytical Jacobian as evidenced by Figure 1, and it behaves well numerically.

Consistent Estimation of Population Eigenvalues
Once the truth of Theorem 2.1 has been established, it becomes tempting to construct an estimator of population covariance matrix eigenvalues simply by minimizing the expression on the left-hand side of equation (2.18) over all possible sets of population eigenvectors. This is exactly what we do in the following theorem. where λ n . . = (λ n,1 , . . . , λ n,p ) ′ are the sample covariance matrix eigenvalues, and Q n,p (t) . . = q 1 n,p (t), . . . , q p n,p (t) ′ is the nonrandom QuEST function defined in equations (2.11)-(2.14); both τ n and λ n are assumed sorted in increasing order. Let τ n,i denote the ith entry of τ n (i = 1, . . . , p), and let τ n . . = (τ n,1 , . . . , τ n,p ) ′ denote the population covariance matrix eigenvalues sorted in increasing order. Then Theorem 2.2 shows that the estimated eigenvalues converge individually to the population eigenvalues, in the same sense as above, using the dimension-normalized Euclidian distance.
Remark 2.2. Mathematically speaking, equation (2.19) performs two tasks: it projects λ n onto the image of the QuEST function, and then inverts the QuEST function. Since the image of the QuEST function is a strict subset of [0, ∞) p , λ n will generally be outside of it. It is the first of these two tasks that gets around any potential ill-posedness by regularizing the set of observed sample eigenvalues.
Practically speaking, both tasks are performed simultaneously by a nonlinear optimizer. We use a standard off-the-shelf commercial software called SNOPT TM (Version 7.4), see Gill et al. (2002), but other choices may work well too.

Comparison with Other Approaches
El Karoui (2008) also attempts to discretize equation (2.5) and invert it, but he opts for a completely opposite method of discretization which does not exploit the natural discreteness of the population spectral distribution for finite p. If the population spectral distribution H n is approximated by a convex linear combination of step functions then in the optimization problem (2.19), we keep the weights w i (i = 1, . . . , p) fixed at 1/p while varying the location parameters t i (i = 1, . . . , p). In contrast, El Karoui (2008) does exactly the reverse: he keeps the location parameters t i fixed on a grid while varying the weights w i . Thus, El Karoui (2008) projects the population spectral distribution onto a "dictionary". Furthermore, instead of matching population eigenvalues to sample eigenvalues on R, he matches a function of m H to a function of m Fn on C + , which makes his algorithm relatively complicated; see Ledoit and Wolf (2012, pages 1043-1044. Despite our best efforts, we were unable to replicate his convergence results in Monte Carlo simulations: in our implementation, his estimator performs poorly overall and does not even appear to be consistent; see Section 5.1.1. Unless someone circulates an implementation of the algorithm described in El Karoui (2008) that works, we have to write off this approach as impractical. Another related paper is the one by Ledoit and Wolf (2012). They use the same discretization strategy as El Karoui (2008) (fix location parameters and vary weights) but, as we do here, match population eigenvalues to sample eigenvalues on the real line. Unlike we do here, they measure closeness by a sup-distance rather than by the Euclidean distance. Ledoit and Wolf (2012) only consider the case p < n. Unfortunately, their nonlinear optimizer no longer converges reliably in the case p > n, as we found out in subsequent experiments; this necessitated the development of the alternative discretization strategy described above, as well as the change from sup-distance to Euclidean distance to measure closeness.
Furthermore, Ledoit and Wolf (2012) are not directly interested in estimating the population eigenvalues; it is just an intermediary step towards their ultimate objective, which is the estimation of the covariance matrix itself. Therefore they do not report any Monte Carlo simulations on the finite-sample behavior of their estimator of the population eigenvalues.
In any case, the aim of the present paper is to develop an estimator of the population eigenvalues that works also when p > n, so the approach of Ledoit and Wolf (2012) is ruled out. The different discretization strategy that we employ here, together with the alternative distance measure, enables us to construct an estimator of τ n that works across both cases p < n and p > n. It is important to point out that the new estimator of population eigenvalues is not only more general, in the sense that it also works for the case p > n, but it also works better for the case p < n; see Section 5.2.
As for the papers of Mestre (2008) and Yao et al. (2012), their methods are based on contour integration of analytic functions in the complex plane. They can only extract a finite numberM of functionals of H n , such as the locations of high-multiplicity eigenvalue clusters, or the trace of powers of Σ n . The main difference with our method is that we extract many more items of information: namely, p population eigenvalues. This distinction is crucial because the ratioM /p vanishes asymptotically. It explains why we are able to recover the whole population spectrum in the general case, whereas they are not.

Covariance Matrix Estimation
The estimation of the covariance matrix Σ n is already considered by Ledoit and Wolf (2012), but only for the case p < n. In particular, they propose a nonlinear shrinkage approach, which we will now extend to the case p > n. To save space, the reader is referred to their paper for a more detailed discussion of the nonlinear shrinkage methodology and a comparison to other estimation strategies of large-dimensional covariance matrices, such as the linear shrinkage estimator of Ledoit and Wolf (2004).

Oracle Shrinkage
The starting point is to restrict attention to rotation-equivariant estimators of Σ n . To be more specific, let W be an arbitrary p-dimensional rotation matrix. Let Σ n . . = Σ n (Y n ) be an estimator of Σ n . Then the estimator is said to be rotation-equivariant if it satisfies Σ n (Y n W ) = W ′ Σ n (Y n )W . In other words, the estimate based on the rotated data equals the rotation of the estimate based on the original data. In the absence of any a priori knowledge about the structure of Σ n , such as sparseness or a factor model, it is natural to consider only estimators of Σ n that are rotation-equivariant.
The class of rotation-equivariant estimators of the covariance that are a function of the sample covariance matrix is constituted of all the estimators that have the same eigenvectors as the sample covariance matrix; for example, see Perlman (2007, Section 5.4). Every such rotation-equivariant estimator is thus of the form and where U n is the matrix whose ith column is the sample eigenvector u i . . = u n,i . This is the class of rotation-equivariant estimators already studied by Stein (1975Stein ( , 1986. We can rewrite the expression for such a rotation-equivariant estimator as This alternative expression shows that any such rotation-equivariant estimator is a linear combination of p rank-1 matrices u i u ′ i (i = 1, . . . , p). But since the {u i } form an orthonormal basis in R p , the resulting estimator is still of full rank p, provided that all the weights d i (i = 1, . . . , p) are strictly positive.
Remark 3.1 (Rotation-equivariant Estimators versus Structured Estimators). By construction, the class (3.1) of rotation-invariant estimators have the same eigenvectors as the sample covariance matrix. In particular, consistent estimation of the covariance matrix is not possible under large-dimensional asymptotics. Another approach would be to impose additional structure on the estimation problem, such as sparseness (Bickel and Levina, 2008), a graph model (Rajaratnam et al., 2008), or an (approximate) factor model (Fan et al., 2008). 1 The advantage of doing so is that, under suitable regularity conditions, consistent estimation of the covariance matrix is possible. The disadvantage is that if the assumed structure is misspecified, the estimator of the covariance matrix can be arbitrarily bad; and whether the structure is correctly specified can never be verified from the data alone.
Rotation-equivariant estimators are widely and successfully used in practice in situations where knowledge on additional structure is not available (or doubtful). This is evidenced by the many citations to Ledoit and Wolf (2004) who propose a linear shrinkage estimator that also belongs to the class (3.1); for example, see the beginning of Section 5.2. Therefore, developing a new, nonlinear shrinkage estimator that outperforms this previous proposal will be of substantial interest to applied researchers in our opinion.
The first objective is to find the matrix in the class (3.1) of rotation-equivariant estimators that is closest to Σ n . To measure distance, we choose the Frobenius norm defined as ||A|| F . . = Tr(AA ′ )/r for any matrix A of dimension r × m .
(3.3) (Dividing by the dimension of the square matrix AA ′ inside the root is not standard, but we do this for asymptotic purposes so that the Frobenius norm remains constant equal to one for the identity matrix regardless of the dimension; see Ledoit and Wolf (2004).) As a result, we end up with the following minimization problem: Elementary matrix algebra shows that its solution is Let y ∈ R p be a random vector with covariance matrix Σ n , drawn independently from the sample covariance matrix S n . We can think of y as an out-of-sample observation. Then d * i is recognized as the variance of the linear combination u ′ i y, conditional on S n . In view of the expression (3.2), it makes intuitive sense that the matrices u i u ′ i whose associated linear combination u ′ i y have higher out-of-sample variance should receive higher weight in computing the estimator of Σ n .
The finite-sample optimal estimator is thus given by Clearly S * n is not a feasible estimator because it depends on knowing the population covariance matrix. By generalizing the Marčenko-Pastur equation (2.5), Ledoit and Péché (2011) show that d * i can be approximated by the asymptotic quantities from which they deduce their oracle estimator The key difference between D * n and D or n is that the former depends on the unobservable population covariance matrix, whereas the latter depends on the limiting distribution of sample eigenvalues, F , which makes it amenable to consistent estimation. It turns out that this estimation problem is solved if a consistent estimator of the population eigenvalues τ n is available.

The Case p < n
We start with the case p < n, which was already considered by Ledoit and Wolf (2012). Silverstein and Choi (1995) show how the support of F , denoted by Supp(F ), is determined; also see Section 2.3 of Ledoit and Wolf (2012). Supp(F ) is seen to be the union of a finite number of disjoint compact intervals, bounded away from zero. To simplify the discussion, we will assume from here on that Supp(F ) is a single compact interval, bounded away from zero, with F ′ > 0 in the interior of this interval. But if Supp(F ) is the union of a finite number of such intervals, the arguments presented in this section as well as in the remainder of the paper apply separately to each interval. In particular, our consistency results presented below can be easily extended to this more general case.
Recall that, for any t . . = (t 1 , . . . , t p ) ′ ∈ [0, +∞) p , equations (2.13)-(2.14) define m t n,p as the Stieltjes transform of F t n,p , the limiting distribution of sample eigenvalues corresponding to the population spectral distribution +∞) . Its domain is the strict upper half of the complex plane, but it can be extended to the real line since Silverstein and Choi (1995) prove that ∀λ ∈ R − {0} lim z∈C + →λ m t n,p (z) = . .m t n,p (λ) exists. Ledoit and Wolf (2012) show how a consistent estimator ofm F can be derived from a consistent estimator of τ n , such as τ n defined in Theorem 2.2. Their Proposition 4.3 establishes thatm τ n n,p (λ) →m F (λ) uniformly in λ ∈ Supp(F ), except for two arbitrarily small regions at the lower and upper end of Supp(F ). Replacingm F withm τ n n,p and c with p/n in Ledoit and Péché (2011)'s oracle quantities d or i of (3.6) yields for i = 1, . . . , p .
(3.8) (Note here that in the case p < n, all sample eigenvalues λ i are positive almost surely, for n large enough, by the results of Bai and Silverstein (1998).) In turn, the bona fide nonlinear shrinkage estimator of Σ n is obtained as: (3.9)

The Case p > n
We move on to the case p > n, which was not considered by Ledoit and Wolf (2012). In this case, F is a mixture distribution with a discrete part and a continuous part. The discrete part is a point mass at zero with mass (c − 1)/c. The continuous part has total mass 1/c and its support is the union of a finite number of disjoint intervals, bounded away from zero; again, see Silverstein and Choi (1995). It can be seen from equations (2.6)-(2.9) that F corresponds to the continuous part of F , scaled to be a proper distribution (function): lim t→∞ F (t) = 1. Consequently, Supp(F ) = {0} ∪ Supp(F ). To simplify the discussion, we will assume from here on that Supp(F ) is a single compact interval, bounded away from zero, with F ′ > 0 in the interior of this interval. But if Supp(F ) is the union of a finite number of such intervals, the arguments presented in this section as well as in the remainder of the paper apply separately to each interval. In particular, our consistency results presented below can be easily extended to this more general case.
Finally, the bona fide nonlinear shrinkage estimator of Σ n is obtained as (3.9) but now with

Strong Consistency
The following theorem establishes that our nonlinear shrinkage estimator, based on the estimator τ n of Theorem 2.2, is strongly consistent for the oracle estimator across both cases p < n and p > n.
Theorem 3.1. Let τ n be an estimator of the eigenvalues of the population covariance matrix satisfying p −1 p i=1 [ τ n,i − τ n,i ] 2 a.s. −→ 0. Define the nonlinear shrinkage estimator S n as in (3.9), where the d i are as in (3.8) in the case p < n and as in (3.11) in the case p > n.
Remark 3.2. We have to rule out the case c = 1 (or p = n) for mathematical reasons. First, we need Supp(F ) to be bounded away from zero to establish various consistency results. But when c = 1, then Supp(F ) can start at zero, that is, there exists u > 0 such that F ′ (λ) > 0 for all λ ∈ (0, u). This was already established by Marčenko and Pastur (1967) for the special case when H is a point mass at one. In particular, the resulting (standard) Marčenko-Pastur distribution F has density function for example, see Bai and Silverstein (2010, Section 3.3.1). Second, we also need the assumption c = 1 'directly' in the proof of Theorem 3.1 to demonstrate that the summand D 1 in (A.17) converges to zero.
Although the case c = 1 is not covered by the mathematical treatment, we can still address it in Monte Carlo simulations; see Section 5.2.

Principal Component Analysis
Principal component analysis (PCA) is one of the oldest and best-known techniques of multivariate analysis, dating back to Pearson (1901) and Hotelling (1933); for a comprehensive treatment, see Jolliffe (2002).

The Central Idea and the Common Practice
The central idea of PCA is to reduce the dimensionality of a data set consisting of a large number of interrelated variables, while retaining as much as possible of the variation present in the data set. This is achieved by transforming the original variables to a new set of uncorrelated variables, the principal components, which are ordered so that the 'largest' few retain most of the variation present in all of the original variables.
More specifically, let y ∈ R p be a random vector with covariance matrix Σ; in this section, it will be convenient to drop the subscript n from the covariance matrix and related quantities. Let ((τ 1 , . . . , τ p ); (v 1 , . . . , v p )) denote a system of eigenvalues and eigenvectors of Σ. To be consistent with our former notation, we assume that the eigenvalues τ i are sorted in increasing order. Then the principal components of y are given by v ′ 1 y, . . . , v ′ p y. Since the eigenvalues τ i are sorted in increasing order, the principal component with the largest variance is v ′ p y and the principal component with the smallest variance is v ′ 1 y. The eigenvector v i is called the vector of coefficients or loadings for the ith principal component (i = 1, . . . , p).
Two brief remarks are in order. First, some authors use the term principal components for the eigenvectors v i ; but we agree with Jolliffe (2002, Section 1.1) that this usage is confusing and that it is preferable to reserve the term for the derived variables v ′ i y. Second, in the PCA literature, in contrast to the bulk of the multivariate statistics literature, the eigenvalues τ i are generally sorted in decreasing order so that v ′ 1 y is the 'largest' principal component (that is, the principal component with the largest variance). This is understandable when the goal is expressed as capturing most of the total variation in the first few principal components. But to avoid confusion with other sections of this paper, we keep the convention of eigenvalues being sorted in increasing order, and then express the goal as capturing most of the total variation in the largest few principal components.
The k largest principal components in our notation are thus given by v ′ p y, . . . , v ′ p−k+1 y (k = 1, . . . , p). Their (cumulative) fraction captured of the total variation contained in y, denoted by f k (Σ), is given by (4.1) The most common rule in deciding how many principal components to retain is to decide on a given fraction of the total variation that one wants to capture, denoted by f target , and to then retain the largest k principal components, where k is the smallest integer satisfying f k (Σ) ≥ f target . Commonly chosen values of f target are 70%, 80%, 90%, depending on the context. For obvious reasons, this rule is known as the cumulative-percentage-of-total-variation rule.
There exist a host of other rules, either analytical or graphical, such as Kaiser's rule or the scree plot; see Jolliffe (2002, Section 6.1). The vast majority of these rules are also solely based on the eigenvalues (τ 1 , . . . , τ p ).
The problem is that generally the covariance matrix Σ is unknown. Thus, neither the (population) principal components v ′ i y nor their cumulative percentages of total variation f k (Σ) can be used in practice.
The common solution is to replace Σ with the sample covariance matrix S, computed from a random sample y 1 , . . . , y n , independent of y. Let ((λ 1 , . . . , λ p ); (u 1 , . . . , u p )) denote a system of eigenvalues and eigenvectors of S; it is assumed again that the eigenvalues λ i are sorted in increasing order. Then the (sample) principal components of y are given by u ′ 1 y, . . . , u ′ p y. The various rules in deciding how many (sample) principal components to retain are now based on the sample eigenvalues λ i . For example, the cumulative-percentage-of-total-variation rule retains the largest k principal components, where k is the smallest integer satisfying (4.2) The pitfall in doing so, unless p ≪ n, is that λ i is not good estimator of the variance of the ith principal component. Indeed, the variance of the ith principal component, u ′ i y, is given by u ′ i Σu i rather than by λ i = u ′ i Su i . By design, for large values of i, the estimator λ i is upward biased for the true variance u ′ i Σu i , whereas for small values of i, it is downward biased. In other words, the variances of the large principal components are overestimated whereas the variances of the small principal components are underestimated. The unfortunate consequence is that most rules in deciding how many principal components to retain, such as the cumulative-percentage-of-total-variation rule, generally retain fewer principal components than really needed.

Previous Approaches under Large-Dimensional Asymptotics
All the previous approaches under large-dimensional asymptotics that we are aware of impose some additional structure on the estimation problem.
Most works assume a sparseness conditions on the eigenvectors v i or on the covariance matrix Σ; see Amini (2011) for a comprehensive review.
Mestre (2008), on the other hand and as discussed before, assumes that Σ has onlyM ≪ p distinct eigenvalues and further that the multiplicity of each of theM distinct eigenvalues is known (which implies that the numberM is known as well). Furthermore, he needs spectral separation. In this restrictive setting, he is able to construct a consistent estimator of every distinct eigenvalue and its associated eigenspace (that is, the space spanned by all eigenvectors corresponding to a specific distinct eigenvalue).

Alternative Approach Based on Nonlinear Shrinkage
Unlike previous approaches under large-dimensional asymptotics, we do not wish to impose additional structure on the estimation problem. In such a setting, improved estimators of the eigenvectors v i are not available and one must indeed use the sample eigenvectors u i as loadings. Therefore, as in common practice, the principal components used are the u ′ i y. Ideally, the rules in deciding how many principal components to retain should be based on the variances of the principal components given by u ′ i Σu i . It is important to note that even if the population eigenvalues τ i were known, the rules should not be based on them. This is because the population eigenvalues τ i = v ′ i Σv i describe the variances of the v ′ i y, which are not available and thus not used. It seems that this important point has not been realized so far. Indeed, various authors have used PCA as a motivational example in the estimation of the population eigenvalues τ i ; for example, see El Karoui (2008), Mestre (2008), and Yao et al. (2012). But unless the population eigenvectors v i are available as well, using the τ i is misleading.
Although, in the absence of additional structure, it is not possible to construct improved principal components, it is possible to accurately estimate the variances of the commonly-used principal components. This is because the variance of the ith principal component is nothing else than the finite-sample-optimal nonlinear shrinkage constant d * i ; see equation (3.4). Its oracle counterpart d or i is given in equation (3.6) and the bona fide counterpart d i is given in equation (3.8) in the case p < n and in equation (3.11) in the case p > n.
Our solution then is to base the various rules in deciding how many principal components to retain on the d i in place of the unavailable d * i = u ′ i Σu i . For example, the cumulativepercentage-of-total-variation rule retains the k largest principal components, where k is the smallest integer satisfying f k ( S) ≥ f target , with (4.3) Remark 4.1. We have taken the total variation to be p m=1 d * m , and the variation attributable to the k largest principal components to be k j=1 d * p−j+1 . In general, the sample principal components u ′ i y are not uncorrelated (unlike the population principal components v ′ i y). This means that u ′ i Σu j can be non-zero for i = j. Nonetheless, even in this case, the variation attributable to the k largest principal components is still equal to k j=1 d * p−j+1 , as explained in Appendix B.
While most applications of PCA seek the principal components with the largest variances, there are also some applications of PCA that seek the principal components with the smallest variances; see Jolliffe (2002, Section 3.4). In the case p > n, a certain number of the λ i will be equal to zero, falsely giving the impression that a certain number of the smallest principal components have variance zero. Such applications also highlight the use of replacing the λ i with our nonlinear shrinkage constants d i , which are always greater than zero.

Monte Carlo Simulations
In this section, we study the finite-sample performance of various estimators in different settings.

Estimation of Population Eigenvalues
We first focus on estimating the eigenvalues of the population covariance matrix, τ n . Of major interest to us is the case where all the eigenvalues are or can be distinct; but we also consider the case where they are known or assumed to be grouped into a small number of high-multiplicity clusters.

All Distinct Eigenvalues
We consider the following estimators of τ n .
• Sample: The sample eigenvalues λ n,i .

• Lawley:
The bias-corrected sample eigenvalues using the formula of Lawley (1956, Section 4). This transformation may not be monotonic in finite samples. Therefore, we post-process it with an isotonic regression.
• El Karoui: The estimator of El Karoui (2008). It provides an estimator of H n , not τ n , so we derive estimates of the population eigenvalues using 'smoothed' quantiles in the spirit of equations (2.17)-(2.16). 2 • LW: Our estimator τ n of Theorem 2.2.
It should be pointed out that the estimator of Lawley (1956) is designed to reduce the finitesample bias of the sample eigenvalues λ n,i ; it is not necessarily designed for consistent estimation of τ n under large-dimensional asymptotics. Let τ n,i denote a generic estimator of τ n,i . The evaluation criterion is the dimensionnormalized Euclidian distance between estimated eigenvalues τ n and population eigenvalues averaged over 1,000 Monte Carlo simulations in each scenario.

Convergence
In the first design, the ith population eigenvalue is equal to τ n,i . . = H −1 ((i − 0.5)/p) (i = 1, . . . , p), where H is given by the distribution of 1 + 10W , and W ∼ Beta(1, 10); this distribution is right-skewed and resembles in shape an exponential distribution. The distribution of the random variates comprising the n × p data matrix X n is real Gaussian. We fix the concentration at p/n = 0.5 and vary the dimension from p = 30 to p = 1, 000. The results are displayed in Figure 2. It can be seen that the empirical mean squared error for LW converges to zero, which is in agreement with the proven consistency of Theorem 2.2. For all the other estimators, the average distance from τ n appears bounded away from zero. This simulation also shows that dividing by p is indeed the appropriate normalization for the Euclidian norm in equation (2.18), as it drives a wedge between estimators such as the sample eigenvalues that are not consistent and τ n , which is consistent.
The second design is similar to the first design, except that we fix the concentration at p/n = 2 and now vary the sample size from n = 30 to n = 1, 000 instead of the dimension. In this design, the sample covariance matrix is always singular. The results are displayed in Figure 3 and are qualitatively similar. Again, LW is the only estimator that appears to be consistent. Notice the vertical scale: the difference between El Karoui and LW is of the same order of magnitude as in Figure 2, but Sample and Lawley are vastly more erroneous now.

Condition Number
In the third design, the focus is on the condition number. The ith population eigenvalue is still τ n,i . . = H −1 ((i − 0.5)/p) (i = 1, . . . , p), but H is now given by the distribution of a + 10W , where W ∼ Beta(1, 10), and a ∈ [0, 7]. As a result, the smallest eigenvalue approaches a, and the previously-used distribution for H is included as a special case when a = 1. The condition number decreases in a from approximately 10, 000 to 2.4.
We use n = 1, 600 and p = 800, so that p/n = 0.5. The random variates are still real Gaussian. The results are displayed in Figure 4. It can be seen that Sample and Lawley perform quite well for values of a near zero (that is, for very large condition numbers) but their performance gets worse as a increases (that is, as the condition number decreases). On the other hand, the performance of El Karoui is more stable across all values of a, though relatively bad. The performance of LW is uniformly the best and also stable across a.
We use n = 1, 600 and p = 800, so that p/n = 0.5. The random variates are real Gaussian. The results are presented in Table 1. It can be seen that LW is uniformly best and Sample is uniformly worst. There is no clear-cut ranking for the remaining two estimators. On average, Lawley is second best, followed by El Karoui.

Heavy Tails
So far, the variates making up the data matrix X n always had a Gaussian distribution. It is also of interest to consider a heavy-tailed distribution instead. We return to the first design with n = 1600 and p = 800, so that p/n = 0.5. In addition to the Gaussian distribution, which can be viewed as a t-distribution with infinite degrees of freedom, we also consider a the t-distribution with three degrees of freedom (scaled to have unit variance). The results are presented in Table 2. It can be seen that all estimators perform worse when the degrees of freedom are changed from infinity to three, but LW is still by far the best.  Table 2: Mean squared error between estimated and population eigenvalues.

Clustered Eigenvalues
We are mainly interested in the case where the population eigenvalues are or can be distinct, but it is also worthwhile seeing how (an adapted version of) our estimator of τ n compares to the one of Mestre (2008) in the setting where the population eigenvalues are known or assumed to be grouped into a small number of high-multiplicity clusters. Let γ 1 < γ 2 < · · · < γM denote the set of pairwise different eigenvalues of the population covariance matrix Σ, whereM is the number of distinct population eigenvalues (1 ≤M < p). Each of the eigenvalues γ j has known multiplicity K j (j = 1, . . . ,M ), so that p = M j=1 K j . (Knowing the multiplicities of the eigenvalues γ j comes from knowing their masses m j in the limiting spectral distribution H, as assumed in Mestre (2008): Then the optimization problem in Theorem 2.1 becomes: subject to: t = (γ 1 , . . . , γ 1 We consider the following estimators of τ n . • Traditional: γ j is estimated by the average of all corresponding sample eigenvalues λ n,i ; under the condition of spectral separation assumed in Mestre (2008), it is known which γ j corresponds to which λ n,i .
The mean squared error criterion (5.1) specializes in this setting tō We report the average MSE over 1,000 Monte Carlo simulations in each scenario.

Convergence
The first design is based on Tables I and II of Mestre (2008). The distinct population eigenvalues are (γ 1 , γ 2 , γ 3 , γ 4 ) = (1, 7, 15, 25) with respective multiplicities (K 1 , K 2 , K 3 , K 4 ) = (p/2, p/4, p/8, p/8). The distribution of the random variates comprising the n × p data matrix X n is circular symmetric complex Gaussian, as in Mestre (2008). We fix the concentration at p/n = 0.32 and vary the dimension from p = 8 to p = 1, 000; the lower end p = 8 corresponds to Table I in Mestre (2008), while the upper end p = 1, 000 corresponds to Table II in Mestre (2008). The results are displayed in Figure 5. It can be seen that the average MSE of both Mestre and LW converges to zero, and that the performance of the two estimators is nearly indistinguishable. On the other hand, Traditional is seen to be inconsistent, as its MSE remains bounded away from zero.

Performance When One Eigenvalue Is Isolated
The second design is based on Table III of Mestre (2008). The distinct population eigenvalues are (γ 1 , γ 2 , γ 3 , γ 4 ) = (1, 7, 15, 25) with multiplicities (K 1 , K 2 , K 3 , K 4 ) = (160, 80, 79, 1). There is a single 'isolated' large eigenvalue. The distribution of the random variates comprising the n × p data matrix X n is circular symmetric complex Gaussian, as in Mestre (2008). We use n = 1, 000 and p = 320, so that p/n = 0.32. The averages and the standard deviations of the estimates γ j over 10,000 Monte Carlos simulations are presented in Table 3; note here that the numbers for Traditional and Mestre have been directly copied from Table III of Mestre (2008). The inconsistency of Traditional is again apparent. In terms of estimating (γ 1 , γ 2 , γ 3 ), the performance of Mestre and LW is nearly indistinguishable. In terms of estimating γ 4 , Mestre has a smaller bias (in absolute value) while LW has a smaller standard deviation; combining the two criteria yields a mean squared error of (25 − 24.9892) 2 + 1.0713 2 = 1.1478 for Mestre and a mean squared error of (25 − 24.9238) 2 + 0.8898 2 = 0.7976 for LW.

Traditional
Mestre LW   (2008), sample eigenvalues, and the proposed estimator. The first six columns are copied from Table III of Mestre (2008). Results are based on 10, 000 Monte Carlo simulations with circularly symmetric complex Gaussian random variates.

Covariance Matrix Estimation
As detailed in Section 3.1, the finite-sample optimal estimator in the class of rotation-equivariant estimators is given by S * n as defined in (3.5). As the benchmark, we use the linear shrinkage estimator of Ledoit and Wolf (2004) instead of the sample covariance matrix. We do this because the linear shrinkage estimator has become the de facto standard among leading researchers because of its simplicity, accuracy, and good conditioning properties. It has been used in several fields of statistics, such as linear regression with a large number of regressors (Anatolyev, 2012), linear discriminant analysis (Pedro Duarte Silva, 2011), factor analysis (Lin and Bentler, 2012), unit root tests (Demetrescu and Hanck, 2012), and vector autoregressive models (Huang and Schneider, 2011), among others. Beyond pure statistics, the linear shrinkage estimator has been applied in finance for portfolio selection (Tsagaris et al., 2012) and tests of asset pricing models (Khan, 2008); in signal processing for cellular phone transmission (Nguyen et al., 2011) and radar detection (Wei et al., 2011); and in biology for neuroimaging (Varoquaux et al., 2010), genetics , cancer research (Pyeon et al., 2007), and psychiatry (Markon, 2010). It has also been used in such varied applications as physics (Pirkl et al., 2012), chemistry (Guo et al., 2012), climatology (Ribes et al., 2009), oil exploration (Saetrom et al., 2012), road safety research (Haufe et al., 2011), etc. In summary, the comparatively poor performance of the sample covariance matrix and the popularity of the linear shrinkage estimator justify taking the latter as the benchmark.
The improvement of the nonlinear shrinkage estimator S n over the linear shrinkage estimator of Ledoit and Wolf (2004), denoted by S n , will be measured by how closely this estimator approximates the finite-sample optimal estimator S * n relative to S n . More specifically, we report the Percentage Relative Improvement in Average Loss (PRIAL), which is defined as where Σ n is an arbitrary estimator of Σ n . By definition, the PRIAL of S n is 0% while the PRIAL of S * n is 100%. We consider the following estimators of Σ n .
• LW (2012) Estimator: The nonlinear shrinkage estimator of Ledoit and Wolf (2012); this version only works for the case p < n.
• New Nonlinear Shrinkage Estimator: The new nonlinear shrinkage estimator of Section 3.2; this version works across both cases p < n and p > n.

Convergence
In our design, 20% of the population eigenvalues are equal to 1, 40% are equal to 3, and 40% are equal to 10. This is a particularly interesting and difficult example introduced and analyzed in detail by Bai and Silverstein (1998); it has also been used in previous Monte Carlo simulations by Ledoit and Wolf (2012). The distribution of the random variates comprising the n × p data matrix X n is real Gaussian. We study convergence of the various estimators by keeping the concentration p/n fixed while increasing the sample size n. We consider the three cases p/n = 0.5, 1, 2; as discussed in Remark 3.2, the case p/n = 1 is not covered by the mathematical treatment. The results are displayed in Figure 6, which shows empirical PRIAL's across 1,000 Monte Carlo simulations (one panel for each case p/n = 0.5, 1, 2). It can be seen that the new nonlinear shrinkage estimator always outperforms linear shrinkage with its PRIAL converging to 100%, though slower than the oracle estimator. As expected, the relative improvement over the linear shrinkage estimator is inversely related to the concentration ratio; also see Figure 4 of Ledoit and Wolf (2012). In the case p < n, it can also be seen that the new nonlinear shrinkage estimator slightly outperforms the earlier nonlinear shrinkage estimator of Ledoit and Wolf (2012). Last but not least, although the case p = n is not covered by the mathematical treatment, it is also dealt with successfully in practice by the new nonlinear shrinkage estimator.

Principal Component Analysis
Recall that in Section 4 on principal component analysis we dropped the first subscript n always, and so the same will be done in this section.
In our design, the ith population eigenvalue is equal to τ i = H −1 ((i − 0.5)/p) (i = 1, . . . , p), where H is given by the distribution of 1 + 10W , and W ∼ Beta(1, 10) The distribution of the random variates comprising the n × p data matrix X n is Gaussian. We consider the two cases (n = 200, p = 100) and (n = 100, p = 200), so the concentration is p/n = 0.5 or p/n = 2.
Let y ∈ R p be a random vector with covariance matrix Σ, drawn independently from the sample covariance matrix S. The out-of-sample variance of the ith (sample) principal component, u ′ i y, is given by d * i . . = u ′ i Σu i ; see (3.4). By our convention, the d * i are sorted in increasing order.
We consider the following estimators of d * i .
• Sample: The estimator of d * i is the ith sample eigenvalue, λ i . • Population: The estimator of d * i is the ith population eigenvalue, τ i ; this estimator is not feasible but is included for educational purposes nevertheless.
• LW: The estimator of d * i is the nonlinear shrinkage quantity d i as given in equation (3.8) in the case p < n and in equation (3.11) in the case p > n.
Let d i be a generic estimator of d * i . First, we are plotting as a function of k, averaged over 1,000 Monte Carlo simulations. The quantity f k serves as an estimator of f k , the fraction of the total variation in y that is explained by the k largest principal components: The results are displayed in Figure 7 (one panel for each case p/n = 0.5, 2.) The upward bias of Sample is apparent, while LW is very close to the Truth. Moreover, Population is also upward biased (though not as much as Sample): the important message is that even if the population eigenvalues were known, they should not be used to judge the variances of the (sample) principal components. As expected, the differences between Sample and LW increase with the concentration p/n; the same is true for the differences between Population and LW.  Figure 7 shows how close the estimator f k is to the truth f k on average. But it does not necessarily answer how close the cumulative-percentage-of-total-variation rule based on f k is to the rule based on f k . For a given percentage (q × 100)%, with q ∈ (0, 1), let k(q) . . = min{k : f k ≥ q} and k(q) . . = min{k : f k ≥ q} .
In words, k(q) is the (smallest) number of the largest principal components that must be retained to explain (q × 100)% of the total variation and k(q) is an estimator of this quantity. We are then also interested in the Root Mean Squared Error (RMSE) of k(q), defined as for the values of q most commonly used in practice, namely q = 0.7, 0.8, 0.9. We compute empirical RMSEs across 1,000 Monte Carlo simulations. The results are presented in Table 4.
It can be seen that in each scenario, Sample has the largest RMSE and LW has the smallest RMSE; in particular LW is highly accurate not only relative to Sample but also in an absolute sense. (The quantity k(q) is a random variable, since it depends on the sample eigenvalues u i ; but the general magnitude of k(q) for the various scenarios can be judged from Figure 7)

Empirical Application
As an empirical application, we study principal component analysis (PCA) in the context of stock return data. Principal components of a return vector of a cross section of p stocks are used for risk analysis and portfolio selection by finance practitioners; for example, see Roll and Ross (1980) and Connor and Korajczyk (1993). We use the p = 30, 60, 240, 480 largest stocks, as measured by their market value at the beginning of 2011, that have a complete return history from January 2002 until December 2011. As is customary in many financial applications, such as portfolio selection, we use monthly data. Consequently, the sample size for the ten-year history is n = 120 and the concentration is p/n = 0.25, 0.5, 2, 4.
It is of crucial interest how much of the total variation in the p-dimensional return vector is explained by the k largest principal components. We compare the approach based on the sample covariance matrix, defined in equation (4.2) and denoted by Sample, to that of nonlinear shrinkage, defined in equation 4.3 and denoted by LW. The results are displayed in Figure 8. It can be seen that Sample is overly optimistic compared to LW and, as expected, the differences between the two methods increase with the concentration p/n. In addition to the visual analysis, the differences can also be presented via the cumulativepercentage-of-total-variation rule to decide how many of the largest principal components to retain; see Section 4. The results are presented in Table 5. It can be seen again that Sample is overly optimistic compared to LW and retains much fewer principal components. Again, as expected, the differences between the two methods increase with the concentration p/n.

Conclusion
The analysis of large-dimensional data sets is becoming more and more common. For many statistical problems, the classic textbook methods no longer work well in such settings. Two cases in point are covariance matrix estimation and principal component analysis, both cornerstones of multivariate analysis. The classic estimator of the covariance matrix is the sample covariance matrix. It is unbiased and the maximum-likelihood estimator under normality. But when the dimension is not small compared to the sample size, the sample covariance matrix contains too much estimation error and is ill-conditioned; when the dimension is larger than the sample size, it is not even invertible anymore.
The variances of the principal components (which are obtained from the sample eigenvectors) are no longer accurately estimated by the sample eigenvalues. In particular, the sample eigenvalues overestimate the variances of the 'large' principal components. As a result the common rules in determining how many principal components to retain generally select too few of them.
In the absence of strong structural assumption on the true covariance matrix, such as sparseness or a factor model, a common remedy for both statistical problems is nonlinear shrinkage of the sample eigenvalues. The optimal shrinkage formula delivers an estimator of the covariance matrix that is finite-sample optimal with respect to the Frobenius norm in the class of rotation-equivariant estimators. The same shrinkage formula also gives the variances of the (sample) principal components. It is noteworthy that the optimal shrinkage formula is different from the population eigenvalues: even if they were available, one should not use them for the ends of covariance matrix estimation and PCA. Unsurprisingly, the optimal nonlinear shrinkage formula is not available, since it depends on population quantities. But an asymptotic counterpart, denoted oracle shrinkage, can be estimated consistently. In this way, bona fide nonlinear shrinkage estimation of covariance matrices and improved PCA result.
The key to the consistent estimation of the oracle shrinkage is the consistent estimation of the population eigenvalues. This problem is challenging and interesting in its own right, solving a host of additional statistical problems. Our proposal to this end is the first one that does not make strong assumptions on the distribution of the population eigenvalues, has proven consistency, and also works well in practice.
Extensive Monte Carlo simulations have established that our methods have desirable finitesample properties and outperform methods that have been previously suggested in the literature.

A Mathematical Proofs
Lemma A.1. Let {G n } and G be c.d.f.'s on the real line and assume that there exists a compact interval that contains the support of G as well as the support of G n for all n large enough. For 0 < α < 1, let G −1 n (α) denote an α quantile of G n and let G −1 (α) denote an α quantile of G. Let {K n } be a sequence of integers with K n → ∞; further, let t n,k . . = (k − 0.5)/K n (k = 1, . . . , K n ).
Then G n ⇒ G if and only if 1 K n Kn k=1 G −1 n (t n,k ) − G −1 (t n,k ) 2 → 0 . (A.1) recalling that λ n,i is a (i − 0.5)/p quantile of F n ; see Remark 2.1. The additional fact that 1 p p i=1 q i n,p (τ n ) − F −1 ((i − 0.5)/p) In the case c < 1, combining Corollary A.1(i) with Lemma A.2 of Ledoit and Wolf (2012) shows that A → 0 almost surely. In the case c > 1, one can also show that A → 0 almost surely: Lemma A.2 of Ledoit and Wolf (2012) implies that F H, cn ⇒ F almost surely; then use equation (2.7) together with the fact that c n → c to deduce that also F H, cn ⇒ F almost surely; finally apply Corollary A.1(i). Combining Corollary A.1(i) with the fact that F n ⇒ F almost surely (Silverstein, 1995) shows that B → 0 almost surely in addition to A → 0 almost surely. Therefore, ||F H n , cn − F || p → 0 almost surely. Using Corollary A.1(i) again shows that F Hn, cn ⇒ F almost surely. A feasible estimator of H is given by The fact that the minimization over a finite but dense family of probability measures, instead of all probability measures on the nonnegative real line, does not affect the strong consistency of the estimator of F follows by arguments similar to those used in the proof of Corollary 5.1(i) of Ledoit and Wolf (2012). Therefore, it also holds that F Hn, cn ⇒ F almost surely.
Having established that F Hn, cn ⇒ F almost surely, it follows that also H n ⇒ H almost surely; see the proof of Theorem 5.1(ii) of Ledoit and Wolf (2012). Since H n is recognized as the empirical distribution (function) of the τ n,i (i = 1, . . . , p), τ n,i is a (i − 0.5)/p quantile of H n ; see Remark 2.1. Therefore, it follows from Corollary A.1(i) that [ τ n,i − τ n,i ] 2 a.s. −→ 0 , which is the statement to be proven.
Proof of Theorem 3.1.
The claim for the case p < n follows immediately from Proposition 4.3(ii) of Ledoit and Wolf (2012).