Global Convergence of Modified Multiplicative
Updates for Nonnegative Matrix Factorization
Author(s): Norikazu Takahashi and Ryota Hibi
Journal: Computational Optimization and Applications
Volume: 57
Pages: 417–440
Month: March
Year: 2014
DOI: 10.1007/s10589-013-9593-0
Published Version: http://link.springer.com/article/10.1007%2Fs10589-013-9593-0

Computational Optimization and Applications manuscript No.
(will be inserted by the editor)

Global convergence of modified multiplicative updates fornonnegative matrix factorization
Norikazu Takahashi · Ryota Hibi

Received: date / Accepted: date

Abstract Nonnegative matrix factorization (NMF) is the problem of approximating
a given nonnegative matrix by the product of two nonnegative matrices. The multiplicative updates proposed by Lee and Seung are widely used as efficient computational methods for NMF. However, the global convergence of these updates is notformally guaranteed because they are not defined for all pairs of nonnegative matrices. In this paper, we consider slightly modified versions of the original multiplicative
updates and study their global convergence properties. The only difference between
the modified updates and the original ones is that the former do not allow variables
to take values less than a user-specified positive constant. Using Zangwill’s global
convergence theorem, we prove that any sequence of solutions generated by either of
those modified updates has at least one convergent subsequence and the limit of any
convergent subsequence is a stationary point of the corresponding optimization problem. Furthermore, we propose algorithms based on the modified updates that always
stop within a finite number of iterations.
Keywords Nonnegative matrix factorization · Multiplicative update · Global
convergence · Finite termination
Part of this paper was presented at 18th International Conference on Neural Information Processing,
Shanghai, China, in November 2011 [14].
N. Takahashi
Department of Informatics, 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
Present address: Graduate School of Natural Science and Technology, Okayama University, 3-1-1
Tsushima-naka, Kita-ku, Okayama 700-8530 Japan
E-mail: takahashi@cs.okayama-u.ac.jp
R. Hibi
Department of Informatics, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395 Japan
Present address: Mitsubishi UFJ Morgan Stanley Securities Co., Ltd.

2

Norikazu Takahashi, Ryota Hibi

1 Introduction
Nonnegative matrix factorization (NMF) [16,17] is the problem of approximating agiven large nonnegative matrix V by the product W H of two flat nonnegative matrices W and H. If we consider the columns of V as data vectors, the columns of
W and those of H are interpreted as a set of nonnegative basis vectors and a set of
nonnegative coefficient vectors, respectively. Each data vector is thus reproduced approximately by a linear combination of the basis vectors with coefficients stored in
the corresponding column of H. In this sense, NMF can generate a reduced representation of the original data. Moreover, the basis vectors often represent parts of theobject because of the nonnegativity constraints [16]. This is a significant difference
between NMF and other factorization methods such as principal component analysis.
So far, NMF has been successfully applied to various problems in machine learning,
signal processing and so on [2, 5,7,15,16,20,22, 25].
Usually, NMF is formulated as a constrained optimization problem in which theapproximation error has to be minimized with respect to W and H subject to thenonnegativity of these matrices. Lee and Seung [17] considered the cases where
the approximation error is measured by the Euclidean distance and the I-divergence,and proposed iterative methods called the multiplicative updates. These updates arewidely used as simple and efficient computational methods for NMF because of thefollowing three advantages. First, the updates do not contain parameters like the step
size in gradient decent methods, and therefore parameter tuning is not needed. Second, nonnegativity of the matrices W k and H k , the solution after k iterations, is automatically satisfied if the initial matrices W 0 and H 0 are chosen to be positive. Third,
implementation is easy because the update formulae are very simple.
However, the multiplicative updates of Lee and Seung have a serious drawbackthat their global convergence is not guaranteed theoretically. By global convergence,
we mean that, for any initial solution, the sequence of solutions contains at least one
convergent subsequence and the limit of any convergent subsequence is a stationary point of the corresponding optimization problem. The main difficulty in proving
global convergence is that the updates, which are expressed in the form of a fraction,
are not defined for all pairs of nonnegative matrices. Hence the convergence analysis of the multiplicative updates and their variants is an important research issue in
NMF, and many authors have addressed this problem so far [1, 10, 12, 18]. Finesso
and Spreij [10] studied convergence properties of the multiplicative update based on
the I-divergence minimization and proved, under the assumption that W k is normalized after each update so that its Frobenius norm becomes one, that the sequencesof W k and W k H k always converge. However, their result does not guarantee convergence of the sequence of H k . Lin [18] considered the case of the Euclidean distanceminimization and showed that some modifications to the original multiplicative update can make it well-defined and globally convergent. However, since Lin’s modified
update is not multiplicative but additive in some cases, this result cannot be directly
applied to the original update. Recently, Badeau et al. [1] studied local stability of ageneralized multiplicative update, which includes the multiplicative updates of Lee
and Seung as special cases, using Lyapunov’s stability theory and showed that theGlobal convergence of modified multiplicative updates for nonnegative matrix factorization3

local optimal solution of the corresponding optimization problem is asymptoticallystable if one of two matrices W k and H k is fixed for all k.
The objective of this paper is to show that a slight modification can guarantee
global convergence of the multiplicative updates of Lee and Seung [17]. Our attention is focused on the modification proposed by Gillis and Glineur [12]. Their
update, which is a modified version of the Euclidean distance-based multiplicative
update of Lee and Seung [17], returns a user-specified positive constant if the original update returns a value less than the constant. Note that unlike the updates of
Lin [18] and Finesso and Spreij [10], normalization procedure is not involved. Gillis
and Glineur proved that their modified multiplicative update decreases the objectivefunction monotonically and that if a sequence of solutions generated by the update has
a limit point then it is necessarily a stationary point of the corresponding optimization
problem [12]. However, this does not imply global convergence of the update.
In this paper, we consider not only the Euclidean distance-based multiplicative
update but also the I-divergence-based one, and prove that their global convergenceis guaranteed if they are modified as described by Gillis and Glineur [12]. Our proof
is based on Zangwill’s global convergence theorem [28, p.91] which is a fundamental result in optimization theory and has played important roles in the convergence
analysis of many algorithms in machine learning [21,23,26]. We also propose twoalgorithms based on the modified updates. They always stop within a finite numberof iterations after finding an approximate stationary point of the optimization problem.
There are many other approaches that attempt to solve NMF optimization problems. For example, some authors modified the multiplicative updates of Lee and Seung by adding a small positive constant to the denominators so that they are definedfor all nonnegative matrices [3, 18]. Also, some authors proposed to apply different
optimization techniques to NMF optimization problems [3,7, 19]. Furthermore, someauthors derived a variety of multiplicative updates by considering various types ofdivergence between V and W H [1,6,9,27]. Although these updates are potentiallysuperior in some cases, we will not consider them in this paper.
The rest of this paper is organized as follows. In Section 2, we introduce briefly
the NMF optimization problems and the multiplicative updates of Lee and Seung. In
Section 3, the modified multiplicative updates based on the idea of Gillis and Glineur
are first introduced and then convergence theorems for these updates are presented. In
addition, algorithms based on the modified multiplicative updates are proposed andtheir finite termination is proved. In Section 4, the convergence theorems in Section
3 are proved using Zangwill’s global convergence theorem. Finally, in Section 5, weconclude the paper with a brief summary.
Part of this paper (Theorem 1 in Section 3) was presented in the authors’ conference paper [14]. However, no rigorous proof was given there because of space
limitation. In this paper, not only Theorem 1 but also some new results (Theorems 2,
3 and 4) are presented with their complete proofs.

4Norikazu Takahashi, Ryota Hibi

