A Randomized Rounding Algorithm for Sparse PCA

arXiv:1508.03337v5 [cs.DS] 22 Nov 2016

Kimon Fountoulakis∗

Abhisek Kundu†

Eugenia-Maria Kontopoulou‡

Petros Drineas§

Abstract
We present and analyze a simple, two-step algorithm to approximate the optimal solution of the
sparse PCA problem. Our approach first solves an ℓ1 -penalized version of the NP-hard sparse PCA optimization problem and then uses a randomized rounding strategy to sparsify the resulting dense solution.
Our main theoretical result guarantees an additive error approximation and provides a tradeoff between
sparsity and accuracy. Our experimental evaluation indicates that our approach is competitive in practice,
even compared to state-of-the-art toolboxes such as Spasm.

Keywords.

Sparce pca, rounding randomized algorithm.

1 Introduction
Large matrices are a common way of representing modern, massive datasets, since an m × n real-valued
matrix X provides a natural structure for encoding information about m objects, each of which is described
by n features. Principal Components Analysis (PCA) and the Singular Value Decomposition (SVD) are
fundamental data analysis tools, expressing a data matrix in terms of a sequence of orthogonal vectors of
decreasing importance. While these vectors enjoy strong optimality properties and are often interpreted
as fundamental latent factors that underlie the observed data, they are linear combinations of up to all the
data points and features. As a result, they are notoriously difficult to interpret in terms of the underlying
processes generating the data [MD09].
The seminal work of [dGJ07] introduced the concept of Sparse PCA, where sparsity constraints are
enforced on the singular vectors in order to improve interpretability. As noted in [dGJ07, MD09, PDK13],
an example where sparsity implies interpretability is document analysis, where sparse principal components
can be mapped to specific topics by inspecting the (few) keywords in their support. Formally, Sparse PCA
can be defined as1 (see eqn. (1) in [PDK13]):
Zopt =

max

x∈Rn , kxk2 ≤1

xT Ax,

s.t.

kxk0 ≤ k.

(1)

In the above formulation, the parameter k controls the sparsity of the resulting vector and is part of the input;
A = XT X ∈ Rn×n is the symmetric positive semidefinite (PSD) covariance matrix that represents all
∗

International Computer Science Institute, Department of Statistics, University of California Berkeley, Berkeley, CA, USA,
kfount@berkeley.edu
†
Speech and Biomedical Analytics Group, Xerox Research Centre India, Bangalore, India, abhisekkundu@gmail.com
‡
Department of Computer Science, Purdue University, West Lafayette IN, USA, ekontopo@purdue.edu
§
Department of Computer Science, Purdue University, West Lafayette IN, USA, pdrineas@purdue.edu
P
p
1
Recall that the ℓp norm of a vector x ∈ Rn is defined as kxkpp = n
i=1 |xi | when 0 < p < ∞; we will only use p = 1, 2
and 3. Furthermore, note that kxk0 is the ℓ0 “norm”, which is the number of non-zero entries of the input vector x.

1

pairwise feature similarities for the data matrix X and xopt denotes a vector that achieves the optimal value
Zopt in the above formulation. In words, the above optimization problem seeks a sparse, unit norm vector
xopt that maximizes the data variance xTopt Axopt . The following two facts are well-known: first, solving
the above optimization problem is NP-hard [MWA12] and, second, this hardness is due to the sparsity
constraint. Indeed, if the sparsity constraint was removed, then the resulting optimization problem can be
easily solved: the optimal vector is the top (left or right) singular vector of A and the related maximal value
is equal to the top singular value of A. Finally, we note that more than one sparse singular vectors can be
constructed using a simple deflation procedure; see Section 3.1 for details. Following the lines of [PDK13],
we will only focus on the formulation of eqn. (1).
Related work. The simplest sparse PCA approaches are to either rotate [Jol95] or threshold [CJ95] the top
singular vector of the matrix A. Such simple methods are computationally efficient and tend to perform very
well in practice (see Section 3). However, there exist cases where they fail (see [CJ95] and Section 3 here).
An alternative line of research focused on solving relaxations of eqn. (1). For example, an ℓ1 relaxation of
eqn. (1) was first used in SCoTLASS [JTU03]. Another possible relaxation is a regression-type approximation [ZHT06], which was implemented in [SCLE12]. (We will compare our approach to this method.)
Finally, efficient optimization methods have been developed for the sparse PCA problem. For example, the
generalized power method was proposed in [JNRS10]: this method calculates stationary points for penalized
versions of eqn. (1).
Despite the many approaches that were developed for sparse PCA, only a handful of them provide any
type of theoretical guarantees regarding the quality of the obtained (approximate) solution. For example,
the semidefinite relaxation of [dGJ07] was analyzed in [AW08], albeit for the special case where A is a
spiked covariance matrix with a sparse maximal singular vector. Briefly, [AW08] studies conditions for the
dimensions m and n of the data matrix X, and the sparsity parameter k, so that the semidefinite relaxation
of [dGJ07] recovers the sparsity pattern of the optimal solution of eqn. (1). Other attempts for provable
results include the work of [dBG08], which was later analyzed in [dBG14]. In the latter paper, the authors
show bounds for the semidefinite relaxation of [dBG08], in the special case that the data points are sampled
using Gaussian models with a single sparse leading singular vector. Strong compressed-sensing-type conditions were used in [YZ13] to guarantee recovery of the optimal solution of eqn. (1) using a truncated power
method. However, [YZ13] requires that the optimal solution is approximately sparse and also that the noise
matrix has sparse submatrices with small spectral norm. Finally, [PDK13] describes a greedy combinatorial
approach for sparse PCA and provides relative-error bounds for the resulting solution under the assumption
that the covariance matrix A has a decaying spectrum. It is important to note that in all the above papers
special assumptions are necessary regarding the input data in order to guarantee the theoretical bounds. Our
work here does not make any assumptions on the input data, but we do pay the cost of increased sparsity
in the final solution (see Theorem 1 for details). There are also connections between sparse approximations
and subspace learning methods, which are widely used in machine learning and data mining. Recently,
methods that enforce sparsity have been developed for Human Activity Recognition [LZW+ 16] and Image
Classification [LWT+ 16]. Moreover, some of these methods have been applied to multidimensional data
[DXXM07].
Our algorithm. We present and analyze a simple, two-step algorithm to approximate the optimal solution
of the problem of eqn. (1). Our approach first finds a stationary point of an ℓ1 -penalized version of problem
(1). Then, a randomized rounding strategy is employed to sparsify the resulting dense solution of the ℓ1 penalized problem. More precisely, we first solve:
√
(2)
Z̃opt =
max
xT Ax, s.t. kxk1 ≤ k.
x∈Rn , kxk2 ≤1

2

