arXiv:1607.08647v1 [] 28 Jul 2016

ASYMPTOTIC PROPERTIES OF PRINCIPAL
COMPONENT ANALYSIS AND SHRINKAGE-BIAS
ADJUSTMENT UNDER THE GENERALIZED SPIKED
POPULATION MODEL
By Rounak Dey∗ and Seunggeun Lee∗
Department of Biostatistics, University of Michigan
With the development of high-throughput technologies, principal component analysis (PCA) in the high-dimensional regime is of
great interest. Most of the existing theoretical and methodological
results for high-dimensional PCA are based on the spiked population model in which all the population eigenvalues are equal except
for a few large ones. Due to the presence of local correlation among
features, however, this assumption may not be satisfied in many realworld datasets. To address this issue, we investigated the asymptotic
behaviors of PCA under the generalized spiked population model.
Based on the theoretical results, we proposed a series of methods for
the consistent estimation of population eigenvalues, angles between
the sample and population eigenvectors, correlation coefficients between the sample and population principal component (PC) scores,
and the shrinkage bias adjustment for the predicted PC scores. Using numerical experiments and real data examples from the genetics
literature, we showed that our methods can greatly reduce bias and
improve prediction accuracy.

1. Introduction. Principal component analysis (PCA) is a very popular tool for analyzing high-dimensional biomedical data, where the number
of features (p) is often substantially larger than the number of observations
(n). PCA is widely used to adjust for population stratification in genomewide association studies (Price et al., 2006) and to identify overall expression patterns in transcriptome analysis (Storey et al., 2005). However, the
asymptotic properties of PCA in high-dimensional data are profoundly different from the properties in low-dimensional (p finite, n → ∞) settings. In
high-dimensional settings, the sample eigenvalues and eigenvectors are not
consistent estimators of the population eigenvalues and eigenvectors (Johnstone and Lu, 2009; Paul, 2007), and the predicted principal component
(PC) scores based on the sample eigenvectors can be systematically biased
toward zero (Lee, Zou and Wright, 2010).
∗

Supported by NIH Grants R00HL113164 and R01HG008773
Primary 62H25; secondary 15A18
Keywords and phrases: PCA, PC Scores, Random Matrix, High-dimensional data

1

2

R. DEY AND S. LEE

There has been extensive effort to investigate the asymptotic behaviors
of PCA in high-dimensional settings. To provide a statistical framework
for PCA in these settings, Johnston introduced a spiked population model,
which assumes that all the eigenvalues are equal except for finitely many
large ones (called the spikes). A spiked population covariance matrix is basically a finite rank perturbation of a scalar multiple of the identity matrix.
A typical example of a spiked population with two spikes is shown in Figure
1a. This two-spike eigenvalue structure arises if the population consists of
three sub-populations, and the features are largely independent with equal
variances. Under this model, convergence of sample eigenvalues, eigenvectors and PC scores have been extensively studied (Johnstone, 2001; Baik
and Silverstein, 2006; Paul, 2007; Lee, Zou and Wright, 2010).
In many biomedical data, however, the assumption of the equality of
non-spiked eigenvalues can be violated due to the presence of local correlation among features. In genome-wide association studies, for example,
the genetic variants are locally correlated due to linkage disequilibrium. In
gene-expression data, since genes in the same pathway are often expressed
together, their expression measurements are often correlated. These local
correlations can cause substantial differences in non-spiked eigenvalues. To
illustrate this phenomenon, we obtained eigenvalues with an autoregressive
within-group correlation structure rather than the independent structure
of the previous example. Figure 1b shows that the equality assumption is
clearly violated. Thus, if methods developed under the equality assumption
are applied to these types of data, we will obtain biased results.
The generalized spiked population model (Bai and Yao, 2012) has been
proposed to address this problem. The condition that the non-spikes have to
be equal is removed in this generalization. In this model the set of population
eigenvalues consists of finitely many large eigenvalues called the generalized
spikes, which are well separated from infinitely many small eigenvalues. Although the generalized spiked population model has a great potential to
provide more accurate inference in high-dimensional biomedical data, only
limited literature is available on the asymptotic properties of PCA under
this model and their application to real data. Bai and Yao (2012) and Ding
(2015) provided results regarding convergence of eigenvalues and eigenvectors. However, their work remained largely theoretical. Moreover, to the best
of our knowledge, no method has been developed for estimating the correlations between the sample and population PC scores, and adjusting biases
in the predicted PC scores under the generalized spiked population model.
In this paper, we systematically investigate the asymptotic behaviors of
PCA under the generalized spiked population model, and develop methods

PCA IN HIGH-DIMENSIONAL DATA

3

to estimate the population eigenvalues and adjust for the bias in the predicted PC scores. We first propose two different approaches to consistently
estimate the population eigenvalues, the angles between the sample and
population eigenvectors, and the correlation coefficients between the sample
and population PC scores. We compare these two methods and show the
asymptotic equivalence of the estimators across them. Finally, we propose a
method to reduce the bias in the predicted PC scores based on the estimated
population eigenvalues.
The paper is organized as follows. We begin in Section 2 by providing
the definition of the generalized spiked population model and present existing theoretical results. We develop our methods to consistently estimate the
population spikes in Section 3. In Section 4, we construct consistent estimators of the angles between the sample and population eigenvectors, and the
correlation coefficients between the sample and population PC scores. We
also propose the bias-reduction technique for the predicted PC scores. Section 5 presents the algorithm (El Karoui, 2008) to estimate the population
limiting spectral distribution and the non-spiked eigenvalues. In Section 6,
we present results from simulation studies and an example from the Hapmap
project to demonstrate the improved performance of our method over the
existing one. Finally, we conclude the paper with a discussion.
2. Generalized spiked population model. In order to formally define generalized spiked population model, we require the concept of spectral
distribution. In random matrix literature, it is natural to associate a probability measure to the set of eigenvalues as the dimension (p) goes to ∞.
More explicitly, if a hermitian matrix Σp has eigenvalues λ1 , λ2 , . . . , λp , we
can define the empirical spectral distribution (ESD) of Σp to be Hp based
on the probability measure
p

1X
dHp (x) =
δλi (x),
p
i=1

where δλi (x) is unity when x = λi , and otherwise zero. Now, for a sequence
{Σp } of covariance matrices, if the corresponding sequence {Hp } of ESDs
converge weakly to a non-random probability distribution H as p → ∞,
then we define H as the limiting spectral distribution (LSD) of the sequence
{Σp }.
The generalized spiked population model (Bai and Yao, 2012) is defined as
follows. Suppose, Hp is the ESD corresponding to the population covariance
matrix Σp and it converges weakly to a non-random probability distribution
H. Let ΓH be the support of H and d(x, A) := inf y∈A |x − y| be the distance

4

R. DEY AND S. LEE

metric from a point x to a set A. Then the set of eigenvalues of Σp comprises
of two subsets of eigenvalues where,
• Non-spike: an eigenvalue β is called a non-spike if d(β, ΓH ) = p → 0,
• Generalized spike: an eigenvalue α is called a generalized spike if ∃δ > 0
such that lim inf d(α, ΓH ) > δ, and α > sup(ΓH ).
It is obvious from the definition that the generalized spikes are measure zero
points of the population LSD. For Johnstone’s spiked population model
(Johnstone, 2001), the population LSD is H = δ{1} , indicating ΓH = {1}.
From the definition above, all eigenvalues larger than one are spikes. Hence,
Johnstone’s spiked population model is a special case of the generalized
spiked population model.
Suppose that the population covariance matrix Σp has eigenvalues λ1 ≥
λ2 ≥ . . . ≥ λp , and the sample covariance matrix Sp = X T X/n has eigenvalues d1 ≥ d2 ≥ . . . ≥ dp , where X is an n × p data matrix. Further, we
will assume the following throughout the paper:
(A) p → ∞, n → ∞, p/n → γ < ∞.
(B) The population eigenvalues follow the generalized spiked population
(GSP) model. The population ESD Hp converges weakly to a nonrandom probability distribution H with support ΓH . Moreover, the
sequence {kΣp k} of spectral norms is bounded.
1/2
(C) The n × p data matrix X = Y Σp where Y is an n × p random matrix
with i.i.d. elements such that E(Yij ) = 0, E(|Yij |2 ) = 1, E(|Yij |4 ) < ∞.
From the Marc̆enko–Pastur theorem (Marc̆enko and Pastur, 1967), the
sample ESD Fp converges weakly to a non-random probability distribution
F with support ΓF . For α > sup ΓH and x > 0, we define the following two
functions
Z
λdH(λ)
x
(2.1)
ψ(α) := α + γα
, fF (x) :=
R τ dF (τ ) .
α−λ
1+γ
x−τ

