Statistica Sinica Preprint No: SS-2022-0319
Title Automatic Sparse PCA for High-Dimensional Data
Manuscript ID SS-2022-0319
URL http://www.stat.sinica.edu.tw/statistica/
DOI 10.5705/ss.202022.0319
Complete List of Authors Kazuyoshi Yata and
Makoto Aoshima
Corresponding Authors Kazuyoshi Yata
E-mails yata@math.tsukuba.ac.jp
Notice: Accepted version subject to English editing.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)Statistica Sinica

AUTOMATIC SPARSE PCA
FOR HIGH-DIMENSIONAL DATA
Kazuyoshi Yata1 and Makoto Aoshima1
1

University of Tsukuba

Abstract: Sparse principal component analysis (SPCA) methods have proven toeﬃciently analyze high-dimensional data. Among them, threshold-based SPCA(TSPCA) is computationally more cost-eﬀective than regularized SPCA, basedon L1 penalties. We herein present an investigation of the eﬃcacy of TSPCA for
high-dimensional data settings and illustrate that, for a suitable threshold value,
TSPCA achieves satisfactory performance for high-dimensional data. Thus, theperformance of the TSPCA depends heavily on the selected threshold value. To
this end, we propose a novel thresholding estimator to obtain the principal component (PC) directions using a customized noise-reduction methodology. The
proposed technique is consistent under mild conditions, unaﬀected by thresholdvalues, and therefore yields more accurate results quickly at a lower computational cost. Furthermore, we explore the shrinkage PC directions and their application in clustering high-dimensional data. Finally, we evaluate the performanceof the estimated shrinkage PC directions in actual data analyses.

Key words and phrases: Clustering, Large p small n, PCA consistency, Shrinkage

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

PC directions, Thresholding.

1.

Introduction

High-dimensional, low-sample-size (HDLSS) data scenarios exist in many
areas of modern science including genomics, medical imaging, text recognition, and ﬁnance. In recent years, substantial work has been conducted on
HDLSS asymptotic theory, wherein the sample size n is ﬁxed or n/d → 0
is used as the data dimension d → ∞. For principal component analysis(PCA), Jung and Marron (2009) and Yata and Aoshima (2009) investigatedinconsistency properties for both the eigenvalues and principal component(PC) directions in a sample covariance matrix. Yata and Aoshima (2012)
developed a new PCA method called the noise-reduction methodology andreported consistent estimators for both eigenvalues and PC directions inaddition with the PC scores using this method. Sparse PCA (SPCA) methods have been investigated in several studies. For example, Zou and Hastie(2006), Shen and Huang (2008), and Lee, Huang and Hu (2010) considered a
regularized SPCA (RSPCA) based on L1 penalties under high-dimensionalsettings. Johnstone and Lu (2009) proposed a thresholded SPCA (TSPCA)
and presented a consistency property of the TSPCA when n/d → 0. Further, Shen, Shen and Marron (2013) showed that the PC directions obtainedStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)by RSPCA and TSPCA are consistent when d → ∞ while n is ﬁxed. In addition, Paul and Johnstone (2007) developed an augmented SPCA method
and Ma (2013) proposed an iterative thresholding procedure for PC directions. In this study, we focused on TSPCA under high-dimensional settings.
Suppose that we have a d × n data matrix X = (x1 , ..., xn ), where xi =(xi(1) , ..., xi(d) )T , i = 1, ..., n, are independent and identically distributed(i.i.d.) as a d-dimensional distribution with mean µ and (non-negative
deﬁnite) covariance matrix Σ. We express the eigen-decomposition of Σas Σ = HΛH T , where Λ represents a diagonal matrix of the eigenvalues,
λ1 ≥ · · · ≥ λd (≥ 0), and H = (h1 , ..., hd ) is an orthogonal matrix ofthe corresponding eigenvectors. The sample covariance matrix is given by
S = (n − 1)−1 (X − X)(X − X)T = (n − 1)−1x̄ = n−1∑n∑nTi=1 (xi − x̄)(xi − x̄) , whereTT
i=1 xi and X = x̄1n with 1n = (1, ..., 1) .

Let λ̂1 ≥ · · · ≥ λ̂d ≥0 be the eigenvalues of S, and let ĥj , j = 1, ..., d be the correspondingTeigenvectors ĥj ĥj ′ = δjj ′ , where δjj ′ is the Kronecker delta. Thus, theeigen-decomposition of S is S =∑dTT
s=1 λ̂s ĥs ĥs . We assume that hj ĥj ≥ 0

for all j without the loss of generality. We now consider the n × n dualsample covariance matrix deﬁned by S D = (n − 1)−1 (X − X)T (X − X).
Here, S and S D share nonzero eigenvalues. Let the eigen-decomposition of
S D be S D =∑n−1T
j=1 λ̂j ûj ûj , where ûj denotes an eigenvector correspondingStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)to λ̂j and ûTj ûj ′ = δjj ′ . Furthermore, ĥj can be calculated as ĥj = {(n −1)λ̂j }−1/2 (X − X)ûj .
Johnstone (2001), Baik and Silverstein (2006), and Paul (2007) considered a spiked model for the eigenvalues.

λj (> κ), j = 1, ..., m, are ﬁxed (not depending on d)and λm+1 = · · · = λd = κ.

(1.1)

Here, m represents a ﬁxed positive integer and κ (> 0) represents a ﬁxed
constant. Under (1.1), the asymptotic behavior of the eigenvalues of S
was studied when both d and n increased at the same rate, that is, fromn/d → γ > 0. Details under the Gaussian assumptions were reported
by Johnstone (2001), Johnstone and Lu (2009), and Paul (2007). Further,
Baik and Silverstein (2006) and Lee, Zou and Wright (2010) reported details
under the non-Gaussian but i.i.d. assumptions as in (2.4). For review, theauthors direct readers to Paul and Aue (2014). The condition λm+1 =· · · = λd = κ is strict for the latter part of (1.1). Without assuming that
λm+1 = · · · = λd = κ in (1.1), Bai and Ding (2012) estimated the forward
eigenvalues. However, the former part of (1.1) is a strict condition because
the eigenvalues depend on d, and it is probable that λj → ∞ for the ﬁrst

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)several js when d → ∞. Details were provided by Fan, Liao and Mincheva(2013), Jung and Marron (2009), Onatski (2012), Shen et al. (2016), Wang
and Fan (2017), and Yata and Aoshima (2012, 2013). They considered
spiked models such as

λj = κj dαj (j = 1, ..., m) and λj = κj (j = m + 1, ..., d).

(1.2)

Here, κj (> 0) and αj (α1 ≥ · · · ≥ αm > 0) are ﬁxed constants preserving
the order that λ1 ≥ · · · ≥ λd . For example, Cai, Han and Pan (2020), Shenet al. (2016), Wang and Fan (2017), and Yata and Aoshima (2012) showedthat

λ̂j /λj = 1 + δ/λj + oP (1) as d → ∞ and n → ∞ for j = 1, ..., munder (1.2), d1−2αm /n = o(1), and (2.4), where δ =∑d

s=m+1 λs /(n − 1).

Further details are provided in Appendix E in the online supplementarymaterial. Here, δ = O(d/n) for (1.2); if δ/λj → ∞, λ̂j is strongly inconsistent in the sense that λj /λ̂j = oP (1). Jung and Marron (2009) and Yataand Aoshima (2009) have reported on the concept of strong inconsistency.
Yata and Aoshima (2012) proposed a noise-reduction (NR) methodologythat uses the geometric representation of high-dimensional eigenspaces toStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)overcome the curse of dimensionality. If the NR method is used, λj s can be
estimated by∑tr(S D ) − js=1 λ̂s
λ̃j = λ̂j −n−j−1(j = 1, ..., n − 1).

(1.3)

