Sparse PCA via Bipartite Matchings

Introduction

Principal Component Analysis (PCA) reduces the dimensionality of a data set by projecting it onto principal subspaces spanned by the leading eigenvectors of the sample covariance matrix. Sparse PCA is a useful variant that offers higher data interpretability , a property that is sometimes desired even at the cost of statistical fidelity . Furthermore, when the obtained features are used in subsequent learning tasks, sparsity potentially leads to better generalization error .

Given a real ${{n}\times d}$ data matrix $\mathbf{S}$ representing ${n}$ centered data points supported on $d$ features, the leading sparse principal component of the data set is the sparse vector that maximizes the explained variance:

where $\mathbf{A}=\nicefrac{{1}}{{{n}}}\cdot\mathbf{S}^{{{\top}}}\mathbf{S}$ is the $d\times d$ empirical covariance matrix. The sparsity constraint makes the problem NP-hard and hence computationally intractable in general, and hard to approximate within some small constant . A significant volume of prior work has focused on algorithms that approximately solve the optimization problem , while a large volume of theoretical results has been established under planted statistical models .

In most practical settings, we tend to go beyond computing a single sparse PC. Contrary to the single-component problem, there has been limited work on computing multiple components. The scarcity is partially attributed to conventional PCA wisdom: multiple components can be computed one-by-one, repeatedly, by solving the single-component sparse PCA problem (1) and deflating the input data to remove information captured by previously extracted components . In fact, the multi-component version of sparse PCA is not uniquely defined in the literature. Different deflation-based approaches can lead to different outputs: extracted components may or may not be orthogonal, while they may have disjoint or overlapping supports . In the statistics literature, where the objective is typically to recover a “true” principal subspace, a branch of work has focused on the “subspace row sparsity” , an assumption that leads to sparse components all supported on the same set of variables. While in , the authors discuss an alternative perspective on the fundamental objective of the sparse PCA problem.

In this work, we develop a novel algorithm for the multi-component sparse PCA problem with disjoint supports. Formally, we are interested in finding $k$ components that are $s$-sparse, have disjoint supports, and jointly maximize the explained variance:

with $\mathbf{X}^{j}$ denoting the $j$th column of $\mathbf{X}$. The number $k$ of the desired components is a user defined parameter and we consider it to be a small constant.

Contrary to the greedy sequential approach that repeatedly uses deflation, our algorithm jointly computes all the vectors in $\mathbf{X}$, and comes with theoretical approximation guarantees. We note that even if one could solve each single-component sparse PCA problem (1) exactly, greedy deflation can be highly suboptimal. We show this through a simple example in Section 7.

We develop an algorithm that provably approximates the solution to the sparse PCA problem (2) within a multiplicative factor arbitrarily close to $1$. To the best of our knowledge, this is the first algorithm that jointly optimizes multiple components with disjoint supports, provably. Our algorithm is combinatorial; it recasts sparse PCA as multiple instances of bipartite maximum weight matching on graphs determined by the input data.

The computational complexity of our algorithm grows as a low order polynomial in the ambient dimension $d$, but is exponential in the intrinsic dimension of the input data, i.e., the rank of $\mathbf{A}$. To alleviate the impact of this dependence, our algorithm can be applied on a low-dimensional sketch of the input data to obtain an approximate solution to (2). This extra level of approximation introduces an additional penalty in our theoretical approximation guarantees, which naturally depends on the quality of the sketch and, in turn, the spectral decay of $\mathbf{A}$. We show how these bounds further translate to an additive PTAS (polynomial-time approximation scheme) for sparse PCA. Our additive PTAS outputs an approximate solution with explained variance of at least ${{\textsf{{\mbox{OPT}}}}}-\epsilon\cdot s$, for any sparsity $s\in\{1,\ldots,n\}$, any constant error $\epsilon>0$ and any $k=O(1)$ number of orthogonal components.Here, OPT is the explained variance captured by the optimal set of $k$ components that are $s$ sparse and have disjoint supports.

We empirically evaluate our algorithm on real datasets, and compare it against state-of-the-art methods for the single-component sparse PCA problem (1) in conjunction with the appropriate deflation step. In many cases, our algorithm—as a result of jointly optimizing over multiple components—leads to significantly improved results, and outperforms deflation-based approaches.

Sparse PCA through Bipartite Matchings

Our algorithm approximately solves the constrained maximization (2) on a $d\times d$ rank-$r$ PSD matrix ${\mathbf{A}}$ within a multiplicative factor arbitrarily close to $1$. It operates by recasting the maximization into multiple instances of the bipartite maximum weight matching problem. Each instance ultimately yields a feasible solution: a set of $k$ components that are $s$-sparse and have disjoint supports. The algorithm examines these solutions, and outputs the one that maximizes the explained variance, i.e., the quadratic objective in (2).

The computational complexity of our algorithm grows as a low order polynomial in the ambient dimension $d$ of the input, but exponentially in its rank $r$. Despite the unfavorable dependence on the rank, it is unlikely that a substantial improvement can be achieved in general . However, decoupling the dependence on the ambient and the intrinsic dimension of the input has an interesting ramification; instead of the original input ${\mathbf{A}}$, our algorithm can be applied on a low-rank surrogate to obtain an approximate solution, alleviating the dependence on $r$. We discuss this in Section 3, and present the approximation bound that this allows us to obtain.

Under the variational characterization of the trace objective in (4), the sparse PCA problem (2) can be re-written as a joint maximization over the variables $\mathbf{X}$ and $\mathbf{C}$ as follows:

