{'fraction_non_alphanumeric': 0.03927573166721473, 'fraction_numerical': 0.01237257481091746, 'mean_word_length': 4.374088800530153, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Homogeneity and Sub-homogeneity Pursuit: Iterative Complement Clustering PCA

Principal component analysis (PCA), the most popular dimension-reduction technique, has been used to analyze high-dimensional data in many areas. It discovers the homogeneity within the data and creates a reduced feature space to capture as much information as possible from the original data. However, in the presence of a group structure of the data, PCA often fails to identify the group-specific pattern, which is known as sub-homogeneity in this study. Group-specific information that is missed can result in an unsatisfactory representation of the data from a particular group. It is important to capture both homogeneity and sub-homogeneity in high-dimensional data analysis, but this poses a great challenge. In this study, we propose a novel iterative complement-clustering principal component analysis (CPCA) to iteratively estimate the homogeneity and sub-homogeneity. A principal component regression based clustering method is also introduced to provide reliable information about clusters. Theoretically, this study shows that our proposed clustering approach can correctly identify the cluster membership under certain conditions. The simulation study and real analysis of the stock return data confirm the superior performance of our proposed methods.


Introduction
Since its introduction, principal component analysis (PCA) (Jolliffe, 2002;Anderson, 2003) has become one of the most popular statistical tools for data analysis in a wide range of areas. Literature on PCA dates back to the early twentieth century (e.g., see Pearson (1901) and Hotelling (1933)). However, with an increasing dimension of data, PCA has now been reconsidered and widely discussed for the purpose of analyzing high-dimensional data. Jolliffe and Cadima (2016) and Fan et al. (2018) reviewed recent developments in PCA for statistical analysis on high-dimensional data, including its sparsity and robustness.
However, are more data together really benefiting statistical analysis? Boivin and Ng (2006) provided a negative answer at the aspect of forecasting. This is due to an increased complexity when more data from different populations are grouped together, as a proportion of data can exhibit a different pattern compared with the rest. More specifically, increased complexity refers to heterogeneity when more data are collected from different populations.
In this study, we consider that information from a particular group of data collected from one population can be divided into two categories. One is shared with other groups and forms the homogeneity within the entire data, while the other is group-specific and is the main source of heterogeneity. This type of heterogeneity can be treated as subhomogeneity, which refers to the homogeneity for a particular group of data.
However, sub-homogeneity may not be identified using traditional estimation methods such as PCA. One reason for this is that the sub-homogeneity for a particular group of data can be relatively small compared with the homogeneity within the entire data because it usually contributes less than the homogeneity to the total variance of all of the data. In such a situation, traditional PCA may regard the sub-homogeneity as negligible compared with homogeneity and ignore this group-specific pattern. Moreover, from the interpretation perspective, principal components that are produced using traditional PCA on all of the data do not target a specific group (e.g., information on which component corresponds to which group of data is not known). In previous studies, the discussion of PCA has mainly focused on the large eigenvalues, which correspond to the homogeneity in this study. For example, Johnstone (2001) studied the asymptotic distribution of the largest eigenvalue of PCA. Recently, Cai et al. (2017) extended the discussion to the asymptotic distribution of the spiked eigenvalues, while Morales-Jimenez et al. (2018) studied the asymptotics for the leading eigenvalues and eigenvectors of the correlation matrix in the class of spiked models.
None of them have considered the existence of sub-homogeneity within the data.
To the best of our knowledge, the sub-homogeneity has not been well discussed, but it can be very important in high-dimensional data analysis. The following example of analyzing stock returns from different industries is used to explain the importance of finding sub-homogeneity. As stated in Fama and French (1997), stock returns from various industries can have varying performance over time, although the homogeneity (e.g., market return) can be deemed to be driven by some common economic variables. In a situation in which a vast number of individual stock returns are collected together for statistical analysis, traditional dimension-reduction techniques such as PCA may be able to capture the market effect but can fail to identify the sub-homogeneity within each industry (industry-specific pattern). This is because some industry-specific components may have much smaller variance than the market component and are highly likely to be omitted by PCA. However, these industry-specific components may be very important in capturing the movement of the stock returns within the industry. This loss of information may result in a very poor forecast of stock returns for some companies in which sub-homogeneity exists. In addition, although it would be interesting to study which industry has a larger industry-specific effect on stock returns, traditional PCA performed on the entire data does not allow us to draw such conclusions. Therefore, this study aims at identifying both homogeneity and sub-homogeneity in high-dimensional data analysis.
The sub-homogeneity in several parts of the data can be identified and estimated by dividing the whole data set into several groups. However, in most situations, the group structure is not known in advance. Therefore, a clustering method can be used to group the data at the first stage, and PCA can then be applied to identify the sub-homogeneity within each group. Similarly, Liu et al. (2002) performed PCA in different blocks of variables, while Tan et al. (2015) penalized the variables in each group differently when using the graphical lasso. Both studies used hierarchical clustering to group the variables to take into account the heterogeneity. However, when homogeneity exists across different groups, the sub-homogeneity in each group cannot be successfully identified. This is mainly because the homogeneity shared by most groups can dominate the sub-homogeneity so that the data from different groups all seem to be highly correlated. In this situation, the subhomogeneity is masked by the homogeneity, and these clustering methods tend to group all individuals (variables) into one cluster. Therefore, the homogeneity must be correctly identified and removed before we can successfully group the individuals and discover the sub-homogeneity of each group.
On the other hand, the estimation of homogeneity can often be problematic in a highdimensional setting. This is because of the inconsistency of PCA estimates when the number of variables p is comparable with or greater than the sample size n, as discussed in Johnstone and Lu (2009). Motivated by the group structure of the data, we suggest first clustering the individuals into groups and then performing a traditional approach such as PCA on the level of groups, followed by a second layer of PCA that extracts common information shared by each group. This approach can improve the estimation accuracy of the homogeneity because the components are now extracted from groups in which the dimension has already been reduced (to the number of variables in a group), while the traditional dimension-reduction method (e.g., PCA) is performed on the full dimension p.
This can be viewed as an effective way to alleviate the potential problem caused by the curse of dimensionality. Other studies have modified the traditional approaches to deal with the curse of dimensionality. For example, Johnstone and Lu (2009) suggest that some initial dimension-reduction work is necessary before using PCA, as long as a sparse basis exists. Further, Zou et al. (2006) introduced sparse PCA, which uses the lasso (elastic net) to modify the principal components so that sparse loadings can be achieved. In addition, a review of some sparse versions of PCA can be found in Tibshirani et al. (2015), while some general discussions about the blessing and curse of dimensionality can be found in Donoho et al. (2000) and Fan et al. (2014).
To conclude, homogeneity must be removed to correctly identify the cluster structure and sub-homogeneity, but the cluster structure of the data is used to accurately find the ho-mogeneity. Therefore, we introduce a novel "iterative complement-clustering principal component analysis" (CPCA) to iteratively estimate the homogeneity and sub-homogeneity.
Details of CPCA are provided in Section 4.
The contributions of this study can be summarized as follows. First, we propose CPCA to identify homogeneity and sub-homogeneity and handle the interaction between them when the whole data set exhibits a group structure. Second, our proposed estimation method not only correctly captures the sub-homogeneity, but also provides very reliable cluster information (e.g., which part of the data is from the same group), which can be useful in understanding and explaining the data structure. Third, inspired by Chiou and Li (2007), we develop a leave-one-out principal component regression (PCR) clustering method that can outperform the hierarchical clustering method used in previous studies.
In addition, we theoretically illustrate that if the sub-homogeneity from different clusters is distinct, our proposed clustering procedure can effectively separate the variables from different clusters. This is also numerically confirmed by the simulation study and real data analysis. Details of the proposed clustering method are provided in Section 4.
The rest of this study is organized as follows. A low-rank representation of the data that captures both homogeneity and sub-homogeneity is introduced in Section 2, and some related PCA methods are discussed in Section 3. In Section 4, a novel estimation method called CPCA is proposed, followed by a discussion of more details of the algorithm. Section 5 demonstrates and explains the effectiveness of the proposed clustering method. Extensive simulations, along with two applications (PCR and covariance estimation) of our proposed method are provided in Sections 6 and 7. Section 8 analyzes a stock return dataset using our method. Lastly, Section 9 concludes the study.

Homogeneity and Sub-homogeneity
Considering the data x i = (x 1i , . . . , x pi ) ∈ R p from i th observation with p variables, the singular value decomposition (SVD) of the data x i can be written as: where g ik is defined as k th principal component score and φ k denotes a p × 1 eigenvector for g ik . Traditional PCA summarizes the data using the first k principal components, and it treats u i as noise because g ik for k = m + 1, . . . , min (n, p) has smaller variance. However, under certain conditions, some of the information contained in u i can be useful in prediction or forecasting problems, particularly for one or more specific groups of data. Therefore, when the data exhibit a group structure, we propose the following low-rank representation for the data to capture both homogeneity and sub-homogeneity: . . , n, j = 1, . . . , J. (2) The first line in (2) measures the homogeneity among all variables from all groups, where g ik , k = 1, . . . , r c is k th principal component, which we call the common component, and φ k is its corresponding eigenvector. The first part of u i in (2) consists of J cluster-specific components from which the sub-homogeneity of the data is derived. Assuming p variables can be split into J clusters, f ih , h = 1, . . . , r j measures the within-cluster principal components for cluster j. The within-cluster eigenvector with dimension p × 1 has the form of h defines a p j × 1 vector and p j is the number of variables in cluster j so that J j=1 p j = p. That is, the values of γ (j) h for variables that do not belong to cluster j are zero. This implies that after removing the effect of the common principal components, variables from different clusters are uncorrelated (e.g., u i has a block-diagonal covariance structure). In the second part of u i , (j) i is simply a p jdimensional error that has variance σ (j) 2 , and I (j) is a diagonal matrix, but the diagonals for variables that do not belong to cluster j are zero.
We further define g i = (g i1 , . . . , g irc ) and f ir j ) , j = 1, . . . , J as the vector forms of the common components and cluster-specific components, and the eigenvector matrices as Φ (p×rc) = (φ 1 , . . . , φ rc ) and Γ (j) r j ), j = 1, . . . , J, respectively. Without loss of generality, we assume E( i , j = 1, . . . , J are mutually uncorrelated. In addition, we also impose the following canonical condition for the identification of both the homogeneity and sub-homogeneity, Under the data structure given in (2), the population covariance matrix is given by a low-rank plus block-diagonal representation If we denote the data X (n×p) = (x 1 , . . . , x n ) , the common components and clusterspecific components in a matrix form G (n×rc) = (g 1 , . . . , g n ) and F (j) n ) , respectively, data representation (2) can also be presented in the following matrix form In general, there are two goals in using representation (2). In the presence of an unknown cluster structure, the first goal is to cluster the variables (i.e., determine which variables are in the same group). From the interpretation perspective, it is interesting to know how variables are clustered and which variables belong to the same cluster. The second goal is to correctly estimate both the common components g i and the cluster-specific components f (j) i . These components serve as a low-rank representation of the data and can therefore be used for further applications. One obvious application is to estimate the covariance matrix Σ and the precision matrix Σ −1 of X. The use of both g i and f (j) i captures both the low rank and the block-diagonal representation of Σ, which results in a more efficient estimation compared with using g i only. Another important application is PCR. In some situations, it is a cluster-specific component f to explain which cluster of variables has a greater effect in predicting the response. We will discuss more about these two applications using simulated data in Section 7.