Yata and Aoshima (2012, 2013) showed the consistency as “λ̃j /λj = 1 +oP (1)” even when δ/λj → ∞. Section 2.3 provides further details. For
PC direction ĥj , under (1.2) and some regularity conditions, Yata andAoshima (2012, 2013) and Shen et al. (2016) showed that Angle(ĥj , hj ) =Arccos{(1 + δ/λj )−1/2 } + oP (1) as d → ∞ and n → ∞ for j = 1, ..., m. If
δ/λj → ∞, ĥj is strongly inconsistent in that Angle(ĥj , hj ) = π/2 + oP (1).
To overcome this inconvenience, Shen, Shen and Marron (2013) showedthat estimators of h1 given by TSPCA and RSPCA have a consistency
property and that they perform equivalently for high-dimensional data.
However, TSPCA is easier to handle as compared to RSPCA. Appendix C
in the online supplementary material provides further details in this regard.
TSPCA can be summarized as follows: Let ĥj = (ĥj(1) , ..., ĥj(d) )T for all j.
Given a sequence of threshold values ζ > 0, we can deﬁne the thresholded

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)entries asĥj(j ′ )ĥj∗(j ′ ) =




0

if |ĥj(j ′ ) | ≥ min{ζ, ĥj max },for j ′ = 1, ..., d,(1.4)

otherwise

where ĥj max = maxs=1,...,d |ĥj(s) |. Let ĥj∗ = (ĥj∗(1) , ..., ĥj∗(d) )T . Then, thethresholded estimator of hj is deﬁned asĥj(ζ) = ĥj∗ /∥ĥj∗ ∥,(1.5)

where ∥ · ∥ denotes the Euclidean norm. Shen, Shen and Marron (2013)
showed the consistency property

Angle(ĥ1(ζ) , h1 ) = oP (1) as d → ∞ for some ζunder (1.2), the Gaussian assumption, and some regularity conditions. Further, they showed consistency even when δ/λj → ∞. However, this estimator depends heavily on the choice of ζ.
We analyzed the microarray data given by Chiaretti et al. (2004) wherein
the dataset comprises 12625 (= d) genes. The dataset consists of two tumor cellular subtypes: B-cell (33 samples) and T-cell (33 samples). B-cell

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

PC2

PC2
Π1
Π2

50

PC2
15

Π1
Π2

50

Π1
Π2

10

5
-50

50

PC1

-60

-40

-20

20

40

60

PC1
-20

-50
-50

-10

10

20

PC1

-5

-10
-100

(i) Conventional PCA(ii) TSPCA with ζ = 0.01

(iii) TSPCA with ζ = 0.05.

Figure 1: Scatter plots of the ﬁrst two PC scores for (i) Conventional PCA,(ii) TSPCA with ζ = 0.01, and (iii) TSPCA with ζ = 0.05.

originally contained 95 samples. We used only the ﬁrst 33 samples to maintain balance in the sample sizes with the T-cell. For the 66 (= n) samples, the scatter plots of the ﬁrst two PC scores are shown in Fig. 1 forthree PCAs: (i) Conventional PCA, (ii) TSPCA with ζ = 0.01, and (iii)
TSPCA with ζ = 0.05. We have that π/2−Angle(ĥ1(0.01) , ĥ2(0.01) ) = 0.154,
and π/2−Angle(ĥ1(0.05) , ĥ2(0.05) ) = 0.266. Further, we observed that theconventional PCA eﬀectively classiﬁes the dataset into two groups using
the ﬁrst PC score. The theoretical clariﬁcation was provided by Yataand Aoshima (2020). The TSPCA with ζ = 0.05 separates them moreclearly than that when using the conventional PCA. However, this may()
not hold consistency as “Angle ĥj(ζ) , hj = oP (1) for j = 1, 2” because
π/2−Angle(ĥ1(0.05) , ĥ2(0.05) ) = 0.266. Thus, TSPCA provides preferable
performance if one chooses a suitable ζ, while it may be inconsistent.
In this study, we investigated TSPCA under high-dimensional settings.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

The contributions of this study are as follows. (I) We propose a new thresholding estimator for PC directions and present its consistency property that
holds freely from threshold values. (II) We propose a shrinkage PC direction
and apply it to clustering.
The remainder of this paper is organized as follows. In Section 2, we
present the asymptotic properties of the NR method under certain conditions. In Section 3, we modify the estimator of the PC directions derived
using the NR method and propose a new TSPCA method. In Section 4, weinvestigate the performance of the proposed TSPCA through simulations.
In Section 5, we present the estimation of the shrinkage PC directions andtheir application to clustering. In Section 6, we investigate the performanceof the estimated shrinkage PC directions in actual data analyses. In the
online supplementary material, we apply the proposed TSPCA to the estimation of the intrinsic component of Σ.

2.

Preliminary

In this section, we lay out basic conditions, assumptions, and asymptotics
for the construction of our PCA methods.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

2.1 Strongly spiked eigenstructures
2.1

Strongly spiked eigenstructures

Aoshima and Yata (2018) provided two disjointed high-dimensional models:
the strongly spiked eigenvalue (SSE) model deﬁned aslim infd→∞

λ21>0
tr(Σ2 )

(2.1)

and the non-SSE (NSSE) model deﬁned byλ21→ 0 as d → ∞.
tr(Σ2 )

(2.2)

Notably, (2.2) is equivalent to “tr(Σ4 )/tr(Σ2 )2 → 0 as d → ∞”. If α1 ≥ 0.5
in (1.2), then the SSE model (2.1) holds. In contrast, if α1 < 0.5 in (1.2),
the NSSE model (2.2) holds. We provide additional examples of the SSE
model in Remark 6 in Section 5.2 and Appendix D in the online supplementary material. The two models are essential for statistical inference ofhigh-dimensional data. We emphasize that it is not possible to ensure theaccuracy of high-dimensional statistical inferences using the SSE model.
The work by Aoshima and Yata (2018) can be referred to for further details. Aoshima and Yata (2018, 2019) proposed data-transformation methods based on the strongly spiked eigenstructures to overcome this inconve-

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

2.1 Strongly spiked eigenstructuresnience. Further, Yata and Aoshima (2020) considered clustering under the
SSE model. The key to these references is the estimation of the strongly
spiked eigenstructures. In this study, we focused on estimating the strongly
spiked eigenstructures.
Let Σ = Σ1 +Σ2 , where Σ1 =

∑m

∑d
T
T
s=1 λs hs hs and Σ2 =
s=m+1 λs hs hs

with an unknown and positive ﬁxed integer m (independent of d). Here,Σ1 is considered an intrinsic part and Σ2 , a noise component. We assume
the following model:(C-i) λ1 , ..., λm are distinct in that lim inf d→∞ (λj /λj ′ − 1) > 0 for 1 ≤ j <j ′ ≤ m when m ≥ 2 and λm and λm+1 satisfy

lim inf
d→∞

λ2m+1
λ2m
>
0
and
→ 0 as d → ∞.
tr(Σ22 )tr(Σ22 )Remark 1. (C-i) is an SSE model. The spiked model (1.2) with αm ≥ 0.5
and κj ̸= κj ′ for 1 ≤ j ̸= j ′ ≤ m satisﬁes (C-i). Aoshima and Yata (2018)provided a method to check whether the SSE model holds or not. Aoshima
and Yata (2018) also provided a consistent estimator of m in (C-i).
Next, we consider a bounded condition for diagonal elements. Let σ(j) =Var(xi(j) ) for all j. Here, σ(j) represents the j-th diagonal element of Σ. Let
A1 =

∑d
T
T
s=m+1 hs hs , where I d represents
s=1 hs hs and A2 = I d − A1 =∑m

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

2.2 Assumptions of high-dimensional distributions
the d-dimensional identity matrix. Let xi,1 = (xi(1),1 , ..., xi(d),1 )T = A1 xi
and xi,2 = (xi(1),2 , ..., xi(d),2 )T = A2 xi for all i, Var(xi,s ) = Σs for s = 1, 2.
Let σ(j),s = Var(xi(j),s ), s = 1, 2, for all j. Notably, σ(j),2 ≤ σ(j),1 + σ(j),2 =σ(j) for all j. We assume the following bounded condition, as necessary:(C-ii) lim inf σ(j),2 > 0 and lim sup σ(j) < ∞ for all j.
d→∞

d→∞