Fig. 1 Nonnegative matrix factorization. A given nonnegative matrix V is approximated by the product of
two nonnegative matrices W and H. The j-th column of V is approximated by the linear combination of
the columns of W where coefficients are the elements of the j-th column of H.

2 Nonnegative matrix factorization and multiplicative updates
Given a nonnegative matrix V ∈ Rn×m
where R+ denotes the set of nonnegative real+
numbers, and a positive integer r ≤ min(n, m), NMF is the problem of finding two
r×m
nonnegative matrices W ∈ Rn×r
such that V is approximately equal to+ and H ∈ R+
W H (see Fig.1). Throughout this paper, we assume the following.
Assumption 1 Each row and column of V has at least one nonzero element.
Let us consider each column of V as a data vector. If the value of r is sufficiently
small, a compact expression for the original data can be obtained through NMF because the total number of elements in the factor matrices W and H is less than that of
the original matrix V . Moreover, the columns of W are regarded as a kind of basis forthe space spanned by the columns of V because each data vector can be approximatedby a linear combination of the columns of W (see Fig. 1).
Lee and Seung [17] employed the Euclidean distance and the I-divergence forthe approximation error between V and W H, and formulated NMF as two types of
optimization problems. In the former case, the problem is expressed asminimize fE (W, H) = ∥V −W H∥2
subject to W ≥ 0, H ≥ 0 ,(1)

where ∥ · ∥ represents the Frobenius norm, that is,nm∥V −W H∥2 = ∑ ∑ (Vi j − (W H)i j )2 ,i=1 j=1

and the inequality W ≥ 0 (resp. H ≥ 0) means that all elements of the matrix W (resp.
H) are nonnegative. In the latter case, the problem is expressed asminimize fD (W, H) = D(V ∥W H)
subject to W ≥ 0, H ≥ 0 ,(2)

where D(·∥·) is defined byn m {D(V ∥W H) = ∑ ∑ Vi j log
i=1 j=1

Vi j−Vi j + (W H)i j(W H)i j}.

Global convergence of modified multiplicative updates for nonnegative matrix factorization5

It is difficult in both cases to find a global optimal solution because the objectivefunctions fE (W, H) and fD (W, H) are not convex. In fact, NP-hardness of NMF wasproved by Vavasis [24]. Therefore, we have to take the second best way, that is, we
try to find a local optimal solution instead of a global one. For this purpose, Lee andSeung [17] proposed the update rulekHak+1
j = Ha j((W k )T V )a j,((W k )T W k H k )a j

Wiak+1 = Wiak

(V (H k+1 )T )ia,k(W H k+1 (H k+1 )T )ia(3)

for the optimization problem (1), and the update rulekHak+1
j = Ha j