Relationship with Existing PCA Methods
Our proposed data representation (2) should be used to perform dimension reduction of the data with a cluster structure because it captures both the common effect (homogeneity) and the cluster-specific effect (sub-homogeneity). Consequently, this method can be viewed as an extension of many other widely used dimension-reduction methods. Three of these methods are discussed below.
• Case 1: When there is no cluster-specific effect, for example, f (j) i = 0, j = 1, . . . , J and σ (j)2 ≡ σ 2 , representation (2) simply reduces to the well-known PCA where g ik can be found by the principal component with k th largest eigenvalue and φ k is the corresponding eigenvector. However, in the presence of a small cluster-specific effect, traditional PCA generally only identifies the large common effect and treats the rest as noise, disregarding the sub-homogeneity within the data.
• Case 2: When the common components do not exist (e.g., g i = 0), representation (2) demonstrates a well-studied block-diagonal covariance structure of the data, as discussed in Liu et al. (2002) and Tan et al. (2015). Thus, representation (2) can be reduced to In this case, hierarchical clustering is usually used to group the variables with high correlations, and PCA is then performed on each group of variables to estimate the cluster-specific components f (j) i . In many situations, this cluster structure is masked by a dominant common effect. Ignoring this common effect will result in a nonidentifiable cluster structure.
• Case 3: Fan et al. (2013) and Li et al. (2018) consider the following low-rank representation, in which the error covariance matrix is assumed to be cross-sectional dependent after the common components have been taken out and where they assume that the covariance of u i , Σ u , is sparse and propose a method called Principal Orthogonal complEment Thresholding (POET) to explore such a high-dimensional structure with sparsity. Li et al. (2018) used weighted PCA to find a more efficient estimator of common components g i . Hong et al. (2018) and Deville and Malinvaud (1983) also applied weighted PCA to estimate the covariance matrix when heteroscedastic noise of samples or variables exists. However, in this study, we consider u i following a cluster structure with a block-diagonal covariance and aim to find its low-rank presentation, which we call sub-homogeneity. In addition, we propose iteratively estimating the common components and cluster-specific components.