The following result by Bai and Yao (2012) provides the almost sure limits of
the sample eigenvalues corresponding to the population generalized spikes.
Result 1 (Bai and Yao (2012)). Suppose assumptions (A)–(C) hold. Let
λk be a generalized spike of multiplicity one and the corresponding sample
eigenvalue is dk . Moreover, let ψ 0 denote the first derivative of the function
ψ. Then,
• If ψ 0 (λk ) > 0, then the sample eigenvalue dk converges almost surely
to ψ(λk ), i.e.
a.s.
|dk − ψ(λk )| −−→ 0.

PCA IN HIGH-DIMENSIONAL DATA

5

• If ψ 0 (λk ) ≤ 0, then let (uk , vk ) ⊂ ΓcH be the maximal interval on which
ψ 0 > 0. Then the sample eigenvalue dk converges almost surely to
ψ(w) where w is a boundary of [uk , vk ] that is nearest to λk . If no such
interval exists, then dk converges almost surely to the γ th quantile of
F.
Since ψ 0 (α) is a strictly increasing function for α > sup ΓH , if a generalized spike λk is large enough such that ψ 0 (α) > 0, according to Result 1
the corresponding sample eigenvalue will converge almost surely to ψ(λk ).
However if the generalized spike lies close enough, i.e. ψ 0 (λk ) ≤ 0, to the set
of non-spikes then the convergence of the corresponding sample eigenvalue
is given by the second part of the result. We will denote a generalized spike
λk as a “distant spike” if ψ 0 (λk ) > 0, otherwise we will call it a “close spike”.
3. Consistent estimation of the generalized spikes. The following
theorem provides two different consistent estimators of the distant spikes.
Theorem 1. Let λk be a distant spike of multiplicity one and the corresponding sample eigenvalue is dk . If the assumptions (A)–(C) hold, then,
p

|ψ −1 (dk ) − λk | →
− 0,
where ψ −1 is the left inverse of ψ. Also,
p

|fF (dk ) − λk | →
− 0.
This theorem shows that for any distant spike λk we have two consistent
estimators ψ −1 (dk ) and fF (dk ). Notice that the function fF depends only
on the sample LSD which can be approximated by the sample ESD. Thus,
fF (dk ) can be approximated directly using the sample eigenvalues. More
explicitly, if there are m unique distant spikes in the population, fF (dk ) can
be closely approximated as
fF (dk ) ≈

1+

γ
p−m

dk
Pp

di
i=m+1 dk −di

.

In contrast, the ψ function depends on the population LSD which is unknown. We can estimate the ψ function using the algorithm described in
Section 5 and then find the inverse function ψ −1 using a Newton-Raphson
type algorithm. For the close spikes, no consistent estimator exists.

6

R. DEY AND S. LEE

4. Consistent estimators of the asymptotic shrinkage in predicting the PC scores. In this section, we investigate the convergence of
sample eigenvectors, PC scores, and shrinkage factors in predicting the PC
scores. Let ei and Ei to be the ith sample and population eigenvectors,
respectively. In addition to assumptions (A)–(C), we further assume that
the distant spikes are of multiplicity one. This assumption is to restrict the
dimension of the corresponding eigenspaces to one, as otherwise the angle
between sample and population eigenvectors, or shrinkage in predicted PC
scores cannot be well defined.
4.1. Angle between sample and population eigenvectors. We first present
the following theorem on the convergence of the quadratic forms of the
sample eigenvectors.
Theorem 2. Let λk be a distant spike of multiplicity one, and the assumptions (A)–(C) hold. Consider the following quadratic form
ηˆk = sT1 ek eTk s2 ,
where s1 and s2 are non-random vectors with uniformly bounded norm for
all p. Then,
a.s.
|ηˆk − ηk | −−→ 0,
where
ηk =

λk ψ 0 (λk ) T
s Ek Ek T s2
ψ(λk ) 1

Mestre (2008a) showed similar asymptotic properties of the quadratic
forms under the assumption that the number of spikes increases with the
dimension. Theorem 2 shows the convergence of the angle between sample
and population eigenvectors. Suppose s1 = s2 = Ek , and then,
ηˆk = EkT ek eTk Ek = hek , Ek i2 ,

ηk =

λk ψ 0 (λk )
.
ψ(λk )

Combining them, we can show
(4.1)

hek , Ek i2 −

λk ψ 0 (λk ) a.s.
−−→ 0.
ψ(λk )

Therefore, {λk ψ 0 (λk )/ψ(λk )}1/2 is a consistent estimator of the cosine of the
angle, i.e. the absolute value of the inner product, between the k th sample
and population eigenvectors. In order to obtain this estimator we first need
to estimate the ψ function using the algorithm described in Section 5.

7

PCA IN HIGH-DIMENSIONAL DATA

The following result by Ding (2015) provides another consistent estimator
for the angle between the k th sample and population eigenvectors. The proof
of the asymptotic equivalence of these two estimators is given in Appendix
A.
Result 2. Let λk be a distant spike of multiplicity one, and dk be the
corresponding sample eigenvalue. Assume that (A)–(C) hold. Define,


Z
τ dF (τ ) −1
gF (x) := 1 + γfF (x)
.
(x − τ )2
Then,

p

|hek , Ek i2 − gF (dk )| →
− 0.
Hence gF (dk )1/2 also works as a consistent estimator of |hek , Ek i|. Since
the function gF depends only on sample LSD, it can be approximated directly using sample eigenvalues. More explicitly, if there are m spikes in the
population, the function gF can be closely approximated as
p
γfF (dk ) X
di
gF (dk ) ≈ 1 +
p−m
(dk − di )2

"

#−1
.

i=m+1

The above equation can be used to estimate the angle between the sample
and population eigenvectors.
4.2. Correlation between sample and population PC scores. The sample
and population PC scores are the projections of the data on the sample and
population eigenvectors respectively. The correlation between them can be
perceived as a measure of accuracy of the PCA. The squared correlation can
also be interpreted as the proportion of variance in the population PC scores
that can be explained by corresponding sample PC scores. The following
theorem provides the consistent estimators of the correlation between the
sample and population PC scores corresponding to a distant spike.
Theorem 3. Suppose λk is a distant spike of multiplicity one, dk is
the corresponding sample eigenvalue, and the assumptions (A)–(C) hold.
Let the normalized k th population PC score is Pk = XEk /(nλk )1/2 and the
normalized k th sample PC score is pk = Xek /(ndk )1/2 . Then,
p

hPk , pk i2 − ψ 0 (λk ) →
− 0,

8

R. DEY AND S. LEE

and,
hPk , pk i2 −

dk gF (dk ) p
→
− 0,
fF (dk )

where the function gF is as defined in Result 2.
Since Pk and pk are normalized random vectors, the absolute value of the
inner product hPk , pk i is identical to the absolute value of their correlation
coefficient. Since correlation is scale invariant, this is also the correlation
between k th sample and population PC scores. Therefore we can consider
both ψ 0 (λk )1/2 and {dk gF (dk )/fF (dk )}1/2 to be consistent estimators of the
correlation between the k th sample and population PC scores.
4.3. Asymptotic shrinkage factor. Suppose λk is a distant spike. Let the
sample PC score for the j th observation xj be pkj = xTj ek , and the k th
predicted PC score for a new observation xnew be qk = xTnew ek . Then the
quantity ρk = limp→∞ {E(qk2 )/E(p2kj )}1/2 describes the asymptotic shrinkage in the k th predicted PC score for a new observation. As both pkj and
qk are centered, i.e. E(pkj ) = E(qk ) = 0, ρk represents the limiting ratio
of the standard deviations of the predicted PC scores and the sample PC
scores. Therefore, if we can estimate ρk , then the shrinkage bias in the k th
predicted PC scores can be easily adjusted by rescaling the predicted scores
by the factor ρ−1
k . The following theorem provides the consistent estimator
of the asymptotic shrinkage factor ρk .
k th