Wiak+1 = Wiak∑ni=1 Wiak Vi j /(W k H k )i j,
∑ni=1 Wiakk k+1 )
∑mj=1 Hak+1
ij
j Vi j /(W H∑mj=1 Hak+1
j

(4)
,

for the optimization problem (2), where k represents the iteration count1 . The updates
like (3) and (4) are called the multiplicative updates because the new estimate isgiven by the product of the current estimate and some factor. An advantage of these
multiplicative updates is that, unlike conventional gradient descent methods, there
are no parameters to tune. Another advantage is that positiveness of W k and H k is
guaranteed for all k under Assumption 1 if the initial matrices W 0 and H 0 are chosen
to be positive [19]. For these reasons, the multiplicative updates (3) and (4) are widelyused as simple and effective methods for finding local optimal solutions of (1) and(2).
3 Modified multiplicative updates and their global convergence
The most serious drawback of the multiplicative update rules described by (3) and (4)is that the right-hand sides are not defined for all nonnegative matrices W k and H k (orH k+1 ). For example, in the case of Euclidean distance, we cannot obtain H k+1 by the
update rule (3) when H k = 0, because the denominator of the first equation vanishes.
As mentioned in Section 2, W k and H k are positive for all k if the initial matrices0
W and H 0 are chosen to be positive. Hence the updates can be performed infinitely
many times. However, even though the sequence {(W k , H k )}∞k=0 converges, it is notguaranteed that both limk→∞ W k and limk→∞ H k are positive. This means that the
update rules may not be defined at limk→∞ (W k , H k ), which makes it difficult to provetheir global convergence using known results such as Zangwill’s global convergence
theorem [28, p.91].
In this section, we introduce slightly modified versions of the update rules (3) and(4) which are based on the idea of Gillis and Glineur [12], and present convergence
1 Although it is not explicitly written in their original paper [17] which of H k and H k+1 is used for thecomputation of W k+1 , we consider the latter case throughout this paper as in [18].

6Norikazu Takahashi, Ryota Hibi

theorems. We also propose two algorithms based on the modified updates and provetheir finite termination.
3.1 Euclidean distanceIn order to prevent elements of matrices W k and H k from vanishing, Gillis andGlineur [12] have proposed to modify the update rule (3) as()((W k )T V )a j
k
,
Hak+1
=
max
H
ε
,
aj
j
((W k )T W k H k )a j
(5)()
(V (H k+1 )T )iakk+1
Wia = max Wia k k+1 k+1 T , ε ,(W H (H ) )ia
where ε is any positive constant specified by the user. Each update in (5) returns the
positive constant ε if the corresponding original update in (3) returns a value less thanε . This is the only difference between these two update rules. With the modificationof the update rule from (3) to (5), we have to modify also the optimization problem(1) as follows:
minimize fE (W, H) = ∥V −W H∥2(6)subject to Wia ≥ ε , Ha j ≥ ε , ∀i, a, j .
The feasible region of this optimization problem is denoted by X, that is,
X = {(W, H) |Wia ≥ ε , Ha j ≥ ε , ∀i, a, j} .
Karush-Kuhn-Tucker (KKT) conditions [4] for the problem (6) are expressed as follows:2
Wia ≥ ε ,∀i, a ,

(7)

Ha j ≥ ε , ∀a, j ,(∇W fE (W, H))ia ≥ 0, ∀i, a ,(∇H fE (W, H))a j ≥ 0, ∀a, j ,(∇W fE (W, H))ia (ε −Wia ) = 0,(∇H fE (W, H))a j (ε − Ha j ) = 0,∀i, a ,
∀a, j ,

(8)
(9)
(10)
(11)
(12)

where
∇W fE (W, H) = 2(W H −V )H T ,∇H fE (W, H) = 2W T (W H −V ) .
Therefore, a necessary condition for a point (W, H) to be a local optimal solutionof (6) is that the conditions (7)–(12) are satisfied. Hereafter, we call a point (W, H) a
stationary point of (6) if it satisfies (7)–(12), and denote the set of all stationary pointsof (6) by SE .
The global convergence property of the modified update rule (5) is stated as follows.
2 The conditions (7)–(12) are derived by eliminating Lagrange multipliers in the original KKT conditions.

Global convergence of modified multiplicative updates for nonnegative matrix factorization

7

Theorem 1 Let {(W k , H k )}∞k=0 be any sequence generated by the modified update
rule (5) with the initial point (W 0 , H 0 ) ∈ X. Then (W k , H k ) ∈ X holds for all positive
integers k. Moreover, the sequence has at least one convergent subsequence and thelimit of any convergent subsequence is a stationary point of the optimization problem(6).
Proof of Theorem 1 will be given in the next section.
By making use of Theorem 1, we can immediately construct an algorithm thatterminates within a finite number of iterations. To do so, we relax the conditions(9)–(12) as(∇W fE (W, H))ia ≥ −δ1 , ∀i, a ,(∇H fE (W, H))a j ≥ −δ1 , ∀ j, a ,Wia − ε ≤ δ2 if (∇W fE (W, H))ia > δ1 , ∀i, a ,Ha j − ε ≤ δ2 if (∇H fE (W, H))a j > δ1 , ∀a, j ,

(13)
(14)
(15)
(16)

where δ1 and δ2 are any positive constants specified by the user, and employ theserelaxed conditions as a stopping criterion. Let S̄E be the set of all (W, H) ∈ X satisfying (13)–(16). Then the proposed algorithm is described as follows:
Algorithm 1
n×m
Input: V ∈ R+
, r ∈ N, ε > 0, δ1 > 0, δ2 > 0
Step 1: Choose (W 0 , H 0 ) ∈ X and set k = 0.
Step 2: Find (W k+1 , H k+1 ) by the update rule (5).
Step 3: If (W k+1 , H k+1 ) ∈ S̄E then return W k+1 and H k+1 . Otherwise add 1 to k andgo to Step 2.

Theorem 2 For any positive constants ε , δ1 and δ2 , Algorithm 1 stops within a finite
number of iterations.
kk ∞Proof Let {(W kl , H kl )}∞l=1 be any convergent subsequence of the sequence {(W , H )}k=0
generated by the modified update rule (5), and (W̄ , H̄) ∈ X be the limit of the subsequence. Then, by Theorem 1, (W̄ , H̄) satisfies(∇W fE (W̄ , H̄))ia ≥ 0,(∇H fE (W̄ , H̄))a j ≥ 0,∀i, a ,
∀a, j ,(∇W fE (W̄ , H̄))ia (ε − W̄ia ) = 0,(∇H fE (W̄ , H̄))a j (ε − H̄a j ) = 0,∀i, a ,
∀a, j .

Recall that ∇W fE (W, H) and ∇H fE (W, H) are continuous for all (W, H) ∈ X. For all
1 such that(i, a) such that (∇W fE (W̄ , H̄))ia = 0, there exists a positive integer Lia|(∇W fE (W kl , H kl ))ia | ≤ δ1 ,1∀l ≥ Lia.

1 suchFor all (i, a) such that (∇W fE (W̄ , H̄))ia > 0, there exists a positive integer Lia
thatk1(∇W fE (W kl , H kl ))ia ≥ −δ1 and Wial − ε ≤ δ2 , ∀l ≥ Lia,

8Norikazu Takahashi, Ryota Hibi

because W̄ia = ε holds. For all (a, j) such that (∇H fE (W̄ , H̄))a j = 0, there exists a
positive integer La2 j such that|(∇W fE (W kl , H kl ))ia | ≤ δ1 ,2
∀l ≥ Lia.

For all (a, j) such that (∇H fE (W̄ , H̄))a j > 0, there exists a positive integer La2 j suchthatk(∇H fE (W kl , H kl ))a j ≥ −δ1 and Ha lj − ε ≤ δ2 , ∀l ≥ La2 j ,because H̄a j = ε holds. From these considerations, we immediately see that there
exists a positive integer L such that the stopping criterion of Algorithm 1 is satisfied
for all (W kl , H kl ) with l ≥ L. This means that Algorithm 1 always stops within a finite
number of iteration.
⊓⊔
3.2 I-divergence
As in the case of Euclidean distance, we modify the update rule (4) as()
n
W k V k /(W k H k )i j
k+1
k ∑i=1 ia i j
Ha j = max Ha j,ε ,∑nµ =1 Wµka()
m
H k+1Vi j /(W k H k+1 )i j
k+1
k ∑ j=1 a j
Wia = max Wia
,ε ,
k
∑m
ν =1 Haν(17)

where ε is any positive constant specified by the user. The modified update rule corresponds to modifying the optimization problem (2) as follows:
minimize fD (W, H) = D(V ∥W H)
subject to Wia ≥ ε , Ha j ≥ ε , ∀i, a, j .

(18)

The feasible region of this optimization problem is X as in the case of (6). KKT
conditions for the problem (18) are expressed as follows:3
Wia ≥ ε ,∀i, a ,

(19)

Ha j ≥ ε , ∀a, j ,(∇W fD (W, H))ia ≥ 0, ∀i, a ,(∇H fD (W, H))a j ≥ 0, ∀a, j ,(∇W fD (W, H))ia (ε −Wia ) = 0,(∇H fD (W, H))a j (ε − Ha j ) = 0,
where

m {

(∇W fD (W, H))ia = ∑

j=1

(20)
(21)
(22)

∀i, a ,
∀a, j ,

Vi j Ha j
Ha j −(W H)i j

(23)
(24)

}
,

3 The conditions (19)–(24) are derived by eliminating Lagrange multipliers in the original KKT conditions.

Global convergence of modified multiplicative updates for nonnegative matrix factorizationn {(∇H fD (W, H))a j = ∑i=1

Vi jWia
Wia −(W H)i j

9}.

Therefore, a necessary condition for a point (W, H) to be a local optimal solution of(18) is that the conditions (19)–(24) are satisfied. Hereafter, we call a point (W, H)
a stationary point of (18) if it satisfies (19)–(24), and denote the set of all stationary
points of (18) by SD .
The global convergence property of the modified update rule (17) is stated asfollows.
Theorem 3 Let {(W k , H k )}∞k=0 be any sequence generated by the modified update
rule (17) with the initial point (W 0 , H 0 ) ∈ X. Then (W k , H k ) ∈ X holds for all positive
integers k. Moreover, the sequence has at least one convergent subsequence and thelimit of any convergent subsequence is a stationary point of the optimization problem(18).
Proof of Theorem 3 will be given in the next section.
By making use of Theorem 3, we can easily construct an algorithm that terminates
within a finite number of iterations. To do so, we relax the conditions (21)–(24) as(∇W fD (W, H))ia ≥ −δ1 , ∀i, a ,(∇H fD (W, H))a j ≥ −δ1 , ∀a, j ,Wia − ε ≤ δ2 if (∇W fD (W, H))ia > δ1 , ∀i, a ,

(25)
(26)
(27)

Ha j − ε ≤ δ2 if (∇H fD (W, H))a j > δ1 ,(28)

∀a, j ,

where δ1 and δ2 are any positive constants specified by the user, and employ theserelaxed conditions as a stopping criterion. Let S̄D be the set of all (W, H) ∈ X satisfying (25)–(28). Then the proposed algorithm is described as follows:
Algorithm 2
n×m
Input: V ∈ R+
, r ∈ N, ε > 0, δ1 > 0, δ2 > 0
Step 1: Choose (W 0 , H 0 ) ∈ X and set k = 0.
Step 2: Find (W k+1 , H k+1 ) by the update rule (17).
Step 3: If (W k+1 , H k+1 ) ∈ S̄D then return W k+1 and H k+1 . Otherwise add 1 to k andgo to Step 2.

Theorem 4 For any positive constants ε , δ1 and δ2 , Algorithm 2 stops within a finite
number of iterations.
We omit the proof of Theorem 4 because it is almost same as the proof of Theorem 2.

3.3 Related works
The modified update rule (5) was first proposed by Gillis and Glineur [12], as statedabove. They proved not only that fE (W k , H k ) is nonincreasing under (5) but also that

10

Norikazu Takahashi, Ryota Hibi

if a sequence of solutions generated by (5) has a limit point then it is necessarily astationary point of the optimization problem (6), but these facts are not sufficient toprove global convergence of (5). As a matter of fact, we cannot rule out, for examkkple, the existence of a sequence {(W k , H k )}∞k=0 such that f E (W , H ) takes the samevalue for all k and the sequence visits a finite number of distinct points periodically.
However, on the other hand, in another paper [13], they showed through numericalexperiments that (5) works better than the original update rule (3) in some cases. This
indicates that (5) is important not only from a theoretical point of view but also inpractice.
Lin [18] proposed a modified version of (3) and proved that any sequence {(W k , H k )}∞k=0
generated by the modified rule has at least one convergent subsequence and their limits are stationary points of the optimization problem (1). However, Lin’s update ruleconsiderably differs from the original one because of many extra operations. In particular, in the case where (∇H fE (W k , H k ))a j is negative and Hak j is less than a userspecified small positive constant, Lin’s update rule is not multiplicative but additive.
Also, the matrix W k must be normalized after each update in order to guarantee thatthe sequence {(W k , H k )}∞k=0 is in a bounded set. In contrast, the normalization is
not required in the modified update rule (5). Nevertheless, the boundedness of thesequence {(W k , H k )}∞k=0 generated by (5) is guaranteed as shown in the next section.
Finesso and Spreij [10] studied the convergence properties of the multiplicative
update (4) by interpreting it as an alternating minimization procedure [8]. Under theassumption that the matrix W k is normalized after each update, they proved that anyksequence {(W k , H k )}∞k=0 generated by (4) satisfies the following properties: 1) Wkk
converges to a nonnegative matrix. 2) For each triple (i, a, j), Wia Ha j converges toa nonnegative number. 3) For each pair (a, j), Hak j converges to a nonnegative number if limk→∞ ∑ni=1 Wiak > 0 [10, Theorem 6.1]. However, they said nothing about theconvergence of Hak j for the case where limk→∞ ∑ni=1 Wiak = 0.
Badeau et al. [1] studied the local stability of a generalized multiplicative update,
which includes (3) and (4) as special cases, using Lyapunov’s stability theory andshowed that the local optimal solution of the corresponding optimization problem is
asymptotically stable if one of two factor matrices W and H is fixed to a nonnegative
constant matrix.