Estimation Methods
Correctly estimating g i and f (j) i poses many challenges. First, PCA, which is widely used to estimate the common component g i , often performs very poorly when p is much larger than n. The group structure of the variables motivates us to separate the data according to the clusters and then extract the common information from each cluster to determine the common components. This is less influenced by the curse of dimensionality because each cluster of the data results in a lower dimension of variable p j . However, to accurately identify the clusters, we must remove the common effect and then perform the clustering method based on the complement (x i − rc k=1 g ik φ k ), but the common components and its eigenvectors are not known in advance. This inspires us to propose a novel iterative method, CPCA, to cluster the variables and estimate the components simultaneously. The details of the method are described in Algorithm 1 and a flowchart that summarizes our estimation method is presented in Figure 1.

Initial
Step: (a) Perform PCA directly on the entire data X and select the number of components using a ratio based estimator defined by Lam and Yao (2012). A brief discussion on selecting the number of components is in Remark 1. The resulting eigenvectors Φ 0 and principal components G 0 are served as the initial estimates of Φ and G. Then, find the initial complement X c 0 = X − G 0 Φ 0 (b) Perform hierarchical clustering for X c 0 based on a similarity matrix given by the absolute value of the empirical correlation matrix X c 0 . The obtained clusters C (j) 0 , j = 1, . . . , J 0 are served as the initial clusters.

