Convergence and prediction of principal component scores in high-dimensional settings

Introduction

Principal component analysis (PCA) jolliffe2002pca is one of the leading statistical tools for analyzing multivariate data. It is especially popular in genetics/genomics, medical imaging and chemometrics studies where high-dimensional data is common. PCA is typically used as a dimension reduction tool. A small number of top ranked principal component (PC) scores are computed by projecting data onto spaces spanned by the eigenvectors of sample covariance matrix, and are used to summarize data characteristics that contribute most to data variation. These PC scores can be subsequently used for data exploration and/or model predictions. For example, in genome-wide association studies (GWAS), PC scores are used to estimate ancestries of study subjects and as covariates to adjust for population stratification price2006pca; patterson2006psa. In gene expression microarray studies, PC scores are used as synthetic “eigen-genes” or “meta-genes” intended to represent and discover gene expression patterns that might not be discernible from single-gene analysis wall2003svd.

Although PCA is widely applied in a number of settings, much of our theoretical understanding rests on a relatively small body of literature. Girshick girshick1936principal introduced the idea that the eigenvectors of sample covariance matrix are maximum likelihood estimators. Here, a key concept in a population view of PCA is that the data arise as $p$-variate values from a distinct set of $n$ independent samples. Later, the asymptotic distribution of eigenvalues and eigenvectors of the sample covariance matrix (i.e., the sample eigenvalues and eigenvectors) were derived for the situation where $n$ goes to infinity and $p$ is fixed girshick1939sampling; Anderson1. With the development of modern high-throughput technologies, it is not uncommon to have data where $p$ is comparable in size to $n$, or substantially larger. Under the assumption that $p$ and $n$ grow at the same rate, that is $p/n\rightarrow\gamma>0$, there has been considerable effort to establish convergence results for sample eigenvalues and eigenvectors (see review bai1999msa). The convergence of the sample eigenvalues and eigenvectors under the “spiked population” model proposed by Johnstone johnstone2001 has also been established Baik2006; paul2007; nadler2008fsa. For this model, it is well known that the sample eigenvectors are not consistent estimators of the eigenvectors of population covariance (i.e., the population eigenvectors) johnstone2007spc; paul2007; nadler2008fsa. Furthermore, Paul paul2007 has derived the degree of discrepancy in terms of the angle between the sample and population eigenvectors, under Gaussian assumptions for $0<\gamma<1$. More recently, Nadler nadler2008fsa has extended the same result to the more general $\gamma>0$ using a matrix perturbation approach.

These results have considerable potential practical utility in understanding the behavior of PC analysis and prediction in modern datasets, for which $p$ may be large. The practical goals of this paper focus primarily on the prediction of PC scores for samples which were not included in the original PC analysis. For example, gene expression data of new breast cancer patients may be collected, and we might want to estimate their PC scores in order to classify their cancer sub-type. The recalculation of PCs using both new and old data might not be practical. For example, if the application of PCs from gene expression is used as a diagnostic tool in clinical applications. For GWAS analysis, it is known that PC analysis which includes related individuals tends to generate spurious PC scores which do not reflect the true underlying population substructures. To overcome this problem, it is common practice to include only one individual per family/sibship in the initial PC analysis. Another example arises in cross-validation for PC regression, in which PC scores for the test set might be derived using PCA performed on the training set jackson2005user. For all of these applications, the predicted PC scores for a new sample are usually estimated in the “naive” fashion, in which the data vector of the new sample is multiplied by the sample eigenvectors from the original PC analysis. Indeed, there appears to be relatively little recognition in the genetics or data mining literature that this approach may lead to misleading conclusions.

For low-dimensional data, where $p$ is fixed as $n$ increases or otherwise much smaller than $n$, the predicted PC scores are nearly unbiased and well-behaved. However, for high-dimensional data, particularly with $p>n$, they tend to be biased and shrunken toward 0. The following simple example of a stratified population with three strata illustrates the shrinkage phenomenon for predicted PC scores. We generated a training data set with $n=100$ and $p=5000$. Among the 100 samples, 50 are from stratum 1, 30 are from stratum 2 and the rest from stratum 3. For each stratum, we first created a $p$-dimensional mean vector $\bm{\mu}_{k}$ $(k=1,2,3)$. Each element of each mean vector was created by drawing randomly with replacement from $\{-0.3,0,0.3\}$, and thereafter considered a fixed property of the stratum. Then for each sample from the $k$th stratum, its $p$ covariates were simulated from the multivariate normal distribution $\operatorname{MVN}(\bm{\mu}_{k},4\mathbf{I})$, where $\mathbf{I}$ is the $p\times p$ identity matrix. A test dataset with the same sample size and $\bm{\mu}_{k}$ vectors was also simulated. Figure 1 shows that the predicted PC scores for the test data are much closer to 0 compared to the scores from the training data. This shrinkage phenomenon may create a serious problem if the predicted PC scores are used to classify new test samples, perhaps by similarity to previous apparent clusters in the original data. In addition, the predicted PC scores may produce incorrect results if used for downstream analyses (e.g., as covariates in association analyses).