ALGORITHM 1: Vector sparsification
Input: x ∈ Rn , integer s.
Output: x̂ ∈ Rn , E ( kx̂k0 ) ≤ s.
for i = 1, . . . , n do
Set:
(
x̂i =

1
pi xi ,

0,

with probability pi = min
otherwise

n

s|xi |
kxk1 , 1

o

end

Notice that we replaced the constraint on the ℓ0 norm of the vector x by a (tighter) constraint on the ℓ1 norm
of x. It is important to mention that problem (2) is difficult and all we can hope in practice is to calculate
a stationary point. However, one should not discount the quality of stationary points. In Section 3 we
show that by calculating stationary points we capture as much of the variance as computationally expensive
convex relaxations. Having said that, the theoretical analysis which follows assumes that we work with the
globally optimal solution of problem (2).
Let x̃opt be a vector that achieves the optimal value Z̃opt for problem (2). Clearly, x̃opt is not necessarily
sparse. Therefore, we employ a randomized rounding strategy to sparsify it by keeping larger entries of x̃opt
with higher probability. Specifically, we employ Algorithm 1 on the vector x̃opt to get a sparse vector x̂opt
that is our approximate solution to the sparse PCA formulation of eqn. (1).
It is obvious that, in expectation, the vector x̂opt has at most s non-zero entries. We will discuss the
appropriate choice for s in Theorem 1 below.
Our results: theory. Surprisingly, this simple randomized rounding approach has not been analyzed in
prior work on Sparse PCA. Theorem 1 is our main theoretical result and guarantees an additive error approximation to the NP-hard problem of eqn. (1). For simplicity of presentation, we will assume that the
rows (and therefore columns) of the matrix A have at most unit norms2 .
Theorem 1 Let xopt be the optimal solution of the Sparse PCA problem of eqn. (1) satisfying kxopt k2 = 1
and kxopt k0 ≤ k. Let x̂opt be the vector returned when Algorithm 1 is applied on the optimal solution x̃opt
of the optimization problem of eqn. (2), with s = 200k/ǫ2 , where ǫ ∈ (0, 1] is an accuracy parameter. Then,
x̂opt has the following properties:
1. E kx̂opt k0 ≤ s.
2. With probability at least 3/4,
kx̂opt k2 ≤ 1 + 0.15ǫ.
3. With probability at least 3/4,
x̂Topt Ax̂opt ≥ xTopt Axopt − ǫ.

(3)

In words, the above theorem states that our sparse vector x̂opt is almost as good as the optimal vector xopt
in terms of capturing (with constant probability) almost as much of the spectrum of A as xopt does. This
comes at a penalty: the sparsity of xopt , which is equal to k, has to be relaxed to O k/ǫ2 . This provides an
elegant trade-off between sparsity and accuracy3 . It is worth emphasizing that one should not worry about
2
We can relax this assumption by increasing s – our sampling factor – by a factor that depends on the upper bound of the row
and column norms of A.
3
A less important artifact of our approach is the fact that the Euclidean norm of the vector x̂opt is slightly larger than one.

3

the small success probability of our approach: by repeating the rounding t times and keeping the vector
x̂opt that satisfies the second bound and maximizes x̂Topt Ax̂opt , we can immediately guarantee that we will
achieve both bounds with probability at least 1 − 2−t .
Our results:experiments. We empirically evaluate our approach on real and synthetic data. We chose to
compare our algorithm with the state-of-the-art Spasm toolbox of [SCLE12, ZHT06]. We also compare
our solution with the simple MaxComp heuristic, which computes the top singular vector of matrix A and
returns a sparse vector by greedily keeping the top k largest components (in absolute value) and setting the
remaining ones to zero. Our empirical evaluation indicates that our simple, provably accurate approach is
competitive in practice.

2 Proof of Theorem 1.1
2.1 Preliminaries
Consider the indicator random variables δi for all i = 1 . . . n which take the following values:
 1
pi , with probability pi
δi =
0 , otherwise
It is easy to see that x̂i = δi x̃i for all i = 1 . . . n. The following trivial properties hold for all i and will
be used repeatedly in the proofs: Eδi = 1, E (1 − δi ) = 0, and E (1 − δi )2 = p−1
i − 1. For simplicity of
notation we will drop the the subscript opt from xopt , x̃opt , and x̂opt in all proofs in this section.

2.2 A bound for kx̂opt k0
Recall that we will use the sampling probabilities pi as defined in Algorithm 1. By definition, pi ≤
s |x̃i | / kx̃k1 , therefore
n
n
X
X
s |x̃i |
= s.
(4)
pi ≤
E ( kx̂k0 ) =
kx̃k1
i=1

i=1

which proves the first bound in Theorem 1.

2.3 A bound for kx̂opt k2
The following lemma immediately implies the second bound of Theorem 1 by setting s = 200k/ǫ2 .
Lemma 1 Given our notation, with probability at least 3/4,
r
kx̂opt k2 ≤ 1 + 2

4

k
.
s

Proof: It is more intuitive to provide a bound on the expectation of kx̂ − x̃k22 and then leverage the triangle
inequality in order to bound kx̂k2 . Using the indicator variables δi and linearity of expectation, we get
E kx̂ − x̃k22 = E
=
=

n
X
i=1

n
X

x̃2i E (1 − δi )2

i=1
n 
X
i=1

(1 − δi )2 x̃2i


1
− 1 x̃2i .
pi

We will now prove the following inequality, which will be quite useful in later proofs:
n 
X
1
i=1


k
− 1 x̃2i ≤ .
pi
s

(5)

Towards that end, we will split the set of indices {1 . . . n} in two subsets: the set I =1 corresponding to
indices i such that pi = 1 and the set I <1 corresponding to indices i such that pi < 1. Note that for all
i ∈ I <1 it must be the case that
s |x̃i |
.
pi =
kx̃k1
We now proceed as follows:


X 1
X 1
− 1 x̃2i +
− 1 x̃2i
pi
pi
i∈I <1
i∈I =1


X 1
=
− 1 x̃2i
p
i
i∈I <1
X 1
x̃2i
≤
p
i
<1

E kx̂ − x̃k22 =

i∈I

=

X kx̃k
1 2
x̃i
s
|x̃
i|
<1

i∈I

≤

n
kx̃k1 X
|x̃i |
s
i=1

kx̃k21

k
.
s
s
√
For the last inequality we used the fact that kx̃k1 ≤ k. We now use Markov’s inequality to conclude that,
with probability at least 3/4,
4k
.
(6)
kx̂ − x̃k22 ≤
s
To conclude the proof note that, from the triangle inequality,
=

≤

kx̂ − x̃k2 ≥ | kx̂k2 − kx̃k2 |
5

and thus
kx̂k2 ≤ kx̃k2 + kx̂ − x̃k2 ≤ 1 + kx̂ − x̃k2 .
Combining with eqn. (6), after taking the square root of both sides, concludes the proof of the lemma.

⋄

2.4 Proving eqn. (3)
The following lemma states that the solution of the convex relaxation of the Sparse PCA problem is at least
as good as the solution of the Sparse PCA problem.
Lemma 2 Given our notation, xopt is a feasible solution of the relaxed Sparse PCA formulation of eqn. (2).
Thus, Z̃opt = x̃Topt Ax̃Topt ≥ xTopt Axopt .
Proof: Recall that xopt is a unit norm vector whose zero norm is at most k. Then, if we let sgn(x) denote
the vector of signs for x (with the additional convention that if xi is equal to zero then the i-th entry of
sgn(x) is also set to zero), we get
√
kxopt k1 = sgn(xopt )T xopt ≤ ksgn(xopt )k2 kxopt k2 ≤ k.
The second inequality follows since sgn(xopt ) has at most k non-zero entries. Thus, xopt is feasible for the
optimization problem of eqn. (2) and the conclusion of the lemma follows immediately.
⋄
T
Getting a lower bound for x̂ Ax̂ is the toughest part of Theorem 1. Towards that end, the next lemma
bounds the error x̃Topt Ax̃opt − x̂Topt Ax̂opt as a function of two other quantities which will be easier to
bound.
Lemma 3 Given our notation,
x̃Topt Ax̃opt − x̂opt Ax̂opt ≤ 2 x̃Topt A (x̃opt − x̂opt ) + | (x̂opt − x̃opt )T A (x̃opt − x̂opt ) |.
Proof: We start with
x̃T Ax̃ − x̂T Ax̂

x̃T Ax̃ − x̂T Ax̃ + x̂T Ax̃ − x̂T Ax̂

=

(x̃ − x̂)T Ax̃ + x̂T A (x̃ − x̂) .

≤

(7)

Next,
(x̂ − x̃)T A (x̃ − x̂)

=
≥

x̂T A (x̃ − x̂) − x̃T A (x̃ − x̂)
x̂T A (x̃ − x̂) − x̃T A (x̃ − x̂) ,

(8)

which implies
x̂T A (x̃ − x̂)

≤

x̃T A (x̃ − x̂) + (x̂ − x̃)T A (x̃ − x̂) .

(9)

We now combine eqns. (7) and (8) to get
x̃T Ax̃ − x̂T Ax̂

≤

(x̃ − x̂)T Ax̃ + x̃T A (x̃ − x̂) + (x̂ − x̃)T A (x̃ − x̂)

= 2 x̃T A (x̃ − x̂) + (x̂ − x̃)T A (x̃ − x̂) .
Our next lemma will provide a bound for the first of the two quantities of interest in Lemma 3.
6

⋄

Lemma 4 Given our notation, with probability at least 7/8,
x̃Topt A (x̃opt − x̂opt ) ≤

p

8k/s.

Proof: Let D ∈ Rn×n be a diagonal matrix with entries Dii = δi for all i = 1 . . . n. Hence, we can write
x̂ = Dx̃. We have that


(1 − δ1 )A1∗ x̃
n
 (1 − δ2 )A2∗ x̃  X


(x̃ − x̂)T Ax̃ = x̃T (In − D) Ax̃ = x̃T 
(1 − δi ) x̃i Ai∗ x̃,
=
..


.
(1 − δn )An∗ x̃

i=1

where Ai∗ is the i-th row of the matrix A as a row vector. Squaring the above expression, we get
n X
n
2 X

(x̃ − x̂) Ax̃
x̃i ATi∗ x̃ (x̃j Aj∗ x̃) (1 − δi ) (1 − δj ) .
T

(10)

i=1 j=1

Recall that E (1 − δi ) = 0 for all i; thus, for all i 6= j, 1 − δi and 1 − δj are independent random variables
and therefore the expectation of their product is equal to zero. Thus, we can simplify the above expression
as follows:
n

2
X
E (1 − δi )2 (x̃i Ai∗ x̃)2
E (x̃ − x̂)T Ax̃
=
i=1
n 
X


1
=
− 1 x̃2i (Ai∗ x̃)2
pi
i=1

n 
X
1
− 1 x̃2i .
≤
pi
i=1

In the last inequality we used |Ai∗ x̃| ≤ kAi∗ k2 kx̃k2 ≤ 1. The last term in the above derivation can be
bounded as shown in eqn. (5), and thus we conclude

2
E (x̃ − x̂)T Ax̃ ≤ k/s.

Markov’s inequality now implies that with probability at least 7/8


2
(x̃ − x̂)T Ax̃ ≤ 8k/s.

⋄
Our next lemma will provide a bound for the second of the two quantities of interest in Lemma 3. The
proof of the lemma is tedious and a number of cases need to be considered.
Lemma 5 Given our notation, with probability at least 7/8,

1/2
√
(x̃opt − x̂opt )T A (x̃opt − x̂opt ) ≤ 24k2 /s2 + 6k2 /s3 + 54 k/s
.
7

Proof: Using the indicator variables δi and linearity of expectation, we get
2

n

2
X
E (x̃ − x̂)T A (x̃ − x̂) = E 
Aij (1 − δi ) x̃i (1 − δj ) x̃j  =
i,j=1

n
X

i1 ,j1 ,i2 ,j2 =1

Ai1 j1 Ai2 j2 x̃i1 x̃i2 x̃j1 x̃j2 × E ((1 − δi1 ) (1 − δi2 ) (1 − δj1 ) (1 − δj2 )) .

(11)

Recall that E (1 − δi ) = 0 for all i. Notice that if any of the four indices i1 , i2 , j1 , j2 appears only once, then
the expectation E ((1 − δi1 ) (1 − δi2 ) (1 − δj1 ) (1 − δj2 )) corresponding to those indices equals zero. This
expectation is non-zero if the four indices are paired in couples or if all four are equal. That is, non-zero
expectation happens if
(A) : i1 = i2 6= j1 = j2
(B) : i1 = j1 6= i2 = j2

(C) : i1 = j2 6= i2 = j1

(D) : i1 = i2 = j1 = j2

(n2 − n terms)

(n2 − n terms)

(n2 − n terms)

(n terms).

For case (A), let i1 = i2 = k and let j1 = j2 = ℓ, in which case the corresponding terms in eqn. (11)
become (notice that δk and δℓ are independent random variables since k 6= ℓ):
n


X
=
A2kℓ x̃2k x̃2ℓ E (1 − δk )2 (1 − δℓ )2
k,ℓ=1, k6=ℓ
n
X



1
1
−1
−1
≤
pk
pℓ
k,ℓ=1, k6=ℓ
X



n
n
X
1
1
2
2
x̃ℓ
−1
−1
=
x̃k
pk
pℓ
ℓ=1
k=1
 !2

n
X
k2
1
≤ 2.
−1
x̃2k
pk
s
A2kℓ x̃2k x̃2ℓ



(12)

k=1

In the first inequality we used |Akℓ | ≤ 1 for all k and ℓ and added extra positive terms (corresponding to
k = ℓ), which reinforce the inequality. The last inequality follows from eqn. (5).
For case (B), let i1 = j1 = k and let i2 = j2 = ℓ, in which case the corresponding terms in eqn. (11)
become (notice that δk and δℓ are independent random variables since k 6= ℓ):
n


X
Akk Aℓℓ x̃2k x̃2ℓ E (1 − δk )2 (1 − δℓ )2
=
k,ℓ=1, k6=ℓ
n
X



1
1
−1
−1
≤
pk
pℓ
k,ℓ=1, k6=ℓ
X



n
n
X
1
1
2
2
x̃ℓ
−1
−1
=
x̃k
pk
pℓ
ℓ=1
k=1
!2

n
X
k2
1
−1
≤ 2.
x̃2k
pk
s
Akk Aℓℓ x̃2k x̃2ℓ



k=1

8

(13)

In the first inequality we used Akk ≤ 1 for all k and the fact that the diagonal entries of a symmetric positive
definite matrix are non-negative, which allows us to add extra positive terms (corresponding to k = ℓ) to
reinforce the inequality. The remainder of the derivation is identical to case (A).
For case (C), let i1 = j2 = k and let i2 = j1 = ℓ, in which case the corresponding terms in eqn. (11)
become (notice that δk and δℓ are independent random variables since k 6= ℓ):
n
X

k,ℓ=1, k6=ℓ
n
X



Akℓ Aℓk x̃2k x̃2ℓ E (1 − δk )2 (1 − δℓ )2 =
A2kℓ x̃2k x̃2ℓ

k,ℓ=1, k6=ℓ





1
k2
1
−1
−1 ≤ 2.
pk
pℓ
s

In the first equality we used the fact that A is symmetric; the remainder of the derivation is identical to
case (A).
Finally, for case (D), let i1 = i2 = j1 = j2 = k, in which case the corresponding terms in eqn. (11)
become:
n
X
k=1

A2kk x̃4k E (1 − δk )4

=
≤
=
≤

n
X


4
1
6
− 2 + 3 −3
pk
pk pk
k=1


n
X
4
6
1
− 2 + 3 −3
x̃4k
pk
pk pk
k=1


X
1
4
6
4
+ 3 −3
−
x̃k
pk p2k
pk
k∈I <1


X
6
1
x̃4k
+ 3 .
pk pk
<1
A2kk x̃4k



k∈I

In the above derivation, we used |Akk | ≤ 1. We also split the set of indices {1 . . . n} in two subsets: the
set I =1 corresponding to indices k such that pk = 1 and the set I <1 corresponding to indices k such that
pk < 1. Note that for all k ∈ I <1 it must be the case that pk = s |x̃i | / kx̃k1 . Thus,
!
n
3
X
X
kx̃k
6
kx̃k
1
1
+
A2kk x̃4k E (1 − δk )4 ≤
x̃4k
3 |x̃ |3
s
|x̃
|
s
k
k
k=1
k∈I <1
!
n
3
X
kx̃k
6
kx̃k
1
1
+
≤
x̃4k
s |x̃k |
s3 |x̃k |3
k=1

=

n
n
kx̃k31 X
6 kx̃k1 X
3
|x̃k |
|x̃k | + 3
s
s
k=1

k=1

=

6 kx̃k1 kx̃k33
s

+

kx̃k41
.
s3

Recall that kx̃kq ≤ kx̃kp for any 1 ≤ p ≤ q ≤ ∞ and use kx̃k2 ≤ 1 and kx̃k1 ≤
9

√
k to get

n
X
k=1

A2kk x̃4k E (1 − δk )4 ≤
≤

6 kx̃k1 kx̃k32
kx̃k4
+ 31
s
s
√
6 k k2
+ 3.
s
s

(14)

Combining all four cases (i.e., eqns. (12), (13), (2.4), and (14)), we get

2
√
E (x̃ − x̂)T A (x̃ − x̂) ≤ 3k2 /s2 + k2 /s3 + 6 k/s.

Using Markov’s inequality and taking square roots concludes the proof of the lemma.
⋄
We can now complete the proof of the lower bound for x̂Topt Ax̂opt . First, combine Lemmas 3, 4, and 5
to get

1/2
√
p
x̃Topt Ax̃opt − x̂opt Ax̂opt ≤ 2 8k/s + 24k2 /s2 + 6k2 /s3 + 54 k/s
.
Since each of Lemmas 4 and 5 fail with probability at most 1/8, it follows from the union bound that
the above inequality holds with probability at least 1 − 2(1/8) = 3/4. Therefore, setting s = 200k/ǫ2
guarantees (after some algebra) that with probability at least 3/4,
x̃Topt Ax̃opt − x̂Topt Ax̂opt ≤ ǫ.
Using the triangle inequality and Lemma 2 we get that with probability at least 3/4
x̂Topt Ax̂opt ≥ x̃Topt Ax̃opt − ǫ ≥ xTopt Axopt − ǫ,
which concludes the proof of eqn. (3).

3 Experiments
We perform experiments on both real and synthetic datasets. We chose to compare our algorithm with
the solution returned by the state-of-the-art Spasm toolbox of [SCLE12], which implements the approach
proposed in [ZHT06]. We also compare our solution with the simple MaxComp heuristic [CJ95]: this
method computes the top singular vector of matrix A and returns a sparse vector by keeping the top k
largest components (in absolute value) and setting the remaining ones to zero.

3.1 Experimental setup
Let X ∈ Rm×n with n ≫ m denote the object-feature data matrix, where every column has zero mean, and
recall that A = XT X in eqn. (2). We use the following function to measure the quality of an approximate
solution x ∈ Rn to the sparse PCA problem:
f (x) = xT Ax/ kAk22 .

(15)

Notice that 0 ≤ f (x) ≤ 1 for all x which satisfy kxk2 ≤ 1. The closer f (x) is to one the more the vector x
is parallel to the top singular vector of A, or, equivalently, the closer f (x) is to one the more x captures the
variance of the matrix A that corresponds to its top singular value. Our goal is to calculate sparse vectors x
with f (x) ≈ 1.
10

ALGORITHM 2: Computing k sparse principal components
Input: X ∈ Rm×n , integer k.
Output: U = {u1 , . . . , uk }, V = {v1 , . . . , vk }.
Y = X⊤ v1 = rspca(X) and u1 = rspca(Y)
for i = 2, . . . , k do
⊤
X = X − Xvi−1 vi−1
and Y = Y − Yui−1 u⊤
i−1
vi = rspca(X) and ui = rspca(Y)
end

Our approach first finds a stationary point of the optimization problem of eqn. (2) and then uses Algorithm 1 (with s = k) to obtain a sparse solution vector x̂opt ∈ Rn . We note that the chosen value of s is
much smaller than the O k/ǫ2 value stipulated by Theorem 1. Also, in practice, our choice of s works
very well and results in solutions that are comparable or better than competing methods in our data.
We note that in our use of Spasm, we used soft thresholding by setting the STOP parameter to − kx̂opt k0
and δ = −∞ (both STOP and δ are parameters of Spasm toolbox). This implies that the solutions obtained
by Spasm and our approach will have the same number of non-zeros, thus making the comparison fair.
Similarly, for MaxComp, after computing the top singular vector of A, we select the kx̂opt k0 largest (in
absolute value) coordinates to form the sparse solution. A final technical note is that the solutions obtained
using either our method or MaxComp may result in vectors with non-unit Euclidean norms. In order to
achieve a fair comparison in terms of eqn. (15), there are two options. The first one (naive approach) is to
simply normalize the resulting vectors. However, a better approach (SVD-based) is possible: given a sparse
solution vector x̂opt , we could keep the rows and columns of A that correspond to the non-zero entries in
x̂opt and compute the top singular vector of the induced matrix. Note that the induced matrix would be a
kx̂opt k0 × kx̂opt k0 matrix and its top singular vector would be a kx̂opt k0 -dimensional vector. Obviously,
this latter vector would be a unit norm vector and it could be padded with zeros to derive a unit norm vector
in Rn with the same sparsity pattern as x̂opt . It is straight-forward to argue that this vector would improve
the value of f compared to the naive normalized vector x̂opt / kx̂opt k2 . In our experimental evaluation, we
will evaluate both the naive and the SVD-based normalization methods.
We also compare the methods based on their wall-clock running times. All experiments were run on a
Intel(R) Core(TM) i7-6700 machine running at 3.40GHz, with 16 GB of RAM.
In some of our experiments we will need to extract more than one sparse singular vectors. Towards that
end, we can use a deflation approach (Algorithm 2) to construct more than one sparse singular vectors by
first projecting the residual matrix into the space that is perpendicular to the top sparse singular vector and
then computing the top sparse singular vector of the residual matrix. In Algorithm 2, rspca refers to the
solution of the optimization problem of eqn. (2) followed by Algorithm 1).

