Approximation Algorithms for Sparse Principal Component Analysis
Authors: Agniva Chowdhury
Petros Drineas
David Woodruff
Samson Zhou
Venue:  (2021)
Keywords: Sparse PCA
Principal component analysis
Randomized linear algebra
Singular value decomposition

Abstract: Principal component analysis (PCA) is a widely used dimension reduction technique in machine learning and multivariate statistics. To improve the interpretability of PCA, various approaches to obtain sparse principal direction loadings have been proposed, which are termed Sparse Principal Component Analysis (SPCA). In this paper, we present three provably accurate, polynomial time, approximation algorithms for the SPCA problem, without imposing any restrictive assumptions on the input covariance matrix. The first algorithm is based on randomized matrix multiplication; the second algorithm is based on a novel deterministic thresholding scheme; and the third algorithm is based on a semidefinite programming relaxation of SPCA. All algorithms come with provable guarantees and run in low-degree polynomial time. Our empirical evaluations confirm our theoretical findings.

Under review as a conference paper at ICLR 2021
APPROXIMATION ALGORITHMS FOR SPARSE PRINCI- PAL COMPONENT ANALYSIS
Anonymous authors Paper under double-blind review
ABSTRACT
Principal component analysis (PCA) is a widely used dimension reduction tech- nique in machine learning and multivariate statistics. To improve the interpretability of PCA, various approaches to obtain sparse principal direction loadings have been proposed, which are termed Sparse Principal Component Analysis (SPCA). In this paper, we present three provably accurate, polynomial time, approximation algorithms for the SPCA problem, without imposing any restrictive assumptions on the input covariance matrix. The first algorithm is based on randomized matrix multiplication; the second algorithm is based on novel deterministic thresholding scheme; and the third algorithm is based on semidefinite programming relaxation of SPCA. All algorithms come with provable guarantees and run in low-degree polynomial time. Our empirical evaluations confirm our theoretical findings.
INTRODUCTION
Principal Component Analysis (PCA) and the related Singular Value Decomposition (SVD) are fundamental data analysis and dimension reduction tools in a wide range of areas including machine learning, multivariate statistics and many others. They return set of orthogonal vectors of decreasing importance that are often interpreted as fundamental latent factors that underlie the observed data. Even though the vectors returned by PCA and SVD) have strong optimality properties, they notoriously difficult to interpret in terms of the underlying processes generating the data (Mahoney & Drineas, 2009), since they are linear combinations of all available data points or all available features. The concept of Sparse Principal Components Analysis (SPCA) was introduced in the seminal work of (d' Aspremont al., 2007), where sparsity constraints were enforced on the singular vectors in order to improve interpretability. A prominent example where sparsity improves interpretability document analysis, where sparse principal components can be mapped to specific topics by inspecting the (few) keywords in their support (d' Aspremont et al., 2007; Mahoney & Drineas, 2009; Papailiopoulos et al., 2013).
Formally, given positive semidefinite (PSD) matrix A E R SPCA can be defined as follows:!
Z max Ix x Ax,
subject to ||x||o < k.
(I)
In the above formulation, A a covariance matrix representing, for example, all pairwise feature object similarities for underlying data matrix. Therefore, SPCA can be applied for either the object or feature space of the data matrix, while the parameter i controls the sparsity of the resulting vector and is part of the input. Let x" denote vector that achieves the optimal value 2" in the above formulation. Then intuitively, the optimization problem of eqn. (1) seeks a sparse, unit norm vector x" that maximizes the data variance.
is well-known that solving the above optimization problem is NP-hard (Moghaddam et al., 2006a) and that its hardness due to the sparsity constraint. Indeed, if the sparsity constraint was removed, then the resulting optimization problem can be easily solved by computing the top left or right singular vector of A and its maximal value 2 equal the top singular value of A.
Notation. We use bold letters denote matrices and vectors. For matrix A E guxx we denote its (-th entry by A its i-th row by Ais and its j-th column by A.j; its 2-norm by
Recaal that the p-th power of the Ep norm of vector E R is defined as ||x||p = Σn, |x|P for 0 < p co. For p- 0, [|z||o is semi-norm denoting the number of non-zero entries of

---

Under review as a conference paper at ICLR 2021
|| maxxeR" ||x||==1 |||x| and its (squared) Frobenius norm by = A We use the notation A č 0 to denote that the matrix A is symmetric positive semidefiitte (PSI) and Tr(A) = E. a to denote its trace, which is also equal to the sum of its singular values. Given a PSD matrix E R=**, its Singular Value Decomposition is given by A = UEU*, where the matrix of left/right singular vectors and Σ is the diagonal matrix of singular values.
1.1 OUR CONTRIBUTIONS
We present three algorithms for SPCA and associated quality-of-apprxximato results (Theo- rems 2.2, 3.1, and 4.1). All three algorithms are simple, intuitive, and run in O (n..a) or less time. They return vector that is provably sparse and, when applied to the input covariance matrix A, provably captures a fraction of the optimal solution Z We note that in all three algorithms, the output vector has sparsity that depends on l (the target sparsity of the original SPCA problem of eqn. (1)) and (an accuracy parameter between zero and one).
The first algorithm is based randomized, approximate matrix multiplication: it randomly (but non-uniformly) selects a subset of O (k/e²) columns of A 1/22 (the square root of the PSD matrix A) and computes its top right singular vector. The output of this algorithm is precisely this singular vector, padded with zeros to become vector in R". It turns out that this simple algorithm, which, surprisingly has not been analyzed in prior work, returns an O (k/e²) sparse vector y E R that satisfies (with constant probability that can be amplified desired, see Section 2 for details):
1 y Ay 2 -'
Tr(A) k
Notice that the above bound depends on both 3d and square root and therefore is not a relative error bound. The second term scales as function of the trace of A divided by k which depends on the properties of the matrix A and the target sparsity.
The second algorithm is deterministic thresholding scheme. It computes small number of the top singular vectors of the matrix A and then applies deterministic thresholding scheme on those singular vectors to (eventually) construct sparse vector z E R" that satisfies z Az 2 (+/2))2--(3))) (3/2)e Tr(A), ). Our analysis provides unconditional guarantees for the accuracy of the solution of this simple thresholding scheme. To the best of our knowledge, (see Section 1.2 for details). The eиror bound of the second algorithm is weaker than the one provided by the first algorithm, but the second algorithm is deterministic and does not need to compute the square root (i.e., all singular vectors and singular values) of the matrix A.
such analyses have appeared in prior work
Our third algorithm provides novel bounds for the following standard convex relaxation of the problem of eqn. (1).
Tr(AZ) s.t. Tr(Z) S and Σ Zr.al S k.
max
ZERAX ZZ
is well-known that the optimal solution to eqn. (2) least the optimal solution to eqn. (1). We present novel, two-step rounding scheme that converts the optimal solution matrix Z E R to vector z E IR" that has expected sparsity² Ô (1²/²²) and satisfies z Az 2 72 -) 2 Here, 7z is a constant that precisely depends on the top singular value of Z, the condition number Z, and the extent to which the SDP relaxation of eqn. (2) is able to capture the original problem (see Theorem 4. I and the following discussion for details). To the best of our knowledge, this is the first analysis of a rounding scheme for the convex relaxation of eqn. (2) that does not assume specific model for the covariance matrix A.
Applications to Sparse Kernel PCA. Our algorithms have immediate applications sparse kernll PCA (SKPCA), where the input matrix A E TR x is instead implicitly given as a kemne matrix whose entry (i,j) is the value k(i, j) (o(X,.), i..,, o (X),)) for some kernel function & that implicitly maps an observation vector into some high-dimensional feature space. Although A is not explicit,
FFor simplicity of presentation and following the lines of (Fountoulakis et al., 2017). assume that the columns of the matrix have unit norm; this assumption was not necessary for previous two algorithms and can be removed as in (Fountoulakis et al., 2017). We are also hiding a poly-logarithmic factor for simplicity, hence the Ó(:) notation. See Theorem 4.1 for detailed statement.
2

---

Under review as a conference paper at ICLR 2021
we can query all O (n²) entries of A using O (n²) time assuming an oracle that computes the kernel function k. We can then subsequently apply our SPCA algorithms and achieve polynomial time runtime with the same approximation guarantees.
1.2 PRIOR WORK
SPCA was formally introduced by (d' Aspremont et al., 2007); however, previously studied PCA approaches based on rotating (Jolliffe, 1995) or thresholding (Cadima & Jolliffe, 1995) the top singular vector of the input matrix seemed to work well, at least in practice, given sparsity constraints. Following (d' Aspremont et al., 2007), there has been an abundance of interest in SPCA. (Jolliffe et al., 2003) considered LASSO (SCoTLASS) on an l relaxation of the problem, while (Zou & Hastie, 2005) considered a non-convex regression-type approximation, penalized similar to LASSO. Additional heuristics based on LASSO (Ando et al., 2009) and non-convex & regularizations (Zou & Hastie, 2005: Zou et al., 2006; Sriperumbudur al., 2007: Shen & Huang, 2008) have also been explored. Random sampling approaches based on non-convex l relaxations (Fountoulakis et al., 2017) have also been studied; we highlight that unlike our approach. (Fountoulakis et al., 2017) solved a non-convex relaxation of the SPCA problem and thus perhaps relied on locally optimal solutions. Additionally, (Moghaddam et al., 2006b) considered a branch-and-boum heuristic motivated by greedy spectral ideas. (Journée et al., 2010; Papailiopoulos et al., 2013; Kuleshov, 2013; Yuan & Zhang, 2013) further explored other spectral approaches based on iterative methods similar to the power method. (Yuan & Zhang, 2013) specifically designed a sparse PCA algorithm with early stopping for the power method, based on the target sparsity.
Another line of work focused on using semidefinite programming (SDP) relaxations (d' Aspremont et al., 2007: d' d'Aspeemont et al., 2008; Amini & Wainwright, 2009). Notably, (Amini & Wain- wright, 2009) achieved provable theoretical guarantees regarding the SDP and thresholding approach of (d' Aspremont et al., 2007) in a specific, high-dimensional spiked covariance model, in which a base matrix is perturbed by adding sparse maximal eigenvector. In other words, the input matrix is the identity matrix plus a "spike", i.c., a sparse rank-one matrix.
Despite the variety of heuristic-based sparse PCA approaches, very few theoretical guarantees have results.
been provided for SPCA; this is partially explained by line of The sparse PCA problem is well-knows to be NP-Hard (Moghaddam et al., 2006a). (Magdon-Ismaii, 2017) shows that if the input matrix is not PSD., then even the sign of the optimal value cannot be determined in polynomial time unless P = NP, ruling out any multiplicative approximation algorithm. In the case where the input matrix is PSD, (Chan et al., 2016) shows that it is NP-hard to approximate the optimal value up to multiplicative (1 + c) error, ruling out any polynomial- time approximation scheme (PTAS). Moreover, they show Small-Set Expansion hardness for any polynomial-time constant factor approximation algorithm and also that the standard SDP relaxation might have an exponential gap.
We conclude by summarizing prior work that offers provable guarantees (beyond the work of (Amini & Wainwright, 2009)), typically given some assumptions about the input matrix. (d' Aspremont et al., 2014) showed that the SDP relaxation can be used to find provable bounds when the covariance input matrix is formed by number of data points sampled from Gaussian models with single sparse singular vector. (Papailiopoulos et al., 2013) presented combinatorial algorithm that analyzed a specific set of vectors in a low-dimensional eigenspace of the input matrix and presented relative error guarantees for the optimal objective, given the assumption that the input covariance matrix has decaying spectrum. (Asteris et al., 2011) gave a polynomial-time PCA exactly for input matrices of constant rank. (Chan et al., 2016) showed that sparse PCA can be approximated in polynomial time within a factor of n 11/3 and also highlighted an additive PTAS of (Asteris et al., 2015) based the idea of finding multiple disjoint components and solving bipartite maximum weight matching problems. This PTAS needs time poly(1)))) whereas all of our algorithms have running times that are low-degree polynomial in n.
algorithm that solves sparse
2 SPARSE PCA VIA RANDOMIZED MATRIX MULTIPLICATION
Our first algorithm for SPCA leverages primitives and ideas from randomized matrix multiplica- tion (Drineas & Kannan, 2001; Drineas et al., 2006; Drineas & Mahoney, 2016; 2018; Woodruff, 2014). Let P E Rm and Q E RR and recall that their product equals PQ EE P., Qi.
3

---

Under review as a conference paper at ICLR 2021
Recall that P.; denotes the i-th column of P and Qix the i-th row of Q. A well-known approach to approximate the product PQ is to sample a subset of columns of P (we will do this without replacement) and the corresponding rows of Q (Drineas & Mahoney, 2018). Formally, let the random variables Z Ad Bernouli(p;). = 1....M n, denote whether the i-th column of P and the i-th row of Q are sampled. Define the diagonal sampling-and-rescallm and-rescaling matrix S R"x" by 2n/vpn). The sampling probabilities {pi) do not have sum up to one. The number of sampled column/row pairs, denoted by satisfies E Σi=1 Pi- (See Algorithm for details.) The next lemma (see Appendix A.1 for its proof) presents accuracy bounds
& diag (Zt//PT,...
Algorithm Construct sampling-andee and-rescaling matrix
Input: Probabilities jı: i = 1... n and integer < n.. Output: Diagonal sampling-and-resene 1: for i = to n do
matrix S E Raxx
S 1/VF, with probability ps min sp., 1}
2:
0.
otherwise
when Algorithm 1 is used approximate matrix multiplication.
Lemma 2.1 Given matrices P E Rmxn and Q E RRXP,, let S E Raxa be constructed using Algorithm 1 with P4 P., 8/IP| 2 for I...n. Then,
E PS
PQI S
P QQI
Our SPCA algorithm uses the above primitive to approximate the product of the (square root) of the input matrix A and its top right singular vector v. Thus, the proposed SPCA algorithm sparsifies the top right singular vector of of A without losing too much of the variance that is captured by v. Interestingly, this conceptually simple algorithm has not been formally analyzed in prior work. Algorithm details approach.
Algorithm 2 SPCA via randomized matrix multiplication
Input: A E R=x", sparsity parameter k. accuracy parameter c E (0, 1).
Output: y E R satisfying E [||y||a] and E [[yy|lo [lo] < k/e2.
1. x + 1/2
2. Use Algorithm I to construct R with p. = [X_.: 1Xxi1///111 and = 1k/2:
3. Let E RR be the top right singular vector of XS;
4. y + Sv:
Theorem 2.2 Let k be the sparsity parameter and e E (0, 1] be the accuracy parameter. Let Raxa be the sampling matrix of Lemma 2.1 with s = 2,/2. Then, Algorithm returns a vector y with expected sparsity at most (i.e., TE. [Ily]la]] S 8) and expected two norm at most one (i.e.,
E || < 1) such that with probability at least 1/4, we have
y Ay 2 1/22" - VTTAA/E.
(4)
See Appendix A. for a proof of the above theorem. We note that the success probability of Algo- rithm 2 can be trivially amplified by repeating the algorithm t times and keeping the vector y that maximizes y Ay. Then, the failure probability of the overall approach diminishes exponentially fast as a function of t to at most (3/4). Finally, the running time of Algorithm is dominated by the computation of a square root of the matrix A the first step, which takes O (n²) time via the computation of the SVD of A.
3 SPCA VIA THRESHOLDING
Our second algorithm is based on a thresholding scheme using the top right singular vectors of the PSD matrix A. Given A and accuracy parameter e, our approach first computes Σ, E cx
4

---

present provable bounds for Under review as a conference paper at ICLR 2021
(the diagonal matrix of the top E singular values of A) and U2 E R (the matrix of the top l right singular vectors of A), for & 1/c. Then, it deterministically selects a subset of O (k/e²) columns of E using a simple thresholding scheme based on the norms of the columns of _ U (Recall that & is the sparsity parameter of the SPCA problem.) In the last step, it returns the top right singular vector of the matrix consisting of the chosen columns of 2 Notice that this right singular vector is an O (k/e")-dimensionall -dimensional vector, which is finally expanded vector in R by appropriate padding with zeros. This sparse vector is our approximate solution to the SPCA problem of eqn. (1).
This simple algorithm is somewhat reminiscent of prior thresholding approaches for SPCA. However, to the best of our knowledge. no provable a priori bounds were known for such algorithms without strong assumptions the input matrix. This might be due to the fact that prior approaches focused on thresholding only the top right singular vector of A, whereas our approach thresholds the top f = 1/c right singular vectors of A. This slight relaxation allows the proposed algorithm.
In more detail, let the SVD of A b UEUT Let Eg E Rex be the diagonal matrix of the top l singular values and let U, E Ruxe be the matrix of the top l right (or left) singular vectors. Let R = (i1,...,,41n } be the set of indices of rows of U. that have squared norms at least e/k and let R be its complement. Here |t denotes the cardinality of the set and R = (I,...,,)) n). Let R E Ra be sampling matrix that selects² the columns of U2 whose indices are in the set R. Given this notation, we are now ready to state Algorithm 3.
Algorithm 3 SPCA via thresholding
Input: A R"xn, sparsity k error parameter > 0. Output: y E I2" such that ||y||2 = 1 and lyyloo = k/e².
1: ( 1/c;
2: Compute U. E TR (top & left singular vectors of A) and Eg E Rixt (square roots the top &
singular values of A);
3: Let R = {i1,.... ((8)) be the set of rows of U, with squared norms at least e/k and let
R. E |nx|R| be the associated sampling matrix (see text for details);
4: y E R + argmax |x E.U, RxIl;;
5: return z Ry E R":
Notice that Ry satisfies |Ryl| = ||yl|2 = 1 (since R has orthogonal columns) and |Ryll = [|] Since R. is the set of rows of U2 with squared norms at least e/k and || = = = I/c. it follows that R S k/e2. Thus, the vector returned by Algorithm has k/e² sparsity and unit norm.
Theorem 3.1 Let k be the sparsity parameter and e E (0, 1] be the accuracy parameter. Then, the vector E R (the output of Algorithum 3) has sparsity k/e,, unit norm, and satisfies z Az Z (1/2)2* - (9/2))T((AA)
We defer the proof of Theorem 3.1 to Appendix A.2. The running time of Algorithm 3 is dominated by the computation of the top singular vectors and singular values the matrix A. In practice, any iterative method, such as subspace iteration using random initial subspace or the Krylov subspace of the matrix, can be used towards this end. However, our current analysis does not account for the inevitable approximation error incurred by such methods, which run in O (nnz( A JE) time. One could always the SVD of the full matrix A (O (n²) time) compute the top singular vectors and singular values of A. Finally, we highlight that, as an intermediate step in the proof of Theorem 3.1, we need to prove the following lemma (see Appendix A.2 for its proof):
Lemma 3.2 Let A Raxa be PSD matrix and Σ E Raxn (respectively, De (R(xE) be the diagonal matrix of all (respectively, top E) singular values and let Raxn (respectively. U, E R be the matrix of all (respectively, top () singular vectors. Then, for all unit vectors x E R",
- (T(A).
Fach column of R. has single non-zero entry (set to onc), corresponding one of the |R] selected columns. Formally, Ra,. for 1...., R all other entries of R. to zero
5

---

Under review as a conference paper at ICLR 2021
The above simple lemma is very much at the heart of our proof of Theorem 3.1 and, unlike prior work, allows us to provide provably accurate bounds for the thresholding Algorithm
Using an approximate SVD solution. The guarantees of Theorem 3.1 1 in Algorithm 3 uses an exact SVD computation, which could take time O (n²)). We can further improve the running time by using an approximate SVD algorithm such as the randomized block Krylov method of Musco & Musco (2015), which runs in nearly input spansity runtime. Our analysis uses the relationships
1/2
Tr(A) and o (Ex) < Tr(A). The randomized block Krylov method of Musco & Musco e
(2015) recovers these guarantees up to a multiplicative (1 +)), using O log
nnz((A) runtime.
Thus by rescaling e, we recover the same guarantees of Theorem 3.1 by using approximate SVD with nearly input sparsity runtime.
4 SPCA VIA A SEMIDEFINIE PROGRAMMING RELAXATION
Our third algorithm is based on the SDP relaxation of eqn. (2). Recall that solving eqn. (2) returns a PSD matrix Z E IR"x* that, by the definition of the semidefinite programming relaxation, satisfies Tr(AZ') ) 2 Z", where 2" the true optimal solution of SPCA in eqn. (2).
We now need to convert the matrix Z" E R into sparse vector that will be the output of our approximate SPCA algorithm and will satisfy certain accuracy guarantees. Towards that end, we employ novel two-step rounding procedure. First, critical observation is that generating random Gaussian vector g E IR" and computing the vector Z" gggg E R results in an unbiased estimator for the trace of (Z*) AZ' in the following sense: E [g (Z*) AZ'g g Tr((Z") TZZ))
Using von Neumann's 's trace inequality, we can prove that
E g (Z')
AZ'gl Tr((Z") AZ') 2 22. Tr(AZ") 2 72 2"
(5)
Here 72- is a constant that precisely depends on the top singular value of Z", the condition number of Zz, and the extent which the SDP relaxation of eqn. (2) able to capture the original problem (see Theorem 4.1 for the exact expression of 72-). The above inequality implies that, at least in expectation, we could use the vector Z g as "rounding" of the output of the semidefinite programming relaxation. However, there is absolutely no guarantee that the vector Z" g is sparse. Thus, in order to sparsify Z" g. we employ a separate sparsification procedure, where each entry of Z' g is kept (and rescaled) with probability proportional to its magnitude. This procedure is similar to the one proposed in (Fountoulakis et al., 2017) and guarantees that larger entries Z gg are more likely to be kept, while smaller entries of Z" g are more likely to be set to zero, without too much loss in accuracy. We also note that to ensure a sufficiently high probability of success for the overall approach, generate multiple Gaussian vectors and keep the one that maximizes the quantity (Z*)* AZ'g. g. See Algorithm 4 and Algorithm 5 for a detailed presentation of our approach.
Algorithm SPARSIFY
Input: y E _ and sparsity parameter 8.
Output: z E 12" with E [|| 5 s.
for 1. 1... do
with probability p; min (1, otherwise.
Yi:
2:
2.
0
The running time of the algorithm is dominated by the time needed to solve the semidefinite pro- gramming relaxation of eqn. (2), which, our setting, O (n.55) (Alizadeh, 1995). We do note that SDP solvers such as the one in (Alizadeh, 1995) return additive error approximation to the optimal solution. However, the running time dependence of SDP solvers on the additive error y is logarithimic in 1/1y and thus highly accurate approximations can be derived without significant increase in the number of iterations of the solver. Thus, for the sake of clarity, we initially omit this additive error from the analysis and address the approximate solution the end.
Our main ofality-of-ppproxmmatio result for Algorithm 5 is Theorem 4.1. For simplicity of presen- tation, and following the lines of (Fountoulakis et al., 2017), we assume that all rows and columns

---

Under review as a conference paper at ICLR 2021
Algorithm 5 Rounding-based SPCA
Input: PSD matrix A R error tolerance
0, and sparsity parameter k.
Output: x E I= with E
[|x [xl]=== =s.
Let Z be the optimal solution the relaxed SPCA problem of eqn. (2);
O k2 (141)
2: M <- 80/e² and
PSee Theorem A.13 for the exact value of
3: Generate M random Gaussian vectors 9 91.... 9M in I";
4: y Z'gg, where j < =r.. gi (Z") AZ" g::
5: Z - (PARSIFY(Y,s))
of A have been normalized to have unit norms. This assumption can be relaxed as in (Fountoulakis al., 2017). In the statement of the theorem, will use the notation Z₁ to denote the best rank-one approximation to the matrix Z.
Theorem 4.1 Given a PSD matrix A E Rn a sparsity parameter k. and an error tolerance e > let Z be an optimal solution the relaxed SPCA problem of eqn. (2). Assume that
0,
Tr(AZ) & a Tr(AZ))
(6)
for some constant a 2 1. Then, Algorithm 5 outputs a vector z E R" that, with probability at least 5/8, satisfies E |||[||| O (Ie log (2/e/22)))) |z|12 O (vlog and
Az 2 72 E) 2 - z
Here 7i = (1--11- x(Z) (1- (Z) with 01 (Z) and (Z) being the top singular value and condition number of ZZ respectively.
Similar Theorem 2.2, the probability of success be boosted 1 by repeating the algorithm O (t) times in parallel. Moreover by using Markov's inequality, we can also guarantee a vector with sparsity O (k* log"/((\/=)/22)) with probability 1 ô, rather than just in expectation.
now discuss the condition of eqn. (6) and the constant 72. Our assumption simply says that much of the trace of the matrix AZ should be captured by the trace of AZ1, as quantified by the constant were a rank-one matrix, then the assumption would hold with a = 1. As the trace of AZZ fails to approximate the trace of AZ (which intuitively implies that the SDP relaxation of eqn. (2) did not sufficiently capture the original problem) the constant & increases and the quality of the approximation decreases. More precisely, first notice that the constant 72 is upper bounded by one, because 01 (Z) < by the SDP relaxation. Second, the quality of the approximation increases as 7z approaches one. This happens if either the condition number of is close one if the constant is close to one; at the same time, 01 (Z) also needs to be close to one. Clearly, these conditions are satisfied if Z is well approximated by Z1 In our experiments, we indeed observed that is close to one and that the top singular value of Z is close to one, which imply that 7z is also close to one (Appendix, Table 6). The proof of Theorem
a. For example,
delegated to Appendix A.3 (as Theorem A. A3), we outline here summary of statements that lead to the final bound.
Lemma 4,2 Let y and z. be defined as in Algorithm 5. If ||y||, S a and ||y||2 S B, then
ly Ay z. Az < y A(y z)1 ( z) A(y z)).
Moreover, with probability at least 7/8, we have ly A(y 2)1 - 408/ s and ((y T) 2)
64o
96
Lemma 4.3 Let M- = 80/2, a = k(1 2 2og M), and 3 = 20/oo M. If the sparsity parameter s is set to 450g /²/2, then with probability at least 3/4, we have y Ay < z Az +
Letting y Z g, now conclude the proof by combining eqn. (5) with the above lemma to bound (at least in expectation) the accuracy of Algorithm 5. To get high probability bound, we leverage a result by (Avron & Toledo, 2011) 1) on estimating the trace of PSD matrices. This approach allows us
7

---

Under review as a conference paper at ICLR 2021
to properly analyze step 4 of Algorithm 5, which uses multiple random Gaussian vectors to achieve measure concentration (see Appendix A.3 for details).
Finally, we can bound the & norm of the vector of Algorithm 5 by proving that, with probability at least 4/4, |zll| = O (/e) Notice that this slightly relaxes the requirement that has unit norm; however, even for accuracy close machine precision, Vlog 1/e is small constant.
Using approximate SDP solution. The guarantees of Theorem 4.1 in Algorithm 5 uses optimal solution Z to the SDP relaxation in eqn. (2). In practice, we will only obtain approximate solution to eqn. (2) using any standard SDP solver, e.g. (Alizadeh, 1995), such that Tr(AZ) > Tr Tr(AZ))
after O (log] )) iterations. Since our amalysis only uses the relationship (x*) ) Ax" I Tr(AZ) then the additive guarantee Thus, we recover the same guarantees of Theorem 4.1 by using an approximate solution to the SDP relaxation in eqn. (2).
be absorbed into the other factors in the guarantees of Theorem 4.1.
5 EXPERIMENTS
05
155
0.4
34
3.0
0.2
01.
32-
0.7
Seciity 155
Apaincy 1993
Spancdy
(a) CHR I, 37, 493
(b) CHR 2, n 40, 844
(c) Gene expression data
Fig. 1: Experimental results on real data: f(y) vs. sparsity.
lpplian
lagg
most
tedtig
ovang
Gaggop
riggls
bowmrss
Toudia shocls
Jeer
kross
dikkst
PVE
2"
spca-d (Alppesten 23
c2)
622
-0.915
apca-ado (Aleseihen.)
-0.124
(1)
-0.28
-0.00
-0.315
-0.50
-0.310
44:247
3,55M
dec (hin clal. 2019)
-0 423
01.33)
-0.08
-0.053
-0.315
-0.50
-6.350
0
96.240%
3.504
espea (Blckk Vildwonnt 20100
-0311
370
-0.350
0
74%
3.946
(Psuulllopetto
(Poe 2006)
-0.157
-0.190 170
-0.990
-0.311
01
-0.450
3.641
depoa (et Asporomn al. 1007)
-0.00
-0.507
0
-B007
-0.357
-0.931
-0.087
-0.130
0
43.3995
3.521
Table 1: Loadings, of variance explained (PVE), and the objective function value for the first pe of the Pit Props data.
We compare the outputs our algorithms with that of the state-of-the-art sparse PCA solvers such as the coordinate-we wise optimization algorithm of Beck & Vaisbourd (2016) (cwpea) and the block decomposition algorithm of Yuan et al. (2019) (dec), along with other standard methods such Papailiopoulos al. (2013) (spca-lorrank) d' Asprrmmmnt et al. (2007) (dspca), and Zou et al. (2006) (spca). For the implementation of dec, we use coordinate descent method with workset: (6, 6) and for ewpca, we use greedy coordinate wise (GCW) method.
First, in order to explore the sparsity patterns of the outputs of our algorithms and how close they compared to standard methods, we first apply our methods on the pil props data which was introduced in (Jeffers, 1967) and is a benchmark example used to test sparse PCA. It is 13 x 13 correlation matrix, originally calculated from 180 observations with 13 explanatory variables. While the existing algorithms aimed to extract top 6 principal components, we restrict only to the top principal component. In partieular, we are interested to apply our algorithms on the Pit Props matrix in view to extract the top pc having a sparsity pattern similar to that of Beck & Vaisbourd (2016). Yuan et al. (2019), (Zou et al., 2006), Papailiopoulos et al. (2013), and (d' Aspremont et al., 2007). is known that the decomposition method of Yuan et al. (2019) can find the global optimal solution. We take A = 7 and Table shows that while spea-d (Algorithm 3) and spca-r (Algorithm 2) perform very similar to spca that of dspca with % of variace explained (PVE) uniformly better than spea, our based method spca-sdp (Algorithm 5) exactly recovers the optimal solution and the output matches with both dec and cwpca and very close spca- aoccalloranti We also
8

---

Under review as a conference paper at ICLR 2021
apply our algorithms another benchmark artificial example from (Zou et al., 2006), please Appendix B.3 for detailed discussion.
Tightness of our bounds. Now, we also verify the tight- ness of the theoretical lower bounds of our results with the guarantee of (Papailiopoulos et al., 2013) on the pir props data. We take e 0.1 and found that the lower bound of spoa-sca (Algorithm 5) (dashed red line the left) is indeed very close to that of Papailiopoulos et al. (2013) with d=3 (dashed grey line on the left) Nevertheless, the accuracy parameter of Papailiopoulos et al. (2013) typi- cally relies the spectrum of A i.e., for highly accurate output. can be much smaller depending on the structure of A, which case the difference between the lower bounds of our Algorithm 5 and Papailiopoulos et al. (2013) becomes
30
12
Sperdty: 19
smaller.
Next, we further demonstrate the empirical performance of our algorithms on larger real-world datasets as well as on synthetic dataset, similar to (Fountoulakis al., 2017) (see Appendix B). We use human genetic data from the HGDP (Consortium, 2007) and HAPMAP (Li et al., 2008) (22 matrices, one for each chromosome). In addition, we also use a lung cancer gene expression dataset (107 x 22, 215) from (Landi et al., 2008) and sparse document-tent matrix (2. 858 x 12, 427) created using the Text-to-Matrix Generator (TMG) (Zeimpekis Gallopoulos, 2006) (see Appendix B).
Comparisons and metries. We compare our Algorithm 2 (spca-r), Algorithm 3 (spoa-d)) and Algorithm (spoa sdp) with the solutions returned by spca (Zou et al., 2006), as well as the simple MaxComp heuristic (Cadima & Jolliffe, 1995). We define the quantity f(y) = Av/ A to measure the quality of an approximate solution y E IR" to the SPCA problem. Notice that As f(y) gets closer to one, the vector y captures more of
0 f(y) < for all y with ||y|| the variance of the matrix A that corresponds to its top singular value and corresponding singular vector. Our goal is to identify sparse vector y with f(y) s Since the outputs of our Algorithm and Algorithm may have norms that slightly exceed one and in order to have a fair comparison between different methods, normalize our outputs in the same way as in (Fountoulakis al., 2017) (Appendix B).
Results. In our experiments, for spca-d, spca-r, and spca-sdp, we the sparsity s to be equal to k. that all algorithms return the same number of non-zero elements. In Figures la-le Ic evaluate the performance the different SPCA algorithms by plotting f(y) against |y i.e., the sparsity of the output vector, on data from chromosome 1, chromosome 2, and the gene expression data. Note that performance of our SDP-based method (spca-sdp) is indeed comparable with the state-of-te oo-the--rt dec, owpca, and spca-lowrank, while both spca-d and spca-r are better than or al least comparable to both spca and maxcomp. However, in practice, the running time of the SDP relaxation is substantially higher than our other methods, which highlights the interesting trade-offs between the accuracy and computation discussed in (Amini & Wainwright, 2009). See Apeendix B for more experimental results.
CONCLUSION AND OPEN PROBLEMS
We present three provably accurate, polynomial time, approximation algorithms for SPCA, without imposing restrictive assumptions on the input covariance matrix. Future directions include: (i) extend the proposed algorithms handle more than one sparse singular vector by deflation or other strategies and (ii) explore matching lower bounds and/or improve the guarantees of Theorems 2.2, 3.1, and 4.1.

---

Under review as a conference paper at ICLR 2021
REFERENCES
Farid Alizadeh. Interior Point Methods in Semidefinite Programming with Applications to Combina-
torial Optimization. SIAM Journal on Optimization, 5(1):13-51, 1995. 6,
Arash A. Amini and Martin J. Wainwright. High-dimensional Analysis of Semidefinite Relaxations for Sparse Principal Components. Annals of Statistics, 37:2877-2921, 2009. 3,9
Ei Ando, Toshio Nakata, and Masafumi Yamashita. Approximating the Longest Path Length of Stochastic DAG by Normal Distribution in Linear Time. Journal of Discrete Algorithms, 7(4): 420-438, 2009. 3
Megasthenis Asteris, Dimitris Papailiopoulos, and George N Karystinos, Sparse Principal Component of Rank-deficient Matrix. In 2011 IEEE International Symposium on Information Theory Proceedings, pp. 673-677, 2011.
Megasthenis Asteris, Dimitris Papailiopoulos. Anastasios Kyrillidis, and Alexandros G Dimakis. Sparse PCA via Bipartite Matchings. In Advances in Neural Information Processing Systems, pp. 766-774, 2015. 3
Haim Avron and Sivan Toledo. Randomized Algorithms for Estimating the Trace of an Implicit Symmetric Positive Semi-definite Matrix. Journal of the ACM, 58(2):8:1-8:34, 2011. 7,
Amir Beck and Yakov Vaisbourd. The Sparse Principal Component Analysis Problem: Optimality Conditions and Algorithms. Journal of Optimization Theory and Applications, 19-0)1199-14,, 2016. 8, 23
Jorge Cadima and Ian T. Jolliffe. Loading and Correlations the Interpretation of Principal Components. Journal of Applied Statistics, 22(2):203-22) 1995. 3,9
Siu On Chan, Dimitris Papailliopoulos, and Aviad Rubinstein. On the Approximability of Sparse
PCA. In Proceedings of the 29th Conference on Learning Theory, pp. 623-646, 2016.
International HapMap Consortium. A Second Generation Human Haplotype Map of over 3.1 Million SNPs. Nature, 449(71 (15)5))) 2007. 9, 22
Alexandre d' Aspremont, Laurent El Ghaoui, Michael I. Jordan, and Gert R. G. Lanckriet. A Direct Formulation for Sparse PCA using Semidefinite Programming. SIAM Review, 49(3):434 448, 2007. 1, 3, 8, 23, 24
Petros Drineas and Ravi Kannan. Fast Monte-Carlo Algorithms for Approximate Matrix Multipli- cation. In Proceedings of the 42nd IEEE Symposium on Foundations of Computer Science, pp. 452-461, 2001. 3
Petros Drineas and Michael W. Mahoney. RandNLA: Randomized Numerical Linear Algebra. Communications of the ACM, 59(6);00990 2016. 3
Petros Drineas and Michael W. Mahoney, Lectures on Randomized Numerical Linear Algebra, volume 25 of The Mathematics of Data, IAS/Park City Mathematics Series. American Mathematical Society, 2018. 3,4 4
Petros Drineas, Ravi Kannan, and Michael W. Mahoney, Fast Monte Carlo Algorithms for Matrices II: Computing Low-Rank Approximation to Matrix. SIAM Journal Computing, 36(1): 158-183, 2006.
Alexandre d' Asprrmonn,, Francis Bach, and Laurent El Ghaoui. Optimal Solutions for Sparse Principal Component Analysis. Journal of Machine Learning Research, 9(ul): 2008.
Alexandre d' Aspremont, Francis Bach, and Laurent El Ghaoui. Approximation Bounds for Sparse Principal Component Analysis. Mathematical Programming, 148(1-2):9--110, 2014. 3
Kimon Fountoulakis, Abhisek Kundu, Eugenia-Maria Kontopoulou, and Petros Drineas. A Random- ized Rounding Algorithm for Sparse PCA. ACM Transactions on Knowledge Discovery from Data 11(3):38:1338-22 2017. 2, 6, 7, 9, 16, 17, 23
10