In this paper, we investigate the degree of shrinkage bias associated with the predicted PC scores, and then propose new bias-adjusted PC score estimates. As the shrinkage phenomenon is largely related to the limiting behavior of the sample eigenvectors, our first step is to describe the discrepancy between the sample and population eigenvectors. To achieve this purpose, we follow the assumption that $p$ and $n$ both are large and grow at the same rate. By applying and extending results from random matrix theory, we establish the convergence of the sample eigenvalues and eigenvectors under the spiked population model. We generalize Theorem 4 of Paul paul2007, which describes the asymptotic angle between sample and population eigenvectors, to non-Gaussian random variables for any $\gamma>0$. We further derive the asymptotic angle between PC scores from sample eigenvectors and population eigenvectors, and the asymptotic shrinkage factor of the PC score predictions. Finally, we construct estimators of the angles and the shrinkage factor. The theoretical results are presented in Section 2.

In Section 3, we report simulations to assess the finite sample accuracy of the proposed asymptotic angle and shrinkage factor estimators. We also show the potential improvements in prediction accuracy for PC regression by using the bias-adjusted PC scores. In Section 4, we apply our PC analysis to a real genome-wide association study, which demonstrates that the shrinkage phenomenon occurs in real studies and that adjustment is needed.

Method

Define the $p\times n$ data matrix, $\mathbf{X}$ as $[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]$, where $\mathbf{x}_{j}$ is the $p$-dimensional vector corresponding to the $j$th sample. For the remainder of the paper, we assume the following.

$\mathbf{X}=\mathbf{E}\bm{\Lambda}^{1/2}\mathbf{Z}$, where $\mathbf{Z}=\{z_{ij}\}$ is a $p\times n$ matrix whose elements $z_{ij}$’s are i.i.d. random variables with $E(z_{ij})=0,E(z_{ij}^{2})=1$ and $E(z_{ij}^{4})<\infty$.

Although the $z_{ij}$’s are i.i.d., Assumption 1 allows for very flexible covariance structures for $\mathbf{X}$, and thus the results of this paper are quite general. The population covariance matrix of $\mathbf{X}$ is $\bm{\Sigma}=\mathbf{E}\bm{\Lambda}\mathbf{E}^{T}$. The sample covariance matrix $\mathbf{S}$ equals

The $\lambda_{k}$’s are the underlying population eigenvalues. The spiked population model defined in johnstone2001 assumes that all the population eigenvalues are 1, except the first $m$ eigenvalues. That is, $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{m}>\lambda_{m+1}=\cdots=\lambda_{p}=1$. The spectral decomposition of the sample covariance matrix is

2 Sample eigenvalues and eigenvectors

Under the classical setting of fixed $p$, it is well known that the sample eigenvalues and eigenvectors are consistent estimators of the corresponding population eigenvalues and eigenvectors anderson. Under the “large $p$, large $n$” framework, however, the consistency is not guaranteed. The following two lemmas summarize and extend some known convergence results.

Let $p/n\rightarrow\gamma\geq 0$ as $n\rightarrow\infty$.

where $k$ is the number of $\lambda_{v}$ greater than $1+\sqrt{\gamma}$, and $\rho(x)=x(1+\gamma/(x-1))$.

The result in (ii) is due to Baik and Silverstein Baik2006, while the proof of (i) can be found in Section 6.3. The result in (i) shows that when $\gamma=0$, the sample eigenvalues converge to the corresponding population eigenvalues, which is consistent with the classical PC result where $p$ is fixed. The result in (ii) shows that for any nonzero $\gamma$, $d_{v}$ is no longer a consistent estimator of $\lambda_{v}$. However, a consistent estimator of $\lambda_{v}$ can be constructed from (2). Define

Then $\rho^{-1}(d_{v})$ is a consistent estimator of $\lambda_{v}$ when $\lambda_{v}>1+\sqrt{\gamma}$. Furthermore, Baik, Ben Arous and Péché baik2005ptl have shown the $\sqrt{n}$-consistency of $d_{v}$ to $\rho(\lambda_{v})$, and Bai and Yao bai2008clt have shown that $d_{v}$ is asymptotically normal.

