2014 International Symposium on Nonlinear Theory and its Applications
NOLTA2014, Luzern, Switzerland, September 14-18, 2014

New Update Rules based on Kullback-Leibler, Gamma, and Rényi Divergences
for Nonnegative Matrix Factorization
Masato Seki† and Norikazu Takahashi†,‡† Graduate School of Natural Science and Technology, Okayama University
3–1–1 Tsushima-naka, Kita-ku, Okayama 700–8530, Japan‡ Institute of Systems, Information Technologies and Nanotechnologies
2–1–22 Momochihama, Sawara-ku, Fukuoka 814–0001 Japan
Email: seki@momo.cs.okayama-u.ac.jp, takahashi@cs.okayama-u.ac.jp
Abstract—Three new update rules based on KullbackLeibler divergence, γ-divergence and Rényi divergence fornonnegative matrix factorization are presented in this paper. An important advantage of these update rules is thatthey are globally convergent in the sense that any sequenceof solutions contains at least one convergent subsequenceand the limit of any convergent subsequence is a stationary
point of the corresponding optimization problem.
1. Introduction
Nonnegative matrix factorization (NMF) [1] is a technique to decompose a given nonnegative matrix X ∈ Rm×n+into two nonnegative matrices W ∈ Rm×rand H ∈ Rr×n++such thatX ≈ WH
(1)
where R+ denotes the set of nonnegative numbers and r
is a positive integer less than min{m, n}. NMF is usuallyformulated as an optimization problem to minimize an error function subject to the nonnegativity constraint on eachvariable. As a simple and eﬃcient method for finding local
optimal solutions of such problems, multiplicative updates
developed by Lee and Seung [1, 2] are widely used.
Recently, Yang and Oja [3] proposed a unified methodfor constructing multiplicative update rules for NMF, andderived eleven update rules from eleven error functions.
However, all of these rules have a common serious problem that the global convergence is not guaranteed. By theglobal convergence, we mean that any sequence of solutions has at least one convergent subsequence and the limitof any convergent subsequence is a stationary point of thecorresponding optimization problem [4]. One of the reasons for this is that the update rules are not well-defined [5].
In order to avoid this problem, Takahashi et al. [6] applied a simple modification technique proposed by Gillis
and Glineur [7] to the eleven multiplicative update rules
and studied their global convergence. As a result, theyproved that eight among the eleven update rules are globally convergent. However, they did not prove the global
convergence of the remaining three update rules, that are
based on Kullback-Leibler divergence, γ-divergence andRéyni divergence (see Table 1).
In this paper, we propose new update rules based on theabove-mentioned three divergences, and prove their global
convergence. The key idea is not to use each divergence
directly as an error function but to add a penalty term to it.
This is based on the observation that if each divergence isdirectly used as an error function then we can increase thevalues of variables as much as we want while keeping thevalue of the error function fixed.
2. Multiplicative Updates for NMF
Suppose that we are given a nonnegative matrix X ∈and a positive integer r < min{m, n}. We will assume
Rm×n+
throughout this paper that every row and column of X has
at least one nonzero entry. Then, the problem of finding W
and H in (1) is formulated asmaximize D(W , H)
subject to W ≥ Om×r , H ≥ Or×n(2)where D(W , H) is an error function and Om×r (Or×n ,resp.) is the m×r (r×n, resp.) zero matrix. Throughout thispaper, matrix inequalities are understood componentwise.
Multiplicative updates are eﬃcient iterative methods forfinding local optimal solutions of the constrained optimization problem (2). Lee and Seung [1, 2] first derived two
multiplicative updates based on Euclidean distance and Idivergence by making use of auxiliary functions. This approach was recently generalized by Yang and Oja [3] sothat it can be applied to many error functions. In fact, they
derived eleven multiplicative updates including the ones ofLee and Seung by using the method summarized below.

This work was partially supported by KAKENHI 24560076 and23310104, and by Proactive Response Against Cyber-attacks ThroughInternational Collaborative Exchange (PRACTICE), Ministry of Internal
Aﬀairs and Communications, Japan.

- 48 -

1. If the error function D(W , H) contains a logarithmic
function, rewrite it in a generalized polynomial form
by using the relationship:zδ − 1.
δ→0+δ

ln z = limTable 1: Error functions based on Kullback-Leibler, gamma, and Rényi divergences [3].
Error function
Kullback-Leibler divergence
γ-divergenceRényi divergence