Iterative
Step: for s = 1, 2, . . . s ) and perform a further step of PCA on Ψ s . Define the principal components as G s and their corresponding eigenvectors as Φ s . Then, compute the updated complement X c s = X − G s Φ s . Details for finding the eigenvectors Φ 1 can be found in Appendix A.
(b) Perform leave-one-out clustering for variables in X c s using PCR (more details of this clustering method are discussed in Remark 3): (i) For k, k = 1, . . . , p, leave k th variable out of X c s .
(ii) Group the rest of the variables of X c s based on C (iii) Fit J s−1 PCRs by using k th variable of X c s as the response and F (j) s−1 as the predictor for each j = 1, . . . , J s−1 , respectively.
(iv) Compute the sum of squared residuals (SSR) for each J s−1 PCR model.
Assign k th variable to the cluster with the minimum SSR. Update the cluster (v) Repeat (i)∼(iv) for each k and denote the updated clusters as C

Final
Step: (a) Repeat (a) in the Iterative Step, but using cluster C and denote the obtained cluster-specific components as F Therefore, CPCA produces these required outputs: the final clusters of p variables . . , J f , the final estimate of the common components G f and the clusterspecific components F (j) f , j = 1, . . . , J f , along with their eigenvectors. Some details and discussions of CPCA are provided in the following remarks.

Remark 1 [Selecting the Number of Components]:
In this work, we use a ratio based estimator (Lam and Yao, 2012) to select the number of common components and estimate the homogeneity. The ratio based estimator of the number of principal components r c is defined as where λ 1 ≥ · · · ≥ λ p are eigenvalues of the sample covariance matrix of data X and r ≤ R ≤ p is a constant. In practice, many researchers may prefer to select the number of principal components by either setting a threshold on the proportion of total variance explained, which may accidentally include some errors into the principal components, or choosing the numbers where the largest drop of eigenvalues exists (i.e. based on the scree plot). Nonetheless, when the number of common components is greater than 1 (e.g., r c > 1), the discrepancy of the eigenvalues (from the sample covariance) among the common components can be accidentally greater than the one between the common components and the cluster-specific component, especially when some groups have a relatively larger cluster-specific effect (e.g., variables in one group are on a large scale or highly correlated).
In this situation, selecting the number of components using the largest drop in eigenvalues in Initial Step (a) and Iterative Step (a) may lead to underestimating the number of common components, while the ratio-based estimator r c is less affected by the scales of variables. In addition, to consider a broader class of sub-homogeneities, we also allow for sub-homogeneities within a specific cluster to have different order of variances. To correctly estimate the number of principal components that may have different orders of variances in the same group, in Iterative Step 2 (a) and (b), we use the iterative estimation procedure purposed by Peng et al. (2020), which is a generalization of the original ratiobased estimator. The exact procedure is lengthy and is not the focus of the present work, hence is omitted here. More details of the iterative estimation procedure of the ratio-based estimator can be found in Section 2 of Peng et al. (2020).

Remark 2 [Iterative Step]:
One of the key contributions of our algorithm is iteratively estimating the common components and cluster-specific components. Directly applying PCA on the entire data X generally leads to a very poor estimate of the common components because a large p may blur the spike structure of the sample covariance matrix. Motivated by the group structure of the data, we utilize Iterative Step (a) to first cluster variables and then perform PCA on the level of clusters, followed by a second layer of PCA that extracts common information shared by each group. Iterative Step (b) is then implemented to update and improve the clusters information. This is more effective in estimating the common components because PCA is performed in a smaller p case, which is also numerically shown in Figure 2. The data X used here are generated by simulation Example 2. Figure 2 (a) presents the correlation plot for the original data X. We observe that the cluster structure is masked by the common effect. Figure 2 (b) demonstrates that after the common effects estimated in the Initial Step are removed, it is still difficult to perceive the cluster structure. However, if we remove the common effects estimated in the Final Step, the block-diagonal structure is clear, implying a prominent heterogeneity within the data. Therefore, iteratively estimating the common components is advantageous. Initial Step (b), an ordinary hierarchical clustering method using average-linkage as dissim-ilarity between clusters is performed. As discussed in Bühlmann et al. (2013), we choose the number of clusters using the largest increment in height between iterations proceeded in a agglomerative way, such that J 0 = argmax j (h j+1 − h j ).
In Iterative Step (b), we propose a new clustering method called "leave-one-out PCR clustering" (LOO-PCR clustering), which is another key contribution of CPCA. The underlying idea of this new clustering method is that if one variable belongs to a cluster, it should be well-predicted by the principal components extracted from that cluster. This is more aligned with our model set-up than hierarchical clustering based on correlations. A similar idea of leaving one out is used for functional clustering, as discussed in (Chiou and Li, 2007). In addition, hierarchical clustering can perform poorly when the dimension of variables p is larger than n. In step (iv), we assign k th variable to the cluster that achieves the minimum SSR. However, when the minimum SSR is larger than a threshold τ , we treat the k th variable itself as a cluster, because in such a situation, this variable cannot be wellpredicted by any clusters of variables. Hence, the number of clusters is driven by the data and can vary in each iteration. This suggests that our proposed method is more flexible than hierarchical clustering. We set τ = 0.95, and Example 4 in the simulation studies demonstrates that our proposed LOO-PCR clustering, served as a clustering method itself, outperforms hierarchical clustering in a large p small n case. Identifiability of LOO-PCR clustering is also demonstrated theoretically in Section 5.
In Iterative Step (c), we stop the iteration when the clusters C (j) s converge. In our algorithm, we adopt the adjusted rand index (ARI) (Rand, 1971;Hubert and Arabie, 1985) between the clusters C (j) s and the one in the previous iteration C (j) s−1 as the stopping criterion. The ARI is a corrected version of the rand index that measures the similarity between two clusters of the same data using the proportions of agreements between these two partitions. The correction is achieved by subtracting the rand index by its expected value. The ARI can yield a maximum value of 1, and a high value implies high similarity.
When the ARI is above a certain threshold η, we stop the iteration. In this study, we use η = 0.95.