Suppose $p/n\rightarrow\gamma\geq 0$ as $n\rightarrow\infty$. Let $\langle\cdot,\cdot\rangle$ be an inner product between two vectors. Under the assumption of multiplicity one:

if $0<\gamma<1$, and the $z_{ij}$’s follow the standard normal distribution, then

removing the normal assumption on the $z_{ij}$’s, the following weaker convergence result holds for all $\gamma\geq 0$:

Here $\phi(x)=\sqrt{(1-\frac{\gamma}{(x-1)^{2}})/(1+\frac{\gamma}{x-1})}$.

The inner product between unit vectors is the cosine angle between these two. Thus, Lemma 2 shows the convergence of the angle between population and sample eigenvectors. For (i), Paul paul2007 proved it for $\gamma<1$; while Nadler nadler2008fsa obtained the same conclusion for $\gamma>0$ using the matrix perturbation approach under the Gaussian random noise model. We relax the Gaussian assumption on $z$ and prove (ii) for $\gamma\geq 0$ in Section 6.4. The result of (ii) is general enough for the application of PCA to, for example, genome-wide association mapping, where each entry of $\mathbf{X}$ is a standardized variable of SNP genotypes, which are typically coded as $\{0,1,2\}$, corresponding to discrete genotypes.

3 Sample and predicted PC scores

In this section, we first discuss convergence of the sample PC scores, which forms the basis for the investigation of the shrinkage phenomenon of the predicted PC scores. For the sample PC scores, we have the following theorem.

Before we formally derive the asymptotic shrinkage factor for the predicted PC scores, we first describe in mathematical terms the shrinkage phenomenon that was demonstrated in the Introduction. Note that the first population eigenvector $\mathbf{e}_{1}$ satisfies

for a random vector $\mathbf{x}$ that follows the same distribution of the $\mathbf{x}_{j}$’s. For the data matrix $\mathbf{X}$, its first sample eigenvector $\mathbf{u}_{1}$ satisfies

Since the $\mathbf{u}_{1}^{T}\mathbf{x}_{j}$’s $(j=1,\ldots,n)$ follow the same distribution,

which demonstrates the shrinkage feature of the predicted PC scores. The amount of the shrinkage, or the asymptotic shrinkage factor, is given by the following theorem.

Suppose $p/n\rightarrow\gamma\geq 0$ as $n\rightarrow\infty$, $\lambda_{v}>1+\sqrt{\gamma}$. Under the multiplicity one assumption,

where $p_{vj}$ is the $j$th element of $\mathbf{p}_{v}$.

The proof is given in Section 6.8. We call $(\lambda_{v}-1)/(\lambda_{v}+\gamma-1)$, the (asymptotic) shrinkage factor for a new subject. As shown, the shrinkage factor is smaller than $1$ if $\gamma>0$. Quite sensibly, it is a decreasing function of $\gamma$ and an increasing function of $\lambda_{v}$. The bias of the predicted PC score can be potentially large for those high-dimensional data where $p$ is substantially greater than $n$, and/or for the data with relatively minor underlying structures where $\lambda_{v}$ is small.

4 Rescaling of sample eigenvalues

The previous theorems are based on the assumption that all except the top $m$ eigenvalues are equal to 1. Even under the spiked eigenvalue model, some rescaling of the sample eigenvalues may be necessary with real data.

For a given data, let its ordered population eigenvalues $\bm{\Lambda}^{*}=\{\zeta\lambda_{1},\ldots,\zeta\lambda_{m}$, $\zeta,\ldots,\zeta\}$, where $\zeta\neq 1$, and its corresponding sample eigenvalues $\mathbf{D}^{*}=\{d_{1}^{*},\ldots,\penalty d_{n}^{*}\}$. We can show that (4), (8) and (5) still hold under such circumstances. However, $\rho^{-1}(d_{v}^{*})$ is no longer a consistent estimator of $\lambda_{v}$, because

To address this issue, Baik and Silverstein Baik2006 have proposed a simple approach to estimate $\zeta$. In their method, the top significant large sample eigenvalues are first separated from the other grouped sample eigenvalues. Then $\zeta$ is estimated as the ratio between the average of the grouped sample eigenvalues and the mean determined by the Marčenko–Pastur law marvcenko1967. To separate the eigenvalues, they have suggested to use a screeplot of the percent variance versus component number. However, for real data, we may not be able to clearly separate the sample eigenvalues in such a manner and readily apply the approach. Thus, we need an automated method which does not require a clear separation of the sample eigenvalues.