The alternative formulation of the sparse PCA problem in (5) takes a step towards decoupling the dependence of the optimization on the ambient and intrinsic dimensions $d$ and ${r}$, respectively. The motivation behind the introduction of the auxiliary variable $\mathbf{C}$ will become clear in the sequel.

For a given $\mathbf{C}$, the value of $\mathbf{X}\in\mathcal{X}_{k}$ that maximizes the objective in (5) for that $\mathbf{C}$ is

where $\mathbf{W}{\triangleq}{\mathbf{U}}{\mathbf{\Lambda}}^{1/2}\mathbf{C}$ is a real $d\times{k}$ matrix. The constrained, non-convex maximization (6) plays a central role in our developments. We will later describe a combinatorial ${O\mathopen{}(d\cdot({s}\cdot{k})^{2})}$ procedure to efficiently compute $\widehat{\mathbf{X}}$ , reducing the maximization to an instance of the bipartite maximum weight matching problem. For now, however, let us assume that such a procedure exists.

For any real $d\times d$ rank-$r$ PSD matrix ${\mathbf{A}}$, desired number of components ${k}$, number ${s}$ of nonzero entries per component, and accuracy parameter $\epsilon\in(0,1)$, Algorithm 1 outputs $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\in\mathcal{X}_{k}$ such that

where ${\mathbf{X}}_{\star}{\triangleq}\operatorname*{arg\,max}_{\mathbf{X}\in\mathcal{X}_{k}}\text{{Tr}}\mathopen{}\left(\mathbf{X}^{{\top}}{{\mathbf{A}}}\mathbf{X}\right),$ in time T_{\texttt{SVD}}(r)+O\mathopen{}\bigl{(}\bigl{(}\tfrac{4}{\epsilon}\bigr{)}^{r\cdot{{k}}}\cdot d\cdot({s}\cdot{k})^{2}\bigr{)}.

Algorithm 1 is the first nontrivial algorithm that provably approximates the solution of the sparse PCA problem (2). According to Theorem 1, it achieves an objective value that lies within a multiplicative factor from the optimal, arbitrarily close to $1$. Its complexity grows as a low-order polynomial in the dimension $d$ of the input, but exponentially in the intrinsic dimension $r$. Note, however, that it can be exponentially faster compared to the ${O(d^{s\cdot k})}$ brute force approach that exhaustively considers all candidate supports for the $k$ sparse components. The complexity of our algorithm follows from the cardinality of the net and the complexity of Algorithm 2, the subroutine that solves the constrained maximization (6). The latter is a key ingredient of our algorithm, and is discussed in detail in the next subsection. A formal proof of Theorem 1 is provided in Section 9.2.

In the core of Algorithm 1 lies Algorithm 2, a procedure that solves the constrained maximization in (6). The algorithm breaks down the maximization into two stages. First, it identifies the support of the optimal solution $\widehat{\mathbf{X}}$ . Determining the support reduces to an instance of the maximum matching problem on a weighted bipartite graph $G$. Then, it recovers the exact values of the nonzero entries in $\widehat{\mathbf{X}}$ based on the Cauchy-Schwarz inequality. In the sequel, we provide a brief description of Algorithm 2, leading up to its guarantees in Lemma 2.27.

Let $\mathcal{I}_{j}{\triangleq}\text{supp}(\widehat{\mathbf{X}}^{j})$ be the support of the $j$th column of $\widehat{\mathbf{X}}$, ${j=1,\ldots,{k}}$. The objective in (6) becomes

The last inequality is an application of the Cauchy-Schwarz Inequality and the constraint ${\|\mathbf{X}^{j}\|_{2}=1}$ $\forall\,j\in\{{1,\ldots,{k}}\}$. In fact, if an oracle reveals the supports $\mathcal{I}_{j}$, ${j=1,\ldots,{k}}$, the upper bound in (7) can always be achieved by setting the nonzero entries of $\widehat{\mathbf{X}}$ as in Algorithm 2 (Line $6$). Therefore, the key in solving (6) is determining the collection of supports to maximize the right-hand side of (7).

By constraint, the sets $\mathcal{I}_{j}$ must be pairwise disjoint, each with cardinality ${s}$. Consider a weighted bipartite graph {G=\bigl{(}{U=\{U_{1},\ldots,U_{k}\}},V,E\bigr{)}} constructed as followsThe construction is formally outlined in Algorithm 4 in Section 8. (Fig. 1):

$V$ is a set of $d$ vertices $v_{1},\ldots,v_{d}$, corresponding to the $d$ variables, i.e., the $d$ rows of $\widehat{\mathbf{X}}$ .

$U$ is a set of ${k}\cdot{s}$ vertices, conceptually partitioned into ${k}$ disjoint subsets $U_{1},\ldots,U_{{k}}$, each of cardinality ${s}$. The $j$th subset, $U_{j}$, is associated with the support $\mathcal{I}_{j}$; the ${s}$ vertices $u^{{(j)}}_{\alpha}$, ${\alpha=1,\ldots,s}$ in $U_{j}$ serve as placeholders for the variables/indices in $\mathcal{I}_{j}$.

Finally, the edge set is ${E=U\times V}$. The edge weights are determined by the $d\times{k}$ matrix $\mathbf{W}$ in (6). In particular, the weight of edge $(u^{{(j)}}_{\alpha},v_{i})$ is equal to $W_{ij}^{2}$. Note that all vertices in $U_{j}$ are effectively identical; they all share a common neighborhood and edge weights.

