advertisement

Sparse Principal Component Analysis Hui Zou∗, Trevor Hastie†, Robert Tibshirani‡ April 26, 2004 Abstract Principal component analysis (PCA) is widely used in data processing and dimensionality reduction. However, PCA suffers from the fact that each principal component is a linear combination of all the original variables, thus it is often difficult to interpret the results. We introduce a new method called sparse principal component analysis (SPCA) using the lasso (elastic net) to produce modified principal components with sparse loadings. We show that PCA can be formulated as a regression-type optimization problem, then sparse loadings are obtained by imposing the lasso (elastic net) constraint on the regression coefficients. Efficient algorithms are proposed to realize SPCA for both regular multivariate data and gene expression arrays. We also give a new formula to compute the total variance of modified principal components. As illustrations, SPCA is applied to real and simulated data, and the results are encouraging. Keywords: multivariate analysis, gene expression arrays, elastic net, lasso, singular value decomposition, thresholding ∗ Hui Zou is a Ph.D student in the Department of Statistics at Stanford University, Stanford, CA 94305. Email: hzou@stat.stanford.edu. † Trevor Hastie is Professor, Department of Statistics and Department of Health Research & Policy, Stanford University, Stanford, CA 94305. Email: hastie@stat.stanford.edu. ‡ Robert Tibshirani is Professor, Department of Health Research & Policy and Department of Statistics, Stanford University, Stanford, CA 94305. Email: tibs@stat.stanford.edu. 1 1 Introduction Principal component analysis (PCA) (Jolliffe 1986) is a popular data processing and dimension reduction technique . As an un-supervised learning method, PCA has numerous applications such as handwritten zip code classification (Hastie et al. 2001) and human face recognition (Hancock et al. 1996). Recently PCA has been used in gene expression data analysis (Misra et al. 2002). Hastie et al. (2000) propose the so-called Gene Shaving techniques using PCA to cluster high variable and coherent genes in microarray data. PCA seeks the linear combinations of the original variables such that the derived variables capture maximal variance. PCA can be done via the singular value decomposition (SVD) of the data matrix. In detail, let the data X be a n × p matrix, where n and p are the number of observations and the number of variables, respectively. Without loss of generality, assume the column means of X are all 0. Suppose we have the SVD of X as X = UDVT (1) where T means transpose. U are the principal components (PCs) of unit length, and the columns of V are the corresponding loadings of the principal components. The variance of the ith PC is 2 . In gene expression data the PCs U are called the eigen-arrays and V are the eigen-genes Di,i (Alter et al. 2000). Usually the first q (q " p) PCs are chosen to represent the data, thus a great dimensionality reduction is achieved. The success of PCA is due to the following two important optimal properties: 1. principal components sequentially capture the maximum variability among X, thus guaranteeing minimal information loss; 2. principal components are uncorrelated, so we can talk about one principal component without 2 referring to others. However, PCA also has an obvious drawback, i.e., each PC is a linear combination of all p variables and the loadings are typically nonzero. This makes it often difficult to interpret the derived PCs. Rotation techniques are commonly used to help practitioners to interpret principal components (Jolliffe 1995). Vines (2000) considered simple principal components by restricting the loadings to take values from a small set of allowable integers such as 0, 1 and -1. We feel it is desirable not only to achieve the dimensionality reduction but also to reduce the size of explicitly used variables. An ad hoc way is to artificially set the loadings with absolute values smaller than a threshold to zero. This informal thresholding approach is frequently used in practice but can be potentially misleading in various respects (Cadima & Jolliffe 1995). McCabe (1984) presented an alternative to PCA which found a subset of principal variables. Jolliffe & Uddin (2003) introduced SCoTLASS to get modified principal components with possible zero loadings. Recall the same interpretation issue arising in multiple linear regression, where the response is predicted by a linear combination of the predictors. Interpretable models are obtained via variable selection. The lasso (Tibshirani 1996) is a promising variable selection technique, simultaneously producing accurate and sparse models. Zou & Hastie (2003) propose the elastic net, a generalization of the lasso, to further improve upon the lasso. In this paper we introduce a new approach to get modified PCs with sparse loadings, which we call sparse principal component analysis (SPCA). SPCA is built on the fact that PCA can be written as a regression-type optimization problem, thus the lasso (elastic net) can be directly integrated into the regression criterion such that the resulting modified PCA produces sparse loadings. In the next section we briefly review the lasso and the elastic net. The method details of SPCA are presented in Section 3. We first discuss a direct sparse approximation approach via the elastic net, which is a useful exploratory tool. We then show that finding the loadings of principal 3 components can be reformulated as estimating coefficients in a regression-type optimization problem. Thus by imposing the lasso (elastic net) constraint on the coefficients, we derive the modified principal components with sparse loadings. An efficient algorithm is proposed to realize SPCA. We also give a new formula, which justifies the correlation effects, to calculate the total variance of modified principal components. In Section 4 we consider a special case of the SPCA algorithm to efficiently handle gene expression arrays. The proposed methodology is illustrated by using real data and simulation examples in Section 5. Discussions are in Section 6. The paper ends up with an appendix summarizing technical details. 2 The Lasso and The Elastic Net Consider the linear regression model. Suppose the data set has n observations with p predictors. Let Y = (y1 , . . . , yn )T be the response and Xj = (x1j , . . . , xnj )T , i = 1, . . . , p are the predictors. After a location transformation we can assume all Xj and Y are centered. The lasso is a penalized least squares method, imposing a constraint on the L 1 norm of the regression coefficients. Thus the lasso estimates β̂lasso are obtained by minimizing the lasso criterion β̂lasso ! !2 ! ! p p " " ! ! ! ! = arg min !Y − Xj βj ! + λ |βj | , β ! ! j=1 j=1 where λ is a non-negative value. (2) The lasso was originally solved by quadratic programming (Tibshirani 1996). Efron et al. (2004) proved that the lasso estimates as a function of λ are piecewise linear, and proposed an algorithm called LARS to efficiently solve the whole lasso solution path in the same order of computations as a single least squares fit. The lasso continuously shrinks the coefficients toward zero, thus gaining its prediction accuracy via the bias variance trade-off. Moreover, due to the nature of the L1 penalty, some coefficients 4 will be shrunk to exact zero if λ1 is large enough. Therefore the lasso simultaneously produces an accurate and sparse model, which makes it a favorable variable selection method. However, the lasso has several limitations as pointed out in Zou & Hastie (2003). The most relevant one to this work is that the number of selected variables by the lasso is limited by the number of observations. For example, if applied to the microarray data where there are thousands of predictors (genes) (p > 1000) with less than 100 samples (n < 100), the lasso can only select at most n genes, which is clearly unsatisfactory. The elastic net (Zou & Hastie 2003) generalizes the lasso to overcome its drawbacks, while enjoying the similar optimal properties. For any non-negative λ1 and λ2 , the elastic net estimates β̂en are given as follows β̂en ! !2 ! ! p p p " " " ! ! 2 ! ! = (1 + λ2 ) arg min !Y − Xj βj ! + λ2 |βj | + λ1 |βj | . β ! ! j=1 j=1 j=1 (3) Hence the elastic net penalty is a convex combination of ridge penalty and the lasso penalty . Obviously, the lasso is a special case of the elastic net with λ2 = 0. Given a fixed λ2 , the LARSEN algorithm (Zou & Hastie 2003) efficiently solves the elastic net problem for all λ 1 with the computation cost as a single least squares fit. When p > n, we choose some λ 2 > 0. Then the elastic net can potentially include all variables in the fitted model, so the limitation of the lasso is removed. An additional benefit offered by the elastic net is its grouping effect, that is, the elastic net tends to select a group of highly correlated variables once one variable among them is selected. In contrast, the lasso tends to select only one out of the grouped variables and does not care which one is in the final model. Zou & Hastie (2003) compare the elastic net with the lasso and discuss the application of the elastic net as a gene selection method in microarray analysis. 5 3 Motivation and Method Details In both lasso and elastic net, the sparse coefficients are a direct consequence of the L 1 penalty, not depending on the squared error loss function. Jolliffe & Uddin (2003) proposed SCoTLASS by directly putting the L1 constraint in PCA to get sparse loadings. SCoTLASS successively maximizes the variance aTk (XT X)ak (4) subject to aTk ak = 1 and (for k ≥ 2) aTh ak = 0, h < k; (5) and the extra constraints p " j=1 |ak,j | ≤ t (6) for some tuning parameter t. Although sufficiently small t yields some exact zero loadings, SCoTLASS seems to lack of a guidance to choose an appropriate t value. One might try several t values, but the high computational cost of SCoTLASS makes it an impractical solution. The high computational cost is due to the fact that SCoTLASS is not a convex optimization problem. Moreover, the examples in Jolliffe & Uddin (2003) show that the obtained loadings by SCoTLASS are not sparse enough when requiring a high percentage of explained variance. We consider a different approach to modify PCA, which can more directly make good use of the lasso. In light of the success of the lasso (elastic net) in regression, we state our strategy We seek a regression optimization framework in which PCA is done exactly. In addition, the regression framework should allow a direct modification by using the lasso (elastic net) penalty such that the derived loadings are sparse. 6 3.1 Direct sparse approximations We first discuss a simple regression approach to PCA. Observe that each PC is a linear combination of the p variables, thus its loadings can be recovered by regressing the PC on the p variables. Theorem 1 ∀i, denote Yi = Ui Di . Yi is the i-th principal component. ∀ λ > 0, suppose β̂ridge is the ridge estimates given by β̂ridge = arg min |Yi − Xβ|2 + λ |β|2 . β Let v̂ = β̂ridge |β̂ridge | (7) , then v̂ = Vi . The theme of this simple theorem is to show the connection between PCA and a regression method is possible. Regressing PCs on variables was discussed in Cadima & Jolliffe (1995), where they focused on approximating PCs by a subset of k variables. We extend it to a more general ridge regression in order to handle all kinds of data, especially the gene expression data. Obviously when n > p and X is a full rank matrix, the theorem does not require a positive λ. Note that if p > n and λ = 0, ordinary multiple regression has no unique solution that is exactly V i . The same story happens when n > p and X is not a full rank matrix. However, PCA always gives a unique solution in all situations. As shown in theorem 1, this discrepancy is eliminated by the positive ridge penalty (λ |β|2 ). Note that after normalization the coefficients are independent of λ, therefore the ridge penalty is not used to penalize the regression coefficients but to ensure the reconstruction of principal components. Hence we keep the ridge penalty term throughout this paper. Now let us add the L1 penalty to (7) and consider the following optimization problem β̂ = arg min |Yi − Xβ|2 + λ |β|2 + λ1 |β|1 . β 7 (8) We call V̂i = β̂ |β̂ | an approximation to Vi , and XV̂i the ith approximated principal component. (8) is called naive elastic net (Zou & Hastie 2003) which differs from the elastic net by a scaling factor (1 + λ). Since we are using the normalized fitted coefficients, the scaling factor does not affect V̂i . Clearly, large enough λ1 gives a sparse β̂, hence a sparse V̂i . Given a fixed λ, (8) is efficiently solved for all λ1 by using the LARS-EN algorithm (Zou & Hastie 2003). Thus we can flexibly choose a sparse approximation to the ith principal component. 3.2 Sparse principal components based on SPCA criterion Theorem 1 depends on the results of PCA, so it is not a genuine alternative. However, it can be used in a two-stage exploratory analysis: first perform PCA, then use (8) to find suitable sparse approximations. We now present a “self-contained” regression-type criterion to derive PCs. We first consider the leading principal component. Theorem 2 Let Xi denote the ith row vector of the matrix X. For any λ > 0, let n " ! ! !Xi − αβ T Xi !2 + λ |β|2 (α̂, β̂) = arg min α,β (9) i=1 subject to |α|2 = 1. Then β̂ ∝ V1 . The next theorem extends theorem 2 to derive the whole sequence of PCs. Theorem 3 Suppose we are considering the first k principal components. Let α and β be p × k matrices. Xi denote the i-th row vector of the matrix X. For any λ > 0, let (α̂, β̂) = arg min α,β n k " " ! ! !Xi − αβ T Xi !2 + λ |βj |2 i=1 j=1 8 (10) subject to α T α = Ik . Then β̂i ∝ Vi for i = 1, 2, . . . , k. Theorem 3 effectively transforms the PCA problem to a regression-type problem. The critical element is the object function !2 !2 #n !! #n !! T T ! ! i=1 Xi − αβ Xi . If we restrict β = α, then i=1 Xi − αβ Xi = !2 #n !! Xi − ααT Xi ! , whose minimizer under the orthonormal constraint on α is exactly the first k i=1 loading vectors of ordinary PCA. This is actually an alternative derivation of PCA other than the maximizing variance approach, e.g. Hastie et al. (2001). Theorem 3 shows that we can still have exact PCA while relaxing the restriction β = α and adding the ridge penalty term. As can be seen later, these generalizations enable us to flexibly modify PCA. To obtain sparse loadings, we add the lasso penalty into the criterion (10) and consider the following optimization problem (α̂, β̂) = arg min α,β n k k " " " ! ! !Xi − αβ T Xi !2 + λ |βj |2 + λ1,j |βj | i=1 j=1 1 j=1 (11) subject to αT α = Ik . Whereas the same λ is used for all k components, different λ1,j s are allowed for penalizing the loadings of different principal components. Again, if p > n, a positive λ is required in order to get exact PCA when the sparsity constraint (the lasso penalty) vanishes (λ1,j = 0). (11) is called the SPCA criterion hereafter. 9 3.3 Numerical solution We propose an alternatively minimization algorithm to minimize the SPCA criterion. From the proof of theorem 3 (see appendix for details) we get !2 #n !! # Xi − αβ T Xi ! + λ k j=1 |βj | i=1 = TrXT X + #k j=1 $ 2 + #k j=1 λ1,j |βj |1 % βjT (XT X + λ)βj − 2αjT XT Xβj + λ1,j |βj |1 . (12) Hence if given α, it amounts to solve k independent elastic net problems to get β̂j for j = 1, 2, . . . , k. On the other hand, we also have (details in appendix) !2 #n !! # Xi − αβ T Xi ! + λ k j=1 |βj | i=1 2 + #k j=1 λ1,j = TrXT X − 2TrαT XT Xβ + Trβ T (XT X + λ)β + |βj |1 #k j=1 λ1,j |βj |1 . (13) Thus if β is fixed, we should maximize TrαT (XT X)β subject to αT α = Ik , whose solution is given by the following theorem. Theorem 4 Let α and β be m × k matrices and β has rank k. Consider the constrained maximization problem & ' α̂ = arg max Tr αT β subject to α α T α = Ik . Suppose the SVD of β is β = U DV T , then α̂ = U V T . Here are the steps of our numerical algorithm to derive the first k sparse PCs. General SPCA Algorithm 1. Let α start at V[, 1 : k], the loadings of first k ordinary principal components. 10 (14) 2. Given fixed α, solve the following naive elastic net problem for j = 1, 2, . . . , k βj = arg min β ∗T (XT X + λ)β ∗ − 2αjT XT Xβ ∗ + λ1,j |β ∗ |1 . ∗ β (15) 3. For each fixed β, do the SVD of XT Xβ = U DV T , then update α = U V T . 4. Repeat steps 2-3, until β converges. 5. Normalization: V̂j = βj |βj | , j = 1, . . . , k. Some remarks: 1. Empirical evidence indicates that the outputs of the above algorithm vary slowly as λ changes. For n > p data, the default choice of λ can be zero. Practically λ is a small positive number to overcome potential collinearity problems of X. Section 4 discusses the default choice of λ for the data with thousands of variables, such as gene expression arrays. 2. In principle, we can try several combinations of {λ1,j } to figure out a good choice of the tunning parameters, since the above algorithm converges quite fast. There is a shortcut provided by the direct sparse approximation (8). The LARS-EN algorithm efficiently deliver a whole sequence of sparse approximations for each PC and the corresponding values of λ 1,j . Hence we can pick a λ1,j which gives a good compromise between variance and sparsity. In this selection, variance has a higher priority than sparsity, thus we tend to be conservative in pursuing sparsity. 3. Both PCA and SPCA depend on X only through XT X. Note that XT X n is actually the sample covariance matrix of variables (Xi ). Therefore if Σ, the covariance matrix of (Xi ), is known, we can replace XT X with Σ and have a population version of PCA or SPCA. If X is 11 standardized beforehand, then PCA or SPCA uses the (sample) correlation matrix, which is preferred when the scales of the variables are different. 3.4 Adjusted total variance The ordinary principal components are uncorrelated and their loadings are orthogonal. Let Σ̂ = XT X, then VT V = Ik and VT Σ̂V is diagonal. It is easy to check that only the loadings of ordinary principal components can satisfy both conditions. In Jolliffe & Uddin (2003) the loadings were forced to be orthogonal, so the uncorrelated property was sacrificed. SPCA does not explicitly impose the uncorrelated components condition too. Let Û be the modified PCs. Usually the total variance explained by Û is calculated by trace(ÛT Û). This is unquestionable when Û are uncorrelated. However, if they are correlated, the computed total variance is too optimistic. Here we propose a new formula to compute the total variance explained by Û, which takes into account the correlations among Û. Suppose (Ûi , i = 1, 2, . . . , k) are the first k modified PCs by any method. Denote Ûj·1,...,j−1 the reminder of Ûj after adjusting the effects of Û1 , . . . , Ûj−1 , that is Ûj·1,...,j−1 = Ûj − H1,...,j−1 Ûj , (16) where H1,...,j−1 is the projection matrix on Ûi i = 1, 2, . . . , j − 1. Then the adjusted variance !2 !2 ! #k !! ! ! ! of Ûj is !Ûj·1,...,j−1 ! , and the total explained variance is given by j=1 !Ûj·1,...,j−1 ! . When the modified PCs Û are uncorrelated, then the new formula agrees with trace(ÛT Û). Note that the above computations depend on the order of Ûi . However, since we have a natural order in PCA, ordering is not an issue here. Using the QR decomposition, we can easily compute the adjusted variance. Suppose Û = QR, 12 where Q is orthonormal and R is upper triangular. Then it is straightforward to see that ! !2 ! ! !Ûj·1,...,j−1 ! = R2j,j . Hence the explained total variance is equal to 3.5 (17) #k 2 j=1 Rj,j . Computation complexity PCA is computationally efficient for both n > p or p ( n data. We separately discuss the computational cost of the general SPCA algorithm for n > p and p ( n. 1. n > p. Traditional multivariate data fit in this category. Note that although the SPCA criterion is defined using X, it only depends on X via XT X. A trick is to first compute the p × p matrix Σ̂ = XT X once for all, which requires np2 operations. Then the same Σ̂ is used at each step within the loop. Computing XT Xβ costs p2 k and the SVD of XT Xβ is of order O(pk 2 ). Each elastic net solution requires at most O(p3 ) operations. Since k ≤ p, the total computation cost is at most np2 + mO(p3 ), where m is the number of iterations before convergence. Therefore the SPCA algorithm is able to efficiently handle data with huge n, as long as p is small (say p < 100). 2. p ( n. Gene expression arrays are typical examples of this p ( n category. The trick of Σ̂ is no longer applicable, because Σ̂ is a huge matrix (p × p) in this case. The most consuming step is solving each elastic net, whose cost is of order O(pJ 2 ) for a positive finite λ, where J is the number of nonzero coefficients. Generally speaking the total cost is of order mO(pJ 2 k), which is expensive for a large J. Fortunately, as shown in Section 4, there exits a special SPCA algorithm for efficiently dealing with p ( n data. 13 4 SPCA for p ! n and Gene Expression Arrays Gene expression arrays are a new type of data where the number of variables (genes) are much bigger than the number of samples. Our general SPCA algorithm still fits this situation using a positive λ. However the computation cost is expensive when requiring a large number of nonzero loadings. It is desirable to simplify the general SPCA algorithm to boost the computation. Observe that theorem 3 is valid for all λ > 0, so in principle we can use any positive λ. It turns out that a thrifty solution emerges if λ → ∞. Precisely, we have the following theorem. Theorem 5 Let V̂i (λ) = β̂i |β̂i | be the loadings derived from criterion (11). Define (α̂∗ , β̂ ∗ ) as the solution of the optimization problem ∗ ∗ (α̂ , β̂ ) = arg min −2Trα X Xβ + T α,β subject to When λ → ∞, V̂i (λ) → T α T α = Ik . k " j=1 βj2 + k " j=1 λ1,j |βj |1 (18) β̂i∗ . |β̂i∗ | By the same statements in Section 3.3, criterion (18) is solved by the following algorithm, which is a special case of the general SPCA algorithm with λ = ∞. Gene Expression Arrays SPCA Algorithm Replacing step 2 in the general SPCA algorithm with Step 2∗ : Given fixed α, for j = 1, 2, . . . , k ( ) ! T T ! λ1,j ! ! βj = αj X X − Sign(αjT XT X). 2 + (19) The operation in (19) is called soft-thresholding. Figure 1 gives an illustration of how the soft-thresholding rule operates. Recently soft-thresholding has become increasingly popular in 14 y Δ x (0,0) Figure 1: An illstration of soft-thresholding rule y = (|x| − ∆)+ Sign(x) with ∆ = 1. the literature. For example, nearest shrunken centroids (Tibshirani et al. 2002) adopts the softthresholding rule to simultaneously classify samples and select important genes in microarrays. 5 Examples 5.1 Pitprops data The pitprops data first introduced in Jeffers (1967) has 180 observations and 13 measured variables. It is the classic example showing the difficulty of interpreting principal components. Jeffers (1967) tried to interpret the first 6 PCs. Jolliffe & Uddin (2003) used their SCoTLASS to find the modified PCs. Table 1 presents the results of PCA, while Table 2 presents the modified PCs loadings by SCoTLASS and the adjusted variance computed using (17). As a demonstration, we also considered the first 6 principal components. Since this is a usual 15 n ( p data set, we set λ = 0. λ1 = (0.06, 0.16, 0.1, 0.5, 0.5, 0.5) were chosen according to Figure 2 such that each sparse approximation explained almost the same amount of variance as the ordinary PC did. Table 3 shows the obtained sparse loadings and the corresponding adjusted variance. Compared with the modified PCs by SCoTLASS, PCs by SPCA account for nearly the same amount of variance (75.8% vs. 78.2%) but with a much sparser loading structure. The important variables associated with the 6 PCs do not overlap, which further makes the interpretations easier and clearer. It is interesting to note that in Table 3 even though the variance does not strictly monotonously decrease, the adjusted variance follows the right order. However, Table 2 shows this is not true in SCoTLASS. It is also worthy to mention that the whole computation of SPCA was done in seconds in R, while the implementation of SCoTLASS for each t was expensive (Jolliffe & Uddin 2003). Optimizing SCoTLASS over several values of t is even a more difficult computational challenge. Although the informal thresholding method, which is referred to as simple thresholding henceforth, has various drawbacks, it may serve as the benchmark for testing sparse PCs methods. An variant of simple thresholding is soft-thresholding. We found that used in PCA, soft-thresholding performs very similarly to simple thresholding. Thus we omitted the results of soft-thresholding in this paper. Both SCoTLASS and SPCA were compared with simple thresholding. Table 4 presents the loadings and the corresponding explained variance by simple thresholding. To make fair comparisons, we let the numbers of nonzero loadings by simple thresholding match the results of SCoTLASS and SPCA. In terms of variance, it seems that simple thresholding is better than SCoTLASS and worse than SPCA. Moreover, the variables with non-zero loadings by SPCA are very different to that chosen by simple thresholding for the first three PCs; while SCoTLASS seems to create a similar sparseness pattern as simple thresholding does, especially in the leading PC. 16 PC 2 0.00 0.00 0.05 0.05 0.10 PEV 0.15 0.10 PEV 0.20 0.25 0.15 0.30 PC 1 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.5 1.0 1.5 λ1 λ1 PC 3 PC 4 2.0 2.5 PEV 0.04 0.00 0.00 0.02 0.05 PEV 0.10 0.06 0.08 0.15 0.0 0.5 1.0 1.5 0.0 1.0 PC 5 PC 6 1.5 0.06 0.04 0.00 0.00 0.02 PEV 0.04 0.06 0.08 λ1 0.02 PEV 0.5 λ1 0.08 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 λ1 0.2 0.4 0.6 0.8 1.0 λ1 Figure 2: Pitprops data: The sequences of sparse approximations to the first 6 principal components. Plots show the percentage of explained variance (PEV) as a function of λ1 . 17 Table 1: Pitprops data: loadings of the Variable PC1 PC2 topdiam -0.404 0.218 length -0.406 0.186 moist -0.124 0.541 testsg -0.173 0.456 ovensg -0.057 -0.170 ringtop -0.284 -0.014 ringbut -0.400 -0.190 bowmax -0.294 -0.189 bowdist -0.357 0.017 whorls -0.379 -0.248 clear 0.011 0.205 knots 0.115 0.343 diaknot 0.113 0.309 Variance (%) 32.4 18.3 Cumulative Variance (%) 32.4 50.7 Table 2: Pitprops data: loadings of the t = 1.75 Variable PC1 topdiam 0.546 length 0.568 moist 0.000 testsg 0.000 ovensg 0.000 ringtop 0.000 ringbut 0.279 bowmax 0.132 bowdist 0.376 whorls 0.376 clear 0.000 knots 0.000 diaknot 0.000 Number of nonzero loadings 6 Variance (%) 27.2 Adjusted Variance (%) 27.2 Cumulative Adjusted Variance (%) 27.2 18 first 6 principal components PC3 PC4 PC5 PC6 -0.207 0.091 -0.083 0.120 -0.235 0.103 -0.113 0.163 0.141 -0.078 0.350 -0.276 0.352 -0.055 0.356 -0.054 0.481 -0.049 0.176 0.626 0.475 0.063 -0.316 0.052 0.253 0.065 -0.215 0.003 -0.243 -0.286 0.185 -0.055 -0.208 -0.097 -0.106 0.034 -0.119 0.205 0.156 -0.173 -0.070 -0.804 -0.343 0.175 0.092 0.301 -0.600 -0.170 -0.326 0.303 0.080 0.626 14.4 8.5 7.0 6.3 65.1 73.6 80.6 86.9 first 6 modified PCs by SCoTLASS PC2 0.047 0.000 0.641 0.641 0.000 0.356 0.000 -0.007 0.000 -0.065 0.000 0.206 0.000 7 16.4 15.3 42.5 PC3 -0.087 -0.076 -0.187 0.000 0.457 0.348 0.325 0.000 0.000 0.000 0.000 0.000 -0.718 7 14.8 14.4 56.9 PC4 0.066 0.117 -0.127 -0.139 0.000 0.000 0.000 -0.589 0.000 -0.067 0.000 0.771 0.013 8 9.4 7.1 64.0 PC5 -0.046 -0.081 0.009 0.000 -0.614 0.000 0.000 0.000 0.000 0.189 -0.659 0.040 -0.379 8 7.1 6.7 70.7 PC6 0.000 0.000 0.017 0.000 -0.562 -0.045 0.000 0.000 0.065 -0.065 0.725 0.003 -0.384 8 7.9 7.5 78.2 Table 3: Pitprops data: loadings of the Variable PC1 topdiam -0.477 length -0.476 moist 0.000 testsg 0.000 ovensg 0.177 ringtop 0.000 ringbut -0.250 bowmax -0.344 bowdist -0.416 whorls -0.400 clear 0.000 knots 0.000 diaknot 0.000 Number of nonzero loadings 7 Variance (%) 28.0 Adjusted Variance (%) 28.0 Cumulative Adjusted Variance (%) 28.0 5.2 first 6 sparse PCs by SPCA PC2 PC3 PC4 PC5 0.000 0.000 0 0 0.000 0.000 0 0 0.785 0.000 0 0 0.620 0.000 0 0 0.000 0.640 0 0 0.000 0.589 0 0 0.000 0.492 0 0 -0.021 0.000 0 0 0.000 0.000 0 0 0.000 0.000 0 0 0.000 0.000 -1 0 0.013 0.000 0 -1 0.000 -0.015 0 0 4 4 1 1 14.4 15.0 7.7 7.7 14.0 13.3 7.4 6.8 42.0 55.3 62.7 69.5 A simulation example We first created three hidden factors V1 ∼ N (0, 290), V3 = −0.3V1 + 0.925V2 + $, V1 , V2 and $ V2 ∼ N (0, 300) $ ∼ N (0, 1) are independent. Then 10 observed variables were generated as the follows Xi = V1 + $1i , $1i ∼ N (0, 1), i = 1, 2, 3, 4, Xi = V2 + $2i , $2i ∼ N (0, 1), i = 5, 6, 7, 8, Xi = V3 + $3i , $3i ∼ N (0, 1), i = 9, 10, {$ji } are independent, j = 1, 2, 3 19 i = 1, · · · , 10. PC6 0 0 0 0 0 0 0 0 0 0 0 0 1 1 7.7 6.2 75.8 Table 4: Pitprops data: loadings of Variable topdiam length moist testsg ovensg ringtop ringbut bowmax bowdist whorls clear knots diaknot Number of nonzero loadings Variance (%) Adjusted Variance (%) Cumulative Adjusted Variance (%) Variable topdiam length moist testsg ovensg ringtop ringbut bowmax bowdist whorls clear knots diaknot Number of nonzero loadings Variance (%) Adjusted Variance (%) Cumulative Adjusted Variance (%) the first PC1 -0.439 -0.441 0.000 0.000 0.000 0.000 -0.435 -0.319 -0.388 -0.412 0.000 0.000 0.000 6 28.9 28.9 28.9 PC1 -0.420 -0.422 0.000 0.000 0.000 -0.296 -0.416 -0.305 -0.370 -0.394 0.000 0.000 0.000 7 30.7 30.7 30.7 20 6 modified PCs PC2 PC3 0.234 0.000 0.000 -0.253 0.582 0.000 0.490 0.379 0.000 0.517 0.000 0.511 0.000 0.272 0.000 -0.261 0.000 0.000 -0.267 0.000 0.221 0.000 0.369 0.000 0.332 -0.350 7 7 16.6 14.2 16.5 14.0 45.4 59.4 PC2 PC3 0.000 0.000 0.000 0.000 0.640 0.000 0.540 0.425 0.000 0.580 0.000 0.573 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.406 0.000 0.365 -0.393 4 4 14.8 13.6 14.7 11.1 45.4 56.5 by simple thresholding PC4 PC5 PC6 0.092 0.000 0.120 0.104 0.000 0.164 0.000 0.361 -0.277 0.000 0.367 0.000 0.000 0.182 0.629 0.000 -0.326 0.000 0.000 -0.222 0.000 -0.288 0.191 0.000 -0.098 0.000 0.000 0.207 0.000 -0.174 -0.812 -0.354 0.176 0.304 -0.620 -0.171 0.306 0.000 0.629 8 8 8 8.6 6.9 6.3 8.5 6.7 6.2 67.9 74.6 80.8 PC4 PC5 PC6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0 1 1 1 1 7.7 7.7 7.7 7.6 5.2 3.6 64.1 68.3 71.9 To avoid the simulation randomness, we used the exact covariance matrix of (X 1 , . . . , X10 ) to perform PCA, SPCA and simple thresholding. In other words, we compared their performances using an infinity amount of data generated from the above model. The variance of the three underlying factors is 290, 300 and 283.8, respectively. The numbers of variables associated with the three factors are 4, 4 and 2. Therefore V 2 and V1 are almost equally important, and they are much more important than V3 . The first two PCs together explain 99.6% of the total variance. These facts suggest that we only need to consider two derived variables with ‘right’ sparse representations. Ideally, the first derived variable should recover the factor V 2 only using (X5 , X6 , X7 , X8 ), and the second derived variable should recover the factor V1 only using (X1 , X2 , X3 , X4 ). In fact, if we sequentially maximize the variance of the first two derived variables under the orthonormal constraint, while restricting the numbers of nonzero loadings to four, then the first derived variable uniformly assigns nonzero loadings on (X5 , X6 , X7 , X8 ); and the second derived variable uniformly assigns nonzero loadings on (X1 , X2 , X3 , X4 ). Both SPCA (λ = 0) and simple thresholding were carried out by using the oracle information that the ideal sparse representations use only four variables. Table 5 summarizes the comparison results. Clearly, SPCA correctly identifies the sets of important variables. As a matter of fact, SPCA delivers the ideal sparse representations of the first two principal components. Mathematically, it is easy to show that if t = 2 is used, SCoTLASS is also able to find the same sparse solution. In this example, both SPCA and SCoTLASS produce the ideal sparse PCs, which may be explained by the fact that both methods explicitly use the lasso penalty. In contrast, simple thresholding wrongly includes X9 , X10 in the most important variables. The explained variance by simple thresholding is also lower than that by SPCA, although the relative difference is small (less than 5%). Due to the high correlation between V2 and V3 , variables X9 , X10 gain loadings which are even higher than that of the true important varaibles (X 5 , X6 , X7 , X8 ). Thus 21 Table 5: Results of the PCA PC1 PC2 X1 0.116 -0.478 X2 0.116 -0.478 X3 0.116 -0.478 X4 0.116 -0.478 X5 -0.395 -0.145 X6 -0.395 -0.145 X7 -0.395 -0.145 X8 -0.395 -0.145 X9 -0.401 0.010 X10 -0.401 0.010 Adjusted Variance (%) 60.0 39.6 simulation example: loadings and variance SPCA (λ = 0) Simple Thresholding PC3 PC1 PC2 PC1 PC2 -0.087 0.0 0.5 0.000 -0.5 -0.087 0.0 0.5 0.000 -0.5 -0.087 0.0 0.5 0.000 -0.5 -0.087 0.0 0.5 0.000 -0.5 0.270 0.5 0.0 0.000 0.0 0.270 0.5 0.0 0.000 0.0 0.270 0.5 0.0 -0.497 0.0 0.270 0.5 0.0 -0.497 0.0 -0.582 0.0 0.0 -0.503 0.0 -0.582 0.0 0.0 -0.503 0.0 0.08 40.9 39.5 38.8 38.6 the truth is disguised by the high correlation. On the other hand, simple thresholding correctly discovers the second factor, because V1 has a low correlation with V3 . 5.3 Ramaswamy data Ramaswamy data (Ramaswamy et al. 2001) has 16063 (p = 16063) genes and 144 (n = 144) samples. Its first principal component explains 46% of the total variance. In a typical microarray data like this, it appears that SCoTLASS cannot be practically useful. We applied SPCA (λ = ∞) to find the sparse leading PC. A sequence of λ1 were used such that the number of nonzero loadings varied in a rather wide range. As displayed in Figure 3, the percentage of explained variance decreases at a slow rate, as the sparsity increase. As few as 2.5% of these 16063 genes can sufficiently construct the leading principal component with little loss of explained variance (from 46% to 40%). Simple thresholding was also applied to this data. It seems that when using the same number of genes, simple thresholding always explains slightly higher variance than SPCA does. Among the same number of selected genes by SPCA and simple thresholding, there are about 2% different genes, and this difference rate is quite consistent. 22 PEV 0.40 0.41 0.42 0.43 0.44 0.45 0.46 Ramaswamy data 0.37 SPCA simple thresholding 200 1600 4000 8000 16063 number of nonzero loadings Figure 3: The sparse leading principal component: percentage of explained variance versus sparsity. Simple thresholding and SPCA have similar performances. However, there still exists consistent difference in the selected genes (the ones with nonzero loadings). 6 Discussion It has been a long standing interest to have a formal approach to derive principal components with sparse loadings. From a practical point of view, a good method to achieve the sparseness goal should (at least) possess the following properties. • Without any sparsity constraint, the method should reduce to PCA. • It should be computationally efficient for both small p and big p data. • It should avoid mis-identifying the important variables. The frequently used simple thresholding is not criterion based. However, this informal ad hoc method seems to have the first two of the good properties listed above. If the explained variance and sparsity are the only concerns, simple thresholding is not such a bad choice, and it is extremely convenient. We have shown that simple thresholding can work pretty well in gene expression 23 arrays. The serious problem with simple thresholding is that it can mis-identify the real important variables. Nevertheless, simple thresholding is regarded as a benchmark for any potentially better method. Using the lasso constraint in PCA, SCoTLASS successfully derives sparse loadings. However, SCoTLASS is not computationally efficient, and it lacks a good rule to pick its tunning parameter. In addition, it is not feasible to apply SCoTLASS to gene expression arrays, while in which PCA is a quite popular tool. In this work we have developed SPCA using the SPCA criterion. The new SPCA criterion gives exact PCA results when its sparsity (lasso) penalty term vanishes. SPCA allows a quite flexible control on the sparse structure of the resulting loadings. Unified efficient algorithms have been proposed to realize SPCA for both regular multivariate data and gene expression arrays. As a principled procedure, SPCA enjoys advantages in several aspects, including computational efficiency, high explained variance and ability of identifying important variables. 7 Appendix: proofs Theorem 1 proof: Using XT X = VD2 VT and VT V = I, we have & T '−1 T X (XVi ) X X + λI ( ) D2 = V V T Vi D2 + λI D2 = Vi 2 i . Di + λ β̂ridge = (20) ! 24 Theorem 2 proof: Note that n " ! ! !Xi − αβ T Xi !2 = #n T i=1 TrXi (I i=1 #n = i=1 Tr(I − βαT )(I − αβ T )Xi − βαT )(I − αβ T )Xi XTi # Tr(I − βαT )(I − αβ T )( ni=1 Xi XTi ) = = Tr(I − βαT − αβ T + βαT αβ T )XT X = TrXT X + Trβ T XT Xβ − 2TrαT XT Xβ. (21) Since αT XT Xβ and β T XT Xβ are both scalars, we get !2 #n !! Xi − αβ T Xi ! + λ |β|2 i=1 = TrXT X − 2αT XT Xβ + β T (XT X + λ)β. (22) For a fixed α, the above quantity is minimized at & '−1 T β = XT X + λ X Xα. (23) Substituting (23) into (22) gives !2 #n !! Xi − αβ T Xi ! + λ |β|2 (24) α̂ = arg max αT XT X(XT X + λ)−1 XT Xα (25) i=1 = TrXT X − 2αT XT X(XT X + λ)−1 XT Xα. Therefore α 25 subject to αT α = 1. & '−1 T And β̂ = XT X + λ X Xα̂. By X = UDVT , we have XT X(XT X + λ)−1 XT X = V D4 VT . D2 + λ (26) D2 1 Hence α̂ = sV1 with s=1 or -1. Then β̂ = s D2 +λ V1 . 1 ! Theorem 3 proof: By the same steps in the proof of theorem 2 we derive (22) as long as α T α = Ik . Hence we have !2 # #n !! Xi − αβ T Xi ! + λ k j=1 |βj | i=1 = 2 TrXT X − 2TrαT XT Xβ + Trβ T (XT X + λ)β = TrXT X + #k j=1 $ (27) βjT (XT X + λ)βj − 2αjT XT Xβj Thus given a fixed α, the above quantity is minimized at βj = j = 1, 2, . . . , k; or equivalently & & '−1 T β = XT X + λ X Xα. % XT X + λ (28) '−1 XT Xαj for (29) Therefore α̂ = arg max TrαT XT X(XT X + λ)−1 XT Xα α (30) subject to αT α = Ik . This is an eigen-analysis problem whose solution is α̂j = sj Vj with sj =1 or -1 for j = 1, 2, . . . , k, 26 D2 j Vj for because the eigenvectors of XT X(XT X + λ)−1 XT X are V. Hence (29) gives β̂j = sj D2 +λ j j = 1, 2, . . . , k. ! Theorem 4 proof: By assumption β = U DV T with U T U = Ik and V V T = V T V = Ik . The constraint αT α = Ik is equivalent to k(k+1) 2 constraints αiT αi = 1, αiT αj = 0, i = 1, 2 . . . , k (31) j > i. (32) Using Lagrangian multipliers method, we define L=− Setting ∂L ∂αi k " i=1 βiT αi + k " 1 i=1 2 λi,i (αiT αi − 1) + k " λi,j (αiT αj ). (33) j>i = 0 gives βi = λi,i α̂i + λi,j α̂j ; or in a matrix form β = α̂Λ, where Λi,j = λj,i . Both β and α are full rank, so Λ is invertible and α = βΛ−1 . We have & ' Trα̂T β = TrΛ−1 β T β = Tr Λ−1,T V D2 V T , Ik = α̂T α̂ = Λ−1,T β T βΛ−1 = Λ−1,T V D2 V T Λ−1 . (34) (35) Let A = V T Λ−1 V , observe k " & ' & ' 2 Tr Λ−1 V D2 V T = Tr V T Λ−1 V D2 = TrAT D2 = Ajj Djj , (36) j=1 AT D2 A = Ik . 27 (37) 2 ≤ 1, Since A2jj Djj k " j=1 2 Ajj Djj ≤ k " (38) Djj . j=1 −1 The “=” is taken if only if A is diagonal and Ajj = Djj . Therefore Λ−1 = V AV T = V D−1 V T , and α̂ = βΛ = U DV T V D−1 V T = U V T . ! Theorem 5 proof: Let β̂ ∗ = (1 + λ)β̂, then we observe V̂i (λ) = means β̂i∗ . |β̂i∗ | On the other hand, β̂ = !2 ! ! ! k n ! k ! " " ! βj ! ! ! ! βj !2 " βT ! ! ! ! ! ! (α̂, β̂ ) = arg min λ1,j ! !X i − α 1 + λ X i ! + λ !1 + λ! + α,β 1 + λ! β̂ ∗ 1+λ ∗ i=1 j=1 j=1 1 (39) subject to αT α = Ik . Then by (12), we have = !2 ! ! #n !! #k !! βj !!2 #k ! ! βj ! βT X − α X + λ + λ ! ! ! ! ! i i=1 j=1 1+λ j=1 1,j 1+λ ! 1+λ i TrXT X + = TrXT X + 1 1+λ 1 1+λ $# k j=1 $# k j=1 $ T 1 X+λ βjT X 1+λ βj − 2αjT XT Xβj + λ1,j |βj |1 %% % % $ T X+λ βjT X 1+λ βj + λ1,j |βj |1 − 2TrαT XT Xβ . (α̂, β̂ ∗ ) = arg min −2TrαT XT Xβ + α,β subject to α α = Ik . T k " j=1 (40) k βjT " XT X + λ βj + λ1,j |βj |1 1+λ (41) j=1 As λ → ∞, (41) approaches (18). Thus the conclusion follows. ! 28 References Alter, O., Brown, P. & Botstein, D. (2000), ‘Singular value decomposition for genome-wide expression data processing and modeling’, Proceedings of the National Academy of Sciences 97, 10101–10106. Cadima, J. & Jolliffe, I. (1995), ‘Loadings and correlations in the interpretation of principal components’, Journal of Applied Statistics 22, 203–214. Efron, B., Hastie, T., Johnstone, I. & Tibshirani, R. (2004), ‘Least angle regression’, The Annals of Statistics 32, In press. Hancock, P., Burton, A. & Bruce, V. (1996), ‘Face processing: human perception and principal components analysis’, Memory and Cognition 24, 26–40. Hastie, T., Tibshirani, R., Eisen, M., Brown, P., Ross, D., Scherf, U., Weinstein, J., Alizadeh, A., Staudt, L. & Botstein, D. (2000), ‘’gene shaving’ as a method for identifying distinct sets of genes with similar expression patterns’, Genome Biology 1, 1–21. Hastie, T., Tibshirani, R. & Friedman, J. (2001), The Elements of Statistical Learning; Data mining, Inference and Prediction, Springer Verlag, New York. Jeffers, J. (1967), ‘Two case studies in the application of principal component’, Applied Statistics 16, 225–236. Jolliffe, I. (1986), Principal component analysis, Springer Verlag, New York. Jolliffe, I. (1995), ‘Rotation of principal components: choice of normalization constraints’, Journal of Applied Statistics 22, 29–35. 29 Jolliffe, I. T. & Uddin, M. (2003), ‘A modified principal component technique based on the lasso’, Journal of Computational and Graphical Statistics 12, 531–547. McCabe, G. (1984), ‘Principal variables’, Technometrics 26, 137–144. Misra, J., Schmitt, W., Hwang, D., Hsiao, L., Gullans, S., Stephanopoulos, G. & Stephanopoulos, G. (2002), ‘Interactive exploration of microarray gene expression patterns in a reduced dimensional space’, Genome Research 12, 1112–1120. Ramaswamy, S., Tamayo, P., Rifkin, R., Mukheriee, S., Yeang, C., Angelo, M., Ladd, C., Reich, M., Latulippe, E., Mesirov, J., Poggio, T., Gerald, W., Loda, M., Lander, E. & Golub, T. (2001), ‘Multiclass cancer diagnosis using tumor gene expression signature’, Proceedings of the National Academy of Sciences 98, 15149–15154. Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society, Series B 58, 267–288. Tibshirani, R., Hastie, T., Narasimhan, B. & Chu, G. (2002), ‘Diagnosis of multiple cancer types by shrunken centroids of gene’, Proceedings of the National Academy of Sciences 99, 6567–6572. Vines, S. (2000), ‘Simple principal components’, Applied Statistics 49, 441–451. Zou, H. & Hastie, T. (2003), Regression shrinkage and selection via the elastic net, with applications to microarrays, Technical report, Department of Statistics, Stanford University. Available at http://www-stat.stanford.edu/˜hastie/pub.htm. 30