Identifiability of LOO-PCR Clustering Approach
Achieving an accurate sub-homogeneity pursuit relies on the effectiveness of the clustering procedure. In this section, we investigate the proposed LOO-PCR clustering method theoretically. More specifically, we examine the identifiability of cluster membership in our proposed method.
Consider one random variable X mi that belongs to cluster l. Based on the variable's structure in cluster l, there exists m ∈ {1, 2, . . . , p}, such that: We now consider another cluster d with components g ik , k = 1, . . . , r c and f (d) ih , h = 1, . . . , r d . When applying our proposed LOO-PCR clustering method, we regress u ih , h = 1, 2, . . . , r d and measure its goodness of fit to determine whether the random variable X mi belongs to cluster d, such that: where β 1 , . . . , β r d are coefficients and ζ (d) mi is the error. Naturally, we expect that features (e.g., principal components) from cluster d cannot explain u (l) mi sufficiently. Based on (10) and (11), intuitively, a large discrepancy between principal components f Theorem 1 For any cluster d and any random variable X mi from another cluster l, we have the following evaluation: where ih , and G is an n × r c matrix with (i, k)-th element being g ik . Here, M F (d) and M G are M F (d) and M G , but with F (d) and G replaced by F (d) and G, respectively.
A brief proof of Theorem 1 is provided in Appendix B.
Remark 4 It is expected that ζ (d) m will be large when cluster d and cluster l are different. When l = d, the first term on the right-hand side of (12) is equal to In contrast, when l = d (e.g., F (l) 2 and F (d) 2 are distinct), the first term on the righthand side of (12) Meanwhile, the second term on the right-hand side of (12) is determined by the estimation accuracy of F (d) and G. It is related to the dimension of cluster d (e.g., p d ), the sample size n and the total number of clusters J. Bai and Ng (2002) provided the rate of convergence for the projection matrices of principal components. Fan et al. (2013) and Fan et al. (2018) also studied the properties of principal components for high-dimensional data. In this study, we do not pursue the exact expression of α np d and γ nJ . However, given that the homogeneity M G and sub-homogeneity F (d) can be estimated accurately, the first term on the right-hand side of (12) is expected to dominate ζ (d) m 2 . This implies that two different clusters d and l are identifiable, given that F (l) 2 has a higher order than ζ (d) m 2 . As a general discussion, Theorem 1 shows that when the sub-homogeneity from the different clusters is distinct, our clustering procedure can effectively separate the variables from the different clusters.