Theorem 4. Suppose λk is a distant spike of multiplicity one, dk is the
corresponding sample eigenvalue, and the assumptions (A)–(C) hold. Let pkj
and qk be as defined above. Then,
s
E(qk2 )
λk p
−
→
− 0.
2
E(pkj ) dk
This is a surprising result in which the asymptotic shrinkage factor is expressed as a simple ratio of the population and sample eigenvalues. Recall
that we already constructed the consistent estimators for population eigenvalues in the previous sections. Using these results, the asymptotic shrinkage
factor ρk can be consistently estimated by λ̂k /dk where λ̂k is any consistent
estimator of λk .
4.4. Comparison of the two different estimators. For each of the quantities discussed above, we proposed two asymptotically equivalent estimators.

PCA IN HIGH-DIMENSIONAL DATA

9

In terms of practical applications they have their own advantages and disadvantages. One of them can be approximated directly based only on the
sample eigenvalues, while the other one requires to estimate the LSD of the
population eigenvalues to obtain the ψ function. For ease of discourse we will
call the former “d-estimator” and the later “λ-estimator”. If the number of
spikes is known, estimating the d-estimator is computationally more efficient
than estimating the λ-estimator as it does not involve estimating the population LSD. However, by estimating the population LSD the λ-estimation
procedure can verify whether an estimated eigenvalue is actually a distant
spike by checking if ψ 0 > 0. Thus it can be used to estimate the number of
distant spikes when it is unknown (see Section 5). On the other hand, the
d-estimation procedure provides no information on the population LSD and
thus cannot distinguish among distant spikes, close spikes and non-spikes.
To summarize, when the number of spikes is known or we only want to estimate few of the largest eigenvalues which are known to be distant spikes,
then the d-estimation procedure has the advantage of a faster computation,
while the λ-estimation procedure is more useful when the number of spikes
is unknown or the distribution of the non-spikes is of interest.
4.5. Comparison of Generalized Spiked Population (GSP) model and Spiked
Population (SP) model. As mentioned before, the SP model (Johnstone,
2001) is a special case of the GSP model. It is easy to verify that when the
population eigenvalues follow the SP model, our consistent estimators for
the spiked eigenvalues, the angles between the eigenvectors, the correlation
coefficients between the PC scores and the shrinkage factors conform to the
consistent estimators derived by Baik and Silverstein (2006); Paul (2007);
Lee, Zou and Wright (2010). For an SP model where all the non-spikes are
equal to one, the LSD H is a degenerate distribution at one, and


γ
γ
ψ(α) = α 1 +
; ψ 0 (α) = 1 −
.
α−1
(α − 1)2
Now, ψ 0 (α) > 0 if and only if α > 1 + γ 1/2 . If α > 1 + γ 1/2 and d is the
corresponding sample eigenvalue, then the consistent estimator of α is given
by ψ −1 (d), and
γ
1 − (α−1)
2
αψ 0 (α)
=
;
γ
ψ(α)
1 + α−1

α
α−1
=
,
ψ(α)
α+γ−1

which show that all our results match with the results from Lee, Zou and
Wright (2010).

10

R. DEY AND S. LEE

It is of interest to investigate how closely methods developed under the
SP model can approximate the consistent estimators for the distant spikes
when the population eigenvalues actually follow a GSP model. Suppose the
population eigenvalues λ1 ≥ λ2 ≥ . . . ≥ λp follow the GSP model with
m distant spikes. The sample eigenvalues are d1 ≥ d2 ≥ . . . ≥ dp . Let
λk be a distant spike with multiplicity one, and the corresponding sample
eigenvalue is dk . Then according to Result 1, dk → ψ(λk ) almost surely.
From the definition of ψ,


Z
Z
λdH(λ)
λdH(λ)
ψ(λk ) = λk 1 + γ
= λk + γ
.
λk − λ
1 − λ/λk
If H is almost degenerate, i.e., the non-spikes are nearly identical, then,
(4.2)

ψ(λk ) ≈ λk +

γ λ̄
,
1 − λ̄/λk

R
where λ̄ = λdH(λ) is the mean of the population LSD which can be closely
approximated by the mean of the non-spikes. On the other hand, if the spike
λk is very large compared to all the non-spikes such that λ/λk ≈ 0 for any
λ ∈ ΓH , then
(4.3)

ψ(λk ) ≈ λk + γ λ̄.

Now, suppose instead of using the GSP assumption, we use the SP assumption to estimate the distant spikes. We assume that under the SP model
the population covariance matrix is scaled by a factor ζ and the population
eigenvalues are β1 ≥ β2 ≥ . . . ≥ βm > ζ = ζ = . . . ζ. If βk is the population
eigenvalue corresponding to dk , then dk → ψ(βk ) almost surely where,


ζ
ψ(βk ) = βk 1 + γ
βk − ζ
γζ
= βk +
.
1 − ζ/βk
Here ζ is estimated as the mean of the non-spikes as they are all assumed
to be equal to ζ. Notice that this expression is approximately equal to the
expression in (4.2) with βk = λk and ζ = λ̄. Therefore, the asymptotic
limit of dk under both the GSP and the SP model are approximately equal
when the non-spikes are nearly identical. On the other hand, when the spike
βk is very large compared to all the non-spikes such that ζ/βk ≈ 0, then
ψ(βk ) ≈ βk + γζ. In this case also, the asymptotic limit of dk under both
the GSP and the SP model are approximately equal with βk = λk and

11

PCA IN HIGH-DIMENSIONAL DATA

ζ = λ̄. Therefore if a generalized spike is very far away from the support of
the population LSD, then the estimate of the spike based on an SP model
will closely approximate the estimate based on a GSP model. However the
SP model will provide potentially biased estimates if the non-spikes are not
similar and the ratio between the largest non-spike and the spike of interest
is substantially larger than zero.
5. Estimation of the population LSD. The λ-estimators rely on
ψ, that is a function of the unknown population LSD H. To use the λestimators, it is thus required to estimate H. Using the Stieltjes transformation and the Marc̆enko–Pastur theorem, El Karoui (2008) developed a
general algorithm to estimate the population LSD from the sample ESD,
Fp . We propose to use Karoui’s method to estimate the population LSD H
and then use it to estimate ψ.
5.1. Karoui’s algorithm. Suppose vFp is the Stieltjes transformation of
the set of eigenvalues in the sample covariance matrix in which
n

vFp (z) =

1X 1
n
di − z
i=1

for any z ∈ C+ , C+ = {x ∈ C, Im(x) > 0}. According to the Marc̆enko–
Pastur theorem (Marc̆enko and Pastur, 1967), when assumptions (A)–(C)
hold, vFp converges pointwise almost surely to a non-random limit vF , which
uniquely satisfies the following equation

Z
vF (z) = − z − γ

λdH(λ)
1 + λvF (z)

−1
.

Karoui’s method first calculates vFp for a grid of values {zj }Jj=1 , and then
finds Ĥ as a solution to minimize the following objective function
J !

Z
1
p
λdH(λ)
,
Ĥ = argH min L
+ zj −
vFp (zj )
n
1 + λvFp (zj ) j=1
where L is any pre-defined convex loss function. In order to approximate
the integral inside of the loss function, the algorithm discretizes H in the
following way,
K
X
dH(λ) '
wk δtk (λ),
k=1

12

R. DEY AND S. LEE

P
where δtk (λ) = 1 if λ = tk and 0 otherwise, K
k=1 wk = 1 with wk > 0 for
all k, and {tk }K
is
a
grid
of
points
on
the
support
of H. This is basically
k=1
approximating H by a discrete distribution with support {tk }K
k=1 . Then the
integral is approximated by
Z

K

X
λdH(λ)
tk
'
,
wk
1 + λvF (z)
1 + tk vFp (zj )
k=1

and the minimization problem transforms into,
(
)J 
K
X
p
tk
1
.
+ zj −
wk
(5.1) Ĥ = argH min L 
vFp (zj )
n
1 + tk vFp (zj )
k=1

j=1

El Karoui (2008) has shown the weak convergence of Ĥ to H, i.e Ĥ → H.
Some examples of the convex loss function L can be,


• L∞ {ej }Jj=1 = maxj max {|Re(ej )|, |Im(ej )|}
 P

• L1 {ej }Jj=1 = Jj=1 |ej |

 P