4 Proofs of Theorems 1 and 3
We will prove Theorems 1 and 3 in this section. The first parts of these theorems apparently follow from the update rules (5) and (17). In order to prove the second parts,
we make use of Zangwill’s global convergence theorem [28, p.91], which is a fundamental result in optimization theory. Let A be a point-to-point mapping4 from X into
itself and S be a subset of X. Then Zangwill’s global convergence theorem claims thefollowing: if the mapping A satisfies the following three conditions then, for any initial point (W 0 , H 0 ) ∈ X, the sequence {(W k , H k )}∞k=0 generated by A contains at least
4 Although A is assumed to be a point-to-set mapping in the original version of Zangwill’s global
convergence theorem, we consider in this paper its special case where A is a point-to-point mapping.

Global convergence of modified multiplicative updates for nonnegative matrix factorization11

one convergent subsequence and the limit of any convergent subsequence belongs toS.
1. All points in the sequence {(W k , H k )}∞k=0 belong to a compact set in X.
2. There is a function z : X → R satisfying the following two conditions.
(a) If (W, H) ̸∈ S then z(A(W, H)) < z(W, H).
(b) If (W, H) ∈ S then z(A(W, H)) ≤ z(W, H).
3. The mapping A is continuous in X \ S.
In the following, we will first prove Theorem 1 by showing that these conditions
are satisfied when the mapping A is defined by (5) and S is set to SE . We will next
prove Theorem 2 by showing that these conditions are satisfied when the mapping A
is defined by (17) and S is set to SD .
4.1 Proof of Theorem 1
Let us rewrite (5) asH k+1 = A1 (W k , H k ) ,W k+1 = A2 (W k , H k+1 ) ,or, more simply,(W k+1 , H k+1 ) = A(W k , H k ) ,where the mapping A is defined byA(W, H) = (A2 (W, A1 (W, H)), A1 (W, H)) .
Let us also set S = SE . Since the mapping A is continuous in X, the third condition of
Zangwill’s global convergence theorem is satisfied. We will thus show in the following that A also satisfies the remaining two conditions.
The following lemma guarantees that the first condition is satisfied.
Lemma 1 For any initial point (W 0 , H 0 ) ∈ X, the sequence {(W k , H k )}∞k=0 generated by the mapping A belongs to a compact set in X.
Proof Let (W, H) be any point in X. Then we have
Ha j

(W T V )a j
∑n WiaVi j= Ha j r i=1 T
T(W W H)a j∑l=1 (W W )al Hl j∑ni=1 WiaVi j= Ha j r∑l=1 (∑ni=1 WiaWil )Hl j∑ni=1 WiaVi j= Ha j n2∑i=1 Wia Ha j + ∑rl=1,l̸=a (∑ni=1 WiaWil )Hl j∑ni=1 WiaVi j= nr∑i=1 Wia2 + ∑l=1,l̸=a (∑ni=1 WiaWil )(Hl j /Ha j )∑n WiaVi j< i=1∑ni=1 Wia2

12

Norikazu Takahashi, Ryota Hibi()√
∑ni=1 Wia / ∑nµ =1 Wµ2a Vi j√=
∑ni=1 Wia2√
∑ni=1 Vi2j
√
≤
,
ε n
from which the inequality√∑ni=1 Vi2j√
(A1 (W, H))a j ≤ max , εε nholds for any pair (a, j). Note that the right-hand side is a constant which depends onneither W nor H. Similarly, we have√T
∑mj=1 Vi2j
(V H )ia√Wia<,(W HH T )ia
ε mfrom which the inequality√∑mj=1 Vi2j√
(A2 (W, H))ia ≤ max , εε mholds for any pair (i, a). Note that the right-hand side is a constant which depends
on neither W nor H. Hence A(W, H) belongs to a compact set in X. This means that,for any initial point (W 0 , H 0 ) ∈ X, the sequence {(W k , H k )}∞k=0 generated by themapping A belongs to a compact set in X.
⊓⊔
The last step is to prove that the second condition of Zangwill’s global convergence theorem is also satisfied. To do this, we first need to introduce two auxiliaryr×m ×
functions for fE . Let (Ŵ , Ĥ) be any point in X. Let the function gŴ
E : [ε , ∞)r×m
[ε , ∞)→ R be defined by
r

m

′
′
Ŵ
′
gŴ
E (H, H ) = f E (Ŵ , H ) + ∑ ∑ gEa j (Ha j , H ) ,a=1 j=1

r×m → R is defined bywhere the function gŴ
Ea j : [ε , ∞) × [ε , ∞)
′
T
′
′
gŴ
Ea j (Ha j , H ) = 2(Ŵ (Ŵ H −V ))a j (Ha j − Ha j ) +(Ŵ T Ŵ H ′ )a j
(Ha j − Ha′ j )2 . (29)Ha′ j

n×r × [ε , ∞)n×r → R be defined by
Similarly, let the function hĤ
E : [ε , ∞)
n

r

′
Ĥ
′
hĤ
E (W,W ) = f E (W, Ĥ) + ∑ ∑ hEia (Wia ,W ) ,i=1 a=1(30)

Global convergence of modified multiplicative updates for nonnegative matrix factorization

13

n×r → R is defined bywhere the function hĤ
Eia : [ε , ∞) × [ε , ∞)
′
′
T
′
hĤ
Eia (Wia ,W ) = 2((W Ĥ −V )Ĥ )ia (Wia −Wia ) +(W ′ Ĥ Ĥ T )ia(Wia −Wia′ )2 .
Wia′

ĤThe functions gŴ
E and hE are essentially the same as the auxiliary functions considered by Lee and Seung [17], though mathematical expressions are slightly different.
r×m × [ε , ∞)r×m
Ĥ
However, note that the domains of gŴ
E and hE are restricted to [ε , ∞)n×r
n×r
and [ε , ∞) × [ε , ∞) , respectively, in the present paper. This is an important difference between our functions and theirs.
In the following, we give five lemmas which are needed to prove that the second condition of Zangwill’s global convergence theorem is satisfied. Although someof them can be immediately obtained from some of the results given by Lee andSeung [17], we will provide proofs for all lemmas in order to make this paper selfcontained.