Any feasible support $\{\mathcal{I}_{j}\}_{j=1}^{k}$ corresponds to a perfect matching in $G$ and vice-versa. Recall that a matching is a subset of the edges containing no two edges incident to the same vertex, while a perfect matching, in the case of an unbalanced bipartite graph ${G=(U,V,E)}$ with $|U|\leq|V|$, is a matching that contains at least one incident edge for each vertex in $U$. Given a perfect matching $\mathcal{M}\subseteq E$, the disjoint neighborhoods of $U_{j}$s under $\mathcal{M}$ yield a support $\{\mathcal{I}_{j}\}_{j=1}^{k}$ . Conversely, any valid support yields a unique perfect matching in $G$ (taking into account that all vertices in $U_{j}$ are isomorphic). Moreover, due to the choice of weights in $G$, the right-hand side of (7) for a given support $\{\mathcal{I}_{j}\}_{j=1}^{k}$ is equal to the weight of the matching $\mathcal{M}$ in $G$ induced by the former, i.e., $\sum_{j=1}^{{k}}\sum_{i\in\mathcal{I}_{j}}W_{ij}^{2}$ $=$ $\sum_{(u,v)\in\mathcal{M}}w(u,v)$ . It follows that determining the support of the solution in (6), reduces to solving the maximum weight matching problem on the bipartite graph $G$.

A more formal analysis and proof of Lemma 2.27 is available in Section 9.1. With Algorithm 2 and Lemma 2.27 in place, we complete the description of our sparse PCA algorithm (Algorithm 1) and the proof sketch of Theorem 1.

Sparse PCA on Low-Dimensional Sketches

Algorithm 1 approximately solves the sparse PCA problem (2) on a $d\times d$ rank-$r$ PSD matrix $\mathbf{A}$, in time that grows as a low-order polynomial in the ambient dimension $d$, but depends exponentially on $r$. This dependence can be prohibitive in practice. To mitigate its effect, instead of the original input, we can apply our sparse PCA algorithm on a low-rank approximation of $\mathbf{A}$. Intuitively, the quality of the extracted components should depend on how well that low-rank surrogate approximates the original input.

More formally, let $\mathbf{S}$ be the real ${n}\times d$ data matrix representing $n$ (potentially centered) datapoints in $d$ variables, and $\mathbf{A}$ the corresponding $d\times d$ covariance matrix. Further, let $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ be a low-dimensional sketch of the original data; an ${n}\times d$ matrix whose rows lie in an $r$-dimensional subspace, with $r$ being an accuracy parameter. Such a sketch can be obtained in several ways, including for example exact or approximate SVD, or online sketching methods . Finally, let \macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}=\scalebox{0.85}{\nicefrac{{1}}{{n}}\cdot{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}^{{{\top}}}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}} be the covariance matrix of the sketched data. Then, instead of $\mathbf{A}$, we can approximately solve the sparse PCA problem by applying Algorithm 1 on the low-rank surrogate $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$. The above are formally outlined in Algorithm 3. We note that the covariance matrix $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ does not need to be explicitly computed; Algorithm 1 can operate directly on the (sketched) input data matrix.

For any $n\times d$ input data matrix $\mathbf{S}$, with corresponding empirical covariance matrix $\mathbf{A}=\nicefrac{{1}}{{n}}\cdot\mathbf{S}^{{{\top}}}\mathbf{S}$, any desired number of components $k$, and accuracy parameters $\epsilon\in(0,1)$ and $r$, Algorithm 3 outputs ${\mathbf{X}_{{(r)}}\in\mathcal{X}_{k}}$ such that

in time T_{\texttt{SKETCH}}(r)+T_{\texttt{SVD}}(r)+O\mathopen{}\bigl{(}\bigl{(}\tfrac{4}{\epsilon}\bigr{)}^{r\cdot{{k}}}\cdot{d}\cdot({s}\cdot{k})^{2}\bigr{)}. Here, $\mathbf{X}_{\star}{\triangleq}\operatorname*{arg\,max}_{\mathbf{X}\in\mathcal{X}_{k}}\text{{Tr}}\mathopen{}\left(\mathbf{X}^{{{\top}}}\mathbf{A}\mathbf{X}\right)$, and $\lambda_{1,s}(\mathbf{A})$ denotes the sparse eigenvalue, i.e., the eigenvalue that corresponds to the principal $s$-sparse eigenvector of $\mathbf{A}$.

The error $\lambda_{1,s}(\mathbf{A}-\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{})$ and in turn the tightness of the approximation guarantees hinges on the quality of the sketch $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$. Higher values of the parameter $r$ (the rank of the sketch) can allow for a more accurate solution and tighter guarantees. That is the case, for example, when the sketch is obtained through exact SVD. In that sense, Theorem 2 establishes a natural trade-off between the running time of Algorithm 3 and the quality of the approximation guarantees. A formal proof of Theorem 2 is provided in Section 9.3. Observe that the error term itself is a sparse eigenvalue that is hard to approximate, however even loose bounds provide tight conditional approximation results, as we see next.

Using the main matrix approximation result of , the next theorem establishes that Algorithm 3 can be turned into an additive PTAS.

Let $\mathbf{A}$ be a $d\times d$ positive semidefinite matrix with entries in $$,$\mathbf{V}$be a$d\times d$matrix such that$\mathbf{A}={\bf V}{\bf V}^{\top}$. Further, let$\mathbf{R}$be a random$d\times r$matrix with entries drawn i.i.d. according to$\mathcal{N}(0,1/r)$, and define