Simulation Studies
In this section, we conduct various simulation studies to investigate the performance of our proposed CPCA compared with traditional PCA under different simulation settings.
In terms of the cluster-specific components, we set J = 5, and for each of the five clusters, we consider r j ≡ 2 and generate λ (j) h ∼ N (θ j,h , 1), where θ j,h represents the mean level of each component's variance in each cluster. In specific, we set θ 1,1 = θ 1,2 = 25 for the first cluster, θ 5,1 = θ 5,2 = 5 for the last cluster, while θ j,1 = 25 and θ j,2 = 5 for j = 2, 3, 4, i.e., the rest three clusters. That is, each of the five clusters has the same number of cluster-specific components, while the components from different clusters are allowed to have different levels of variances. Among those five clusters, the components in the first cluster have the largest variances, those in the last cluster have the smallest variances, while the two components in the rest clusters are designed to have different levels of variances. Then, by setting p j ≡ 20 for each cluster, we end up with a total number of variables p = 100. According to (2), Z (j) ih for h = 1, 2 and (j) i are simulated independently for each cluster j, with σ (j) = 0.5 for j = 1, 2, ..., 5 so that 2n observations are constructed based on (2). Then, the first n observations are served as the training sample X train = (x 1 , . . . , x n ) and the last n observations are designated as the testing sample X test = (x n+1 , . . . , x 2n ). In this example, n = 50.
• Example 2 (n = 30, p = 100). This example is identical to Example 1, except that the sample size is decreased from 50 to 30. We will show that our CPCA is more reliable than other competing methods when the sample size is relatively small.
• Example 3 (n = 50, p = 200). In this example, we consider data consisting more clusters, that is, we consider J = 10 clusters compared with five in Example 1, leading to a number of variables p = 200 given p j ≡ 20. We maintain r j ≡ 2 but double the number of clusters with exactly the same setting of λ (j) h ∼ N (θ j,h , 1). In other words, components in the first two clusters have the largest variances, while those in the last two clusters have the smallest variances. And similarly, the rest six clusters have two levels of variances on their components as described in Example 1. Accordingly, (j) i is simulated independently for each of the 10 clusters, with σ (j) = 0.5 for j = 1, . . . , 10.
The settings for the common components remain the same as those in Example 1.
• Example 4 (n = 30, p = 100). This example is identical to Example 2, but without the common effect, which is a special case of our proposed data representation as in (7). Correspondingly, the response is generated as: In this example, we do not estimate the common effect, and the clustering method is directly applied to the data. Therefore, the purpose of considering this example is to demonstrate that our proposed LOO-PCR clustering (Iterative Step (b) and (c) in CPCA), served as a clustering method itself, can outperform the hierarchical clustering.
For all four examples, we repeat the simulation procedure 100 times and investigate the clustering and recovering accuracy of the following methods: • PCA: The classical PCA method is performed directly on the whole data set without clustering variables. In addition, clustering accuracy is not considered for this method. The number of principal components is determined according to the ratiobased estimator (9) in Remark 1.
• CPCA I: This procedure is simply CPCA, but without the Iterative Step. First, the complement X c 0 = X − G 0 Φ 0 is found, and then clusters C (j) 0 are obtained using hierarchical clustering for X c 0 . Variables of X c 0 are clustered, and the cluster-specific components F (j) 0 are finally obtained using PCA on each cluster.
• CPCA F: This is exactly the estimation procedure described in CPCA. Different to CPCA I, the clusters and components are estimated iteratively.
The main purposes of this comparison are to demonstrate that: 1) traditional PCA fails to capture the cluster-specific components while our proposed CPCA can successfully identify and capture the cluster-specific components; 2) our proposed iterative estimation of the clusters and components outperforms the initial estimation.
In this study, the recovering accuracy is measured by the mean squared recovering error (MSRE). In this study, the MSRE is defined as: where X test is the recovered testing sample computed using the testing sample and the eigenvectors (Φ, Λ (j) ) estimated from the training sample X train . We prefer a lower MSRE because it indicates a better recovering of the data.

Simulation Results
The boxplots of the total number of principal components selected (No.PCs) and MSRE for Example 1,2 and 3 are presented in Figure 3 and Table 1. As shown, CPCA I (blue boxes) achieves a consistently higher number of principal components than its iterative counterpart (red box) and the true number of principal components (dotted line), which implies the advantages of using the iterative estimation. In the initial clustering, the common effects are estimated without partitioning the variables, while the iterative clustering partitions the variables in the previous step to estimate the common components, resulting in a more accurate estimation. This is also stated in Remark 2.
From the recovering accuracy perspective, it is worth noting from Figure 3 and Ta Lastly, we investigate Example 4, which does not take into account the common effect.
In this example, CPCA I and CPCA F simply reduce to performing PCA on each cluster of variables based on hierarchical clustering and LOO-PCR clustering, respectively. From the last row of Figure 3 and the last column of Table 1, we observe that applying PCA after clustering the variables demonstrates better performance than applying PCA directly to the data. Compared with CPCA I, CPCA F achieves a better estimate of the number of PCs, implying that our proposed LOO-PCR clustering outperforms hierarchical clustering. We consider a PCR model with an univariate response y i for i th observation as: f (j) i β j + e i , i = 1, . . . , n, j = 1, . . . , J. where α is a r c ×1 vector representing the regression coefficients of the common components g i , β j is a r j ×1 vector denoting the regression coefficients of j th cluster-specific components f (j) i and e i is simply the error term with mean 0 and variance θ 2 .
Therefore, after the the common components g i and cluster-specific components f (j) i are estimated via CPCA, we can fit a PCR to predict the response. However, in many situations, only a few groups are useful in predicting the response. To investigate which clusters of variable have an impact on predicting the response, we utilize the group lasso (Yuan and Lin, 2006) using the components produced by CPCA as the coviarates to estimate the regression coefficients as in (15), where λ is the tuning parameter that controls the sparsity of the regression coefficients.
Using this group lasso penalty, β j for some j will be shrunk to zero exactly so that we can identify which clusters of variables are important in predicting the response.
To better demonstrate the performance of our method, we conduct similar simulation studies as in Section 6. We consider the following three examples: • Example 1 (n = 50, p = 100). We generate X using the same settings as Example 1 in Section 6.1. Then, we simulate the response y i according to (15) by setting regression coefficients α = (1, 1, 1), β = (β 1 , . . . , β 10 ) = (0, 0, . . . , 0, 25, 25) and standard deviation θ = 1. That is, only the common effect and the cluster-specific effects for the last cluster are important in predicting the response. This regression setting is interesting because the cluster-specific components with smaller variance (e.g. the last cluster) are highly likely to be omitted or estimated poorly in traditional PCA, but they can sometimes be very important in predicting the response.
• Example 2 (n = 30, p = 100). This example is identical to Example 1 above, except that the sample size is decreased from 50 to 30.
• Example 3 (n = 50, p = 200). We generate X using the same settings as Example 3 in Section 6.1. That is, we consider J = 10 clusters in this example. In regard to the PCR, we set regression coefficients α = (1, 1, 1) again and consider β = For all three examples, we repeat the simulation procedure 100 times and investigate prediction accuracy of the methods mentioned in Section 6.1. The MSPE is used to determine whether the selected components can accurately predict the response: where y i,test is computed using (15) with the regression coefficients (γ, β) and the eigenvectors Φ, Λ (j) estimated from the training sample.
The first panel of Table 2