Lemma 2 For any Ŵ ∈ [ε , ∞)n×r , the function gŴ
E satisfies the following two conditions:gŴ
E (H, H) = f E (Ŵ , H),′gŴ
E (H, H ) ≥ f E (Ŵ , H),∀H ∈ [ε , ∞)r×m ,′∀H, H ∈ [ε , ∞)

(31)

r×m

.

(32)

Also, for any Ĥ ∈ [ε , ∞)r×m , the function hĤ
E satisfies the following two conditions:hĤ
E (W,W ) = f E (W, Ĥ),∀W ∈ [ε , ∞)n×r ,′hĤE (W,W ) ≥ f E (W, Ĥ),

∀W,W ′ ∈ [ε , ∞)n×r .

Proof We prove only the first part because the second one can be proved in the same
r×m and indices a and j, the firstway. Since gŴ
Ea j (Ha j , H) = 0 holds for all H ∈ [ε , ∞)condition (31) is satisfied. To see that the second condition (32) is also satisfied, we
first rewrite fE (Ŵ , H) using the Taylor series expansion asmrfE (Ŵ , H) = fE (Ŵ , H ′ ) + ∑ ∑ 2(Ŵ T (Ŵ H ′ −V ))a j (Ha j − Ha′ j )a=1 j=1rrm+ ∑ ∑ ∑ (Ŵ T Ŵ )ab (Ha j − Ha′ j )(Hb j − Hb′ j ) .
a=1 b=1 j=1

Then we have′gŴ
E (H, H ) − f E (Ŵ , H) =

m

r

r

∑∑∑

j=1 a=1 b=1

where

(
( j)
Mab

Ha j − Ha′ j
Ha′ j

)(

Hb j − Hb′ j
Hb′ j

Mab = δab (Ŵ T Ŵ H ′ )a j Hb′ j − (Ŵ T Ŵ )ab Ha′ j Hb′ j ,
( j)

)
,

(33)

(34)

and δab represents the Kronecker’s delta. We next show that the matrices M ( j) =( j)
[Mab ] ( j = 1, 2, . . . , m) are positive semi-definite for all Ŵ ∈ [ε , ∞)n×r and H ′ ∈14Norikazu Takahashi, Ryota Hibi

[ε , ∞)r×m . If this is true, the right-hand side of (33) is nonnegative for all H, H ′ ∈[ε , ∞)r×m . Since the right-hand side of (34) can be rewritten asr

Mab = δab ∑ (Ŵ T Ŵ )al Hl′ j Hb′ j − (Ŵ T Ŵ )ab Ha′ j Hb′ j
( j)

{
=

l=1
∑rl=1,l̸=a (Ŵ T Ŵ )al Hl′ j Ha′ j , if a = b−(Ŵ T Ŵ )ab Ha′ j Hb′ j ,if a ̸= b ,the matrix M ( j) satisfies
( j)

r

Maa =

∑ |Mal |,( j)

a = 1, 2, . . . , r ,l=1,l̸=awhich means that M ( j) is real, symmetric and diagonally dominant with positive diagonal elements. Therefore, M ( j) ( j = 1, 2, . . . , m) are positive semi-definite.
⊓⊔
Lemma 3 Let (Ŵ , Ĥ) be any point in X. Then gŴ
E (H, Ĥ), which is considered as ar×mfunction of H, is strictly convex in [ε , ∞) . Also, hĤ
E (W, Ŵ ), which is considered asa function of W , is strictly convex in [ε , ∞)n×r .
Proof The second-order partial derivatives of gŴ
E (H, Ĥ) are given by∂ 2 gŴ
E (H, Ĥ)=∂ Ha j ∂ Ha′ j′{(Ŵ T Ŵ Ĥ)a j, if (a, j) = (a′ , j′ )
Ĥa j

0,

otherwise ,where (a, j), (a′ , j′ ) ∈ {1, 2, . . . , r}×{1, 2, . . . , m}. Since (Ŵ T Ŵ Ĥ)a j /Ĥa j is a positive
r×m . The second part can be proved inconstant, gŴ
E (H, Ĥ) is strictly convex in [ε , ∞)the same way.
⊓⊔
Lemma 4 Let (Ŵ , Ĥ) be any point in X. The optimization problem
minimize gŴ
E (H, Ĥ)subject to Ha j ≥ ε , ∀a, j(35)

has a unique optimal solution which is given by A1 (Ŵ , Ĥ). Also, the optimization
problemminimize hĤ
E (W, Ŵ )subject to Wia ≥ ε , ∀i, a(36)

has a unique optimal solution which is given by A2 (Ŵ , Ĥ).
Proof It suffices for us to show that for any pair (a, j), the optimization problem
minimize gŴ
Ea j (Ha j , Ĥ)subject to Ha j ≥ ε(37)

Global convergence of modified multiplicative updates for nonnegative matrix factorization

15

has a unique optimal solution which is given by (A1 (Ŵ , Ĥ))a j and that for any pair(i, a), the optimization problem
minimize hĤ
Eia (Wia , Ŵ )subject to Wia ≥ ε(38)

has a unique optimal solution which is given by (A2 (Ŵ , Ĥ))a j . In the following, we
consider only the first part because the second part can be proved similarly. SinceŴgŴ
Ea j (Ha j , Ĥ) is strictly convex in [ε , ∞), the equation dgEa j (Ha j , Ĥ)/dHa j = 0 has at
most one solution in [ε , ∞). By solving this equation, we have
Ha j = Ĥa j(Ŵ T V )a j.
(Ŵ T Ŵ Ĥ)a j

Let the right-side hand be denoted by Ha∗j which is a nonnegative number. If Ha∗j ≥ εthen Ha∗j is apparently the optimal solution of (37). If Ha∗j < ε then ε is the optimal
solution of (37) because gŴ
Ea j (Ha j , Ĥ) is strictly monotone increasing in [ε , ∞). Therefore the optimal solution of (37) is identical with (A1 (Ŵ , Ĥ))a j .
⊓⊔Lemma 5 The inequality fE (A(Ŵ , Ĥ)) ≤ fE (Ŵ , Ĥ) holds for all (Ŵ , Ĥ) ∈ X.
Proof By Lemmas 2 and 4, we haveŴfE (Ŵ , A1 (Ŵ , Ĥ)) ≤ gŴ
E (A1 (Ŵ , Ĥ), Ĥ) ≤ gE (Ĥ, Ĥ) = f E (Ŵ , Ĥ), ∀(Ŵ , Ĥ) ∈ X

andĤ
fE (A2 (Ŵ , Ĥ), Ĥ) ≤ hĤ
E (A2 (Ŵ , Ĥ), Ŵ ) ≤ hE (Ŵ , Ŵ ) = f E (Ŵ , Ĥ), ∀(Ŵ , Ĥ) ∈ X .

From these two inequalities, we havefE (A(Ŵ , Ĥ)) = fE (A2 (Ŵ , A1 (Ŵ , Ĥ)), A1 (Ŵ , Ĥ)) ≤ fE (Ŵ , A1 (Ŵ , Ĥ)) ≤ fE (Ŵ , Ĥ)(39)which completes the proof.
⊓⊔
Lemma 6 (Ŵ , Ĥ) ∈ SE if and only if Ĥ and Ŵ are the optimal solutions of (35) and(36), respectively.
Proof It suffices for us to show that (Ŵ , Ĥ) ∈ SE if and only if Ĥa j is the optimal
solution of (37) for any pair (a, j) and Ŵia is the optimal solution of (38) for any pair(i, a). By the definition (29) of gŴ
Ea j (Ha j , Ĥ), we have∂ gŴEa j (Ha j , Ĥ)∂ Ha j= (∇H fE (Ŵ , Ĥ))a j .
Ha j =Ĥa j

Since gŴEa j (Ha j , Ĥ) is strictly convex in [ε , ∞), the necessary and sufficient condition
for Ĥa j to be the optimal solution of (37) is given by{= 0, if Ĥa j > ε(∇H fE (Ŵ , Ĥ))a j≥ 0, if Ĥa j = ε ,16Norikazu Takahashi, Ryota Hibi

which is equivalent to the set of conditions (8) and (10). By the definition (30) ofhĤ
Eia (Wia , Ŵ ), we have

∂ hĤEia (Wia , Ŵ )∂ Wia

= (∇W fE (Ŵ , Ĥ))ia ,∀i, a .

Wia =Ŵia

Hence the necessary and sufficient condition for Ŵ to be the optimal solution of (36)
is given by{= 0, if Ŵia > ε(∇W fE (Ŵ , Ĥ))ia∀i, a ,
≥ 0, if Ŵia = εwhich is equivalent to the set of conditions (7) and (9).

⊓⊔From Lemmas 4–6, we derive the following lemma which claims that the second
condition of Zangwill’s global convergence theorem is satisfied by setting z = fE .
Lemma 7 Let (Ŵ , Ĥ) be any point in X. If (Ŵ , Ĥ) ∈ SE then A(Ŵ , Ĥ) = (Ŵ , Ĥ).
Otherwise fE (A(Ŵ , Ĥ)) < fE (Ŵ , Ĥ). That is, SE is identical with the set of fixed
points of the mapping A.
Proof We first consider the case where (Ŵ , Ĥ) ∈ SE . By Lemma 6, Ĥ and Ŵ are
unique optimal solutions of (35) and (36), respectively. By Lemma 4, this implies
A1 (Ŵ , Ĥ) = Ĥ and A2 (Ŵ , Ĥ) = Ŵ . Therefore, we haveA(Ŵ , Ĥ) = (A2 (Ŵ , A1 (Ŵ , Ĥ)), A1 (Ŵ , Ĥ)) = (A2 (Ŵ , Ĥ), Ĥ) = (Ŵ , Ĥ) .
We next consider the case where (Ŵ , Ĥ) ̸∈ SE . In this case, by Lemma 6, at least oneof the following statements must be false: 1) Ĥ is the unique optimal solution of (35).
2) Ŵ is the unique optimal solution of (36). Suppose that the statement 1) does notŴhold true. Then, by Lemma 4, we have gŴ
E (A1 (Ŵ , Ĥ), Ĥ) < gE (Ĥ, Ĥ) which implies
that the second inequality of (39) holds as a strict inequality. Therefore, fE (A(Ŵ , Ĥ))is strictly less than fE (Ŵ , Ĥ). Suppose next that the statement 1) holds true but 2) doesĤnot. Then, by Lemma 4, we have A1 (Ŵ , Ĥ) = Ĥ and hĤ
E (A2 (Ŵ , Ĥ), Ŵ ) < hE (Ŵ , Ŵ ).
From these facts and (39), we havefE (A(Ŵ , Ĥ)) = fE (A2 (Ŵ , A1 (Ŵ , Ĥ)), A1 (Ŵ , Ĥ)) = fE (A2 (Ŵ , Ĥ), Ĥ) < fE (Ŵ , Ĥ) .
Therefore, fE (A(Ŵ , Ĥ)) is strictly less than fE (Ŵ , Ĥ).

