Sparse PCA

Statistical analysis technique

Sparse principal component analysis (SPCA or sparse PCA) is a technique used in statistical analysis and, in particular, in the analysis of multivariate data sets. It extends the classic method of principal component analysis (PCA) for the reduction of dimensionality of data by introducing sparsity structures to the input variables.

A particular disadvantage of ordinary PCA is that the principal components are usually linear combinations of all input variables. SPCA overcomes this disadvantage by finding components that are linear combinations of just a few input variables (SPCs). This means that some of the coefficients of the linear combinations defining the SPCs, called loadings,[note 1] are equal to zero. The number of nonzero loadings is called the cardinality of the SPC.

Mathematical formulation

Consider a data matrix, X {\displaystyle X} , where each of the p {\displaystyle p} columns represent an input variable, and each of the n {\displaystyle n} rows represents an independent sample from data population. One assumes each column of X {\displaystyle X} has mean zero, otherwise one can subtract column-wise mean from each element of X {\displaystyle X} . Let Σ = 1 n 1 X X {\displaystyle \Sigma ={\frac {1}{n-1}}X^{\top }X} be the empirical covariance matrix of X {\displaystyle X} , which has dimension p × p {\displaystyle p\times p} .

Given an integer k {\displaystyle k} with 1 k p {\displaystyle 1\leq k\leq p} , the sparse PCA problem can be formulated as maximizing the variance along a direction represented by vector v R p {\displaystyle v\in \mathbb {R} ^{p}} while constraining its cardinality:

max v T Σ v subject to v 2 = 1 v 0 k . {\displaystyle {\begin{aligned}\max \quad &v^{T}\Sigma v\\{\text{subject to}}\quad &\left\Vert v\right\Vert _{2}=1\\&\left\Vert v\right\Vert _{0}\leq k.\end{aligned}}} Eq. 1

The first constraint specifies that v is a unit vector. In the second constraint, v 0 {\displaystyle \left\Vert v\right\Vert _{0}} represents the 0 {\displaystyle \ell _{0}} pseudo-norm of v, which is defined as the number of its non-zero components. So the second constraint specifies that the number of non-zero components in v is less than or equal to k, which is typically an integer that is much smaller than dimension p. The optimal value of Eq. 1 is known as the k-sparse largest eigenvalue.

If one takes k=p, the problem reduces to the ordinary PCA, and the optimal value becomes the largest eigenvalue of covariance matrix Σ.

After finding the optimal solution v, one deflates Σ to obtain a new matrix

Σ 1 = Σ ( v T Σ v ) v v T , {\displaystyle \Sigma _{1}=\Sigma -(v^{T}\Sigma v)vv^{T},}

and iterate this process to obtain further principal components. However, unlike PCA, sparse PCA cannot guarantee that different principal components are orthogonal. In order to achieve orthogonality, additional constraints must be enforced.

The following equivalent definition is in matrix form. Let V {\displaystyle V} be a p×p symmetric matrix, one can rewrite the sparse PCA problem as

max T r ( Σ V ) subject to T r ( V ) = 1 V 0 k 2 R a n k ( V ) = 1 , V 0. {\displaystyle {\begin{aligned}\max \quad &Tr(\Sigma V)\\{\text{subject to}}\quad &Tr(V)=1\\&\Vert V\Vert _{0}\leq k^{2}\\&Rank(V)=1,V\succeq 0.\end{aligned}}} Eq. 2

Tr is the matrix trace, and V 0 {\displaystyle \Vert V\Vert _{0}} represents the non-zero elements in matrix V. The last line specifies that V has matrix rank one and is positive semidefinite. The last line means that one has V = v v T {\displaystyle V=vv^{T}} , so Eq. 2 is equivalent to Eq. 1.

Moreover, the rank constraint in this formulation is actually redundant, and therefore sparse PCA can be cast as the following mixed-integer semidefinite program[1]

max T r ( Σ V ) subject to T r ( V ) = 1 | V i , i | z i , i { 1 , . . . , p } , | V i , j | 1 2 z i , i , j { 1 , . . . , p } : i j , V 0 , z { 0 , 1 } p , i z i k {\displaystyle {\begin{aligned}\max \quad &Tr(\Sigma V)\\{\text{subject to}}\quad &Tr(V)=1\\&\vert V_{i,i}\vert \leq z_{i},\forall i\in \{1,...,p\},\vert V_{i,j}\vert \leq {\frac {1}{2}}z_{i},\forall i,j\in \{1,...,p\}:i\neq j,\\&V\succeq 0,z\in \{0,1\}^{p},\sum _{i}z_{i}\leq k\end{aligned}}} Eq. 3

Because of the cardinality constraint, the maximization problem is hard to solve exactly, especially when dimension p is high. In fact, the sparse PCA problem in Eq. 1 is NP-hard in the strong sense.[2]

Computational considerations

As most sparse problems, variable selection in SPCA is a computationally intractable nonconvex NP-hard problem,[3] therefore greedy sub-optimal algorithms are required to find solutions.

Algorithms for SPCA