Covariance and Precision Matrix Estimation
Another useful application of CPCA is to estimate Σ as in (4). Using g i and f (j) i leads to a better estimation of Σ compared with using g i only, because g i and f (j) i can capture both the low rank and the block-diagonal representation of Σ. Therefore, we estimate Σ using clusters, the common components, and the cluster specific components, along with their associated eigenvectors estimated via CPCA. In this section, we numerically illustrate the advantage of using CPCA to estimate the covariance matrix in comparison with other traditional PCA methods. Again, we conduct three simulation examples mentioned in Section 7.1 and utilize the Euclidean distance (ED) between the estimated covariance and population covariance, Σ−Σ 2 F , to measure the performances of different methods. Recall Σ is generated from (4): In this section, we add another prevalent covariance estimation method POET (Fan et al., 2013) as aforementioned into our comparison.
The second panel of Table 2 and the second column of Figure 4 demonstrate the mean and standard errors and boxplot of Σ − Σ 2 F for all three examples, computed using methods discussed before. From these results, we see that POET performs better than PCA but worse than CPCA I and CPCA F, indicating that the sparsity structure implemented in POET can partly capture the sub-homogeneity, while not in a very efficient way. Among these methods, both CPCA I and CPCA F can achieve the lowest mean and standard error of Σ − Σ 2 F because it can best identify both homogeneity and sub-homogeneity to accurately estimate Σ.

Real Data Analysis
For further illustration, we analyze a stock return data set using the proposed CPCA and compare it with traditional PCA. In particular, we are interested in the performance of our proposed CPCA on clustering stock returns of companies from various industries.
Moreover, we also use CPCA to obtain an estimate of the covariance matrix and use it to construct a minimum variance portfolio (MVP).

Clustering Stock Returns
The data are collected from the Center for Research in Security Prices (CRSP) and include the daily stock returns of 160 companies from 1st January, 2014 to 31st December, 2014, with 252 trading days. The 160 stocks are selected from eight different industries according to Fama and French's 48-industry classification (Fama and French, 1997), namely, Candy and Soda, Tobacco Products, Apparel, Aircraft, Shipbuilding and Railroad Equipment, Petroleum and Natural Gas, Measuring and Control Equipment, and Shipping Containers, with 20 stocks from each industry. The data for the first 126 trading days are treated as the training sample, and the rest are the testing sample. Thus, the training data have the dimensions n = 126 and p = 160. In this example, after removing the common effect from the data, we aim to identify the clusters that consist of companies from the same industry.
This cluster structure of stock returns is also discussed in Fan et al. (2013).  vague cluster structure, but the off-diagonals are clearly non-zero, implying that the stocks from the different industries selected in our study tend to be positively correlated. We apply hierarchical clustering directly to the original data and find that most stocks from different industries are clustered in the same group, as displayed in Figure 6 (b). Only stocks from the 5 th industry stand out as another cluster, because they have a relatively lower correlation with other stocks. If we consider the industries as true clusters, the hierarchical clustering of the original data results in an ARI of only 0.09. This indicates that common components may exist and conceal the cluster structure within the data.
Next, we apply our proposed CPCA to the stock data, and one common component is determined using the ratio-based estimator (9). Figure 5 (b) presents the correlation plot of the data after the common effect is removed. Comparatively, the cluster structure is more apparent. One can interpret the common effect as the market effect and the rest of the components as industry-specific effects. As demonstrated in Figure 6 (c), after the initial step of CPCA, a few clusters of stocks from the 1 st , 2 nd , and 7 th industries are correctly identified, while stocks from other industries (the 3 rd , 4 th , 5 th , 6 th , and 8 th ) cannot be correctly separated from two large clusters. The ARI after the initial step of CPCA is about 0.44 which is a significant improvement to 0.09. Figure 6 (d) depicts the final clusters after the iterative steps of CPCA where most stocks from the aforementioned two large clusters can be correctly separated according to industries with a minor cost of a few miss-classifications in the 1 st and 7 th cluster. The ARI is also further improved from 0.44 to 0.72. The final clusters mainly capture the industry information, but a few stocks are not clustered into their own industries. This is to be expected because some stocks from the same industry are not highly correlated, as shown in Figure 5 (e.g., some stocks from the first, third, sixth and seventh industries), which can be common in reality due to diversification of enterprises.
From the prediction perspective, we compute the MSRE for PCA, CPCA I, and CPCA F as 3.56, 2.23, and 1.67 (in unit of 10 −4 ), respectively. In other words, our CPCA can not only capture the market effect of stock returns as the homogeneity but also include some sub-homogeneities which can be interpreted as industrial effects. All of these findings firmly support that our proposed CPCA is appropriate for analyzing these stock return data.