4.2 Proof of Theorem 3
As in the proof of Theorem 1, let us rewrite (17) asH k+1 = A1 (W k , H k ) ,W k+1 = A2 (W k , H k+1 ) ,⊓⊔Global convergence of modified multiplicative updates for nonnegative matrix factorization

17

or, more simply,(W k+1 , H k+1 ) = A(W k , H k ) ,where the mapping A is defined byA(W, H) = (A2 (W, A1 (W, H)), A1 (W, H)) .
Let us also set S = SD . Since the mapping A is continuous in X, the third condition of Zangwill’s global convergence theorem is satisfied. We will thus show in thefollowing that A also satisfies the remaining two conditions.
The following lemma guarantees that the first condition is satisfied.
Lemma 8 For any initial point (W 0 , H 0 ) ∈ X, the sequence {(W k , H k )}∞k=0 generated by the mapping A belongs to a compact set in X.
Proof Let (W, H) be any point in X. Then we have
Ha jn
WiaVi j∑ni=1 WiaVi j /(W H)i j= Ha j ∑ nrnW∑i=1 Wia∑i=1 µ =1 µ a ∑l=1 Wil Hl jn

WiaVi jnr
i=1 ∑µ =1 Wµ a (Wia Ha j + ∑l=1,l̸=a Wil Hl j )= Ha j ∑n

WiaVi j
r
n
W
{W
+
∑
ia ∑l=1,l̸=a Wil (Hl j /Ha j )}i=1 µ =1 µ a=∑nWiaVi jni=1 (∑µ =1 Wµ a )Wia
n
Vi j
=
n
i=1 ∑µ =1 Wµ a∑ni=1 Vi j

<∑

∑

≤

εn

,

from which the inequality((A1 (W, H))a j ≤ max∑ni=1 Vi j
,ε
εn

)

holds for any pair (a, j). Note that the right-hand side is a constant which depends onneither W nor H. Similarly, we haveWia

∑mj=1 Ha jVi j /(W H)i j∑mj=1 Vi j<,mεm
∑ j=1 Ha jfrom which the inequality(A2 (W, H))ia ≤ max(∑mj=1 Vi j
,ε
εm

)

holds for any pair (i, a). Note that the right-hand side is a constant which depends
on neither W nor H. Hence A(W, H) belongs to a compact set in X. This means that,for any initial point (W 0 , H 0 ) ∈ X, the sequence {(W k , H k )}∞k=0 generated by themapping A belongs to a compact set in X.
⊓
⊔

18

Norikazu Takahashi, Ryota Hibi