• L2 {ej }Jj=1 = Jj=1 |ej |2
For the convex loss functions described above, the estimation of H in (5.1)
reduces to a convex optimization problem (Boyd and Vandenberghe, 2004).
Karoui also provided a translation of this problem into a linear programming
problem when L∞ loss function is used. Further details can be found in
El Karoui (2008).
5.2. Implementing Karoui’s algorithm when the number of spikes is known.
Since the generalized spikes fall outside the support of the population LSD,
Karoui’s algorithm cannot be directly applied to estimate the spikes. Furthermore, Bai and Silverstein (1998) showed that the probability of a sample
eigenvalue falling outside the support of the sample LSD will go to zero as
p increases, which implies that the sample eigenvalues corresponding to the
population generalized spikes will be measure zero points in the sample LSD.
Since the spikes behave like measure zero points (or outliers) when we are
concerned about estimating the population LSD, we can exclude the sample
eigenvalues corresponding to the population generalized spikes while calculating vFp and that will lead to a more robust estimation of H. Therefore,
we will apply Karoui’s algorithm in the following way,
(I) Suppose the population covariance matrix possesses m generalized
spikes. We exclude the top m sample eigenvalues while calculating

PCA IN HIGH-DIMENSIONAL DATA

vFp ,

13

n

vFp (z) =

1 X 1
.
n−m
di − z
i=1

(II) Apply Karoui’s algorithm to obtain Ĥ. Since it is reasonable to assume
that the true population LSD is a continuous or piecewise continuous
distribution function, kernel smoothing algorithm can be used on Ĥ
to obtain a more continuous approximation of H.
(III) The quantiles of Ĥ can be considered as the estimators of the nonspikes.
(IV) Suppose, λ̂m+1 , λ̂m+2 , . . . , λ̂p are the estimated non-spikes. Then the
ψ function is estimated by,
ψ̂(α) = α +

p
X
γα
λ̂i
.
p−m
α
−
λ̂
i
i=m+1

Due to the weak convergence Ĥ → H, ψ̂ will also converge to ψ point-wise.
Thus, all the estimates provided in Section 3 and 4 will still be consistent if
we replace ψ with ψ̂.
5.3. Estimating the number of spikes. Our application of Karoui’s algorithm to the GSP model depends on the number of spikes m, which is
usually unknown. If we have some knowledge of the underlying structure of
the data, we can use it to estimate m roughly. Suppose we know that the
data are coming from a mixture of K subpopulations, and within each subpopulation the observations are i.i.d.. Since the spikes represent the between
group differences, the number of spikes should be the same as the rank of
the between group covariance matrix which is (K − 1). However in real data,
it is often hard to accurately assess the number of such homogeneous subpopulations. In those cases we can use the following algorithm to estimate
m.
(I) Start with a reasonable finite upper bound mmax of the number of
spikes. The upper bound can be selected based on prior information
on the subpopulations, or by examining the sample eigenvalues. Set
m = mmax .
(II) Use Karoui’s algorithm to estimate the population LSD and the nonspikes. Suppose the estimated non-spikes are λ̂m+1 ≥ λ̂m+2 ≥ . . . ≥
λ̂p , and the ψ function is estimated by,
ψ̂(α) = α +

p
X
γα
λ̂i
.
p−m
i=m+1 α − λ̂i

14

R. DEY AND S. LEE

(III) Find Sψ > λm+1 using Newton-Raphson algorithm such that
p
X
γ
0
ψ̂ (Sψ ) = 1 −
p−m

i=m+1

λ̂i
Sψ − λ̂i

!2
= 0.

(IV) Since any distant spike must be larger than Sψ , and ψ̂, ψ̂ 0 are both
continuous and strictly increasing functions on (Sψ , ∞), the equation
ψ̂(λ) − dk = 0 has a root in (Sψ , ∞) if and only if ψ̂(Sψ ) − dk < 0.
Therefore, find the smallest index i∗ in 1, 2, . . . , m such that di∗ ≤
ψ̂(Sψ ). If all d1 , d2 , . . . , dm are larger than ψ̂(Sψ ) then stop and select
m as the number of distant spikes. Otherwise, set m = i∗ − 1 and
repeat step (II)-(IV).
Note that the close spikes occur so close to the support of the population
LSD that they cannot be distinguished separately from the non-spikes when
the number of spikes is unknown.
6. Simulation studies and real data example.
6.1. Simulation studies. In this section we will present simulation studies
of five different scenarios to compare the performances of the proposed GSPbased methods and the existing SP-based method proposed by Lee, Zou and
Wright (2010). For each of the first four studies, we simulated a training
dataset with n = 500 individuals and p = 5000 features. The data were
generated from three subpopulations with sample sizes 100, 150 and 250.
For each subpopulation we first selected a mean vector µi by drawing its
elements randomly with replacement from {−0.3, 0, 0.3}. Then samples in
the ith subpopulation were drawn from Np (µi , V ) where V is the AR(1)
covariance matrix with variance σ 2 and autocorrelation ρ. The (σ 2 , ρ) pairs
used for the four studies were (4, 0.8), (1, 0.7), (7.5, 0.8) and (4, 0). In study
5, we considered an ultra high-dimensional setting where n = 500 samples
and p = 50000 features were simulated from a population with two spikes
at 300 and 250, and the non-spikes generated from the AR(1) covariance
structure with (σ 2 , ρ) = (4, 0.8). The population eigenvalue plots for all the
studies are shown in Figure 2.
We also generated test datasets for each study with the same settings
as the training datasets. Then we applied our GSP-based methods and the
existing SP-based method to estimate the population spikes, the angles between the sample and population eigenvectors, the correlations between the
sample and population PC scores and the asymptotic shrinkage factors. For

PCA IN HIGH-DIMENSIONAL DATA

15

all of the studies, we used the upper bound mmax = 5 to estimate the number of distant spikes using the algorithm described in 5.3. We simulated each
study 200 times to calculate the empirical biases and standard errors of the
estimates. The results are presented in Table 1.
It is clear from Table 1 that for Study 1, 2, 3 and 5 our methods reduced
the bias in all the estimates while having similar standard errors as the
existing method. The positive empirical biases in all the SP estimates suggest
that the SP method tends to overestimate all the quantities. In Study 4, since
the underlying population satisfied the SP assumption, all methods provided
very similar and almost unbiased estimates (< 1%). The results also verify
that the λ-estimates and d-estimates are asymptotically equivalent. The
performances of the λ-estimates and the d-estimates are nearly identical in
all the simulation studies.
In Study 1, the ratio of the largest non-spike with the two spikes are 0.29
and 0.48, which are substantially larger than zero. Thus according to the
discussion in Section 4 the SP model does not closely approximate the GSP
model. The results support this assertion as the SP model-based estimates
are highly biased whereas the estimates based on our methods have very
little empirical bias. On the other hand, in Study 2 the largest non-spike is
very small compared to the smallest spike (ratio 0.08). Thus the estimates
based on the SP model closely approximate the estimates based on the GSP
model, and we find very little empirical bias (< 1%) in all of the SP modelbased estimates. In Study 3, even though there were two spikes present, only
the largest population eigenvalue was a distant spike. So we presented only
the estimates corresponding to the largest population eigenvalue. Since the
ratio of the largest non-spike and the largest spike is substantially larger
than zero (0.53) in this study, we observe very high empirical bias in the SP
model-based estimates. However, our methods provided negligible empirical
biases even in the presence of a close spike. We also presented the estimated
number of distant spikes in each of the simulation studies in Supplement A.
Note that in some cases our algorithm over-estimates the number of distant
spikes. However, as the over-estimation is finite, the estimates of the distant
spikes still remain consistent.
Study 5 shows the performances of our methods when p >> n. In this
setting, p/n = 100 is comparable to n = 500. Therefore this asymptotic
setting falls into the category of ultra-high dimensional data where p/n →
∞. Although our asymptotic results were derived based on the assumption
that p/n → γ < ∞, the results from study 5 show that even for this ultrahigh dimensional setting our methods provided almost unbiased estimates.

16

R. DEY AND S. LEE