Several alternative approaches (of Eq. 1) have been proposed, including

  • a regression framework,[4]
  • a penalized matrix decomposition framework,[5]
  • a convex relaxation/semidefinite programming framework,[6]
  • a generalized power method framework[7]
  • an alternating maximization framework[8]
  • forward-backward greedy search and exact methods using branch-and-bound techniques,[9]
  • a certifiably optimal branch-and-bound approach[10]
  • Bayesian formulation framework.[11]
  • A certifiably optimal mixed-integer semidefinite branch-and-cut approach [1]

The methodological and theoretical developments of Sparse PCA as well as its applications in scientific studies are recently reviewed in a survey paper.[12]

Notes on Semidefinite Programming Relaxation

It has been proposed that sparse PCA can be approximated by semidefinite programming (SDP).[6] If one drops the rank constraint and relaxes the cardinality constraint by a 1-norm convex constraint, one gets a semidefinite programming relaxation, which can be solved efficiently in polynomial time:

max T r ( Σ V ) subject to T r ( V ) = 1 1 T | V | 1 k V 0. {\displaystyle {\begin{aligned}\max \quad &Tr(\Sigma V)\\{\text{subject to}}\quad &Tr(V)=1\\&\mathbf {1} ^{T}|V|\mathbf {1} \leq k\\&V\succeq 0.\end{aligned}}} Eq. 3

In the second constraint, 1 {\displaystyle \mathbf {1} } is a p×1 vector of ones, and |V| is the matrix whose elements are the absolute values of the elements of V.

The optimal solution V {\displaystyle V} to the relaxed problem Eq. 3 is not guaranteed to have rank one. In that case, V {\displaystyle V} can be truncated to retain only the dominant eigenvector.

While the semidefinite program does not scale beyond n=300 covariates, it has been shown that a second-order cone relaxation of the semidefinite relaxation is almost as tight and successfully solves problems with n=1000s of covariates [13]

Applications

Financial Data Analysis

Suppose ordinary PCA is applied to a dataset where each input variable represents a different asset, it may generate principal components that are weighted combination of all the assets. In contrast, sparse PCA would produce principal components that are weighted combination of only a few input assets, so one can easily interpret its meaning. Furthermore, if one uses a trading strategy based on these principal components, fewer assets imply less transaction costs.

Biology

Consider a dataset where each input variable corresponds to a specific gene. Sparse PCA can produce a principal component that involves only a few genes, so researchers can focus on these specific genes for further analysis.

High-dimensional Hypothesis Testing

Contemporary datasets often have the number of input variables ( p {\displaystyle p} ) comparable with or even much larger than the number of samples ( n {\displaystyle n} ). It has been shown that if p / n {\displaystyle p/n} does not converge to zero, the classical PCA is not consistent. In other words, if we let k = p {\displaystyle k=p} in Eq. 1, then the optimal value does not converge to the largest eigenvalue of data population when the sample size n {\displaystyle n\rightarrow \infty } , and the optimal solution does not converge to the direction of maximum variance. But sparse PCA can retain consistency even if p n . {\displaystyle p\gg n.}

The k-sparse largest eigenvalue (the optimal value of Eq. 1) can be used to discriminate an isometric model, where every direction has the same variance, from a spiked covariance model in high-dimensional setting.[14] Consider a hypothesis test where the null hypothesis specifies that data X {\displaystyle X} are generated from a multivariate normal distribution with mean 0 and covariance equal to an identity matrix, and the alternative hypothesis specifies that data X {\displaystyle X} is generated from a spiked model with signal strength θ {\displaystyle \theta } :

H 0 : X N ( 0 , I p ) , H 1 : X N ( 0 , I p + θ v v T ) , {\displaystyle H_{0}:X\sim N(0,I_{p}),\quad H_{1}:X\sim N(0,I_{p}+\theta vv^{T}),}

where v R p {\displaystyle v\in \mathbb {R} ^{p}} has only k non-zero coordinates. The largest k-sparse eigenvalue can discriminate the two hypotheses if and only if θ > Θ ( k log ( p ) / n ) {\displaystyle \theta >\Theta ({\sqrt {k\log(p)/n}})} .

Since computing k-sparse eigenvalue is NP-hard, one can approximate it by the optimal value of semidefinite programming relaxation (Eq. 3). If that case, we can discriminate the two hypotheses if θ > Θ ( k 2 log ( p ) / n ) {\displaystyle \theta >\Theta ({\sqrt {k^{2}\log(p)/n}})} . The additional k {\displaystyle {\sqrt {k}}} term cannot be improved by any other polynomial time algorithm if the planted clique conjecture holds.

Software/source code

  • amanpg - R package for Sparse PCA using the Alternating Manifold Proximal Gradient Method [15]
  • elasticnet – R package for Sparse Estimation and Sparse PCA using Elastic-Nets [16]
  • epca – R package for exploratory principal component analysis for large-scale dataset, including sparse principal component analysis and sparse matrix approximation. [17]
  • nsprcomp - R package for sparse and/or non-negative PCA based on thresholded power iterations[18]
  • scikit-learn – Python library for machine learning which contains Sparse PCA and other techniques in the decomposition module.[19]