The diagonal elements are typically bounded. Thus, (C-ii) generally
holds. Under (C-ii), σ(j),2 ∈ (0, ∞) as d → ∞ for all j. Here, for function f (·), “f (d) ∈ (0, ∞) as d → ∞” implies that lim inf d→∞ f (d) > 0
and lim supd→∞ f (d) < ∞. Then, tr(Σ22 )/d ≥∑d2
s=1 σ(s),2 /d ∈ (0, ∞) andtr(Σ)/d ∈ (0, ∞) as d → ∞ under (C-ii). Therefore, under (C-i) and (C-ii),

lim sup λj /d < ∞ and lim inf λj /d1/2 > 0 for j = 1, ..., m.
p→∞

2.2

p→∞

(2.3)

Assumptions of high-dimensional distributions

Let X − µ1Tn = HΛ1/2 Z. Then, Z represents a d × n sphered data matrix
obtained from a distribution with an identity covariance matrix. Here, wewrite Z = (z 1 , ..., z d )T and z j = (z1j , ..., znj )T , j = 1, ..., d. E(zij zij ′ ) =0 (j ̸= j ′ ) and Var(z j ) = I n . For convenience, when λj = 0 for some j, we
assume Var(z j ) = I n . Let Mj = Var(zij2 ), j = 1, ..., d. If X is Gaussian,zij s are i.i.d. as the standard normal distribution N (0, 1) and Mj = 2 forStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

2.2 Assumptions of high-dimensional distributions
all j. We consider the following assumptions:(A-i) E(zij2 1 zij2 2 ) = 1, E(zij1 zij2 zij3 ) = 0, E(zij1 zij2 zij3 zij4 ) = 0 for all j1 ̸=
j2 , j3 , j4 ; and lim sup Mj < ∞ for all j.
d→∞Here, (A-i) naturally holds when X is Gaussian. Another example satisfying (A-i) is the case when X has a skew normal distribution (Remark
S4.1 in the work by Aoshima and Yata (2018)). This assumption was provided by Bai and Saranadasa (1996), Chen and Qin (2010), and Aoshima
and Yata (2011). In contrast, Baik and Silverstein (2006), Lee, Zou andWright (2010), and Yata and Aoshima (2012) considered the following assumption:
zi1 , ..., zid are independent (or i.i.d.).

(2.4)

The ﬁrst assumption in (A-i) is milder than (2.4).
Next, we consider a sub-exponential distribution. Let hj = (hj(1) , ..., hj(d) )Tfor all j. From xi(j),2 =∑d1/2
s=m+1 λs hs(j) zis , we note that Var(xi(j),2 zij ′ ) =σ(j),2 under (A-i) for j = 1, ..., d; j ′ = 1, ..., m. Under (C-ii), we consider
the following assumption:

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

2.3 Estimation of eigenvalues and PC directions using the NR method(A-ii) There exist positive (ﬁxed) constants t1 and t2 such that{}lim sup E exp(tx2i(j),2 ) < ∞ for |t| ≤ t1 and all j; andd→∞{}lim sup E exp(txi(j),2 zij ′ ) < ∞ for |t| ≤ t2 , all j and j ′ = 1, ..., m.
d→∞

1/2

If X is Gaussian, xi(j),2 /σ(j),2 follows N (0, 1), and xi(j),2 and zij ′ are independent for all j and j ′ = 1, ..., m, so that E{exp(txi(j),2 zij ′ )} ≤ E{exp(|t|x2i(j),2 +|t|zij2 ′ )} = E{exp(|σ(j),2 t|x2i(j),2 /σ(j),2 )}E{exp(|t|zij2 ′ )} = (1−2|σ(j),2 t|)−1/2 (1−2|t|)−1/2 if |t| < min{1/2, 1/(2σ(j),2 )} for all j and j ′ = 1, ..., m. Thus, when
X is Gaussian, (A-ii) holds under (C-ii). In Appendix B in the online
supplementary material, we consider a milder assumption than (A-ii).

2.3

Estimation of eigenvalues and PC directions using the NR
methodEigenvalue estimation using the NR method is given by (1.3). The second
term in (1.3) is the estimator of δ.
Proposition 1 (Aoshima and Yata, 2018; Yata and Aoshima, 2013). Assume (A-i) and (C-i). Then, it holds for j = 1, ..., m that λ̃j /λj = 1 +Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

2.3 Estimation of eigenvalues and PC directions using the NR method
OP (n−1/2 ) and()(n/Mj )1/2 λ̃j /λj − 1 ⇒ N (0, 1) if lim inf Mj > 0d→∞

as d → ∞ and n → ∞.
Notably, λ̃j is a consistent estimator of λj even when δ/λj → ∞. Thus,
we recommend using λ̃j instead of λ̂j for high-dimensional data. Wang
and Fan (2017) applied λ̃j s in “Shrinkage Principal Orthogonal complEment Thresholding” when d1/2 = o(λm ). Anderson (1963) showed thatn1/2 (λ̂j /λj − 1) ⇒ N (0, 2) as n → ∞ for Gaussian data when d is ﬁxed.
Thus, λ̃j has the same limiting distribution as in Proposition 1 for theGaussian data.
When applying the NR method to the PC direction vector, we obtainh̃j = {(n − 1)λ̃j }−1/2 (X − X)ûj for j = 1, ..., n − 1.

Notably, h̃j = (λ̂j /λ̃j )1/2 ĥj . From Aoshima and Yata (2018), we obtain thefollowing result.
Proposition 2 (Aoshima and Yata, 2018). Assume (A-i) and (C-i). Then,T

it holds for j = 1, ..., m that h̃j hj = 1 + OP (n−1 ) as d → ∞ and n → ∞.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

2.3 Estimation of eigenvalues and PC directions using the NR method
Aoshima and Yata (2018, 2019) used h̃j to develop a data transformation method under the SSE model in the two-sample test and highdimensional classiﬁcation. Here, h̃j is not a unit vector and ∥h̃j ∥2 = λ̂j /λ̃j .
Further, Angle(h̃j , ĥj ) = 0, and therefore, Angle(h̃j , hj ) = Angle(ĥj , hj ).
From Propositions 2, E.1 in Appendix E, and (G.1) in Appendix G in the online supplementary material, under (A-i) and (C-i), it holds for j = 1, ..., mthat∥h̃j ∥2 = 1 + (δ/λj ){1 + OP (n−1/2 )} + OP (n−1 ) and(2.5)∥h̃j − hj ∥2 = (δ/λj ){1 + OP (n−1/2 )} + OP (n−1 )as d → ∞ and n → ∞. Considering Remark E.1 in Appendix E by notingthat 2{1 − (1 + δ/λj )−1/2 } < δ/λj , the norm loss of h̃j is larger than thatof ĥj . However, from Proposition 2, h̃j represents a consistent estimatorof hj considering the inner product even when δ/λj → ∞. Further, from
Propositions 2 and (2.5), under (A-i) and (C-i), there exists a random ddimensional vector ṽ j such that hTj ṽ j = 0,h̃j = {1 + OP (n−1 )}hj + ṽ j and ∥ṽ j ∥2 = (δ/λj ){1 + OP (n−1/2 )} + OP (n−1 )(2.6)

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

2.3 Estimation of eigenvalues and PC directions using the NR method
for j = 1, ..., m. The coeﬃcient of hj in h̃j is asymptotically 1. In contrast,from (2.6), Proposition 1, and Preposition E.1, it holds thatĥj = (1 + δ/λj )−1/2 {1 + oP (1)}h̃j = (1 + δ/λj )−1/2 {1 + oP (1)}(hj + ṽ j ).
(2.7)

The coeﬃcient of hj in ĥj depends on noise δ. The NR estimator h̃j has
an advantage over ĥj by applying the two steps described in (1.4) and (1.5)
to the threshold estimation of PC directions.
We consider the following divergence conditions for d and n:
(⋆)

log d
= o(1) as d → ∞ and n → ∞.
n

Notably, (⋆) holds even when d/n → ∞. Let h̃j = (h̃j(1) , ..., h̃j(d) )T for allj.
Lemma 1. Assume (A-i), (A-ii), (C-i), and (C-ii). Under (⋆), it holds for
j = 1, ..., m and j ′ = 1, ..., d that(h̃j(j ′ ) = hj(j ′ ) + OP−1
(λ−1log d)1/2
j n)

as d → ∞ and n → ∞.

