Adachi, Kohei

Trendafilov, Nickolay

English

Open Research Online 

Open Research OnlineThe Open University’s repository of research publicationsand other research outputsPenalty-free sparse PCAConference or Workshop ItemHow to cite:Trendafilov, Nickolay and Adachi, Kohei (2014). Penalty-free sparse PCA. In: Proceedings of COMPSTAT2014: 21st International Conference on Computational Statistics, pp. 197–203.For guidance on citations see FAQs.c© [not recorded]Version: Accepted ManuscriptCopyright and Moral Rights for the articles on this site are retained by the individual authors and/or other copyrightowners. For more information on Open Research Online’s data policy on reuse of materials please consult the policiespage.oro.open.ac.ukPenalty-free sparse PCAKohei Adachi, Osaka University, adachi@hus.osaka-u.ac.jpNickolay Trendalov, Open University, Nickolay.Trendafilov@open.ac.ukAbstract. A drawback of the sparse principal component analysis (PCA) procedures usingpenalty functions is that the number of zeros in the matrix of component loadings as a wholecannot be specied in advance. We thus propose a new sparse PCA procedure in which theleast squares PCA loss function is minimized subject to a pre-specied number of zeros in theloading matrix. The procedure is called unpenalized sparse matrix PCA (USMPCA), as it doesnot use a penalty function and obtains component loadings matrix-wise, i.e., simultaneouslyrather than sequentially. The key point of USMPCA is to use the fact that the PCA lossfunction can be decomposed into sum of two terms, one of them irrelevant to loadings, andanother one being a function easily minimized under the considered cardinality constraint. Thisdecomposition makes it possible to construct an ecient alternate least squares algorithm forUSMPCA. Another useful feature is that the PC score matrix is column-orthonormal, whichhelps to dene naturally the percentage of explained variance by the sparse PCs. USMPCA isillustrated with real data examples.Keywords. Sparse component loadings, loss function decomposition, constrained matrix com-plexity.1 IntroductionFor an n observations  p variables column-centered data matrix X, principal componentanalysis (PCA) can be formulated as minimizingf(F;A) = kX   FA>k2 (1)over an nm PC score matrix F and a pm component loading matrix A, with k k2 indicatingthe squared Frobenius norm and the number of components m  min(n; p). The resultingsolution is interpreted by noting the loadings in A which quantify the relationships betweenthe p variables and m components. It is desired for A to be sparse, i.e., to have a number ofzero elements, since a sparse matrix is easily interpreted by focusing only on the variables andcomponents linked with nonzero elements. However, such sparse A cannot be obtained by thestandard PCA. For this reason, a number of modied PCA procedures have been proposed inthe last decade, which produce sparse solutions [8]. Such procedures are called sparse PCA.2Almost all existing sparse PCA procedures are using penalized approaches: they are formu-lated by combining a PCA objective function with penalty functions that penalize A to havenonzero elements. Such examples are SCoTLASS [3], SPCA [10], and sPCA-rSVD [6], wherethe relative importance of penalty functions is controlled by tuning parameters. That is, theycontrol the number of nonzero elements, which is called cardinality. Though a number of otherpenalized procedures have been developed for improving the preceding ones [4, 7, 8], they areformulated by the same format.A common drawback of the penalized sparse PCA is that the appropriate value of the tun-ing parameter which corresponds to the desired cardinality is not obvious. Thus, the penalizedsparse PCA is not convenient for users who wish to have a loading matrix with a speciednumber of zero elements. The procedures studied in [1] and [5] avoid such a diculty. Theirauthors presented ecient heuristic algorithms called "greedy" search to nd component load-ings sequentially with direct cardinality constraint. In this paper, we also propose a directlyconstrained cardinality procedure without using a penalty function. However, our proposedprocedure diers from the "greedy" search approaches in that all component are extracted si-multaneously (not sequentially), i.e., F and A are obtained matrix-wise (not column-wise). We,thus, refer to our proposed procedure as unpenalized sparse matrix PCA (USMPCA). More-over, the resulting PC scores are uncorrelated, which helps to dene naturally the percentage ofexplained variance as described in Section 4.2 Unpenalized Sparse Matrix PCAIn USMPCA, the PCA loss function (1) is minimized subject to the column-orthnormalitycondition for n 1=2F and the constraint on card(A) which denotes the cardinality of A. Thatis, USMPCA is formulated asminF;Af(F;A) = kX   FA>k2 ; subject to 1nF>F = Im and card(A) = c (2)with Im denoting the mm identity matrix and c being a specied integer.The key point of USMPCA is to use the fact that the orthonormality 1nF>F = Im allowsthe loss function (1) to be decomposed askX   FA>k2 = kX   FB> + FB>   FA>k2 = kX   FB>k2 + nkB  Ak2 ; (3)with B being the cross-product matrix of p variables m components:B =1nX>F : (4)The decomposition (3), which is derived from (X   FB>)(FB>   FA>) being the zeromatrix, shows that a simple function kB   Ak2 is only relevant to A, which allows us to easilyattain the cardinality constrained minimization of (1) as found in the next section.3 AlgorithmThe USMPCA problem (4) can be solved by alternately performing the two steps:COMPSTAT 2014 ProceedingsKohei Adachi and Nickolay T. Trendalov 3A-step minimizing (1) over A subject to card(A) = c with F being kept xed;F-step minimizing (1) over F subject to 1nF>F = Im with A kept xed.First, let us consider the A-step, which is equivalent to minimizing g(A) = kB  Ak2 undercard(A) = c, since of (3). Using A = (aij) and B = (bij), we can rewrite g(A) asg(A) = kB  Ak2 =X(i;j)2Ob2ij +X(i;j)2O?(aij   bij)2 X(i;j)2Ob2ij : (5)Here, O denotes the set of the q = pm  c indexes (i; j)'s indicating the locations of the loadingsaij to be zero, while the complement set O? contains the c (i; j)'s of nonzero aij . The inequalityin (5) shows that g(A) attains its lower limitP(i;j)2O b2ij when the non-zero loadings aij with(i; j) 2 O? are taken equal to the corresponding bij . Moreover, the limitP(i;j)2O b2ij is minimal,when O contains the indexes for the q smallest b2ij among all squared elements of B. Thus, g(A)is minimized for A = (aij) beingaij =0 if b2ij  b2[q]bij otherwise; (6)with b2[q] the qth smallest value among all b2ij .Next, let us consider the minimization in F-step. It is attained forF =pnKL> = XAL 1L>; (7)where K and L are given by the singular value decomposition (SVD) of XA dened as1pnXA = KL> (8)with K>K = L>L = Ip and  a diagonal matrix. However, it is shown in the next paragraphthat the update of F by (7) can be skipped.Using 1nF>F = Im and (4), the loss function (1) can be expanded asf(F;A) = trX>X + trAF>FA>   2trX>FA = ntrS + ntrA>A  2ntrB>A ; (9)with S = 1nX>X. Noting that (9) is a function of B and the use of (7) in (4) leads toB =1nX>XAL 1L> = SAL 1L> ; (10)we can nd that (1) or (9) is minimized for B given by (10) and this B is also used for (6): Fmay not be obtained in F-step. Moreover, the original data matrix X may not be available andonly the sample covariance matrix S suces for minimizing (1), since L 1L> in (10) can beobtained through the eigenvalue decomposition (EVD)A>SA = L2L> ; (11)following from (8): X is found to vanish in (9), (10), and (11).@ COMPSTAT 20144It should be noted that the A resulting in (6) satises trA>A = trB>A. We can use it in(9) to nd that the value of loss function (1) after the update (6) is expressed asf(A) = ntrS   ntrA>A = ntrS  fN (A) : (12)Here, fN (A) = 1  trA>A=trS is normalized so as to take a value within [0; 1], thus convenientfor checking convergence. Thus, the USMPCA algorithm can be formed as follows:1. Initialize A.2. Perform EVD (11) to obtain B with (10).3. Obtain A with (6).4. Finish if fN (A)  " ; otherwise go back to 2.Here, fN (A) denotes the change in fN (A) from the previous round. In this paper, " = 0:17 andthe algorithm is repeated fty times with random initialization. Among the resulting solutions,we select the one with the lowest fN (A) value as the optimal solution, in order to avoid localminimizers. After those procedures, F can be obtained using (7).4 Percentages of Explained VariancesThe loss function value (12) allows us to dene the goodness of the resulting A asPEV = 100trA>A=trS ; (13)with trA>A = 1nkFA>k2 following from 1nF>F = Im. The statistic (13) can be called totalpercentage of explained variance (PEV), since trS in (13) is the total variance of the variables,while trA>A = 1nkFA>k2 is the total variance of FA>, since (7) shows that F is column-centeredas X is so.The total PEV (13) can be decomposed as the sum ofPEV(j) = 100a>j aj=trS ; (14)over j = 1; :::;m. It serves as the PEV index for each component. On the other hand, the PEVfor each variable is derived from the fact that (12) can be rewritten as nPpi=1(sii   keaik2) =nPpi=1 sii(1   keaik2=sii)  0, with ea>i the ith row of A and sii the variance of variable i. Itgives the percentage of keaik2 = 1nkFea>i k2 to sii,PEV[i] = 100keaik2=sii : (15)In the same forms as (13), (14), and (15), PEV indices are dened for the standard PCA,which is formulated as minimizing (1) with 1nF>F = Im and A>A being a diagonal matrix. Thesame forms of denitions facilitate the comparison of solutions between USMPCA and PCAin goodness-of-t. Since PCA is the best rank m approximation of X, the value of the totalPEV (13) for USMPCS cannot exceed the one for PCA. However, if the former value is notsubstantially less than the latter, the USMPCS solution can be considered to be acceptable.It should be noted that USMPCS can be superior to PCA in (14) and (15), as illustrated inSection 6.1.COMPSTAT 2014 ProceedingsKohei Adachi and Nickolay T. Trendalov 55 Nonzero Loadings as CovariancesThe matrix B dened in (4) contains the covariances of p variables to m components, sinceX and F are column-centered. By taking this fact into account in (6), the nonzero loadingsin A are found to equal the corresponding covariances in B: nonzero aij equals the covariancebetween variable i and component j. It implies that the nonzero loadings equal the correlationcoecients of variables to components, when the columns of X have unit variances or S is acorrelation matrix, since of 1nF>F = Im.6 Two ExamplesThe rst example is the Pitprop data set [2] given as the correlation matrix obtained from a180 13 data matrix. We set m = 6 following the previous studies to perform USMPCA. Thesolution subject to card(A) = pm=2 = 39 is shown left in Table 1 with blank indicating zeroloadings. There, the total PEV 86.7 is found to be almost equivalent to the PEV 87.0 for PCA:USMPCA approximated the data as well as PCA with a half of loadings vanishing in the formersolution. We further performed USMPCA with decreasing card(A) one by one, to nd that thePEV for the solution with card(A) = 17 nearly exceeded 80, a benchmark percentage not beingvery lower than 87.0 for PCA. That solution is shown right in Table 1. Bold font is used for thePEV for variables and components which exceed the corresponding ones for PCA. One notesthat the USMPCA components with j = 4; 5; 6 explain more variance than the PCA ones.Vars USMPCA: card(A) = 39 USMPCA: card(A) = 17 PCA1 2 3 4 5 6 PEV 1 2 3 4 5 6 PEV PEVtopdiam .86 .40 90.4 .89 79.2 90.9length .90 .33 92.0 .91 82.9 92.5moist .98 -.10 97.5 .96 92.4 97.8testsg .90 -.40 97.5 .94 88.6 97.5ovensg -.17 -.93 88.7 .81 65.0 86.8ringtop .32 .19 .59 -.55 0.29 87.1 .37 .79 76.6 86.4ringbut .61 .61 -.41 -.14 93.1 .67 .62 83.4 92.7bowmax .54 .15 -0.60 67.8 .61 .51 63.4 68.4bowdist .75 .15 -0.20 62.9 .80 63.5 64.0whoris .66 .33 -0.38 -.42 86.8 .75 .44 75.1 87.2clear .15 .97 96.9 -.98 95.3 95.9knots -.11 .25 .25 0.80 77.5 -.92 85.5 80.4diaknot .15 .00 -.87 .10 0.31 88.6 -.96 91.6 90.7PEV 25.8 17.1 12.5 12.0 10.5 8.8 86.7 29.1 13.9 12.8 8.8 8.6 7.0 80.2 87.0PEVPCA 32.5 18.3 14.5 8.5 7.0 6.3 87.0 32.5 18.3 14.5 8.5 7.0 6.3 87.0Table 1. USMPCA solutions for Pitprop data with PCA's PEV in the nal row and column.The variables are well clustered with every variable loading only one or two components. Itmakes sense to compare the USMPCA solutions with the classic (subjective) interpretation ofthe Pitprop component loadings [2], which is summarized in Table 2. The adopted notationsmean that the rst component is determined by topdiam, length, ringbut, bowmax, bowdistand whoris, the second { by moist and testsg, and etc. The ringbut value for component four in[2, Table 4, p.229] seems incorrect, by inspecting the corresponding eigenvalue. The corrected"classic" interpretation is given in [8], where ringbut is dropped o the fourth component.Clearly, the USMPCA solution with card(A) = 17 suggests identical interpretation of therst three components as the one given in [2, p.230]. The fourth component is, indeed, acontrast, but between clear and whoris. The fth component is also a contrast between knots@ COMPSTAT 20146and bowmax, and the sixth component is a direct measure of diaknot (the average diameter ofthe knots in inches).Vars 1 2 3 4 5 6topdiam xlength xmoist xtestsg x xovensg x xringtop x xringbut x xbowmax xbowdist xwhoris xclear xknots xdiaknot xTable 2. Classic interpretation of the Pitprop component loadings [2, p.229-30].The second example concerns the gene expression data matrix of n = 17 time points byp = 384 genes presented by [9] and available at http://faculty.washington.edu/kayee/pca.The 384 genes are categorized into ve phases of cell cycles, with each phase containing 67,135, 75, 52, and 55 genes, respectively. It suggests m = 5, but this choice yielded one trivialcomponent in preliminary trials. We thus reduced m to 4. For card(A), we rst used theinteger nearest to the one-third of pm, then increase card(A) one by one to nd that the totalPEV of the solution with card(A) = 538 = 0:35pm nearly exceeds a benchmark 70, which is notconsiderably lower than the PEV 81.2 for PCA withm = 4. The resulting A with card(A) = 538are presented block-wise in Figure 1. There, the blocks correspond to the ve phases, with theblock for the second one divided into two, and positive/negative nonzero loadings represented aslled squares/triangles, respectively. The solution is considered to be reasonable, as each phasehas a unique feature of loadings: [a] Phases 1, 2, and 4 are characterized by positive loadingsfor Components 1, 2, and 3, respectively; [b] Phases 5 are characterized by positive loadings forComponent 4 and negative ones for 2; [c] Phases 3 consists of the genes positively loaded byComponent 2 or 3 and by both.7 Final RemarksIn this paper, we proposed the penalty-free sparse PCA procedure USMSPCA and presented itsalternate least squares algorithm. An advantage of USMPCA over the penalized sparse PCA isthat the cardinality of loadings can be set to a specied integer in advance. For that integer wecan use the one conceived easily such as a half or the one-third of the number of loadings, whichcan be exibly changed for nding a better solution, as illustrated in the examples. There, it wasalso illustrated that a solution obtained can be validated by comparing the PEV value with thecorresponding one for the standard PCA. The reasonableness of this PEV comparison followsfrom that the PEV indices for USMPCA are dened in the same manner as in PCA.AcknowledgementThis work is supported by a grant RPG-2013-211 from The Leverhulme Trust, UK.COMPSTAT 2014 ProceedingsKohei Adachi and Nickolay T. Trendalov 7Figure 1. USMPCA solution for gene expression data with blank indicating zeroBibliography[1] d'Aspremont, A., Bach, F., & Ghaoui, L. E. (2008) Optimal solutions for sparse principalcomponent analysis, Journal of Machine Learning Research, 9, 1269-1294.[2] Jeers, J. N. R. (1967). Two case studies in the application of principal component analysis.Applied Statistics, 16, 225-236.[3] Jollie, I. T., Trendalov, N. T., & Uddin, M. (2003). A modied principal componenttechnique based on the LASSO. Journal of Computational and Graphical Statistics, 12,531-547.[4] Journe M, Nesterov Y, Richtrik P, & Sepulchre R (2010) Generalized power method forsparse principal component analysis, Journal of Machine Learning Research, 11, 517-553.@ COMPSTAT 20148[5] Moghaddam, B., Weiss, Y., & Avidan, S. (2006). Spectral bounds for sparse PCA: Exactand greedy algorithms. Advances in Neural Information Processing System, 18, 915-922.[6] Shen, H. & Huang, J. Z. (2008). Sparse principal component analysis via regularized lowrank matrix approximation. Journal of Multivariate Analysis, 99, 1015 - 1034.[7] Qi, X., Luo, R., & Zhao, H. (2013). Sparse principal component analysis by choice of norm.Journal of Multivariate Analysis, 114, 127-160.[8] Trendalov, N. T. (2013). From simple structure to sparse components: a review. Compu-tational Statistics, published on line DOI:10.1007/s00180-013-0434-5.[9] Yeung, K. Y., & Ruzzo, W. L. (2001). Principal component analysis for clustering geneexpression data. Bioinformatics, 17, 763-774.[10] Zou, D. M., Hastie, T., & Tibshirani, R. (2006). Sparse principal component analysis.Journal of Computational and Graphical Statistics, 15, 265-286.COMPSTAT 2014 Proceedings

Penalty-free sparse PCA

http://oro.open.ac.uk/41228/1/FreeKANT.pdf

Penalty-free sparse PCA

Abstract

Similar works

Full text

Available Versions

Open Research Online