For any constant $\epsilon\in(0,1]$, let $r=O(\epsilon^{-2}\log d)$. Then, for any desired sparsity $s$, and number of components $k={O}(1)$, Algorithm 1 with input argument $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ and accuracy parameter $\epsilon$, outputs ${\mathbf{X}_{{(r)}}\in\mathcal{X}_{k}}$ such that

with probability at least $1-1/\text{poly}(d)$, in time $n^{O(\log(1/\epsilon)/\epsilon^{2}))}$.

Note that $\lambda_{1}(\mathbf{A}-\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{})$ serves as another elementary upper bound on $\lambda_{1,s}(\mathbf{A}-\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{})$. If $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ is a the rank-$d$ SVD approximation of $\mathbf{A}$, then—similar to —we can obtain a multiplicative PTAS for sparse PCA, under the assumption of a decaying spectrum (e.g., under a power-law decay), and for $s=\Omega(n)$.

Related Work

We are not aware of any algorithm with provable guarantees for sparse PCA with disjoint supports. Multiple components can be extracted by repeatedly solving (1) using one of the aforementioned methods. To ensure disjoint supports, variables “selected” by a component are removed from the dataset. This greedy approach, however, can result in highly suboptimal objective value (See example in Sec. 7).

Parallel to the algorithmic and optimization perspective, there is large line of statistical analysis for sparse PCA that focuses on guarantees pertaining to planted models and the recovery of a “true” sparse component .

There has been some work on the explicit estimation of principal subspaces or multiple components under sparsity constraints. Non-deflation-based algorithms include extensions of the diagonal thresholding algorithm and iterative thresholding approaches , while and propose methods that rely on the “row sparsity for subspaces” assumption of . These methods yield components supported on a common set of variables, and hence solve a problem different from (2). Magdon-Ismail and Boutsidis discuss the multiple component Sparse PCA problem, propose an alternative objective function and for that problem obtain interesting theoretical guarantees. Finally, develops a framework for sparse matrix factorizaiton problems, based on a novel atomic norm. That framework captures sparse PCA – although not explicitly the constraint of disjoint supports – but the resulting optimization problem, albeit convex, is NP-hard.

Experiments

We evaluate our algorithm on a series of real datasets, and compare it to deflation-based approaches for sparse PCA using TPower , EM , and SpanSPCA . The latter are representative of the state of the art for the single-component sparse PCA problem (1). Multiple components are computed one by one. To ensure disjoint supports, the deflation step effectively amounts to removing from the dataset all variables used by previously extracted components. For algorithms that are randomly initialized, we depict best results over multiple random restarts. Additional experimental results are listed in Section 11 of the appendix.

Our experiments are conducted in a Matlab environment. Due to its nature, our algorithm is easily parallelizable; its prototypical implementation utilizes the Parallel Pool Matlab feature to exploit multicore (or distributed cluster) capabilities. Recall that our algorithm operates on a low-rank approximation of the input data. Unless otherwise specified, it is configured for a rank-$4$ approximation obtained via truncated SVD. Finally, we put a time barrier in the execution of our algorithm, at the cost of the theoretical approximation guarantees; the algorithm returns best results at the time of termination. This “early termination” can only hurt the performance of our algorithm.

Leukemia Dataset. We evaluate our algorithm on the Leukemia dataset . The dataset comprises $72$ samples, each consisting of expression values for $12582$ probe sets. We extract ${k=5}$ sparse components, each active on $s=50$ features. In Fig. 2(a), we plot the cumulative explained variance versus the number of components. Deflation-based approaches are greedy: the leading components capture high values of variance, but subsequent ones contribute less. On the contrary, our algorithm jointly optimizes the ${k=5}$ components and achieves higher total cumulative variance; one cannot identify a top component. We repeat the experiment for multiple values of $k$. Fig. 2(b) depicts the total cumulative variance capture by each method, for each value of $k$.

Additional Datasets. We repeat the experiment on multiple datasets, arbitrarily selected from . Table 1 lists the total cumulative variance captured by ${k=5}$ components, each with ${s=40}$ nonzero entries, extracted using the four methods. Our algorithm achieves the highest values in most cases.

Bag of Words (BoW) Dataset. This is a collection of text corpora stored under the “bag-of-words” model. For each text corpus, a vocabulary of ${d}$ words is extracted upon tokenization, and the removal of stopwords and words appearing fewer than ten times in total. Each document is then represented as a vector in that $d$-dimensional space, with the $i$th entry corresponding to the number of appearances of the $i$th vocabulary entry in the document.

We solve the sparse PCA problem (2) on the word-by-word cooccurrence matrix, and extract ${k=8}$ sparse components, each with cardinality ${s=10}$. We note that the latter is not explicitly constructed; our algorithm can operate directly on the input word-by-document matrix. Table 2 lists the variance captured by each method; our algorithm consistently outperforms the other approaches.

Finally, note that here each sparse component effectively selects a small set of words. In turn, the $k$ extracted components can be interpreted as a set of well-separated topics. In Table 3, we list the topics extracted from the NY Times corpus (part of the Bag of Words dataset). The corpus consists of ${3\cdot 10^{5}}$ news articles and a vocabulary of ${d=102660}$ words.

Conclusions