6.2. Application on Hapmap III data. For this demonstration we used
genetic data from the Hapmap Phase III project (http://hapmap.ncbi.
nlm.nih.gov/). Our sample consisted of unrelated individuals sampled from
two different populations: a) Utah residents with Northern and Western European ancestry (CEU) and b) Toscans in Italy (TSI). We only included
genomic markers that are on chromosome 1-22, have less than 5% missing
values, and those with minor allele frequency more than 0.05. We also excluded 2 samples (both from CEU) with outlier PC scores (more than six
standard deviations away from the mean PC score corresponding to at least
one distant spike). We then mean-centered and variance-standardized the
data for each marker. The final sample consisted of 198 individuals (110
from CEU and 88 from TSI). Total number of markers selected across chromosome 1-22 was 1389511.
To evaluate the performance of the proposed methods with different p,
we performed PCA on each chromosome separately. The number of markers
varied from 19331 (chromosome 21) to 116582 (chromosome 2). The distribution of the number of markers across different chromosomes are presented
in Supplement A. We first estimated the number of distant spikes using the
algorithm described in Section 5.3. We found no distant spike in chromosome
22 and only one distant spike in chromosome 2. Then we applied our GSPbased methods and the existing SP-based method (Lee, Zou and Wright,
2010) to estimate the asymptotic shrinkage factors corresponding to the
distant spikes. Figure 3a, 3b compares the estimated asymptotic shrinkage
factors for the first two PCs across different chromosomes. The plots show
that for all the chromosomes, λ-GSP and d-GSP methods provided almost
equal estimates while the SP estimates are larger than both the GSP estimates. This suggests that the SP method would over-estimate the shrinkage
factors when the population eigenvalues deviate from the assumption that
the non-spiked eigenvalues are the same.
To investigate whether the proposed shrinkage-bias adjustment can improve the prediction accuracy, we performed a leave-one-out cross-validation.
In each iteration we removed one individual (test sample) and performed
PCA on the remaining individuals (training samples) to predict the PC score
of the test sample. For each predicted PC score, we adjusted the shrinkagebias using the GSP-based and SP-based shrinkage factor estimates. One
important issue with this cross-validation is that the exclusion of one individual can substantially change the PC-coordinates, in which the PC score
plots from the training sample-based and complete sample-based PCA can
be substantially different. In order to circumvent this problem, in each iteration we first rescaled the PC scores based on their corresponding sample

PCA IN HIGH-DIMENSIONAL DATA

17

eigenvalues to make the PCs comparable. In addition, we obtained the mean
squared difference of the training sample PC1-2 scores with and without the
exclusion of the test sample (for chromosome 2, only PC1 is used), and
excluded the test sample from the prediction error estimation if the mean
squared difference was above a threshold . We used four different values
0.5, 1, 5 and 10 for the threshold parameter , and for each value of  we
calculated the mean squared errors (MSE) of the unadjusted and adjusted
(using the GSP and SP-based methods) predicted PC scores of the test samples. The sample sizes of the test samples that were finally included in the
prediction error estimation for different values of  are shown in Supplement
A. Figure 3c shows the estimated MSEs for  = 1. It is clear that both the
λ-GSP and d-GSP methods have much smaller MSEs than the SP method.
As expected, the unadjusted predicted PC scores have substantially larger
MSE than all the proposed adjustments. The plots are very similar for the
other values of , and they can be found in Supplement A.
Figure 4 illustrates the shrinkage-bias adjustment for the PC1 and PC2
scores of an individual based on the markers on chromosome 7. The plot
clearly shows that the bias-adjusted PC score based on the SP model is still
biased towards zero, whereas the bias-adjusted PC score based on the GSP
model is very close to the original sample PC score. We only showed the
d-GSP adjusted score in the plot as the d-GSP and λ-GSP adjusted scores
were almost equal.
7. Conclusions and discussion. In this paper, we investigated the
asymptotic properties of PCA under the Generalized Spiked Population
model and derived estimators of the population eigenvalues, the angles between the sample and population eigenvectors, and the correlation coefficients between the sample and population PC scores. We also proposed
methods to adjust the shrinkage bias in the predicted PC scores. Since the
proposed methods do not require the equality of the non-spiked eigenvalues,
they can be widely used in high-dimensional biomedical data analysis. We
also implemented all our algorithms in the R package hdpca.
We note that Mestre (2008b,a) proposed an asymptotic setting similar
to the generalized spiked population model but with a different assumption
on the number of spikes in which the number of spikes increases with the
dimension. Under this assumption, he provided asymptotic properties of
sample eigenvalues and eigenvectors. However, in many biomedical data,
the number of spikes is usually finite as the spikes represent the difference
between finitely many underlying subpopulations. Therefore we believe that
the generalized spiked population model is more appropriate in such cases.

18

R. DEY AND S. LEE

Since our results were derived under the assumption that p and n grow
at the same rate, one open question is whether the method can be applied
for the ultra high-dimensional data where p is greatly larger than n. In
our simulation studies, however, we showed that our method can provide
almost unbiased estimators when p/n is as large as 100. We note that, there
exist asymptotic regimes which are suitable for the ultra high-dimensional
data. For example, the high-dimension low sample size regime, in which
p → ∞ but n is finite, was discussed in Hall, Marron and Neeman (2005);
Ahn et al. (2007); Jung and Marron (2009). Lee, Zou and Wright (2014)
considered an alternative ultra high-dimensional regime, where p and n grow
to infinity with p/n → ∞. They also showed that methods developed under
the assumption that p/n is finite can be applied to the ultra-high dimensional
data. In future research, we will use this asymptotic regime to investigate
behaviors of our methods in the ultra-high dimensional regime.
In some special cases, even though the features exhibit strong local correlation, one can use the spiked population model based methods after some
suitable data manipulation. In genome-wide association studies, SNP pruning (Anderson et al., 2010) can be used to remove locally correlated SNPs
to satisfy the spiked population model. For example, Lee, Zou and Wright
(2010) reported good performance of the spiked population model-based
methods with the SNP-pruned Hapmap III dataset. This approach, however,
can lead to a considerable loss of information; the SNP-pruning in Hapmap
III data removed nearly 90% of the SNPs. Since the proposed approach does
not require this additional step, it can use most of the information present
in the data.
APPENDIX A: PROOFS
Proof of Theorem 1. The first part of the proof follows directly from
Result 1 along with the fact that on the domain of the distant spikes, the ψ
function is strictly increasing, and hence is left invertible. Since ψ 00 (α) > 0
for any α > sup ΓH , ψ 0 (α) is a strictly increasing function for α > sup ΓH .
Let Sψ > sup ΓH be a solution for ψ 0 (α) = 0. Then for any α > sup ΓH ,
ψ 0 (α) > 0 if and only if α > Sψ . Therefore the interval (Sψ , ∞) is the domain
of the distance spikes, and ψ is a strictly increasing function on this interval.
The second part follows from Lemma 2.
Proof of Theorem 2. The proof closely follows the proof of Theorem 2
in Mestre (2008b). However, contrary to Mestre (2008b), we do not assume
that the population LSD contains the generalized spikes. Thus, some of the
derivation steps and results are substantially different from Mestre (2008b).

19

PCA IN HIGH-DIMENSIONAL DATA

We start the derivation by first noting that the quadratic forms η̂k can be
expressed as contour integrals of a special class of Stieltjes transforms of the
sample covariance matrix. Let us define,
m̂p (z) := sT1 (Sp − zIp )−1 s2 =

p
X
sT1 ej eTj s2
j=1

dj − z

;

∀z ∈ C+

where s1 and s2 are non-random vectors with uniformly bounded norms.
Girko (1996) and Mestre (2006) showed that under the assumption that the
population LSD contains the generalized spikes,
a.s.

|m̂p (z) − mp (z)| −−→ 0;

(A.1)
where

mp (z) =

sT1

−1

[w(z)Σp − zIp ]

∀z ∈ C+

p
X
sT1 Ej EjT s2
s2 =
.
w(z)λj − z
j=1

w(z) is defined as w(z) = 1 − γ − γzbF (z) where bF (z) =
RThe function
(τ − z)−1 dF (τ ) is the Stieltjes transform of the sample LSD. It is easy to
check, by the same arguments provided in Mestre (2006), that the result still
holds when the generalized spikes are considered lying outside the support
of the population LSD. The functions m̂p , mp and bF can be extended to
C− = {z ∈ C : Im(z) < 0} by defining m̂p (z) = m̂∗p (z ∗ ), mp (z) = m∗p (z ∗ )
and bF (z) = b∗F (z ∗ ) for z ∈ C− where z ∗ is the complex conjugate of z. With
a.s.
this definition, |m̂p (z) − mp (z)| −−→ 0 even when z ∈ C− . Now η̂k can be
expressed as an integral of m̂p ,
I
1
η̂k =
m̂p (z)dz,
2πi ∂ R̂−
y (k)
√
where i = −1, y > 0 and ∂ R̂−
y (k) is the negatively (clockwise) oriented
boundary of the rectangle R̂y (k) = {z ∈ C : â1 ≤ Re(z) ≤ â2 , |Im(z)| ≤ y}.
â1 and â2 can be arbitrarily chosen provided that R̂y (k) contains only the
sample eigenvalue dk and no other sample eigenvalue. Then the following
lemma gives the almost sure limit of η̂k .
Lemma 1.
1
2πi

