Block Coordinate Descent for Sparse NMF
Authors: Vamsi Potluru, Sergey M. Plis, Jonathan Le Roux, Barak A. Pearlmutter, Vince D. Calhoun, Thomas P. Hayes
ICLR 2013

Abstract: Nonnegative matrix factorization (NMF) has become a ubiquitous tool for data analysis. An important variant is the sparse NMF problem which arises when we explicitly require the learnt features to be sparse. A natural measure of sparsity is the L$_0$ norm, however its optimization is NP-hard. Mixed norms, such as L$_1$/L$_2$ measure, have been shown to model sparsity robustly, based on intuitive attributes that such measures need to satisfy. This is in contrast to computationally cheaper alternatives such as the plain L$_1$ norm. However, present algorithms designed for optimizing the mixed norm L$_1$/L$_2$ are slow and other formulations for sparse NMF have been proposed such as those based on L$_1$ and L$_0$ norms. Our proposed algorithm allows us to solve the mixed norm sparsity constraints while not sacrificing computation time. We present experimental evidence on real-world datasets that shows our new algorithm performs an order of magnitude faster compared to the current state-of-the-art solvers optimizing the mixed norm and is suitable for large-scale datasets.

Block Coordinate Descent for Sparse NMF
Vamsi K. Potluru
Department of Computer Science,
University of New Mexico
ismav@cs.unm.edu
Sergey M. Plis
Mind Research Network,
splis@mrn.org
Jonathan Le Roux
Mitsubishi Electric Research Labs
leroux@merl.com
Barak A. Pearlmutter
Department of Computer Science,
National University of Ireland Maynooth
barak@cs.nuim.ie
Vince D. Calhoun
Electrical and Computer Engineering, UNM and
Mind Research Network
vcalhoun@mrn.org
Thomas P. Hayes
Department of Computer Science,
University of New Mexico
hayes@cs.unm.edu
Abstract
Nonnegative matrix factorization (NMF) has become a ubiquitous tool for
data analysis. An important variant is the sparse NMF problem which
arises when we explicitly require the learnt features to be sparse. A natural
measure of sparsity is the L 0 norm, however its optimization is NP-hard.
Mixed norms, such as L 1/L2 measure, have been shown to model sparsity
robustly, based on intuitive attributes that such measures need to satisfy.
This is in contrast to computationally cheaper alternatives such as the plain
L1 norm. However, present algorithms designed for optimizing the mixed
norm L1/L2 are slow and other formulations for sparse NMF have been pro-
posed such as those based on L 1 and L0 norms. Our proposed algorithm
allows us to solve the mixed norm sparsity constraints while not sacriﬁcing
computation time. We present experimental evidence on real-world datasets
that shows our new algorithm performs an order of magnitude faster com-
pared to the current state-of-the-art solvers optimizing the mixed norm and
is suitable for large-scale datasets.
1 Introduction
Matrix factorization arises in a wide range of application domains and is useful for extracting
the latent features in the dataset (Figure 1). In particular, we are interested in matrix
factorizations which impose the following requirements:
•nonnegativity
•low-rankedness
•sparsity
Nonnegativity is a natural constraint when modeling data with physical constraints such as
chemical concentrations in solutions, pixel intensities in images and radiation dosages for
cancer treatment. Low-rankedness is useful for learning a lower dimensionality represen-
tation. Sparsity is useful for modeling the conciseness of the representation or that of the
1
arXiv:1301.3527v2  [cs.LG]  18 Mar 2013
latent features. Imposing all these requirements on our matrix factorization leads to the
sparse nonnegative matrix factorization (SNMF) problem.
SNMF enjoys quite a few formulations [2, 14, 13, 11, 24, 17, 25, 26] with successful applica-
tions to single-channel speech separation [27] and micro-array data analysis [17, 25].
However, algorithms [14, 11] for solving SNMF which utilize the mixed norm of L 1/L2 as
their sparsity measure are slow and do not scale well to large datasets. Thus, we develop
an eﬃcient algorithm to solve this problem and has the following ingredients:
• A theoretically eﬃcient projection operator (O(mlog m)) to enforce the user-deﬁned
sparsity where m is the dimensionality of the feature vector as opposed to the
previous approach [14].
• Novel sequential updates which provide the bulk of our speedup compared to the
previously employed batch methods [14, 11].
Figure 1: (Left) Features learned from the ORL dataset 2with various matrix factorization
methods such as principal component analysis (PCA), independent component analysis
(ICA), and dictionary learning. The relative merit of the various matrix factorizations
depends on both the signal domain and the target application of interest. (Right) Features
learned under the sparse NMF formulation where roughly half the features were constrained
to lie in the interval [0 .2,0.4] and the rest are ﬁxed to sparsity value 0 .7. This illustrates
the ﬂexibility that the user has in ﬁne tuning the feature sparsity based on prior domain
knowledge. White pixels in this ﬁgure correspond to the zeros in the features.
2 Preliminaries and Previous Work
In this section, we give an introduction to the nonnegative matrix factorization (NMF) and
SNMF problems. Also, we discuss some widely used algorithms from the literature to solve
them.
Both these problems share the following problem and solution structure. At a high-level,
given a nonnegative matrix X of size m×n, we want to approximate it with a product of
two nonnegative matrices W,H of sizes m×r and r×n, respectively:
X ≈WH.(1)
The nonnegative constraint on matrix H makes the representation a conical combination
of features given by the columns of matrix W. In particular, NMF can result in sparse
representations, or a parts-based representation, unlike other factorization techniques such
as principal component analysis (PCA) and vector quantization (VQ). A common theme in
2Scikit-learn package was used in generating the ﬁgure.
2
the algorithms proposed for solving these problems is the use of alternating updates to the
matrix factors, which is natural because the objective function to be minimized is convex
in W and in H, separately, but not in both together. Much eﬀort has been focused on
optimizing the eﬃciency of the core step of updating one of W,H while the other stays
ﬁxed.
2.1 Nonnegative Matrix Factorization
Factoring a matrix, all of whose entries are nonnegative, as a product of two low-rank non-
negative factors is a fundamental algorithmic challenge. This has arisen naturally in diverse
areas such as image analysis [20], micro-array data analysis [17], document clustering [31],
chemometrics [19], information retrieval [12] and biology applications [4]. For further appli-
cations, see the references in the following papers [1, 7].
We will consider the following version of the NMF problem, which measures the reconstruc-
tion error using the Frobenius norm [21]:
min
W,H
1
2∥X −WH∥2
F s.t. W ≥0, H ≥0, ∥Wj∥2 = 1, ∀j ∈{1,··· ,r}(2)
where ≥is element-wise. We use subscripts to denote column elements. Simple multiplica-
tive updates were proposed by Lee and Seung to solve the NMF problem. This is attractive
for the following reasons:
• Unlike additive gradient descent methods, there is no arbitrary learning rate pa-
rameter that needs to be set.
• The nonnegativity constraint is satisﬁed automatically, without any additional pro-
jection step.
• The objective function converges to a limit point and the values are non-increasing
across the updates, as shown by Lee and Seung [21].
Algorithm 1 is an example of the kind of multiplicative update procedure used, for instance,
by Lee and Seung [21]. The algorithm alternates between updating the matrices W and H
(we have only shown the updates for H—those for W are analogous).
Algorithm 1nnls-mult(X,W,H)
1: repeat
2: H = H ⊙ W⊤X
W⊤WH.
3: until convergence
4: Output: Matrix H.
Here, ⊙indicates element-wise (Hadamard) product and matrix division is also element-
wise. To remove the scaling ambiguity, the norm of columns of matrix W are set to unity.
Also, a small constant, say 10 −9, is added to the denominator in the updates to avoid
division by zero.
Besides multiplicative updates, other algorithms have been proposed to solve the NMF
problem based on projected gradient [22], block pivoting [18], sequential constrained opti-
mization [6] and greedy coordinate-descent [15].
2.2 Sparse Nonnegative Matrix Factorization
The nonnegative decomposition is in general not unique [9]. Furthermore, the features may
not be parts-based if the data resides well inside the positive orthant. To address these
issues, sparseness constraints have been imposed on the NMF problem.
Sparse NMF can be formulated in many diﬀerent ways. From a user point of view, we can
split them into two classes of formulations: explicit and implicit. In explicit versions of
SNMF [14, 11], one can set the sparsities of the matrix factors W,Hdirectly. On the other
3
hand, in implicit versions of SNMF [17, 25], the sparsity is controlled via a regularization
parameter and is often hard to tune to speciﬁed sparsity values a priori. However, the
algorithms for implicit versions tend to be faster compared to the explicit versions of SNMF.
In this paper, we consider the explicit sparse NMF formulation proposed by Hoyer [14].
To make the presentation easier to follow, we ﬁrst consider the case where the sparsity is
imposed on one of the matrix factors, namely the feature matrix W—the analysis for the
symmetric case where the sparsity is instead set on the other matrix factor H is analogous.
The case where sparsity requirements are imposed on both the matrix factors is dealt with in
the Appendix. The sparse NMF problem formulated by Hoyer [14] with sparsity on matrix
W is as follows:
min
W,H
f(W,H) =1
2∥X −WH∥2
F s.t. W ≥0,H ≥0,
∥Wj∥2 = 1, sp(Wj) = α, ∀j ∈{1,··· ,r}(3)
Sparsity measure for a d-dimensional vector xis given by:
sp(x) =
√
d−∥x∥1/∥x∥2√
d−1
(4)
The sparsity measure (4) deﬁned above has many appealing qualities. Some of which are as
follows:
• The measure closely models the intuitive notion of sparsity as captured by the L0
norm. So, it easy for the user to specify sparsity constraints from prior knowledge
of the application domain.
• Simultaneously, it is able to avoid the pitfalls associated with directly optimizing
the L0 norm. Desirable properties for sparsity measures have been previously ex-
plored [16] and it satisﬁes all of these properties for our problem formulation. The
properties can be brieﬂy summarized as: (a) Robin Hood — Spreading the energy
from larger coordinates to smaller ones decreases sparsity, (b) Scaling — Sparsity
is invariant to scaling, (c) Rising tide — Adding a constant to the coordinates
decreases sparsity, (d) Cloning — Sparsity is invariant to cloning, (e) Bill Gates
— One big coordinate can increase sparsity, (f) Babies — coordinates with zeros
increase sparsity.
• The above sparsity measure enables one to limit the sparsity for each feature to lie
in a given range by changing the equality constraints in the SNMF formulation (3)
to inequality constraints [11]. This could be useful in scenarios like fMRI brain
analysis, where one would like to model the prior knowledge such as sizes of artifacts
are diﬀerent from that of the brain signals. A sample illustration on a face dataset
is shown in Figure 1 (Right). The features are now evenly split into two groups of
local and global features by choosing two diﬀerent intervals of sparsity.
A gradient descent-based algorithm called Nonnegative Matrix Factorization with Sparse-
ness Constraints (NMFSC) to solve SNMF was proposed [14]. Multiplicative updates
were used for optimizing the matrix factor which did not have sparsity constraints spec-
iﬁed. Heiler and Schn¨ orr[11] proposed two new algorithms which also solved this prob-
lem by sequential cone programming and utilized general purpose solvers like MOSEK
(http://www.mosek.com). We will consider the faster one of these called tangent-plane
constraint (TPC) algorithm. However, both these algorithms, namely NMFSC and TPC,
solve for the whole matrix of coeﬃcients at once. In contrast, we propose a block coordinate-
descent strategy which considers a sequence of vector problems where each one can be solved
in closed form eﬃciently.
3 The Sequential Sparse NMF Algorithm
We present our algorithm which we call Sequential Sparse NMF (SSNMF) to solve the
SNMF problem as follows:
4
First, we consider a problem of special form which is the building block (Algorithm 2) of our
SSNMF algorithm and give an eﬃcient, as well as exact, algorithm to solve it. Second, we
describe our sequential approach (Algorithm 3) to solve the subproblem of SNMF. This uses
the routine we developed in the previous step. Finally, we combine our routines developed
in the previous two steps along with standard solvers (for instance Algorithm 1) to complete
the SSNMF Algorithm (Algorithm 4).
3.1 Sparse-opt
Sparse-opt routine solves the following subproblem which arises when solving problem (3):
maxy≥0
b⊤ys.t. ∥y∥1 = k,∥y∥2 = 1 (5)
where vector b is of size m. This problem has been previously considered [14], and an
algorithm to solve it was proposed which we will henceforth refer to as the Projection-
Hoyer. Similar projection problems have been recently considered in the literature and
solved eﬃciently [10, 5].
Observation 1. For anyi,j, we have that ifbi ≥bj, then yi ≥yj.
Let us ﬁrst consider the case when the vector bis sorted. Then by the previous observation,
we have a transition point p that separates the zeros of the solution vector from the rest.
Observation 2. By applying the Cauchy-Schwarz inequality ony and the all ones vector,
we get p≥k2.
The Lagrangian of the problem (5) is :
L(y,µ,λ, γ) = b⊤y+ µ
(m∑
i=1
yi −k
)
+ λ
2
(m∑
i=1
y2
i −1
)
+ γ⊤y
Setting the partial derivatives of the Lagrangian to zero, we get by observation 1:
m∑
i
yi = k,
m∑
i
y2
i = 1
bi + µ(p) + λ(p)yi = 0,∀i∈{1,2,··· ,p}
γi = 0,∀i∈{1,··· ,p}
yi = 0,∀i∈{p+ 1,··· ,m}
where we account for the dependence of the Lagrange parameters λ,µ on the transition
point p. We compute the objective value of problem (5) for all transition points p in
the range from k2 to m and select the one with the highest value. In the case, where
the vector b is not sorted, we just simply sort it and note down the sorting permutation
vector. The complete algorithm is given in Algorithm 2. The dominant contribution to the
running time of Algorithm 2 is the sorting of vector b and therefore can be implemented
in O(mlog m) time3. Contrast this with the running time of Projection-Hoyer whose worst
case is O(m2) [14, 28].
3.2 Sequential Approach —Block Coordinate Descent
Previous approaches for solving SNMF [14, 11] use batch methods to solve for sparsity
constraints. That is, the whole matrix is updated at once and projected to satisfy the con-
straints. We take a diﬀerent approach of updating a column vector at a time. This gives us
the beneﬁt of being able to solve the subproblem (column) eﬃciently and exactly. Subse-
quent updates can beneﬁt from the newly updated columns resulting in faster convergence
as seen in the experiments.
3This can be further reduced to linear time by noting that we do not need to fully sort the input
in order to ﬁnd p∗.
5
Algorithm 2Sparse-opt(b,k)
1: Set a= sort(b) and p∗= m. Get a mapping π such that ai = bπ(i) and aj ≥aj+1 for
all valid i,j.
2: Compute values of µ(p),λ(p) as follows:
3: for p∈{⌈k2⌉,m}do
4: λ(p) = −
√
p∑p
i=1 a2
i −(
∑p
i=1 ai)
2
(p−k2)
5: µ(p) = −
∑p
i=1 ai
p −k
pλ(p)
6: if a(p) <−µ(p) then
7: p∗= p−1
8: break
9: end if
10: end for
11: Set xi = −ai+µ(p∗)
λ(p∗) ,∀i∈{1,··· ,p∗}and to zero otherwise.
12: Output: Solution vector ywhere yπ(i) = xi.
In particular, consider the optimization problem (3) for a column j of the matrix W while
ﬁxing the rest of the elements of matrices W,H:
min
Wj≥0
˜f(Wj) = 1
2g∥Wj∥2
2 + u⊤Wj s.t. ∥Wj∥2 = 1, ∥Wj∥1 = k
where g= H⊤
j Hj and u= −XH⊤
j +∑
i̸=jWi(HH⊤)ij. This reduces to the problem (5) for
which we have proposed an exact algorithm (Algorithm 2). We update the columns of the
matrix factor W sequentially as shown in Algorithm 3. We call it sequential for we update
the columns one at a time. Note that this approach can be seen as an instance of block
coordinate descent methods by mapping features to blocks and the Sparse-opt projection
operator to a descent step.
Algorithm 3sequential-pass(X,W,H)
1: C = −XH⊤+ WHH⊤
2: G = HH⊤
3: repeat
4: for j = 1 to r (randomly) do
5: Uj = Cj −WjGjj
6: t= Sparse-opt(−Uj,k).
7: C = C + (t−Wj) G⊤
j
8: Wj = t.
9: end for
10: until convergence
11: Output: Matrix W.
3.3 SSNMF Algorithm for Sparse NMF
We are now in a position to present our complete Sequential Sparse NMF (SSNMF) algo-
rithm. By combining Algorithms 1, 2 and 3, we obtain SSNMF (Algorithm 4).
Algorithm 4ssnmf(X,W,H)
1: repeat
2: W = sequential-pass(X,W,H)
3: H = nnls-mult(X,W,H)
4: until convergence
5: Output: Matrices W,H.
6
4 Implementation Issues
For clarity of exposition, we presented the plain vanilla version of our SSNMF Algorithm 4.
We now describe some of the actual implementation details.
• Initialization: Generate a positive random vector v of size m and obtain z =
Sparse-opt(v,k) where k = √m−α√m−1 (from equation (4)). Use the solu-
tion zand its random permutations to initialize matrix W. Initialize the matrix H
to uniform random entries in [0 ,1].
•Incorporating faster solvers: We use multiplicative updates for a fair comparison
with NMFSC and TPC. However, we can use other NNLS solvers [22, 18, 6, 15]
to solve for matrix H. Empirical results (not reported here) show that this further
speeds up the SSNMF algorithm.
•Termination: In our experiments, we ﬁx the number of alternate updates or equiv-
alently the number of times we update matrix W. Other approaches include spec-
ifying total running time, relative change in objective value between iterations or
approximate satisfaction of KKT conditions.
•Sparsity constraints: We have primarly considered the sparse NMF model as for-
mulated by Hoyer [14]. This has been generalized by Heiler and Schn¨ orr [11] by
relaxing the sparsity constraints to lie in user-deﬁned intervals. Note that, we can
handle this formulation [11] by making a trivial change to Algorithm 3.
5 Experiments and Discussion
In this section, we compare the performance of our algorithm with the state-of-the-art
NMFSC and TPC algorithms [14, 11]. Running times for the algorithms are presented when
applied to one synthetic and three real-world datasets. Experiments report reconstruction
error (∥X−WH∥F) instead of objective value for convenience of display. For all experiments
on the datasets, we ensure that our ﬁnal reconstruction error is always better than that of
the other two algorithms. Our algorithm was implemented in MATLAB ( http://www.
mathworks.com) similar to NMFSC and TPC. All of our experiments were run on a 3 .2Ghz
Intel machine with 24GB of RAM and the number of threads set to one.
5.1 Datasets
For comparing the performance of SSNMF with NMFSC and TPC, we consider the following
synthetic and three real-world datasets :
•Synthetic: 200 images of size 9 ×9 as provided by Heiler and Schn¨ orr [11] in their
code implementation.
•CBCL face dataset consists of 2429 images of size 19 ×19 and can be obtained at
http://cbcl.mit.edu/cbcl/software-datasets/FaceData2.html.
•ORL: Face dataset that consists of 400 images of size 112 ×92 and can be ob-
tained at http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.
html. http://www.kyb.tuebingen.mpg.de/bethge/vanhateren/iml/. We use
400 of these images in our experiments.
•sMRI: Structural MRI scans of 269 subjects taken at the John Hopkins University
were obtained. The scans were taken on a single 1.5T scanner with the imaging
parameters set to 35mm TR, 5ms TE, matrix size of 256 ×256. We segment
these images into gray matter, white matter and cerebral spinal ﬂuid images, using
the software program SPM5 (http://www.ﬁl.ion.ucl.ac.uk/spm/software/spm5/),
followed by spatial smoothing with a Gaussian kernel of 10 ×10 ×10 mm. This
results in images which are of size 105 ×127 ×46
7
Problem size (omitted scale factor = 100,000)
oyer
Figure 2: Mean running times for Sparse-opt and the Projection-Hoyer are presented for
random problems. The x-axis plots the dimension of the problem while the y-axis has the
running time in seconds. Each of the subﬁgures corresponds to a single sparsity value in
{0.2,0.4,0.6,0.8}. Each datapoint corresponds to the mean running time averaged over 40
runs for random problems of the same ﬁxed dimension.
0.5
1.0
1.5
2.0reconstruction error
0.2
SSNMF
NMFSC
TPC
0.3
 0.4
 0.5
10-2
10-1
100
101
time (seconds)
0.5
1.0
1.5
2.0reconstruction error
0.6
10-2
10-1
100
101
time (seconds)
0.7
10-2
10-1
100
101
time (seconds)
0.8
10-2
10-1
100
101
time (seconds)
0.9
Figure 3: Running times for SSNMF and NMFSC and TPC algorithms on the synthetic
dataset where the sparsity values range from 0 .2 to 0 .8 and number of features is 5. Note
that SSNMF and NMFSC are over an order of magnitude faster than TPC.
5.2 Comparing Performances of Core Updates
We compare our Sparse-opt (Algorithm 2) routine with the competing Projection-
Hoyer [14]. In particular, we generate 40 random problems for each sparsity constraint
in {0.2,0.4,0.6,0.8}and a ﬁxed problem size. The problems are of size 2 i ×100 where
i takes integer values from 0 to 12. Input coeﬃcients are generated by drawing samples
uniformly at random from [0 ,1]. The mean values of the running times for Sparse-opt and
the Projection-Hoyer for each dimension and corresponding sparsity value are plotted in
Figure 2.
We compare SSNMF with SSNMF+Proj on the CBCL dataset. The algorithms were run
with rank set to 49. The running times are shown in Figure 6. We see that in low-dimensional
datasets, the diﬀerence in running times are very small.
5.3 Comparing Overall Performances
SSNMF versus NMFSC and TPC: We plot the performance of SSNMF against
NMFSC and TPC on the synthetic dataset provided by Heiler and Schn¨ orr [11] in Fig-
ure 3. We used the default settings for both NMFSC and TPC using the software provided
8
5
6
7
8
9
10reconstruction error
0.1
SSNMF
NMFSC
 0.2
 0.3
 0.4
101
102
103
104
time (seconds)
5
6
7
8
9
10reconstruction error
0.5
101
102
103
104
time (seconds)
0.6
101
102
103
104
time (seconds)
0.7
101
102
103
104
time (seconds)
0.8
Figure 4: Convergence plots for the ORL dataset with sparsity from [0.1,0.8] for the NMFSC
and SSNMF algorithms. Note that we are an order of magnitude faster, especially when
the sparsity is higher.
300
400
500
600
700
800
900reconstruction error
0.1
SSNMF
NMFSC
 0.2
 0.3
 0.4
103
104
time (seconds)
300
400
500
600
700
800
900reconstruction error
0.5
103
104
time (seconds)
0.6
103
104
time (seconds)
0.7
103
104
time (seconds)
0.8
Figure 5: Running times for SSNMF and NMFSC algorithms for the sMRI dataset with
rank set to 40 and sparsity values of αfrom 0.1 to 0.8. Note that for higher sparsity values
we converged to a lower reconstruction error and are also noticeably faster than the NMFSC
algorithm.
by the authors. Our experience with TPC was not encouraging on bigger datasets and hence
we show its performance only on the synthetic dataset. It is possible that the performance
of TPC can be improved by changing the default settings but we found it non-trivial to do
so.
SSNMF versus NMFSC: To ensure fairness, we removed logging information from
NMFSC code [14] and only computed the objective for equivalent number of matrix updates
as SSNMF. We do not plot the objective values at the ﬁrst iteration for convenience of
display. However, they are the same for both algorithms because of the shared initialization
. We ran the SSNMF and NMFSC on the ORL face dataset. The rank was ﬁxed at 25
in both the algorithms. Also, the plots of running times versus objective values are shown
in Figure 4 corresponding to sparsity values ranging from 0 .1 to 0 .7. Additionally, we ran
our SSNMF algorithm and NMFSC algorithm on a large-scale dataset consisting of the
structural MRI images by setting the rank to 40. The running times are shown in Figure 5.
9
40
50
60
70
80
90
100reconstruction error
0.2
SSNMF
SSNMF+Proj
0.3
 0.4
 0.5
100
101
102
time (seconds)
40
50
60
70
80
90
100reconstruction error
0.6
100
101
102
time (seconds)
0.7
100
101
102
time (seconds)
0.8
100
101
102
time (seconds)
0.9
Figure 6: Running times for SSNMF and SSNMF+Proj algorithms for the CBCL face
dataset with rank set to 49 and sparsity values ranging from 0 .2 to 0.9
5.4 Main Results
We compared the running times of our Sparse-opt routine versus the Projection-Hoyer and
found that on the synthetically generated datasets we are faster on average.
Our results on switching the Sparse-opt routine with the Projection-Hoyer did not slow
down our SSNMF solver signiﬁcantly for the datasets we considered. So, we conclude that
the speedup is mainly due to the sequential nature of the updates (Algorithm 3).
Also, we converge faster than NMFSC for fewer number of matrix updates. This can be
seen by noting that the plotted points in Figures 4 and 5 are such that the number of matrix
updates are the same for both SSNMF and NMFSC. For some datasets, we noted a speedup
of an order of magnitude making our approach attractive for computation purposes.
Finally, we note that we recover a parts-based representation as shown by Hoyer [14]. An
example of the obtained features by NMFSC and ours is shown in Figure 7.
(a) sparsity 0.5
 (b) sparsity 0.6
 (c) sparsity 0.75
Figure 7: Feature sets from NMFSC algorithm (Left) and SSNMF algorithm (Right) using
the ORL face dataset for each sparsity value of α in {0.5,0.6,0.75}. Note that SSNMF
algorithm gives a parts-based representation similar to the one recovered by NMFSC.
6 Connections to Related Work
Other SNMF formulations have been considered by Hoyer [13], Mørup et al. [24] , Kim
and Park [17], Pascual-Montano et al. [25] (nsNMF) and Peharz and Pernkopf [26]. SNMF
formulations using similar sparsity measures as used in this paper have been considered for
applications in speech and audio recordings [30, 29].
We note that our sparsity measure has all the desirable properties, extensively discussed
by Hurley and Rickard [16], except for one (“cloning”). Cloning property is satisﬁed when
10
two vectors of same sparsity when concatenated maintain their sparsity value. Dimensions
in our optimization problem are ﬁxed and thus violating the cloning property is not an
issue. Compare this with the L 1 norm that satisﬁes only one of these properties (namely
“rising tide”). Rising tide is the property where adding a constant to the elements of a
vector decreases the sparsity of the vector. Nevertheless, the measure used in Kim and Park
is based on the L1 norm. The properties satisﬁed by the measure in Pascual-Montano et al.
are unclear because of the implicit nature of the sparsity formulation.
Pascual-Montano et al. [25] claim that the SNMF formulation of Hoyer, as given by prob-
lem (3) does not capture the variance in the data. However, some transformation of the
sparsity values is required to properly compare the two formulations [14, 25]. Preliminary
results show that the formulation given by Hoyer [14] is able to capture the variance in the
data if the sparsity parameters are set appropriately. Peharz and Pernkopf [26] propose
to tackle the L 0 norm constrained NMF directly by projecting from intermediate uncon-
strained solutions to the required L 0 constraint. This leads to the well-known problem of
getting stuck in local minima. Indeed, the authors re-initialize their feature matrix with
an NNLS solver to recover from the local suboptimum. Our formulation avoids the local
minima associated with L 0 norm by using a smooth surrogate.
7 Conclusions
We have proposed a new eﬃcient algorithm to solve the sparse NMF problem. Experiments
demonstrate the eﬀectiveness of our approach on real datasets of practical interest. Our
algorithm is faster over a range of sparsity values and generally performs better when the
sparsity is higher. The speed up is mainly because of the sequential nature of the updates
in contrast to the previously employed batch updates of Hoyer. Also, we presented an exact
and eﬃcient algorithm to solve the problem of maximizing a linear objective with a sparsity
constraint, which is an improvement over the heuristic approach in Hoyer.
Our approach can be extended to other NMF variants [13]. Another possible application is
the sparse version of nonnegative tensor factorization. A diﬀerent research direction would
be to scale our algorithm to handle large datasets by chunking [23] and/or take advantage
of distributed/parallel computational settings [3].
Acknowledgement
The ﬁrst author would like to acknowledge the support from NIBIB grants 1 R01 EB 000840
and 1 R01 EB 005846. The second author was supported by NIMH grant 1 R01 MH076282-
01. The latter two grants were funded as part of the NSF/NIH Collaborative Research in
Computational Neuroscience Program.
References
[1] Sanjeev Arora, Rong Ge, Ravindran Kannan, and Ankur Moitra. Computing a nonneg-
ative matrix factorization – provably. In Proceedings of the 44th symposium on Theory
of Computing, STOC ’12, pages 145–162, New York, NY, USA, 2012. ACM.
[2] Michael W Berry, Murray Browne, Amy N Langville, V Paul Pauca, and Robert J
Plemmons. Algorithms and applications for approximate nonnegative matrix factor-
ization. Computational Statistics & Data Analysis, 52(1):155–173, 2007.
[3] Joseph K. Bradley, Aapo Kyrola, Danny Bickson, and Carlos Guestrin. Parallel coor-
dinate descent for L1-regularized loss minimization. In ICML, pages 321–328, 2011.
[4] G. Buchsbaum and O. Bloch. Color categories revealed by non-negative matrix factor-
ization of munsell color spectra. Vision research, 42(5):559–563, 2002.
[5] Yunmei Chen and Xiaojing Ye. Projection onto a simplex. arXiv preprint
arXiv:1101.6081, 2011.
11
[6] A. Cichocki and A. H. Phan. Fast local algorithms for large scale nonnegative matrix
and tensor factorizations. IEICE Transactions on Fundamentals of Electronics, 92:
708–721, 2009.
[7] J. E. Cohen and U. G. Rothblum. Nonnegative ranks, decompositions, and factor-
izations of nonnegative matrices. Linear Algebra and its Applications, 190:149–168,
1993.
[8] Chris Ding, Tao Li, Wei Peng, and Haesun Park. Orthogonal nonnegative matrix t-
factorizations for clustering. In Proceedings of the 12th ACM SIGKDD international
conference on Knowledge discovery and data mining, KDD ’06, pages 126–135, New
York, NY, USA, 2006. ACM.
[9] David Donoho and Victoria Stodden. When does non-negative matrix factorization give
a correct decomposition into parts? In Sebastian Thrun, Lawrence Saul, and Bernhard
Sch¨ olkopf, editors,Advances in Neural Information Processing Systems 16. MIT Press,
Cambridge, MA, 2004.
[10] John Duchi, Shai Shalev-Shwartz, Yoram Singer, and Tushar Chandra. Eﬃcient pro-
jections onto the l1-ball for learning in high dimensions. In Proceedings of the 25th
international conference on Machine learning, pages 272–279, 2008.
[11] Matthias Heiler and Christoph Schn¨ orr. Learning sparse representations by non-
negative matrix factorization and sequential cone programming. The Journal of Ma-
chine Learning Research, 7:2006, 2006.
[12] T. Hofmann. Unsupervised learning by probabilistic latent semantic analysis. Machine
Learning, 42(1):177–196, 2001.
[13] P. O. Hoyer. Non-negative sparse coding. In Neural Networks for Signal Processing,
2002. Proceedings of the 2002 12th IEEE Workshop on, pages 557–565, 2002.
[14] Patrik O. Hoyer. Non-negative matrix factorization with sparseness constraints. J.
Mach. Learn. Res., 5:1457–1469, December 2004.
[15] C. J. Hsieh and I. Dhillon. Fast coordinate descent methods with variable selection for
non-negative matrix factorization. ACM SIGKDD Internation Conference on Knowl-
edge Discovery and Data Mining, page xx, 2011.
[16] Niall Hurley and Scott Rickard. Comparing measures of sparsity. IEEE Trans. Inf.
Theor., 55:4723–4741, October 2009.
[17] Hyunsoo Kim and Haesun Park. Sparse non-negative matrix factorizations via alter-
nating non-negativity-constrained least squares for microarray data analysis. Bioinfor-
matics, 23(12):1495–1502, 2007.
[18] Jingu Kim and Haesun Park. Toward faster nonnegative matrix factorization: A new
algorithm and comparisons. Data Mining, IEEE International Conference on, 0:353–
362, 2008.
[19] W. H. Lawton and E. A. Sylvestre. Self modeling curve resolution. Technometrics,
pages 617–633, 1971.
[20] D. D. Lee and H. S. Seung. Learning the parts of objects by non-negative matrix
factorization. Nature, 401(6755):788–791, October 1999.
[21] Daniel D. Lee and Sebastian H. Seung. Algorithms for non-negative matrix factoriza-
tion. In NIPS, pages 556–562, 2000.
[22] Chih-Jen Lin. Projected gradient methods for nonnegative matrix factorization. Neural
Comp., 19(10):2756–2779, October 2007.
[23] Julien Mairal, Francis Bach, Jean Ponce, and Guillermo Sapiro. Online learning for
matrix factorization and sparse coding. The Journal of Machine Learning Research,
11:19–60, 2010.
12
[24] Morten Mørup, Kristoﬀer Hougaard Madsen, and Lars Kai Hansen. Approximate L0
constrained non-negative matrix and tensor factorization. In ISCAS, pages 1328–1331,
2008.
[25] A. Pascual-Montano, J. M. Carazo, K. Kochi, D. Lehmann, and R. D. Pascual-Marqui.
Nonsmooth nonnegative matrix factorization (nsNMF). Pattern Analysis and Machine
Intelligence, IEEE Transactions on, 28(3):403–415, March 2006.
[26] R. Peharz and F. Pernkopf. Sparse nonnegative matrix factorization withl0-constraints.
Neurocomputing, 2011.
[27] M. N. Schmidt and R. K. Olsson. Single-channel speech separation using sparse non-
negative matrix factorization. In International Conference on Spoken Language Pro-
cessing (INTERSPEECH), volume 2, page 1. Citeseer, 2006.
[28] Fabian J Theis, Kurt Stadlthanner, and Toshihisa Tanaka. First results on uniqueness
of sparse non-negative matrix factorization. In Proceedings of the 13th European Signal
Processing Conference (EUSIPCO05), 2005.
[29] Tuomas Virtanen. Monaural sound source separation by nonnegative matrix factor-
ization with temporal continuity and sparseness criteria. Audio, Speech, and Language
Processing, IEEE Transactions on, 15(3):1066–1074, 2007.
[30] Felix Weninger, Jordi Feliu, and Bjorn Schuller. Supervised and semi-supervised sup-
pression of background music in monaural speech recordings. In Acoustics, Speech and
Signal Processing (ICASSP), 2012 IEEE International Conference on, pages 61–64.
IEEE, 2012.
[31] W. Xu, X. Liu, and Y. Gong. Document clustering based on non-negative matrix
factorization. In Proceedings of the 26th annual international ACM SIGIR conference
on Research and development in informaion retrieval, pages 267–273. ACM, 2003.
Appendix
Bi-Sparse NMF
In some applications, it is desirable to set the sparsity on both matrix factors. However,
this can lead to the situation where the variance in the data is poorly captured [25]. To
ameliorate this condition, we formulate it as the following optimization problem and call it
as Bi-Sparse NMF:
min
W,H,D
1
2∥X −WDH∥2
F
s.t.W ≥0,H ≥0,D ≥0
∥Wj∥2 = 1,sp(Wj) = α,∀j ∈{1,··· ,r}
∥Hi∥2 = 1,sp(Hi) = β,∀i∈{1,··· ,r}(6)
where D is a r×r matrix. In the above formulation, we constrain the L 2 norms of the
columns of matrix W to unity. Similarly, we constrain the L 2 norms of rows of matrix H
to be unity. This scaling is absorbed by the matrix D. Note that this formulation with the
matrix D constrained to be diagonal is equivalent to the one proposed in Hoyer when both
the matrix factors have their sparsity speciﬁed.
We can solve for the matrix D with any NNLS solver. A concrete algorithm is the one
presented in Ding et al. and is reproduced here for convenience (Algorithm 5). If D is a
diagonal matrix, we only update the diagonal terms and maintain the rest at zero. Algo-
rithms 1 and 5 can be sped up by pre-computing the matrix products which are unchanged
during the iterations.
Also, the matrix D captures the variance of the dataset when we have sparsity set on both
the matrices W,H.
13
Algorithm 5Diag-mult(X,W,H,D)
repeat
D = D ⊙ W⊤XH
W⊤WDHH⊤
until convergence
Output: Matrix D.
14


Peer review:
[Review]:
This paper proposes new algorithms to minimize the Non-negative Matrix Factorization (NMF) reconstruction error in Frobenius norm subject to additional sparseness constraints (NMFSC) as originally proposed by [R1]. The original method from [R1] to minimize the reconstruction error is a projected gradient descent. While in [R1] a geometrically inspired method is used to compute the projection onto the sparseness constraints, this paper proposes to use Lagrange multipliers instead. To solve the NMFSC problem, the authors propose to update the basis vectors one at a time (therefore their method is called Sequential Sparse NMF or SSNMF), while in ordinary NMF/NMFSC the entire matrix with the basis vectors is updated at once. Experiments are reported that show that SSNMF is one order of magnitude faster compared to the algorithm of [R1].

The paper may only propose more efficient algorithms to solve a known optimization problem instead of proposing new learnable representations, but the approach is interesting and the results are promising. There are however some major issues with the paper:

(1) The sparseness projection of [R1] is essentially a Euclidean projection onto the intersection of a scaled probabilistic simplex (L1 sphere intersected with positive cone) and the scaled unit sphere (in L2 norm). The method of [R1] to compute this projection is an alternating projection algorithm (similar to the Dykstra algorithm for convex sets). The method was proven correct by [R2], and additionally it was shown that the projection is unique almost everywhere. Therefore, the method of [R1] and Algorithm 2 of the paper (Sparse-opt) should almost always compute the same result. In the paper, however, the sparseness projection of [R1] is denoted the 'projection-heuristic' while Sparse-opt is called 'exact', and when the projection of [R1] is used in the SSNMF algorithm instead of Sparse-opt the reconstruction error is no more monotonically decreasing as optimization proceeds. As both projection algorithms should compute the same, the plot should be identical for them when using the same starting points. Section 5.2 of the paper should be enhanced to verify whether both algorithms actually compute the same result and to find the bug that causes this contradiction.

(2) The proposed Algorithm 2 can be considered a (non-trivial) extension of the projection onto a scaled probabilistic simplex as described by [R3] and is a valuable contribution. In the paper, there is however a bug in the execution (which may explain the discrepancies described in Issue (1)): There are no multipliers that enforce the entries of the projection to be non-negative, as would be required by Problem (5) in the paper. Analogously, in Algorithm 2 there is no check in the loop of Line 2 to guarantee the values for lambda and mu produce a feasible (that is non-negative) solution. I implemented the algorithm in Matlab and compared it to the sparseness projection of [R1] (which is freely available on the author's homepage). In the algorithm as given in the paper, p_star always equals m after line 3 and no correct solution to Problem (5) is found in general. If I add the check for a feasible solution, both Sparse-opt and the sparseness projection of [R1] compute numerically equal results. I first suspected there was a typo in the manuscript, but that still would not explain the contradictory results from Section 5.2 of the paper.

On the positive side, I did check the expressions for lambda, mu and obj as given in Algorithm 2, and found them correct. Further, the algorithm is empirically faster than that of [R1], and its run-time is guaranteed theoretically to be at most quasilinear.

Based on the bugfix, I realized that the method from [R4] could be adapted to Sparse-opt to further enhance its run-time efficiency: Set p_star to m before the for loop of line 2 (in case all elements of the projection will be non-zero). Then, after computation of lambda and mu (obj does not need to be computed anymore with this modification), check if a_p < -mu(p) holds. If it does, set p_star to p - 1 and break the for loop. Line 3 of the algorithm should then be omitted. This modification fixes the algorithm, and additionally obj is not needed, and for lambda and mu simple scalars are sufficient to store at most two values of each.

(3) As noted by Paul Shearer and confirmed by the first author of the paper (see public comments), the proof of Theorem 3 is flawed as the arguments there would only apply if the sparseness constraints would induce a convex set (which they don't). I wouldn't have any objections if Theorem 3 and its proof were withdrawn and removed from the manuscript.

Moreover, I verified Algorithm 3 from the paper and found no obvious bugs. I implemented all algorithms and ran them on the ORL face data set and found that SSNMF computes a sparse representation. I did not check what happens without the bugfix for Algorithm 2, though. The authors should definitely fix the major issues and repeat the experiments before publication (judging from the run-time given in Figure 3 and Figure 4 this shouldn't take too long).

There are some minor issues too:
- It should be briefly discussed whether SSNMF could benefit from a multi-threaded implementation as NMF/NMFSC do (in the experiments, the number of threads was set to one).
- Figures should be enlarged and improved such that a difference between the plots is also noticeable when printed in black and white on letter size paper.
- The references should be polished to achieve a consistent style (remove the URLs and ISSNs, don't use two different styles for JMLR ([7] and [10]) and NIPS ([6] and [17]), fix the page of [11], add volume and issue to [15] and [23], add the journal version of [21] unless that citation is withdrawn with Theorem 3, etc.).
- Always cite using numbers in the main text ('cite{}') instead of using only the author names ('citeauthor{}') (e.g. Hoyer, Kim and Park, etc.), because now [9], [10] and [13], [14] could be confused.
- The termination criteria should be described more elaborately for Algorithms 1, 3, and 4.
- Page 2, just after Expression (1): This is only a convex combination if the rows of H are normed (wrt. L1), otherwise it's a conical combination.
- Page 2, just after Expression (2): We use *subscripts* to denote... (missing s). Please also define what H_j^T would mean (is it (H_j)^T or (H^T)_j or something else?).
- It would be nice to add line numbers to all algorithms (some have ones, some don't).
- In Algorithm 3, Line 7: This should probably read G_j^T, as i is not defined here?
- Mention the number of images for the sMRI data set in Section 5.1, and use '' or a footnote for the URL there.
- Cite [7] in the third bullet point in Section 2.2.

References:
[R1] Hoyer. Non-negative Matrix Factorization with Sparseness Constraints. JMLR, 2004, vol. 5, pp. 1457-1469.
[R2] Theis et al. First results on uniqueness of sparse non-negative matrix factorization. EUSIPCO, 2005, vol. 3, pp. 1672-1675.
[R3] Duchi et al. Efficient Projections onto the l1-Ball for Learning in High Dimensions. ICML, 2008, pp. 272-279.
[R4] Chen & Ye. Projection Onto A Simplex. arXiv:1101.6081v2, 2011.