We considered the sparse PCA problem for multiple components with disjoint supports. Existing methods for the single component problem can be used along with an appropriate deflation step to compute multiple components one by one, leading to potentially suboptimal results. We presented a novel algorithm for jointly optimizing multiple sparse and disjoint components with provable approximation guarantees. Our algorithm is combinatorial and exploits interesting connections between the sparse PCA and the bipartite maximum weight matching problems. It runs in time that grows as a low-order polynomial in the ambient dimension of the input data, but depends exponentially on its rank. To alleviate this dependency, we can apply the algorithm on a low-dimensional sketch of the input, at the cost of an additional error in our theoretical approximation guarantees. Empirical evaluation of our algorithm demonstrated that in many cases it outperforms deflation-based approaches.

Acknowledgments

DP is generously supported by NSF awards CCF-1217058 and CCF-1116404 and MURI AFOSR grant 556016. This research has been supported by NSF Grants CCF 1344179, 1344364, 1407278, 1422549 and ARO YIP W911NF-14-1-0258.

References

Supplemental Material

On the sub-optimality of deflation – An example

We provide a simple example demonstrating the sub-optimality of deflation based approaches for computing multiple sparse components with disjoint supports. Consider the real $4\times 4$ matrix

with $\epsilon,\delta>0$ such that ${\epsilon+\delta<1}$. Note that $\mathbf{A}$ is PSD; $\mathbf{A}=\mathbf{B}^{{\top}}\mathbf{B}$ for

We seek two $2$-sparse components with disjoint supports, i.e., the solution to

Iterative computation with deflation. Following an iterative, greedy procedure with a deflation step, we compute one component at the time. The first component is

Recall that for any unit norm vector $\mathbf{x}$ with support $I=\text{supp}(\mathbf{x})$,

where $\mathbf{A}_{I,I}$ denotes the principal submatrix of $\mathbf{A}$ formed by the rows and columns indexed by $I$. Equality can be achieved in (10) for $\mathbf{x}$ equal to the leading eigenvector of $\mathbf{A}_{I,I}$. Hence, it suffices to determine the optimal support for $\mathbf{x}_{1}$. Due to the small size of the example, it is easy to determine that the set $I_{1}=\{1,4\}$ maximizes the objective in (10) over all sets of two indices, achieving value

Since subsequent components must have disjoint supports, it follows that the support of the second $2$-sparse component $\mathbf{x}_{2}$ is $I_{2}=\{2,3\}$, and $\mathbf{x}_{2}$ achieves value

In total, the objective value in (8) achieved by the greedy computation with a deflation step is

The sub-optimality of deflation. Consider an alternative pair of $2$-sparse components ${\mathbf{x}_{1}^{\prime}}$ and ${\mathbf{x}_{2}^{\prime}}$ with support sets ${I_{1}^{\prime}=\{1,2\}}$ and ${I_{2}^{\prime}=\{3,4\}}$, respectively. Based on the above, such a pair achieves objective value in (8) equal to

which clearly outperforms the objective value in (13) (under the assumption $\epsilon+\delta<1$), demonstrating the sub-optimality of the $\mathbf{x}_{1}$, $\mathbf{x}_{2}$ pair computed by the deflation-based approach. In fact, for small $\epsilon,\delta$ the objective value in the second case is larger than the former by almost a factor of two.

Construction of Bipartite Graph

The following algorithm formally outlines the steps for generating the bipartite graph {G=\bigl{(}\{{U_{j}}\}_{j=1}^{{k}},V,E\bigr{)}} given a weight $d\times{k}$ matrix $\mathbf{W}$.

Proofs

For any real $d\times{k}$ matrix $\mathbf{W}$, and Algorithm 2 outputs

in time $O\mathopen{}\left(d\cdot({s}\cdot{k})^{2}\right)$.

Consider a matrix $\mathbf{X}\in\mathcal{X}_{k}$ and let $I_{j}$, $j=1\ldots,{k}$ denote the support sets of its columns. By the constraints in $\mathcal{X}_{k}$, those sets are disjoint, i.e., $I_{j_{1}}\cap I_{j_{2}}=\emptyset$ $\forall j_{1},j_{2}\in\{{1,\ldots,{k}}\},j_{1}\neq j_{2}$, and

The last inequality is due to Cauchy-Schwarz and the fact that $\|\mathbf{X}^{j}\|_{2}\leq 1$, $\forall\,j\in\{{1,\ldots,{k}}\}$. In fact, if the supports sets $I_{j}$, ${j=1,\ldots,{k}}$ were known, the upper bound in (15) would be achieved by setting ${\mathbf{X}}_{I_{j}}^{j}=\mathbf{W}_{I_{j}}^{j}/\|\mathbf{W}_{I_{j}}^{j}\|_{2}$, i.e., setting the nonzero subvector of the $j$th column of ${\mathbf{X}}$ colinear to the corresponding subvector of the $j$th column of $\mathbf{W}$. Hence, the key step towards computing the optimal solution $\widetilde{\mathbf{X}}$ is to determine the support sets $I_{j}$, $j=1,\ldots,{k}$ of its columns.

The set represents all possible supports for the members of $\mathcal{X}_{k}$. Taking into account the previous discussion, the maximization in (14) can be written with respect to ${\mathbf{Z}\in\mathcal{Z}}$:

$V$ is a set of $d$ vertices $v_{1},\ldots,v_{d}$, corresponding to the $d$ variables, i.e., the $d$ rows of $\widehat{\mathbf{X}}$ .

