計測自動制御学会東北支部 第 284 回研究集会 (2013.11.8)
資料番号 284-11

Sparse Representation by Nonnegative Matrix
Factorization with the Itakura-Saito Divergence andSparsity Constraint
Zhenni Li∗ Shuxue Ding∗ Yujie Li∗ and Zunyi Tang∗∗* School of Computer Science and Engineering, The University of Aizu,**Faculty of Biomedical Engineering, Osaka Electro-Communication University

キーワード：Nonnegative matrix factorization, Sparse representation, Dictionary learning,
Itakura-Saito divergence, Multiplicative updates
連絡先: 〒965-8580 348E Room, The University of Aizu, Tsuruga, Ikki- Machi,
Aizu-Wakamatsu City,Fukushima, Japan
Zhenni Li，Tel.: (024)237-2799，Fax.: (024)237-2799，Email:lizhenni2012@gmail.com____________________________________
Abstract—In this paper, we propose a novel and efficientdictionary learning method for sparse representation of signals.
The proposed algorithm is based on the nonnegative matrix
factorization (NMF). We adopt the Itakura-Saito (IS) divergenceas the error function and impose ℓ1 -norm as the sparsity
constraint. The error function is quite different from conventional
dictionary learning methods using Euclidean (EUC) distance aserror function. Numerical experiments show that the proposed
algorithm performs better than the other three compared algorithms which all use Euclidean distance as the error function.

1. INTRODUCTION

model of sparse representation of signals, using an overcomplete dictionary matrix W ∈ Rm×r that contains r atoms ofsize m × 1 for columns, signals Y ∈ Rm×n can be describedby sparse and efficient linear combinations of few atoms.
An overcomplete dictionary can lead to sparse representation
where overcomplete means m < r. Y = WH or Y ≈ WH
are two ways to represent Y, where the H ∈ Rr×n is termed
as the coefficient matrix which contains the coefficients forrepresentation of signals Y. The popular EUC distance can bechosen as the error function to measure the approximation error
between Y and WH, namely satisfying ∥ Y − WH ∥2 ≤ ε.

Amazing nonnegative matrix factorization (NMF)[1, 2],which is a method for dimensionality reduction and data analysis [3], has attracted great attention over the past few years.
The NMF has found a wide variety of application, including
machine learning, image processing, blind source separation,m×netc. Given such a representation Y ∈ R+
, the NMF consists
in finding a factorization of the form Y ≈ WH, whereW ∈ Rm×r
is termed as the base matrix, and H ∈ Rr×n
is++termed as the coefficient matrix. The factorization is obtainedby minimizing error function D(Y, W H) between Y and
WH. For D(Y, W H), the Euclidean (EUC) distance and thegeneralized Kullbank-Leibler (KL) divergence are often chosen
in the literature. People note that nonnegativity constraint could
lead to sparse representation which is attracting more and more
attention in recent years. The sparse representations, wheremost coefficients are zero, are proven to be an interesting andpowerful tool for analysis and processing of signals. In the

It is obvious and interesting that dictionary learning, building a dictionary consisting of atoms or subspaces so that aclass of signals can be efficiently and sparsely represented interms of the atoms, is an important topic. In recent practices[4, 5], a learned dictionary has been proved to be criticalfor achieving superior results in the field of signal and image processing. On the other hand, in some applications thenonnegativity of the signals and the dictionary are required,such as the multilateral data analysis [6] and the nonnegative
factorization for recognition [7, 8]. These requirements call forthe dictionary learning method imposed with the nonnegativity,
namely, so-called nonnegative dictionary learning. Extensive
research in this field concentrates mainly on the study of
pursuit dictionary learning algorithms. To some extent, sparse
representation of nonnegative signals is similar to NMF in
the sense that it can increase interpretability and admitting
nonnegative combinations of the coefficient matrix can lead-1-