I

1
m̂p (z)dz −
−
2πi
∂ R̂y (k)

I
∂R−
y (k)

a.s.

mp (z)dz −−→ 0,

where y > 0 and ∂R−
y (k) is the negatively (clockwise) oriented boundary
of the rectangle Ry (k) = {z ∈ C : a1 ≤ Re(z) ≤ a2 , |Im(z)| ≤ y}. a1 and a2

20

R. DEY AND S. LEE

can be arbitrarily chosen so that ψ(λk ) ∈ [a1 , a2 ] and [a1 , a2 ] ⊂ ψ (Sψ , ∞)
where Sψ > sup ΓH , ψ 0 (Sψ ) = 0. ψ (Sψ , ∞) denotes the image of the interval
(Sψ , ∞) under ψ.
Lemma 1 implies
(A.2)

ηˆk −

p
X
j=1

1
2πi

I
∂R−
y (k)

dz
w(z)λj − z

!
a.s.

sT1 Ej EjT s2 −−→ 0.

Now we need to evaluate the integral in (A.2) in order to get the almost sure
limit of the random variable η̂k . First, we extend the ψ function to Ry (k) as
follows,


Z
λdH(λ)
ψ(z) := z 1 + γ
, ∀z ∈ Ry (k).
z−λ
According to Marc̆enko and Pastur (1967), for all z ∈ C+ , bF (z) = b is the
unique solution to the following equation
Z
dH(λ)
(A.3)
b=
λ(1 − γ − γzb) − z
in the set {b ∈ C : γb − (1 − γ)/z ∈ C+ }. It is easy to see that bF also satisfies (A.3) when z ∈ C− . Now we formally define the fF function introduced
in (2.1),
(A.4)

fF (z) :=

z
z
=
,
w(z)
1 − γ − γzbF (z)

∀z ∈ C \ R.

Then bF can be expressed in terms of fF as,
bF (z) =

(1 − γ)fF (z) − z
.
γzfF (z)

By replacing b with [(1 − γ)f − z] /γzf in (A.3),


Z
λdH(λ)
(A.5)
f 1+γ
= z.
f −λ
It is easy to see that bF is a solution to (A.3) if and only if fF is a solution
to (A.5). Therefore, for all z ∈ C+ (similarly for z ∈ C− ), fF (z) = f is the
unique solution to (A.5) on C+ (respectively, C− ). This implies ψ (fF (z)) = z
for all z ∈ Ry (k) \ [a1 , a2 ].
Now we focus on the case when z ∈ R\{0}. According to Silverstein and
Choi (1995), we can extend bF to R\{0} by defining bF (z) = limy→0+ bF (z + iy)

PCA IN HIGH-DIMENSIONAL DATA

21

for any z ∈ R\{0}. The definition of fF can also be extended in a similar
fashion. In Lemma 2 we have shown that fF is the inverse function of ψ
on (Sψ , ∞), and there exists Mf > sup ΓF for which ψ (Sψ , ∞) = (Mf , ∞).
Thus, [a1 , a2 ] ⊂ ψ (Sψ , ∞) implies ψ (fF (z)) = z for all z ∈ Ry (k). Furthermore, the function ψ is continuous and differentiable on Ry (k), and the
derivative is given by,
Z 

0

ψ (z) = 1 − γ

λ
z−λ

2
dH(λ).

Then the integral in (A.2) can be expressed in terms of ψ and fF as follows,
I
I
1
1
dz
dz
=
z
−
−
2πi ∂Ry (k) w(z)λj − z
2πi ∂Ry (k) fF (z) λj − z
(A.6)
I
1
fF (z)
1
=
.
dz.
−
2πi ∂Ry (k) λj − fF (z) ψ (fF (z))
The integrand in the final expression is holomorphic on R−
y (k) when j 6= k
and possesses a simple pole ψ(λk ) when j = k. Therefore, when j 6= k the
integral in (A.6) is zero. When j = k, Applying the residue theorem on the
final integral,
I
1
ψ(λk ) − z
fF (z)
dz
= lim
.
−
2πi ∂Ry (k) w(z)λk − z z→ψ(λk ) λk − fF (z) ψ (fF (z))
ψ(λk ) − ψ (fF (z)) fF (z)
.
λk − fF (z)
ψ (fF (z))
z→ψ(λk )
λk ψ 0 (λk )
=
.
ψ(λk )
=

lim

This implies,
1
2πi

I
∂R−
y (k)

dz
=
w(z)λj − z

(λ

0 (λ )
k
ψ(λk )

kψ

0

j=k
j 6= k

and the proof is complete.
Proof of Lemma 1. First, we show that â1 , â2 , a1 , a2 can be chosen
satisfying â1 → a1 and â2 → a2 . This is possible due to the fact that
a.s.
dk −−→ ψ(λk ) and ψ(λk ) ⊂ ψ (Sψ , ∞) = (Mf , ∞) where Mf > sup ΓF .
Therefore, we can choose a neighborhood [a1 , a2 ] around ψ(λk ) so that
[a1 , a2 ] ⊂ (Mf , ∞). Moreover, as Mf is bounded away from the support

22

R. DEY AND S. LEE
a.s.

of the sample LSD F and dk −−→ ψ(λk ), we can select a neighborhood
[â1 , â2 ] around dk which does not contain any other eigenvalue for which
â1 → a1 , â2 → a2 . Then,
(A.7)

I
1
m̂p (z)dz −
mp (z)dz
2πi ∂R−
∂ R̂−
y (k)
y (k)
)
(
1
|m̂p (z)| (|â1 − a1 )| + |â2 − a2 )|)
≤
sup
2π z∈∂ R̂−
+
(k)
(k)∪∂R
y
y
I
1
+
|m̂p (z)dz − mp (z)| |dz|.
2π ∂R−
y (k)

1
2πi

I

From Cauchy-Schwartz inequality, we can obtain the following upper bound
for m̂p ;
ks1 kks2 k
,
|m̂p (z)| ≤
d z, ΓFp

where d z, ΓFp = inf y∈ΓFp |z − y|. Since Fp → F point-wise and [a1 , a2 ] is

bounded away from ΓF , d z, ΓFp is bounded away from zero with probability one for large enough p and n. Therefore |m̂p (z)| is finite for z ∈ Ry (k) with
probability one for large enough p and n. Moreover, since [â1 , â2 ] → [a1 , a2 ],
the interval [â1 , â2 ] will eventually be bounded away from ΓF . Thus, eventually the upper bound for |m̂p (z)| will also be finite for z ∈ R̂y (k). Therefore, the first term on the right hand side of (A.7) will go to zero as
â1 → a1 , â2 → a2 .
Now, as m̂p (z) and mp (z) are holomorphic functions on the compact set
∂R−
y (k),
sup |m̂p (z) − mp (z)| < ∞.
z∈∂R−
y (k)
a.s.

Also from (A.1), |m̂p (z) − mp (z)| −−→ 0 point-wise for all z ∈ C \ R. Therefore, by dominated convergence theorem the second term on the right hand
side of (A.7) also converges to zero almost surely.
We can show the asymptotic equivalence of the limits derived in Theorem
2 and Result 2 as a direct application of the following lemma.
Lemma 2. Suppose assumptions (A)–(C) hold. If λk is a distant spike
with multiplicity one, and dk is the corresponding sample eigenvalue, then
p

fF (dk ) →
− λk ;

dk gF (dk ) p 0
→
− ψ (λk ).
fF (dk )

PCA IN HIGH-DIMENSIONAL DATA

23

Proof. We have already established in the proof of Theorem 2 that for
all z ∈ C+ (similarly for z ∈ C− ), fF (z) = f is the unique solution to (A.5)
−
on C+ (respectively,
z is restricted to C \ R, using (A.4) and the
R C ). When
−1
fact that bF (z) = (τ − z) dF (τ ) we can write,
fF (z) =

1+γ

z
R τ dF (τ ) .
z−τ