$U$ is a set of ${k}\cdot{s}$ vertices, conceptually partitioned into ${k}$ disjoint subsets $U_{1},\ldots,U_{{k}}$, each of cardinality ${s}$. The $j$th subset, $U_{j}$, is associated with the support $\mathcal{I}_{j}$; the ${s}$ vertices $u^{{(j)}}_{\alpha}$, ${\alpha=1,\ldots,s}$ in $U_{j}$ serve as placeholders for the variables/indices in $\mathcal{I}_{j}$.

Finally, the edge set is ${E=U\times V}$. The edge weights are determined by the $d\times{k}$ matrix $\mathbf{W}$ in (6). In particular, the weight of edge $(u^{{(j)}}_{\alpha},v_{i})$ is equal to $W_{ij}^{2}$. Note that all vertices in $U_{j}$ are effectively identical; they all share a common neighborhood and edge weights.

It is straightforward to verify that any $\mathbf{Z}\in\mathcal{Z}$ corresponds to a perfect matching in $G$ and vice versa; $Z_{ij}=1$ if and only if vertex $v_{i}\in V$ is matched with a vertex in $U_{j}$ (all vertices in $U_{j}$ are equivalent with respect to their neighborhood). Further, the objective value in (16) for a given $\mathbf{Z}\in\mathcal{Z}$ is equal to the weight of the corresponding matching in $G$. More formally,

Given a perfect matching $\mathcal{M}$, the support $I_{j}$ of the $j$th column of $\mathbf{Z}$ is determined by the neighborhood of $U_{j}$ in the matching:

Note that the sets $I_{j}$, $j=1,\ldots,{k}$ are indeed disjoint, and each has cardinality equal to ${s}$. The weight of the matching $\mathcal{M}$ is

which is equal to the objective function in (16).

Conversely, given an indicator matrix $\mathbf{Z}\in\mathcal{Z}$, let $I_{j}{\triangleq}\text{supp}(\mathbf{Z}^{j})$, and let $I_{j}(\alpha)$ denote the $\alpha$th element in the set, $\alpha=1,\ldots,{s}$ (with an arbitrary ordering). Then,

is a perfect matching in $G$. The objective value achieved by $\mathbf{Z}$ is equal to the weight of $\mathcal{M}$:

It follows from (18) and (19) that to determine $\widetilde{\mathbf{Z}}$, it suffices to compute a maximum weight perfect matching in $G$. The desired support is then obtained as described in (17) (lines 4-7 of Algorithm 2). This complete the proof of correctness of Algorithm 2 which proceeds in the steps described above to determine the support of $\widetilde{\mathbf{X}}$.

The weighted bipartite graph $G$ is generated in $O(d\cdot({s}\cdot{k}))$. The running time of Algorithm 2 is dominated by computing the maximum weight matching of $G$. For the case of unbalanced bipartite graph with $|U|={{s}\cdot{k}}<d=\lvert{V}\rvert$ the Hungarian algorithm can be modified to compute the maximum weight bipartite matching in time O\mathopen{}\left(|E||U|+|U|^{2}\log{|U|}\right)={O\mathopen{}\bigl{(}d\cdot({s}\cdot{k})^{2}\bigr{)}}. This completes the proof. ∎

2 Guarantees of Algorithm 1 – Proof of Theorem 1

We first prove a more general version of Theorem 1 for arbitrary constraint sets. Combining that with the guarantees of Algorithm 2, we prove the Theorem 1.

in time T_{\texttt{SVD}}(r)+O\mathopen{}\bigl{(}\bigl{(}\tfrac{4}{\epsilon}\bigr{)}^{r\cdot{{k}}}\cdot\bigl{(}T_{\mathcal{X}}+{{k}}{d}\bigr{)}\bigr{)}, where $T_{\mathcal{X}}$ is the time required to compute $P_{\mathcal{X}}(\cdot)$ and $T_{\texttt{SVD}}(r)$ the time required to compute the truncated SVD of $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$.

In fact, equality in (20) is achieved for $\mathbf{c}$ colinear to ${\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}^{1/2}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\mathbf{x}$, and hence,

Recall that $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}$ is the optimal solution of the trace maximization on $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$, i.e.,

Let $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}$ be the maximizing value of $\mathbf{C}$ in (22) for $\mathbf{X}=\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}$, i.e., $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}$ is an $r\times{k}$ matrix with unit-norm columns such that for all $j\in\{1,\ldots,{k}\}$,

Based on the above, for all $j\in\{1,\ldots,{k}\}$,

The first step follows by the definition of $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}$, the second by the linearity of the inner product, the third by the triangle inequality, the fourth by Cauchy-Schwarz inequality and the last by (24). Rearranging the terms in (25),

Summing the terms in (26) over all $j\in\{1,\ldots,{k}\}$,

Let $\mathbf{X}_{\sharp}\in\mathcal{X}$ be the candidate solution produced by the algorithm at $\mathbf{C}_{\sharp}$, i.e.,

where $(\alpha)$ follows from the observation in (22), $(\beta)$ from the sub-optimality of $\mathbf{C}_{\sharp}$, $(\gamma)$ by the definition of $\mathbf{X}_{\sharp}$ in (28), while $(\delta)$ follows from (27). According to (29), at least one of the candidate solutions produced by Algorithm 1, namely $\mathbf{X}_{\sharp}$, achieves an objective value within a multiplicative factor $(1-\epsilon)$ from the optimal, implying the guarantees of the lemma.