to sparse results. Unfortunately, the degree of sparseness can
not be controlled and the results may not be enough if only
using the nonnegativity. In order to render a sparser representation, some kinds of sparsity-constrained NMF methods, withdifferent constraints imposed on the matrix factors, have been
proposed. Expansion speaking, ℓ0 -norm, ℓ1 -norm and ℓ2 -norm
constraint have been usually used as the sparsity constraints[9–11]. Since the ℓ0 -norm optimization problem is generally
NP-hard, fortunately it can be replaced by ℓ1 -norm for theconvenience of optimization in the real-world applications.
Some authors also impose sparsity constraints by using ℓ2 norm constraint[12], because of the particularity of the sparse
NMF. The classical methods include NMFSC Hoyer proposed[10] and NMFℓ0 -H [13]. And other classical algorithms is NNKSVD (nonnegative variant of K-SVD) algorithm [14, 17].
These algorithms all use the EUC distance as the measure of
approximation error. In addition, the KL divergence is also
used for the measure error[15].
In this paper, we propose using the IS divergence as themeasure of approximation error. This divergence is obtainedby Itakura-Satio [16] from the maximum likelihood (ML)estimation. The divergence has good properties, in particularscaleinvariant meaning that low energy components of input
signal bear the same relative importance as high energy ones.
And we use the ℓ1 -norm as sparsity constraint to reach thesparsest representation for the coefficient matrix under thecondition of nonnegativity. Then we combine sparsity constraint with IS divergence into a novel algorithm. We used
multiplicative updating to develop an algorithm termed as IS
divergence with ℓ1 -norm sparsity constraint (ISSC). Gradientmultiplicative updating is an efficient way in this case, since it
can easily preserve nonnegativity constraint at each iteration.
Numerical experiments show that the proposed algorithm canrecover almost all aimed dictionary atoms from training data,
which is superior over the other algorithms including NMFSC,
NN-KSVD, and NMFℓ0 -H.
The remaining part of paper is organized as follows. In
section II, we review the relationships between the IS divergence and other error functions used in NMF, especially the
Euclidean distance and the generalized Kullback-Leibler (KL)divergence. In section III, we describe proposed algorithm
ISSC. In Section IV we present the results of numericalexperiments for the proposed algorithm and compare theseresults with those of several other algorithms. Finally, we givethe conclusion and discuss the future work.

2. PROBLEM FORMULATIONThe NMF problem is formulated as follows. Given an input
matrix Y ∈ Rm×n
, where each element is nonnegative. NMF+
aims to find nonnegative matrices W ∈ Rm×rand H ∈ Rr×n++
satisfying the condition of Y = WH or Y ≈ WH, wherer ≪ min(m, n). This problem can be formulated as theminimization of an objective function,minsubject toD(Y|WH)(1)

W ∈ Rm×r, H ∈ Rr×n
+
+

-2-