---

Under review as a conference paper at ICLR 2021
John NR Jeffers. Two case studies in the application of principal component analysis. Journal of the
Royal Statistical Society: Series C (Applied Statistics), 16(3):225- 236, 1967.
lan T. Jolliffe. Rotation of principal components: Choice of Normalization Constraints. Journal of Applied Statistics, 22(1):29-35, 1995. 3
Ian T. Jolliffe, Nickolay Trendafilov, and Mudassir Uddin. A Modified Principal Component Technique Based on the LASSO. Journal of Computational and Graphical Statistics, 12(3): 531-547, 2003.
Michel Journée, Yurii Nesterov, Peter Richtárik, and Rodolphe Sepulchre, Generalized Power Method for Sparse Principal Component Analysis. Journal Machine Learning Research, 11(2), 2010.
Volodymyr Kuleshov.
Fast Algorithms for Sparse Principal Component Analysis Based Rayleigh Quotient Iteration. In Proceedings of the 30th International Conference on Machine Learning, pp.
1418-1425, 2013.
Maria Teresa Landi, Tatiana Dracheva, Melissa Rotunno, Jonine D. Figueroa, Huaitian Liu, Abhijit Dasgupta, Felecia E. Mann, Junya Fukuoka, Megan Hames, Andrew W. Bergen, et al. Gene Expression Signature Cigarette Smoking and Its Role Lung Adenocarcinoma Development and Survival. PloS one, 3(2), 2008. 9, 22
Jun Z. Li, Devin M. Absher, Hua Tang, Audrey M. Southwick, Amanda M. Casto, Sohini Ra- machandran, Howard M. Cann, Gregory S. Barsh, Marcus Feldman, Luigi Cavalli-Sforza, et al. Worldwide Human Relationships Inferred from Genome-Win Patterns of Variation. Science, 319 (0066):1100-1100 104, 2008. 9, 22
Malik Magdon-Ismail. NP-Hardness and Inapproximability of Sparse PCA. Information Processing Letters, 126:35-38, 2017.
Michael W. Mahoney and Drineas. CUR Matrix Decompositions for Improved Data Analysis. In Proceedings of the National Academy of Sciences, pp. 697-702, 106 (3), 2009. 1
Baback Moghaddam, Yair Weiss, and Shai Avidan. Generalized Spectral Bounds for Sparse LDA. In Proceedings of the 23rd International Conference on Machine learning, pp. 641-648, 2006a. 1, 3
Baback Moghaddam, Yair Weiss, and Shai Avidan. Spectral Bounds for Sparse PCA: Exact and Greedy Algorithms. In Advances in Neural Information Processing Systems. pp. 915-922, 2006b.
Cameron Musco and Christopher Musco. Randomized block krylov methods for stronger and faster approximate singular value decomposition. In Advances Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems, pp. 1396-1404, 2015. 6
Dimitris Papailiopoulos. Alexandros Dimakis, and Stavros Korokythakis. Sparse PCA through Low-rank Approximations. In Proceedings of the 30th International Conference on Machine Learning, pp. 747-755, 2013. 3, 8, 9,2,33
Haipeng Shen and Jianhua Z. Huang. Sparse Principal Component Analysis via Regularized Low Rank Matrix Approximation. Journal of Multivariate Analysis, 1915)1011 2008. 3
Bharath K. Sriperumbudur, David A. Torres, and Gert R.G. Lanckriet. Sparse Eigen Methods by D.C. Programming. In Proceedings the 24th International Conference on Machine Learning,
pp. 831-838, 2007. 3
Martin J. Wainwright. Lecture Notes, 2015. https::/www.stt stat wew.teee..eeeeeee "m hwtpjainnstat21////o/m http://cpppp.iii.ooooaaaaaaaaaaa pdf. 19
David Woodruff. Sketching as Tool for Numerical Linear Algebra. Foundations and Trends in
Theoretical Computer Science, 10(1-2), 2014.
Ganzhao Yuan, Li Shen, and Wei-Shi Zheng. A Decomposition Algorithm for the Sparse Generalized Eigenvalue Problem. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2019. 8, 23
11