Thus, under the conditions in Lemma 1, we have consistency in thesense that h̃j(j ′ ) = hj(j ′ ) {1+oP (1)} for j ′ satisfying (nλj h2j(j ′ ) )−1 log d = o(1).

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)From (2.7), the elements in ĥj are not consistent unless δ/λj = o(1).

3.

Automatic sparse PCA methodology

In this section, we propose a thresholded estimator of the PC directions
by using the NR method. We emphasize that consistency properties of the
SPCA methods heavily depend on threshold (tuning) values. To overcome
the inconvenience, we propose an SPCA method from (2.5) and (2.6).

3.1

Thresholded estimator of PC directions using the NR method

For the PC direction hj = (hj(1) , ..., hj(d) )T , we arrange the elements hj(1) , ..., hj(d)in the descending order as|hoj(1) | ≥ · · · ≥ |hoj(d) |.

Further,∑d2
s=1 hoj(s) = 1.

(3.1)

Under (C-i), we assume the following condition

for the PC direction as necessary:(C-iii) For j = 1, ..., m, there exists an integer kj∗ (which may depend on

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

3.1 Thresholded estimator of PC directions using the NR methodd), such thatkj∗∑

h2oj(s) → 1 as d → ∞,s=1

lim supd→∞lim inf λj h2oj(kj∗ ) > 0;
d→∞

and

|hoj(kj∗ +1) |< 1 when kj∗ ≤ d − 1.
|hoj(kj∗ ) |Remark 2. When Σ = Γd , h1 is 1d /d1/2 , where Γd is given by (D.1) inAppendix D in the online supplementary material. Then, (C-iii) holds withk1∗ = d and m = 1.
Remark 6 in Section 5.2 and Appendix D present some models thatsatisfy (C-i) and (C-iii). From (G.15) in Appendix G in the online supplementary material, under (C-i) through (C-iii), we note that kj∗ → ∞ andkj∗ /λj ∈ (0, ∞) as d → ∞ for j = 1, ..., m. Now, we consider removing ṽ jfrom h̃j in (2.6) for each j (= 1, ..., m). Notably, ∥ṽ j ∥ ≈ δ/λj from (2.5).
For h̃j = (h̃j(1) , ..., h̃j(d) )T by the NR method, we arrange the elementsh̃j(1) , ..., h̃j(d) in the descending order as h̃oj(1) , ..., h̃oj(d) , where|h̃oj(1) | ≥ · · · ≥ |h̃oj(d) |.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

3.1 Thresholded estimator of PC directions using the NR method∑
Further, ds=1 h̃2oj(s) = ∥h̃j ∥2 . We consider the following thresholded estimator for an integer k ∈ [1, d] ash̃oj (k) = (h̃oj(1) , ..., h̃oj(k) , 0, ..., 0)T

whose last d − k elements are zero; that is, we replace h̃oj(k+1) , ..., h̃oj(d) in(h̃oj(1) , ..., h̃oj(d) )T with 0. Here, we provide the optimal integer k from (2.5)
and (2.6). From ∥h̃j ∥ ≥ 1 w.p.1, there exists a unique integer k̃j ∈ [1, d]such thatk̃j −1∑h̃2oj(s) < 1

and

s=1

For a sequence {As }, we deﬁnek̃j∑h̃2oj(s) ≥ 1.

(3.2)

s=1

∑0

s=1 As

= 0 for convenience, so thatk̃j = 1 if h̃2oj(1) ≥ 1. Notably, ∥h̃oj (k̃j )∥ ≈ 1 and∑d2
s=k̃j +1 h̃oj(s)≈ δ/λjfor a suﬃciently large d. Thus, we can remove ṽ j from h̃j in (2.6) to obtain the following result. We write the elements of h̃oj (k̃j ) as h̃oj (k̃j ) =(h̃oj∗(1) , ..., h̃oj∗(d) )T . Further, we adjust the subscript j ′ of h̃oj∗(j ′ ) as h̃j∗(j ′′ )corresponding to the order of elements in h̃j . Let h̃j∗ = (h̃j∗(1) , ...., h̃j∗(d) )Tfor all j. Let ηj =d → ∞ in (C-iii).

∑d2
s=kj∗ +1 hoj(s)for j = 1, ..., m. Notably, ηj → 0 as

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

3.1 Thresholded estimator of PC directions using the NR method
Theorem 1. Assume (A-i), (A-ii), and (C-i) through (C-iii). Under (⋆),
it holds for j = 1, ..., m that ∥h̃j∗ ∥2 = 1 + OP (ηj + n−1 ) = 1 + oP (1) andTh̃j∗ hj = 1 + OP (ηj + n−1 ) = 1 + oP (1) as d → ∞ and n → ∞.
From Theorem 1, we have the following result.
Corollary 1. Assume (A-i), (A-ii), and (C-i) through (C-iii). Under (⋆),
it holds that()∥h̃j∗ − hj ∥2 = OP ηj + n−1 = oP (1) for j = 1, ..., m;
( 1/2
T
1/2
1/2 )
h̃j∗ hj ′ = OP ηj + n−1/2 min{1, λj /λj ′ } = oP (1) and
( 1/2
)
T
1/2
h̃j∗ h̃j ′ ∗ = OP ηj + ηj ′ + n−1/2 = oP (1)when m ≥ 2 for j, j ′ ≤ m; j ̸= j ′as d → ∞ and n → ∞.
Thus, h̃j∗ has the aforementioned consistency properties even whenδ/λj → ∞ without threshold (tuning) values such as ζ in (1.4). From Theorem 1, we have Angle(h̃j∗ , hj ) = oP (1) under the conditions in Corollary
1.
T

Proposition 3. Assume (A-i) and (C-i). Then, it holds that h̃j hj ′ =OP (n−1/2 min{1, λj /λj ′ }) as d → ∞ and n → ∞ when m ≥ 2 for j, j ′ ≤
1/2

1/2

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

3.2 Automatic SPCA′

m; j ̸= j .
T

T

From Proposition 3, h̃j∗ hj ′ and h̃j hj ′ are of the same order for j, j ′ ≤m; j ̸= j ′ if ηj = O(n−1/2 min{1, λj /λj ′ }) as d → ∞ and n → ∞.
1/2

1/2

Remark 3. We assume that |h̃oj(k̃j ) | > |h̃oj(k̃j +1) |, j = 1, ..., m for the sakeof simplicity. Then, we can simply obtain h̃j∗(j ′ ) by




h̃

if |h̃j(j ′ ) | ≥ |h̃oj(k̃j ) | ,




0

otherwisej(j ′ )h̃j∗(j ′ ) =3.2

for j ′ = 1, ..., d.

Automatic SPCA

The computational cost of the SPCA method is high because it heavily
depends on threshold (tuning) values determined by some cross-validation
or information criteria. One can automatically yield the threshold estimation with a low computational cost because h̃j∗ does not depend on any
threshold (tuning) values such as ζ in (1.4).
We call the new PCA method that uses λ̃j s and h̃j∗ s as the “automatic
SPCA (A-SPCA)”. Proposition 1 provides the details of λ̃j . λ̃j s have theconsistency as “λ̃j /λj = 1 + oP (1)” even when n/d → 0. The A-SPCA
algorithm is given as[Automatic SPCA (A-SPCA)]

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

3.2 Automatic SPCA
(Step 1) Set the dual-sample covariance matrix deﬁned by S D = (n −1)−1 (X − X)T (X − X).
(Step 2) Find the eigenvalues λ̂j , and the eigenvectors ûj of S D .
(Step 3) Calculate λ̃j = λ̂j − {tr(S D ) −∑j

s=1 λ̂s }/(n − j − 1) and h̃j ={(n − 1)λ̃j }−1/2 (X − X)ûj for each j. Estimate the j-th eigenvalueby λ̃j .
(Step 4) Arrange the elements of h̃j = (h̃j(1) , ..., h̃j(d) )T in the descending order as |h̃oj(1) | ≥ · · · ≥ |h̃oj(d) |. Find the integer k̃j such that∑k̃j −1s=1h̃2oj(s) < 1 and∑k̃j2
s=1 h̃oj(s) ≥ 1.