The last step is to prove that the second condition of Zangwill’s global convergence theorem is also satisfied. To do this, we first need to introduce two auxilr×m ×iary functions fD . Let (Ŵ , Ĥ) be any point in X. Let the function gŴ
D : [ε , ∞)r×m
→ R be defined by[ε , ∞)}{
n mr
Ŵia Ha′ j
Vi j
′
Ŵ
′
gD (H, H ) = ∑ ∑ Vi j logVi j −Vi j +∑ Ŵia Ha j log (Ŵ H ′ )i j
(Ŵ H ′ )i j a=1
i=1 j=1rm′+ ∑ ∑ gŴ
Da j (Ha j , H ) ,a=1 j=1

r×m is defined bywhere the function gŴ
Da j : [ε , ∞) × [ε , ∞)}nn {ŴiaVi jŴ′′
gDa j (Ha j , H ) = Ha j ∑ Ŵia − Ha j ∑log(Ŵia Ha j ) .
′i=1
i=1 (Ŵ H )i j
n×r × [ε , ∞)n×r → R be defined by
Similarly, let the function gĤ
D : [ε , ∞)
{
}
r
n m
Vi j
Wia′ Ĥa j′
Ĥ
′hD (W,W ) = ∑ ∑ Vi j logVi j −Vi j +∑ Wia Ĥa j log (W ′ Ĥ)i j(W ′ Ĥ)i j a=1
i=1 j=1nr′+ ∑ ∑ hĤDia (Wia ,W ) ,i=1 a=1

n×r → R is defined bywhere the function hĤ
Dia : [ε , ∞) × [ε , ∞)′hĤDia (Wia ,W ) =

m {

∑

j=1

}
Wia′ Ĥa j
Wia Ĥa j −Vi j ′log(Wia Ĥa j ) .
(W Ĥ)i j

ĤThe functions gŴ
D and hD are essentially the same as the auxiliary functions considered by Lee and Seung [17], though mathematical expressions are slightly different. In the following, we give five lemmas which are needed to prove that the second
condition of Zangwill’s global convergence theorem is satisfied. Although Lemmas
9 and 10 below can be immediately obtained from some of the results given by Lee
and Seung [17], we will provide proofs for these lemmas in order to make this paper self-contained. On the other hand, as for Lemmas 11, 12 and 13, we omit proofs
because they are similar to those for Lemmas 4, 5 and 6.

Lemma 9 Let (Ŵ , Ĥ) be any point in X. The function gŴ
D satisfies the following two
conditions:gŴ
D (H, H) = f D (Ŵ , H),′gŴ
D (H, H ) ≥ f D (Ŵ , H),∀H ∈ [ε , ∞)r×m ,′∀H, H ∈ [ε , ∞)

r×m

(40)
.

Also, the function hĤ
D satisfies the following two conditions:hĤ
D (W,W ) = f D (W, Ĥ),∀W ∈ [ε , ∞)n×r ,′hĤ
D (W,W ) ≥ f D (W, Ĥ),

∀W,W ′ ∈ [ε , ∞)n×r .

(41)

Global convergence of modified multiplicative updates for nonnegative matrix factorization

19

Proof We prove only the first part because the second part can be proved in the sameway. For any Ŵ ∈ [ε , ∞)n×r and H ∈ [ε , ∞)r×m , gŴ
D (H, H) can be transformed asgŴ
D (H, H) =

{

r
Vi j
Ŵia Ha j
V
logV
−V
+
∑ ∑ i j i j i j (Ŵ H)i j ∑ Ŵia Ha j log (Ŵ H)i j
i=1 j=1
a=1
{
}
r m
n
n
ŴiaVi j+ ∑ ∑ Ha j ∑ Ŵia − Ha j ∑log(Ŵia Ha j )a=1 j=1i=1i=1 (Ŵ H)i j
n

m

n

m

n

m

n

}

m

= ∑ ∑ (Vi j logVi j −Vi j ) − ∑ ∑ Vi j log(Ŵ H)i j + ∑ ∑ (Ŵ H)i j
i=1 j=1n m {i=1 j=1= ∑ ∑ Vi j log
i=1 j=1

Vi j−Vi j + (Ŵ H)i j(Ŵ H)i j

i=1 j=1}= fD (Ŵ , H) .
Thus the condition (40) holds true. In order to show (41), we consider}Ŵia Ha′ j
Ha j′∑ ∑ Vi j log(Ŵ H)i j − log H ′ − (Ŵ H ′ )i j log(Ŵ H )i j .
aj
i=1 j=1(42)
From the concavity of the log function,′gŴ
D (H, H )− f D (Ŵ , H) =

n

{

m

(
log(Ŵ H)i j = log)r∑ Ŵia Ha j

a=1

r

(

Ŵia Ha j≥ ∑ µa log
µa
a=1

)
(43)

for any H ∈ [ε , ∞)r×m and any set of positive numbers µ1 , µ2 , . . . , µr such that ∑ra=1 µa =1. By substituting µa = (Ŵia Ha′ j )/(Ŵ H ′ )i j for a = 1, 2, . . . , r into (43), we have(log(Ŵ H)i j = logr)∑ Ŵia Ha j

a=1)
(Ŵ H ′ )i j
log Ŵia Ha j ·≥ ∑′Ŵia Ha′ j
a=1 (Ŵ H )i j{}
r Ŵ H ′Ha jia a j′≥ ∑log(Ŵ H )i j + log ′,′Ha j
a=1 (Ŵ H )i jr

