rs-Sparse principal component analysis: A mixed integer nonlinear programming approach with VNS

rs-Sparse principal component analysis: A mixed integer nonlinear programming approach with VNS $ Emilio Carrizosa, Vanesa Guerrero n Instituto de Matemáticas de la Universidad de Sevilla, Fac. de Matemáticas, Avda Reina Mercedes s/n, 41012 Sevilla, Spain article info Available online 3 May 2013 Keywords: Sparse principal component analysis Variable neighborhood search Nonlinear mixed integer programming abstract Principal component analysis is a popular data analysis dimensionality reduction technique, aiming to project with minimum error for a given dataset into a subspace of smaller number of dimensions. In order to improve interpretability, different variants of the method have been proposed in the literature, in which, besides error minimization, sparsity is sought. In this paper we formulate as a mixed integer nonlinear program the problem of finding a subspace with a sparse basis minimizing the sum of squares of distances between the points and their projections. Contrary to other attempts in the literature, with our model the user can fix the level of sparseness of the resulting basis vectors. Variable neighborhood search is proposed to solve the problem obtained this way. Our numerical experience on test sets shows that our procedure outperforms benchmark methods in the literature, both in terms of sparsity and errors. &2013 Elsevier Ltd. All rights reserved. 1. Introduction and literature review Principal component analysis (PCA) was first introduced by [20] as a method for projecting a set of points u 1 ;…;u p ∈R n onto a lower dimensional space in such a way that the distances between the points and their projections are minimized, see e.g., [14,9]. Given kvectors c 1 ;…;c k ∈R n , let π fc 1 ;…;c k g denote the projection onto the linear space Lðfc 1 ;…;c k gÞ spanned by the vectors c 1 ;…;c k π fc 1 ;…;c k g ðuÞ¼arg minf∥u−z∥:z∈Lðfc 1 ;…;c k gÞg: The aim of PCA is to find a set of k≤northonormal vectors c 1 ;…;c k such that the mean squared distance between the points in the dataset fu 1 ;…;u p gand their projections onto the vector space Lðfc 1 ;…;c k gÞ generated by fc 1 ;…;c k gis minimized. In other words, the following mixed integer nonlinear program (MINLP) is to be solved min c 1 ;…;c k :orthonormal 1 p∑ p h¼1 ∥u h −π fc 1 ;…;c k g ðu h Þ∥ 2 The optimal solutions, c n ¼ðc n 1 ;…;c n k Þ,arecalledprincipal components (PC). The main drawback of this dimensionality-reduction technique is interpretability: interpreting the projections is usually quite difficult due to the fact that most of the original variables are present in each c n i ,i¼1,…,k, i.e., the PC are not sparse. Interpretability is improved if some loadings (the coefficients of PCs) are zero, and this has been pursued in different papers. Afirst attempt for achieving this is rounding, by considering all loadings smaller than some threshold absolute value as zero. However, this procedure has been shown to be misleading, see [4].Varimax rotation, see [15], has also been proposed, but it hardly ever achieves the aim of an easier interpretation despite losing orthogonality of the loadings and uncorrelation of the components. Some authors relate the notion of simplicity to the fact that loadings belong to the set of integers. Such idea was developed in [26] and later in [23], who called their method simple component analysis (SCA). SCA allows the user, under his or her criterion, to modify a simple system of components in order to improve interpretability. However, SCA does not yield either orthogonal or uncorrelated components. An exploratory approach to SCA was presented in [1] for achieving orthogonality. Also in [24], SCA is modified by using genetic algorithms. Another way of obtaining sparsity is by constraining the number of non-zero loadings in each PC. In this line, Ref. [5] proposes a convex relaxation method based on semidefinite programming, which does not preserve orthogonality or uncorrelation, as while [6] a branch-and-bound approach lets the user choose between keeping orthogonality or uncorrelation. A related approach is presented in [11], in which a bound on the sum of the absolute values of the loadings is added, combining this way the lasso paradigm, [25],withPCA. In [28], PCA is formulated as a regression-type optimization problem. Sparse loadings are achieved by imposing elasticnet constraint, [27],ageneralizationoflasso, on the regression coefficients. The sparse PCs obtained from sparse principal component analysis (SPCA) are neither orthogonal nor uncorrelated. See another lasso-based model in [21]. Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/caor Computers & Operations Research 0305-0548/$ - see front matter &2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cor.2013.04.012 ☆ This research is supported by Grants MTM2012-36136 (Ministerio de Economía y Competitividad, Spain) and FQM329 (Junta de Andalucía), both co-funded by EU ERD funds. n Corresponding author. Tel.: +34 687153536. E-mail addresses: [email protected] (E. Carrizosa), [email protected] (V. Guerrero). Computers & Operations Research 52 (2014) 349–354 Some other authors address sparsity in PCA with tools of feature selection. This is the case of [17], who called principal variables those which are considered to be relevant. Rejection methods are introduced and tested in artificial and real datasets in [12,13]. Also the well-known variable neighborhood search (VNS), [18,8], is adapted in [3] for developing two heuristics for feature selection in PCA. More recently, an exact method (branch-and- bound) is presented in [19]. In this paper we present a new procedure for obtaining a set of orthonormal vectors which are as sparse as desired, and minimizing the sum of squared distances between the data points and their projections. This is achieved by heuristically solving a MINLP with VNS. Our procedure has been tested in five datasets and compared with some benchmark methods in the literature, namely, VARIMAX, SCA and SPCA. In all cases we have outperformed the results given by these other procedures, both in terms of sparsity and errors. The remainder of the paper is structured as follows. In Section 2, our problem is stated. Section 3 analyses an important particular case, namely, the case in which sparsity is obtained by forcing each original variable to appear with nonzero value in at most one component. The model is extended in Section 4 to give more freedom to the user to control sparsity. Section 5 reports our numerical experience, in which our method is compared against state-of-the-art procedures for sparse PCA, showing the advantages of our procedure. The paper ends with Section 6,inwhichsome concluding remarks and future lines of research are outlined. 2. Problem statement The special feature of our sparse PCA method against the classical one is the fact that we force principal components to have the most of their coordinates equal to zero. Our proposal, as in the classical PCA, tries to minimize the sum of the squares of the distances between the dataset and its projection on the vector space generated by the principal components. However, those k principal components will be forced to satisfy some sparsity constraints, which are basically that each variable is nonzero in at most rcomponents, and each component has at most snonzero elements. Since we have two parameters to control the model and make principal components sparse, we call our proposal rs-sparse principal component analysis,rs-SPCA. So, we formulate the sparse problem as a MINLP by defining the following set of variables. First, let c 1 ;…;c k denote the k vectors in R n which are sought as principal components. We call c il the l-th coordinate of the i-th principal component c i , where i¼1…kand l¼1…n. Binary variables z il are defined as z il ¼1ifc il ≠0 0 else i¼1;…;k;l¼1;…;n:  This assignment will let us control the sparsity of the principal components c 1 ;…;c k . With the previous notation, our aim is to solve the following MINLP: min 1 p∑ p h¼1 ∥u h −π fc 1 ;…;c k g ðu h Þ∥ 2 s:t: c ⊤ i c j ¼δ ij ∀i;j¼1;…;k jc il j≤z il ∀i¼1;…;k;l¼1;…;n ∑ k i¼1 z il ≤r∀l¼1;…;n ∑ n l¼1 z il ≤s∀i¼1;…;k z il ∈f0;1g∀i¼1;…;k;l¼1;…;n 8 > > > > > > > > > > > < > > > > > > > > > > > : where δ ij ¼1ifi¼jand 0 otherwise. Every constraint plays an important role in the achievement of the sparsity. The first one is given by the classical PCA, which forces orthonormality. The second constraint is obtained by imposing that, if z il ¼0, then c il ¼0, i¼1,…,k,l¼1,…,n. Finally, the constraints which really control the sparsity are the last two. In the third one, each variable is forced to appear in at most r components, with the fourth one every component guarantees to have no more than snonzero loadings. Denote the variance–covariance matrix (or correlation matrix) by V. The previous problem is equivalent to max ∑ k i¼1 c ⊤ i Vc i s:t: c ⊤ i c j ¼δ ij ∀i;j¼1;…:k jc il j≤z il ∀i¼1;…;k;l¼1;…;n ∑ k i¼1 z il ≤r∀l¼1;…;n ∑ n l¼1 z il ≤s∀i¼1;…;k z il ∈f0;1g∀i¼1;…;k;l¼1;…;n 8 > > > > > > > > > > > < > > > > > > > > > > > : ð1Þ We must point out that principal components given by (1) will be orthonormal but no longer uncorrelated. In order to compare the results obtained with (1) against other procedures, a goodness criterion is introduced. We suggest proceeding in the same way as in the classical problem, which means to calculate the percentage of total variance explained by, in this case, rs-sparse principal components, as follows: f¼1 trðVÞ∑ k i¼1 c ⊤ i Vc i 100;ð2Þ where trðVÞis the trace of the variance–covariance or correlation matrix V. Other related measures can be found in [7]. In the following sections two special cases of (1), namely, r¼1, and the general case, r≥1, will be analyzed. 3. The case r¼1 Imposing in (1) the condition r¼1 we are forcing that each variable appears just once. This induces a partition in the set of variables, i.e., those assigned to each component representing one cluster. Hence, interpretability will be high, though perhaps at expense of a decrease in variance explained. In this case we have directly the orthonogonality of the components, so the constraint c ⊤ i c j ¼δ ij becomes in c ⊤ i c i ¼1, ∀i¼1…k and (1) is reduced to max ∑ k i¼1 c ⊤ i Vc i s:t: c ⊤ i c i ¼1∀i¼1;…:k jc il j≤z il ∀i¼1;…;k;l¼1;…;n ∑ k i¼1 z il ≤1∀l¼1;…;n ∑ n l¼1 z il ≤s∀i¼1;…;k z il ∈f0;1g∀i¼1;…;k;l¼1;…;n 8 > > > > > > > > > > > < > > > > > > > > > > > : ð3Þ However, if z il could be fixed by any procedure, the resulting problem, P(z), turns into a very familiar one, namely max ∑ k i¼1 c ⊤ i Vc i s:t: c ⊤ i c i ¼1∀i¼1;…;k c il ¼0∀i¼1;…;k;l¼1;…;n:z il ¼0 (PðzÞ E. Carrizosa, V. Guerrero / Computers & Operations Research 52 (2014) 349–354350 P(z) can be expressed in a much more manageable form. Indeed, the problem is separable and can thus be split into kindependent subproblems. In other words, P(z) can be written as ∑ k i¼1 fmax c ⊤ i Vc i g s:t: c ⊤ i c i ¼1∀i¼1;…;k c il ¼0∀i;l:z il ¼0 (ð4Þ For z¼ðz 1 ;…;z k Þgiven, denote by V z i i ,i¼1,…,k, the symmetric matrix obtained from Vby deleting all rows and columns lsuch that z il ¼0, l¼1,…,n. Let d z i i denote an optimal solution max d ⊤ i V z i i d i s:t:fd ⊤ i d i ¼1ð5Þ Moreover, define the vector c z i i as c z i il ¼d z i il if z il ≠0 0 else ( It immediately follows that the so-constructed c z i i ,i¼1,…,k solve (4), and thus P(z). In other words, P(z) is solved by solving kproblems of type (5), which, since V z i i is symmetric, its optimal solution corresponds to any unit eigenvector d z i i associated with the highest eigenvalue λ z i i of V z i i . Nevertheless, we must find a way for fixing those binary variables which make our initial problem turn into P(z). 3.1. Problem resolution: fixing an initial feasible solution Given V, the correlation or variance–covariance matrix, classical principal components c n ¼ðc n 1 ;…;c n k Þcan be easily computed. We will find z il ,i¼1,…,k,l¼1,…,n, such that the sum of the absolute values of loadings of principal components are maximized by imposing constraints which are required in (3). A constraint forcing at least one variable appears in each principal component will be also included in this formulation. So, we want z il to satisfy max ∑ k i¼1 ∑ n l¼1 jc n il jz il s:t: ∑ k i¼1 z il ≤1∀l¼1;…;n ∑ n l¼1 z il ≤s∀i¼1;…;k ∑ n l¼1 z il ≥1∀i¼1;…;k z il ≥0∀i¼1;…;k;l¼1;…;n 8 > > > > > > > > > > > < > > > > > > > > > > > : ð6Þ Problem (6) has a flow problem constraint structure, and it attains its optimal value at some point z n with all coordinates z n il ∈f0;1g. Such z n will be used as values for zin P(z). Combining problems (6) and P(z), a solution (c,z), feasible in terms of sparsity and orthonormality, is achieved. However, it may not be optimal for (3), but it could be a good starting point for a search procedure such as VNS. 3.2. Problem resolution: improving the solution via a VNS algorithm VNS is a metaheuristic for avoiding been trapped in a local optimum. It does not assure reaching the global one, but it sometimes achieves a good improvement on the solution initially given. Now we describe how VNS is customized for our problem (3). As we have seen above, problem (3) is a MINLP in the variables (c,z). For z¼ðz 1 ;…;z k Þgiven, the optimal solution to P(z), cðzÞ¼ðc 1 ðz 1 Þ;…;c k ðz k ÞÞ, is obtained by calculating eigenvectors of the matrices V z i i defined above. Our algorithm will perform a neighborhood search in the z-space, yielding ðcðzÞ;zÞas solution to problem (3). Since zconsists of the binary variables z il , a natural neighborhood structure in the z-space is defined as follows: given z, satisfying the constraints in (6), for each radius ρ¼1;…;n, the neighborhood N ρ ðzÞof zis defined as the set of all z′fulfilling the constraints in (6) and obtained form zby exchanging ρrows. For instance, consider the case n¼4, k¼2, r¼1, s¼2, and z given by. z¼ 10 00 01 01 0 B B B @ 1 C C C A ð7Þ The elements in the neighborhood N 1 ðzÞof radius 1 of zare the five matrices z′given in Table 1. Observe that each element zinduces a partition of the features space 1;2;…;ninto kþ1 clusters: For l¼1,…,n, variable lis associated with cluster i(i¼1,…,k)ifz il ¼1; if z il ¼0 for all i¼1, …,k, then variable lis associated with the extra cluster k+1. Hence, the elements z′in a neighborhood of radius ρof z correspond to those partitions of the features obtained by exchanging ρfeatures form the partition in z. For instance, zgiven in (7) induces the clustering {1}, {3,4}, {2} for the features {1,2,3,4}; the partitions in the neighborhood of radius 1 are shown in Table 2. As a stopping criterion, let the algorithm looks for no more than Q max solutions in each neighborhood and impose that the total running time does not exceed a given T max . Assuming these conditions, the implementation of the algorithm is done as is shown in Fig. 1. 4. General r Let us now consider problem (1). We propose to proceed as we did in the case r¼1. Firstly, we obtain the principal components c n ¼ðc n 1 ;…;c n k Þand then we obtain an approximation to c n which is feasible (in terms of the number of nonzero elements) and as close as possible to c n . To do that, we find an initial zby solving the following flow problem (which is identical to problem (6), except that ris an arbitrary integer number, r≤kÞ max ∑ k i¼1 ∑ n l¼1 jc n il jz il Table 1 Example N 1 ðzÞ. z′Perturbed row –1 z′ 1 ¼ 10 10 01 01 0 B B B @ 1 C C C A 2 z′ 2 ¼ 10 00 00 01 0 B B B @ 1 C C C A ,z′ 3 ¼ 10 00 10 01 0 B B B @ 1 C C C A 3 z′ 4 ¼ 10 00 01 00 0 B B B @ 1 C C C A ,z′ 5 ¼ 10 00 01 10 0 B B B @ 1 C C C A 4 E. Carrizosa, V. Guerrero / Computers & Operations Research 52 (2014) 349–354 351 s:t: ∑ k i¼1 z il ≤r∀l¼1…n ∑ n l¼1 z il ≤s∀i¼1…k ∑ n l¼1 z il ≥1∀i¼1…k z il ≥0∀i;l 8 > > > > > > > > > > > < > > > > > > > > > > > : ð8Þ Let z n be the optimal solution of the flow linear program (8), and then we find cðz n Þby solving the nonlinear continuous program max ∑ k i¼1 c ⊤ i Vc i s:t: c ⊤ i c j ¼δ ij ∀i;j c il ¼0∀i;l:z n il ¼0 (ð9Þ Orthogonality constraints cannot be simplified anymore, contrary to what we did in the previous section. Hence, a generalpurpose optimization routine is to be used to solve this problem. 5. Computational experiments In this section we analyze the empirical behavior of our rs- SPCA, which is compared against benchmark procedures in the literature, namely rotation methods (varimax criterion) [15], SCA, [23] and SPCA, [28]. Those different procedures for achieving simplicity will be carried out on datasets proposed in [22]. The numerical experience will show that rs-SPCA provides principal components which outperform, in terms of sparsity and error minimization, the ones given by the procedures cited above. 5.1. Data description Hearing loss data, come from [9], who considered the first four principal components in their study. It consists of eight measurements of hearing loss taken on 100 males who had no indication of hearing difficulties. The correlation matrix is available in Rpackage sca. Reflexes data were described in [14], gives measurements of strength of reflexes at 10 sites of the body, taken form 143 individuals. As the previous data, reflexes correlation matrix can also be found in Rpackage sca. Pitprop, data [10], is an example which has been proven the difficult achievement of simplicity in the past. About 13 variables arised from a study on the strength of pitprops cut from home-grown timber, measured in 180 observations. Six principal components were calculated. Correlation matrix can be found in Rpackage elasticnet. Associated Movements data were introduced in [16], where 10 measures of the amount of associated movements observed during some finemotor and grossmotor tasks in 484 children and adolescents were taken. Five principal components are proposed to be calculated. Muscle Strength data consist of 51 maximal isometric muscle strength measurements made on all areas of the body from 569 healthy subjects. The first six principal components will be calculated. See [22]. 5.2. Experiments description and results In order to show the power of rs-SPCA, several computational experiments have been carried out. The measure of goodness used is the percentage of the variance explained by the basis selected, as given by (2). We compare our rs-SPCA (r¼1) against benchmark methods. We have implemented our procedure in Matlab, using the routines linprog to numerically solve the flow problems (6) and (8), while the routine fmincon is used to solve (9). VNS was implemented as described in Section 3.2,with Q max ¼5 trials within each neighborhood, and a time limit Table 2 Partition. z′Partition z′ 1 {1,2}, {3,4}, f∅g z′ 2 {1}, {4}, {2,3} z′ 3 ({1,3}, {4}, {2} z′ 4 {1}, {3}, {2,4} z′ 5 {1,4}, {3}, {2} -Initialization: - find c∗, principal components. - find z0, optimal solution to (6). - find c0, optimal solution to P(z0). -Repeat the following steps until a maximum running time Tmax is achieved. (a) q←1, ρ←1 (b) Repeat the following steps until ρ>n. (i) Shaking: generate a ˜z= (˜z1,..., ˜zk)inNρ(z0) and find c(˜z) = (c1(˜z1),..., ck(˜zk)), optimal solution to P( ˜z). (ii) Neighborhood change: if k i=1 ci(˜zi)Vc i(˜zi)> k i=1 c0 iVc 0 i, then move ( z0←˜z), continue the search in the neighborhood N1(˜z) and set ( q←1); otherwise, set q←q+1 . If q+1 >Qmaxthen set q←1 and ρ←ρ+1. Go to (b). Fig. 1. VNS algorithm. E. Carrizosa, V. Guerrero / Computers & Operations Research 52 (2014) 349–354352 T max of 15, 30, 45, 60 and 100 s respectively for the different datasets. Varimax, SCA and SPCA algorithms are included in stats, sca and elasticnet packages of Rsoftware, respectively. For SCA, two parameters, namely, the number of block and difference components, must be chosen a priori. We have followed [22] to chose such parameters. SPCA also depends on two parameters. For the quadratic penalty, we follow [28], and fix it to 0, since the number of variables is lower than the number of observations in all our datasets. The second parameter has been chosen so that the soobtained principal components have the same non-zero coordinates than the ones given by our method. In other words, we are giving valuable information to the user, since we are already giving a clue on the distribution of non-zeros. The results are given in Table 3. The first column gives the name of the dataset. Columns 2 and 3 give respectively n(the number of variables) and k(the number of components to be obtained). Then, for the different methods tested (PCA, VARIMAX, SCA, SPCA and our rs-SPCA), we represent the sparsity (as the percentage of zeros in the components) and the percentage of variance explained by the components, f, as given by (2). We conclude from Table 3 that, in terms of sparsity, our procedure clearly outperforms PCA as well as the methods seeking sparseness. This happens for all methods except SPCA, for which, thanks to the extra information we provided, gives identical sparsity results to ours, but at the expense of a lower variance explained. Unlike in our proposal, orthogonality is not ensured in SPCA or the remaining sparse procedures. In other words we have proposed a method with a high sparsity and a low loss in terms of the variance explained, as compared with the classical PCA. A second analysis has been done to study how sensitive our procedure is with respect to the choice of kand s. The performance of our method on the different test sets is depicted in Figs. 2–6, which show ffor r¼1andvaryingsand k.Thesefigures correspond to the case r¼1, as analyzed in Section 3. For the case of an arbitrary r, as discussed in Section 4, the local search routine may fail to find a feasible solution. Table 4 shows the percentages of feasible solutions obtained for the different datasets, when ris varying between 1 and k, and svaries between 1 and n. We conclude that failing to obtain feasibility is rare in small-dimensional datasets, while for the larger set (n¼51), a feasible solution was obtained in less than 40% of the instances. 6. Conclusions In this paper we have introduced a new dimensionality reduction method which ensures sparsity of the procedure. We have modeled the problem as a MINLP, heuristically solved by using VNS. The numerical experience reported shows that our procedure outperforms, in terms of error minimization and sparsity, benchmark methods in the literature. Several research lines remain open. Although for high-dimensional problems heuristics such as VNS seem to be the only feasible strategy, it would be interesting to be able to solve the MINLPs obtained with exact methods. Plugging directly the problem into state-of-the-art MINLP branch- Table 3 Percentages of zeros of the components and percentages of explained variance. Datasets nkSparsity PCA VARIMAX SCA SPCA rs- SPCA Hearing Loss 8 4 %zeros 0 0 12.50 75.00 75.00 f87.37 68.40 85.40 84.08 85.37 Reflexes 10 5 %zeros 0 16.00 52.00 80.00 80.00 f97.05 72.20 91.50 96.17 96.70 Pitprop 13 6 %zeros 1.28 7.69 80.77 83.33 83.33 f87.00 78.90 74.80 71.99 76.85 Movements 22 5 %zeros 0 2.27 20.45 75.00 75.00 f55.00 43.20 53.80 49.84 53.60 Musclestrength 51 6 %zeros 1.63 8.50 34.64 80.07 80.07 f70.40 70.39 68.10 60.00 61.39 0 1 2 3 4 5 6 7 8 9 10 20 30 40 50 60 70 80 90 Hearing loss s f k = 1 k = 2 k = 3 k = 4 Fig. 2. fvarying sand kfor Hearingloss dataset. 0 1 2 3 4 5 6 7 8 9 10 11 10 20 30 40 50 60 70 80 90 100 Reflexes s f k = 1 k = 2 k = 3 k = 4 k = 5 Fig. 3. fvarying sand kfor Reflexes dataset. 0 2 4 6 8 10 12 14 10 20 30 40 50 60 70 80 Pitprop s f k = 1 k = 2 k = 3 k = 4 k = 5 k = 6 Fig. 4. fvarying sand kfor Pitprop dataset. E. Carrizosa, V. Guerrero / Computers & Operations Research 52 (2014) 349–354 353 and-bound methods such as Couenne, [2], was unsuccessful, since no solution was provided in reasonable times. Our method, as well as its competitors, seeks orthogonal vectors, but the output is not uncorrelated. Seeking a basis (not necessarily orthogonal) yielding uncorrelated projections should lead to more challenging optimization problems, which are now under study. Since both uncorrelation and orthogonality are of interest, addressing a biobjective problem trading-off components correlation and orthogonality, but maintaining a high variance explained, is also a challenging problem which deserves a deeper analysis. Acknowledgment We are grateful to Professor Valentin Rousson for providing us the muscle strength correlation matrix. References [1] Anaya-Izquierdo K, Critchley F, Vines K. Orthogonal simple component analysis: a new, exploratory approach. The Annals of Applied Statistics 2011; 5(1):486–522. [2] Belotti P. Couenne: a user's manual. Technical Report, Lehigh University; 2009. [3] Brusco M, Singh R, Steinley D. Variable neighborhood search heuristics for selecting a subset of variables in principal component analysis. Psycometrika 2009;74(4):705–26. [4] Cadima J, Jolliffe T. Loading and correlation in the interpretation of principal components. Journal of Applied Statistics 1995;22:203–14. [5] d'Aspremont A, Ghaoui L, Jordan M, G. L. Direct formulation of sparse PCA using semidefinite programming. SIAM Review 2007;49(3):434–48. [6] Farcomeni A. An exact approach to sparse principal component analysis. Computational Statistics 2009;24(583–604). [7] Gervini D, Rousson V. Criteria for evaluating dimension-reducing components for multivariate data. The American Statitician 2004;58:72–6. [8] Hansen P, MladenovićN. Variable neighbourhood search: methods and applications. Annals of Operations Research 2010;175:35–45. [9] Jackson JE. A user's guide to principal components. New York: Wiley; 1991. [10] Jeffers J. Two case studies in the application of principal component analysis. Applied Statistics 1967;16(3):225–36. [11] Jollife IT, Trendafilov NT, Uddin M. A modified principal component technique based on the Lasso. Journal of Computational and Graphical Statistics 2003;12: 531–47. [12] Jolliffe IT. Discarding variables in a principal component analysis I: artificial data. Applied Statistics 1972;21(2):160–73. [13] Jolliffe IT. Discarding variables in a principal component analysis II: artificial data. Applied Statistics 1973;22(1):21–31. [14] Jolliffe IT. Principal component analysis. 2nd ed. New York: Springer; 2002. [15] Kaiser HF. The varimax criterion for analytic rotation in factor analysis. Psychometrica 1958;23(3):187–200. [16] Largo R, Caflisch J, Hug F, Muggli K, Molnar A, Molinari L. Neuromotor development from 5 to 18 years. Part 2: associated movements. Developmental Medicine on Child Neurology 2001;43:444–53. [17] McCabe GP. Principal variables. Technometrics 1984;26(2):137–44. [18] MladenovićN, Hansen P. Variable neighborhood search. Computers and Operations Research 1997;24:1097–100. [19] Pacheco J, Casado S, Porras S. Exact methods for variable selection in principal component analysis: guide functions and pre-selection. Computational Statistics & Data Analysis 2013;57:95–111. [20] Pearson K. On lines and planes of closest fit to systems of points in space. Philosophical Magazine 1901;2:559–72. [21] Qi X, Luo R, Zhao H. Sparse principal component analysis by choice of norm. Journal of Multivariate Analysis 2013;114:127–60. [22] Rousson V, Gasser T. Some case studies of simple component analysis, 〈http:// www.biostat.uzh.ch/research/manuscripts/scacases.pdf〉; 2003. [23] Rousson V, Gasser T. Simple component analysis. Applied Statistics 2004; 53(4):539–55. [24] Sabatier R, Reynès C. Extensions of simple component analysis and simple linear discriminant analysis using genetic algorithms. Computational Statistics & Data Analysis 2008;52:4779–89. [25] Tibshirani R. Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society, Series B 1996;58:267–88. [26] Vines S. Simple principal components. Applied Statistics 2000;49(4):441–51. [27] Zou H, Hastie T. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society Series B 2005;(67):301–20. [28] Zou H, Hastie T, Tibshirani R. Sparse principal component analysis. Journal of Computational and Graphical Statistics 2006;15(2):265–86. 0 5 10 15 20 0 10 20 30 40 50 60 Movements s f k = 1 k = 2 k = 3 k = 4 Fig. 5. fvarying sand kfor Movements dataset. 0 5 10 15 20 25 30 35 40 45 50 0 10 20 30 40 50 60 Muscle strength s f k = 1 k =2 k = 3 k = 4 k = 5 k = 6 Fig. 6. fvarying sand kfor Musclestrength dataset. Table 4 Percentage of feasibility in rs-SPCA general problem. Datasets Feasibility Hearing Loss 100% Reflexes 100% Pitprop 94.87% Movements 92.04% Musclestregth 38.89% E. Carrizosa, V. Guerrero / Computers & Operations Research 52 (2014) 349–354354