Constructing Minimum Variance Portfolio
Next, we consider the problem of constructing a global minimum variance portfolio (MVP) from stock returns. The minimum variance portfolio is built by solving the following optimization problem min w w Σw where w is a p-dimensional vector of weights allocated on each stock, 1 p is a p-dimensional vector of ones and Σ is defined as an estimate of the p by p covariance matrix of stock returns. The analytic solution of w is where we assume short sale is allowed and there are no transaction costs.
Constructing the global MVP is another important application of our proposed CPCA due to its ability of correctly identifying both the homogeneity and the sub-homogeneity when estimating the covariance matrices of stock returns. tradings days before day i (from day (i − 110) to day (i − 1)) are considered, and then compute the optimal weight w i by (18). Consequently, the return of the global MVP on day i can be computed using w i and the return of each stock on day i. Similar to Wang et al.
(2021), we compute the standard deviation (STD), the information ratio (IR), which is the ratio of mean return to STD, and the Sharpe ratio (SR), which is the ratio of mean excess return (the portfolio-wide return minus the risk-free return) to STD. As discussed by Engle et al. (2019), Chen et al. (2019) and Wang et al. (2021), the STD should be the primary measure when comparing the MVPs constructed using different estimates of covariance matrices, while the IR and SR are also desirable but less crucial. The numerical results in Table 3 compare the portfolios constructed with different estimates of the covariance matrices. As depicted, CPCA with both the initial and final clusters can achieve a slightly better estimate of the covariance matrix, since the portfolio corresponding to CPCA with has the smallest STD. In addition, the global MVP constructed using CPCA with the final clustering index can also achieve highest IR and SR across all methods, which is also desirable. As a result, CPCA works well when constructing the global MVP since the sub-homogeneity of industrial effects can be well estimated.

Discussion
To conclude, we introduced a novel CPCA to study the homogeneity and sub-homogeneity of high-dimensional data collected from different populations, where the sub-homogeneity refers to a group-specific feature from a particular population. Our numerical simulations confirmed that traditional PCA can only extract the homogeneity from the data, whereas CPCA not only provides a more accurate estimate of the common components, but also identifies the group-specific features, even for the group with small variations. Under a PCR setting, the features extracted using CPCA can significantly outperform those selected by classical PCA in terms of prediction, especially when n is small but p is large. Our real analysis of the stock return data also demonstrated that, when using CPCA, we can capture the industry information after the common component (i.e., market effect) is removed. All of these findings support the use of our proposed CPCA in dimension-reduction problems.
The applications of CPCA are not limited to producing principal components in PCR and revealing the industry structure from the stock return data. CPCA can also be applied to any other data sets in which a group structure exists but hides in the homogeneity.
Further, it can be used to estimate a covariance matrix and its inverse of a large data set that exhibits a group structure. More applications of CPCA will be explored in our future work. Appendices Appendix A Estimation of G 1 and Φ 1 in CPCA Iterative Step (a) In this part, we describe the estimation of the principal components, along with their eigenvectors, combined from each cluster in Iterative Step (a) of Algorithm 1.
First, we perform PCA in cluster j by: where Ψ (j) is an n × r j matrix of principal components for variables in cluster j and Π (j) is a p j × r j matrix in which each column represents an eigenvector of X (j) X (j) . Here, recall that p j denotes the number of variables in cluster j, and r j is the number of principal components in cluster j. Then, we can combine the principal components Ψ (j) from each cluster as Ψ = (Ψ (1) , ..., Ψ (J 0 ) ) and perform a further step of PCA on Ψ to obtain: where we assume the first r c principal components can be used to summarize the common effects among all clusters. Then, G 1 is a n × r matrix of principal components of Ψ, and H is a ( J 0 j=1 r j ) × r c matrix in which each column represents an eigenvector of Ψ Ψ. Lastly, we can find the complement for cluster j as: where Π (j) = (0 (1) , 0 (2) , ..., Π (j) , ..., 0 (J 0 ) ) is a p j ×( J 0 j=1 r j ) matrix in which Π (j) is a p j ×r j matrix, as we described earlier, and 0 (l) is a p j × r l zero matrix for l = 1, 2, ..., J 0 ; l = j.
Once we obtain the complement from each cluster, we can combine them and present the (20) and (21), our estimator of ζ (d) m can be written as: , we obtain the result in this theorem.