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

A Generalized Suﬃcient Condition for Global Convergence of Modified
Multiplicative Updates for NMF
Norikazu Takahashi†,♯ , Jiro Katayama‡ and Jun’ichi Takeuchi♭,♯† Graduate School of Natural Science and Technology, Okayama University
3–1–1 Tsushima-naka, Kita-ku, Okayama 700–8530, Japan‡ Graduate School of Information Science and Electrical Engineering, Kyushu University
744 Motooka, Nishi-ku, Fukuoka 819–0395 Japan♭ Faculty of Information Science and Electrical Engineering, Kyushu University
744 Motooka, Nishi-ku, Fukuoka 819–0395 Japan♯ Institute of Systems, Information Technologies and Nanotechnologies
2–1–22 Momochihama, Sawara-ku, Fukuoka 814–0001 Japan
Email: takahashi@cs.okayama-u.ac.jp
Abstract—Multiplicative updates are widely used computational methods for nonnegative matrix factorization(NMF). However, the global convergence of the original
updates is not theoretically guaranteed. By the global convergence, we mean that the sequence of solutions contains
at least one convergent subsequence and the limit of any
convergent subsequence is a stationary point of the NMF
optimization problem. In this paper, we consider a modified multiplicative update for a general error function andgive a suﬃcient condition for the global convergence.

proved that this modification makes the updates of Lee andSeung globally convergent [5, 6].
In this paper, we generalize the global convergence analysis presented in Reference [6]. More specifically, we consider a general error function and show that under someconditions on the auxiliary function the modified multiplicative update is globally convergent. We also show thateight among eleven updates presented in Reference [3] satisfy those conditions.
2. NMF and Multiplicative Updates

1. Introduction
Nonnegative matrix factorization (NMF) is a technique
to decompose a given nonnegative matrix into the productof two nonnegative matrices. Since NMF is useful for finding a set of nonnegative basis vectors for the given nonnegative data, it has attracted great attention from researchers
in machine learning, signal processing, pattern classification, statistics, and so on.
Multiplicative updates proposed by Lee and Seung [1, 2]
are widely used as an eﬃcient computational approach to
NMF. They derived two kinds of multiplicative updates
based on Euclidean distance and I-divergence by taking thefollowing two steps. The first step is to construct an auxiliary function for the original error function, and the second
one is to find its unique minimum point. Yang and Oja [3]recently generalized this procedure and derived eleven multiplicative updates including the ones of Lee and Seung.
However, all of the eleven multiplicative updates are notwell-defined because each of them contains a rational function and its denominator can become zero. A simple way
to avoid this problem is to apply the modification proposedby Gillis and Glineur [4]. Moreover, it has recently beenThis work was partially supported by KAKENHI 24560076,
24500018 and 23310104, and by Proactive Response Against Cyberattacks Through International Collaborative Exchange (PRACTICE),
Ministry of Internal Aﬀairs and Communications, Japan.

Given a nonnegative matrix X ∈ R+m×n , where R+ denotes the set of nonnegative numbers, let us consider the
problem of finding two nonnegative matrices W ∈ Rm×r+and H ∈ Rr×n
+ such that
X ≈ WH(1)

where r is a positive integer less than min{m, n}. A variety
of techniques to find W and H in (1) are called nonnegative matrix factorization (NMF). Although it is important
in NMF how to set the value of r, we will not consider thisproblem in this paper: we simply assume that the value ofr is given together with X. Also, we assume throughout
this paper that every row and column of X has at least one
nonzero entry.
The problem of finding W and H in (1) is usually formulated as a mathematical programming problem of thefollowing form:
minimize 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. So far, various
error functions such as Euclidean distance and I-divergencehave been used for NMF.

- 44 -

Because (2) is a nonconvex optimization problem, it isdiﬃcult to find its global optimal solution. As an approach to find local optimal solutions, multiplicative updates [1–3, 7] are widely used. For example, the multiplicative update rule developed by Lee and Seung [1, 2] for∑
Euclidean distance error function D(W , H) = i j (Xi j −(W H)i j )2 is given by
Wik(l+1)=

Hk(l+1)j=(X(H (l) )T )ik,(W (l) H (l) (H (l) )T )ik((W (l+1) )T X)k j,
Hk(l)j((W (l+1) )T W (l+1) H (l) )k j

Wik(l)(3)(4)

