arXiv:1708.05026v1 [] 16 Aug 2017

Adjusting systematic bias in high
dimensional principal component scores
Sungkyu Jung
Department of Statistics, University of Pittsburgh, Pittsburgh, PA 15260, USA
e-mail: sungkyu@pitt.edu
Abstract: Principal component analysis continues to be a powerful tool
in dimension reduction of high dimensional data. We assume a variancediverging model and use the high-dimension, low-sample-size asymptotics
to show that even though the principal component directions are not consistent, the sample and prediction principal component scores can be useful
in revealing the population structure. We further show that these scores are
biased, and the bias is asymptotically decomposed into rotation and scaling
parts. We propose methods of bias-adjustment that are shown to be consistent and work well in the finite but high dimensional situations with small
sample sizes. The potential advantage of bias-adjustment is demonstrated
in a classification setting.
Keywords and phrases: proportional bias, HDLSS, high-dimension, lowsample-size, jackknife, principal component analysis, pervasive factor.

1. Introduction
Principal component analysis is a workhorse method of multivariate analysis,
and has been used in a variety of fields for dimension reduction, visualization
and as exploratory analysis. The standard estimates of principal components,
obtained by either the eigendecomposition of the sample covariance matrix or
the singular value decomposition of the data matrix, are now well-known to be
inconsistent when the number of variables, or the dimension d, is much larger
than the sample size n (Paul, 2007; Johnstone & Lu, 2009; Jung & Marron,
2009). These observations were paralleled with a vast amount of proposals on,
e.g., sparse principal component estimations (cf. most notably, Zou et al., 2006),
which performs better in some models with high dimensions.
However, the standard estimates of principal components continues to be
useful, partly due to available fast computations (see, e.g., Abraham & Inouye,
2014). Many of the sparse estimation methods, unfortunately, do not computationally scale well for large data with hundreds of thousands of variables.
Moreover, the standard estimation has shown to be useful in some application
areas such as imaging, genomics and big-data analysis (Fan et al., 2014). In
these areas, the sample principal component scores (the projection scores of the
data points onto the principal component directions) are often used in the next
stage of analysis, such as regression and classification. The predicted principal
component scores of new observations can be obtained as well, and can be used
as the input to fitted models for prediction.
1

S. Jung/Bias in principal component scores

60

Sample
True

40

40

20

20

PC 1

PC 1

60

0

-20

-40

-40

-60

-60
-20

-10

0

PC 2

10

20

30

Prediction
True

0

-20

-30

2

-30

-20

-10

0

10

20

30

PC 2

Fig 1. Sample and prediction principal component scores connected to their true values. This
toy data set of size pd, nq “ p10000, 50q is generated from the spike model with m “ 2 spikes,
with polynomially-decreasing eigenvalues with β “ 0.3; see Section 4.2 for details.

In this paper, we revisit the standard estimates of principal components in
ultra-high dimensions and reveal that while the component directions and variances are inconsistent, the sample and prediction scores are useful for moderately large sample size. For low sample sizes, the scores are biased. We quantify
the bias, decompose it into two systematic parts, and propose to estimate biasadjustment factors.
As a visual example of the systematic bias, a toy data set with 2 distinguishable principal components is simulated and plotted in Fig. 1. Each observation
in the data set consists of d “ 10, 000 variables. The first two sample principal
component directions are estimated from n “ 50 observations, and are used
to obtain the sample and prediction scores (the latter are computed from 20
new observations). The true principal scores for each of observations are also
plotted and connected to their empirical counterparts. This example visually
reveals that the sample scores are systematically biased, that is, uniformly rotated and stretched. What is more surprising is that the prediction scores are
also uniformly rotated, by the same angle as the sample scores, and uniformly
shrunk.
On the other hand, the third component scores from this example appear to
be quite arbitrary; see Fig. 2. (The estimates for component 3 in this example
is only as good as random guess.) Moreover, unlike the first two components
plotted in Fig. 1, the sample scores of the third component are grossly inflated,
while the prediction scores are much smaller than the sample scores.

S. Jung/Bias in principal component scores

Sample
True

40

40

20

20

0

0

-20

-20

-40

-40

-60

-60

-40

-20

PC 3

0

20

Prediction
True

60

PC 1

PC 1

60

3

-40

-20

0

20

PC 3

Fig 2. Sample and prediction principal component scores connected to their true values.
Models and data are the same as in Fig. 1.

In Section 2, we provide theoretical justification of the phenomenon observed
in Figs. 1 and 2, and asymptotically quantify the two parts of the systematic
bias. We assume m-component models with diverging variances, and use the
high-dimension, low-sample-size asymptotic scenario (i.e. d Ñ 8 while n is
fixed). The models and asymptotics are used in giving the contrasting results
of the sample and prediction scores. The correlation coefficients between the
sample (or prediction) and true scores turn out to be close to 1, for large signals
and large sample sizes, indicating the situations where the principal component
scores are most useful.
Since the bias is asymptotically quantified, the natural next step is to adjust
the bias by estimating the bias-adjustment factor. In Section 3, we propose a
simple, yet consistent, estimator and several variants of estimators based on
the idea of Jackknife. Adjusting these biases improves the performance of prediction modeling, and we demonstrate its potential by an example involving
classification. Results from numerical studies are summarized in Section 4.
There are several related work on the principal component scores in high
dimensions (Lee et al., 2010; Fan et al., 2013; Lee et al., 2014; Sundberg &
Feldmann, 2016; Shen et al., 2016; Hellton & Thoresen, 2017; Jung et al., 2017).
This paper is built upon these previous work. Connections to those are discussed
in Section 5. A survey of high-dimension, low-sample-size asymptotics can be
found in Aoshima et al. (2017).

S. Jung/Bias in principal component scores

4

2. Asymptotic behavior of principal component scores
2.1. Model and assumptions
Let X “ rX1 , . . . , Xn s be a d ˆ n data matrix, where each Xi is mutually
independent and has zero mean and covariance matrix Σd . Population principal components are obtained by the eigendecomposition of Σd “ U ΛU 1 , where
Λ “ diagpλ1 , . . . , λd q is the diagonal matrix of principal component variances
and U “ ru1 , . . . , ud s consists of principal component directions. For a fixed m,
we assume an m-component model, where the first m component variances are
distinguishably larger than the rest. Specifically, the larger variances increasing
at the same rate as the dimension d, i.e. λi — d, which was previously noted as
the “boundary situation” (Jung et al., 2012). This diverging-variance condition
seems to be more realistic than the other simpler cases λi " d and λi ! d (Hellton & Thoresen, 2017; Shen et al., 2016), and is satisfied for high-dimensional
models used in factor analysis (Fan et al., 2013; Li et al., 2017; Sundberg &
Feldmann, 2016).
We assume that the population principal component variances satisfy the
following:
2
(A1) λi “ σi2 d, i “ 1, . . . , m, σ12 ě ¨ ¨ ¨ ě σm
.
řd
2
(A2) limdÑ8 i“m`1 λi {d :“ τ P p0, 8q.
(A3) There exist B ă 8 such that for all i ą m, lim supdÑ8 λi ă B.

The conditions (A2) and (A3) are used to allow λi for i ą m increase as d
increases. All of our results hold when the condition (A3) is relaxed to, e.g., allow
the situation that λi — dα , α ă 1{2. Such a generalization is straightforward,
but invites nonintuitive technicality (see, e.g., Jung et al., 2012, 2017). By
decomposing each independent observation into the first m components and the
remaining term, we write
Xj “

m
ÿ
i“1

1{2

λi ui zij `