Finally, the running time of Algorithm 1 follows immediately from the cost per iteration and the cardinality of the $\nicefrac{{\epsilon}}{{2}}$-net on the unit-sphere. Note that matrix multiplications can exploit the singular value decomposition which is performed once. ∎

For any real $d\times d$ rank-$r$ PSD matrix $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$, desired number of components ${k}$, number ${s}$ of nonzero entries per component, and accuracy parameter $\epsilon\in(0,1)$, Algorithm 1 outputs $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\in\mathcal{X}_{k}$ such that

where ${\mathbf{X}}_{\star}{\triangleq}\operatorname*{arg\,max}_{\mathbf{X}\in\mathcal{X}_{k}}\text{{Tr}}\mathopen{}\left(\mathbf{X}^{{\top}}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}\mathbf{X}\right),$ in time T_{\texttt{SVD}}(r)+O\mathopen{}\bigl{(}\bigl{(}\tfrac{4}{\epsilon}\bigr{)}^{r\cdot{{k}}}\cdot d\cdot({s}\cdot{k})^{2}\bigr{)}. $T_{\texttt{SVD}}(r)$ is the time required to compute the truncated SVD of $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$.

3 Guarantees of Algorithm 3 – Proof of Theorem 2

We prove Theorem 2 with the approximation guarantees of Algorithm 3.

Then, for any $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\in\mathcal{X}$ such that \text{{Tr}}\mathopen{}\bigl{(}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}^{{\top}}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\bigr{)}\geq\gamma\cdot\text{{Tr}}\mathopen{}\bigl{(}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}^{{\top}}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}\bigr{)} for some $0<\gamma<1$,

By the optimality of $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}$ for ${\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}$,

In turn, for any $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\in\mathcal{X}$ such that $\text{{Tr}}\mathopen{}\left({\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}^{{\top}}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\right)\geq\gamma\cdot\text{{Tr}}\mathopen{}\left(\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}^{{\top}}{\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}}\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{\star}\right)$ for some $0<\gamma<1$,

Let $\mathbf{E}{\triangleq}\mathbf{A}-\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$. By the linearity of the trace,

where the first inequality follows from (33) the second from (30), the third from (34), and the last from the fact that $R\geq 0$ and $0<\gamma\leq 1$. This concludes the proof. ∎

This follows from the fact that if $\mathbf{E}{\triangleq}\mathbf{A}-\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ is PSD, then

Further, by Lemma 10.37, the bound in (32) can be improved to

For any $n\times d$ input data matrix $\mathbf{S}$, with corresponding empirical covariance matrix $\mathbf{A}=\nicefrac{{1}}{{n}}\cdot\mathbf{S}^{{{\top}}}\mathbf{S}$, any desired number of components $k$, and accuracy parameters $\epsilon\in(0,1)$ and $r$, Algorithm 3 outputs ${\mathbf{X}_{{(r)}}\in\mathcal{X}_{k}}$ such that

where $\mathbf{X}_{\star}{\triangleq}\operatorname*{arg\,max}_{\mathbf{X}\in\mathcal{X}_{k}}\text{{Tr}}\mathopen{}\left(\mathbf{X}^{{{\top}}}\mathbf{A}\mathbf{X}\right)$, in time T_{\texttt{SKETCH}}(r)+T_{\texttt{SVD}}(r)+O\mathopen{}\bigl{(}\bigl{(}\tfrac{4}{\epsilon}\bigr{)}^{r\cdot{{k}}}\cdot{d}\cdot({s}\cdot{k})^{2}\bigr{)}.

The theorem follows from Lemma 9.29 and the approximation guarantees of Algorithm 1. ∎

4 Proof of Theorem 3

First, we restate and prove the following Lemma by .

for all $i,j$ with probability at least $1-1/d$.

Observe that each element of $\mathbf{A}$ is in $$, hence can be rewritten as an inner product of two unit-norm vectors:

Setting $r=O(\epsilon^{-2}\log d)$ and using the JL lemma and a union bound over all $O(d^{2})$ vector pairs ${\bf V}_{:,i}$, ${\bf V}_{:,j}$ we obtain the desired result. ∎

Next, we provide the proof of Theorem 3 for the simple case of $k=1$; the proof easily generalizes to the multi-component case $k>1$. According to Lemma 9.30, choosing $d=O\left((\delta/6)^{-2}\log{n}\right)=O\left(\delta^{-2}\log{n}\right)$ suffices for all entries of $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ constructed as described in the lemma to satisfiy

with probability at least $1-1/d$. In turn, for any $s$-sparse, unit-norm $\mathbf{x}$,

where the second inequality follows from the fact that $\mathbf{x}$ is $s$-sparse and unit norm.

We run Algorithm 1 (for $k=1$) with input argument the rank-$r$ matrix $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$, desired sparsity ${s}$ and accuracy parameter $\epsilon=\delta/6$. Algorithm 1 outputs a ${s}$-sparse unit-norm vector $\widehat{\mathbf{x}}$ which according to Theorem 1 satisfies

where $\mathbf{x}_{d}$ is the true $s$-sparse principal component of $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$. This, in turn, implies that $\widehat{\mathbf{x}}$ satisfies

where the second inequality follows from the fact that the entries of $\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$ lie in $[-1-\frac{\delta}{6},1+\frac{\delta}{6}]$ and $\widehat{\mathbf{x}}$ is ${s}$-sparse and unit-norm.