3.2 Data Sets
We used 22 matrices emerging from population genetics, namely the 22 matrices (one for each chromosome) that encode all autosomal genotypes that are available in both the Human Genome Diversity
Panel [Con07] and the HAPMAP [LAT+ 08] project. Each of these matrices has approximately 2,500 samples (objects) described with respect to tens of thousands of features (Single Nucleotide Polymorphisms or
SNPs); see [PLJD10] for details. We also used a gene expression dataset (GSE10072, lung cancer gene
expression data) from the NCBI Gene Expression Omnibus database. This dataset includes 107 samples (58
lung tumor cases and 49 normal lung controls) measured using 22,215 probes (features) from the GPL96

11

platform annotation table. Both the population genetics and the gene expression datasets are interesting
in the context of sparse PCA beyond numerical evaluations, since the sparse components can be directly
interpreted to identify small sets of SNPs or genes that capture the data variance.
Our synthetic dataset has been carefully designed in order to highlight a setting where the MaxComp
heuristic will fail. More specifically, the absolute values of the entries of the largest singular vector of a
matrix in this family of matrices is not a good indicator of the importance of the respective entry in a sparse
PCA solution. Instead, the vector that emerges from the optimization problem of eqn. (2) is a much better
indicator. For a detailed description of the synthetic data generator see Section A.1 in the Appendix.
Our third data set comes from the field of text categorization and information retrieval. In such applications, documents are often represented as a “bag of words” and a vector space model is used. We use a
subset of the Classic-34 document collection, which we will call Classic-2. This subset consists of the CISI
(Comités Interministériels pour la Société de l’Information) collection (1,460 information science abstracts)
and the CRANFIELD collection (1,398 aeronautical systems abstracts). We created a document-by-term
matrix using the Text-to-Matrix Generator (TMG) [ZG06]5 ; the final matrix is a sparse 2, 858 × 12, 427 matrix with entries between zero and one, representing the weight of each term in the corresponding document.

