Homogeneity and Sub-homogeneity Pursuit: Iterative Complement Clustering PCA

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 sub-homogeneity, 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 sub-homogeneity 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 high-dimensional 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 homogeneity. 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.

Considering the data $\bm{x}_{i}=(x_{1i},\ldots,x_{pi})^{\top}\in\mathbb{R}^{p}$ from $i^{th}$ observation with $p$ variables, the singular value decomposition (SVD) of the data $\bm{x}_{i}$ can be written as:

where ${g}_{ik}$ is defined as $k^{th}$ principal component score and $\bm{\phi}_{k}$ denotes a $p\times{1}$ eigenvector for ${g}_{ik}$. Traditional PCA summarizes the data using the first $k$ principal components, and it treats $\bm{u}_{i}$ as noise because ${g}_{ik}$ for $k=m+1,\ldots,\min{(n,p)}$ has smaller variance. However, under certain conditions, some of the information contained in $\bm{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:

The first line in (2) measures the homogeneity among all variables from all groups, where ${g}_{ik},k=1,\ldots,r_{c}$ is $k^{th}$ principal component, which we call the common component, and $\bm{\phi}_{k}$ is its corresponding eigenvector. The first part of $\bm{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}^{(j)}_{ih},h=1,\ldots,r_{j}$ measures the within-cluster principal components for cluster $j$. The within-cluster eigenvector with dimension $p\times 1$ has the form of $\bm{\gamma}^{(j)}_{h}=(\bm{0}^{(1)\top},\ldots,\bm{\eta}^{(j)\top}_{h},\ldots,\bm{0}^{(J)\top})^{\top}$, where $\bm{\eta}^{(j)}_{h}$ defines a $p_{j}\times 1$ vector and $p_{j}$ is the number of variables in cluster $j$ so that $\sum_{j=1}^{J}p_{j}=p$. That is, the values of $\bm{\gamma}^{(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., $\bm{u}_{i}$ has a block-diagonal covariance structure). In the second part of $\bm{u}_{i}$, $\bm{\epsilon}^{(j)}_{i}$ is simply a $p_{j}$-dimensional error that has variance ${\sigma^{(j)}}^{2}$, and $\bm{I}^{(j)}$ is a diagonal matrix, but the diagonals for variables that do not belong to cluster $j$ are zero.

We further define $\bm{g}_{i}=({g}_{i1},\ldots,{g}_{ir_{c}})^{\top}$ and $\bm{f}^{(j)}_{i}=({f}^{(j)}_{i1},\ldots,{f}^{(j)}_{ir_{j}})^{\top},j=1,\ldots,J$ as the vector forms of the common components and cluster-specific components, and the eigenvector matrices as $\bm{\Phi}_{(p\times r_{c})}=(\bm{\phi}_{1},\ldots,\bm{\phi}_{r_{c}})$ and $\bm{\Gamma}^{(j)}_{(p\times r_{j})}=(\bm{\gamma}^{(j)}_{1},\ldots,\bm{\gamma}^{(j)}_{r_{j}}),j=1,\ldots,J$, respectively. Without loss of generality, we assume $E(\bm{g}_{i})=E(\bm{f}^{(j)}_{i})=E(\bm{\epsilon}^{(j)}_{i})=\bm{0}$, and $\bm{g}_{i},\bm{f}^{(j)}_{i},\bm{\epsilon}^{(j)}_{i},j=1,\ldots,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 $\bm{X}_{(n\times p)}=(\bm{x}_{1},\ldots,\bm{x}_{n})^{\top}$, the common components and cluster-specific components in a matrix form $\bm{G}_{(n\times{r_{c}})}=(\bm{g}_{1},\ldots,\bm{g}_{n})^{\top}$ and $\bm{F}^{(j)}_{(n\times{r_{j}})}=(\bm{f}^{(j)}_{1},\ldots,\bm{f}^{(j)}_{n})^{\top}$ , 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 $\bm{g}_{i}$ and the cluster-specific components $\bm{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 $\bm{\Sigma}$ and the precision matrix $\bm{\Sigma}^{-1}$ of $\bm{X}$. The use of both $\bm{g}_{i}$ and $\bm{f}^{(j)}_{i}$ captures both the low rank and the block-diagonal representation of $\bm{\Sigma}$, which results in a more efficient estimation compared with using $\bm{g}_{i}$ only. Another important application is PCR. In some situations, it is a cluster-specific component $\bm{f}^{(j)}_{i}$ that contributes most to determining the response in PCR rather than the common component $\bm{g}_{i}$. It is important to identify each cluster-specific component $\bm{f}^{(j)}_{i}$ 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.

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, $\bm{f}^{(j)}_{i}=\bm{0},j=1,\ldots,J$ and ${\sigma^{(j)2}}\equiv{\sigma}^{2}$, representation (2) simply reduces to the well-known PCA

  $\displaystyle\bm{x}_{i}=\sum_{k=1}^{r_{c}}{g}_{ik}\bm{\phi}_{k}+\bm{\epsilon}_{i},\quad i=1,\ldots,n,$

  (6)

  where $\bm{g}_{ik}$ can be found by the principal component with $k^{th}$ largest eigenvalue and $\bm{\phi}_{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., $\bm{g}_{i}=\bm{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

  $\displaystyle\bm{x}_{i}=\sum_{j=1}^{J}\sum_{h=1}^{r_{j}}{f}^{(j)}_{ih}\bm{\gamma}^{(j)}_{h}+\sum_{j=1}^{J}\bm{I}^{(j)}\bm{\epsilon}^{(j)}_{i},\quad i=1,\ldots,n,\quad j=1,\ldots,J.$

  (7)

  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 $\bm{f}^{(j)}_{i}$. In many situations, this cluster structure is masked by a dominant common effect. Ignoring this common effect will result in a non-identifiable 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

  $\displaystyle\bm{x}_{i}=\sum_{k=1}^{r_{c}}{g}_{ik}\bm{\phi}_{k}+\bm{u}_{i},\quad i=1,\ldots,n,$

  (8)

  where they assume that the covariance of $\bm{u}_{i}$, $\Sigma_{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 $\bm{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 $\bm{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.

Correctly estimating $\bm{g}_{i}$ and $\bm{f}^{(j)}_{i}$ poses many challenges. First, PCA, which is widely used to estimate the common component $\bm{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 $\left(\bm{x}_{i}-\sum_{k=1}^{r_{c}}{g}_{ik}\bm{\phi}_{k}\right)$, 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.

Therefore, CPCA produces these required outputs: the final clusters of $p$ variables $C^{(j)}_{f},j=1,\ldots,J_{f}$, the final estimate of the common components $\bm{G}_{f}$ and the cluster-specific components $\bm{F}_{f}^{(j)},j=1,\ldots,J_{f}$, along with their eigenvectors. Some details and discussions of CPCA are provided in the following remarks.

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\in\{1,2,\ldots,p\}$, such that:

We now consider another cluster $d$ with components $g_{ik},k=1,\ldots,r_{c}$ and $f_{ih}^{(d)},h=1,\ldots,r_{d}$. When applying our proposed LOO-PCR clustering method, we regress $u_{mi}^{(l)}$ on $f_{ih}^{(d)},h=1,2,\ldots,r_{d}$ and measure its goodness of fit to determine whether the random variable $X_{mi}$ belongs to cluster $d$, such that:

where $\beta_{1},\ldots,\beta_{r_{d}}$ are coefficients and $\zeta_{mi}^{(d)}$ is the error. Naturally, we expect that features (e.g., principal components) from cluster $d$ cannot explain $u_{mi}^{(l)}$ sufficiently. Based on (10) and (11), intuitively, a large discrepancy between principal components $f_{ih}^{(d)}$ from cluster $d$ and principal components $f_{ih}^{(l)}$ from cluster $l$ will result in a large residual $\zeta_{mi}^{(d)}$ in (11).

The following theorem will show the property of $\widehat{\zeta}_{mi}^{(d)}$, which is an estimator of the error $\zeta_{mi}^{(d)}$ from our proposed clustering approach.

A brief proof of Theorem 1 is provided in Appendix B.

In this section, we conduct various simulation studies to investigate the performance of our proposed CPCA compared with traditional PCA under different simulation settings.

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).

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.

SUPPLEMENTARY MATERIAL

R-codes for CPCA:

R-codes including R functions to perform the iterative complement-clustering principal component analysis (CPCA) described in the article. Also included are the R-codes for simulation studies and real analysis of the stock return data shown in the article. (.R files)

Data of stock returns:

Data set of daily returns of stocks from different industries in the US (2014). (.rds file)