In the following, we bound the difference of the performance of $\widehat{\mathbf{x}}$ on the original matrix $\mathbf{A}$ from the optimal value. Let $\mathbf{x}_{\star}$ denote the ${s}$-sparse principal component of $\mathbf{A}$ and define

Utilizing (35) and the triangle inequality, one can verify that

where $(\alpha)$ follows from (37). Continuing from (38), combining (39) and (40) we obtain

which is the desired result. $\blacksquare$

Auxiliary Technical Lemmata

For any real $d\times n$ matrix $\mathbf{M}$, and any $r,k\leq\min\{d,n\}$,

where $\sigma_{i}(\mathbf{M})$ is the $i$th largest singular value of $\mathbf{M}$.

Note that $\sigma_{r+1}(\mathbf{M}),\ldots,\sigma_{r+k}(\mathbf{M})$ are the ${k}$ smallest among the $r+k$ largest singular values. Hence,

Combining the two inequalities, the desired result follows. ∎

For any real $d\times n$ matrix $\mathbf{M}$ and $k\leq\min\{d,n\}$, $\sigma_{{k}}(\mathbf{M})\leq k^{-1/2}\cdot\|\mathbf{M}\|_{{\textnormal{F}}}$.

It follows immediately from Lemma 10.31. ∎

Let $a_{1},\ldots,a_{n}$ and $b_{1},\ldots,b_{n}$ be $2n$ real numbers and let $p$ and $q$ be two numbers such that ${1/p}+{1/q}=1$ and $p>1$. We have

For any two real matrices $\mathbf{A}$ and $\mathbf{B}$ of appropriate dimensions,

Let $\mathbf{b}_{i}$ denote the $i$th column of $\mathbf{B}$. Then,

Similarly, using the previous inequality,

Combining the two upper bounds, the desired result follows. ∎

The inequality follows from Lemma 10.32 for ${p=q=2}$, treating $\mathbf{A}$ and $\mathbf{B}$ as vectors. ∎

For any real $m\times n$ matrix $\mathbf{A}$, and any $k\leq\min\{{m},\;{n}\}$,

The maximum is attained by $\mathbf{Y}$ coinciding with the ${k}$ leading right singular vectors of $\mathbf{A}$.

Let $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{{\top}}$ be the singular value decomposition of $\mathbf{A}$; $\mathbf{U}$ and $\mathbf{V}$ are $m\times m$ and $n\times n$ unitary matrices respectively, while $\Sigma$ is a diagonal matrix with ${\Sigma_{jj}=\sigma_{j}}$, the $j$th largest singular value of $\mathbf{A}$, $j=1,\ldots,d$, where $d{\triangleq}\min\{{m},{n}\}$. Due to the invariance of the Frobenius norm under unitary multiplication,

Let $z_{j}{\triangleq}\sum_{i=1}^{{k}}\left(\mathbf{v}_{j}^{{\top}}\mathbf{y}_{i}\right)^{2}$, $j=1,\ldots,d$. Note that each individual $z_{j}$ satisfies

where the last inequality follows from the fact that the columns of $\mathbf{Y}$ are orthonormal. Further,

Finally, it is straightforward to verify that if $\mathbf{y}_{i}=\mathbf{v}_{i}$, $i=1,\ldots,{{k}}$, then (42) holds with equality. ∎

For any real $d\times n$ matrix $\mathbf{A}$, and pair of $d\times{{k}}$ matrix $\mathbf{X}$ and $n\times{{k}}$ matrix $\mathbf{Y}$ such that $\mathbf{X}^{{\top}}\mathbf{X}=\mathbf{I}_{{k}}$ and $\mathbf{Y}^{{\top}}\mathbf{Y}=\mathbf{I}_{{k}}$ with ${k}\leq\min\{{d},\;{n}\}$, the following holds:

where the last inequality follows from the fact that $\|\mathbf{X}\|_{{\textnormal{F}}}^{2}=\text{{Tr}}\mathopen{}\left(\mathbf{X}^{{\top}}\mathbf{X}\right)=\text{{Tr}}\mathopen{}\left(\mathbf{I}_{{k}}\right)={k}$. Combining with a bound on $\|\mathbf{A}\mathbf{Y}\|_{{\textnormal{F}}}$ as in Lemma 10.35, completes the proof. ∎

For any real $d\times d$ PSD matrix $\mathbf{A}$, and $k\times d$ matrix $\mathbf{X}$ with $k\leq d$ orthonormal columns,

where $\lambda_{i}(\mathbf{A})$ is the $i$th largest eigenvalue of $\mathbf{A}$. Equality is achieved for $\mathbf{X}$ coinciding with the ${k}$ leading eigenvectors of $\mathbf{A}$.

Let $\mathbf{A}=\mathbf{V}\mathbf{V}^{{\top}}$ be a factorization of the PSD matrix $\mathbf{A}$. Then, $\text{{Tr}}\mathopen{}\left(\mathbf{X}^{{\top}}\mathbf{A}\mathbf{X}\right)=\text{{Tr}}\mathopen{}\left(\mathbf{X}^{{\top}}\mathbf{V}\mathbf{V}^{{\top}}\mathbf{X}\right)=\|\mathbf{V}^{{\top}}\mathbf{X}\|_{{\textnormal{F}}}^{2}$. The desired result follows by Lemma 10.35 and the fact that $\lambda_{i}(\mathbf{A})=\sigma_{i}^{2}(\mathbf{V})$, $i=1,\ldots,d$. ∎

Additional Experimental Results