3.3 Results
First we compare the performance of the different methods on a synthetic dataset, using the data generator
which was described in Section 3.2 with m = 27 , n = 212 , and θ ≈ 0.27π. In Figure 1a we plot (in
the y-axis) the value of the performance ratio f (x) (as defined in eqn. (15)) for our method (rspca), the
Spasm software, and the MaxComp heuristic. We also note that for each of the three methods, we use two
different approaches to normalize the resulting sparse vector: the naive and the SVD-based ones (see the
last paragraph in Section 3.1). As a result, we have a total of six possible methods to create and normalize a
sparse vector for sparse PCA. The x-axis shows the sparsity ratio of the resulting vector, namely kx̂opt k0 /n.
We remark that all six methods produce sparse vectors with exactly the same sparsity in order to have a fair
comparison. Notice that in Figure 1a, the MaxComp heuristic has worse performance when compared to
either our approach or to Spasm: this is expected, since we constructed this family of matrices in order to
precisely guarantee that the largest components of the top singular vector would not be good elements to
retain in the construction of a sparse vector. To further visualize this, we look at the sparse vectors, returned
by the different methods, in Figure 4 in the Appendix. In this figure, we present the resulting sparse vectors
for each of the three methods (normalization, obviously, is not relevant in Figure 4) for a particular choice
of the sparsity ratio ( kx̂opt k0 /n ≈ 0.1). Notice that MaxComp fails to capture the right sparsity pattern,
whereas our method and Spasm succeed.
In the second experiment, we present the performance of the different methods on the real datasets
described in Section 3.2. The results are shown in Figures 1b, 1c, and 1d. (We only show results for the first
two chromosomes of the joint HapMap and HGDP datasets; the other 20 chromosomes behave very similarly
and are shown in Figures 5, 6, 7, and 8 in the Appendix). We note that in the population genetics data,
our method has approximately the same or better performance compared to both MaxComp and Spasm.
Not surprisingly, the naive normalization approach is consistently worse than the SVD-based one. It is
worth noting that our SVD-based normalization approach easily improves the output of Spasm. This is
because Spasm does detect the correct sparsity pattern but fails to compute the appropriate coefficients of
the resulting sparse vectors.
4