d
ÿ

1{2

λi ui zij ,

pj “ 1, . . . , nq,

(1)

i“m`1

where zij is the normalized principal component score.
(A4) For each j “ 1, 2, . . ., pz1j , z2j , . . .q is a sequence of independent random
variables such that for any i, Epzij q “ 0, Varpzij q “ 1, and that the fourth
moment of zij is uniformly bounded.
2.2. Sample and prediction principal component scores
Suppose we have a data matrix X “ rX1 , . . . , Xn s and a vector X˚ , independently drawn from the same population with principal component directions ui .
The principal component analysis is performed for data X and is used to predict
the principal component scores of X˚ .

S. Jung/Bias in principal component scores

5

We define the ith true principal component scores of X as the vector of n
projection scores:
wiT “ uTi X “ pwi1 , . . . , win q, pi “ 1, . . . , dq,
(2)
?
where wij “ uTi Xj “ λi zij . The last equality is given by the decomposition
of Xj in ?
(1). Likewise, the true ith principal component score of X˚ is wi˚ “
uTi X˚ “ λi zi˚ .
The classical estimators of the pair of the ith principal component direction and variance are pûi , λ̂i q, obtained by either the eigendecomposition of the
sample covariance matrix Sd “ n´1 X X T ,
Sd “

n
ÿ

λ̂i ûi ûTi ,

i“1

or by the singular value decomposition of the data matrix,
X “

n b
? ÿ
n
λ̂i ûi v̂iT ,

(3)

i“1

where v̂i is the right singular vector of X . By replacing ui in (2) with its estimator
ûi , we define the ith sample principal component scores of X as
ŵiT “ ûTi X “ pŵi1 , . . . , ŵin q,

pi “ 1, . . . , nq.

(4)

The sample principal component scores
are in fact weighted right singular veca
tors of X ; comparing to (3), ŵi “ nλ̂i v̂iT .
For an independent observation X˚ , the definition (4) gives
ŵi˚ “ ûTi X˚ ,
which is called the ith prediction principal component score for X˚ .
2.3. Main results
Denote W1 “ pσi zij qi,j “ pd´1{2 wij qi,j “ d´1{2 ru1 , . . . , um sT X for the m ˆ n
matrix of the scaled true scores for the first m principal components. The ith
row of W1 is d´1{2 wiT . Similarly, the scaled sample scores for the first m principal
x1 “ d´1{2 rû1 , . . . , ûm sT X .
components are denoted by W
xT “
For a new observation X˚ , write W˚ “ d´1{2 pw1˚ , . . . , wm˚ qT and W
˚
d´1{2 pŵ1˚ , . . . , ŵm˚ qT for the scaled true scores, and prediction scores, respectively, of the fist m principal components.
Denote that W “ W1 W1T for the scaled mˆm sample covariance matrix of the
first m scores. Let tλi pSq, vi pSqu denote the ith largest eigenvalue-eigenvector
pair of a non-negative definite matrix S and vij pSq denote the jth loading of
the vector vi pSq.

S. Jung/Bias in principal component scores

6

Theorem 1. Assume the m-component model under Conditions (A1)–(A4) and
let n ą m ě 0 be fixed and d Ñ 8. Then, the first m sample and prediction
scores are systematically biased:
x1 “ SRT W1 ` Op pd´1{4 q
W
x˚ “ S
W

´1

T

´1{2

R W˚ ` Op pd

(5)
q

where R “ rv1 pWq, . . . , vm pWqs, S “ diagpρ1 , . . . , ρm q, and ρk “
Moreover, for k ą m,
ŵkj “ Op pd1{2 q,
ŵk˚ “ Op p1q.

j “ 1, . . . , n,

(6)
a
1 ` τ 2 {λk pWq.
(7)
(8)

Our main results show that the first m sample and predictions scores are
comparable to the true scores. The asymptotic relation tells that for large d, the
x1 converge to the true scores in W1 , uniformly rotated
first m sample scores in W
and scaled for all data points. It is thus valid to use the first m sample principal
scores for exploration of important data structure, to reduce the dimension of
the data space from d to m in the high-dimension, low-sample-size context.
Theorem 1 explains and quantifies the two parts of the bias, exemplified
in Fig. 1. In particular, the same rotational bias applies to both sample and
prediction scores. The scaling bias factors ρk in the matrix S are all greater
than 1. Thus, while the sample scores are all stretched, the prediction scores
are all shrunk. The second part of the theorem shows that the magnitude of
inflation for the sample scores of the “noise” component (see, e.g., component
3 scores in Fig. 2) is of order d1{2 . On the other hand, the prediction scores of
the noise component do not diverge.
Remark 1. Suppose m “ 1 in Theorem 1. Then the sample and prediction scores
are simply proportionally-biased in the limit: ŵ1j {w1j Ñ ρ1 and ŵ1˚ {w1˚ Ñ ρ´1
1
in probability as d Ñ 8. There is no rotational bias.
The proof of Theorem 1 relies on the asymptotic behavior of the principal
component direction and variance, which is now well-understood. For reference
we restate it here.
Lemma 2. [Theorem S2.1, Jung et al. (2017)] Assume the conditions of Theorem 1. (i) the sample principal component variances converge in probability as
d Ñ 8;
"
λi pWq ` τ 2 ` Op pd´1{2 q, i “ 1, . . . , m;
´1
d nλ̂i “
τ 2 ` Op pd´1{2 q,
i “ m ` 1, . . . , n.
(ii) The inner product between sample and population PC directions converges
in probability as d Ñ 8;
" ´1
ρi vij pWq ` Op pd´1{2 q, i, j “ 1, . . . , m;
T
ûi uj “
Op pd´1{2 q,
otherwise.

S. Jung/Bias in principal component scores

7

This result is abridged later in Section 2.4 for discussion. To handle prediction
scores, we need in addition the following observation, summarized in Lemma 3.
For each k “ 1, . . . , m, the kth projection score ŵk˚ is decomposed into
m
ÿ

ŵk˚ “ ûTk X˚ “

wi˚ ûTk ui ` k˚ ,

(9)

i“1

řd
where k˚ “ i“m`1 wi˚ ûTk ui . In the next lemma, we show that the “error
term,” k˚ , is stochastically bounded.
Lemma 3. Assume the m-factor model with (A1)–(A4) and let n ą m ě 0 be
fixed. For k “ 1, . . . , n, Epk˚ |W1 q “ 0, and
2
lim Varpk˚ | W1 q “ υO
{pλk pWq ` τ 2 q,

lim

dÑ8

dÑ8
n
ÿ

for k ď m;

1
2
Varpk˚ | W1 q “ υO
{τ 2 ,
n ´ m k“m`1

2
where υO
“ limdÑ8 d´1

(10)
(11)

řd

λ2i . As d Ñ 8, k˚ “ Op p1q.
?
Proof of Lemma 3. Fix k “ 1, . . . , n. Let Yi “ λi zi˚ pki , where pki “ ûTk ui .
řd
Then k˚ “ i“m`1 Yi . Since zi˚ and pki are independent, for each i ą m,
EpYi | W1 q “ 0 and
Varp

d
ÿ

i“m`1

Yi | W1 q “ Ep

i“m`1

d
ÿ

2 2
λi zi˚
pki | W1 q “

i“m`1

λi Epp2ki | W1 q,

i“m`1

where we use the fact that Epzi˚ q “ 0,
For k ď m, if the following claim,
Epp2ki | W1 q “ d´1

d
ÿ

2
Epzi˚
q

“ 1.

λi
` Opd´3{2 q,
pλk pWq ` τ 2 q

(12)

is true for any i ą m, then it is easy to check (10).
To show (12), we first post-multiply v̂i to (3) to obtain ûi “ pnλ̂i q´1{2 X v̂i .
´1{2
By writing ziT “ λi wiT “ pzi1 , . . . , zin q from (2), we have
pki “ uTi ûk
“ pnλ̂k q´1{2 uTi X v̂k
1{2

“ pnλ̂k q´1{2 λi ziT v̂k .
Thus,
p2ki “ d´1

λi
nd´1 λ̂k

pziT v̂k q2

λi
pz T v̂k q2
λk pWq ` τ 2 ` Op pd´1{2 q i
λi
“ d´1
pz T v̂k q2 ` Op pd´3{2 q.
λk pWq ` τ 2 i

“ d´1

(13)

S. Jung/Bias in principal component scores

8

In (13), we used Lemma 2(i) and that p1 ` xq´1 “ 1 ` Opxq, and the fact that
2
2
2
|ziT v̂k |2 ď }zi }2 }v̂k }2 “ }zi }2 “ Op p1q.
Write pziT v̂k q2 “ rziT vk pW1T W1 q ` ziT pv̂k ´ vk pW1T W1 qqs2 . Note that W1T W1 is
an n ˆ n matrix, and is different from the m ˆ m matrix W “ W1 W1T . It can
be shown that the right singular vector v̂k converges to vk pW1T W1 q (see, e.g.,
Lemma S1.1 of Jung et al., 2017): For k “ 1, . . . , m,
v̂k “ vk pW1T W1 q ` Op pd´1{2 q.

(14)

Thus we get |ziT pv̂k ´ vk pW1T W1 qq| ď }zi }2 }v̂k ´ vk pW1T W1 qq}2 “ Op pd´1{2 q.
Therefore,
EppziT v̂k q2 | W1 q “ EppziT vk pW1T W1 qq2 | W1 q ` Opd´1{2 q
n
ÿ
2
2
“
qvk`
pW1T W1 q ` Opd´1{2 q
Epzi`
`“1

“ 1 ` Opd´1{2 q.

(15)

Combing (13) and (15), we get (12) for k ď m as desired.
To show (11), note that W “ W1 W1T is of rank m. For k ą m, with λk pWq “
0, (13) holds. Thus,
n
n
d
ÿ
ÿ
ÿ
1
1
Varpk˚ | W1 q “
λi Epp2ki | W1 q
n ´ m k“m`1
n ´ m k“m`1 i“m`1

“

(16)

d
n
ÿ
ÿ
1
λ2i {τ 2
EppziT v̂k q2 | W1 q.
dpn ´ mq i“m`1
k“m`1

To simplify the expression EppziT v̂k q2 | W1 q, one should not naively try (15).
This is because that (15) does not apply for k ą m due to the non-unique kth
eigenvector vk pW1T W1 q of the rank-m matrix W1T W1 . Instead, from
n
ÿ

pziT v̂k q2 “ ziT zi ´

k“m`1

m
ÿ

pziT v̂k q2 ,

k“1

and (15) for k ď m, we get
n
ÿ

EppziT v̂k q2 | W1 q “ n ´ m ` Opd´1{2 q.

(17)

k“m`1

Taking the limit d Ñ 8 to (16), combined with (17), leads to (11).
The last statement, k˚ “ Op p1q, easily follows from the fact limdÑ8 Varpk˚ q ď
2
υO
{τ 2 pn ´ mq ă 8, which is obtained by (10) and (11).
We are now ready to show Theorem 1. Note that the results on the sample
scores,
(5) and (7), can be easily shown, using the decomposition d´1{2 ŵk “
a
d´1 nλ̂k v̂k , together with Lemma 2(i) and (14). We show (6) and (8).

S. Jung/Bias in principal component scores

9

Proof of Theorem 1. Proof of (6). Recall the decomposition (9). Using the notation pki “ ûTk ui , we write ŵk˚ “ ppk1 , . . . , pkm qpw1˚ , . . . , wm˚ qT ` k˚ . Putting
all parts together, we have
¨
˛
p11 ¨ ¨ ¨ p1m
.
.. ‹ W ` ˜ ,
..
x˚ “ d´1{2 pŵ1˚ , . . . , ŵm˚ qT “ ˚
W
˝ ..
k˚
.
. ‚ ˚
pm1

¨¨¨

pmm

where ˜k˚ “ d´1{2 p1˚ , . . . , m˚ qT . By Lemma 3, as d Ñ 8, ˜k˚ “ Op pd´1{2 q.
´1{2
Since pki “ ρ´1
q, by Lemma 2(ii), we have
k vki pWq ` Op pd
x T “ S ´1 RT W T ` Op pd´1{2 q.
W
˚
˚
Proof of (8). Using the decomposition
řm (9), and by the fact k˚ “ Op p1q,
from Lemma 3, it is enough to show i“1 wi˚ pki “ Op p1q. But, since Lemma 2
1
implies d 2 pki “ Op p1q for any pair of pk, iq such that k ą m, i ď m, we have
řm
1
1
2
2
i“1 wi˚ pki “ σi pd pk1 , . . . , d pkm qpz1˚ , . . . , zm˚ q “ Op p1q,
Next result shows that the sample and true scores (or prediction and true
scores) are highly correlated with each other. For this,
a we compute the inner product between standardized sample scores ŵk { ŵkT ŵk and true scores
a
?
wk { wkT wk . Define for a pair px, yq of n-vectors rpx, yq “ xT y{ xT x ¨ y T y,
which is an empirical correlation coefficient between x and y when the mean is
assumed to be zero.
řm 2
řm 2
Theorem 4. Let ζkj “ λk pWq{p `“1 v`j
pWqλ` pWqq and ζ̄kj “ σk2 {p `“1 v`j
pWqσ`2 q.
Under the assumptions of Theorem 1, as d Ñ 8, for k, j “ 1 . . . , m,
1{2

(i) rpŵk , wj q Ñ vkj pWqζkj in probability ;
(ii) limdÑ8 Corrpŵk˚ , wj˚ | W1 q “ vkj pWqζ̄kj 1{2 .
Remark 2. In the special case, m “ 1, both the sample and prediction scores
of the first principal component are perfectly correlated with the true scores, in
the limit. Specifically, Theorem 4 leads that |rpŵ1 , w1 q| Ñ 1 in probability and
|Corrpŵk˚ , wj˚ q| Ñ 1 as d Ñ 8.
1{2

Remark 3. The somewhat complex limiting quantity vkj pWqζkj is an artifact
of fixed sample size. To simplify the expression for the case k “ j, write
´
¯2
ÿ
λ` pWq
1
1{2
2
, ξk pWq “
v`k
pWq
.
vkk pWqζkk
“
1 ` ξk pWq
λ
k pWq
`‰k
Note that W “ W1 W1T is proportional to the sample covariance matrix of the
first m true scores, and that vkk pWq is the inner product between the kth sample
and theoretical principal component directions of the data set W1 , where the
number of variables, m, is smaller than the sample size n. Therefore, we expect
that |vkk pWq| « 1 and ξk pWq « 0 for large sample size n. Taking the additional
limits n Ñ 8, the results in Theorem 4 become more interpretable:
|rpŵk , wj q| Ñ 1pk“jq in probability, and |Corrpŵk˚ , wj˚ q| Ñ 1pk“jq ,

S. Jung/Bias in principal component scores

10

as d Ñ 8, n Ñ 8 (limits are taken progressively).
Remark 4. What is the correlation coefficient rpŵk , wk q for k ą m in the limit
d Ñ 8? In an attempt to answer this question, we note ŵk “ pnλ̂k q1{2 v̂k ,
řd
v̂k “ vk pX T X q and X T X “ i“1 wi wiT . Thus,
rpŵk , wk q “ wkT vk p

d
ÿ

a
wi wiT q{ λk ,

i“1

and it is natural to guess that the dependence of v̂k on any wi , including the
case i “ k, would diminish as d tends to infinity. In fact, d´1 X T X converges to
the rank-m matrix S0 :“ W1T W1 ` τ 2 In , (Jung et al., 2012), and wk and S0 are
independent. Thus, it is reasonable to conjecture that limdÑ8 Errpŵk , wk qs “ 0,
for k ą m. Unfortunately, in the limit d Ñ 8, the kth, k ą m, eigenvector
of d´1 X T X becomes an arbitrary choice in the left null space of W1 . Due to
this non-unique eigenvector, the inner product wkT vk pS0 q is not defined, and
consequently discussing the convergence of rpŵk , wk q is somewhat demanding.
We numerically confirm the conjecture in Section 4.1.
Proof of Theorem 4. Proof of (i). Write the singular value decomposition of the
m ˆ n matrix of scaled scores W1 as
a
a
(18)
W1 “ Rdiagp λ1 pWq, . . . , λ1 pWqqGT ,
where G “ rg1 , . . . , gm s is the n ˆ m matrix consisting of right singular vectors
of W1 . The left singular vector matrix R “ rv1 pWq, . . . , vm pWqs is exactly the
matrix R appearing in Theorem 1. Since
m a
ÿ
λ` pWqv` pWqg`T ,

W1 “

`“1

the jth row of W1 is, for j ď m,
1

d´ 2 wjT “

m a
ÿ

λ` pWqv`j pWqg`T .

`“1
´1{2
For the scaled sample
ŵa
k , k ď m, we obtain from Theorem 1 and
ascore d
x
(18) that Wa
λ1 pWqqGT ` Op pd´1{4 q and its kth row
1 “ Sdiagp λ1 pWq, . . . ,
´1{2
´1{4
2
d
ŵk “ λk pWq ` τ gk ` Op pd
q. Since g` ’s are orthonormal,
1

}d´ 2 ŵk }2 “

a
λk pWq ` τ 2 ` Op pd´1{4 q,

and
d´1 ŵkT wj “ pd´1{2 ŵk qT pd´1{2 wj q
a
a
“ λk pWq λk pWq ` τ 2 vkj pWq ` Op pd´1{4 q.

S. Jung/Bias in principal component scores

Since d´1 wjT wj “

řm
`“1

11

2
v`j
pWqλ` pWq, we have

rpŵk , wj q “

d´1 ŵkT wj
1{2
Ñ vkj pWqζkj
}d´1{2 ŵk }2 ¨ }d´1{2 wj }2

in probability, as d Ñ 8.
Proof of (ii). From Theorem 1, write
d´1{2 ŵk˚ “ ρ´1
k

m
ÿ

vk` pWqd´1{2 w`˚ ` Op pd´1{2 q,

(19)

`“1

and note that Epwk˚ q “ Epŵk˚ q “ 0. Then for k “ 1, . . . , m, we have
Varpd´1{2 wk˚ q “ d´1 Epwk˚ q2 “ σk2 Epzk˚ q2 “ σk2 ,
and, by (19),
Varpd´1{2 ŵk˚ | W1 q “ ρ´2
k

m
ÿ

2

pvk` pWqq σ`2 ` Opd´1{2 q.

`“1

The independence of w`˚ and wk˚ for k ‰ ` and (19) give
Covpd´1{2 ŵk˚ , d´1{2 wj˚ | W1 q “ Epd´1 ŵk˚ wj˚ | W1 q
2
´1{2
“ ρ´1
q,
k vkj pWqσj ` Opd

which in turn leads to
corrpŵk˚ , wj˚ | W1 q “ `

Covpd´1{2 ŵk˚ , d´1{2 wj˚ | W1 q

˘1{2
Varpd´1{2 wj˚ qVarpd´1{2 ŵk˚ | W1 q
σj
´1{2
“ vkj pWq ”
q.
ı1{2 ` Opd
řm
2 2
pv
pWqq
σ
k`
`“1
`

2.4. Inconsistency of the direction and variance estimators
The findings in the previous subsection may be summarized as that the first m
principal component scores convey about the same visual information as the true
values when displayed. (The information is further honed by bias adjustment in
Section 3.) In a practical point of view, the scores and their graph matter the
most.
On the other hand, a quite different conclusion about the standard principal
component analysis is made when the standard estimator of the principal component direction ûi is of interest. The asymptotic behavior of the direction ûi

S. Jung/Bias in principal component scores

12

as well as the variance estimator λ̂i are obtained as a special case of Lemma 2.
Under our model,
" ´1
pρi vii pWq, λi pWq ` τ 2 q, i “ 1, . . . , m;
T
´1
pûi ui , d nλ̂i q Ñ
(20)
p0, τ 2 q,
i “ m ` 1, . . . , n.
in probability as d Ñ 8 (n is fixed).
The variance estimator λ̂i , for i ď m, is asymptotically proportionally-biased.
Specifically, λ̂i {λi Ñ pλi pWq ` τ 2 q{pnσi2 q in probability as d Ñ 8. Thus by
using a classical result on the expansion of the eigenvalues of W for large n,
«
ff
m
σj2
τ2
1 ÿ
` 2 ` Opn´2 q,
Epλ̂i {λi q Ñ 1 `
n j‰i σi2 ´ σj2
σi
as d Ñ 8. Note that even when m “ 1, the bias is still of order n´1 . This
proportional bias may be empirically adjusted, using good estimates of σi2 and
τ 2 . We do not pursue it here. Note that all empirical principal component variances, for i ą m, converge to τ 2 {n, when scaled by d, and thus do not reflect
any information of the population.
The result (20) also shows that the direction estimator ûi is inconsistent and
asymptotically-biased, compared to ui . The estimator ûi is closer to ui when
ρ´1
i |vii pWq| is closer to 1. This is impossible to achieve since for finite n, both
|vii pWq| and ρ´1
are strictly less than 1. Although the “angle” between ûi and
i
ui is quantified in (20), the theorem itself is useless in adjusting the bias. This is
because that the direction to which ûi moves away from ui is not specified in the
theorem, and is conjectured to be random (i.e. uniformly distributed). See Jung
et al. (2012) and Shen et al. (2016) for related discussions on the inconsistent
principal component directions.
In short, while the bias in the principal component direction is challenging to
remove, the bias in sample and prediction scores can be quantified and removed.
3. Bias-adjusted scores
In this section, we describe and compare several choices in estimation of the biasadjustment factor ρi , of which the matrix S in Theorem 1 consists. Since both
sample and prediction scores are rotated by the same direction and amount,
specified in the matrix R, there is little practical advantage in estimating R. We
focus on adjusting the scores by estimating ρi .
Suppose that the number of effective principal components, m, is prespecified
or estimated inaadvance. Our first estimator is obtained by replacing τ 2 and
λi pWq in ρi “ 1 ` τ 2 {λi pWq with reasonable estimators. In particular, we set
řn
λ̂i n
2
τ̃ “ i“m`1
, λ̃i pWq “ d´1 nλ̂i ´ τ 2 ,
(21)
n´m d
and
ρ̃i “

b
1 ` τ̃ 2 {λ̃i pWq,

pi “ 1, . . . , mq.

(22)

S. Jung/Bias in principal component scores

Sample (Adj)
True

40

13

Prediction (Adj)
True

50
40

30
30
20
20
10

PC 1

PC 1

10

0

0
-10

-10

-20
-20
-30
-30
-40
-40
-50
-20

-10

0

10

20

-20

PC 2

-10

0

10

20

30

PC 2

Fig 3. Bias-adjusted sample and prediction scores using (23) for the toy data introduced in
Fig. 1. The estimates (22) are pρ̃1 , ρ̃2 q “ p1.385, 1.546q, very close to the theoretical values
pρ1 , ρ2 q “ p1.385, 1.557q. Both sample and prediction scores are simultaneously rotated about
16 degrees clockwise.

This simple estimator ρ̃i is in fact consistent.
Corollary 5. Suppose the assumptions of Lemma 2 are satisfied. Let d Ñ 8.
For i “ 1, . . . , m, conditional to W1 , τ̃ 2 , λ̃i pWq and ρ̃i are consistent estimators
of τ 2 , λi pWq and ρi , respectively.
Proof. Lemma 2 is used to show that τ̃ 2 and λ̃i pWq converge in probability
to τ 2 and λi pWq as d Ñ 8, respectively. By continuous mapping theorem, ρ̃i
converges in probability to ρi .
padjq

Using (22), the bias-adjusted sample and prediction scores are ŵi
“ ρ̃´1
i ŵi
padjq
and ŵi˚ “ ρ̃i ŵi˚ for i “ 1, . . . , m. The sample and prediction scores matrices
in (5) and (6) are then adjusted to, using S̃ “ diagpρ̃‘ , . . . , ρ̃m q,
x padjq “ S̃ ´1 W
x1 ,
W
1

x˚padjq “ S̃ W
x˚ .
W

(23)

An application of the above bias-adjustment procedure is exemplified in
Fig. 3. There, the magnitudes of the sample and prediction scores are welladjusted. Adjustment for the ‘rotation’ part is not needed, since both sample
and prediction scores are simultaneously rotated.
Our next proposed estimators are motivated by the well-known jackknife
bias adjustment procedures and also by the leave-one-out cross-validation. For
simplicity, assume m “ 1. The bias-adjustment factor we aim to estimate is

S. Jung/Bias in principal component scores

14

ρ1 “ p1 ` τ 2 {}ω1 }22 q1{2 , where ω1 “ d´1{2 w1 “ σ1 pz11 , . . . , z1n qT is the scaled
true scores for the first principal component.
Write, for each j “ 1, . . . , n, the jth scaled sample score as ω̂1j “ d´1{2 ûT1 Xj ,
and the jth scaled prediction score as
ω̂1pjq “ d´1{2 ûT1p´jq Xj ,
where û1p´jq is the first principal component direction, computed from Xp´jq ,
i.e., the data except the jth observation.
Then from Theorem 1, ρ1 is the asymptotic bias-adjustment factor for ω̂1 ;
ω̂1j “ ρ1 ω1j ` Op pd´1{4 q. For ω̂1pjq , again applying Theorem 1, we get ω̂1pjq “
´1{2
ρ´1
q, where ρ1p´jq “ p1 ` τ 2 {}ω1p´jq }22 q1{2 is the bias-adjustment
1p´jq ω1j `Op pd
factor computed from Xp´jq , using ω1p´jq “ σ1 pz11 , . . . , z1,j´1 , z1,j`1 , . . . z1n qT .
To simplify terms, Taylor expansion is used to expand ρ1p´jq as a function of
2
ω1j
{n, resulting in
˜
ρ1p´jq “

τ 2 {n
1`
2 {n
}ω1 }22 {n ´ ω1j

¸1{2
“ ρ1 `

2
1 }ω1 }22 {n ω1j
1
` Op p 2 q. (24)
2
2ρ1 τ
n
n

Using the approximation
ρ1 ρ1p´jq « ρ21 `

2
}ω1 }22 ω1j
2τ 2 n2

given by (24), we write the ratio of the sample and prediction scores to cancel
out the unknown true score ω1j as follows:
ˆ

ŵ1j
ŵ1pjq

˙1{2

ˆ
“

ω̂1j
ω̂1pjq

˙1{2
« ρ1 .

Based on above heuristic, we define the following estimators of the bias-adjustment
factors:
˙1{2
n ˆ
1 ÿ ŵij
p1q
ρ̂i “
,
(25)
n j“1 ŵipjq
˜ řn
¸1{2
j“1 ŵij
p2q
ρ̂i “ řn
,
(26)
j“1 ŵipjq
˜ řn
¸1{4
2
j“1 ŵij
p3q
ρ̂i “ řn
.
(27)
2
j“1 ŵipjq
In implementing the above estimators, we used absolute values of the sample
and predicted scores.
p1q
p2q
p3q
The estimators ρ̂i , ρ̂i , and ρ̂i tend to overestimate ρ for small sample size
n, as expected from (24). In our numerical experiments, these three estimators
perform similarly.

S. Jung/Bias in principal component scores

15

4. Numerical studies
4.1. Simulations to confirm the asymptotic bias and near-perfect
correlations
In this section, we compare the theoretical asymptotic quantities derived in
Section 2.3 with their finite-dimensional empirical counterparts.
First, the theoretical values of the scaling bias ρi and the rotation matrix
R in Theorem 1 are compared with their empirical counterparts. The empirical counterparts of the two matrices R, S are defined as the minimizer of the
Procrustes problem
›2
›
›
x T S ´1 R0 ›› ,
(28)
min ›W1 ´ W
1

0

F

with the constraint that S0 is a diagonal matrix with positive entries and R0 is an
orthogonal matrix. The solutions are denoted by Sq “ diagpρ̌1 pW1 q, . . . , ρ̌m pW1 qq
q For simplicity, we consider the m “ 2 case, and parameterize R by the
and R.
q by θ̌R “ cos´1 pR
q1,1 q. We compare
rotation angle, θR “ cos´1 pR1,1 q, and R
θR with θ̌R and ρi pW1 q with ρ̌i pW1 q, from a 2-component model with pn, dq “
p50, 5000q (precisely, the spike model with m “ 2 and β “ 0.3 in Section 4.2).
Note that both the theoretical values and the best-fitted values depend on the
true scores W1 . To capture the natural variation given by W1 , the experiment
is repeated for 100 times. The results, summarized in the top row of Fig. 4,
confirm that the asymptotic statements in Theorem 1 approximately hold for
q are very close
finite dimensions. In particular, the rotation matrices R and R
to each other. The Procrustes-fitted, or “best”, ρ̌i tends to be larger than the
asymptotic, or theoretical, ρi , especially for i “ 2 (shown as
in Fig. 4) and for
larger values of ρ2 . This is not unexpected. Larger values of ρ2 are from smaller
λ2 pWq. Take an extreme case where λ2 pWq “ 0, then by (7) in Theorem 1,
the sample scores are of magnitude d1{2 compared to the true scores. Thus,
as λ2 pWq decreases to 0, the Procrustes scaler ρ̌2 empirically interpolates the
finite-scaling case (5) to the diverging case (7) of Theorem 1.
Second, we compare the limit of correlation coefficients in Theorem 4 with
finite-dimensional empirical correlations, rpŵk , wk q, for k “ 1, 2. For the correlation coefficient of the prediction scores, we use the sample correlation coefficient
between pŵk˚ , wk˚ q, as an estimate of Corrpŵk˚ , wk˚ | W1 q. The simulated results are shown in the bottom row of Fig. 4. The empirical correlation coefficients
tend to be smaller than the theoretical counterparts. Note that the approximation of rpŵk , wk q by its limit is of better quality if nσk2 “ Epλk pWqq is larger.
Moreover, the correlation coefficient tends to be larger for larger nσk2 , which
represents “signal strength” of the kth principal component.
Third, from the same simulations, it can be checked that the kth, where
k ą m, sample scores are diverging, while the prediction scores are stable, as
indicated in (7) and (8). For this, we choose k “ 3 and for each experiment, comy ŵ3 q, the sample variance of the sample scores, and an approximation
pute Varp
of Varpŵ3˚ q. The results are shown in Table 1. As expected, the sample scores

S. Jung/Bias in principal component scores

Angles (in degree) in R

30

16

Bias-adjustment factors in S
2.5

25

Best

Best

20
15

2

10
1.5

5
0
0

5

10

15

20

25

30

1.5

2

Theory

1

Correlation coefficient (sample)

1

Correlation coefficient (prediction)

0.95

Empirical

0.95

Empirical

2.5

Theory

0.9

0.85

0.9

0.85

0.8

0.8
0.8

0.85

0.9

Theory

0.95

1

0.8

0.85

0.9

0.95

1

Theory

Fig 4. (Top row) Theoretical rotation angles θR and bias-adjustment factors ρ1 (ˆ), ρ2 ( ),
compared with the best-fitting Procrustes counterparts (θ̌R , ρ̌i pW1 q). (Bottom row) Empirical
correlations compared with their limits in Theorem 4.

are grossly inflated, while the prediction scores are stable. Finally, the conjecture in Remark 4 is also empirically checked; Table 1 also shows that for large
d, the sample (or prediction) and true scores for the kth, k ą m, component are
nearly uncorrelated.
4.2. Numerical performance of the bias-adjustment factor
estimation
We now test our estimators of the bias-adjustment factor ρi , using the following
data-generating models with m “ 2.
The first one is called a spike model. We sample from the d-dimensional
zero-mean normal distribution where the first two largest eigenvalues of the
covariance matrix are λi “ σi2 d, for i “ 1, 2, where pσ12 , σ22 q “ p0.02, 0.01q.
´β
The rest
řd of eigenvalues are slowly-decreasing. In particular, λi “ τ i , where
τ “ r i“3 i´β {pd ´ 2qs´1 . We set β “ 0.3 or 0.5. This spike model has more
than two unique principal components, for each fixed dimension, but in the limit
d Ñ 8, only the first two principal components are useful.

S. Jung/Bias in principal component scores

17

Sample scores
Prediction scores
120.7(4.4)
1.38(0.2)
-0.0024(0.2)
-0.004(0.15)
Table 1
The case k ą m. The 3rd sample and prediction scores (unadjusted). Shown are the mean
(standard deviation) of the variances and correlation coefficients to true scores, from 100
repetitions. The true variance is λ3 “ Varpw3˚ q « 6.5.
Variance
Corr. Coef.

The second model is a mixture model. Let µg (g “ 1, 2, 3) be d-dimensional
vectors, the elements of which are randomly drawn from t´a, 0, au with replacement for a given a ą 0, then assumed as fixed quantities. Given µg ’s we
sample from the mixture model X | G “ g „ N pµg , Id q, P pG “ gq “ pg ą 0,
ř3
g“1 pg “ 1. We set pp1 , p2 , p3 q “ p0.5, 0.3, 0.2q. It can be shown that the mean
EpXq is non-zero, and covariance matrix satisfies the assumption of 2-component
model in (A1)–(A4).
For various cases of high-dimension, low-sample-size situations, ranging d “
5, 000 to 20, 000 and n “ 50 to 100, samples from each of these models are
obtained. For each of samples, the theoretical quantity ρi “ ρi pW1 q and the
best-fitted Procrustes scaler ρ̌i “ ρ̌i pW1 q are computed. These quantities depend
on the m ˆ n random matrix W1 . The mean and standard deviation of ρi (from
100 repetitions) are shown in the first column of Table 2. As expected, the
theoretical value ρi depends on the sample size n; large sample size decreases
the bias, Epρi q, and also decreases the variance Varpρi q.
The mean of the best-fitted scaler ρ̌i (i “ 1, 2) is displayed in the second
column of the table. While they are quite close to the theoretical counterpart,
ρ̌i s are significantly larger for the mixture model, whose signal-to-noise ratio is
smaller than the spike model, and for the not-so-large dimension d “ 5, 000. This
is not unexpected, since the theoretical values are also based on the increasingdimension asymptotic arguments.
We further compute the proposed estimators of ρi , given in (22), (25)–(27).
We also compute the estimator derived from Lee et al. (2010), which is the
square-root of the reciprocal of the shrinkage factor, obtained by numerical
iterations, denoted dˆν in Lee et al. (2010). (The relation of Lee et al. (2010) to
our work is further discussed in Section 5.) Table 2 summarizes the results from
our simulation study. All of the methods considered provide accurate estimates
of the theoretical quantity ρi . We omit the numerical results from the estimators
(26) and (27), as their performances are very close to that from (25).
4.3. Bias-adjustment improves classification
Our last simulation study is an application of the bias-adjustment procedure to
classification. Our training and testing data, each with sample size 100, are sampled from the mixture model with three groups, as described in Section 4.2. As
frequently used in practice, dimension reduction by the standard principal component analysis is performed first, then our classifier, a support vector machine
(SVM, Cristianini & Shawe-Taylor, 2000), is trained on the sample principal

S. Jung/Bias in principal component scores

Spike model
β “ 0.3

Spike model
β “ 0.5

Mixture model
a “ 0.15

d
5000
10000
10000
20000
5000
10000
10000
20000
5000
10000
10000
20000

n
50
50
100
100
50
50
100
100
50
50
100
100

d
5000
10000
10000
20000
5000
10000
10000
20000
5000
10000
10000
20000

n
50
50
100
100
50
50
100
100
50
50
100
100

Theory
1.41 (0.07)
1.42 (0.06)
1.23 (0.03)
1.23 (0.02)
1.42 (0.08)
1.43 (0.07)
1.22 (0.02)
1.23 (0.02)
2.06 (0.06)
2.09 (0.06)
1.63 (0.02)
1.64 (0.02)

Best
1.42
1.43
1.23
1.23
1.45
1.45
1.23
1.23
2.22
2.17
1.67
1.66

ρ1
Asymp.
1.40
1.42
1.23
1.23
1.41
1.43
1.22
1.23
1.92
1.98
1.61
1.62
ρ2
Asymp.
1.75
1.77
1.43
1.43
1.72
1.76
1.43
1.42
2.20
2.35
1.90
1.93

18

Jackknife
1.43
1.44
1.24
1.24
1.45
1.46
1.23
1.24
2.14
2.14
1.65
1.66

LZW
1.41
1.42
1.23
1.23
1.40
1.42
1.21
1.22
2.00
2.02
1.63
1.63

Theory
Best
Jackknife
LZW
1.79 (0.11)
1.86
1.78
1.79
Spike model
1.79 (0.11)
1.82
1.77
1.79
β “ 0.3
1.43 (0.06)
1.44
1.42
1.43
1.43 (0.05)
1.44
1.42
1.43
1.79 (0.11)
1.99
1.81
1.71
Spike model
1.80 (0.11)
1.88
1.79
1.74
β “ 0.5
1.44 (0.05)
1.47
1.44
1.41
1.42 (0.05)
1.44
1.41
1.40
2.62 (0.21)
5.44
2.68
2.46
Mixture model
2.68 (0.19)
3.20
2.68
2.50
a “ 0.15
2.00 (0.09)
2.13
2.00
1.99
1.99 (0.10)
2.05
1.97
1.97
Table 2
Simulation results from 100 repetitions. “Theory” is mean (standard deviation) of ρi ;
p1q
“Best” is ρ̌i (28); “Asymp.” is ρ̃i (22); “Jackknife” is ρ̂i (25); “LZW” is from Lee et al.
(2010). Averages are shown for the latter four columns. The standard errors of the
quantities in estimation of ρi are at most 0.04.

component scores. In this simulation, we fix m “ 2 and d “ 5000. We compare the training and testing missclassification error rates (estimated by 100
repetitions) of SVMs trained (and tested) either on the unadjusted sample and
x1 and W
x‹ , or on the bias-adjusted sample and prediction
prediction scores, W
padjq
padjq
x
x
scores, W1
and W‹
in (23). The estimated error rates are shown in Table 3. It is clear that using the bias-adjusted scores, proposed in (23), greatly
improves the performance of classification.
To better understand the huge improvement of classification performances,
we plot the sample and prediction scores that are inputs of the classifier. In
Fig. 5, the classifier is estimated from the the sample scores (symbol ) and is
used to classify future observations, i.e. the prediction scores (symbol ˆ). Due
to the scaling bias, the unadjusted sample and prediction scores are of different
scales (shown in the left panel), and classification is bound to fail. On the other
hand, the proposed bias-adjustment, shown in the right panel, works well for
this data set, leading to a better classification performance.

S. Jung/Bias in principal component scores

19

Unadjusted scores
Bias-adjusted scores
0.04(0.02)
0.07(0.03)
21.4(1.33)
1.98(0.23)
Table 3
Means (standard errors) of Missclassification error rates (in percent).
Training Error
Testing Error

Unadjusted scores

Bias-adjusted scores
10

PC 1 Scores

PC 1 Scores

10
5
0
-5
-10

5
0
-5
-10

-10

0

10

PC 2 Scores

20

-10

0

10

20

PC 2 Scores

Fig 5. Bias-adjusted scores from the mixture models greatly improve the classification performance. Different colors correspond to different groups. Symbol
represents the sample scores
(unadjusted in the left, adjusted in the right); symbol ˆ represents the prediction scores.

5. Discussion
The standard principal component analysis is shown to be useful in the dimension reduction of data from the m-component models with diverging variances.
In particular, in the high-dimension, low-sample-size asymptotic scenario we
reveal that the sample and prediction scores have systematic bias that can be
consistently adjusted. We propose several estimators of the scaling bias, while
there is no compelling reason to adjust rotational bias. The amount of bias is
large when the sample size is small and when the variance of accumulated noise
is large compared to the variances of the first m components.
This work is built upon several previous findings on the principal component scores in high dimensions. The decomposition of the bias into rotation and scaling parts is also found in Hellton & Thoresen (2017) and Jung
et al. (2017). While our current work is a continuation of the latter, Hellton
& Thoresen (2017) focused on the standardized sample scores ŵi { }ŵi }. Using
the standardized scores complicates the interpretation since the resulting scalpHTq
ing bias, say ρi
, can be any positive value. For example, if m “ 1, then
a
pHTq
1{2
“ n ρ1 { }ŵ1 } “
ρ1
n{λ1 pWq « 1{σ1 , and the the sample scores are
stretched if σ1 ă 1, and shrunk if σ1 ą 1, for any value of error variance τ 2 . In
comparison, the unstandardized sample scores ŵi are always stretched by ρi ą 1.
The identification of the systematic bias is the key in revealing the bias of prediction scores, which is not available in the previous work (Hellton & Thoresen,
2017; Jung et al., 2017).
Shen et al. (2016) used a stronger model, with eigenvalues diverging faster

S. Jung/Bias in principal component scores

20

than the dimension, i.e., λ1 {d Ñ 8, to show that there is an asymptotic scaling
bias for standardized sample scores. They, however, did not consider either the
multiple component model or the prediction scores. Moreover, their result, when
applied to unstandardized scores, is less useful as, following Shen et al. (2016),
the scores are asymptotically unbiased. For example, under their assumptions,
for a 1-component model with λ1 " d, |ŵ11 {w11 | Ñ 1 as d Ñ 8 for any fixed
sample size n. This is an artifact of assuming fast-diverging variance; our scaling
factor ρi also approaches 1 when λi pWq « nλ1 {d increases.
Lee et al. (2010) discussed adjusting bias in the prediction of principal components, based on random matrix theory and the asymptotic scenario of d{n Ñ γ P
p0, 8q, n Ñ 8. They showed that the prediction scores tend to be smaller than
the sample scores, and the ratio of the shrinkage is asymptotically sdpŵi1 q{sdpŵi˚ q «
pLZWq
pLZWq
i ´1
. This “shrinkage factor” ρi
corresponds to the squared
ρi
“ λiλ`γ´1
´2
reciprocal of our scaling bias, ρi . Specifically, these two quantities, for finite
d and n, d " n, are in fact close to each other; replacing λi by σi2 d and γ by
pLZWq
´2
i pWq
i ´n{d
« nλ
d{n results in ρi
« λnλ
by order
nλi `τ 2 , which approximates ρi
i `1´n{d
n´1{2 when τ “ 1. Our work can be thought of as an extension of Lee et al.
(2010) from the asymptotic regime d — n to the high-dimension, low-sample-size
situations (see also Lee et al. (2014)). Finally, we note that in the asymptotic
scenario of Lee et al. (2010) and Lee et al. (2014) there is no rotational bias. This
is because in their limit the sample size is infinite. We show that the rotational
bias is universal to both sample and prediction scores and is of order n´1{2 .
Principal component analysis is often thought of as a special case of factor
model or its estimation method. Our covariance model with diverging variances
frequently appears in recent investigations of high-dimensional factor models
(e.g., Fan et al., 2013; Li et al., 2017; Sundberg & Feldmann, 2016). In particular, Sundberg & Feldmann (2016) also investigated the high-dimension, lowsample-size asymptotic scenario for factor analysis, and reached a similar conclusion to this work. Our work echoes the message that “(principal component or
factor) scores are useful in high dimensions” and further provides interpretable
decomposition of bias and methods of bias-adjustment.
References
Abraham, G. & Inouye, M. (2014). Fast principal component analysis of
large-scale genome-wide data. PloS one 9, e93766.
Aoshima, M., Shen, D., Shen, H., Yata, K., hui Zhou4, Y. & Marron,
J. S. (2017). A survey of high dimension low sample size asymptotics. Aust.
N. Z. J. Stat. to appear.
Cristianini, N. & Shawe-Taylor, J. (2000). An Introduction to Support
Vector Machines. Cambridge University Press.
Fan, J., Han, F. & Liu, H. (2014). Challenges of big data analysis. Natl. Sci.
Rev. 1, 293–314.
Fan, J., Liao, Y. & Mincheva, M. (2013). Large covariance estimation by

S. Jung/Bias in principal component scores

21

thresholding principal orthogonal complements. J. R. Stat. Soc. B 75, 603–
680.
Hellton, K. H. & Thoresen, M. (2017). When and why are principal component scores a good tool for visualizing high-dimensional data? Scand. J.
Stat , to appear.
Johnstone, I. M. & Lu, A. Y. (2009). On Consistency and Sparsity for
Principal Components Analysis in High Dimensions. J. Am. Stat. Assoc.
104, 682–693.
Jung, S., Ahn, J. & Lee, M. H. (2017). On the number of principal components in high dimensions. submitted to Biometrika .
Jung, S. & Marron, J. S. (2009). PCA consistency in high dimension, low
sample size context. Ann. Stat. 37, 4104–4130.
Jung, S., Sen, A. & Marron, J. (2012). Boundary behavior in High Dimension, Low Sample Size asymptotics of PCA. J. Multivar. Anal. 109, 190–203.
Lee, S., Zou, F. & Wright, F. A. (2010). Convergence and prediction of
principal component scores in high-dimensional settings. Ann. Stat. 38, 3605.
Lee, S., Zou, F. & Wright, F. A. (2014). Convergence of sample eigenvalues,
eigenvectors, and principal component scores for ultra-high dimensional data.
Biometrika 101, 484.
Li, Q., Cheng, G., Fan, J. & Wang, Y. (2017). Embracing the blessing of
dimensionality in factor models. J. Am. Stat. Assoc. To appear.
Paul, D. (2007). Asymptotics of sample eigenstructure for a large dimensional
spiked covariance model. Stat. Sin. 17, 1617–1642.
Shen, D., Shen, H., Zhu, H. & Marron, J. (2016). The statistics and
mathematics of high dimension low sample size asymptotics. Stat. Sin. 26,
1747.
Sundberg, R. & Feldmann, U. (2016). Exploratory factor analysisparameter
estimation and scores prediction with high-dimensional data. J. Multivar.
Anal. 148, 49–59.
Zou, H., Hastie, T. & Tibshirani, R. (2006). Sparse Principal Component
Analysis. J. Comp. Graph. Stat. 15, 265–286.