D(W , H)
∑Xi j∑
Xi j ln(WH)/ijab (W H)abij




 ∑

∑


 ∑1γ
1+γ 1+γ ln 

Xi j (W H)i j  (γ , −1, 0)
X  + γ ln  (W H)i j  − (1 + γ) ln γ(1 + γ)   i j i j ijij ()ρ ()1−ρ ∑X(WH)1ijij∑∑ln  (ρ > 0, ρ , 1)
ρ−1 Xab(W H)abij

ababf , H)f :2. Construct an auxiliary function D̄(W , H, W
r×n
r×n
m×rm×rR++ × R++ × R++ × R++ → R of the error function D(W , H), where R++ denotes the set of positive
numbers, that satisfies
f > Om×r , H
f > Or×n ,∀W > Om×r , H > Or×n , Wf , H)f ≥ D(W , H)
D̄(W , H, Wand
∀W > Om×r , H > Or×n ,D̄(W , H, W , H) = D(W , H) .
f and Hf be any positive matrices and let W ∗3. Let W
be the unique minimizer of the function ϕ(W ) =fWf , H).
f Then W ∗ is expressed as a funcD̄(W , H,ikf and Hf which gives a multiplicative updatetion of W
f and Hf be any positive matrirule for Wik . Also, let W∗
ces and let H be the unique minimizer of the function
f , H, Wf , H).
f Then H ∗ is expressed as
ψ(H) = D̄(Wkjf and H
f which gives a multiplicative
a function of W
update rule for Hk j .
Table 2 shows three multiplicative update rules obtainedby applying this method to the error functions in Table 1.
3. Modified Updates and Their Global Convergence
Every multiplicative update rules presented in [3] is notdefined for all pairs of nonnegative matrices W and H.
This is because each update rule contains a rational function of W and H and the denominator can become zero(see Table 2). A simple way to avoid this problem is to use
the modified update rule [7] which can be expressed as()Wik(l+1) = max ϵ, fik (W (l) , H (l) ) ,(3)()(l+1)
Hk j= max ϵ, gk j (W (l+1) , H (l) ) ,(4)where fik and gk j represent the original updates for Wik and
Hk j , respectively, and ϵ is a small positive constant specified by the user. In the following, we will focus our attention on the modified update rule. Note that, in so doing, weshould not consider the original optimization problem (2)but its modified version:
maximize D(W , H)
(5)subject to W ≥ ϵ1m×r , H ≥ ϵ1r×n
where 1m×r (1r×n , resp.) is the m × r (r × n, resp.) matrix
consisting of all ones. The feasible region of (5) will bedenoted by Fϵ .
Takahashi et al. [6] have recently presented a suﬃcientcondition on the auxiliary function for the modified update rule to be globally convergent in the sense that any
sequence {(W (l) , H (l) )}∞l=0 generated by (3) and (4) with(W (0) , H (0) ) ∈ Fϵ contains at least one convergent subsequence and the limit of any convergent subsequence isa stationary point of (5). Also, they have pointed out thatthe suﬃcient condition is not satisfied for the three multiplicative update rules shown in Table 2. The main reason
for this is that the boundedness of the sequence of solutions
generated by (3) and (4) is not guaranteed [5].
4. New Updates based on Three Divergences
In this section, we propose to modify error functions inTable 1 and then derive fik and gk j by using the method ofYang and Oja [3] so that the modified update rule expressedby (3) and (4) is globally convergent.
By looking at the formulas in Table 1, we can easily seethat the value of D(W , H) does not change if W and H
are multiplied by any positive scalars. In other words, wecan increase the values of nonzero entries of W and H asmuch as we want, while keeping the value of D(W , H)
fixed. This may be the reason why boundedness of solutions is not guaranteed. As a simple way to keep W andH bounded, we propose to add a penalty term2∑C ∑
Xi j −(W H)i j 2 ijij
to each error function, where C is a positive constant. Forexample, our new error function based on Kullback-Leibler
divergence is given by∑
Xi j∑
D(W , H) =Xi j ln(WH)/ijab (W H)abij- 49 -

Table 2: Multiplicative update rules [3] for error functions in Table 1. Z = (Zi j ) is defined by Zi j = Xi j /(W H)i j .
Error function
Kullback-Leibler divergence
γ-divergenceRényi divergence

Multiplicative updates for Wik(ZH T )ik ∑Wiknew = Wik ∑(W H)abj Hk jab∑∑1+γ η1
 j Xi j (W H)γ−1
i j Hk j 1+γ , if γ > 0ab (W H)ab
new


,
η=·Wik = Wik  ∑∑γγ 
1

(W H)i j Hk j, if γ < 0 and γ , −1ab Xab (W H)ab1−γ
 ∑ jρη∑
1ZHkj,ifρ>1j(WH)ijabab , η =  ρ
Wiknew = Wik  ∑·∑
 1, if 0 < ρ < 1X ρ (W H)1−ρ
j Hk jab2∑C ∑+ 
Xi j −(W H)i j  .
2 ij
ij

ab

previous section.

Let us now apply the unified method of Yang and Oja to thiserror function. In the first step, we rewrite the error function in the form of a generalized polynomial as follows:D(W , H)=ab

Theorem 1 Let ϵ and C be any positive numbers. For any
initial matrices W (0) ≥ ϵ1m×r and H (0) ≥ ϵ1r×n , the sequence {(W (l) , H (l) )}∞l=0 generated by any of the three update rules given in the previous section is bounded.
Proof We prove this by using Lemma 1 of [5]. To do so,
it suﬃces for us to show that there exist constants cik (> 0)and νik (< 1) such thatδ[ ∑∑1 ∑
Xi j (W H)δi j
Xab  (W H)i j  −lim
δ→0+ δij
ij
ab
2 ]

∑
δC ∑
+
Xi j −(W H)i j  + constant2 ijij∀W ≥ ϵ1m×r , ∀H ≥ ϵ1r×n , fik (W , H) ≤ cik Wikνik

where “constant” represents terms independent of W andH. In the second step, we obtain an auxiliary function:f , H)f
D̄(W , H, Wδ−1[ ∑ ∑ Wik2 Hk2 j∑1 δff
= lim
Xab  (W H)ab 
δ→0+ δ 2e ei jk Wik Hk jabab∑f H)f δ−1 (We ik Hek j )1−δ W δ H δ−
Xi j (Wijik k jfor all pairs of i and k and there exist constants dk j (> 0)and µk j (< 1) such thatµ∀W ≥ ϵ1m×r , ∀H ≥ ϵ1r×n , gk j (W , H) ≤ dk j Hk kj j (7)for all pairs of k and j. Let us consider the new update rule
based on Kullback-Leibler divergence. Since
fik (W , H)∑∑
 12∑
 j Xi j (W H)−1j Hk j i j Hk j + C ab Xab∑∑
< Wik C ab (W H)ab j Hk j
 1 ∑
1
∑∑
 rϵ 2 j Xi j Hk j + C ab Xab j Hk j  2∑∑< Wik C ab (W H)ab j Hk j
1∑∑ 12 Xi j Hk j + C ab Xab Hk j  2rϵ∑∑= Wik C ab (W H)ab l Hkl j1∑∑ 12 Xi j Hk j + C ab Xab Hk j  2rϵ< Wik CWHnϵikkjj
1∑ 12 Xi j + C ∑ab Xab  21 rϵ ,= Wik2 Cnϵi jk2 2
δC ∑ f f ∑ Wik Hk j+(W H)abe e
2 abi jk Wi j Hk j−C∑ab(6)]∑1−δ δ δeeXab(Wi j Hk j ) Wik Hk j + constant .
i jkIn the last step, we find the unique minimizer of the function ϕ(W ) = D̄(W , H (l) , W (l) , H (l) ), from which we obtain fik (W (l) , H (l) ). Similarly, we find the unique minimizer of the function ψ(H) = D̄(W (l+1) , H, W (l+1) , H (l) ),from which we obtain gk j (W (l+1) , H (l) ).
By using the same procedure as above, we can obtain fikin (3) and gk j in (4) for each of the error functions basedon γ-divergence and Rényi divergence. Explicit formulasof fik and gk j are given in Table 3.
5. Global Convergence of New Updates
Let us consider the boundedness of the sequence of solutions generated by the three update rules obtained in thejthe condition (6) is satisfied. In the same way, we can showthat the condition (7) is also satisfied. As for the update
rules based on γ-divergence and Rényi divergence, we canalso prove that both (6) and (7) are satisfied. However, weomit the details due to lack of space.
□- 50 -

Table 3: fik and gk j in new update rules.
Error function
Kullback-Leibler divergence

γ-divergenceRényi divergence

fik (W , H) and gk j (W , H)
 21∑
∑
∑
−1

j Xi j (W H)i j Hk j + Cab Xabj Hk j

∑∑
∑
∑
fik (W , H) = Wik  ∑−1X((WH))
H+C(WH)
Hababkjabkjababjabj

 21
∑
∑
∑−1

i Xi j (W H)i j Wik + Cab Xabi Wik
∑∑
∑
∑
gk j (W , H) = Hk j  ∑−1X((WH)
W+C(WH)
Wabikabikab ababiabi ∑η∑∑∑ ( ab Xab (W H)γab )−1 j Xi j (W H)γ−1ab Xabj Hk j i j Hk j + C
fik (W , H) = Wik  ∑∑∑1+γ −1 ∑γ((W H)ab )j (W H)i j Hk j + Cab (W H)abj Hk j
 ∑ ab
∑
∑
∑
η
γγ−1
 ( ab Xab (W H)ab )−1 i Xi j (W H)i j Wik + C ab Xab i Wik 
gk j (W , H) = Hk j  ∑∑∑∑γ−1( ab (W H)1+γi (W H)i j Wik + C
ab (W H)abi Wikab )
1/(1 + γ), if γ > 1where η = 1/2,if 0 < γ < 1
 1/(2 − γ), if γ < 0 and γ , −1η ∑∑
∑ρ
ρ
−ρ
−1 ∑
 ( ab Xab(W H)1−ρj Xi j (W H)i j Hk j + Cab Xabj Hk j ab )

∑
∑
∑
∑
fik (W , H) = Wik ( ab (W H)ab )−1 j Hk j + C ab (W H)ab j Hk jη ∑∑∑ρ−ρρ−1 ∑ ( ab Xab(W H)1−ρ
i Wik ab Xabi Xi j (W H)i j Wik + Cab )
∑
∑
∑
∑
gk j (W , H) = Hk j ( ab (W H)ab )−1 i Wik + C ab (W H)ab i Wik
 1/(1 + ρ), if ρ > 1
where η = 
 1/2,if 0 < ρ < 1

We finally consider the global convergence of the new
update rules. Because all auxiliary functions used for obtaining the three update rules satisfy the condition in Theorem 1 of [6], we have the following result.
Theorem 2 Let ϵ and C be any positive numbers. For any
initial matrices W (0) ≥ ϵ1m×r and H (0) ≥ ϵ1r×n , the sequence {(W (l) , H (l) )}∞l=0 generated by any of the three update rules given in the previous section has at least one convergent subsequence and the limit of any convergent subsequence is a stationary point of the problem (5).
6. Conclusion
We have proposed three new error functions based on
Kullback-Leibler divergence, γ-divergence and Rényi divergence for NMF, and derived three new update rules. We
have also shown that all of these update rules have a global
convergence property. However, it is not clear how the solutions obtained by the proposed update rules depends on
two parameters ϵ and C. We leave this as a future problem.

Information Processing Systems (T. K. Leen, T. G. Dietterich, and V. Tresp, eds.), vol. 13, pp. 556–562,
2001.
[3] Z. Yang and E. Oja, “Unified development of multiplicative algorithm for linear and quadratic nonnegative matrix factorization,” IEEE Transactions on Neural Networks, vol. 22, pp. 1878–1891, December 2011.
[4] W. I. Zangwill, Nonlinear programming: A unified approach. Englewood Cliﬀs, NJ: Prentice-Hall, 1969.
[5] J. Katayama, N. Takahashi, and J. Takeuchi, “Boundedness of modified multiplicative updates for nonnegative matrix factorization,” in Proceedings of the Fifth
International Workshop on Computational Advances inMulti-Sensor Adaptive Processing, pp. 252–255, December 2013.

References[6] N. Takahashi, J. Katayama, and J. Takeuchi, “A generalized suﬃcient condition for global convergence ofmodified multiplicative updates for NMF.” to appear
in the Proceedings of 2014 International Symposiumon Nonlinear Theory and its Applications.

[1] D. D. Lee and H. S. Seung, “Learning the parts ofobjects by non-negative matrix factorization,” Nature,
vol. 401, pp. 788–792, 1999.

[7] N. Gillis and F. Glineur, “Nonnegative factorizationand the maximum edge biclique problem,” ArXiv eprints, October 2008.

[2] D. D. Lee and H. S. Seung, “Algorithms for nonnegative matrix factorization,” in Advances in Neural- 51 -