---

Under review as a conference paper at ICLR 2021
Xiao-Tome Tong Yuan and Tong Zhang. Truncated Power Method for Sparse Eigenvalue Problems. Journal of Machine Learning Research, (4(Ap)):899-92)) 2013.
Dimitrios Zeimpekis and Efstratios Gallopoulos. TMG: A MATLAB Toolbox for Generating Term- Document Matrices from Text Collections. In Grouping Multidimensional Data, pp. 187-210. Springer, 2006. 9, 23
Hui Zou and Trevor Hastie. Regularization and Variable Selection via the Elastic Net. Journal of the Royal Statistical Society: Series B, 67(2):301-3200 2005. 3
Hui Zou, Trevor Hastic, and Robert Tibshirani. Sparse Principal Component Analysis. Journal of Computational and Graphical Statistics, (5(2):665-286 2006. 3, 8, 9, 23, 24, 25
12

---

Under review as a conference paper at ICLR 2021
Appendix to Approximation Algorithms for Sparse Principal Component Analysis
A APPENDIX
A.1 SPCA VIA RANDOMIZED MATRIX MULTIPLICATION. PROOFS
First, we prove two lemmas that are crucial in proving Lemma 2.1.
Lemma A.1 Given matrices P E Rmxx and Q E IRaX,, let S E R be constructed using Algorithm Then,
E [ = =(PQ);
P Q
var (PS" =
Σ
Pl
Oi:
Pk
k=1
for any indices j, E (1,...,n).
Proof:
For any i,j E [1,...,,... n), we have that
n
Z
n Σ
P Qki
B Q) =
Σ
2 RK
Qki
=
Σ
P sk
Qkj
n1 Pa Quj (PQ)
since Z3 _ Zg ~ Ber(p) and thus E [
= E [Zx]
d = denotes equality in distribution. = Pk. where
By the independence the Zy s, and noting that var (Zx) Pa (1 - -PR), we have that
Z
(PS
P ck
Qki
Pk
Rm
P Qui
1 Pk
var
Z3
P Q
Pk
Pk
k=
P
Q
12
P Q
pk
Lemma A.2 Given matrices P E Rmxx and Q E Raxp, let E Ruxn be constructed using
Algorithm 1. Then.
Q
E [ -
QQ.
Pi
1=1
Here P.; and Qi. the i-th column of P and i-th of respectively.
Proof:
Using Lemma A.1, we have that
E
PPQ - PS Q =
E ((PQ) -PPS
var (PS*Q)
P
P2
Q
ik
Σ
p2 Q
Pk
13

---

Under review as a conference paper at ICLR 2021
m
Σ
P
Σ 2i
Pkc
P Pakil 22 Qu: l1
Σ
[P Ae 11
||
12.
pk
km
Proof of Lemma 2.1: From Lemma A.2,
[P -1QQ:1112
7
P
PQ
P QQ-11
Pr
P
|Q 12
P.1. ||Q
Σ
P
12
1Q:112
spi
{ispe> 1/3
P.il 11
sp:
sp:
We conclude the proof by setting p0: 1P..1E//IP13..
Proof of Theorem 2.2. Proof: In Lemma 2.1, let P X and Q xx to get
E [|XSxx Xx
[ X
x
x
X
The last inequality follows from ||x||2 I Moreover, by Markov's inequality, with probability at least 8/4,
2 h
Xx'
X
X
Let Sx". Taking square roots of both sides of the above inequality, and applying the triangle inequality on the left hand side of the above inequality, we get:
X|x 1/2 |XSX|12|= S
|
VK
= XXSX| 2 2 V
|| r
Vk
2
2
|XSX 2 2 ~
X
2/2
|
(10)
k
Vi
Note that 2" |Xx" I|.. Ignoring the non-negative term in eqn. (10), we conclude
2c
XSX 2 Z
Vz
|X|FF
(11)
V
Next, using sub-multiplicativity on the left hand side of eqn. (11),
XSX 1/= <
X
|
|
(12)
where v E R is the top right singular vector of XS. Letting x be the i-th entry of x", we have,
2
E E[|
( E
x
x=1
(13)
Pi
Pi
since E [Z?]] =P;- Using Markov's 's inequality, with probability least 1/2,
||*|| 2.
(14)
Conditioning on this event, we can rewrite eqn. (12) follows:
XS* 2 2|[Sv|2 (Sv) x X(Sv) = =2y Tx XXX.
(15)
14

---

Under review as a conference paper at ICLR 2021
Combining eqns. (11) and (15), conclude
Xy
Vz* |
Vk
Using X A and Tr(A) |X| p concludes the proof of eqn. (4). Finally, following the lines of eqn. (13), we can prove
E(lly ( = B( [Sv] 1.
To conclude the proof of the theorem, notice that the failure probability at most 1/4 + 1/2 3/4 from union bound on the failure probabilities of eqns. (9) and (14).
0
A.2 SPCA VIA THRESHOLDING: PROOFS
We will use the notation of Section 3. For notational convenience, let ,ly...;0 On. entries of the matrix Σ E Rexn, i.e., the singular values of A.
be the diagonal
Proof of Lemma 3.2: Let U,,1 E (n-(n-)) be matrix whose columns form a basis for the subspace
perpendicular to the subspace spanned by the columns of Ug. Similarly, let Σ IR (=( be the diagonal matrix of the bottom n. l singular values of A. Notice that U = [Up U2 and
E Σg 0; 0 Σe,1]: thus,
U
/ + Ur,1 E/
By the Pythagorean theorem,
U. U
x
U2
2
Using invariance properties of the vector two-norm and sub- sub-multiplicativiiyy we get
Σ 2 U x
Σ
|
2
| x
We conclude the proof by noting that
1
x
x11
= UUUT x= x
T
Ax and
1/2
1
Σ = Tr(A)
0e+1
-
8
i=1
The inequality above follows since 00 2 de 2 01 2 on. We conclude the proof by setting e = 1/c.
Proof of Theorem 3.1. Let R fi fi...... 1 | } be the set of indices of rows of U. (columns of U that have squared norms at least c/k and let R be its complement. Here [R] denotes the cardinality of R be the sampling matrix that selects the columns of U, whose indices are in the set R and let R. E (x be the sampling matrix that selects the columns of U. whose indices are in the set R. Thus, each column of R (respectively R2) has a single non-zero entry, equal to one, corresponding to one of the [ (respectively [R]) selected columns. Formally, R,,, 1 for all
the set R and R R (l..... n). Let R E gnx
1, ... RR). while all other entries of R (respectively R_) are set to zero; R. can be defined analogously. The following properties are easy to prove: [ R I; R | R = 0. Recall that x" is the optimal solution to the SPCA problem from eqn. (1). We proceed as follows:
R... + R I,; R R I; R
1
U. x
2
U. (RR T
+ R
R
=
2
2 +2
2
2
U. RR. xx
Σ
R. R x
2
2
2 U
1
RR T xx
| +
20
U
R
R
(16)
The above inequalities follow from the Pythagorean theorem and sub-muliiplicativity. We now bound the second term in the higd side of the above inequality.
R x 2 =
15

---

Under review as a conference paper at ICLR 2021
I(U, R.) Ro)lll22 ((. S
(RIx)))
& V Ve.
| x
R
(17)
In the above derivations we use standard properties of norms and the fact that the columns of U that have indices in the set R have squared norm at most e/k. The last inequality follows from [R] [RIx S |xx 1| S Vž, since x" has at most k non-zero entries and Euclidean norm at most one.
Recall that the vector y of Algorithm maximizes [ dimensions (including Rx") and thus
2U RxI12 over all vectors x of appropriate
E
E Ryl22
U
R. x
(18)
Combining eqns. (16), (17), and (18), we get
Σソ U. ||
U 212 Tr(A).
(19)
In the above we used z = Ry (as in Algorithm 3) and 01 S Tr()). Notice that
U, U U E , 2= / 2
and use the Pythagorean theorem to get
| 2 | Σ
1 1.,
U 1 2 1
= | / 2122
Using the unitary invariance of the two norm and dropping non-negative term, we get the bound
E D U 2 <
Σ
UT
(20)
Combining eqns. (20) and (19), we conclude
1
Σ
< | U
21 + TT(A).
(21)
We now apply Lemma 3.2 to the optimal vector x to get
- Tr(A) <
Σ
/2 U.
xt
Combining with eqn. (21) we get
1 2" Tr(A).
3
z Az
2
In the above we used 1 2 Z Az and
|ز x (xx) Ax* =
A.3 SPCA VIA A SEMIDEFINITE PROGRAMMING RELAXATION: PROOFS
We start with lemma arguing that the sparsification procedure of Algorithm 5 does not significantly distort the Eg norm of the input/output vectors.
Lemma A.3 Let y and be defined in Algorithm 5. Iflly 1| S a, then, with probability at least 15/16.
16a²
||
y11 <
Proof: Notice that
n
n
n
Iyy
E v12
Σ
1
y/
Pi
Pi
which is at most a² from our assumptions. The lemma follows by Markov's inequality.
To prove Lemma 4.3, we start with the following consequence of the triangle inequality, e.g., Lemma 3 from (Fountoulakis et al., 2017).
16

---

Under review as a conference paper at ICLR 2021
Lemma A.4 Let y and be defined in Algorithm 5. Then
T Ay 2. Az| S y ( x)] + ((y - AAy A(y - -2)l-
Proof: This Lemma from (Fountoulakis et al., 2017).
We now proceed to upper bound the two terms in the above lemma separately. The following lemma bounds the first term.
Lemma A.5 Let y and z be defined as in Algorithm 5. If lly I and Iylle S B, then, with probability at least 15/16, ly A(y z)| S 408/v7.
Proof: Recall that Then,
wi/p. with probability pt and zero otherwise, for all 1...n.
E [y T - Σ -1) ve
( (A..)))
Using |A,-y| I
1yyll2 S 33 (from our assumption on the & norm of y as well our assumption that the rows/columns of A have unit norm), it follows that
V
8 ly
2
E yy A(y 2))
pi
The lemma follows from Markov's inequality.
The next lemma provides upper bound for the second term in the right hand side of Lemma A.4.
Lemma A.6 Let y and z be defined as in Algorithm If llyll, I a and ||yl| S B, then, with probability at least 16/16,
64a
96a,33
((y - 2)
AAyy A(y 2))
Proof: Let Ce Pi
with probability ps and zero otherwise for all i 1...n. Then,
E ((y-2) T
v)) A(zz-v)) Σ
Anc [[[ (a) (1) 6.) Ca)l a.b.c,d
We immediately have E - G.] 0. Thus, if any of the indices b, c, d appears only once in the above summation, then
E [11 (1 G1) (1) - -Ca)]
Let
B, = Σ
2 2. E [11 (1 - -,,
axb
B₂ = Σ
A 0.a
1 2 [(1 21 G.
azpb
B3 Σ
Aa. [[1 21 -
aybb
B4 = Σ
2 a.a 2
ME [11 (1--60 Co)
a=1
now follows that
E [(y (
-)- z) A(z y)) = Σ B
1=1
Using |Ars < for all i,j. we can bound B1, Bz, and Bg by
max [B,} 112
< Σ 2 (1 Sa) Σ
v E (1 - (6)
==1
17

---

Under review as a conference paper at ICLR 2021
Using E [1 [
for all we get
Pi
2
||
2
a²
max (B,})= 1
Σ
1,233
Pa
s
where the inequality follows by |yyll S a. To bound B4, use |yy 1|
S a
and
|yyl
S 3, to get
1
4
6
B1
A
[[1
Sa)
E
[1
Sa)
Σ
Vn
+
4 Pa
Pd
pa
Pa
q=1
6 ||yll
I s3
6 lll
lyyl2
a²
60,33
Vu
ya
s3
8
The last inequality follows from properties of norms, namely that |yl1 S 11
Thus,
a
3a
6a:33
4a
+ 60,33
E
[(y - 2) A(z y)) A(z--))
33
82
Using Markov's inequality, we conclude
4
6
Pr (z - y) A(z
x)) A(z-
4V
+
g2
16
Observe that Lemma 4.2 follows immediately from Lemma A.4, Lemma A.5, Lemma A.6, and union bound.
Proof of Lemma 4.2:
D
The next two lemmas bound the & and é norms of the vectors y; for all i = 1... M. We will bound the norm of single vector y, (we will drop the index) and then apply a union bound on all M vectors.
Lemma A.7 Let y be defined as in Algorithm 5. Then,
1 M
Pr ||yl| 2 2
2 log M
Proof: Let Z UΣV the singular value decomposition of and let o, = Σ,,, for 1... n, be the singular values of Z. Since Tr(Z) = 1, it follows that >x 1 and also o, S 1 for all i 1...n. n. Additionally,
E i=1
0: 11
(22)
1=1
Then.
gE'Vig-/EV 1ly1 - = z Zg - EV gl.
The rotational invariance of the Gaussian distribution implies that ~ h2, where h is random vector whose i-th entry h; satisfies h; N/0 of). Hence,
E
1
E b
Σ 1=1
of 11
Now, from Markov's inequality, for any C > O,
I|yll2 2 + log M
Pr Cllillla 2 LogClloMMi)
E [exll]] ca
Pr
MC/t
Then.
c Clylla = =E [1bbIl2]
18

---

Under review as a conference paper at ICLR 2021
2
cx
C dr
2TT
- -vtr, (vt,-- Cry )g))² dr
2
0
i=1
2
c
0
e -1- dl
II
0 exp
Σ c² c
i=1
Using eqn. (22), we get
E [[alllla]
I 19
Setting C = 2f and S 1, we get
log M
e
Pr
|yll2 +
< M
ecM
Setting & Vlog M concludes the proof.
Prior to bounding the & norm of y. we present measure concentration result that will be useful in our proof. First, recall the definition of L-Lipschitz functions.
Definition A.8 Let) f R - R be any function. If lf(x)) - f(y)lla S L L|Ix - for all x, y E IR", then is L-Lipschitz
Theorem A.9 (Gaussian Lipschitz Concentration) (Wainwright, 2015) Let f an L-Lipschitz function and let g E R be vector of i.i.d. Gaussians. Then f(x) is mab-Gaussian with variance L and. for all 2 0,
Pr - E[f(x)) 2 2 2e -f/222
Lemma A.10 Let y be defined as in Algorithm Then,
| 2 k(1 + 2 /lgg M)] S M
Pr
Proof: Since gj N0, 1) for all 1. .... n, the 2-stability of the Gaussian distribution implies that
2 Z
2 |
E
Zg
Z
Σ
112
12122
1=1
i=1
Let f(x) = ||Zx|| The triangle inequality implies that
"
11
[Zx||]
|Zyll 1
Σ
Zux
Z,,yl Σ
Zin (x y)
Thus, by Cauchy-Sean Schhy-
Z ||
|Zy
S
1 11 x
and f(x) is |||||, Lipschitz". Using Theorem A.9,
Pr
|
2
2
< 2e- -,/1
1.2
(> |Aco 1)),)))8)
'Recall that the Lp., norm of A is ||Al| pm
e.g. | ||
19

---

Under review as a conference paper at ICLR 2021
for all 2 0. Setting = 2/log M and noting that | 1,2 $
[ZI11,1 I k. get
1
Pr |ylli 2 k(1 + log
1og M)
M
Proof of Lemma 4.3: Using Lemma A.7 and Lemma A.10, we conclude that ||yll I and
l|yll2 S 3 both hold with probability at least
2 Using Lemma A.4, we get M
ly Ay z Az| S 2 A
z)1 (ly 2) A(y =)|-
Since ly T A(y z)] S
408 vs
with probability at least 15 (by Lemma A.5) and
16
64a g2
96a,33
(ly z) 2) A(y 2)
450n
with probability at least
16 (by Lemma A.6), setting 8
we get
Ay < x Az +
with probability at least 14/16 - 2/M 2 3/4, since M 2 16.
We now prove an inequality that was used eqn. (5) compare y T Ay and Tr(AZ).
Lemma A.11 Let Z, E R=x* be PSD matrices and Tr(AZ) I Tr(A)) for some 2 1, where Z₁ is the best rank-1 approximation of ZZ. Then,
Tr(ZAZ) 2 72. AZZ)
Here 72 = (1-(1- 0 (Z) with 01 (Z) and kZ)) being the top singular value and condition number of ZZ respectively.
Proof: For simplicity of exposition, assume that rank(Z)) = n. Let Z = UΣU be the SVD of Z.
Suppose = -U1 U1.1) and E = / 0 such that we have Z = -UE,UU and Z =
Σ1,1
UL E UT,1- As Z is the best rank-1 approximation of Z, we have Z, 7, = Z Z = Using this, we rewrite Tr(ZAZ) as the following
Tr(ZAZ) = Tr ((Z) + (Z, + ZI.1))
Tr(Z,AZ,) + Tr(Z,, (Z + AZ( + IZ
(T(Z)AZ,) + Tr(Z,, AZ + Tr(AZ, Zi) + (I(AAiZ,,.. TT1Z,AZ,) +
(23)
where the third equality follows from the invariance of matrix trace under cyclic permutations and the last step due to Z, Z1, Z
Next, we rewrite TrAZZ) as
Tr(AZ) Tr(AZ)) TT(AZ,,)) Z + Tr(AZ,,)) = ZHZZZ,AZ, + Tr(AZ,,)) o1 (Z)) (Z, TI(AZI,))
(24)
where Z is the pseudo-inverse of Z₁ and we have used the fact Z, Z Z, Z, and the last inequality follows from the von Neumann's 's trace inequality. Now, noting that a (Z)) o1( along with the fact that o AZ, < AZ,Z, applying eqn. (23), we have
Tr(AZ) (Z) (TiZAA)) Tr(Z, AZ) Tr(AZ),.)
(25)
20

---

Under review as a conference paper at ICLR 2021
Next, we will show that 1
o (Z)-T((()),, First, note that Σ1, > on (Z) I n-1-
as on ()) S 0, (Z) for all i = 2,..., n. Therefore, pre- and post-multipyinng both sides by Σ 1, 1 we
further have Σ > (Z) Σ E Again, pre- and post-multillying both sides by 11,1 and U 1,1+
we have:
7 UI EJ UT + (Z) U E1, UT (Z) ZI.1.
(26)
As the matrix A is also PSD, it has PSD square-root A such that A = A A/2,, Now, pre- and post-multillingg both sides of eqn. (26) by A we have
A 1 Z
A > on (Z) A Z A 1/2
(27)
Next, we rewrite Tr AZ follows
Tr(Zi, AZI, ) = Tr(AZ) = Tr (ZA A'") = Σ e A A A Ave. i=1
2 (Z) Σ e A A 1=1
A ee = on (Z) Tr A A 1
= o (Z) Tr (AZ
(28)
In the above, e,, ..,..., en E Rn are the canonical basis vectors and we have used the invariance
property of matrix trace under cyelic permutations. Finally, the inequality in eqn. (28) directly follows
from eqn. (27), as eqn. (27) boils down to x A A A 12x 2 on (Z) x A A 1/xx for any vector x 7 0.
Next, combine eqns. (25) and (28) and replacing (Z)
oz(Z) (Z) to get
Tr(ZAZ)
1
Tr(AZ)
Tr(AZ).))
(Z)
(Z)
Tr(ZAZ)
1
1
(Tr(AZ) Tr())))
01 (Z)
k(()
Tr(ZAZ)
k(Z) ) (1-2) Tr(AZ),
1
(29)
(Z)
where the equality is holds as Tr(AZ) = Tr(AZ,) + Tr(AZ,,)) and the last inequality is due to assumption that Tr(AZ) < TrrAZ.). This concludes the proof.
D
To finalize our proof, we use the following result of (Avron & Toledo, 2011) for estimating the trace of PSD matrices.
Theorem A.12 (Avron & Toledo, 2011) Given a PSD matrix A E Raxn, let M 80/22, Lel gi (for i == -=...MM M) be standard Gaussian random vectors. Then with probability at least 7/8,
1 M
Tr(A)
Σ g.
Ag S. Tr(A).
M
We are now ready to prove the correctness of Algorithm 5 by establishing Theorem A.
Theorem A.13 Let Z be an optimal solution to the relaxed SPCA problem eqn. (2) and Assume that Tr (Z) < a Tr( (AZ)) for some constant a 2 1, where Z₁ is the best rank-1 approximation of Z. Then, there exists an algorithu that takes as input a PSD matrix A E gaxn, an approximation parameter > 0, and a parameter k. and outputs a vector z such that with probability at least 5/8,
E [|2]1oo] s 450c² 33
1z/12 S
2 A 2 72 -2)
2
Vs
In the above, a = k(1 + 2/ M), 3 = 2 M, and M = 80/2, and zL (1 number of respectively.
R (1- a (Z) with 01 (Z) and zZ) being the top singular value and condition
21

---

Under review as a conference paper at ICLR 2021
Proof: Consider Algorithm 5 and let Z" be optimal solution to the SPCA semidefinite relaxation of eqn. (2). Then, as already discussed, (x*) I "Tr(AZ)). where x" is the optimal solution to the SPCA problem of eqn. (I). Then, using Lemma A.11, I, it follows that
7z. Tr(AZ*) x Tr ((Z)) ((Z))AZZ)
Applying Theorem A.12 to the matrix (Z*) "Z"" and using our choice of in Algorithm 4, get
M M
y Ay
2
ZZ')
AZ' gi 2 - c) ) AZZ),
1=1
with probability at least 7/s. By Lemma 4.3, have y Ay S z Az e with probability least 8/1. Thus, with probability at least 8/s,
(1- =" (1 c)72 (x') Ax' Sz Az+
To conclude the proof, we need to bound the & norm of the solution vector z. Let be the event that [[ 1 S k(1 2ol M) and || |Z - 2 M for all = 1... From Lemma A.7 and Lemma A0 and the union bound, we have Pr [8] > 1 - Conditioned on E, Lemma A.3 implies that, with probability at least 15/16,
16 (1 + 2og My²
lyy - 2
Therefore, with probability at least inequality gets
15 16
2
2 (since M 2 16), an application of the triangle
M
4k(1 2o M)
| S
y
yll S 2 log M
vs
Using our chosen values for a and B concludes the proof.
B ADDITIONAL NOTES ON EXPERIMENTS
In addition, we normalize the outputs of Algorithm 2 and Algorithm 5 by keeping the rows and columns A corresponding to the nonzero elements the output vectors and then getting the top singular vector of the induced matrix and padding it with zeros. The above two considerations make our comparisons fair in terms of function f(y) (see Section 5 for the definition of f(y)). For Algorithm 3, fix the threshold parameter to 30 for human genetic data, as well for the text data; we set L = 10 for the gene expression data. Finally, for Algorithm 5, we fix M (the number of random Gaussian vectors) to 300 and experiments were implemented on a single-core Intel(R) Xeon(R) Gold 6126 CPU @ 2.60GHz.
Python's 's cvxpy package to solve eqn. (2). All the
B. REAL DATA
Population geneties data. We use population genetics data from the Human Genome Diversity Panel (Consortium, 2007) and the HAPMAP (Li al., 2008). In particular, we use the 22 matrices (one for each chromosome) that encode all autosomal genotypes. Each matrix contains 2,240 rows and a varying number of columns that is equal the number of single nucleotide polymorphisms (SNPs, well-known biallelic loci of genetic variation across the human genome) in the respective chromosome. The columns of each matrix were mean-centered preprocessing step. See Table for summary statistics.
Gene expression data. We also use lung cancer gene expression dataset (GSE10072) from from the NCBI Gene Expression Omnibus database (Landi et al., 2008). This dataset contains 107 samples (58 cases and 49 controls) and 22,215 features. Both the population genetics and the gene expression datasets interesting in the context of sparse PCA beyond numerical evaluations, since the sparse components can be directly interpreted to identify small sets SNPs or genes that capture the data variance.
22

---

Under review as a conference paper at ICLR 2021
Text classification data. We also evaluate algorithms on text classification dataset used in (Fountoulakis et al., 2017). This consists of two publicly available standard test collections for ad hoc information retrieval system evaluation: the Cranfield collection that contains 398 abstracts of aerodynamics journal articles and the CISI (Centre for Inventions and Scientific Information) data that contains 1, 460 information science abstracts. Finally, using these two collections, sparse, 2, 858 x 12, 427 document-term matrix was created using the Text-to-Matrix Generator (TMG) (Zeimpekis & Gallopoulos, 2006), with the entries representing the weight of each term in the corresponding document. See Table 5 for summary statistics.
B.2 SYNTHETIC DATA
We also use a synthetic dataset generated using the same mechanism as in (Fountoulakis et al., 2017). Specifically, we construct the m X matrix X such that X = UEV + E,. Here, E, is noise matrix, containing i.i.d. Gaussian elements with zero mean and we set a 10 10-3; U E R
a Hadamard matrix with normalized columns; Σ = Σ 0) E Rmx" such that Σ E R a diagonal matrix with Σn = 100 and E; = e- for = 2,..., m, V E R such that = (0),, where V E Ra** is also a Hadamard matrix with normalized columns and (o) = G( in + 1.0) G i2 + 1,
G f//1 11 + 1,0).
composition of 2 Givens rotation matrices with is う + 2k- for k 1, ..... 4. (, E n be Givens rotation matrix, which rotates the plane i j by an angle 0. For 0 0.27T and = 212, the matrix (0) rotates the bottom 3/2 components of the columns of V, making half of them almost zero and the rest half larger. Figure 2 shows the absolute values of the elements of the first column of the matrices V and V.
Here
B.3 ADDITIONAL EXPERIMENTS
X
X
A
Xa
X
Xg
Xs
Xn
10
PVE
2
spea-r (Algorithm 2)
0
0
0.51
0.51
0.50
0
0
0.48
39.0%
1164.2
spea (Algorithm
0.50
0.50
0.50
0.50
0
0
0
0
39.5%
1161.0
apoa-adp (A (Algorithm 5)
0
0
0
0.50
0.50
0.50
0.50
0
0
40.9%
1201.0 .0
dec (Yuan al., 2019)
0
0.50
0.50
0.50
0.50
0
0
40.9%
1201.0
expeca (Beck Vaisbourd, 2016)
0.50
0.50
0.50
0. 50
0
D
0
01
39.5%
1161.0
spea lowzank (Papaillopoulos ot 2013)
0
0.50
0.50
0.50
0.50
0
0
40.9%
1200.9
apoa (Zou al., 2006)
0.50
0.50
0.50
40.9%
1201.0
dapco (d' Aspremont al., 2007)
0
0.50
0.50
0.50
0.50
0
0
40.9%
1201.0
Table 2: Loadings, % variance explained (PVE), and the objective function value for the first principal component of the artificial data.
Additionally, in order to further explore the sparsity patterns of the outputs of our algorithms and how close they compared standard methods, further apply our methods on simulation example proposed by (Zou et al., 2006). We describe them below:
Artificial Data of (Zou et al., 2006). In this example, three hidden factors Vi, V2, and V/3 created in the following way:
V1 (00,200)) V/2 300,3000)
V3 0.3 - 0.925 V2 te, E~ 10, and
V1, V2, and are independent.
Next, we create 10 observable variables X1, X2.... X100 in the following way:
X; Vi + En, En NN(0 1), i 1
X; V/2 + 6:2, 6n2 ~ .V(0,)) 5=5,6,7,8,
X; V/3 + Ea, 6a .N(0,)) 1), 1=9,10,
Eig are independent i= 1,2, ..., 10; j = 1,2,3.
23

---

Under review as a conference paper at ICLR 2021
We take A to be the exact covariance matrix of (X1 X2 X10) to compute the top principal component. As the first two factors i.e, V1 and V/2 are associated with four variables while the last one L.., V3 is associated with only two variables and noting that all the three factors V,, V2. and V3 roughly have the same variance, V1 and V/2 are almost equally important, and they are both significantly more important than V3. Therefore, for the first sparse principal component, the ideal solution would be to use either (X1, X2, X3, or (Xa, Xo. X7, Xs). Using the true covariance matrix and the oracle knowledge that the ideal sparsity (Zou et al., 2006) as well as the SDP-based algorithm of (d' Aspremont et al., 2007). We found that while two of our methods, namely, spea-d (Algorithm 3) and spca-sdp (Algorithm 5) are able to identify the correct sparsity pattern of the optimal solution, spca-r (Algorithm 2) wrongly includes variable X instead of Xs. possibly due to high correlation between 1/2 and V/3 (see Table 2 for details). However, the output of the spca -d is much more interpretable, even though it has slightly lower PVE than spca-r.
4, we apply our algorithms and compare it with spca of
In additional experiments the large datasets, Figure 2b shows the performance of various SPCA algorithms on the synthetic data. Notice that the performance of the max maxcomp heuristic is worse than spca well synthetic data. In particular, turning the bottom elements of V into large values guarantees that these would not be good elements to retain the construction of the output vector in maxcomp, as they fail to capture the right sparsity pattern. On the other hand, algorithms perform better than comparable to spca. Similar to the real data, performance of spca-sdp closely matches with that of dec, ewpca, and spca lowrank. In Figure 3, we demonstrate how our algorithms perform on CUR 3 and CIIR 4 of the population genetics data. We see similar behavior as observed for CIR and CHR Figures la-lb. Table 3, we report the variance fyy captured by the output vectors of different methods for the text data, which again validates the accuracy of our algorithms.
algorithms. This is quite evident from the way constructed the
0.020
0.6
0.5
dec spoaloerark
0.4
0.010
0.005
rotation without rotation
1000
2000
3000
4000
0
250
500
750
1000
1250
1500
1750
index
Spansily: leel
(b)
Fig. 2: Experimental results on synthetic data with m 2 and 12; (a) the red and the blue lines are the sorted absolute values of the elements of the first column of matrices and V respectively. (b) f(y) vs. sparsity ratio.
PC3
5pc3
naxcomp
2pca-d
Spca-F
Spc3-5dp
CNSCa
doc
Spico-JoMrans
100
0.497
0.0174
00009
0.1385
000297
01144
0.1319
0.177
0.1566
100
0.497
1.884
(2140)
0.2948
011681
0.354
0.337
0.31I
03323
0.280X
04214
5,00
0.4397
0.4257
0.4411
0.4422
0.4276
0.4413
0,4558
0.4221
0.4560
10,000
0.1397
0.9962
0.4512
0.3505
0.3412
0.4985
0,4559
0.4343
0:4571
Table 3: Text data: f(y) vs. the sparsity parameter k for various SPCA algorithms.
24

---

Under review as a conference paper at ICLR 2021
spcad
spcad
spca
Spc0
0.6
0.5
spca
0.4
0.4
0.3
0.1
0.0
0.0
1000
2000
3000
4000
5000
6000
7000
1000
2000
3000
4000
5000
6000
SparsRy: (bÍ
SparsRy: lbillo
(a) CHR 34, 258
(b) CHR 4. 30, 328
Fig. 3: Experimental results on real data: f(y) vs. sparsity ratio.
Table 4: Statistics the population genetics data.
Dataset
Rows (m²
Columns (n)
Density
CHR
2,240
37,493
0.986
CHR 2
2,240
40,844
0.987
CHR
2,240
34,258
0.986
CHR 4
2,240
30,328
0.986
Table Statistics of gene expression and text data.
Dataset
# Rows (m)
# Columns (n)
Density
Gene expression
107
22,215 12,427
0.999
Text classification 2,858
0.004
Table 6: Values of 01 (Z), (z(), a, and 7/2 for various datasets.
Dataset
(Z)
x(Z)
7'2
Pit Props with k = 7
0.9999
3.95 x 1000
1.0000
0.9999
Data of Zou et al. (2006) with
0.9999
5.14 x 1040
1.0000
0.9999
CHR with l: 10, 000
0.9985
7.12 106
1.0017
0.9968
CHR 2 with k 10, 000
0.9996
6.94 X 106
1.0003
0.9982
CHR with 10, 000
0.9995
9.47 X 107
1.0005
0.9989
CIIR 4 with
10, 000
0.9998
1.27 107
1.0002
0.9991
Gene expression with & = 500
0.9913
2.05 x 106
1.0001
0.9914
Text Classification with & 5,0
0.9997
5.78 x 106
1.0001
0.9987
25


Reviews:
[{"review_id": "xAXBrLyyvks", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposed three simple algorithms for sparse principal component analysis (SPCA): a) randomized matrix multiplication; b) deterministic thresholding scheme; and c) semidefinite programming relaxation. All of the proposed algorithms look like native combinations of existing techniques and simple sparsification steps. However, it is somewhat interesting to have novel theoretical guarantees for these simple strategies whose error bounds depend on the properties of the input matrix and the target sparsity.\n\nThis paper is mathematically sound and the theoretical bound of existing SPCA tricks is also interesting to know. But the significance of the work is not clearly written. Novelty is marginal as these tricks were already there it seems. It needs clarity on improvement on the bound in terms of assumptions and tightness of existing bounds related to SPCA provided by d’Aspremont et al., 2014, Papailiopoulosetal.,2013\n\nThe experiment section is weak as the author did not compare to recent methods like d’Aspremont et al., 2014, Papailiopoulosetal.,2013, The result achieved by artificial data is not demonstrating the significance of the proposed method in comparison with the existing methods. (supplementary material)\n\nNo experimental validation is given on the bound. It will be good to verify how much tight these theoretical bounds are. \n\nThe analysis of Algorithm 2 and Algorithm 3 are based on that of full SVD whose complexity is  O(n^3).  It will be nice to know what is the running time for the proposed algorithm on a moderately large dataset.\n\nI recommended a reject as it has marginal novelty and the claims of the proposed method are not experimentally well supported in the paper.\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "wuR1dDSpbhK", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper presents three approximation algorithms for the sparse PCA problem. They are (i) randomized matrix multiplication, (ii) deterministic thresholding scheme, (iii) semidefinite programming relaxation. All the three algorithms have provable theoretical guarantees and low-degree polynomial time. \n\nI have the following comments.\n1. The authors propose a set of algorithms. It makes me feel that if any of the algorithms really works, the authors do not necessarily propose three algorithms.\n\n2. The algorithm only output a solution with expected k-sparsity. However, its variance could be large. This significantly limits it practical use since we often need a exactly k sparse solution. \n\n3. The experimental comparisons are not sufficient. The authors should compare their methods with state-of-the-art sparse PCA solvers such as coordinate-wise optimization method (The sparse principal component analysis problem: Optimality conditions and algorithms, JOTA 2016) and block decomposition method (A Decomposition Algorithm for the Sparse Generalized Eigenvalue Problem, CVPR 2019). In addition, the dimension of the data set is too small to show the effectiveness of the algorithm. For the pit props data set which only contains 13 dimensions, it is shown that the block decomposition method can find the global optimal solution.\n\n4. This paper is of theoretical research interest only. ICLR  emphasises a lot on experiments and empirical evaluation, STOC/FOCS/SODA could be a better place for the present paper. \n\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "h0Oa1Nl8dz5", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Overall, I vote for weak rejecting. The theoretical findings in this paper is presented clearly and looks solid. However, my major concern is the meaningfulness of the problem (in particular, the setting that imposes no assumption on the input matrix) studied in this paper (see cons below). Hopefully the authors can address my concern in the rebuttal period. \n\n \n##########################################################################Pros: \n\n \n1. The paper studies the sparse PCA problem. While no assumptions is imposed on the input covariance matrix, the authors give theoretical guarantees on the optimality gap (i.e. the gap between the optimal value and the objective value achieved by the solution) as well as the sparsity level of the solution. \n\n2. The authors provide numerical experiments to illustrate their theoretical findings for the three presented algorithms.\n\n \n##########################################################################\n\nCons: \n\n \n1. The sparse PCA problem studied in prior literatures aims to estimate the leading eigenvectors of a covariance matrix under the high dimensional setting, when the leading eigenvectors are assumed to be sparse. This paper imposes no constraint on the input covariance matrix. This leads to my major concerns as follows:\n\n(i) Is the problem still meaningful when we do not assume that the leading eigenvectors of the input covariance matrix are not sparse? If this is the case, I think it is not reasonable to seek for a sparse estimate of the leading eigenvector.\n\n(ii) In the case when the top eigenvector is sparse, the theoretical findings in this paper does not provide guarantees that the output of the three algorithms are reliable approximation of the top eigenvector. For Algorithm 2 and Algorithm 3, the results in Theorem 2.2 and 3.1 indicates that there is a nonnegligible optimality gap, which can be as large as half the optimal value. For Algorithm 3, the results in Theorem 4.1 also relies on the properties of the SDP solution (through constants $\\alpha$ and $\\kappa(Z)$), which also might result in a large optimality gap. \n\nI understand that the main purpose of this paper is to study SPCA problem under no assumptions on the covariance matrix, so it might not be necessary to address (ii). In that case, I suggest the authors provide more explanations for the meaning of solving SPCA problem (equation (1)) for general input covariance matrix, as well as the implications of the theoretical guarantees (e.g. the reason why the outputs of the three algorithms are reliable/meaningful estimates, given the fact that they satisfy the theoretical guarantees shown in the theorems) to help address (i).\n\n2. In Theorem 4.1, the exceptional probability is 3/8. It might be better to set the exceptional probability to be $\\delta$, and show the dependency of other quantities on $\\delta$.", "extra_scores": {}, "extra_text": {}}, {"review_id": "PDimeKVgz7F", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposed three simple algorithms for sparse principal component analysis (SPCA): a) randomized matrix multiplication; b) deterministic thresholding scheme; and c) semidefinite programming relaxation. \nAll of the proposed algorithms look like native combinations of existing techniques and simple sparsification steps. However, it is somewhat surprising that such simple strategies have reasonable theoretical guarantees whose error bounds depend on the properties of input matrix and the target sparsity. \n\nI have some comments as follows:\n\n1. The analysis of Algorithm 2 and 3 are based on full SVD whose complexity is $O(n^3)$. The authors also mentioned that we can use iterative methods to replace full SVD in practice.  Hence, I think it is necessary to provide the theoretical analysis for the proposed frameworks with inexact SVD and establish the error bounds contain the error from approximation of SVD. Similarly, it is also desired to give the results for Algorithm 5 when we only obtain an approximate solution of problem (2).\n\n2. It is prefer to test the proposed algorithms on large scale datasets, which could make the paper be more convincing to machine learning community. Additionally, since  $O(n^3)$ is unacceptable for large $n$, it is necessary to implement the ideas of this paper with inexact SVD and report the corresponding empirical results.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes three algorithms for the sparse PCA problem, where one imposes the additional constraint that the vectors have a small number of non-zero entries. The proposed algorithms run in polynomial time and achieve provable approximation guarantees on the accuracy and sparsity. The reviewers identified the following strengths of the contributions: the algorithms are simple and have different strengths; the theoretical results are sound and perhaps even surprising; the presentation is clear. The reviewers identified the following weaknesses of the contributions: the running times of the proposed algorithms are high and they may not scale to large datasets, which significantly limits their application to machine learning datasets; the experimental evaluation is insufficient and it does not compare with some of the state of the art algorithms; the algorithmic novelty is limited. After weighing these strengths and weaknesses as well as evaluating the paper relative to other ICLR submissions, I recommend reject.
Decision: Reject