(Step 5) Deﬁne h̃j∗ = (h̃j∗(1) , ..., h̃j∗(d) )T by




h̃

if |h̃j(j ′ ) | ≥ |h̃oj(k̃j ) | ,




0

otherwisej(j ′ )h̃j∗(j ′ ) =for j ′ = 1, ..., d.

Estimate the j-th PC direction using h̃j∗ for each j.
In Appendix A in the online supplementary material, we describe theapplication of A-SPCA to estimate Σ1 . Appendix F presents an R code for
A-SPCA.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)Remark 4. Aoshima and Yata (2018) created a data-transformation method
based on the strongly spiked eigenstructures and proposed a high-dimensional
two-sample test using data transformation. Aoshima and Yata (2019); Ishii,
Yata and Aoshima (2022) proposed high-dimensional classiﬁers using datatransformation. The key to these inferences is the estimation of the strongly
spiked eigenstructures. In future, we apply A-SPCA to these inferences forhigh-dimensional data.

4.

Simulation

In this section, we compare the performance of A-SPCA with the conventional PCA and TSPCA given by (1.5).
First, we set d = 2s , s = 6, ..., 12 and n = ⌈d1/2 ⌉, where ⌈x⌉ denotes
the smallest integer ≥ x. We generate xi , i = 1, 2, ... independently fromNd (0, Σ). Then, we consider the following two cases:(S-i) Σ is given by λ1 = d2/3 , λ2 = d1/2 and λ3 = · · · = λd = 1 together
with h1 = (1, 0, ..., 0)T and h2 = (0, 1, 0, ..., 0)T ;(S-ii) Σ is given by (D.2) in Appendix D in the online supplementarymaterial with α = 0.5, β = 2, d1 = ⌈d2/3 ⌉, d2 = ⌈d1/2 ⌉ and Ωd−d1 −d2 =I d−d1 −d2 .

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)(S-i) satisﬁes (C-i) and (C-iii) when m = 2 and k1∗ = k2∗ = 1, and(S-ii) satisﬁes (C-i) and (C-iii) when m = 2, k1∗ = d1 and k2∗ = d2 . (S-i)
does not satisfy (C-ii) because σ(1) ≥ λ1 h21(1) → ∞. We considered (S-i) tohandle a considerably sparse PC direction.
Next, we set n = 10s, s = 2, ..., 10 and d = 1000. We generated xi ,i = 1, 2, ... independently from the mixture model (5.4) with g = 2 andε1 = ε2 = 1/2 as follows:(S-iii) For i = 1, 2, fi (x; µi , Ψi ) is the probability density function of
Nd (µi , Ψi ), where µ1 = (1, ..., 1, 0, ..., 0)T the ﬁrst ⌈d2/3 ⌉ elements are1, µ2 = −µ1 , Ψ1 = (0.3|i−j|

1/3) and Ψ2 = (0.4|i−j|

1/3

).

Further, λ1 ≈ d2/3 and h1 ≈ µ1 /d1/3 in (S-iii). (S-iii) satisﬁes (C-i) and(C-iii) when m = 1. Details are provided in Remark 6 in Section 5.2. (S-iii)
does not satisfy (A-i). Further details can be found in Section 4.1.1 in thework by Qiao et al. (2010).
Finally, we set d = 2s , s = 6, ..., 12 and n = ⌈d1/2 ⌉. We consider aconsiderably non-sparse and non-Gaussian case for the ﬁrst PC direction
as follows:(S-iv) We generated xi = HΛ1/2 z i , i = 1, 2, ... independently from zij =(yij − 5)/101/2 (j = 1, ..., d) where yij s are i.i.d. as the chi-squared

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

distribution with 5 degrees of freedom. We set Σ = Γd where α = 0.5
and β = 1, and Γd is given by (D.1) in Appendix D.
λ1 = 0.5d + 0.5 and h1 = d−1/2 1d in (S-iv). Further, (S-iv) satisﬁes(A-i), (C-i), and (C-iii) when m = 1 and k1∗ = d. However, (S-iv) does notsatisfy (A-ii). Here, h1 appears to be a non-sparse vector in the sense that
all elements of h1 are nonzero.
We set ζ = 0.01, 0.05 and 0.1 for TSPCA using (1.5) in (S-i) to (S-iii).
Further, ζ = 0.01 is a soft threshold and ζ = 0.1 is a hard threshold. In(S-iv), we set ζ = 0.01, 0.015 and 0.02 for TSPCA because all elementsof h1 are nonzero and small. Thus, ĥj(ζ) with a hard threshold results
in a considerably poor performance in (S-iv). We considered the meansquared error MSE(ĥj ) = ∥ĥj − hj ∥2 (j = 1, ..., m) and took the average of
its outcomes from 2000 independent replications. Similar procedures were
performed for MSE(h̃j∗ ) and MSE(ĥj(ζ) ). We summarized the results in
Figs. 2 and 3.
As expected, we observed that A-SPCA achieved preferable performances in (S-i) to (S-iii). For (S-iv), the conventional PCA or TSPCA
with ζ = 0.01 performed better than A-SPCA because all elements of h1
are nonzero. In addition, λ1 is considerably large because λ1 = O(d). λ̂1is consistent in the sense that λ̂1 /λ1 = 1 + oP (1) (see Appendix E in theStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

MSE
0.7
0.6
0.5
0.4
0.3
0.2
0.1
6

̭
h1

MSE
1.2

̭
h2

h͂ 1
̭*
h1(0.01)

1.0
0.8

h͂ 2
̭*
h2(0.01)h1(0.05)0.6

h2(0.05)h1(0.1)0.4

h2(0.1)

̭
̭

̭
̭

0.2
7

8

9

10

11

12

log2d

6

7

8

9

10

11

12

log2d

(S-i) Nd (0, Σ), d = 2s (s = 6, ..., 12), where Σ has λ1 = d2/3 , λ2 = d1/2 , andλ3 = · · · = λd = 1 together with h1 = (1, 0, ..., 0)T and h2 = (0, 1, 0, ..., 0)T .
̭

MSE

h1

1.2
1.0
0.8
0.6
0.4
0.2

h͂ 1
̭*
h1(0.01)̭

h1(0.05)̭

h1(0.1)

6

7

8

9

10

11

12

log2d

̭

MSE

h2

1.2
1.0
0.8
0.6
0.4
0.2

h͂ 2
̭*
h2(0.01)̭h2(0.05)̭h2(0.1)

6

7

8

9

10

11

12

log2d

(S-ii) Nd (0, Σ), d = 2s (s = 6, ..., 12), where Σ has λ1 ≈ d2/3 and λ2 ≈ d1/2
together with h1 = (1, ..., 1, ..., 0)T whose ⌈d2/3 ⌉ elements are 1 andh2 = (0, ..., 0, 1, ..., 1, 0, ..., 0)T , whose ⌈d1/2 ⌉ elements are 1.
Figure 2: Average mean-squared errors for PC directions in (S-i) and (S-ii).

online supplementary material). However, A-SPCA performed preferably
as d increased, even for the non-sparse case. TSPCA with ζ = 0.015 andζ = 0.02 exhibited poor performances when d is large. This is because all
elements of h1 = (1/d1/2 , ..., 1/d1/2 )T in (S-iv) become close to zero with anincrease in d.
Throughout the simulations, TSPCA depended heavily on the threshold value. In contrast, A-SPCA performed preferably without using anyStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

MSE

MSE

̭
h1

0.8
0.6

h͂ 1
̭*
h1(0.01)

0.4

h1(0.05)̭̭

h1(0.1)

0.2
0.0

n
20

40

60

80

100

̭
h1

1.2
1.0
0.8
0.6
0.4
0.2

h͂ 1
̭*
h1(0.01)̭

h1(0.015)̭

h1(0.02)
6

7

8

9

10

11

12

log2d