The expectation of the sum of the sample eigenvalues when $\zeta=1$ is

Thus, the sum of the rescaled eigenvalues is expected to be close to $(\sum_{v=1}^{m}\lambda_{v}+p-m)$. Let $r_{v}=d_{v}^{*}/(\sum_{v=1}^{p}d_{v}^{*})$ and $\hat{d}_{v}$ be a properly rescaled eigenvalue, then $\hat{d}_{v}$ should be very close to $r_{v}(\sum_{v=1}^{m}\lambda_{v}+p-m)$. Note that ${p}/(\sum_{v=1}^{m}\lambda_{v}+p-\penalty m)\rightarrow 1$ for fixed $m$ and $\lambda_{v}$. Thus, $pr_{v}$ is a properly adjusted eigenvalue. However, for finite $n$ and $p$, the difference between $p$ and $(\sum_{v=1}^{m}\lambda_{v}+p-m)$ can be substantial, especially when the first several $\lambda_{v}$’s are considerably larger than $1$. To reduce this difference, we propose a novel method which iteratively estimates the $(\sum_{v=1}^{m}\lambda_{v}+p-m)$ and $\hat{d}_{v}$.

2. For the $l$th iteration, set $\hat{\lambda}_{v,l}=\rho^{-1}(\hat{d}_{v,l-1})$ for $\hat{d}_{v,l-1}>(1+\sqrt{\gamma})^{2}$, and $\hat{\lambda}_{v,l}=1$ for $\hat{d}_{v,l-1}\leq(1+\sqrt{\gamma})^{2}$. Define $k_{l}$ as the number of $\hat{\lambda}_{v,l}$’s that are greater than 1, and let

3. If $\sum_{v=1}^{k_{l}}\hat{\lambda}_{v,l}+p-k_{l}$ converges, let

The consistency of $\hat{d}_{v}$ to $\rho(\lambda_{v})$ is shown in the following theorem.

Let $\hat{d}_{v}$ be the rescaled sample eigenvalue from the proposed algorithm. Then, for $\lambda_{v}>1+\sqrt{\gamma}$ with multiplicity one,

Since $\rho^{-1}(\hat{d}_{v})\stackrel{{\scriptstyle p}}{{\rightarrow}}\lambda_{v}$, $\phi(\rho^{-1}(\hat{d}_{v}))^{2}$ is a consistent estimator of $\phi(\lambda_{v})^{2}$. Combining this fact with Theorems 1 and 2, we can obtain the bias-adjusted PC score $q_{v}^{*}$

Simulation

First, we applied our bias-adjustment process to the simulated data described in the Introduction. Our estimated asymptotic shrinkage factors are 0.465 and 0.329 for the first and second PC scores, respectively. The scatter plot of the top two bias-adjusted PC scores is given in Figure 2. After the bias adjustment,

the predicted PC scores of the test data are comparable to those of the training data. This indicates that our method is effective in correcting for the shrinkage bias.

Next, we conducted a new simulation to check the accuracy of our estimators. For the $j$th sample ($j=1,\ldots,n$), its $i$th variable was generated as

where $\lambda_{1}>\lambda_{2}>1$ and $z_{ij}\sim N(0,2^{2})$. Under this setting, $\lambda_{1}$ and $\lambda_{2}$ are the first and the second population eigenvalues. The first and second population eigenvectors are $e_{1}=\{1,0,\ldots,0\}$ and $e_{2}=\{0,1,0,\ldots,0\}$, respectively. We set the standard deviation of $z_{ij}$ to 2 instead of 1, which allows us to test whether the rescaling procedure works properly. We tried different values of $\gamma$ and $n$, but set $\lambda_{1}$ and $\lambda_{2}$ to $4(1+\sqrt{\gamma})$ and $2(1+\sqrt{\gamma})$, respectively.

We split the simulated samples into test and training sets, each with $n$ samples. We first estimated the asymptotic shrinkage factor based on the training samples. We then calculated the predicted PC scores on the test samples. To assess the accuracy of shrinkage factor estimator for each PC, we empirically estimated the shrinkage factor by the ratio of the mean predicted PC scores of the test samples to the mean PC scores of the training samples. That is, for the $v$th PC, the empirical shrinkage factor is estimated by $\sqrt{\sum_{i=1}^{n}q_{vi}^{2}/\sum_{k=1}^{n}p_{vk}^{2}}$. On the training samples, we also estimated the empirical angle between the sample and (known) population eigenvectors, as well as the empirical angle between PC scores from sample and population eigenvectors. The asymptotic theoretical estimates were also calculated. Tables 1 and 2 summarize the simulation results. Our asymptotic estimators provide accurate estimates for the angles and the shrinkage factor.