Now suppose z = x ∈ R \ {0}. Then both equations (A.3) and (A.5) will
have multiple roots (both real and complex valued depending on x and H).
If we look at (A.5) closely, we can see for real valued x it can be represented
as ψ(f (x)) = x, where the ψ function is as defined in (2.1). As we have seen
in the proof of Theorem 1, ψ is strictly increasing in the interval (Sψ , ∞)
where Sψ > sup ΓH and ψ 0 (Sψ ) = 0. Therefore, any real-valued solution f
of ψ(f (x)) = x in (Sψ , ∞) has to be the inverse of ψ, which is unique due
to the strict monotonicity of ψ on (Sψ , ∞). Now suppose ΓF is the support
of the sample LSD F . We will show that there exists Mf > sup ΓF such
that for any x > Mf , the function fF is real-valued and it is a solution to
(A.5) in the interval (Sψ , ∞). Thus it is also the unique such solution and
the inverse of the ψ function in (Sψ , ∞).
Let x ∈ R, x > sup ΓF and z = x + iy ∈ C+ . Now, as z ∈ C+ , fF (z) is the
unique solution to (A.5) in C+ . Therefore, if we express fF (z) as u(z)+iv(z),
then v(z) > 0. Also, the imaginary part of (A.5) can be written as


Z
λ2
v(z) 1 − γ
= y.
{u(z) − λ}2 + v(z)2
Both v(z) and y being positive implies that
Z
λ2
(A.8)
1−γ
> 0.
{u(z) − λ}2 + v(z)2
Due to the continuity of fF on the set {z ∈ C+ : z = x + iy, x > sup ΓF },
fF (x) = lim

y→0+

x + iy
x
R τ dF (τ ) =
R τ dF (τ ) ,
1 + γ x+iy−τ
1+γ
x−τ

which is real-valued. Thus u(z) → fF (x) and v(z) → 0 as y → 0+ . Therefore
as y → 0+ , the inequality (A.8) becomes
Z
λ2
1−γ
> 0,
{fF (x) − λ}2

24

R. DEY AND S. LEE

which implies ψ 0 (fF (x)) > 0.
We can see that fF (x) attains zero at sup ΓF and it is strictly and unboundedly increasing for x > sup ΓF . This ensures the existence of a threshold MF > sup ΓF such that the function fF maps the interval (MF , ∞) to
(Sψ , ∞). Therefore, fF and ψ are both strictly increasing, continuous and
bijective mappings between the intervals (MF , ∞) and (Sψ , ∞). Since fF (z)
is the unique solution to (A.5) in C+ when z ∈ C+ , fF is also a solution to
(A.5) in (Sψ , ∞) when x > MF due to the continuity of the left hand side
of (A.5) on the set {f ∈ C+ : f = u + iv, u > Sψ }, which further implies
that fF is the inverse function of ψ on (Sψ , ∞).
The first part of this lemma is proved as a corollary to Result 1 as ψ −1 =
fF on the domain of distant spikes, i.e. (Sψ , ∞). For the second part we first
need to derive the expression of fF0 , and then derive the expression of ψ 0 in
terms of fF and F .
Z
f (x)
τ dF (τ )
0
fF (x) =
[1 + γvF (x)] ; vF (x) =
.
x
(x − τ )2
For a distant spike λk , using the expression of fF0 we get,
λk ψ 0 (λk )
λk
1
=
=
R
0
ψ(λk )
ψ(λk )fF (ψ(λk ))
1 + γfF (ψ(λk ))

τ dF (τ )
[ψ(λk )−τ ]2

= gF (ψ(λk )) .

p

As ψ(λk ) > Mf , gF is continuous at ψ(λ). Since dk →
− ψ(λk ),
p

gF (dk ) →
− gF (ψ(λk )) =

λk ψ 0 (λk )
;
ψ(λk )

dk gF (dk ) p 0
→
− ψ (λk ).
fF (dk )

Proof of Theorem 3.

2
T
1
1
2
TX X
hPk , pk i = 2
hXEk , Xek i =
Ek
ek
n λk dk
λk dk
n
"
! #2
p
X
1
dk
=
EkT
di ei eTi ek =
hek , Ek i2 .
λk dk
λk
2

i=1

Using the limits derived in Theorem 2 and Result 1,
dk
p
hek , Ek i2 − ψ 0 (λk ) →
− 0.
λk

PCA IN HIGH-DIMENSIONAL DATA

25

Using Lemma 2,
dk gF (dk ) p
dk
hek , Ek i2 −
→
− 0.
λk
fF (dk )

 
Proof of Theorem 4. We show that the denominator E p2kj converges

to ψ(λk ) and the numerator E qk2 converges to λ2k /ψ(λk ). The proof will
p
be complete using the fact that dk →
− ψ(λk ).
The denominator,
!
!
n
n
X
X

1
1
E p2kj = E
p2ki = E
(xTi ek )2
n
n
i=1
i=1
"
! #


p
T
X
T
T
TX X
ek = E ek
di ei ei ek = E(dk ) → ψ(λk ).
= E ek
n
i=1

The numerator,





E qk2 = E (xTnew ek )2 = E E(xTnew ek )2 |ek




= E V ar(xTnew ek )|ek = E eTk Σp ek .
Now, usingRthe notations in the proof of Theorem 2 and Lemma 2, we
have bF (z) = (τ − z)−1 dF (τ ) as the Stieltjes transform of the sample LSD
and the function fF defined as fF (z) = z [1 − γ − γzbF (z)]−1 . Therefore,
bF (z) =

(1 − γ)fF (z) − z
.
γzfF (z)

The functions bF and fF can be extended to the real axis by defining the
extensions as shown in the proof of Lemma 2. Thus, for the sample eigenvalue
dk corresponding to the distant spike λk we have
bF (dk ) =

(1 − γ)fF (dk ) − dk
.
γdk fF (dk )

According to Theorem 4 in Ledoit and Péché (2010), the limit of eTk Σp ek
is given by dk [1 − γ − γdk bF (dk )]−2 . Replacing the expression of bF (dk ) in
this limit, we get
f 2 (dk ) p
eTk Σp ek − F
→
− 0.
dk

26

R. DEY AND S. LEE
p

Using Result 1 and Lemma 2 we have fF2 (dk )/dk →
− λ2k /ψ(λk ). Therefore,
the limit of the numerator is given by,


E(qk2 ) = E eTk Σp ek →

λ2k
.
ψ(λk )

SUPPLEMENTARY MATERIAL
Supplement A: Supplementary tables and figures for “Asymptotic properties of Principal Component Analysis and shrinkagebias adjustment under the Generalized Spiked Population model”
(SuppA.pdf). Some additional tables and figures relevant to this paper are
provided in this supplementary material.
REFERENCES
Ahn, J., Marron, J. S., Muller, K. M. and Chi, Y.-Y. (2007). The High-Dimension,
Low-Sample-Size Geometric Representation Holds under Mild Conditions. Biometrika
94 760-766.
Anderson, T. W. (1963). Asymptotic Theory for Principal Component Analysis. Ann.
Math. Statist. 34 122–148.
Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed.
Wiley Series in Probability and Statistics, Hoboken, NJ.
Anderson, C. A., Pettersson, F. H., Clarke, G. M., Cardon, L. R., Morris, A. P.
and Zondervan, K. T. (2010). Data quality control in genetic case-control association
studies. Nature protocols 5 1564-1573.
Bai, Z. (1999). Methodologies in spectral analysis of large-dimensional random matrices,
a review. Statist. Sinica 9 611–677.
Bai, Z. D. and Silverstein, J. W. (1998). No eigenvalues outside the support of the
limiting spectral distribution of large-dimensional sample covariance matrices. Ann.
Probab. 26 316–345.
Bai, Z. D. and Silverstein, J. W. (1999). Exact Separation of Eigenvalues of Large
Dimensional Sample Covariance Matrices. Ann. Probab. 27 1536–1555.
Bai, Z. and Silverstein, J. W. (2006). Spectral Analysis of Large Dimensional Random
Matrices. Science Press, Beijing.
Bai, Z. and Yao, J. (2012). On sample eigenvalues in a generalized spiked population
model. Journal of Multivariate Analysis 106 167 - 177.
Baik, J. and Silverstein, J. W. (2006). Eigenvalues of large sample covariance matrices
of spiked population models. Journal of Multivariate Analysis 97 1382 - 1408.
Barrett, J. C., Fry, B., Maller, J. and Daly, M. J. (2005). Haploview: analysis and
visualization of LD and haplotype maps. Bioinformatics 21 263-265.
Boyd, S. P. and Vandenberghe, L. (2004). Convex Optimization. Cambridge University
Press, Cambridge.
Ding, X. (2015). Convergence of Sample Eigenvectors of Spiked Population Model. Communications in Statistics - Theory and Methods 44 3825-3840.
El Karoui, N. (2008). Spectrum estimation for large dimensional covariance matrices
using random matrix theory. Ann. Statist. 36 2757–2790.

PCA IN HIGH-DIMENSIONAL DATA