ftp://ftp.cs.cornell.edu/pub/smart/
The matrix was created using the stoplist provided by TMG, the tf-idf weighting scheme, and document normalization. We
removed numbers and alphanumerics and we didn’t use stemming.
5

12

0.9

0.9

0.8

0.8

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

0.7

0.4

0.4
0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1
0.05

0.1

0.15

0.2

0.25

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.3

0.1

0.2

0.3

kx̂opt k0 /n

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

(a) Synthetic data: n = 212 .

(b) HapMap+HGDP data (chromosome 1): n = 37493.

0.9
0.9
0.8
0.8
0.7
0.7
0.6

f (x)

f (x)

0.6
0.5
0.4

0.5
0.4

0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.9

0.1

kx̂opt k0 /n

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

(c) HapMap+HGDP data (chromosome 2): n = 40844.

(d) Gene Expression: n = 22, 215.

Figure 1: f (x) vs. sparsity ratio kx̂opt k0 /n for various real and synthetic datasets.
In Figure 2 (and Figures 9, 10, 11 at the Appendix) we plot (in the y-axis) the running time for our
method (rspca), the Spasm software, and the MaxComp heuristic. We also note that for each of the
methods we use the naive approach (see the last paragraph in Section 3.1). The x-axis shows the sparsity
ratio of the resulting vector. Our algorithm presents a variable time performance that tends to decrease as the
sparsity ratio increases. This is an artifact of the projected gradient ascent algorithm [JNRS10] that we used
to find stationary points for the problem of eqn. (2). The smaller k is in (2), the harder the problem of (2)
becomes for projected gradient ascent. In practice we observed that the smaller k is the more iterations were
needed for projected gradient ascent to converge since the algorithm balances between the standard PCA
objective and the sparsity constraint. Projected gradient ascent, which solves problem (2) tends to be slower
for small k than the algorithm implemented in Spasm for solving the problem in [ZHT06]. We anticipate
that this issue can be solved by using more sophisticated optimization algorithms, we leave this for future
work. Notice that MaxComp is the fastest method, but based on Figure 1a there are cases where it fails to
find a good solution. The running time of MaxComp and Spasm appears to be constant. This is because
both methods require calculation of the first principal component which dominates the computational cost.
In particular, MaxComp is computing the first principal component which is then sparsified, while Spasm
is computing the first principal component to initialize the algorithm that solves the problem in [ZHT06].
All algorithms take into account the sparsity of the data.
13