Notes

  1. ^ The term loadings is inappropriately used for these vectors, which should be called coefficients, instead. The naming comes from the term loadings used in factor analysis to design the values generating the common covariance matrix. Since in standard PCA the loadings are equal to the coefficients, the term loadings has been used for the coefficients. This is quite unfortunate, because in SPCA the coefficients are not equal to the loadings.


References

  1. ^ a b Dimitris Bertsimas; Ryan Cory-Wright; Jean Pauphilet (2020). "Solving Large-Scale Sparse PCA to Certifiable (Near) Optimality". arXiv:2005.05195 [math.OC].
  2. ^ Andreas M. Tillmann; Marc E. Pfetsch (2013). "The Computational Complexity of the Restricted Isometry Property, the Nullspace Property, and Related Concepts in Compressed Sensing". IEEE Transactions on Information Theory. 60 (2): 1248–1259. arXiv:1205.2081. CiteSeerX 10.1.1.760.2559. doi:10.1109/TIT.2013.2290112. S2CID 2788088.
  3. ^ Moghaddam, Baback; Weiss, Yair; Avidan, Shai (2007-09-07). Schölkopf, Bernhard; Platt, John; Hofmann, Thomas (eds.). Advances in Neural Information Processing Systems 19: Proceedings of the 2006 Conference. The MIT Press. pp. 915–922. doi:10.7551/mitpress/7503.001.0001. ISBN 978-0-262-25691-9.
  4. ^ Hui Zou; Trevor Hastie; Robert Tibshirani (2006). "Sparse principal component analysis" (PDF). Journal of Computational and Graphical Statistics. 15 (2): 262–286. CiteSeerX 10.1.1.62.580. doi:10.1198/106186006x113430. S2CID 5730904.
  5. ^ Fan Chen; Karl Rohe (2021). "A New Basis for Sparse Principal Component Analysis". Journal of Computational and Graphical Statistics. arXiv:2007.00596. doi:10.1080/10618600.2023.2256502.
  6. ^ a b Alexandre d’Aspremont; Laurent El Ghaoui; Michael I. Jordan; Gert R. G. Lanckriet (2007). "A Direct Formulation for Sparse PCA Using Semidefinite Programming" (PDF). SIAM Review. 49 (3): 434–448. arXiv:cs/0406021. doi:10.1137/050645506. S2CID 5490061.
  7. ^ Michel Journee; Yurii Nesterov; Peter Richtarik; Rodolphe Sepulchre (2010). "Generalized Power Method for Sparse Principal Component Analysis" (PDF). Journal of Machine Learning Research. 11: 517–553. arXiv:0811.4724. Bibcode:2008arXiv0811.4724J. CORE Discussion Paper 2008/70.
  8. ^ Peter Richtarik; Majid Jahani; S. Damla Ahipasaoglu; Martin Takac (2021). "Alternating Maximization: Unifying Framework for 8 Sparse PCA Formulations and Efficient Parallel Codes". Optimization and Engineering. 22 (3): 1493–1519. arXiv:1212.4137. doi:10.1007/s11081-020-09562-3. S2CID 2549610.
  9. ^ Baback Moghaddam; Yair Weiss; Shai Avidan (2005). "Spectral Bounds for Sparse PCA: Exact and Greedy Algorithms" (PDF). Advances in Neural Information Processing Systems. Vol. 18. MIT Press.
  10. ^ Lauren Berk; Dimitris Bertsimas (2019). "Certifiably optimal sparse principal component analysis". Mathematical Programming Computation. 11 (3). Springer: 381–420. doi:10.1007/s12532-018-0153-6. hdl:1721.1/131566. S2CID 126998398.
  11. ^ Yue Guan; Jennifer Dy (2009). "Sparse Probabilistic Principal Component Analysis" (PDF). Journal of Machine Learning Research Workshop and Conference Proceedings. 5: 185.
  12. ^ Hui Zou; Lingzhou Xue (2018). "A Selective Overview of Sparse Principal Component Analysis". Proceedings of the IEEE. 106 (8): 1311–1320. doi:10.1109/jproc.2018.2846588.
  13. ^ Dimitris Bertsimas; Ryan Cory-Wright (2020). "On polyhedral and second-order cone decompositions of semidefinite optimization problems". Operations Research Letters. 48 (1). Elsevier: 78–85. arXiv:1910.03143. doi:10.1016/j.orl.2019.12.003.
  14. ^ Quentin Berthet; Philippe Rigollet (2013). "Optimal Detection of Sparse Principal Components in High Dimension". Annals of Statistics. 41 (1): 1780–1815. arXiv:1202.5070. doi:10.1214/13-aos1127. S2CID 7162068.
  15. ^ [1] https://cran.r-project.org/web/packages/amanpg/index.html
  16. ^ [2] https://cran.r-project.org/web/packages/elasticnet/index.html
  17. ^ [3] https://cran.r-project.org/web/packages/epca/index.html
  18. ^ [4] https://cran.r-project.org/web/packages/nsprcomp/index.html
  19. ^ [5] http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.SparsePCA.html