where l represents the number of iterations and Wik(l) (Hk(l)j ,resp.) is the value of Wik (Hk j , resp.) after l updates.
One can easily see from (3) and (4) that if the initialmatrices W (0) and H (0) are positive then W (l) and H (l)are positive for all l ≥ 1. Also, it is known that thesequence {D(W (l) , H (l) )}∞l=1 is monotone decreasing [2],which means that the sequence converges to some constantbecause it is bounded from below. However, this does not
imply that the sequence {(W (l) , H (l) )}∞l=0 converges to a local optimal solution of (2).
The multiplicative update rule given by (3) and (4) is
derived by constructing an auxiliary function for the error∑
function D(W , H) = i j (Xi j − (W H)i j )2 and minimizf , H)f :ing it [2]. In this paper, a function D̄(W , H, W
m×r
r×n
m×r
r×n
R++ × R++ × R++ × R++ → R, where R++ denotes theset of positive numbers, is called an auxiliary function ofD(W , H) if the following two conditions are satisfied:

where ϵ is a small positive constant specified by the user.
This simple modification can be applied to all multiplicative updates. Note that, with this modification, the problem(2) has to be modified asminimize D(W , H)
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. Hibi and Takahashi [5] provedthat the update rule given by (5) and (6) has the global convergence property in the sense of Zangwill [8]. Also, they
have recently shown that this modification can also guarantee the global convergence of the multiplicative update forI-divergence [6].
In the following, we will consider a general error function D(W , H) and give a suﬃcient condition on the auxf , H)
f for the modified update
iliary function D̄(W , H, W
rule to be globally convergent. Before proceeding further,
we need to introduce some notations. The feasible regionof the problem (7) is denoted by Fϵ , that is,Fϵ = {(W , H) | W ≥ ϵ1m×r , H ≥ ϵ1r×n } .
A point (W , H) that satisfies Karush-Kuhn-Tucker (KKT)
conditions for (7):W
H
∂D(W , H)∂Wik∂D(W , H)∂Hk j∂D(W , H)(ϵ − Wik )∂Wik∂D(W , H)(ϵ − Hk j )∂Hk jf > Om×r , H
f > Or×n ,∀W > Om×r , H > Or×n , W
f , H)
f ≥ D(W , H) ,D̄(W , H, W
∀W > Om×r , H > Or×n ,D̄(W , H, W , H) = D(W , H) .
Note that this approach is not restricted to Euclidean distance error function but can be applied to various error
functions. In fact, Yang and Oja [3] proposed a unified approach to develop multiplicative update rules and applied it
to eleven error functions (see Table 1).
3. Global Convergence of Modified Update Rules
The multiplicative update rule given by (3) and (4) andother update rules shown in Table 1 have a common serious
problem that they are not well-defined. In order to solve
this problem, Gillis and Glineur [4] proposed to modify (3)and (4) as()(X(H (l) )T )ik
Wik(l+1) = max ϵ, Wik(l),(5)(W (l) H (l) (H (l) )T )ik()((W (l+1) )T X)k j(l)
Hk(l+1)
=
max
ϵ,
H
, (6)
j
kj
((W (l+1) )T W (l+1) H (l) )k j

(7)

≥
≥

ϵ1m×r ,
ϵ1r×n ,≥

0,∀i, k ,≥

0,∀k, j ,=

0,∀i, k ,=

0,∀k, j ,

is called a stationary point of (7). The set of stationary
points is denoted by S ϵ .
The following assumptions are also needed for later discussions.
m×rf ∈ R++
f ∈ Rr×n
Assumption 1 For any W
and H
++ ,
f
f
D̄(W , H, W , H) is diﬀerentiable with respect to Wik , andsatisfies

f , H)f∂D̄(W , H, W
∂Wik

=
g ,H)
f
(W ,H)=(W∂D(W , H)
.
∂Wik
g ,H)
f
(W ,H)=(WfW
f , H)
f can be expressed asAssumption 2 D̄(W , H,
∑
f
f
f f
ik uik (Wik , W , H) and uik (Wik , W , H) is strictly convex
f , H)
f ∈with respect to Wik on R++ . Also, for each (Wm×r
r×n
R++ × R++ , the problem- 45 -

f , H)fminimize uik (Wik , W
subject to Wik > 0

Table 1: Error functions and multiplicative update rules [3]. Z = (Zi j ) is defined by Zi j = Xi j /(W H)i j .
Error function

Multiplicative updates for Wik

Euclidean distance(XH )ik
Wiknew = Wik (W
HH T )

I-divergence)∑
Wiknew = Wik (ZH
j Hk j(∑)(ln Z )H
Wiknew = Wik exp j ∑ j Hi jk j k j√∑−2
j Xi j (W H)i j Hk j∑
Wiknew = Wik−1
j (W H)i j Hk j( ∑ Z α H ) α1
j ij kj, Wik ∑ j Hk j
new(∑)
Wik = 

(ln
Z
)H

ij kj
j


∑
,
 Wik exp
H
j kj

T

ik

T

Dual I-divergence
Itakura-Saito divergence
α-divergence(∑
β-divergence

Wiknew = Wik

Log-Quad cost

Wiknew = Wik

αβ-Bregman divergence

Wiknew = Wik

Kullback-Leibler divergence

Wiknew = W

γ-divergence(

1−βα(α−1)

∑

(∑

g ,H)
f
(W ,H)=(W(α ≥ 1, 0 < β < 1)

)η
{ 1
∑
1+γ
, γ>0
(W H)ab1+γ· ∑ abX (W H)
, η=
γ
1
ab ab
, γ<0
ab
1−γ
{ 1
( ∑ Zr H
)η
∑
,
r
>
1
k
j
j
(W
H)
i
j
ab
r
Wiknew = Wik ∑ j Hk j · ∑ ab
, η=
r
1−r
ab Xab (W H)ab1, 0 < r < 1(8)

The next theorem is the main result of this paper.
Theorem 1 Let ϵ be any positive number. If an auxiliaryf , H)
f of D(W , H) satisfies Assumpfunction D̄(W , H, W
tions 1–4 then, for any initial matrices W (0) ≥ ϵ1m×r andH (0) ≥ ϵ1r×n , the sequence {(W (l) , H (l) )}∞l=0 generated by
the update rule given by∂D(W , H)
.
∂Hk j
g ,H)
f
(W ,H)=(Wf , H)fvik (Hk j , W
Hk j > 0

has a unique optimal solution that can be explicitly exf , H).
f Furthermore, gk j (Wf , H)
fpressed as Hk∗ j = gk j (Wis continuous and, for each ϵ > 0, there exist dk j > 0 andµk j < 1 such thatf ≥ ϵ1m×r , H
f ≥ ϵ1r×n ,∀W1) α−β+1(W H)

f , H, W
f , H)
f can be expressed asAssumption 4 D̄(W
∑
f
f
f f
k j vk j (Hk j , W , H) and vk j (Hk j , W , H) is strictly convex
f , H)
f ∈with respect to Hk j on R++ . Also, for each (W
r×n
Rm×r
×
R
,
the
problem++++minimizesubject to

β>1
0<β≤1
β<0

γ−1
j Xi j (W H)i j Hk j∑
γj (W H)i j Hk j

f , H)f ≤ cik We νik .
fik (W
ik

=

,

∑

r×n
f ∈ Rm×rf
Assumption 3 For any W++ and H ∈ R++ ,f , H)
f is diﬀerentiable with respect to Hk j , andD̄(W , H, W
satisfies

f , H)f∂D̄(W , H, W
∂Hk j

,

(ZH T +2XH T )ik∑
Tj Hk j +2(W HH )ik

has a unique optimal solution that can be explicitly exf , H).
f Furthermore, fik (Wf , H)
f ispressed as Wik∗ = fik (Wcontinuous and, for each ϵ > 0, there exist cik > 0 andνik < 1 such thatf ≥ ϵ1m×r , H
f ≥ ϵ1r×n ,
∀W

)η

α=0
 1

,



 β
1,
η=



 1

∑
β−2
α−2
j Xi j (W H)i j Hk j +β(1−β) j Xi j (W H)i j Hk j
∑
∑
β−1
α−1
α(α−1) j (W H)i j Hk j +β(1−β) j (W H)i j Hk j
(ZH T )ikab
ik ∑ Habj kj

Wiknew = WikRényi divergence√

β−1
j Xi j (W H)i j Hk j
∑
β
(W
H)
j
i j Hk j

α,0f , H)f ≤ dk j Heµk j . (9)g k j (Wkj

Wik(l+1)= max(ϵ, fik (W (l) , H (l) ))Hk(l+1)j= max(ϵ, gk j (W (l+1) , H (l) ))has at least one convergent subsequence and the limit ofany convergent subsequence is a stationary point of the
problem (7).
Proof Proof sketch will be given in the next section.

□

Let us now consider eleven update rules shown in Table 1. By analyzing the unified method [3] for developing
multiplicative update rules, we can prove that the auxiliary
functions for the first eight error functions (Euclidean distance, I-divergence, Dual I-divergence, Itakura-Saito divergence, α-divergence, β-divergence, Log-Quad cost and αβBregman divergence) satisfy Assumptions 1–4. Therefore,
we can conclude that the modified update rules corresponding to the first eight update rules in Table 1 have a global- 46 -

convergence property. On the other hand, the global convergence of the last three update rules cannot be proved by
Theorem 1 because the inequalities (8) and (9) are not satisfied [9]. This issue will be discussed in [10].
4. Proof Sketch of Theorem 1
Let (5) and (6) be expressed as W (l+1) = A1 (W (l) , H (l) )and H (l+1) = A2 (W (l+1) , H (l) ), respectively. Moreover,
let the update from (W (l) , H (l) ) to (W (l+1) , H (l+1) ) be expressed as A(W (l) , H (l) ) = (W (l+1) , H (l+1) ). Then we have==A(W , H )(A1 (W (l) , H (l) ), A2 (A1 (W (l) , H (l) ), H (l) )) .
(l)

We prove Theorem 1 by using Zangwill’s global convergence theorem [8]. To do so, we need to show that thefollowing statements hold true.
1. (Boundedness) For each initial value (W (0) , H (0) ) ∈Fϵ , the sequence {(W (l) , H (l) )}∞l=0 generated by the update rule A belongs to a compact subset of Fϵ .
2. (Monotoneness) There exists a function h : Fϵ → R
such that(W , H) < S ϵ ⇒ h(A(W , H)) < h(W , H) ,(W , H) ∈ S ϵ ⇒ h(A(W , H)) ≤ h(W , H) .
3. (Continuity) A is continuous in Fϵ \ S ϵ .
Among these three statements, the boundedness of A
can be directly proved by Assumptions 2 and 4 and [9,Lemma 1]. The continuity of A is also apparent from Assumptions 2 and 4. So it suﬃces for us to show the monotoneness of A. This can be done by using Assumptions 1–4
and the following lemmas. However, we omit the details
due to lack of space.
f , H)
f be any point in Fϵ . Then
Lemma 1 Let (WfWf , H)fminimize u(W ) = D̄(W , H,
subject to W ≥ ϵ1m×r

f , H, W
f , H)
fv(H) = D̄(WH ≥ ϵ1r×n

References[2] D. D. Lee and H. S. Seung, “Algorithms for nonnegative matrix factorization,” in Advances in Neural Information Processing Systems, 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 Trans. Neural Networks, vol. 22, pp. 1878–1891, December 2011.
[4] N. Gillis and F. Glineur, “Nonnegative factorizationand the maximum edge biclique problem,” arXiv eprints, October 2008.
[5] R. Hibi and N. Takahashi, “A modified multiplicative
update algorithm for Euclidean distance-based nonnegative matrix factorization and its global convergence,” in Proc. ICONIP2011, Part-II, pp. 655–662,
November 2011.
[6] N. Takahashi and R. Hibi, “Global convergence ofmodified multiplicative updates for nonnegative matrix factorization,” Computational Optimization andApplications, vol. 57, pp. 417–440, March 2014.
[7] L. Finesso and P. Spreij, “Nonnegative matrix factorization and I-divergence alternating minimization,” Linear Algebra and its Applications, vol. 416,
pp. 270–287, 2006.

(10)[8] W. I. Zangwill, Nonlinear programming: A unified
approach. Englewood Cliﬀs, NJ: Prentice-Hall, 1969.

f , H).
f Also,has a unique optimal solution W ∗ = A1 (Wminimizesubject to

We have given a suﬃcient condition for the modified
multiplicative update for NMF to be globally convergent.
Because the suﬃcient condition is very general, the resultof this paper can be applied to many multiplicative updates.
In fact, the global convergence of the modified versionsof eight among the eleven multiplicative updates presented
in [3] is guaranteed by Theorem 1.

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

(W (l+1) , H (l+1) )(l)

5. Conclusion(11)f , H).
f
has a unique optimal solution H ∗ = A2 (W
Lemma 2 The necessary and suﬃcient condition for
f , H)
f to be a stationary point of (7) is that W
f is the(Wf is the unique optiunique optimal solution of (10) and H
mal solution of (11).

[9] J. Katayama, N. Takahashi, and J. Takeuchi, “Boundedness of modified multiplicative updates for nonnegative matrix factorization,” in Proc. CAMSAP2013,
pp. 252–255, December 2013.
[10] M. Seki and N. Takahashi, “New updates based on
Kullback-Leibler, gamma, and Rényi divergences fornonnegative matrix factorization,” to appear in Proc.
NOLTA2014.

- 47 -