Finally, we conducted simulation to demonstrate an application of the bias-adjusted PC scores in PC regression. PC regression has been widely used in microarray gene-expression studies bovelstad2007. In this simulation, we let $p=5000$, and our set up is very similar to the first simulation of Bair et al. bair2006prediction. Let $x_{ij}$ denote the gene expression level of the $i$th gene for the $j$th subject. We generated each $x_{ij}$ according to

where $n$ is the number of samples, $g$ is the number of genes that are differentially expressed and associated with the phenotype, $\varepsilon\sim N(0,2^{2})$ and $\varepsilon_{y}\sim N(0,1)$. A total of eight different combinations of $n$ and $g$ were simulated. For the training data, we fit the PC regression with the first PC as the covariate and computed the mean square error (MSE). For the test samples with the same configuration of the training samples, we applied the PC model built on the training data to predict the phenotypes using the unadjusted and adjusted PC scores. The results are presented in Table 3. We see that the MSE of the test

set without bias adjustment is appreciably higher than that of the test set with bias adjustment, and the MSE of the test set with bias adjustment is comparable with the MSE of the training set.

Real data example

Here, we demonstrate that the shrinkage phenomenon appears in real data, and can be adjusted by our method. For this purpose, genetic data on samples from unrelated individuals in the Phase 3 HapMap study (http://hapmap.ncbi.nlm.nih.gov/) were used. HapMap is a dense genotyping study designed to elucidate population genetic differences. The genetic data are discrete, assuming the values 0, 1 or 2 at each genomic marker (also known as SNPs) for each individual. Data from CEU individuals (of northern and western European ancestry) were compared with data from TSI individuals (Toscani individuals from Italy, representing southern European ancestry).

Some initial data trimming steps are standard in genetic analysis. We first removed apparently related samples, and removed genomic markers with more than a 10% missing rate, and those with frequency less than 0.01 for the minor genetic allele. To avoid spurious PC results, we further pruned out SNPs that are in high linkage disequlibrium (LD) fellay2007. Lastly, we excluded $7$ samples with PC scores greater than 6 standard deviations away from the mean of at least one of the top significant PCs [i.e., with Tracy–Widom (TW) Test $p$-value $<$0.01] price2006pca; patterson2006psa. The final dataset contained 178 samples (101 CEU, 77 TSI) and 100,183 markers. We mean-centered and variance-standardized the genotypes for each marker price2006pca. The screeplot of the sample eigenvalues is presented in Figure 3. The first eigenvalue is substantially larger than the rest of the eigenvalues, although the TW test actually identifies two significant PCs. Figure 3 suggests that our data approximately satisfies the spiked eigenvalue assumption.

We estimated the asymptotic shrinkage factor and compared it with the following jackknife-based shrinkage factor estimate. For the first PC, we first computed the scores of all samples. Next, we removed one sample at a time and computed the (unadjusted) predicted PC score. We then calculated the jackknife estimate as the square root of the ratio of the means of the sample PC score and the predicted PC score. The jackknife shrinkage factor estimate is $0.319$, which is close to our asymptotic estimate $0.325$. Figure 4 shows the PC scores from the

whole sample, the predicted PC score of an illustrative excluded sample, and its bias-adjusted predicted score. Clearly, the predicted PC score without adjustment is very biased toward zero, while the bias-adjusted PC score is not.

Discussion and conclusions

In this paper, we have identified and explored the shrinkage phenomenon of the predicted PC scores, and have developed a novel method to adjust these quantities. We also have constructed the asymptotic estimator of correlation coefficient between PC scores from population eigenvectors and sample eigenvectors. In simulation experiments and real data analysis, we have demonstrated the accuracy of our estimates, and the capability to increase prediction accuracy in PC regression by adopting shrinkage bias adjustment. For achieving these, we consider asymptotics in the large $p$, large $n$ framework, under the spiked population model.

We believe that this asymptotic regime applies well to many high-dimensional datasets. It is not, however, the only model paradigm applied to such data. For example, the large $p$ small $n$ paradigm hall2005geometric; ahn2007high, which assumes $p/n\rightarrow\infty$, has also been explored. Under this assumption, Jung and Marron sungkyu have shown that the consistency and the strong inconsistency of the sample eigenvectors to population eigenvectors depend on whether $p$ increases at a slower or faster rate than $\lambda_{v}$. It may be argued that for real data where $p/n$ is “large,” we should follow the paradigm of Hall, Marron and Neeman hall2005geometric, Ahn et al. ahn2007high. However, for any real study, it is unclear how to test whether $p$ increases at a faster rate than $\lambda_{v}$, or vice versa, making the application of Hall, Marron and Neeman hall2005geometric, Ahn et al. ahn2007high difficult in practice. Furthermore, the scenario where $p$ and $\lambda_{v}$ grow at the same rate is scientifically more interesting, for which we are aware of no theoretical results. In contrast, our asymptotic results can be straightforwardly applied. Further, our simulation results indicate that for $p/n$ as large as 500, our asymptotic results still hold well. We believe that the approach we describe here applies to many datasets.

Although the results from the spiked model are useful, it is likely that observed data has more structure than allowed by the model. Recently, several methods have been suggested to estimate population eigenvalues under more general scenarios elkaroui2008sel; rao2008. However, no analogous results are available for the eigenvectors. In data analysis, jackknife estimators, as demonstrated in the real data analysis section, can be used. However, resampling approaches are very computationally intensive, and it remains of interest to establish the asymptotic behavior of eigenvectors in a variety of situations.

We note that inconsistency of the sample eigenvectors does not necessarily imply poor performance of PCA. For example, PCA has been successfully applied in genome-wide association studies for accurate estimation of ethnicity price2006pca, and in PC regression for microarrays ma2006additive. However, for any individual study we cannot rule out the possibility of poor performance of the PC analysis. Our asymptotic result on the correlation coefficient between PC scores from sample and population eigenvectors provides us a measure to quantify the performance of PC analysis.

For the CEU/TSI data, SNP pruning was applied to adjust for strong LD among adjacent SNPs. Such SNP pruning is a common practice in the analysis of GWAS data, and has been implemented in the popular GWAS analysis software Plink purcell2007pts. The primary goal of SNP pruning is to avoid spurious PC results unrelated to population substructures. Technically, our approach does not rely on any independence assumption of the SNPs. However, strong local correlation may affect eigenvalues considerably. Thus, the value in SNP pruning may be viewed as helping the data better accord with the assumptions of the spiked population model. From the CEU/TSI data and our experience in other GWAS data, we have found that the most common pruning procedure implemented in Plink is sufficient for us to then apply our methods.

Proofs

Note that $\mathbf{E}\bm{\Lambda}^{1/2}\mathbf{Z}\mathbf{Z}^{T}\bm{\Lambda}^{1/2}\mathbf{E}^{T}$ and $\bm{\Lambda}^{1/2}\mathbf{Z}\mathbf{Z}^{T}\bm{\Lambda}^{1/2}$ have the same eigenvalues, and $\mathbf{E}^{T}\mathbf{U}$ is the eigenvector matrix of $\bm{\Lambda}^{1/2}\mathbf{Z}\mathbf{Z}^{T}\bm{\Lambda}^{1/2}$. Since eigenvalues and angles between sample and population eigenvectors are what we concerned about, without loss of generality (WLOG), in the sequel, we assume $\bm{\Lambda}$ to be the population covariance matrix.

We largely follow notation in Paul paul2007. We denote $\lambda_{v}(\mathbf{S})$ as the $v$th largest eigenvalue of $\mathbf{S}$. Let suffice $A$ represent the first $m$ coordinates and $B$ represent the remaining coordinates. Then we can partition $\mathbf{S}$ into

We similarly partition the $v$th eigenvector $\mathbf{u}_{v}^{T}$ into $(\mathbf{u}_{A,v},\mathbf{u}_{B,v})$ and $\mathbf{Z}^{T}$ into $[\mathbf{Z}_{A}^{T},\mathbf{Z}_{B}^{T}]$. Define $R_{v}$ as $\|\mathbf{u}_{B,v}\|$ and let $\mathbf{a}_{v}=\mathbf{u}_{A,v}/\sqrt{1-R_{v}^{2}}$, then we get $\|\mathbf{a}_{v}\|=1$.

Applying singular value decomposition (SVD) to $\mathbf{Z}_{B}/\sqrt{n}$, we get

where $\mathbf{M}=\operatorname{diag}(\mu_{1},\ldots,\mu_{p-m})$ is a $(p-m)\times(p-m)$ diagonal matrix of ordered eigenvalues of $\mathbf{S}_{BB}$, $\mathbf{V}$ is a $(p-m)\times(p-m)$ orthogonal matrix and $\mathbf{H}$ is an $n\times(p-m)$ matrix. For $n\geq p-m$, $\mathbf{H}$ has full rank orthogonal columns. When $n<p-m$, $\mathbf{H}$ has more columns than rows, hence it does not have full rank orthogonal columns. For the later case, we make $\mathbf{H}=[\mathbf{H}_{n},0]$ where $\mathbf{H}_{n}$ is an $n\times n$ orthogonal matrix.

2 Propositions

We introduce two propositions for later use. The proofs of the two propositions can be found in Sections 6.5 and 6.6.

Suppose $\mathbf{Y}$ is an $n\times m$ matrix with fixed $m$ and each entry of $\mathbf{Y}$ is i.i.d. random variable which satisfies the moment condition of $z_{ij}$ in Assumption 1. Let $\mathbf{C}$ be an $n\times n$ symmetric nonnegative definite random matrix and independent of $\mathbf{Y}$. Further, assume $\|\mathbf{C}\|=O(1)$. Then

as $n\rightarrow\infty$, where $F_{\gamma}(x)$ is a distribution function of Marčenko–Pastur law with parameter $\gamma$ marvcenko1967.

3 Proof of part (i) of Lemma 1

3.2 When p→∞p\rightarrow\infty

Similarly by the interlacing inequality, we get

Part (i) of Lemma 1 follows by (10) and (11).

4 Proof of part (ii) of Lemma 2

Our proof of Lemma 2(ii) closely follows the arguments in Paul paul2007. From paul2007, it can be shown that

where $\mathbf{e}_{A,v}$ is a vector of the first $m$ coordinates of the $v$th population eigenvector $\mathbf{e}_{v}$, $\rho_{v}$ is $\lambda_{v}(1+\frac{\gamma}{\lambda_{v}-1})$ and $\mathbf{z}_{Av}$ is a

vector of $v$th row of $\mathbf{Z}_{A}$. The proofs can be found in Section 6.4.3. Note that $\mathbf{e}_{v}$ is a vector with $1$ in its $v$th coordinate and $0$ elsewhere. WLOG, we assume that $\langle\mathbf{e}_{v},\mathbf{u}_{v}\rangle\geq 0$. Since $\langle\mathbf{e}_{v},\mathbf{u}_{v}\rangle=\sqrt{1-R_{v}^{2}}\langle\mathbf{e}_{A,v},\mathbf{a}_{v}\rangle$, $\langle\mathbf{e}_{v},\mathbf{u}_{v}\rangle\stackrel{{\scriptstyle p}}{{\rightarrow}}\sqrt{1-R_{v}^{2}}$. By (13) and (15), we can show that

It concludes the proof of the first part of Lemma 2(ii).

4.2 When 1<λv≤1+γ1<\lambda_{v}\leq 1+\sqrt{\gamma}

Here, we only need to consider $\gamma>0$ because no eigenvalue satisfies this condition when $\gamma=0$. We first show that $R_{v}\stackrel{{\scriptstyle p}}{{\rightarrow}}1$, which implies $\mathbf{u}_{A,v}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$, hence $\langle\mathbf{e}_{v},\mathbf{u}_{v}\rangle\stackrel{{\scriptstyle p}}{{\rightarrow}}0$. For any $\varepsilon>0$ and $x\geq 0$, define

where $a=(1-\sqrt{\gamma})^{2}$ and $b=(1+\sqrt{\gamma})^{2}$. Since (21) equals $\infty$ for any $a\leq\rho_{v}\leq b$, we conclude that

Therefore $R_{v}\stackrel{{\scriptstyle p}}{{\rightarrow}}1$, which proves the second part of Lemma 2(ii).

4.3 Proof of (14) and (15)

With the exactly same argument of paul2007, it can be shown that

where $\mathbf{r}_{v}=-(1-\langle\mathbf{e}_{A,v},\mathbf{a}_{v}\rangle)\mathbf{e}_{A,v}-\mathcal{R}_{v}\mathcal{D}_{v}(\mathbf{a}_{v}-\mathbf{e}_{A,v})+(d_{v}-\rho_{v})\mathcal{R}_{v}(\mathbf{a}_{v}-\mathbf{e}_{A,v})$. By Lemma 1 of paul2005, $r_{v}=o_{p}(1)$, if $\alpha_{v}=o_{p}(1)$ and $\beta_{v}=o_{p}(1)$.

When $\gamma=0$, $\mathbf{S}_{AA}-(\rho_{v}/\lambda_{v})\bm{\Lambda}_{A}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$ and the remainder of $\mathcal{D}_{v}$ is

When $\gamma>0$, $\mathcal{D}_{v}$ can be written as

The first term of the right-hand side is $o_{p}(1)$ by the weak law of large number. The second and third terms are $o_{p}(1)$ by Propositions 1 and 2. For the fourth term, $\rho_{v}-d_{v}=o_{p}(1)$ and its remainder part is $O_{p}(1)$. Therefore, $\mathcal{D}_{v}=o_{p}(1)$. By combining the above results and $\mathcal{R}_{v}=O_{p}(1)$ plus $d_{v}-\rho_{v}=o_{p}(1)$, we prove (14).

5 Proof of Proposition 1

Let $\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{n}$ be the ordered eigenvalues of $\mathbf{C}$, and $c_{ij}$ be the $(i,j)$th element of $\mathbf{C}$. Suppose $\mathbf{y}_{s}$ is the $s$th column of $\mathbf{Y}$, and $y_{ij}$ is the $(i,j)$th element of $\mathbf{Y}$. We further define $\psi(s,s)=\frac{1}{n}\mathbf{y}_{s}^{T}\mathbf{C}\mathbf{y}_{s}-\frac{1}{n}\operatorname{trace}(\mathbf{C})$ and $\psi(s,t)=\frac{1}{n}\mathbf{y}_{s}^{T}\mathbf{C}\mathbf{y}_{t}$ for $s\neq t$. The conditional mean of $\psi(s,s)$ given $\mathbf{C}$ is

Thus, $E(\psi(s,s))=E(E(\psi(s,s)|\mathbf{C}))=E(0)=0$.

Next, the conditional variance of $\psi(s,s)$ given $\mathbf{C}$ is

where $\alpha=\max(1,E(y_{is}^{4})-1)$. Since $\|\mathbf{C}\|=O(1)$, $\mu_{i}^{2}\leq\|\mathbf{C}\|^{2}=O(1)$. Therefore, $\operatorname{Var}(\psi(s,s)|\mathbf{C})\leq O(1/n)$ and $\operatorname{Var}(\psi(s,s))=\operatorname{Var}(E(\psi(s,s)|\mathbf{C}))+\penalty E(\operatorname{Var}(\psi(s,s)|\mathbf{C}))\leq 0+O(1/n)\rightarrow 0$ as $n\rightarrow\infty$. By the Chebyshev inequality, we can conclude that

We can similarly show $\psi(s,t)\stackrel{{\scriptstyle p}}{{\rightarrow}}0$, which we omit here.

6 Proof of Proposition 2

We show that both (a) and (b) converge to 0 in probability.

(b): Let $F_{p-m}$ be an empirical spectral distribution of $\mathbf{S}_{BB}$, then

and $\int{f(x)\,dF_{n}(x)}\stackrel{{\scriptstyle p}}{{\rightarrow}}\int{f(x)\,dF_{\gamma}(x)}$ marvcenko1967; bai1999msa. Thus,

Combining (a) and (b), we finish the proof.

7 Proof of Theorem 1

8 Proof of Theorem 2

First, we show the square of the denominator converges to $\rho(\lambda_{v})$. Since $p_{vj}=\mathbf{u}_{v}^{T}\mathbf{x}_{j}$, and $E(p_{vi}^{2})=E(p_{vj}^{2})$ for $i\neq j$,

Next, we show the square of numerator converges to $\phi(\lambda_{v})^{2}(\lambda_{v}-1)+1$. Define $\mathbf{u}_{v}^{\bot}:=\frac{1}{\sqrt{1-(\mathbf{u}_{v}^{T}\mathbf{e}_{v})^{2}}}(I-\mathbf{e}_{v}\mathbf{e}_{v}^{T})\mathbf{u}_{v}$, then $\mathbf{u}_{v}$ can be expressed as

9 Proof of Theorem 3

Since $\rho^{-1}(pr_{v})\rightarrow\lambda_{v}$ for $v\leq k$, WLOG we assume that $k_{0}=k$, where $k$ is the number of $\lambda_{v}$ bigger than $1+\sqrt{\gamma}$. Set

The first and second partial derivatives of $h(x)$ are

Since $d_{v}\rightarrow\rho(\lambda_{v})$ for $v\leq k$,

From the facts that $h(x)$ is a continuous concave function, $\omega>p$, and $h(p)>0$, we conclude that

for $v\leq k$, which concludes the proof.

References