22

25

rspca
Spasm
MaxComp

20
18

rspca
Spasm
MaxComp

20

time

time

16
14

15

12
10
10
8
6
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

kx̂opt k0 /n
(a) HapMap+HGDP data (chromosome 1): n = 37493.

(b) HapMap+HGDP data (chromosome 2): n = 40844.

22
12

rspca
Spasm
MaxComp

20
18

10

time

16

time

rspca
Spasm
MaxComp

11

14
12

9
8
7

10

6
8
5
6
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

kx̂opt k0 /n
(c) HapMap+HGDP data (chromosome 3): n = 34258.

(d) HapMap+HGDP data (chromosome 4): n = 30328.

Figure 2: Running time vs. sparsity ratio kx̂opt k0 /n for real data.
Finally, we evaluate the performance of our algorithm in a text mining application, using the Classic-2
document collection. First, we run Algorithm 2 to obtain two sparse singular vectors and we use the number
of their non-zero entries to compute two sparse singular vectors for MaxComp and Spasm6. This way we
can guarantee the same sparsity levels for all three pairs of singular vectors. We repeat this procedure for
8 different values of k (sparsity parameter in eqn. (1)). Table 1 summarizes, for each k, the variance and
the sparsity (in parenthesis) captured by the top two principal components using PCA, randomized sPCA
(rspca), MaxComp heuristic and Spasm or solving eqn. (2) (cvx).
6

Spasm does not support sparse input matrices.

14

Table 1: Variance and sparsity captured by the principal components. PCA results in dense principal components, while Spasm and MaxComp share the same sparsity with rspca.
k
Top Principal Comp.
Top two Principal Comp.
Top Principal Comp.
Top two Principal Comp.
Top Principal Comp.
Top two Principal Comp.
Top Principal Comp.
Top two Principal Comp.
Top Principal Comp.
Top two Principal Comp.
Top Principal Comp.
Top two Principal Comp.
Top Principal Comp.
Top two Principal Comp.
Top Principal Comp.
Top two Principal Comp.

100
500
1000
1500
2000
4000
8000
10500

pca
0.4351
0.6802
0.4351
0.6802
0.4351
0.6802
0.4351
0.6802
0.4351
0.6802
0.4351
0.6802
0.4351
0.6802
0.4351
0.6802

cvx
0.3077 (99%)
0.4897 (99%)
0.3880 (95%)
0.6073 (95%)
0.4136 (90%)
0.6486 (90%)
0.4242 (84%)
0.6649 (82%)
0.4295 (75%)
0.6730 (70%)
0.4350 (6%)
0.6801 (3%)
0.4351 (0%)
0.6802 (0%)
0.4351 (0%)
0.6802 (0%)

rspca
0.2942 (99%)
0.4680 (99%)
0.3728 (98%)
0.5864 (98%)
0.4005 (95%)
0.6294 (95%)
0.4120 (93%)
0.6470 (93%)
0.4190 (91%)
0.6572 (91%)
0.4278 (81%)
0.6700 (81%)
0.4324 (68%)
0.6764 (69%)
0.4332 (63%)
0.6776 (64%)

MaxComp
0.1955
0.3391
0.3353
0.5399
0.3825
0.6074
0.4013
0.6342
0.4133
0.6503
0.4284
0.6710
0.4326
0.6768
0.4333
0.6778

Spasm
0.2768
0.4227
0.3601
0.5701
0.3912
0.6163
0.4039
0.6361
0.4131
0.6491
0.4271
0.6690
0.4316
0.6752
0.4324
0.6764

It seems that Spasm and MaxComp capture less variance than the randomized sPCA. Furthermore, the
variance captured by randomized sPCA is constantly close to the one captured by the solution of eqn. (2),
but with a sparser component as Table 1 indicates.
Table 2 summarizes the terms with non-zero weights in randomized sPCA principal components with
sparsity parameter k = 100. The terms are ranked in descending order with respect to their weights. Notice
that the first principal component reveals terms that appear more often in the CRANFIELD collection while
the second principal component reveals terms that appear mostly in the CISI collection. CRANFIELD’s
terms are more singular than these of CISI’s and they tend to dominate the singular vectors since they tend
to appear more in the documents associated with the CRANFIELD collection than in the entire CLASSIC2 collection (e.g., the word “boundary” has one appearance in CISI and 459 in CRANFIELD). The exact
opposite happens for terms in CISI: a significant amount of these terms appear with high weights in both collections (e.g., the word “information” has 664 appearances in CISI and 44 in CRANFIELD). This indicates
that these terms will appear in singular vectors that do not separate the two collections.

15

Table 2: Non-zero terms of the randomized sPCA when k = 100
1st Principal Component
boundary
layer
heat
flow
transfer
laminar
plate
shock
turbulent
hypersonic
pressure
temperature
wall
flat
surface
mach
friction
compressible
velocity
reynolds
stagnation
skin
number
stream
equations
transition
solution
viscous
layers
separation
solutions
gradient
injection
dimensional
incompressible
fluid
edge
cylinder
shear
point
interaction
numbers
body
wave
method
leading
region
bodies
supersonic
gas
measurements
approximate
local
equation
case
constant

2nd Principal Component
information
library
retrieval
systems
system
scientific
science
libraries
research
layer
services
boundary
indexing
search
literature
heat
university
documents
users
classification
technology

16

4 Open problems
From a theoretical perspective, it would be interesting to explore whether other relaxations of the sparse
PCA problem of eqn. (1), combined with randomized rounding procedures, could improve our error bounds
in Theorem 1. It would also be interesting to formally analyze the deflation algorithm that computes more
than one sparse singular vectors in a randomized manner (Algorithm 2). Finally, from a complexity theory
perspective, we are not aware of any inapproximability results for the sparse PCA problem; to the best of
our knowledge, it is not known whether relative error approximations are possible (without any assumptions
on the input matrix).

17