(S-iii)(S-iv)(S-iii) Mixture model (5.4) with g = 2, d = 1000, n = 10s (s = 2, ..., 10),λ1 ≈ d2/3 and h1 ≈ d−1/3 (1, ..., 1, 0, ..., 0)T whose ﬁrst ⌈d2/3 ⌉ elements are
d−2/3 .
1/2
(S-iv) xi = HΛ z i ; zij = (yij − 5)/101/2 (j = 1, ..., d) in which yij s are
i.i.d. as the chi-squared distribution with 5 degrees of freedom,d = 2s (s = 6, ..., 12), λ1 = 0.5d + 0.5, and h1 = d−1/2 1d .
Figure 3: Average mean-squared errors for the ﬁrst PC direction in (S-iii)and (S-iv). In the left panel, MSE(ĥ1 (0.1)) is too high to describe when n
is small.
threshold values.

5.

Shrinkage PC directions and its application to clustering

We consider the shrinkage PC directions using A-SPCA, and we apply them
to clustering. Fig. 1 shows that the TSPCA yields a preferable performanceeven though it may not hold the consistency as in Corollary 1. Further, weshow that PCA in the shrinkage PC direction is eﬀective for clustering highdimensional data. We emphasize that the shrinkage PC directions depend
on a parameter unlike in A-SPCA.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

5.1 Shrinkage PC directions
5.1

Shrinkage PC directions

For a given constant ωj ∈ (0, 1] for each j (= 1, ..., m), we consider the PC
shrinkage directions with a cumulative contribution ratio greater than orequal to ωj . For each j (= 1, ..., m), there exists a unique integer kjω ∈ [1, d]such thatkjω −1∑h2oj(s) < ωj

and

s=1

kjω
∑

h2oj(s) ≥ ωj ,

(5.1)

s=1

where hoj(s) are given in (3.1). We assume that |hoj(kjω ) | > |hoj(kjω +1) | forsimplicity. We deﬁne




h

if |hj(j ′ ) | ≥ |hoj(kjω ) | ,




0

otherwisej(j ′ )hjω(j ′ ) =for j ′ = 1, ..., d.

Let hjω = (hjω(1) , ..., hjω(d) )T for j = 1, ..., m. Further, ∥hjω ∥2 =∑kjω2
s=1 hoj(s)(≥ ωj ) if |hoj(kjω ) | > |hoj(kjω +1) |. We seek to estimate the shrinkage PC
direction hjω .
As in (3.2), there exists a unique integer k̃jω ∈ [1, d] such thatk̃jω −1∑s=1h̃2oj(s) < ωj

and

k̃jω
∑
s=1

h̃2oj(s) ≥ ωj .

(5.2)

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

5.1 Shrinkage PC directions
We assume that |h̃oj(k̃jω ) | > |h̃oj(k̃jω +1) | for simplicity. We deﬁne that




h̃

if |h̃j(j ′ ) | ≥ |h̃oj(k̃jω ) | ,




0

otherwisej(j ′ )h̃jω(j ′ ) =for j ′ = 1, ..., d.

Let h̃jω = (h̃jω(1) , ..., h̃jω(d) )T for j = 1, ..., m. Notably, h̃jω = h̃j∗ whenωj = 1. Under (C-i), we assume the following conditions:(C-iii’) For some ﬁxed integers rj (≥ 0), lim supd→∞lim inf λj h2oj(kjω +rj ) > 0,d→∞|hoj(kjω +rj +1) |< 1 and|hoj(kjω +rj ) |and ωj λj → ∞ as d → ∞ for j = 1, ..., m.

From (G.47) in Appendix G in the online supplementary material, under(C-i), (C-ii), and (C-iii’), we note that kjω → ∞ and kjω /(ωj λj ) ∈ (0, ∞)as d → ∞ for j = 1, ..., m.
Theorem 2. Assume (A-i), (A-ii), (C-i), (C-ii), and (C-iii’). Under (⋆),
−1
it holds for j = 1, ..., m that ∥h̃jω − hjω ∥2 = OP {ωj (kjω
+ n−1/2 )} = oP (ωj )as d → ∞ and n → ∞.
Thus, h̃jω is consistent even when δ/λj → ∞. However, we cannot
construct a consistent estimator of hjω using ĥj unless δ/λj = o(1). The
reason is explained in (2.7).
We propose shrinkage PC (SH-PC) scores using the shrinkage PC di-

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

5.2 Application to clustering
rections as(xi − x̄)T h̃jω , i = 1, ..., n(5.3)

for each j. In Section 6, we investigated the performance of SH-PC scoresfor real datasets.
Remark 5. One can use a normalized shrinkage PC direction h̃jω /∥h̃jω ∥in (5.3).

5.2

Application to clustering

For the population distribution of xi , we consider a g (≥ 2)-class mixture
model whose probability density function is given byf (xi ) =g∑

εs fs (xi ; µs , Ψs ),

(5.4)

s=1

where εs ∈ (0, 1), s = 1, ..., g,∑g

s=1 εs = 1 and fs (xi ;µs , Ψs ) represents

the probability density function of a d-variate population Πs with mean µs
and covariance matrix Ψs .
∑g

s=1 εs µs

E(xi ) (= µ) =µs′ )(µs − µs′ )T +∑g

and Var(xi ) =s=1 εs Ψs (= Σ).

∑g−1 ∑gs=1