In order to complete the approximate factorization, we needto define some error functions to measure quality of theapproximation error. Popular choices are the EUC distancewhich can be defined as,
1∥Y − WH∥2F(2)
2
and the generalized KL divergence, also termed as Idivergence, can be defined as)∑(
Yij−Yij +(WH)ij
min DKL (Y|WH) =Yij log(WH)ijij(3)
1 ≤ i ≤ m, 1 ≤ j ≤ n
min DEU C (Y|WH) =Above two error functions are both nonnegative, and if Y =WH the value equal to zero.
In this paper, we consider NMF under other error function
which is the IS divergence, whose expression is given by)∑ ( Yij
Yij
min DIS (Y|WH) =−log−1 (4)(WH)ij(WH)ijij
1 ≤ i ≤ m, 1 ≤ j ≤ n
The three mentioned error functions are the special casesof β-divergence [18–20], which show that β = 2 is the caseof EUC distance, β = 1 is the case of KL divergence, and the
IS divergence is a limiting case when β = 0.

3.

THE PROPOSED ALGORITHM

Although NMF yields sparseness for nonnegativity to someextent, it is believed that much sparser representation canbe learned by imposing sparsity constraints on the matrix
factors. ℓ0 -norm, ℓ1 -norm and ℓ2 -norm constraint have been
usually used as the sparsity constraints. The ℓ0 -norm constraint
is the sparsest of the three. Unfortunately, solving the ℓ0 norm constraint optimization problem is generally NP-hard.
For solving this problem, one can replace ℓ0 -norm by ℓ1 norm constraint. In addition, some authors also impose sparsity
constraints by using ℓ2 -norm, because of the particularity ofthe sparse NMF. However ℓ2 -norm constraint is not sparser
than ℓ0 -norm and ℓ1 -norm constraint. In this paper, we adoptℓ1 -norm as sparsity constraint.
Combining the error function of the IS divergence with ℓ1 norm sparsity constraint, then we propose our problem as thefollowing optimization problem,)∑ ( Yij
Yij
min DIS (Y|WH) =−log−1 +λ ∥ H ∥1(WH)ij(WH)ijij(5)
r×n
subject to W ∈ Rm×r
,
H
∈
R
+
+
∑
where ||H||1 meansi,j |Hi,j | and λ ≥ 0 is a regularization parameter which can be adjusted for controlling thetradeoff between the approximation error and the sparsity of
the coefficient matrix H. One can repeat experiments withdifferent values of λ and determine which value for λ is
optimal according to the output results.

For solving the constrained NMF problem, many algorithms
have been developed and most of them are structured with
iterative strategy, which utilize the fact that the problem canbe reduced into two sequential convex nonnegative problems
about W or H whereas the other of them is regarded as fixedand known. The traditional multiplicative gradient descentapproach [1, 19] consists in updating each parameter bymultiplying its value at the previous iteration by a certain
coefficient, which can be computed by the following formulas,
Wij = Wij − φij∂D(W, H)∂Wij

(6)

1 ≤ i ≤ m, 1 ≤ j ≤ r∂D(W, H)
Hij = Hij − ψij
∂Hij

(7)

Algorithm 1 ISSC
Require: Data Matrix Y ∈ Rm×n+
1) Initialize W ∈ Rm×rand H ∈ Rr×n
as random non++negative matrices;
2) Scale rows of H to sum to one;Normalize columns of W to a unit ℓ2 -norm
3) Iterate until converge or stop;[(WH).(−2) ⊙ Y]HT(WH).(−1) HT + δ
WT [(WH).(−2) ]H←H⊙
WT (WH).(−1) + λ + δW←W⊙δ = 10−9 (overcoming the numerical instabilities)
4) Scale rows of H to sum∑to one and normalize columnsm2
of W to a unit ℓ2 -norm: i=1 wij= 1, ∀j.

1 ≤ i ≤ r, 1 ≤ j ≤ n
The gradients of criterion DIS (Y|WH) with regard to W
4. NUMERICAL EXPERIMENTS
and H can be written as,
Dictionary learning for sparse representation of signals has(() T ) been successfully applied in many areas, e.g. data classification∂D(W, H).(−1) T
.(−2)
=
(WH)
H − (WH)⊙Y H∂Wij
ij and image processing. Especially in the latter area, it has been(8) applied for image denoising [22], image inpainting [23], image(() compression [14, 24], and superresolution reconstruction [25],∂D(W, H)=WT (WH).(−1) − WT (WH).(−2) ⊙ Y
etc. In this section, we made experiments by using ISSC
∂Hij)
algorithm for synthetic signals, to test whether this algorithm+ λcan recover the original dictionary and to compare its resultsij(9) with other algorithms.
where ⊙ and .(−2), .(−1) denote element-wise multiplicationand power. If we choose
φij

=

Wij
(
)
(WH).(−1) HT ij

(10)

ψij

=

Hij
(
)
WT (WH).(−1) + λ ij

(11)