APPENDIX
A.1 Synthetic data generator
In order to generate our synthetic dataset, we generated the following matrix:
X = UΣV⊺ + Eσ ,
where U ∈ Rm×m and V ∈ Rn×n are orthonormal. The matrix Σ ∈ Rm×n has m distinct singular values
σi in its diagonal and the matrix Eσ ∈ Rm×n is a noise matrix parameterized by σ > 0.
We set U to be a Hadamard matrix with normalized columns; we set Σ to have entries σ1 = 100 and
σi = 1/ei for all i = 2, . . . , m. The entries of the matrix Eσ follow a zero-mean normal distribution with
standard deviation σ = 10−3 . We now describe how to construct the matrix V: we set V = Gd (θ)Ṽ, where
Ṽ is also a Hadamard matrix with normalized columns. Here Gn (θ) is a composition of Givens rotations.
In particular, Gn (θ) is a composition of n/4 Givens rotations with the same angle θ for every rotation. More
precisely, let G(i, j, θ) ∈ Rn×n be a Givens rotation matrix, which rotates the plane i-j by an angle θ:


1 ··· 0 ··· 0 ··· 0
 .. . .
.. 
..
..
.
.
.
.
.


0 · · · c1 · · · −c2 · · · 0



.. . .
..
..  ,
G(i, j, θ) =  ...
.
.
.
.


0 · · · c2 · · · c1 · · · 0


 ..
..
..
. . .. 
.
. .
.
.
0 ··· 0 ··· 0 ··· 1
where i, j ∈ {1, . . . , n}, c1 = cos θ and c2 = sin θ. We define the composition as follows:

Gn (θ) = G(i1 , j1 , θ)G(i2 , j2 , θ) · · · · · · , G(ik , jk , θ), · · · , G(in/4 , jn/4 , θ)
with
ik =

n
+ 2k − 1,
2

jk =

n
+ 2k
2

n
for k = 1, . . . , .
4

The matrix Gn (θ) rotates (in a pairwise manner) the bottom n/2 components of the columns of Ṽ. Since
the Hadamard matrix has entries equal to +1 or -1 (up to normalization), we will pick a value of θ that
guarantees that, after rotation, n/4 components of the columns of Ṽ will be almost zero. Thus, the resulting
matrices will have about a quarter of components set at a large value, a quarter of their components set
at roughly zero, and the rest set at a moderate value. For example, let n = 212 and θ ≈ 0.27π; then,
the difference between the first column of matrix V and Ṽ is presented in Figure 3, where we plotted the
(sorted) absolute values of the components of the first column of the matrices V and Ṽ.

18

magnitude

0.02

0.015

0.01

0.005

without rotation
with rotation
0

500

1000

1500

2000

2500

3000

3500

4000

index

Figure 3: The red dashed line corresponds to the sorted absolute values of the components of the first
column of matrix V. Similarly, the blue line corresponds to the first column of Ṽ.

(a) Actual eigenvector

(b) rspca

(c) Spasm

(d) MaxComp

Figure 4: Sparse PCA solutions for a synthetic experiment.

19

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

0.9

0.4

0.4
0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.9

0.1

0.2

0.3

0.5

0.6

0.7

0.8

0.9

(b) HapMap+HGDP data (chromosome 4): n = 30328.

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

(a) HapMap+HGDP data (chromosome 3): n = 34258.

0.4

0.4
0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.9

0.1

0.2

0.3

kx̂opt k0 /n

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

(c) HapMap+HGDP data (chromosome 5): n = 31479.

(d) HapMap+HGDP data (chromosome 6): n = 32800.

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

0.4

kx̂opt k0 /n

kx̂opt k0 /n

0.4

0.4
0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.9

0.1

kx̂opt k0 /n

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

(e) HapMap+HGDP data (chromosome 7): n = 27130.

(f) HapMap+HGDP data (chromosome 8): n = 28508.

Figure 5: Performance of sparse PCA algorithms on additional population genetics data (chromosomes 3-8).

20

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

0.9

0.4

0.4

0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.1

0.9

0.1

kx̂opt k0 /n

0.2

0.3

0.5

0.6

0.7

0.8

0.9

(b) HapMap+HGDP data (chromosome 10): n = 26327.
0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

0.9

0.4

0.4
0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.9

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

kx̂opt k0 /n

(c) HapMap+HGDP data (chromosome 11): n = 24440.

(d) HapMap+HGDP data (chromosome 12): n = 24084.

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

0.4

kx̂opt k0 /n

(a) HapMap+HGDP data (chromosome 9): n = 24070.

f (x)

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2

0.4

0.4
0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.9

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

kx̂opt k0 /n

(e) HapMap+HGDP data (chromosome 13): n = 18958.

(f) HapMap+HGDP data (chromosome 14): n = 16469.

Figure 6: Performance of sparse PCA algorithms on additional population genetics data (chromosomes
9-14).

21

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

0.9

0.4

0.4

0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.9

0.1

0.2

0.3

kx̂opt k0 /n

0.6

0.7

0.8

0.9

(b) HapMap+HGDP data (chromosome 16): n = 15289.

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

0.5

kx̂opt k0 /n

(a) HapMap+HGDP data (chromosome 15): n = 15351.

0.4

0.4
0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.9

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

kx̂opt k0 /n

(c) HapMap+HGDP data (chromosome 17): n = 12945.

(d) HapMap+HGDP data (chromosome 18): n = 15373.

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

0.4

0.4

0.4

0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.9

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

kx̂opt k0 /n

(e) HapMap+HGDP data (chromosome 19): n = 8465.

(f) HapMap+HGDP data (chromosome 20): n = 13015.

Figure 7: Performance of sparse PCA algorithms on additional population genetics data (chromosomes
15-20).

22

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

f (x)

f (x)

0.9

0.4

0.4
0.3

0.3

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

rspca (SVD-based)
Spasm (SVD-based)
MaxComp (SVD-based)
rspca (naive)
Spasm (naive)
MaxComp (naive)

0.2
0.1

0.1

0.9

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

kx̂opt k0 /n

(a) HapMap+HGDP data (chromosome 21): n = 7556.

(b) HapMap+HGDP data (chromosome 22): n = 7334.

Figure 8: Performance of sparse PCA algorithms on additional population genetics data (chromosomes
21-22).

23

16

16

rspca
Spasm
MaxComp

14
12

12

time

time

rspca
Spasm
MaxComp

14

10

10

8

8

6

6

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.3

(a) HapMap+HGDP data (chromosome 5): n =
31479.

10

0.5

0.6

0.7

0.8

0.9

(b) HapMap+HGDP data (chromosome 6): n =
32800.

rspca
Spasm
MaxComp

11

0.4

kx̂opt k0 /n

kx̂opt k0 /n

rspca
Spasm
MaxComp

16
14

time

time

9
8
7

12
10

6

8

5

6

4
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

kx̂opt k0 /n

(c) HapMap+HGDP data (chromosome 7): n =
27130.

(d) HapMap+HGDP data (chromosome 8): n =
28508.

20

rspca
Spasm
MaxComp

14

rspca
Spasm
MaxComp
15

time

time

12
10

10

8
6

5

4
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

kx̂opt k0 /n

(e) HapMap+HGDP data (chromosome 9): n =
24070.

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

(f) HapMap+HGDP data (chromosome 10): n =
26327.

Figure 9: Running time of sparse PCA algorithms on additional population genetics data (chromosomes
5-10).

24

12
12

rspca
Spasm
MaxComp

11
10

10
9

time

9

time

rspca
Spasm
MaxComp

11