s′ =s+1 εs εs′ (µs −Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

5.2 Application to clustering[Two-class mixture model]
We assume g = 2. Let µ12 = µ1 − µ2 and ∆12 = ∥µ12 ∥2 . Then, thecovariance matrix of the mixture model is given by Σ = ε1 ε2 µ12 µT12 +ε1 Ψ1 + ε2 Ψ2 . We assume hT1 µ12 ≥ 0 without loss of generality. If

λmax (Ψs )/∆12 → 0 as d → ∞ for s = 1, 2,(5.5)

it holds that

λ1 /(ε1 ε2 ∆12 ) → 1 and Angle(h1 , µ12 ) → 0 as d → ∞,(5.6)

where λmax (M ) denotes the largest eigenvalue of any positive-semideﬁnite
1/2

matrix M . Thus, we can obtain a threshold estimation of µ12 /∆12 using
A-SPCA. Furthermore, for the ﬁrst (true) PC score, it follows thathT1 (xi − µ) =




1/2
 λ1/2+ oP (1)}1 {(ε2 /ε1 )when xi ∈ Π1 ,1/2
 −λ1/2+ oP (1)}1 {(ε1 /ε2 )when xi ∈ Π2

(5.7)

as d → ∞ under (5.5). Hence, one can cluster xi s into two groups based on
the sign of the estimated ﬁrst PC scores, as shown in Fig. 1 (i). Sections 2
and 3 in the work by Yata and Aoshima (2020) provide further details on

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

5.2 Application to clustering(5.6) and (5.7).
Let µ12 = (µ(1) , ..., µ(d) )T . We assume |µ(1) | ≥ · · · ≥ |µ(d) | without lossof generality.
Remark 6. If it holds thatµ(k+1) = · · · = µ(d) = 0, lim inf |µ(k) | > 0 and k ≥ d1/2 for some integer k,d→∞(C-i) and (C-iii) with m = 1 are satisﬁed when maxs=1,2 tr(Ψ2s ) = O(d) andΨs s have NSSE models.
For an integer k, we write xi(k) = (xi(1) , ..., xi(k) )T for all i, Var(xi(k) ) =Σ(k) , and Var(xi(k) ) = Ψs(k) when xi ∈ Πs . Let h1(k) = (h1(1) , ...., h1(k) )T ,µ12(k) = (µ(1) , ..., µ(k) )T , and ∆12(k) = ∥µ12(k) ∥2 . As in (5.6)–(5.7), under
the condition that

λmax (Ψs(k) )/∆12(k) → 0 as d → ∞ for s = 1, 2,(5.8)

the following holds: Angle(h1(k) , µ12(k) ) → 0 as d → ∞. Furthermore, forStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

5.2 Application to clustering
the ﬁrst PC score, it follows that as d → ∞hT1(k) (xi − µ) =




 {λ

1/2
{(ε2 /ε1 )1/2 + oP (1)}max (Σ(k) )} −{λmax (Σ(k) )}1/2 {(ε1 /ε2 )1/2 + oP (1)}when xi ∈ Π1 ,when xi ∈ Π2 .

Thus, the shrinkage (true) scores were consistent with those under (5.8).
If
lim inf ∆12(k) /∆12 > 0 and λmax (Ψs(k) ) = o{λmax (Ψs )},

(5.9)

d→∞

(5.8) is milder than that in (5.5). The second condition in (5.9) often holds
1/2

when k ≪ d. Further, from (5.6), it holds that h1(k) ≈ µ12(k) /∆12 , so thatk∑∆12(k)
h2oj(s) .
≈ ∥h1(k) ∥2 =∆12
s=1

On the basis of (5.1) and (5.9), if lim inf d→∞ ωj > 0, the SH-PC scores
were consistent under milder conditions than those for the PC scores. The
SH-PC scores by h̃jω eﬀectively cluster xi s into two groups (see Fig. 1 (iii)and Section 6.1).

[g (≥ 3)-class mixture model]
Let µ23 = µ2 − µ3 . When g = 3, from Theorem 2 and Lemma 2 inYata and Aoshima (2020), under certain regularity conditions, it holds that

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

5.2 Application to clustering
Angle(h1 , µ12 ) → 0 and Angle(h2 , µ23 ) → 0 as d → ∞. Furthermore, for(true) PC scores, it follows thathT1 (xi − µ) =


[
]


 λ1/2 {(1 − ε )/ε }1/2 + o (1)when xi ∈ Π1 ,[]
 −λ1/2{ε1 /(1 − ε1 )}1/2 + oP (1)1otherwise

1

1

1

P

(5.10)
and


1/2

oP (λ2 )





]
T
1/2 [
1/2
h2 (xi − µ) =
λ
[ε
/{ε
(1
−
ε
)}]
+
o
(1)
3
2
1
P
2






 −λ1/2 [[ε2 /{ε3 (1 − ε1 )}]1/2 + oP (1)]2

when xi ∈ Π1 ,when xi ∈ Π2 ,when xi ∈ Π3
(5.11)

for i = 1, ..., n. Thus, as in the case of g = 2, the SH-PC scores by h̃jω with
lim inf d→∞ ωj > 0 can eﬀectively cluster xi s into three groups (see Section
6.2).
When g ≥ 4, Yata and Aoshima (2020) provided the consistency properties of the PC scores.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

6.

Example: Clustering

We assess the performance of clustering based on SH-PC scores using (5.3).
We choose ω1 = ω2 (= ω, say) in (5.2) as ω = 0.1, 0.2 and 0.5.

6.1

Two-class mixture model

We used microarray data provided by Chiaretti et al. (2004) with d = 12625(see Fig 1). We consider a dataset consisting of 33 samples from Π1 : B-celland 33 samples from Π2 : T-cell. For the mixed 66 samples, we calculated
the ﬁrst two PC scores using PCA, A-SPCA, and SH-PC with ω = 0.1, 0.2,
and 0.5. We have that (k̃1ω , k̃2ω ) = (17, 50), (k̃1ω , k̃2ω ) = (52, 134) and(k̃1ω , k̃2ω ) = (381, 640) for ω = 0.1, 0.2 and 0.5, respectively. Fig. 4 depicts
the scatter plots of the PC scores. We observed that the 66 samples were
perfectly classiﬁed into two groups based on the sign of the ﬁrst PC scoreseven when the cumulative contribution ratio ω is as small as ω = 0.1.
These data can be clustered using the SH-PC score with only 17 signiﬁcantvariables. Details are provided in Section 5.2.
Next, we consider an unbalanced case. The B-cells originally contained
95 samples. We set Π1 : B-cell with 95 samples and Π2 : T-cell with 33
samples. As in Fig. 4, we present the scatter plots of the PC scores inFig. 5. The following hold: (k̃1ω , k̃2ω ) = (88, 13), (k̃1ω , k̃2ω ) = (229, 35),

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

6.1 Two-class mixture model
PC2

SPC2
Π1
Π2

50

-50

PC1

50

Π1
Π2

50

-50

SPC1

50

-50

-50

-100

(i) Conventional PCA(ii) A-SPCA
SPC2

SPC2

SPC2
Π1
Π2

40

Π1
Π2

15
10

20

Π1
Π2

10

5

5
-30

-20

-10

10

20

30

SPC1
-10

-5

5

10

SPC1

-6

-4

-2

2

4

6

SPC1

-5
-20

-5

-10
-40

(iii) SH-PCs with ω = 0.5

-15
-10

(iv) SH-PCs with ω = 0.2(v) SH-PCs with ω = 0.1

Figure 4: Scatter plots of the ﬁrst two PC scores for the dataset comprisingΠ1 : B-cell with 33 samples and Π2 : T-cell with 33 samples.
and (k̃1ω , k̃2ω ) = (977, 223) for ω = 0.1, 0.2, and 0.5, respectively. The
128 samples were eﬀectively classiﬁed into two groups based on the signof the second SH-PC score; this ﬁgure is remarkably diﬀerent from Fig. 4.
Further, ε1 ε2 (µ1 −µ2 )(µ1 −µ2 )T becomes small in such an imbalanced case,and therefore, λmax (Σ1 ) or λmax (Σ2 ) is probably the largest eigenvalue ofΣ because Σ = ε1 ε2 (µ1 − µ2 )(µ1 − µ2 )T + ε1 Σ1 + ε2 Σ2 . The second PC
direction is probably (µ1 −µ2 )/∆1/2 . Thus, 128 samples were classiﬁed into
two groups according to the sign of the second SH-PC scores. Section 4.3
in the work by Yata and Aoshima (2020) for further details.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

6.2 Three-class mixture model
PC2

SPC2
Π1
Π2

50

Π1
Π2

60
40
20

-100

-50

PC1

50

-100

-50

SPC1

50
-20
-40

-50
-60
-80

(i) Conventional PCA
SPC2(ii) A-SPCA
SPC2

SPC2
Π1
Π2

30
20

Π1
Π2

10

Π1
Π2

4
2

5
10
-60

-40

-20

20

40

SPC1

-20

-10

10

SPC1

-10

-5

5

SPC1

-2

-10

-5
-4

-20

-10

-6

-15

-8

-30
-40

(iii) SH-PCs with ω = 0.5(iv) SH-PCs with ω = 0.2(v) SH-PCs with ω = 0.1

Figure 5: Scatter plots of the ﬁrst two PC scores for the dataset thatcomprises Π1 : B-cell with 95 samples and Π2 : T-cell with 33 samples.
6.2

Three-class mixture model

We analyzed microarray data from Bhattacharjee et al. (2001) in which
the dataset comprised ﬁve lung carcinoma types with d = 3312. We used
only three classes: Π1 : pulmonary carcinoids with 20 samples; Π2 : normallung with 17 samples; and Π3 : squamous cell lung carcinoma with 21
samples. As in Fig. 4, we present the scatter plots of the PC scores inFig. 6. The following hold: (k̃1ω , k̃2ω ) = (17, 16), (k̃1ω , k̃2ω ) = (44, 41), and(k̃1ω , k̃2ω ) = (195, 171) for ω = 0.1, 0.2 and 0.5, respectively.
We observed that all samples were eﬀectively classiﬁed into three groups
based on the ﬁrst two PC scores. Further details are provided in (5.10) andStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

PC2

20

-20

SPC2

Π1
Π2
Π3

40

20

40

60

20

PC1

-20

20

-20

-40

10

-10

10

20

60

SPC1

-40

Π1
Π2
Π3

20

40

-20

(i) Conventional PCA
SPC2

Π1
Π2
Π3

40

30 SPC1(ii) A-SPCA

SPC2
5

-5

SPC2

Π1
Π2
Π3
5

10

SPC1

Π1
Π2
Π3

2

-2

-2

2

4

6

SPC1

-2
-10

-5
-4

-20

(iii) SH-PCs with ω = 0.5(iv) SH-PCs with ω = 0.2(v) SH-PCs with ω = 0.1

Figure 6: Scatter plots of the ﬁrst two PC scores for the dataset thatcomprises Π1 : pulmonary carcinoids with 20 samples; Π2 : normal lungwith 17 samples; and Π3 : squamous cell lung carcinomas with 21 samples.

(5.11). Fig. 6 (v) shows that these data can be clustered using the SH-PC
score with (k̃1ω , k̃2ω ) = (17, 16). Here, the 17 variables in h1ω and the 16
variables in h2ω were all diﬀerent. Therefore, these data can be clusteredusing the ﬁrst two SH-PC scores with only 33 signiﬁcant variables.

7.

Conclusion

In this study, we investigated TSPCA in high-dimensional settings. We
proposed a new TSPCA method called automatic SPCA (A-SPCA) anddemonstrated that it exhibits the consistency property under mild conditions free from any threshold values. Therefore, we can quickly obtain aStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)more accurate result at a lower computational cost. Further, we proposed
shrinkage PC directions and applied them to the clustering. We demonstrated the performance of clustering based on the shrinkage PC directions.
We demonstrated that the datasets could be clustered using a set of significant variables.