27

Girko, V. L. (1996). Strong Law for the eigenvalues and eigenvectors of empirical covariance matrices. Random Operators and Stochastic Equations 4 176–204.
Hall, P., Marron, J. S. and Neeman, A. (2005). Geometric representation of high
dimension, low sample size data. Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 67 427–444.
Johnstone, I. M. (2001). On the distribution of the largest eigenvalue in principal components analysis. Ann. Statist. 29 295–327.
Johnstone, I. M. and Lu, A. Y. (2009). On Consistency and Sparsity for Principal Components Analysis in High Dimensions. Journal of the American Statistical Association
104 682-693. PMID: 20617121.
Jung, S. and Marron, J. S. (2009). PCA consistency in high dimension, low sample size
context. Ann. Statist. 37 4104–4130.
Ledoit, O. and Péché, S. (2010). Eigenvectors of some large sample covariance matrix
ensembles. Probability Theory and Related Fields 151 233–264.
Lee, S., Zou, F. and Wright, F. A. (2010). Convergence and prediction of principal
component scores in high-dimensional settings. Ann. Statist. 38 3605–3629.
Lee, S., Zou, F. and Wright, F. A. (2014). Convergence of sample eigenvalues, eigenvectors, and principal component scores for ultra-high dimensional data. Biometrika.
Marc̆enko, V. A. and Pastur, L. A. (1967). Distribution of eigenvalues for some sets
of random matrices. Mathematics of the USSR-Sbornik 1 457.
Mestre, X. (2006). On the asymptotic behavior of quadratic forms of the resolvent of
certain covariance-type matrices. Technical Report, CTTC/RC/2006-001, Centre Tecnologic de Telecomunicacions de Catalunya.
Mestre, X. (2008a). On the Asymptotic Behavior of the Sample Estimates of Eigenvalues
and Eigenvectors of Covariance Matrices. IEEE Transactions on Signal Processing 56
5353-5368.
Mestre, X. (2008b). Improved Estimation of Eigenvalues and Eigenvectors of Covariance
Matrices Using Their Sample Estimates. IEEE Transactions on Information Theory 54
5113-5129.
Paul, D. (2007). Asymptotics of sample eigenstruture for a large dimensional spiked
covariance model. Statist. Sinica 17 1617–1642.
Price, A. L., Patterson, N. J., Plenge, R. M., Weinblatt, M. E., Shadick, N. A.
and Reich, D. (2006). Principal components analysis corrects for stratification in
genome-wide association studies. Nature Genetics 38 904–909.
Silverstein, J. W. (1995). Strong Convergence of the Empirical Distribution of Eigenvalues of Large Dimensional Random Matrices. Journal of Multivariate Analysis 55 331
- 339.
Silverstein, J. W. and Choi, S. I. (1995). Analysis of the Limiting Spectral Distribution
of Large Dimensional Random Matrices. Journal of Multivariate Analysis 54 295 - 309.
Stein, E. M. and Shakarchi, R. (2013). Complex Analysis. Princeton University Press,
Princeton, NJ.
Storey, J. D., Xiao, W., Leek, J. T., Tompkins, R. G. and Davis, R. W. (2005).
Significance analysis of time course microarray experiments. Proceedings of the National
Academy of Sciences of the United States of America 102 12837-12842.

28

R. DEY AND S. LEE

APPENDIX B: TABLES AND FIGURES
Table 1
Simulation results for GSP-based and SP-based methods for estimating the population
eigenvalues, cosine of the angles between sample and population eigenvectors, correlations
between sample and population PC scores, and the asymptotic shrinkage factors. Each
cell has empirical bias (%) with coefficients of variations (%) in parenthesis.
Settings
No.

Method

1

SP

n = 500
p = 5000
σ2 = 4
ρ = 0.8

λ-GSP
d-GSP

2

n = 500
p = 5000
σ2 = 1
ρ = 0.7

SP
λ-GSP
d-GSP

3

n = 500
p = 5000
σ 2 = 7.5
ρ = 0.8

SP
λ-GSP
d-GSP

4

n = 500
p = 5000
σ2 = 4
ρ=0

SP
λ-GSP
d-GSP

5

n = 500
p = 50000
σ2 = 4
ρ = 0.8

SP
λ-GSP
d-GSP

Eigenvalue
1
5.27
(2.37)
0.43
(2.67)
0.47
(2.67)
0.10
(0.90)
-0.04
(0.90)
-0.004
(0.90)
25.68
(2.54)
2.92
(5.7)
2.45
(5.74)
0.05
(1.58)
0.03
(1.58)
0.16
(1.58)
8.35
(5.75)
0.46
(6.80)
1.35
(6.70)

2
18.27
(3.11)
0.95
(5.27)
0.69
(5.45)
0.46
(1.27)
0.04
(1.28)
0.10
(1.28)
-0.26
(2.35)
-0.35
(2.35)
-0.12
(2.35)
8.78
(5.01)
-3.38
(6.59)
-2.16
(6.45)

Angle
1
6.52
(0.32)
0.53
(0.77)
0.47
(0.77)
0.16
(0.04)
0.01
(0.04)
0.03
(0.04)
64.06
(0.52)
12.62
(11.90)
12.25
(10.52)
0.06
(0.23)
0.02
(0.24)
0.10
(0.23)
13.41
(1.71)
4.25
(2.99)
4.96
(2.90)

2
34.07
(0.60)
3.28
(6.26)
2.48
(6.70)
0.44
(0.08)
0.004
(0.10)
0.03
(0.10)
-0.06
(0.53)
-0.18
(0.54)
-0.03
(0.53)
19.20
(1.64)
3.77
(3.76)
4.84
(3.62)

Correlation
1
2
3.83
23.33
(0.03)
(0.08)
0.33
2.79
(0.31)
(4.69)
0.24
2.11
(0.31)
(5.07)
0.08
0.24
(0.001) (0.003)
0.01
0.01
(0.03)
(0.01)
0.004
0.01
(0.03)
(0.01)
46.50
(0.07)
10.95
(10.13)
10.87
(8.58)
0.03
0.05
(0.02)
(0.08)
0.01
-0.02
(0.02)
(0.09)
0.01
0.02
(0.02)
(0.09)
8.98
13.53
(0.09)
(0.12)
4.05
4.87
(0.85)
(1.46)
4.29
5.29
(0.81)
(1.39)

Shrinkage
1
2
5.32
17.88
(0.60) (1.06)
0.47
0.58
(0.92) (3.31)
0.51
0.31
(0.92) (3.51)
0.18
0.39
(0.08) (0.16)
0.03
-0.03
(0.08) (0.18)
0.07
0.03
(0.08) (0.18)
26.41
(0.90)
3.47
(4.20)
3.00
(4.24)
0.07
-0.22
(0.43) (0.90)
0.04
-0.31
(0.43) (0.91)
0.18
-0.08
(0.42) (0.90)
8.51
9.46
(3.23) (3.04)
0.58
-2.82
(4.28) (4.64)
1.48
-1.59
(4.19) (4.50)

29

PCA IN HIGH-DIMENSIONAL DATA

Fig 1: Eigenvalue structures for SP and GSP models
(a) Spiked population

(b) Generalized spiked population

Fig 2: Eigenvalue structures in simulation studies
(a) Study 1

(b) Study 2

(d) Study 4

(c) Study 3

(e) Study 5

30

R. DEY AND S. LEE

Fig 3: Estimated shrinkage factors for (a) PC1 and (b) PC2 across chromosomes 1-21 based on three different methods. (c) Comparison of the mean
squared errors (MSE) of the unadjusted and adjusted PC scores based on
the d-GSP and SP methods with the adjusted PC scores based on the λGSP method. The ratios of the MSEs are presented for chromosome 1-21
using the threshold  = 1. The Y-Axis is presented in a logarithmic scale.
(a)

(b)

(c)

PCA IN HIGH-DIMENSIONAL DATA

31

Fig 4: PC1 vs PC2 plot of the Hapmap III CEU and TSI samples based
on chromosome 7. The predicted PC scores for the illustrative individual,
and its bias-adjusted PC scores are also presented. Since the d-GSP and the
λ-GSP adjusted scores are nearly the same, the λ-GSP adjusted scores are
not presented.

32

R. DEY AND S. LEE
Department of Biostatistics
University of Michigan School of Public Health
1415 Washington Heights
Ann Arbor, MI 48109-2029, USA
E-mail: deyrnk@umich.edu
leeshawn@umich.edu