Ŵia Ha′ j(which implies that the right-hand side of (42) is nonnegative for all H, H ′ ∈ [ε , ∞)r×m .
⊓⊔Lemma 10 Let (Ŵ , Ĥ) be any point in X. Then gŴ
D (H, Ĥ), which is considered as ar×mfunction of H, is strictly convex in [ε , ∞) . Also, hĤ
D (W, Ŵ ), which is considered asa function of W , is strictly convex in [ε , ∞)n×r .

20

Norikazu Takahashi, Ryota Hibi

Proof The second-order partial derivatives of gŴ
D (H, Ĥ) are given by∂ 2 gŴ
D (H, Ĥ)=∂ Ha j ∂ Ha′ j′

{ Ĥ

aj

Ha2j

0,

ŴiaVi j∑ni=1 (Ŵ Ĥ) , if (a, j) = (a′ , j′ )ijotherwise ,where (a, j), (a′ , j′ ) ∈ {1, 2, . . . , r}×{1, 2, . . . , m}. Note here that (Ĥa j /Ha2j ) ∑ni=1 (ŴiaVi j /(Ŵ Ĥ)i j )is positive for all Ha j ∈ [ε , ∞) because of Assumption 1. Therefore, gŴ
D (H, Ĥ) is
strictly convex in [ε , ∞)r×m . The second part can be proved in the same way.
⊓⊔
Lemma 11 Let (Ŵ , Ĥ) be any point in X. The optimization problem
minimize gŴ
D (H, Ĥ)subject to Ha j ≥ ε , ∀a, j(44)

has a unique optimal solution which is given by A1 (Ŵ , Ĥ). Also, the optimization
problemminimize hĤ
D (W, Ŵ )(45)
subject to Wia ≥ ε , ∀i, a
has a unique optimal solution which is given by A2 (Ŵ , Ĥ).
Lemma 12 The inequality fD (A(Ŵ , Ĥ)) ≤ fD (Ŵ , Ĥ) holds for all (Ŵ , Ĥ) ∈ X.
Lemma 13 (Ŵ , Ĥ) ∈ SD if and only if Ĥ and Ŵ are the optimal solutions of (44) and(45), respectively.
From Lemmas 11–13, we derive the following lemma which claims that the second condition of Zangwill’s global convergence theorem is satisfied by setting z = fD .
The proof is omitted because it is similar to that for Lemma 7.
Lemma 14 Let (Ŵ , Ĥ) be any point in X. If (Ŵ , Ĥ) ∈ SD then A(Ŵ , Ĥ) = (Ŵ , Ĥ).
Otherwise fD (A(Ŵ , Ĥ)) < fD (Ŵ , Ĥ). That is, SD is identical with the set of fixed
points of the mapping A.

5 Conclusion
We have shown that the global convergence of the multiplicative updates proposed by
Lee and Seung is established if they are slightly modified as discussed by Gillis andGlineur. Their idea is just to prevent each variable from becoming smaller than a userspecified positive constant ε , but this slight modification guarantees the boundednessof solutions without normalization. Using Zangwill’s global convergence theorem,
we have proved that any sequence of solutions generated by the modified updates has
at least one convergent subsequence and the limit of any convergent subsequence isa stationary point of the corresponding optimization problem. Furthermore, we havedeveloped two algorithms based on the modified updates which always stop within afinite number of iterations after finding an approximate stationary point.

Global convergence of modified multiplicative updates for nonnegative matrix factorization

21

One may be concerned with the fact that matrices obtained by the modified updates are always dense. However, when sparse matrices are preferable, we only haveto replace all ε in the obtained matrices with zero. If ε is set to a small positive number, this replacement will not affect the results significantly. It is in fact proved thatsetting the entries of W and H equal to ε to zero gives a solution which is O(ε ) close
to a stationary point of the original problem, and that the objective function is affectedby an additive factor of at most O(ε ) [11].
The approach presented in this paper may be applied to various multiplicativealgorithms for NMF or other optimization problems. Developing a unified framework
for the global convergence analysis of multiplicative updates is a topic for future
research.
Acknowledgements This work was partially supported by JSPS KAKENHI Grant Numbers 24560076
and 23310104, and by the project “R&D for cyber-attack predictions and rapid response technology bymeans of international cooperation” of the Ministry of Internal Affairs and Communications, Japan.

References
1. Badeau, R., Bertin, N., Vincent, E.: Stability analysis of multiplicative update algorithms and application to nonnegative matrix factorization. IEEE Transactions on Neural Networks 21(12), 1869–1881(2010)
2. Berry, M.W., Browne, M.: Email surveillance using non-negative matrix factorization. Computationaland Mathematical Organization Theory 11, 249–264 (2005)
3. Berry, M.W., Browne, M., Langville, A.N., Pauca, V.P., Plemmons, R.J.: Algorithms and applications
for approximate nonnegative matrix factorization. Computational Statistics and Data Analysis 52,
155–173 (2007)
4. Boyd, S., Vandenberghe, L.: Convex Optimization. Cambridge University Press, Cambridge, UK(2004)
5. Brunet, J.P., Tamayo, P., Golub, T.R., Mesirov, J.P.: Metagenes and molecular pattern discovery using
matrix factorization. Proceedings of National Academy of Science 101(12), 4164–4169 (2004)
6. Cai, D., He, X., Han, J., Huang, T.S.: Graph regularized nonnegative matrix factorization for data
representation. IEEE Transactions on Pattern Analysis and Machine Intelligence 33(8), 1548–1560(2011)
7. Cichocki, A., Zdunek, R., Phan, A.H., Amari, S.I.: Nonnegative Matrix and Tensor Factorizations.
John Wiley & Sons, West Sussex, U.K. (2009)
8. Csiszár, I., Tusnády, G.: Information geometry and alternating minimization procedures. Statistics
and Decisions, Supplemental Issue pp. 205–237 (1984)
9. Févotte, C., Bertin, N., Durrieu, J.L.: Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis. Neural Computation 21(3), 793–830 (2009)
10. Finesso, L., Spreij, P.: Nonnegative matrix factorization and I-divergence alternating minimization.
Linear Algebra and its Applications 416, 270–287 (2006)
11. Gillis, N.: Nonnegative matrix factorization: Complexity, algorithms and applications. Ph.D. thesis,
Université Catholique de Louvain, Louvain-la-Neuve (2011)
12. Gillis, N., Glineur, F.: Nonnegative factorization and the maximum edge biclique problem. ArXiv
e-prints (2008)
13. Gillis, N., Glineur, F.: Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization. Neural Computation 24(4), 1085–1105 (2012)
14. Hibi, R., Takahashi, N.: A modified multiplicative update algorithm for Euclidean distance-based
nonnegative matrix factorization and its global convergence. In: Proceedings of 18th International
Conference on Neural Information Processing, Part-II, pp. 655–662 (2011)15. Holzapfel, A., Stylianou, Y.: Musical genre classification using nonnegative matrix factorizationbased features. IEEE Transactions on Audio, Speech, and Language Processing 16(2), 424–434(2008)

22

Norikazu Takahashi, Ryota Hibi

16. Lee, D.D., Seung, H.S.: Learning the parts of objects by non-negative matrix factorization. Nature
401, 788–792 (1999)
17. Lee, D.D., Seung, H.S.: Algorithms for non-negative matrix factorization. In: T.K. Leen, T.G. Dietterich, V. Tresp (eds.) Advances in Neural Information Processing Systems, vol. 13, pp. 556–562(2001)
18. Lin, C.J.: On the convergence of multiplicative update algorithms for nonnegative matrix factorization.
IEEE Transactions on Neural Networks 18(6), 1589–1596 (2007)
19. Lin, C.J.: Projected gradient methods for non-negative matrix factorization. Neural Computation
19(10), 2756–2779 (2007)
20. Lu, W., Sun, W., Lu, H.: Robust watermarking based on DWT and nonnegative matrix factorization.
Computers and Electrical Engineering 35, 183–188 (2009)
21. Sha, F., Min, Y., Saul, L.K., Lee, D.D.: Multiplicative updates for nonnegative quadratic programming. Neural Computation 19, 2004–2031 (2007)
22. Shahnaz, F., Berry, M.W., Pauca, V.P., Plemmons, R.J.: Document clustering using nonnegative matrix
factorization. Information Processing and Management 42, 373–386 (2006)
23. Takahashi, N., Nishi, T.: Global convergence of decomposition learning methods for support vector
machines. IEEE Transactions on Neural Networks 17(6), 1362–1369 (2006)
24. Vavasis, S.A.: On the complexity of nonnegative matrix factorization. SIAM Journal on Optimization
20(3), 1364–1377 (2009)
25. Wang, R.S., Zhang, S., Wang, Y., Zhang, X.S., Chen, L.: Clustering complex networks and biologicalnetworks by nonnegative matrix factorization with various similarity measures. Neurocomputing 72,
134–141 (2008)
26. Wu, C.F.J.: On the convergence properties of the EM algorithm. The Ananls of Statistics 11(1),95–103 (1983)
27. Yang, Z., Oja, E.: Unified development of multiplicative algorithm for linear and quadratic nonnegative matrix factorization. IEEE Transactions on Neural Networks 22(12), 1878–1891 (2011)
28. Zangwill, W.I.: Nonlinear programming: A unified approach. Prentice-Hall, Englewood Cliffs, NJ(1969)

Erratum Sheet
Norikazu Takahashi and Ryota Hibi, “Global convergence of modified multiplicative updates fornonnegative matrix factorization”, Computational Optimization and Applications, vol.57, pp.417–
440, 2014.
1. Eq.(42) on Page 437 is not correct. It should be replaced with
{
}
n ∑
m
r
′
∑
∑
Ŵia Haj
HajŴ′′
Vij log(Ŵ H)ij −gD (H, H ) − fD (Ŵ , H) =log ′ − log(Ŵ H )ij
′)
Haj
(
Ŵ
H
ij
a=1
i=1 j=1
which can be derived as follows:ŴgD(H, H ′ ) − fD (Ŵ , H)
{
n ∑
m
∑
=
Vij log Vij − Vij +i=1 j=1+mr ∑∑a=1 j=1

−

n ∑
m
∑

{
Haj

n
∑

{

=

′
Ŵia − Haj

Vij log{

Vij(Ŵ H)ij+

Vij log Vij − Vij + Vij

i=1 j=1−mn ∑∑i=1 j=1=n ∑m∑
i=1 j=1n ∑m∑{Ŵia Haj − Vij

a=1

Vij log{

Vij

1(Ŵ H)ij

log(Ŵ H)ij +{′Ŵia Hajlogn∑Ŵia Vij′
Ŵia Haj}(Ŵ H ′ )ij}

log(Ŵia Haj )′i=1 (Ŵ H )ij}− Vij + (Ŵ H)ij

i=1 j=1

{ r
n ∑
m
∑
∑

r
∑

(Ŵ H ′ )ij a=1

i=1

i=1 j=1
m
n ∑
∑

Vij

r
′
∑
Ŵia Haj

log′
a=1 (Ŵ H )ijr′∑Ŵia Haj′
a=1 (Ŵ H )ij′
Ŵia Haj}(Ŵ H ′ )ij}

log(Ŵia Haj )}− Vij + (Ŵ H)ijr′∑
Ŵia Haj′
a=1 (Ŵ H )ij}(Ŵ H ′ )ij Ŵia Haj′
Haj ∑ Ŵia Haj=
Vij log(Ŵ H)ij −log ′ −
log(Ŵ H ′ )ij
′)
′)
H
(
Ŵ
H
(
Ŵ
H
aj
ij
ij
a=1
a=1
i=1 j=1
}
{
n
m
r
′
∑∑
∑ Ŵia Haj
Haj=
Vij log(Ŵ H)ij −log ′ − log(Ŵ H ′ )ij .
Haj(Ŵ H ′ )ij
i=1 j=1r′∑
Ŵia Hajlog′
Ŵia Haj

r

}

a=1

(Last updated: September 3, 2014)