Supplementary Materials
We give an estimate of the intrinsic part in Σ using A-SPCA, examples ofthe strongly spiked eigenstructures, asymptotic results for A-SPCA under amilder assumption than (A-ii), comparisons between TSPCA and RSPCA,
asymptotic properties of the conventional PCA, an R code for A-SPCA andproofs of the theoretical results in the online supplementary material.

Acknowledgements
We are very grateful to the associate editor and the reviewers for theirconstructive comments. The research of the ﬁrst author was partially supported by Grant-in-Aid for Scientiﬁc Research (C), Japan Society for the
Promotion of Science (JSPS), under Contract Number 22K03412. The research of the second author was partially supported by Grants-in-Aid for
Scientiﬁc Research (A) and Challenging Research (Exploratory), JSPS, un-Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

REFERENCES
der Contract Numbers 20H00576 and 22K19769.

References
Anderson, T.W. (1963). Asymptotic theory for principal component analysis. The Annals of Mathematical Statistics 34, 122-148.
Aoshima, M. and Yata, K. (2011). Two-stage procedures for highdimensional data. Sequential Analysis (Editor’s special invited paper) 30,
356-399.
Aoshima, M. and Yata, K. (2018). Two-sample tests for high-dimension,strongly spiked eigenvalue models. Statistica Sinica 28, 43-62.
Aoshima, M. and Yata, K. (2019). Distance-based classiﬁer by data transformation for high-dimension, strongly spiked eigenvalue models. Annalsof the Institute of Statistical Mathematics 71, 473-503.
Bai, Z. and Saranadasa, H. (1996). Eﬀect of high dimension: By an example
of a two sample problem. Statistica Sinica 6, 311-329.
Bai, Z. and Ding, X. (2012). Estimation of spiked eigenvalues in spiked
models. Random Matrices: Theory and Applications 1, [21 pages].
Baik, J. and Silverstein, J.W. (2006) Eigenvalues of large sample covarianceStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

REFERENCES
matrices of spiked population models. Journal of Multivariate Analysis 97, 1382-1408.
Bhattacharjee A., Richards W.G., Staunton, J., Li, C., Monti, S., Vasa, P.,
Ladd, C., Beheshti, J., Bueno, R., Gillette, M., Loda, M., Weber, G.,
Mark, E.J., Lander, E.S., Wong, W., Johnson, B.E., Golub, T.R., Sugarbaker, D.J., and Meyerson, M. (2001). Classiﬁcation of human lung carcinomas by mRNA expression proﬁling reveals distinct adenocarcinoma
subclasses. Proceedings of the National Academy of Sciences U.S.A. 98,
13790-13795.
Cai. T. T., Han, X., and Pan, G. (2020). Limiting laws for divergent spiked
eigenvalues and largest nonspiked eigenvalue of sample covariance matrices. The Annals of Statistics 48, 1255-1280.
Chen, S. X. and Qin, Y.-L. (2010). A two-sample test for high-dimensional
data with applications to gene-set testing. The Annals of Statistics 38,
808-835.
Chiaretti, S., Li, X., Gentleman, R., Vitale, A., Vignetti, M., Mandelli, F.,
Ritz, J., and Foa, R. (2004). Gene expression proﬁle of adult T-cell acutelymphocytic leukemia identiﬁes distinct subsets of patients with diﬀerent
response to therapy and survival. Blood 103, 2771-2778.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

REFERENCES
Fan, J., Liao, Y., and Mincheva, M. (2013). Large covariance estimationby thresholding principal orthogonal complements. Journal of the Royal
Statistical Society: Series B 75, 603-680.
Ishii, A., Yata, K. and Aoshima, M. (2022). Geometric classiﬁers for highdimensional noisy data. Journal of Multivariate Analysis (Editor’s invited
paper) 188, 104850.
Johnstone, I. (2001). On the distribution of the largest eigenvalue in principal components analysis. The Annals of Statistics 29, 295-327.
Johnstone, I. and Lu, Y. (2009). On consistency and sparsity for principal
components analysis in high dimensions. Journal of the Royal StatisticalSociety: Series B 104, 682-693.
Jung, S. and Marron, J. S. (2009). PCA consistency in high dimension, low
sample size context. The Annals of Statistics 37, 4104-4130.
Lee, S., Huang, J. Z., and Hu, J. (2010). Sparse logistic principal components analysis for binary data. The Annals of Applied Statistics 4,
1579-1601.
Lee, S., Zou, F., and Wright F.A. (2010). Convergence and prediction ofStatistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

REFERENCES
principal component scores in high-dimensional settings. The Annals of
Statistics 38, 3605-3629.
Ma, Z. (2013). Sparse principal component analysis and iterative thresholding. The Annals of Statistics 41, 772-801.
Onatski, A. (2012). Asymptotics of the principal components estimator oflarge factor models with weakly inﬂuential factors. Journal of Econometrics 168, 244-258.
Paul, D. (2007). Asymptotics of sample eigenstructure for a large dimensional spiked covariance model. Statistica Sinica 17, 1617-1642.
Paul, D. and Aue, A. (2014). Random matrix theory in statistics: A review.
Journal of Statistical Planning and Inference 150, 1-29.
Paul, D. and Johnstone, I. (2007). Augmented sparse principal component
analysis for high dimensional data. TechnicalReport, UC Davis.
Qiao, X., Zhang, H.H., Liu, Y., Todd, M.J., and Marron, J.S. (2010).
Weighted distance weighted discrimination and its asymptotic properties. Journal of the Royal Statistical Society: Series B 105, 401-414.
Shen, H. and Huang, J.Z. (2008). Sparse principal component analysis via

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

REFERENCES
regularized low rank matrix approximation Journal of Multivariate Analysis 99, 1015-1034.
Shen, D., Shen, H., and Marron, J.S. (2013). Consistency of sparse PCA inhigh dimension, low sample size contexts. Journal of Multivariate Analysis 115, 317-333.
Shen, D., Shen, H., Zhu, H., and Marron, J.S. (2016). The statistics andmathematics of high dimension low sample size asymptotics. Statistica
Sinica 26, 1747-1770.
Wang, W. and Fan, J. (2017). Asymptotics of empirical eigenstructure forhigh dimensional spiked covariance. The Annals of Statistics 45, 13421374.
Yata, K. and Aoshima, M. (2009). PCA consistency for non-Gaussian datain high dimension, low sample size context. Communications in Statistics- Theory and Methods, Special Issue Honoring Zacks, S. (ed. Mukhopadhyay, N.) 38, 2634-2652.
Yata, K. and Aoshima, M. (2012). Eﬀective PCA for high-dimension,
low-sample-size data with noise reduction via geometric representations.
Journal of Multivariate Analysis 105, 193-215.

Statistica Sinica: Newly accepted Paper(accepted author-version subject to English editing)

REFERENCES
Yata, K. and Aoshima, M. (2013). PCA consistency for the power spiked
model in high-dimensional settings. Journal of Multivariate Analysis 122,334-354.
Yata, K. and Aoshima, M. (2020). Geometric consistency of principal component scores for high-dimensional mixture models and its application.
Scandinavian Journal of Statistics 47, 899-921.
Zou, H. and Hastie, R. (2006). Sparse principal component analysis. Journalof Computational and Graphical Statistics 15, 265-286.

Institute of Mathematics, University of Tsukuba, Ibaraki 305-8571, Japan
E-mail: (yata@math.tsukuba.ac.jp)
Institute of Mathematics, University of Tsukuba, Ibaraki 305-8571, Japan
E-mail: (aoshima@math.tsukuba.ac.jp)