8

8
7

7
6

6

5

5
4

4
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.3

kx̂opt k0 /n

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

(a) HapMap+HGDP data (chromosome 11): n =
24440.

(b) HapMap+HGDP data (chromosome 12): n =
24084.

12

rspca
Spasm
MaxComp

10

rspca
Spasm
MaxComp

9
8

time

time

7
8

6

6
5
4

4
3
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.3

kx̂opt k0 /n

0.6

0.7

0.8

0.9

(d) HapMap+HGDP data (chromosome 14): n =
16469.

rspca
Spasm
MaxComp

7

0.5

kx̂opt k0 /n

(c) HapMap+HGDP data (chromosome 13): n =
18958.

8

0.4

rspca
Spasm
MaxComp

6
5.5

time

time

5
6
5

4.5
4
3.5

4

3

3

2.5
2
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

kx̂opt k0 /n

(e) HapMap+HGDP data (chromosome 15): n =
15351.

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

(f) HapMap+HGDP data (chromosome 16): n =
15289.

Figure 10: Running time of sparse PCA algorithms on additional population genetics data (chromosomes
11-16).

25

5

9

rspca
Spasm
MaxComp

4.5

7

time

4

time

rspca
Spasm
MaxComp

8

3.5

6
5

3

4
2.5
3
2
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.3

kx̂opt k0 /n

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

(a) HapMap+HGDP data (chromosome 17): n =
12945.

(b) HapMap+HGDP data (chromosome 18): n =
15373.

4.5

7

rspca
Spasm
MaxComp

4

0.4

rspca
Spasm
MaxComp

6

time

time

3.5
3

5

4

2.5
2

3

1.5

2
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.3

kx̂opt k0 /n

0.6

0.7

0.8

0.9

(d) HapMap+HGDP data (chromosome 20): n =
13015.

4.5

rspca
Spasm
MaxComp

3

0.5

kx̂opt k0 /n

(c) HapMap+HGDP data (chromosome 19): n =
8465.

3.5

0.4

rspca
Spasm
MaxComp

4

time

time

3.5
2.5

3
2.5

2

2
1.5
1.5
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

kx̂opt k0 /n

(e) HapMap+HGDP data (chromosome 21): n =
7556.

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

kx̂opt k0 /n

(f) HapMap+HGDP data (chromosome 22): n =
7334.

Figure 11: Running time of sparse PCA algorithms on additional population genetics data (chromosomes
17-22).

26

References
[AW08]

A. A. Amini and M. J. Wainwright. High-dimensional Analysis of Semidefinite Relaxations
for Sparse Principal Components. In Information Theory, 2008. ISIT 2008. IEEE International
Symposium on,, pages 2454–2458, 2008.

[CJ95]

J. Cadima and I. T. Jolliffe. Loading and correlations in the interpretation of principle compenents. Journal of Applied Statistics, pages 203–214, 1995.

[Con07]

The International HapMap Consortium. A second generation human haplotype map of over
3.1 million SNPs. Nature, 449:851–861, 2007.

[dBG08]

A. d’Aspremont, F. Bach, and L. El Ghaoui. Optimal Solutions for Sparse Principal Component
Analysis. The Journal of Machine Learning Research, pages 1269–1294, 2008.

[dBG14]

A. d’Aspremont, F. Bach, and L. El Ghaoui. Approximation Bounds for Sparse Principal
Component Analysis. Math. Program., Ser. B, pages 89–110, 2014.

[dGJ07]

A. d’Aspremont, L. El Ghaoui, and M. I. Jordan. A Direct Formulation for Sparse PCA using
Semidefinite Programming. SIAM Review, pages 434–448, 2007.

[DXXM07] T. Dacheng, L. Xuelong, W. Xindong, and S. Maybank. General tensor discriminant analysis
and gabor features for gait recognition. IEEE Trans. Pattern Anal. Mach. Intell., 29(10):1700–
1715, 2007.
[JNRS10]

Michel Journée, Y. Nesterov, P. Richtárik, and R. Sepulchre. Generalized Power Method for
Sparse Principal Component Analysis. Journal of Machine Learning Research, pages 517–553,
2010.

[Jol95]

I. T. Jolliffe. Rotation of principal components: choice of normalization constraints. Journal
of Applied Statistics, pages 29–35, 1995.

[JTU03]

I. T. Jolliffe, N. T. Trendafilov, and M. Uddin. A modified principal component technique
based on the lasso. Journal of Computational and Graphical Statistics, pages 531–547, 2003.

[LAT+ 08]

J. Z. Li, D. M. Absher, H. Tang, A. M. Southwick, A. M. Casto, S. Ramachandran, H. M.
Cann, G. S. Barsh, M. Feldman, L. L. Cavalli-Sforza, and R. M. Myers. Worldwide human
relationships inferred from genome-wide patterns of variation. Science, 319:1100–1104, 2008.

[LWT+ 16] Y. Luo, Y. Wen, D. Tao, J. Gui, and C. Xu. Large margin multi-modal multi-task feature
extraction for image classification. IEEE Trans. on Image Processing, 25(1):414–427, 2016.
[LZW+ 16] W. Liu, Z. Zha, Y. Wang, K. Lu, and D. Tao. p-laplacian regularized sparse coding for human
activity recognition. IEEE Trans. on Industrial Electronics, 63(8):5120–5129, 2016.
[MD09]

M.W. Mahoney and P. Drineas. CUR matrix decompositions for improved data analysis. In
Proceedings of the National Academy of Sciences, pages 697–702, 106 (3), 2009.

[MWA12]

B. Moghaddam, Y. Weiss, and S. Avidan. Generalized Spectral Bounds for Sparse lda. In
Proceedings of the 23rd international conference on Machine learning, pages 641–648, 2012.

27

[PDK13]

D. S. Papailiopoulos, A. G. Dimakis, and S. Korokythakis. Sparse PCA through Low-rank
Approximations. http://arxiv.org/abs/1303.0551, 2013.

[PLJD10]

P. Paschou, J. Lewis, A. Javed, and P. Drineas. Ancestry informative markers for fine-scale
individual assignment to worldwide populations. J. Med. Genet., 47(12):835–847, Dec 2010.

[SCLE12]

K. Sjöstrand, L.H. Clemmensen, R. Larsen, and B. Ersbøll. Spasm: A matlab toolbox for
sparse statistical modeling. In Journal of Statistical Software (Accepted for publication), 2012.

[YZ13]

X.-T. Yuan and T. Zhang. Truncated Power Method for Sparse Eigenvalue Problems. Journal
of Machine Learning Research, pages 899–925, 2013.

[ZG06]

D. Zeimpekis and E. Gallopoulos. TMG: A MATLAB toolbox for generating term document
matrices from text collections. In Grouping Multidimensional Data: Recent Advances in Clustering, pages 187–210. J. Kogan and C. Nicholas and M. Teboulle, 2006.

[ZHT06]

H. Zou, T. Hastie, and R Tibshirani. Sparse Principal Component Analysis. Journal of Computational and Graphical Statistics, pages 265–286, 2006.

28