Substituting (8),(9),(10),(11) into (6),(7), then the algorithm isobtained as the following updates,(())(WH).(−2) ⊙ Y HT
ij()
Wij ←− Wij ·
(12)
.(−1)
T
(WH)
H
ij
(
(
))
T
.(−2)
W (WH))ij
Hij ←− Hij · (WT (WH).(−1) + λ

4.1Experiment on Synthetic Signals Generated with a Dictionary
For our experiments, we begun with generating a stochastic
nonnegative matrix W of size 20×50 with i.i.d. uniformly
distributed entries. Each column was normalized to the unit
ℓ2 -norm. The stochastic matrix was referred to as the true
dictionary W, which was not used in the learning but only forevaluation. Then we synthesized 1500 test signals of dimension
20, each of which was produced by a linear combination of
three different atoms in the true dictionary W. These test
signals are uniformly distributed i.i.d coefficients in random
and independent locations. We then synthesized 1500 test
signals Y of dimension 20, each of which was produced
by a linear combination of three different atoms in the true
dictionary, with three corresponding coefficients in random
and independent positions. The uniformly distributed noise ofvarying signal-to-noise ratio (SNR) for performance analysis
of noise-robustness was added to the experiments.

ij

(13)
Obviously, we can observe that the criterion is still nonincreasing under updates (12) and (13). In addition, the convergence had been theoretically proven by the paper [21] basedon the expectation maximization (EM) algorithm.
According to the analysis above, the proposed IS divergencewith ℓ1 -norm sparsity constraint (ISSC) is summarized in
Algorithm 1.

-3-

4.2 Applying the ISSC
The initialized dictionary was composed of the signals
selected randomly from 1500 test signals. The corresponding
coefficients were initialized with random entries that were
i.i.d. uniformly distributed and nonnegative. The maximum
number of iterations was set to 200. For ISSC, the sparsityof the coefficient matrices was adjusted via the regularization
parameters λ. The parameter λ could be determined off-line
calibrating. We repeated the experiment with different λ and4.3 Comparison with the Other Algorithms
Since NN-KSVD, NMFSC and NMFℓ0 -H were the threestate-in-art algorithms for nonnegative dictionary learning, we
compared our algorithm with these algorithms. We executed
NMFSC, NN-KSVD, NMFℓ0 -H using the same test signals
with ISSC, respectively. Then it is easy to estimate learned W
and evaluate every accuracy by comparing with the true W.
For the three algorithms, the initialized dictionary matrices ofsize 20×50 were composed of the randomly selected parts ofthe test signals. Note that for NMFSC the corresponding coefficient matrices were initialized with i.i.d. uniformly distributed
random nonnegative entries. NN-KSVD and NMFℓ0 -H did not
require a specified coefficient matrix, as they could generate
the corresponding coefficient matrix by sparse coding. The
implementation of NN-KSVD1 algorithm is online available.
We executed the NN-KSVD algorithm for a total number of
200 iterations. Matlab code for NMFSC 2 and NMFℓ0 -H3algorithms are also online available. The learning procedure
with NMFSC was stopped after 3000 iterations because it
converged fairly slower than the other algorithms. And the
maximum number of iterations of NMFℓ0 -H algorithm wasfairly set to 200. It was worth noting that, in the experiment,
NN-KSVD and NMFℓ0 -H needed the specified exact numberof non-zero elements in the coefficient matrix (3/50=0.06 forthe case), while NMFSC was executed with a sparsity factorof 0.85 on the coefficients.
4.4 Results of Experiment
The learned dictionaries were compared with the true generated dictionary. The comparisons were done as describedin [14] by sweeping through the columns of the true andthe learned dictionaries and finding the closest column (inℓ2 -norm distance) between the two dictionaries. A distanceless than 0.01 was considered as a success. All trials wererepeated 15 times. In the experiment, the ISSC algorithm
could recovery averaged 9.07%, 69.2%, 88.13% and 93.07%
atoms under the noise levels of 10 dB, 20 dB, 30 dB, and in
the noiseless case, respectively. For NN-KSVD could obtainaveraged 15.7%, 68.0%, 82.9% and 86.5%. While NMFℓ0 -H
could recover 23.7%, 80.8%, 84.9% and 84.0% atoms under
the same conditions. For NMFSC, it recovered only averaged
0.4%, 13.5%, 38.4% and 49.3% atoms. The detailed results ofthe experiment for these algorithms are shown in the Fig. 1.
The proposed ISSC performed best on dictionary learning inlevels of 30dB and noiseless case.

___________________________________________
1 Online available http://www.cs.technion.ac.il/˜elad/software/
2 Online available http://www.cs.helsinki.fi/u/phoyer/contact.html
3 Online available http://www3.spsc.tugraz.at/people/robert-peharz

-4-

100
NMFSC

The rates of recovery (%)

determined which value for λ was optimal according to the
output results. In the our algorithm ISSC λ was set to 0.015.
Then we executed ISSC on the test signals to estimate learned
W and evaluate its accuracy by comparing with the true W.

NN−KSVD

80

NMFl0−H
ISSC

60
40
20
0

10 dB

20 dB
30 dB
No Noise
Noise level and no noise caseFig. 1. Experiment result when apply to the synthetic signal: for each of
the tested algorithms and for each noise level, 15 trials were performed andtheir results were sorted. The averaged recovery rates of learned atoms andcorresponding deviation of recovery rates are displayed. λ = 0.015

5. CONCLUSIONSIn this paper we have presented a novel and efficientnonnegative dictionary learning algorithm which is obtainedby combining IS divergence as the error function and with
ℓ1 -norm as the sparsity constraint. Using IS divergence aserror function is quite distinguished comparing with the conventional algorithms. Results of experiments on dictionary
recovery show that the ISSC algorithm can correctly learn
an overcomplete, nonnegative dictionary on synthetic signals
and further show that the proposed algorithm performs therobustness against different level of noise in comparison withthe other compared algorithms which all using the EUC
distance as the error function. It implies that different error
functions lead to different results. In future work, we will
employ the proposed algorithm for inpainting, image denoisingand other application.
R EFERENCES[1] D. D. Lee and H. S. Seung, “Learning the parts of objects
by nonnegative matrix factorization,” Nature, vol. 401,
pp. 788–791, 1999.
[2] A. Cichocki and A.-H. Phan, “Fast local algorithms forlarge scale nonnegative matrix and tensor factorizations,”
IEICE Trans. on Fundamentals of Electronics, vol. E92A, no. 3, pp. 708–721, 2009.
[3] V. P. Pauca, J. Piper, and R. J. Plemmons, “Nonnegative
matrix factorization for spectral data analysis,” Linear
Algebra and its Applications, vol. 416, no. 1, pp. 29–47,
2006.
[4] M. Plumbley, T. Blumensath, L. Daudet, R. Gribonval,
and M. Davies, “Sparse representations in audio andmusic: From coding to source separation,” Proceedingsof the IEEE, vol. 98, no. 6, pp. 995–1005, Jun 2010.
[5] M. Elad, M. Figueiredo, and Y. Ma, “On the role of sparseand redundant representations in image processing,” Proceedings of the IEEE, vol. 98, no. 6, pp. 972–982, 2010.

[6] L. Miao and H. Qi, “Endmember extraction from highly
mixed data using minimum volume constrained nonnegative matrix factorization,” IEEE Transactions on Geoscience and Remote Sensing, vol. 45, no. 3, pp. 765–777,
2007.
[7] S. Li, X. Hou, H. Zhang, and Q. Cheng, “Learning
spatially localized, parts-based representation,” in Proceedings of the 2001 IEEE Computer Society Conference
on Computer Vision and Pattern Recognition (CVPR
2001), vol. 1, 2001, pp. 207–212.
[8] I. Kotsia, S. Zafeiriou, and I. Pitas, “A novel discriminant non-negative matrix factorization algorithm withapplications to facial image characterization problems,”
IEEE Transactions on Information Forensics and Security, vol. 2, no. 3, pp. 588–595, Sep 2007.
[9] R. Peharz and F. Pernkopf, “Sparse nonnegative matrix factorization using ℓ0 -constraints,” Neurocomputing,vol. 80, pp. 38–46, Mar 2012.
[10] P. O. Hoyer, “Non-negative matrix factorization withsparseness constraints,” Journal of Machine Learning
Research, vol. 5, pp. 1457–1469, 2004.
[11] V. P. Pauca, F. Shahnaz, M. W. Berry, and R. J. Plemmons, “Text mining using non-negative matrix factorizations,”in Proc. of the Fourth SIAM International Conference on
Data Mining, 2004, pp. 452–456.
[12] F. Shahnaz, M. W. Berry, V. P. Pauca, and R. J. Plemmons, “Document clustering using nonnegative matrix
factorization,” Information Processing & Management,
vol. 42, no. 2, pp. 373–386, 2006.
[13] R. Peharz, M. Stark, and F. Pernkopf, “Sparse nonnegative matrix factorization using ℓ0 -constraints,” pp. 83–88,
Sep. 2010.
[14] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An
algorithm for designing overcomplete dictionaries forsparse representation,” IEEE Trans. on Signal Processing,vol. 54, no. 11, pp. 4311–4322, Nov. 2006.
[15] W. Liu, N. Zheng, and X. Lu, “Non-negative matrix
factorization for visual coding,” in IEEE Conference on
Acoustics, Speech, and Signal Processing, April 2003, pp.
293–296.
[16] F. Itakura and S. Saito, “Analysis synthesis telephonybased on the maximum likelihood method,” In Proc 6th
International Congress on Acoustics, pp. C–17 – C–20,
1968.
[17] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD and its
nonnegative variant for dictionary design,” vol. 5914, pp.
327–339, Jul. 2005.
[18] S. Eguchi and Y. Kano, “Robustifying maximum likelihood estimation,” Technical report, Institute of Statistical
Mathematics, 2001.
[19] A. Cichocki, R. Zdunek, and S. Amari, “Csiszárs divergences for non-negative matrix factorization : Family of
new algorithms,” Proc. Int. Conf. Independent Component Analysis and Blind Signal Separation, p. 32C39,
2006.
[20] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative

matrix factorization with the itakura-saito divergence:With application to music analysis,” Neural Computation,
vol. 21, no. 3, pp. 793–830, 2009.
[21] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono,and S. Sagayama, “Convergence-guaranteed multiplicative algorithms for nonnegative matrix factorization with
β-divergence,” in Machine Learning for Signal Processing (MLSP), 2010 IEEE International Workshop on,
2010, pp. 283–288.
[22] M. Elad and M. Aharon, “Image denoising via sparseand redundant representations over learned dictionaries,”
IEEE Transactions on Image Processing, vol. 15, no. 12,
pp. 3736–3745, Dec 2006.
[23] M. Fadili, J.-L. Starck, and F. Murtagh, “Inpainting andzooming using sparse representations,” The Computer
Journal, vol. 52, no. 1, pp. 64–79, 2009.
[24] J. Zepeda, C. Guillemot, and E. Kijak, “Image compression using sparse representations and the iteration-tuned
and aligned dictionary,” IEEE Journal of Selected Topics
in Signal Processing, vol. 5, no. 5, pp. 1061 –1073, Sept.
2011.
[25] J. Yang, J. Wright, T. Huang, and Y. Ma, “Image superresolution as sparse representation of raw image patches,”in IEEE Conference on Computer Vision and Pattern
Recognition, 2008. CVPR 2008., June 2008, pp. 1–8.

-5-