A Unified Global Convergence Analysis ofMultiplicative Update Rules for Nonnegative Matrix
Factorization
Author(s): Norikazu Takahashi, Jiro Katayama, Masato Seki and Jun’ichi
Takeuchi
Journal: Computational Optimization and Applications
Volume: 71
Number: 1
Pages: 221–250
Month: September
Year: 2018
DOI: 10.1007/s10589-018-9997-y

This is a post-peer-review, pre-copyedit version of an article published in Computational
Optimization and Applications. The final authenticated version is available online at:
http://dx.doi.org/10.1007/s10589-018-9997-y.

Noname manuscript No.
(will be inserted by the editor)

A unified global convergence analysis ofmultiplicative update rules for nonnegative matrix
factorization
Norikazu Takahashi · Jiro Katayama ·Masato Seki · Jun’ichi Takeuchi

Received: date / Accepted: date

Abstract Multiplicative update rules are a well-known computational method
for nonnegative matrix factorization. Depending on the error measure between
two matrices, various types of multiplicative update rules have been proposedso far. However, their convergence properties are not fully understood. This
paper provides a suﬃcient condition for a general multiplicative update ruleto have the global convergence property in the sense that any sequence ofsolutions has at least one convergent subsequence and the limit of any convergent subsequence is a stationary point of the optimization problem. Using thiscondition, it is proved that many of the existing multiplicative update rules
have the global convergence property if they are modified slightly so that all
variables take positive values. This paper also proposes new multiplicative update rules based on Kullback-Leibler, Gamma, and Rényi divergences. It isshown that these three rules have the global convergence property if the same
modification as above is made.
Part of this paper was presented in 2014 International Symposium on Nonlinear Theory andits Applications [34, 37].
N. Takahashi
Graduate School of Natural Science and Technology, Okayama University, 3–1–1 Tsushimanaka, Kita-ku, Okayama 700–8530, Japan
Tel.: +81–86–251–8179
Fax: +81–86–251–8256
E-mail: takahashi@cs.okayama-u.ac.jp
J. Katayama
Department of Informatics, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819–0395,
Japan
M. Seki
Graduate School of Natural Science and Technology, Okayama University, 3–1–1 Tsushimanaka, Kita-ku, Okayama 700–8530, Japan
J. Takeuchi
Department of Informatics, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819–0395,
Japan

2

Norikazu Takahashi et al.

Keywords Nonnegative matrix factorization · Multiplicative update rule ·Global convergence
Mathematics Subject Classification (2010) 90C30 · 90C90 · 68W40 ·15A23
1 Introduction
Nonnegative matrix factorization (NMF), which was first introduced by Paatero
and Tapper [32], is a mathematical operation that decomposes a given nonnegative matrix into two nonnegative low-rank factor matrices. Since the seminal
papers by Lee and Seung [28,29], NMF has attracted a lot of of attention fromresearchers in diﬀerent disciplines, and has been shown to be useful for many
applications in pattern recognition [18], text mining [3], document clustering [35, 43], signal processing [10], music analysis [20, 33], graph analysis [40],cyber security [44], and so on.
NMF is formulated as an optimization problem in which an error between
the given nonnegative matrix and the product of two factor matrices is minimized subject to the nonnegativity constraints on the factor matrices. Unlike
nonnegative rank factorization [2, 4, 5] which has a longer history, the ranksof the factor matrices in NMF are allowed to be less than the maximum.
However, because the objective function is nonconvex with respect to the two
factor matrices, finding a global optimal solution is very diﬃcult in general.
In fact, it has been proved that the NMF optimization problem is NP-hardwhen the error is measured by Euclidean distance [39]. Therefore the goal insolving an NMF optimization problem is to find a local optimal solution.
During the last two decades, many algorithms for NMF have been developed [6–9, 11–13, 15, 17, 19, 21, 23–25, 27–29, 31, 45]. The most well-known andwidely used algorithms are the two multiplicative update rules proposed by
Lee and Seung [28, 29]: one is for NMF with Euclidean distance and the otheris for NMF with the generalized Kullback-Leibler divergence, also known as
I-divergence. The main idea behind their algorithms is to repeat the minimization of an auxiliary function with respect to a subset of the variables
until a certain stopping condition is satisfied. Because the minimization of theauxiliary function does not increase the error, it is guaranteed that the error
decreases monotonically. This idea was later applied to various error measures
such as Bregman divergence [11], a parametric divergence containing Euclidean
distance and I-divergence as special cases [27], α-divergence [7], Itakura-Saito
divergence [12], and β-divergence [13]. Furthermore, a unified method for deriving multiplicative update rules from a wide variety of error measures was
proposed by Yang and Oja [45]. Using this method, they obtained eleven different multiplicative update rules from a collection of error measures including
the ones mentioned above, Kullback-Leibler divergence, γ-divergence, Rényi
divergence, and so on. Like the multiplicative update rules of Lee and Seung,all update rules obtained by this method have the property that the error
decreases monotonically.

A unified global convergence analysis of multiplicative update rules for NMF

3

Although the original multiplicative update rules have some good properties as stated above, they have a serious drawback that the convergence to astationary point is not guaranteed theoretically. The error decreases monotonically, but this does not mean the convergence of the sequence of solutions. In
fact, it was shown by numerical experiments that the multiplicative update
rule based on Euclidean distance does not converge to an local optimal solution
for some problems [16]. Another drawback is that every multiplicative update
is not well-defined in the feasible region of the NMF optimization problem.
For example, if one of the two factor matrices is a zero matrix, the update
rule is not defined because the denominator is zero. Here we should note thatthese two drawbacks are not independent but closely related to each other.
Many authors have studied the convergence property of the multiplicative update rules proposed by Lee and Seung and their variants. Finesso andSpreij [14] considered the I-divergence based multiplicative update rule andproved the convergence of the left factor matrix and the product of the two
factor matrices under the assumption that the right factor matrix is normalized after each update. Lin [30] considered a modified version of the Euclidean
distance based multiplicative update rule and proved that the sequence of solutions converges to a stationary point in some sense. Gillis and Glineur [15]considered another modified version of the Euclidean distance based multiplicative update rule, which returns a user-specified positive constant if theoriginal update rule returns a value less than the constant, and proved that
if the sequence of solutions converges then it is a stationary point. Takahashi
and Hibi [36] proved that the multiplicative update rules based on Euclidean
distance and I-divergence combined with the modification of Gillis and Glineur
have the global convergence property in the sense that any sequence of solutions has at least one convergent subsequence and the limit of any convergent
subsequence is a stationary point. Badeau et al. [1] studied the local stabilityof a generalized multiplicative update rules by Lyapunov’s stability theory andshowed that the local optimal solution is asymptotically stable if one of thetwo factor matrices is fixed to a nonnegative constant matrix. Zhao and Tan
recently performed a unified convergence analysis of the multiplicative update
rules for NMF with ℓ1 regularization [47].
In this paper, we consider a wide class of error measures including Euclidean distance, I-divergence and all other divergences mentioned above, andprovide a general suﬃcient condition on the error measure and the auxiliary
function for the obtained multiplicative update rule combined with the modification of Gillis and Glineur [15] to have the global convergence property
in the same sense as Takahashi and Hibi [36]. We then apply our results tothe eleven error measures considered by Yang and Oja [45], and show that
all of them satisfy the suﬃcient condition for the global convergence except
Kullback-Leibler divergence, γ-divergence and Rényi divergence. Finally, we
propose three new error measures based on these divergences and show thatthe obtained update rules have the global convergence property.
It is often said that the multiplicative update rules are slow. In fact, many
faster algorithms have been proposed for NMF with Euclidean distance [8,4Norikazu Takahashi et al.

9, 15, 17, 21, 23–25, 31]. Nevertheless, the multiplicative update rules are still
important because of its simplicity, easiness to implement, and applicability
to a wide range of error measures. We therefore focus our attention on themultiplicative update rules in this paper, and clarify the global convergencecondition. For the convergence properties of other algorithms, see [6,19,25,26]
and references therein.
The rest of this paper is organized as follows. In Section 2, NMF andthe multiplicative update rule are briefly reviewed. In Section 3, the global
convergence of the modified multiplicative update rule is proved under somemild assumptions on the error measure and the auxiliary function. In Section 4,
the result in Section 3 is applied to the eleven multiplicative updates obtainedby the unified method of Yang and Oja, and their global convergence is proved
except Kullback-Leibler divergence, γ-divergence, and Rényi divergence. In
Section 5, three new multiplicative update rules corresponding to these threeerror measures are derived, and their global convergence is proved. In Section 6
we conclude the paper.
2 NMF and Multiplicative Update Rules
2.1 NMF Optimization Problem
m×n
Given a nonnegative matrix X ∈ R+
, where R+ denotes the set of nonnegative real numbers, we consider the problem of finding two nonnegativer×n
matrices W ∈ Rm×rand H ∈ R+
such that+

X ≈ WH(1)

where r is a positive integer less than m and n. The operation that decomposesa given nonnegative matrix X into two nonnegative factor matrices W and H
as shown in (1) is called nonnegative matrix factorization (NMF). Although
it is important in NMF how to choose the value of r, we do not consider thisissue in this paper. We simply assume that the value of r is given together
with X. Also, we assume throughout this paper that every row and columnof X has at least one nonzero entry. The problem of finding W and H in (1)is formulated as a constrained optimization problem:
minimize D(W , H)
subject to W ≥ Om×r , H ≥ Or×n(2)

where D(W , H) is a function representing the error between X and W H andOm×r (Or×n , resp.) is the m × r (r × n, resp.) zero matrix. The inequality
signs between matrices are to be interpreted componentwise. In the remainderof this paper, we call D(W , H) the error function. Also, we denote the feasible
region of (2) by F0 , that is, F0 = Rm×r
× Rr×n++ .
So far, various kinds of error functions such as the one based on Euclidean
distance:∑
D(W , H) =(Xij − (W H)ij )2(3)
ijA unified global convergence analysis of multiplicative update rules for NMF

5

where (W H)ij denotes the (i, j)-th entry of the matrix W H ∈ Rm×n , andthe one based on I-divergence:)∑(
Xij
D(W , H) =Xij ln− Xij + (W H)ij
(W H)ij
ij
have been used [45]. All of those error functions are in general continuous in
the interior of F0 , which is denoted by int F0 in this paper. However, someof them are not well-defined on the boundary of F0 . Hence, when using suchan error function, we allow it take the value +∞, that is, we regard it as anextended real-valued function defined on F0 .

2.2 Multiplicative Update Rules
Because (2) is not a convex optimization problem, it is diﬃcult to find a
global optimal solution. As an approach to find a local optimal solution, aclass of iterative methods called multiplicative update rules [14, 28, 29, 45] arewidely used. Given an error function D(W , H), a multiplicative update rule
is obtained in the following manner. First, an auxiliary function of the error
function is constructed in some way. The definition of the auxiliary function
is given below. Second, for each variable, a problem of finding a unique minimum point of the auxiliary function under some constraints is solved. If theminimum point is explicitly expressed in terms of the current values of W andH, this expression leads to the multiplicative update rule.
Let us now give a formal definition of the auxiliary function.
Definition 1 (Auxiliary Function) Given an error function D(W , H), any
f , H)
f : int F0 × int F0 → R that satisfies
D̄(W , H, W
f , H)f ∈ int F0 × int F0 , D̄(W , H, W
f , H)
f ≥ D(W , H)∀(W , H, W

(4)∀(W , H) ∈ int F0 , D̄(W , H, W , H) = D(W , H)

(5)

and
is called an auxiliary function of D(W , H).
It is important to note that we use only a single auxiliary function in thispaper, while two auxiliary functions have been used in the literature: one isfor W and the other is for H. The reason for doing this is to simplify the
analysis. In fact, the update rule for the entries of W and the update rule forthe entries of H are obtained from a common auxiliary function. When two
auxiliary functions are given, one can immediately obtain a single auxiliary
function by averaging them. If the obtained single auxiliary function satisfies
Assumptions 1–3 given later, our method can be applied.
f , H)
f be an auxiliary function of D(W , H). Let {(W (l) , H (l) )}∞Let D̄(W , H, W
l=0
be any sequence such that i) (W (0) , H (0) ) ∈ int F0 and ii) for each l ≥ 0,

6Norikazu Takahashi et al.

W (l+1) is obtained from W (l) and H (l) as an optimal solution of the problem:
minimize D̄(W , H (l) , W (l) , H (l) )(6)subject to W > Om×r ,and iii) for each l ≥ 0, H (l+1) is obtained from W (l+1) and H (l) as an optimal
solution of the problem:
minimize D̄(W (l+1) , H, W (l+1) , H (l) )subject to H > Or×n .

(7)

Then the sequence {D(W (l) , H (l) )}∞l=0 is monotone decreasing because
D(W (l+1) , H (l+1) ) ≤ D̄(W (l+1) , H (l+1) , W (l) , H (l) )≤ D̄(W (l+1) , H (l) , W (l) , H (l) )≤ D̄(W (l) , H (l) , W (l) , H (l) )= D(W (l) , H (l) ) .
Here, the first inequality follows from (4), the second inequality follows fromthe fact that H (l+1) is an optimal solution of (7), the third inequality follows
from the fact that W (l+1) is an optimal solution of (6), and the equality follows
from (5).
For example, using Lemma 9 given in Appendix A, we can obtain an auxiliary function of (3) as follows:
2
2
∑
∑
∑ Wik
∑
Hkj
2
f H)
f ij
−2
(W
Xij (W H)ij +
Xij
.
fik H
e kj
W
ij
ij
ij
k
(8)
If (W (l) , H (l) ) ∈ int F0 then the problem (6) with (8) has the unique optimal solution, which is explicitly expressed by W (l) and H (l) . Similarly, if(W (l+1) , H (l) ) ∈ int F0 then the problem (7) with (8) has the unique optimal solution, which is explicitly expressed by W (l+1) and H (l) . From these
expressions, a multiplicative update rule described by

f , H)f =D̄(W , H, W

(l+1)

= Wik

(l+1)

= Hkj

Wik

Hkj

(X(H (l) )T )ik,(W (l) H (l) (H (l) )T )ik(9)((W (l+1) )T X)kj((W (l+1) )T W (l) H (l) )kj

(10)

(l)

(l)

is obtained [28, 29]. Note that the right-hand sides of (9) and (10) are positive
because we have assumed that every row and column of X has at least one
nonzero entry. Therefore, if (W (0) , H (0) ) ∈ int F0 then (W (l) , H (l) ) ∈ int F0
for all l ≥ 1. Note also that the sequence {D(W (l) , H (l) )}∞l=1 converges tosome constant because it is monotone decreasing and 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). In addition, even if the sequence converges
to some point, it is not guaranteed that the limit point belongs to int F0 .

A unified global convergence analysis of multiplicative update rules for NMF

7

The approach described above is not restricted to Euclidean distance but
can be applied to various error functions. In fact, Yang and Oja [45] proposed aunified method to develop multiplicative update rules and applied it to elevenerror functions. Details will be given in Section 4.
2.3 Modified Multiplicative Update Rules
The multiplicative update rule given by (9) and (10) and other multiplicative
update rules [7, 11–13, 27, 45] have a common serious problem that they are
not defined for all points in F0 . For example, if all entries in the k-th row of
H (l) are zero then the denominator of the right-hand side of (9) becomes zero.
One may think this is not a serious issue because if (W (0) , H (0) ) belongs to
int F0 then so does (W (l) , H (l) ) for all l ≥ 0, as mentioned above. However,it may occur that some of the entries of W (l) and H (l) converges to zero as
l → ∞. In this case, some of the entries of W (l) and H (l) may go to infinity
with l. In order to avoid this situation, Gillis and Glineur [15] modified (9)and (10) as()(X(H (l) )T )ik
(l+1)
(l)
Wik
,
(11)
= max ϵ, Wik(W (l) H (l) (H (l) )T )ik)(((W (l+1) )T X)kj
(l+1)
(l)
,
(12)
Hkj
= max ϵ, Hkj((W (l+1) )T W (l) H (l) )kj
respectively, where ϵ is a user-specified positive constant. This simple modification is general and can be applied to all multiplicative update rules.
When a modified multiplicative update rule is used, it is natural to consider
a modified optimization problem:
minimize D(W , H)
subject to W ≥ ϵ1m×r , H ≥ ϵ1r×n(13)instead of the original optimization problem (2), where 1m×r (1r×n , resp.) is
the m × r (r × n, resp.) matrix consisting of all ones. Let the feasible regionof the problem (13) be denoted by Fϵ , that is,Fϵ = {(W , H) | W ≥ ϵ1m×r , H ≥ ϵ1r×n } .
If (W (0) , H (0) ) belongs to Fϵ then the sequence {(W (l) , H (l) )}∞l=0 generatedby (11) and (12) is contained in Fϵ . Moreover, if some entry of W (l) and H (l)converges as l → ∞, its limit point belongs to [ϵ, ∞) because this is a closed
set.
Although the original multiplicative update rule given by (9) and (10) often
produces sparse matrices, the modified update rule given by (11) and (12)always produces dense matrices because ϵ is positive. However, if we replace
all ϵ in the factor matrices produced by the modified update rule with zero, theresulting matrices are expected to be sparse because local optimal solutions of(2) or (13) are often located at the boundary of the feasible region.

8

Norikazu Takahashi et al.

3 Global Convergence of Modified Multiplicative Update Rules
3.1 Problem Setting and Main Result
In this section, we consider a general error function D(W , H), which is definedon F0 as an extended real-valued function and satisfies the following properties:i) it is continuously diﬀerentiable on int F0 (see Assumption 1 given below),ii) D(W , H) ≥ 0 for all (W , H) ∈ F0 , and iii) D(W , H) = 0 if W H = X.
We then give a suﬃcient condition on D(W , H) and the auxiliary function
f , H)
f for the modified update rule to have the global convergence
D̄(W , H, W
property in the sense of Zangwill [46]. For the update rule we are considering,
it is defined as follows.
Definition 2 (Global Convergence) Let ϵ be a user-specified positive constant. An update rule is said to have the global convergence property if, for
any initial point (W (0) , H (0) ) ∈ Fϵ , the sequence {(W (l) , H (l) )}∞l=0 generated
by the update rule has at least one convergent subsequence and the limit ofany convergent subsequence is a stationary point of the optimization problem(13).
A stationary point of (13) is a point (Ŵ , Ĥ) ∈ Fϵ such that the following
conditions are satisfied (see, e.g., [36]):∇W D(Ŵ , Ĥ) ≥ Om×r ,(14)∇H D(Ŵ , Ĥ) ≥ Or×n ,(15)∇W D(Ŵ , Ĥ) ⊙ (ϵ1m×r − Ŵ ) = Om×r ,(16)∇H D(Ŵ , Ĥ) ⊙ (ϵ1r×n − Ĥ) = Or×n ,(17)

where ∇W D(Ŵ , Ĥ) (∇H D(Ŵ , Ĥ), resp.) is the m×r (r ×n, resp.) matrix ofwhich the (i, k)-th ((k, j)-th, resp.) entry is the value of ∂D/∂Wik (∂D/∂Hkj ,resp.) at (Ŵ , Ĥ) and ⊙ denotes the componentwise multiplication of two
matrices of the same dimension. In this paper, the set of stationary points of(13) is denoted by Sϵ .
We now give three assumptions about the error function and the auxiliary
function.
f , H)
f are continuously diﬀerenAssumption 1 D(W , H) and D̄(W , H, W
tiable on int F0 and int F0 × int F0 , respectively.
Assumption 2 For any (Ŵ , Ĥ) ∈ int F0 , the following equalities hold:∇W D̄(Ŵ , Ĥ, Ŵ , Ĥ) = ∇W D(Ŵ , Ĥ) ,∇H D̄(Ŵ , Ĥ, Ŵ , Ĥ) = ∇H D(Ŵ , Ĥ)where ∇W D̄(Ŵ , Ĥ, Ŵ , Ĥ) (∇H D̄(Ŵ , Ĥ, Ŵ , Ĥ), resp.) is the m × r (r × n,
resp.) matrix of which (i, k)-th ((k, j)-th, resp.) entry is the value of ∂ D̄/∂Wik(∂ D̄/∂Hkj , resp.) at (Ŵ , Ĥ, Ŵ , Ĥ).

A unified global convergence analysis of multiplicative update rules for NMF

9f , H)
f satisfies the following conditions.
Assumption 3 D̄(W , H, W
f , H)
f is separable in the sense that it is expressed as
1. D̄(W , H, W∑
1
f , H)
f =f , H)
f + D̄2 (Wf , H)
f + constant .
D̄(W , H, W
D̄ijk
(Wik , Hkj , W
ijk

(18)
2. For each (Ŵ , Ĥ) ∈ int F0 ,uik (Wik ) ≜

∑

1
D̄ijk
(Wik , Ĥkj , Ŵ , Ĥ)j∗is strictly convex in R++ and minimized at Wik∈ R++ which is explicitly
expressed as∗Wik= fik (Ŵ , Ĥ) ,where R++ is the set of positive numbers. Also, for each (Ŵ , Ĥ) ∈ int F0 ,∑1
vkj (Hkj ) ≜D̄ijk(Ŵik , Hkj , Ŵ , Ĥ)i∗is strictly convex in R++ and minimized at Hkj∈ R++ which is explicitly
expressed as∗
Hkj= gkj (Ŵ , Ĥ) .

3. fik (W , H) is continuous in int F0 and, for each ϵ > 0, there exist cik > 0
and ϕik < 1 such that∀(W , H) ∈ Fϵ ,ϕikfik (W , H) ≤ cik Wik.

(19)

Also, gkj (W , H) is continuous in int F0 and, for each ϵ > 0, there exist
dkj > 0 and ψ kj < 1 such that∀(W , H) ∈ Fϵ ,

ψgkj (W , H) ≤ dkj Hkjkj .

(20)

Under Assumptions 1–3, we can prove the global convergence of the modified update rule, as stated in the following theorem. The proof is given in thenext subsection.
Theorem 1 Suppose that an error function D(W , H) and an auxiliary funcf , H)
f of it satisfy Assumptions 1–3. Then, for any positivetion D̄(W , H, W
constant ϵ, the modified multiplicative update rule described by
(l+1)

Wik

= max(ϵ, fik (W (l) , H (l) )) ,

(l+1)
Hkj
= max(ϵ, gkj (W (l+1) , H (l) ))

(21)
(22)

has the global convergence property, where fik and gkj are defined in Assumption 3.

10

Norikazu Takahashi et al.

If a modified update rule has the global convergence property, it is easy
to obtain an algorithm that terminates within a finite number of iterations
after finding an approximate solution. Let us relax the conditions (14)–(17) as
follows:∇W D(W , H) ≥ −δ1 1m×r ,∇H D(W , H) ≥ −δ1 1r×n ,

(23)
(24)

(∇W D(W , H))ik > δ1 ⇒ Wik − ϵ ≤ δ2 ,(∇H D(W , H))kj > δ1 ⇒ Hkj − ϵ ≤ δ2 ,

(25)
(26)

where δ1 and δ2 are positive constants specified by the user. Employing these
conditions as a stopping criterion, we obtain the following algorithm.
Algorithm 1 Modified Multiplicative Update Algorithm
Input: X ∈ Rm×n
, r ∈ N, ϵ > 0, δ1 > 0, δ2 > 0+
Output: A point (W , H) ∈ Fϵ that satisfies (23)–(26)1: Choose (W (0) , H (0) ) ∈ Fϵ and set l = 0.
2: Find (W (l+1) , H (l+1) ) by the update rule described by (21) and (22).
3: If (W , H) = (W (l+1) , H (l+1) ) ∈ Fϵ satisfies (23)–(26) then return (W (l+1) , H (l+1) ).
Otherwise add 1 to l and go to Step 2.

The finite termination of this algorithm is guaranteed by the following
theorem. We omit the proof because it is the same as that of Theorem 2 inReference [36].
Theorem 2 Suppose that an error function D(W , H) and its auxiliary funcf , H)f satisfy Assumptions 1–3. Then, for any positive contion D̄(W , H, W
stants ϵ, δ1 and δ2 , Algorithm 1 stops within a finite number of iterations.

3.2 Relation to Existing Work
There are many results on the convergence properties of NMF algorithms
in the literature [6, 15, 17, 21, 23, 25, 27, 29, 31]. However, most of them are
insuﬃcient to prove that an approximate stationary point is always obtainedby the algorithm. As mentioned before, the monotone decrease of the error
value does not imply the convergence of the sequence of solutions. It is claimed
in some papers [6, 15, 21, 23, 25, 31] that every limit point of the sequence ofsolutions is a stationary point. However, this claim says nothing about theconvergence of the sequence or even the existence of a convergent subsequence.
The convergence of the sequence of solutions generated by the multiplicative update rule was proved by some authors in diﬀerent settings [14, 30, 47].
Finesso and Spreij [14] considered a variant of the I-divergence based multiplicative update rule, and proved the convergence of the right factor matrix
and the product of the two factor matrices. Lin [30] proved the convergenceA unified global convergence analysis of multiplicative update rules for NMF

11

of a variant of the Euclidean distance based multiplicative update rule. A
common key feature of these two algorithms is that one of the two factor
matrices is normalized after each update. This guarantees the boundednessof the corresponding factor matrix. Zhao and Tan [47] recently performed aunified convergence analysis of the multiplicative update rules for NMF withℓ1 regularization. They proved the boundedness of the sequence of solutions
by showing that the level set of the objective function is bounded. However,this result relies on the existence of ℓ1 regularization terms, and cannot bedirectly applied to NMF without regularizations.
The unified global convergence analysis in the present paper is a generalization of the results given by Takahashi and Hibi [36]. They showed that themodification of Gillis and Glineur [15] is suﬃcient to guarantee the global convergence of the multiplicative update rules based on Euclidean distance andI-divergence. However, because these update rules were studied separately, it
was not clear what kind of conditions are required in general for the error
function and the auxiliary function in order for the obtained multiplicative
update rule to have the global convergence property.
Finally, it is important to make a comment on the stopping condition. Someauthors [17, 31] used a stopping condition similar to the one in Algorithm 1.
Roughly speaking, their condition corresponds to (23)–(26) with δ2 = 0. In
this case, the finite termination of the algorithm is not guaranteed.
3.3 Proof of Theorem 1
We prove Theorem 1 by Zangwill’s global convergence theorem [46], which
has been used extensively to prove the convergence of various algorithms (see[38, 42] for example).
Theorem 3 (Zangwill [46]) Let A be a point-to-set mapping defined on a
space V that assigns to every point x ∈ V a subset of V . Let {x(l) }∞l=0 be asequence generated satisfying x(0) ∈ V and x(l+1) ∈ A(x(l) ). Also let a solution
set Ω ⊂ V be given. Suppose that the following three statements hold true.
1. All points x(l) are in a compact set X ⊂ V .
2. There is a continuous function Z : V → R such that i) if x is not asolution then Z(y) < Z(x) for any y ∈ A(x), and ii) if x is a solution thenZ(y) ≤ Z(x) for any y ∈ A(x).
3. The mapping A is closed at points outside Ω.
Then the limit of any convergent subsequence of {x(l) }∞l=0 is a solution.
In the following discussion, we express for simplicity (21) and (22) as
W (l+1) = A1 (W (l) , H (l) ) and H (l+1) = A2 (W (l+1) , H (l) ), respectively. Furthermore, we express the update process from (W (l) , H (l) ) to (W (l+1) , H (l+1) )as (W (l+1) , H (l+1) ) = A(W (l) , H (l) ). Then the mapping A is expressed interms of A1 and A2 as follows:A(W (l) , H (l) ) = (A1 (W (l) , H (l) ), A2 (A1 (W (l) , H (l) ), H (l) )) .

12

Norikazu Takahashi et al.

Note that A is not a point-to-set mapping but a point-to-point mapping definedon Fϵ . In this case, the closedness of A in Theorem 3 reduces to the continuityof A. Therefore, it suﬃces to show that the following statements hold true.
1. (Boundedness) For any initial point (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 Z : Fϵ → R such that
(W , H) ∈/ Sϵ ⇒ Z(A(W , H)) < Z(W , H) ,(W , H) ∈ Sϵ ⇒ Z(A(W , H)) ≤ Z(W , H) ,where Sϵ is the set of stationary points of (13).
3. (Continuity) A is continuous in Fϵ \ Sϵ .
We first consider the first statement. Our proof is based on the following
lemma.
Lemma 1 (Katayama et al. [22]) Let ϵ be any positive constant. Let f bea mapping from [ϵ, ∞) to R. If there exist a positive constant c and a constantϕ less than 1 such that∀x ≥ ϵ, f (x) ≤ cxϕthen any sequence {x(l) }∞l=0 generated by the update rule:x(l+1) = max(ϵ, f (x(l) )),l = 0, 1, 2, . . .

with the initial value x(0) ≥ ϵ is contained in a closed and bounded set.
Some illustrative examples for Lemma 1 are shown in Fig. 1. In Fig. 1(a),two sequences generated by x(l+1) = max(ϵ, f (x(l) )) with ϵ = 0.1 and f (x) =1.5x0.3 are plotted. The sequence starting from x(0) = 0.2 increases monotonically and converges to the unique fixed point. The sequence starting fromx(0) = 3.2 decreases monotonically and converges to the fixed point. In Fig. 1(b),the sequence generated by x(l+1) = max(ϵ, f (x(l) )) with x(0) = 0.7, ϵ = 0.2
and f (x) = 1.5x−0.8 − 0.28 is plotted. We see from the figure that it converges
to periodic sequence with period 2.
In view of Assumption 3 and Lemma 1, we immediately obtain the following
lemma.
Lemma 2 Let ϵ be any positive constant. If Assumption 3 holds then, for any
initial point (W (0) , H (0) ) ∈ Fϵ , the sequence {(W (l) , H (l) )}∞l=0 generated by(21) and (22) is contained in a closed and bounded set.
We next consider the third statement. The continuity of the mapping A isshown as follows.
Lemma 3 Let ϵ be any positive constant. If Assumption 3 holds then the
mapping A : Fϵ → Fϵ is continuous.

A unified global convergence analysis of multiplicative update rules for NMF
3.5

6f(x(l))

3

13

f(x(l))

5

2.5x(l+1)

x(l+1)

4
2
1.5

3
2

1
1

0.5
0

0
0

0.5

1

1.5
x

2

2.5

3

3.5

0

1

(l)

2

3

4

5

6

(l)

x

(a)

(b)

Fig. 1 Illustrative examples for Lemma 1. (a) Two sequences generated by x(l+1) =max(ϵ, f (x(l) )) with ϵ = 0.1 and f (x) = 1.5x0.3 . (b) A sequence generated by x(l+1) =max(ϵ, f (x(l) )) with ϵ = 0.2 and f (x) = 1.5x−0.8 − 0.28.

Proof By Assumption 3, fik and gkj are continuous. Also, for any positive constant ϵ, the mapping max(ϵ, ·) is continuous. So (21), which is the compositionof fik and max(ϵ, ·), and (22), which is the composition of gkj and max(ϵ, ·),are both continuous. Furthermore, because the mapping A is the compositionof (21) and (22), it is continuous.
⊔⊓
We finally consider the second statement. As the first step, we discuss therelationship between the stationary points of (13) and the mapping A in thefollowing two lemmas.
Lemma 4 Let ϵ be any positive constant. Let (Ŵ , Ĥ) be any point in Fϵ . If
Assumption 3 holds then the optimization problem:
minimize D̄(W , Ĥ, Ŵ , Ĥ)subject to W ≥ ϵ1m×r(27)

has a unique optimal solution which is given by A1 (Ŵ , Ĥ). Similarly, if Assumption 3 hold then the optimization problem:
minimize D̄(Ŵ , H, Ŵ , Ĥ)subject to H ≥ ϵ1r×n(28)

has a unique optimal solution which is given by A2 (Ŵ , Ĥ).
Proof We prove only the first part because the second one can be proved in
the same way. Because D̄(W , Ĥ, Ŵ , Ĥ) is separable by Assumption 3, theproblem (27) can be divided into mr independent problems of the form:∑ 1
minimizej D̄ijk (Wik , Ĥkj , Ŵ , Ĥ)(29)subject to Wik ≥ ϵ .

14Norikazu Takahashi et al.

∑ 1
It follows from Assumption 3 that uik (Wik ) ≜ j D̄ijk(Wik , Ĥkj , Ŵ , Ĥ) has
a unique minimum point given by fik (Ŵ , Ĥ) ∈ R++ . If fik (Ŵ , Ĥ) ≥ ϵthen it is apparently the unique optimal solution of (29). Otherwise, ϵ is theunique optimal solution of (29) because, by Assumption 3, uik (Wik ) is strictly
convex in R++ . Therefore, the unique optimal solution of (29) is given bymax(ϵ, fik (Ŵ , Ĥ)) = (A1 (Ŵ , Ĥ))ik .
⊔⊓Lemma 5 Under Assumptions 1–3, the necessary and suﬃcient condition for(Ŵ , Ĥ) ∈ Fϵ to be a stationary point of (13) is that Ŵ and Ĥ are the uniqueoptimal solutions of (27) and (28), respectively.
Proof As shown in the proof of Lemma 4, the necessary and suﬃcient condition
for Ŵ ≥ ϵ1m×r to be the unique optimal solution of (27) is that∀i, k,

u′ik (Ŵik ){= 0, if Ŵik > ϵ ,≥ 0, if Ŵik = ϵ .

(30)

By Assumptions 2 and 3, we have
∂ D̄∂Wik (Ŵ ,Ĥ,Ŵ ,Ĥ)()= ∇W D̄(Ŵ , Ĥ, Ŵ , Ĥ)ik().
= ∇W D(Ŵ , Ĥ)u′ik (Ŵik ) =ik

Therefore, (30) is rewritten as∀i, k,() {= 0, if Ŵik > ϵ ,∇W D(Ŵ , Ĥ)≥ 0, if Ŵik = ϵ ,ikwhich is equivalent to the conjunction of (14) and (16). We can show in thesame way as above that the necessary and suﬃcient condition for Ĥ ≥ ϵ1r×n
to be the unique optimal solution of (28) is equivalent to the conjunction of(15) and (17).
⊔⊓Using Lemmas 4 and 5, we can prove the monotoneness of the error function
as follows.
Lemma 6 Let ϵ be any positive constant. Under Assumptions 1–3, the following propositions hold true:(Ŵ , Ĥ) ∈/ Sϵ ⇒ D(A(Ŵ , Ĥ)) < D(Ŵ , Ĥ) ,(31)(Ŵ , Ĥ) ∈ Sϵ ⇒ D(A(Ŵ , Ĥ)) = D(Ŵ , Ĥ) .

(32)

A unified global convergence analysis of multiplicative update rules for NMF

15

Proof We first prove (32). Let (Ŵ , Ĥ) be any point in Sϵ . Then, it follows
from Lemma 5 that Ŵ and Ĥ are the unique optimal solutions of (27) and(28), respectively. Furthermore, it follows from Lemma 4 that Ŵ = A1 (Ŵ , Ĥ)and Ĥ = A2 (Ŵ , Ĥ). Therefore, we haveA(Ŵ , Ĥ) = (A1 (Ŵ , Ĥ), A2 (A1 (Ŵ , Ĥ), Ĥ))= (A1 (Ŵ , Ĥ), A2 (Ŵ , Ĥ))= (Ŵ , Ĥ)which implies that D(A(Ŵ , Ĥ)) = D(Ŵ , Ĥ). We next prove (31). Note that
if (Ŵ , Ĥ) ∈/ Sϵ then, by Lemma 5, at least one of the following two statements
holds: i) Ŵ is not the unique optimal solution of (27), that is, Ŵ ̸= A1 (Ŵ , Ĥ)and ii) Ĥ is not the unique optimal solution of (28), that is, Ĥ ̸= A2 (Ŵ , Ĥ).
If the first statement holds, we have
D(Ŵ , Ĥ) = D̄(Ŵ , Ĥ, Ŵ , Ĥ)> D̄(A1 (Ŵ , Ĥ), Ĥ, Ŵ , Ĥ)= D(A1 (Ŵ , Ĥ), Ĥ) − D(A1 (Ŵ , Ĥ), Ĥ)+ D̄(A1 (Ŵ , Ĥ), Ĥ, Ŵ , Ĥ)≥ D(A1 (Ŵ , Ĥ), Ĥ)= D̄(A1 (Ŵ , Ĥ), Ĥ, A1 (Ŵ , Ĥ), Ĥ)≥ D̄(A1 (Ŵ , Ĥ), A2 (A1 (Ŵ , Ĥ), Ĥ), A1 (Ŵ , Ĥ), Ĥ)= D(A1 (Ŵ , Ĥ), A2 (A1 (Ŵ , Ĥ), Ĥ))
− D(A1 (Ŵ , Ĥ), A2 (A1 (Ŵ , Ĥ), Ĥ))
+ D̄(A1 (Ŵ , Ĥ), A2 (A1 (Ŵ , Ĥ), Ĥ), A1 (Ŵ , Ĥ), Ĥ)≥ D(A1 (Ŵ , Ĥ), A2 (A1 (Ŵ , Ĥ), Ĥ))
= D(A(Ŵ , Ĥ)) .
If the first statement does not hold but the second one does, we have
D(Ŵ , Ĥ) = D̄(Ŵ , Ĥ, Ŵ , Ĥ)> D̄(Ŵ , A2 (Ŵ , Ĥ), Ŵ , Ĥ)= D(Ŵ , A2 (Ŵ , Ĥ)) − D(Ŵ , A2 (Ŵ , Ĥ))
+ D̄(Ŵ , A2 (Ŵ , Ĥ), Ŵ , Ĥ)≥ D(Ŵ , A2 (Ŵ , Ĥ))
= D(A(Ŵ , Ĥ)) .
Therefore, the inequality D(A(Ŵ , Ĥ)) < D(Ŵ , Ĥ) holds in both cases.

⊓
⊔

16

Norikazu Takahashi et al.

Table 1 Error functions considered by Yang and Oja [45].
Name
Euclidean distanceI-divergence
Dual I-divergence
Itakura-Saito divergence
α-divergence
β-divergence
Log-Quad cost
αβ-Bregman divergence

Kullback-Leibler divergence
γ-divergence

D(W , H)∑2ij (Xij − (W H)ij )
)
∑ (
Xij
ij Xij ln (W H)ij − Xij + (W H)ij)
∑ ((W H)ij− (W H)ij + Xijij (W H)ij ln
Xij
(
(
)
)
∑
Xij
Xij
ij − ln (W H)ij + (W H)ij − 1
(
)
∑
1−α
1
α
(α ̸= 0, 1)
ij αXij + (1 − α)(W H)ij − Xij (W H)ijα(1−α)
)
(
β+1
β+1
β
β
∑
Xij −(W H)ij
Xij −(W H)ij−
(β ̸= 0, −1)
ij Xij
β
β+1
)
∑ (
Xij
2 + X ln
−
X
(X
−
(W
H)
)
ij + (W H)ij
ij
ij
ij
ij
(W H)ij(∑ββαα
ij Xij − Xij − (W H)ij + (W H)ij)−(α(W H)α−1− β(W H)β−1(α ≥ 1, 0 < β < 1)
ij
ij )(Xij − (W H)ij )()∑∑
Xij / pq Xpq
Xij∑ln (W H) / ∑ (W H)ij
pqij
pq Xpq
pq
(
(∑
)
(∑
)
1+γ
1+γ
1
ln
+ γ lnij Xijij (W H)ij
γ(1+γ)− (1 + γ) ln

Rényi divergence

1
ln
ρ−1

(
∑

(∑

(
ij

γ
ij Xij (W H)ij

Xij
pq Xpq

∑

)ρ (

)
)

(γ ̸= 0, −1)(W H)ij
pq (W H)pq∑

)1−ρ )(ρ > 0, ρ ̸= 1)

4 Application of Theorem 1 to Multiplicative Update Rules
Derived by the Unified Method of Yang and Oja
Theorem 1 can be applied to various multiplicative update rules to prove theirglobal convergence. In this section, we apply it to the eleven multiplicative
update rules derived from the error functions shown in Table 11 by the unified method proposed by Yang and Oja [45]. To make the analysis simpler,
we assume for the moment that X is a positive matrix. Then all error functions in Table 1 are well-defined on F0 as extended real-valued functions. For
each of the first eight error functions (Euclidean distance, I-divergence, Dual
I-divergence, Itakura-Saito divergence, α-divergence, β-divergence, Log-Quadcost and αβ-Bregman divergence), we show later that the auxiliary function
derived by the method of Yang and Oja satisfies all conditions in Assumptions 1–3. This means that the modified update rules obtained from theseerror functions have the global convergence property. In contrast, for the last
three error functions, the global convergence of the modified update rule cannot be proved by Theorem 1 because the inequalities (19) and (20) are notsatisfied. This issue is discussed in the next section.

1 The error function based on Kullback-Leibler divergence in Table 1 is slightly diﬀerent∑
from the one in [45]. Instead of assuming thatij Xij = 1, Xij has been replaced with∑Xij / pq Xpq so that the result can be applied to a general nonnegative matrix X.

A unified global convergence analysis of multiplicative update rules for NMF

17

Let D(W , H) be a given error function. The procedure of the unified
method of Yang and Oja to derive a multiplicative update rule from D(W , H)is summarized as follows.
1. If D(W , H) contains one or more logarithmic functions then we replaceeach of them with a generalized polynomial by using the following relationship:xµ − 1.
(33)
ln x = lim+µµ→0
2. Applying Lemmas 8–10 given in Appendix A and then taking the limitµ → 0+ if necessary, we can derive an auxiliary function of D(W , H).
Note that we may need to apply L’Hôpital’s rule when taking the limit.
The derived auxiliary function is expressed in the form of (18) and, for∑ 1
any (Ŵ , Ĥ) ∈ int F0 , both uik (Wik ) ≜j D̄ijk (Wik , Ĥkj , Ŵ , Ĥ) and∑ 1
vkj (Hkj ) ≜ i D̄ijk (Ŵik , Hkj , Ŵ , Ĥ) are strictly convex in R++ .
3. Let (Ŵ , Ĥ) ∈ int F0 be the current solution. Solving the equations u′ik (Wik ) =′(Hkj ) = 0, we obtain multiplicative update formulae for Wik and
0 and vkj
Hkj , respectively.
Carrying out this procedure for each of the eleven error functions shown in
Table 1, we obtain eleven auxiliary functions in the form of (18). Table 2 shows1
f , H)
f of the auxiliary function
the formula for the first term D̄ijk(Wik , Hkj , W
for each error function. It is easily seen that every auxiliary function shownin Table 2 satisfies Assumptions 1, Assumption 2, and the first and second
f , H)
fconditions of Assumption 3. It is also seen that every D̄1 (Wik , Hkj , W
ijk

is expressed as a1ijk (Wik Hkj )b1 + a2ijk (Wik Hkj )b2 where b1 ̸= b2 and a1ijk
and a2ijk are independent of Wik and Hkj . Therefore, letting uik (Wik ) ≜∑ 1
j D̄ijk (Wik , Ĥkj , Ŵ , Ĥ) where (Ŵ , Ĥ) ∈ int F0 and solving the equation′uik (Wik ) = 0, we obtain Wik = fik (Ŵ , Ĥ). The formulae for fik (W , H)
∑ 1
(Ŵik , Hkj ,are shown in Reference [45]2 . Similarly, letting vkj (Hkj ) ≜ i D̄ijk′
Ŵ , Ĥ) and solving the equation vkj (Hkj ) = 0, we obtain Hkj = gkj (Ŵ , Ĥ).
Finally, by simple algebraic manipulations, we can prove that fik (W , H) andgkj (W , H) are upper bounded for the first eight error functions. Upper bounds
for fik (W , H) are explicitly given in Reference [22].
Let us next consider the case where X has a zero entry. For Dual Idivergence and α-divergence with a negative α, there exists a pair (i, k) such
∑ 1
f , H)
f and fik (W , H) are not defined. There(Wik , Hkj , W
that both j D̄ijk
fore, these error functions need a stronger assumption that X is positive. For
Itakura-Saito divergence, β-divergence with β less than −1, γ-divergence withγ less than −1, D(W , H) contains a constant term which takes +∞. However,
if we construct a new error function from the remaining terms, it is well-defined∑ 1
f , H)
f
on F0 as an extended real-valued function. Also, both j D̄ijk(Wik , Hkj , W
2 In the case of Kullback-Leibler divergence, X
ij in fik (W , H) must be replaced with∑Xij / pq Xpq .

18

Norikazu Takahashi et al.

and fik (W , H) are well-defined on R++ × R++ × int F0 and int F0 , respectively, for all (i, k). As for other functions, D(W , H) is well-defined on F0 even
∑ 1
f , H)
f andthough some entries of X are zero. In addition, j D̄ijk(Wik , Hkj , W
fik (W , H) are well-defined on R++ × R++ × int F0 and int F0 , respectively,
for all (i, k).
From these observations, we obtain the following lemma.
Lemma 7 For Euclidean distance, I-divergence, Itakura-Saito divergence, αdivergence with a positive α, β-divergence, Log-Quad cost and αβ-Bregman
divergence, the auxiliary function shown in Table 2 satisfies all conditionsin Assumptions 1–3. For Dual I-divergence and α-divergence with a negative
α, the auxiliary function satisfies all conditions in Assumptions 1–3 if X ispositive.
By Lemma 7 and Theorem 1, we immediately obtain the following theorem.
Theorem 4 For Euclidean distance, I-divergence, Itakura-Saito divergence,
α-divergence with a positive α, β-divergence, Log-Quad cost and αβ-Bregman
divergence, the modified multiplicative update rule obtained by the unified
method of Yang and Oja has the global convergence property for any positive
constant ϵ. For Dual I-divergence and α-divergence with a negative α, themodified multiplicative update rule has the global convergence property forany positive constant ϵ if X is positive.

5 New Multiplicative Update Rules for Kullback-Leibler, Gamma
and Rényi Divergences
In the previous section, we proved the global convergence of the modified
multiplicative update rule described by (21) and (22) for the first eight error
functions shown in Table 1. As for the last three (Kullback-Leibler divergence,γ-divergence and Rényi divergence), the multiplicative update rule can bederived, but the boundedness of the sequence {(W (l) , H (l) )}∞l=0 cannot be
proved. A possible reason is that the value of the error function does notchange even though W and H are multiplied by any positive scalar. In otherwords, we can increase the values of nonzero entries of W and H as much as
we want, while keeping the value of the error function fixed. As a simple way
to avoid this situation, we add a regularization term
(
)2
∑
C ∑
Xij −
(W H)ij
(34)
2 ij
ij
to each error function, where C is any positive constant. Applying Lemmas 8–10, we obtain new auxiliary functions for Kullback-Leibler divergence,γ-divergence and Rényi divergence. How to derive the auxiliary function for
Kullback-Leibler divergence is explained in Appendix B. Although other auxiliary function is omitted due to space constraints, they are obtained in a

A unified global convergence analysis of multiplicative update rules for NMF

19

Table 2 Auxiliary functions obtained from the error functions in Table 1 by the methodof Yang and Oja.
Error function

1 (W , H , W
f , H)
f
D̄ijk
ik
kj

Dual I-divergencef H)
f ij (Wfik He kj )−1 (Wik Hkj )2−2Xij Wik Hkj + (W
−1 f ef
f−Xij (W H)ij Wik Hkj ln(Wik Hkj ) + Wik Hkj(())−1 f ffik H
e kj )−1 − 1 Wik Hkjln(Wik Hkj ) + ln Xij(W H)ij (WItakura-Saito divergence

f H)
f −1 Wik Hkj + Xij (Wf H)
f −2 (Wfik He kj )2 (Wik Hkj )−1
(W
ij
ij

α-divergence

1
1
α (W
f H)
f −α (Wfik He kj )α (Wik Hkj )1−α
Wik Hkj − α(1−α)
Xij
ij
α

β-divergence1) β > 1

f H)
f β−1 Wik Hkj + 1 (Wf H)
f β (Wfik H
e kj )−β (Wik Hkj )β+1
−Xij (Wij
ij
β+1

Euclidean distanceI-divergence(α ̸= 0, 1)

2) 0 < β ≤ 1

f H)
f β−1 (Wfik He kj )1−β (Wik Hkj )β− β1 Xij (Wij3) β < 0, β ̸= −1

1
f H)
f β (Wfik H
e kj )−β (Wik Hkj )β+1
+ β+1
(W
ij
1
f H)
f β−1 (Wf H)
f β Wik Hkj
fik He kj )1−β (Wik Hkj )β + (W− β Xij (Wijij

Log-Quad cost

αβ-Bregman divergence

Kullback-Leibler divergence

(
)
f H)
f ij + 1 (Wfik H
e kj )−1 (Wik Hkj )2
(W
2
(
)
−1
f
f
fik H
e kj ln(Wik Hkj )−Xij (W H)ij + 2 W()β(1−β)
f H)
f α−1 +
f H)
f β−1 (Wfik He kj )1−α (Wik Hkj )α(α − 1)(W
(W
ij
ij
α
)
(
α(α−1) f f α−2
β−2
f
f
f
e kj )2−β (Wik Hkj )β−1
(W H)ij + β(W H)ij
Xij (Wik H+
1−β(α ≥ 1, 0 < β < 1)∑∑
−1 W H−1 X (Wff
f H)
f −1 Wfik H
e kj ln(Wik Hkj )ij
ik kj − ( pq Xpq )pq (W H)pq )ij(

γ-divergence1) γ > 02) γ < 0, γ ̸= −1

1
1+γ

−
(∑

(∑
(∑

1) ρ > 1

2) 0 < ρ < 1

)−1
ff γ
f H)
f γ−1 Wfik H
e kj ln(Wik Hkj )Xij (W
pq Xpq (W H)pq
ij)−1
f f 1+γf H)
f γ Wik Hkj(W
pq (W H)pqij

− γ1Rényi divergence)−1
f f 1+γf H)
f γ (Wfik H
e kj )−γ (Wik Hkj )1+γ
(W
pq (W H)pq
ij

(∑

)−1
(
)1−γ
ff γ
f H)
f γ−1 Wfik H
e kj
Xij (W(Wik Hkj )γ
pq Xpq (W H)pq
ij

(∑
)
(
)ρ
ρ f e 1−ρ −1
ρ f f −ρ f e
1
− 1−ρ
Xij
(W H)ij Wik Hkj (Wik Hkj )1−ρ
pq Xpq (W H)pq
(∑
)−1
ff
+
Wik Hkj
pq (W H)pq(∑)
ρ f f 1−ρ −1
ρ f f −ρ f e−
Xij(W H)ij Wik Hkj ln(Wik Hkj )pq Xpq (W H)pq
(∑
)−1
f
f
+
Wik Hkj
pq (W H)pq

similar way. All of the three auxiliary functions satisfy Assumption 1, Assumption 2, and the first and second conditions of Assumption 3. Also, each1
f , H)
f can be expressed as a1ijk (Wik Hkj )b1 +a2ijk (Wik Hkj )b2
D̄ijk(Wik , Hkj , W
where b1 ̸= b2 and a1ijk and a2ijk are independent of Wik and Hkj . Therefore,
the multiplicative update rules can be obtained from the auxiliary functions.
The formulae for fik (W , H) for the three error functions are shown in Table 3. Furthermore, for each of the three error functions, it is easy to provethat fik (W , H) and gjk (W , H) satisfy the third condition of Assumption 3.
Upper bounds for fik (W , H) on Fϵ are shown in Table 4, and how to derive
the upper bound for Kullback-Leibler divergence is described in Appendix C.
From these results, we obtain the following theorem.

20

Norikazu Takahashi et al.

Table 3 Multiplicative update rules obtained from the last three error functions in Table 1
with regularization term.
Error function

fik (W , H)

Kullback-Leibler divergenceWikγ-divergenceRényi divergence

( ∑
)1
∑
∑
∑
2
( pq Xpq )−1 j Xij (W H)−1
pq Xpqj Hkj
ij Hkj + C
∑
∑
∑
∑
−1
( pq (W H)pq )j Hkj + C
pq (W H)pq
j Hkj
( ∑
)η
∑
∑
∑
( pq Xpq (W H)γpq )−1 j Xij (W H)γ−1
pq Xpqj Hkj
ij Hkj + C
Wik
∑
∑
∑
∑
γ
−1
( pq (W H)1+γ
pq )j (W H)ij Hkj + C
pq (W H)pqj Hkj{
1/(1 + γ), if γ > 0
where η =1/(1 − γ), if γ < 0, γ ̸= −1
( ∑
)η
∑
∑
∑
ρ
−ρ
ρ
−1
(W H)1−ρ( pq Xpq
pq )j Xij (W H)ij Hkj + C
pq Xpq
j Hkj
Wik
∑
∑
∑
∑
( pq (W H)pq )−1 j Hkj + C pq (W H)pq j Hkj{
1/ρ, if ρ > 1
where η =1, if 0 < ρ < 1

Table 4 Upper bounds for fik (W , H) in Table 3 on Fϵ .
Error function

Upper bound(Kullback-Leibler divergence(
γ-divergence(Rényi divergence

1
1 ∑
+ ϵn
pq Xpqϵ3 nrC

)1
2

1 ∑
1
+ ϵn
pq Xpqϵ3 mrC

)η

1
1 ∑
+ ϵn
pq Xpqϵ3 mrC

)η

1
2
Wik

{

1−η
Wik

where η =1−ηWik

where η ={1/(1 + γ), if γ > 0
1/(1 − γ), if γ < 0, γ ̸= −1
1/ρ, if ρ > 1
1, if 0 < ρ < 1

Theorem 5 For Kullback-Leibler divergence, γ-divergence and Rényi divergence with the regularization term (34), the modified multiplicative update
rule obtained by the unified method of Yang and Oja has the global convergence property for any positive constant ϵ.

6 Conclusions
A unified global convergence analysis of the multiplicative update rule for
NMF has been presented. We have given a suﬃcient condition on the error
function and the auxiliary function for a slightly modified version of the multiplicative update rule to have the global convergence property in the sensethat any sequence of solutions contains at least one convergent subsequenceand the limit of any convergent subsequence is a stationary point of the corresponding optimization problem. This result can be applied to a wide variety
of multiplicative update rules to examine their global convergence. In fact, wehave proved the global convergence of eleven diﬀerent multiplicative update
rules.
Recently, many variants of NMF have been proposed (see the review paper
by Wang and Zhang [41] and references therein). Extending the results in thispaper to such variants is an interesting and important direction for futureA unified global convergence analysis of multiplicative update rules for NMF

21

research. Another interesting direction is to develop a method for determining
the global convergence of the modified multiplicative update rule only from theerror function. The motivation of this problem comes from the observation thatthe multiplicative update rule is uniquely determined by the error function if
we assume that the auxiliary function is always derived by the method of Yang
and Oja.

A How to Derive Auxiliary Functions
In the unified method of Yang and Oja [45], an auxiliary function is systematically derived
from a given generalized polynomial by using three rules. They are described by the following lemmas. Because the mathematical expressions diﬀer from those in [45] due to theintroduction of the framework of a single auxiliary function, we provide proofs for the sakeof readers’ convenience.
(∑)dcLemma 8 Suppose that the error function is expressed as D(W , H) = a
ij bij (W H)ij
where a and c are nonzero constants, bij are positive constants, and d is a constant otherthan 0 or 1. If ξ(x) ≜ axd is convex in R++ , letf , H)f =a
D̄(W , H, W

(
∑

)d−1
f H)
fc
bij (W
ij

∑

ij

(
f H)
fc
bij (Wij

ij(W H)ij
f H)
f ij
(W

)cd
.

If ξ(x) is concave in R++ , let(f , H)f =a
D̄(W , H, W

∑

(
×

)d
f H)
fc
bij (W
ij

(
+ ad

ij

∑

)d−1
f H)
fc
bij (Wij

ij∑

bij (W H)cij −

ij

∑

f H)
fc
bij (Wij).

ij

f , H)
f is an auxiliary function of D(W , H), and satisfies the conditions in
Then D̄(W , H, W
Assumptions 1 and 2.
Proof There are two cases to consider: One is that ξ(x) ≜ axd is convex, and the other isthat ξ(x) is concave. In either case, it is easy to see that the following statements hold true:
1. D̄(Ŵ , Ĥ, Ŵ , Ĥ) = D(Ŵ , Ĥ) for all (Ŵ , Ĥ) ∈ F0 ,f , H)
f is diﬀerentiable at any point in int F0 × int F0 , and2. D̄(W , H, W
3. ∇W D̄(Ŵ , Ĥ, Ŵ , Ĥ) = ∇W D(Ŵ , Ĥ) and ∇H D̄(Ŵ , Ĥ, Ŵ , Ĥ) = ∇H D(Ŵ , Ĥ) for all(Ŵ , Ĥ) ∈ int F0 .
Therefore, it suﬃces for us to show thatf , H)
f ≥ D(W , H)
D̄(W , H, W
f , H)f ∈ int F0 × int F0 . Suppose first that ξ(x) = axd is convex in R++ .
for all (W , H, W
Then,∑for any numbers x11 , x12 , . . . , xmn and any positive numbers λ11 , λ12 , . . . , λmn suchthat ij λij = 1, it follows from Jensen’s inequality that

ξ

∑
ij





xij  = ξ 


(
)
∑
xij
xij
≤
λij ξ
.
λij ·
λij
λij
ij
ij

∑

22

Norikazu Takahashi et al.

f H)
fc /
Substituting xij = bij (W H)cij and λij = bij (Wijwe have

d
∑
c 

D(W , H) = abij (W H)ij
ij

≤a

∑

∑

ij


= a

∑

(

f H)
fc
bij (Wij
ff c
pq bpq (W H)pq

d−1

f H)
fc 
bij (W
ij

ij

∑

ff c
pq bpq (W H)pq into this equation,bij (W H)cij
f H)
f c / ∑ bpq (Wf H)
fc
bij (W
pq
pq
ij

∑

(
f H)
fc
bij (Wij

ij(W H)ij
f H)
f ij
(W

)d

)cd

f , H)f .
= D̄(W , H, W
Suppose next that ξ(x) = axd is concave in R++ . Then, for any positive numbers x and x̃,the following inequality holds:
ξ(x) ≤ ξ(x̃) + ξ ′ (x̃)(x − x̃) .
Substituting x =∑cij bij (W H)ij and x̃ =

D(W , H) = a 

≤ a

∑

d

d
f H)
fc 
bij (Wij

ij

+ ad 

ff c
ij bij (W H)ij into this inequality, we havebij (W H)cij 

ij

∑

∑

∑

d−1 
f H)
fc bij (W
ij



ij

∑
ij

bij (W H)cij −

∑


f H)
fc 
bij (Wij

ij

f , H)f= D̄(W , H, W⊓⊔

which completes the proof.
∑bij aij (W H)ijb
where aij are nonzero constants and b is a constant other than 0 or 1. If ξij (x) ≜ aij x is

Lemma 9 Suppose that the error function is expressed as D(W , H) =convex in R++ , letf , H)
f = aij (W
f H)
f b−1
D̄ij (W , H, Wij∑

e kj )1−b (Wik Hkj )b .
fik H
(W

k

If ξij (x) is concave in R++ , let()
f , H)
f = aij (Wf H)
f b + aij b(Wf H)
f b−1 (W H)ij − (Wf H)
f ij .
D̄ij (W , H, W
ijijf , H)
f = ∑ D̄ij (W , H, W
f , H)
f is an auxiliary function of D(W , H), andThen D̄(W , H, W
ijsatisfies the conditions in Assumptions 1 and 2.
Proof Let Dij (W , H) = aij (W H)bij . There are two cases to consider: One is that ξij (x) ≜aij xb is convex and the other is that ξij (x) is concave. In either case, we easily see that thefollowing statements hold true:
1. D̄ij (Ŵ , Ĥ, Ŵ , Ĥ) = Dij (Ŵ , Ĥ) for all (Ŵ , Ĥ) ∈ int F0 ,f , H)
f is diﬀerentiable at any point in int F0 × int F0 , and2. D̄ij (W , H, W

A unified global convergence analysis of multiplicative update rules for NMF

23

3. ∇W D̄ij (Ŵ , Ĥ, Ŵ , Ĥ) = ∇W Dij (Ŵ , Ĥ) and ∇H D̄ij (Ŵ , Ĥ, Ŵ , Ĥ) = ∇H Dij (Ŵ , Ĥ)for all (Ŵ , Ĥ) ∈ int F0 .
Therefore, it suﬃces for us to show thatf , H)
f ≥ Dij (W , H)
D̄ij (W , H, W
f , H)f ∈ int F0 × int F0 . Suppose first that ξij (x) = aij xb is convex in R++ .
for all (W , H, W
Then,
for
any
numbers
x1 , x2 , . . . , xr and any positive numbers λ1 , λ2 , . . . , λr such that∑k λk = 1, it follows from Jensen’s inequality that
ξij

(
∑

)
xk

= ξij

(
∑

k

xk
λk

λk ·

k

)
≤

∑

(
λk ξij

k

xk
λk

)
.

f H)
f ij into this equation, we havefik H
e kj /(WSubstituting xk = Wik Hkj and λk = W
Dij (W , H) = aij (W H)bij
≤

fik H
e kj
∑ W
k

f H)
f ij
(W

(
aij

f H)
f b−1
= aij (Wij

Wik Hkj
f H)
f ij
fik H
e kj /(W
W

)b

∑
fik H
e kj )1−b (Wik Hkj )b(Wkf , H)f .
= D̄ij (W , H, W
Suppose next that ξij (x) = aij xb is concave in R++ . Then, for any positive numbers x andx̃, the following inequality holds:′
ξij (x) ≤ ξij (x̃) + ξij(x̃)(x − x̃) .

f H)
f ij into this inequality, we have
Substituting x = (W H)ij and x̃ = (WDij (W , H) = aij (W H)bij

(
)
f H)
f b + aij b(Wf H)
f b−1 (W H)ij − (W
f H)
f ij
≤ aij (Wijijf , H)f= D̄ij (W , H, W⊓⊔

which completes the proof.
∑

Lemma 10 Suppose that the error function is expressed as D(W , H) = tijk atijk (Wik Hkj )bt
where atijk are nonzero constants and bt are constants, atijk xbt is convex in R++ , and{bt } contains at least two distinct nonzero numbers. Let bmax = max{bt | bt ̸= 0} andf , H)
f on int F0 × int F0 as follows:bmin = min{bt | bt ̸= 0}. Let us define D̄tijk (W , H, W
1. If bt ∈ {bmin , bmax , 0}, letf , H)
f = atijk (Wik Hkj )bt .
D̄tijk (W , H, W
2. If bt ̸∈ {bmin , bmax , 0} and(a) if (bt > 1) ∨ ((bt = 1) ∧ (atijk > 0)), letf , H)f =D̄tijk (W , H, W)batijk bt f e bt −bmax ((Wik Hkj )Wik Hkj max
bmax
(
)
fik H
e kj )bt 1 − bt
+ atijk (W
,
bmax

24

Norikazu Takahashi et al.
(b) if (bt < 1) ∨ ((bt = 1) ∧ (atijk < 0)), letf , H)f =D̄tijk (W , H, W

f , H)
f =Then D̄(W , H, W)batijk bt f e bt −bmin ((Wik Hkj )Wik Hkj min
bmin
(
)
fik H
e kj )bt 1 − bt
+ atijk (W
.
bmin

∑

f f

tijk D̄tijk (W , H, W , H) is an auxiliary function of D(W , H),

f , H)
f satisfies the conditions in Asand strictly convex in int F0 . Furthermore, D̄(W , H, W
sumptions 1 and 2.

Proof Let Dtijk (W , H) = atijk (Wik Hkj )bt . There are three cases to consider dependingf , H)
f ison the values of atijk and bt . In either case, we easily see that D̄tijk (W , H, W
diﬀerentiable at any point in int F0 × int F0 and that
∇W D̄tijk (Ŵ , Ĥ, Ŵ , Ĥ) = ∇W Dtijk (Ŵ , Ĥ) ,∇H D̄tijk (Ŵ , Ĥ, Ŵ , Ĥ) = ∇H Dtijk (Ŵ , Ĥ)hold for all (Ŵ , Ĥ) ∈ F0 . Also, it has already been shown by Yang and Oja [45, Lemma 2]
f , H)
f is an auxiliary function of D(W , H).
that D̄(W , H, W⊓⊔

B Derivation of Auxiliary Function for Kullback-Leibler
Divergence with Regularization Term
We derive an auxiliary function for Kullback-Leibler divergence with the regularization termby using the unified method of Yang and Oja. First of all, we rewrite the error function byusing (33) as follows:)
1(
D1 (W , H) + D2 (W , H) + D3 (W , H) + D4 (W , H)µ→0+ µ
)2
)
(
(
∑
Xij
Xij
C ∑
∑
+
Xij
ln ∑
+
2
pq Xpqpq Xpqij
ij

D(W , H) = lim

where
D1 (W , H) =(∑)µ
(W H)ij

,

ij

∑

Xij
(W H)µij ,pq Xpqij∑∑
D3 (W , H) = −µC
Xij ·(W H)ij ,D2 (W , H) = −µC
D4 (W , H) =
2

(

∑

ij

ij

∑
(W H)ij)2.

ij

Let us assume that µ is a suﬃciently small positive constant. Applying Lemmas 8 and 9 tothese functions, we have the following auxiliary functions:f , H)f =µ
D1 (W , H, W

(
(
)µ−1
)µ
∑
∑
∑
f
f
f
f
(W H)ij
Wik Hkj + (1 − µ)(W H)ij
,
ij

ijk

ij

A unified global convergence analysis of multiplicative update rules for NMF
∑
Xij
f H)
f µ−1
fik H
e kj )1−µ (Wik Hkj )µ ,
(W
(W
ij
X
pq
pq
ij
k
∑
∑
f , H)
f = −µC
D3 (W , H, W
Xij ·Wik Hkj ,f , H)f =−D2 (W , H, W

∑

25

∑

ij

ijk

)−1
∑
∑(
f H)
f ij ·
f , H)
f = µC
fik H
e kj
(W
W
D4 (W , H, W
(Wik Hkj )2 .
2 ij
ijk
The exponents of Wik Hkj in these auxiliary functions are 1, µ, 1 and 2. The minimum is µ
and the maximum is 2. So we apply Lemma 10 to some of these auxiliary functions to obtainan auxiliary function of D(W , H) such that the exponents of Wik Hkj are restricted to µ
and 2. Applying Lemma 10 to D1 , we obtain another auxiliary function of D1 as follows:f , H)f =D1 (W , H, W

(
)µ−1
∑
µ ∑ fffik H
e kj )−1 (Wik Hkj )2(W(W H)ij
2 ij
ijk
(
)µ
(
µ) ∑ ff+ 1−(W H)ij.
2
ij

Applying Lemma 10 to D3 , we obtain another auxiliary function of D3 as follows:f , H)f
D3 (W , H, W
∑
∑
∑
∑
f H)
f ij .
fik H
e kj )1−µ (Wik Hkj )µ + (1 − µ)C
= −C
Xij ·
(W
Xij ·
(W
ij

ij

ijk

ij

As a result, we have the following auxiliary function:(
1f
f
f , H)
f + D2 (W , H, W
f , H)f
D̄(W , H, W , H) = lim
D1 (W , H, Wµ→0+ µ)ff
f
f+ D3 (W , H, W , H) + D4 (W , H, W , H)

+

∑
ij

∑

Xij
ln
pq Xpq(Xij
pq Xpq

∑

)
+

C
2

(
∑

)2
Xij

.

ij

Because
(
lim

µ→0+

f , H)
f + D2 (W , H, W
f , H)f
D1 (W , H, W)
f , H)
f + D4 (W , H, W
f , H)f+ D3 (W , H, W

we apply L’Hôpital’s rule. Then we havef , H)f =D̄(W , H, W)−1(∑
1 ∑ fffik H
e kj )−1 (Wik Hkj )2(W(W H)ij
2 ij
ijk
(
)
∑
1
f H)
f ij
− + ln
(W
2
ij
−

∑
ij

∑

(
)∑
Xij
f H)
f −1 ln (W
f H)
f ij
fik H
e kj
(W
W
ij
pq Xpq
k

= 0,

26

Norikazu Takahashi et al.
)
(
∑
Xij
f H)
f −1
fik H
e kj ln Wik Hkj
W
(W
ij
e kj
fik H
W
pq Xpq
ij
k
(
)
∑
∑
fik H
e kj ln Wik Hkj
−C
Xij ·
W
fik H
e kj
W
ij
ijk
∑
∑
f H)
f ij
(W
−C
Xij ·
−

∑

∑

ij

ij

∑
C ∑ fffik H
e kj )−1 (Wik Hkj )2+(W H)ij ·
(W
2 ij
ijk
(
)
(
)2
∑
Xij
Xij
C ∑
+
∑
ln ∑
+
Xij
2
pq Xpqpq Xpqij
ij
which can be rewritten in the form of (18) with1
f , H)
f = 1
D̄ijk(Wik , Hkj , W
2

(

∑
f H)
f pq
(W

)−1

fik H
e kj )−1 (Wik Hkj )2
(W

pq

Xij
f H)
f −1 Wfik H
e kj ln(Wik Hij )(W
ijpq Xpq∑fik H
e kj ln(Wik Hkj )−C
Xpq · W−∑pq

C ∑ fffik H
e kj )−1 (Wik Hkj )2 .
(W H)pq · (W+2 pq

C Derivation of Upper Bound for Multiplicative Update Rule
Obtained from Kullback-Leibler Divergence With Regularization
Term
For the first multiplicative update rule shown in Table 3, which is obtained from KullbackLeibler divergence with the regularization term, we derive an upper bound for fik (W , H)
on Fϵ . By simple mathematical manipulations, we have the following inequalities:∑fik (W , H) < Wik 

j

1

∑
∑
Xij
2
(W H)−1
pq Xpqj Hkj ij Hkj + C
pq XpqC∑≤ Wik ∑∑pq (W H)pq

∑

j Hkj

1
2
∑
(W H)−1
ij Hkjpq Xpq∑∑+∑ .
C pq (W H)pq j Hkj
pq (W H)pq
j

∑

Xij

pq XpqHere note that
∑
j

∑

Xij

pq Xpq∑−1(W H)ij
Hkj

j Hkj

∑

Hkj
Xij
(W H)−1
ij ∑
X
pq
pq
q Hkq
j
∑
1
j Xij
< 2 ∑
ϵ r pq Xpq
=

<

1
ϵ2 r

∑



A unified global convergence analysis of multiplicative update rules for NMF
and∑1
pq (W H)pq

< ∑1
q Wik Hkq

Therefore we have12

fik (W , H) ≤ Wik

(

=

Wik

1
∑

q Hkq

1
1 ∑
+
Xpq
ϵ3 nrC
ϵn pq

≤

27

1
.
ϵnWik

)1
2

.

Acknowledgements This work was partially supported by JSPS KAKENHI Grant Number JP15K00035.

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. Berman, A., Plemmons, R.: Nonnegative Matrices in the Mathematical Sciences. Academic Press, New York (1979)
3. Berry, M.W., Browne, M.: Email surveillance using non-negative matrix factorization.
Computational and Mathematical Organization Theory 11, 249–264 (2005)
4. Campbell, S.L., Poole, G.D.: Computing nonnegative rank factorizations. Linear Algebra and its Applications 35, 175–182 (1981)
5. Chen, J.C.: The nonnegative rank factorizations of nonnegative matrices. Linear algebra
and its applications 62, 207–217 (1984)
6. Chi, E.C., Kolda, T.G.: On tensors, sparsity, and nonnegative factorizations. SIAM
Journal on Matrix Analysis and Applications 33(4), 1272–1299 (2012)
7. Cichocki, A., Lee, H., Kim, Y.D., Choi, S.: Non-negative matrix factorization with αdivergence. Pattern Recognition Letters 29(9), 1433–1440 (2008)
8. Cichocki, A., Phan, A.H.: Fast local algorithms for large scale nonnegative matrix andtensor factorization. IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences E92-A(3), 708–721 (2009)
9. Cichocki, A., Zdunek, R., Amari, S.I.: Hierarchical ALS algorithms for nonnegative
matrix and 3D tensor factorization. In: Lecture Notes in Computer Science, vol. 4666,
pp. 169–176. Springer (2007)
10. Cichocki, A., Zdunek, R., Phan, A.H., Amari, S.I.: Nonnegative Matrix and Tensor
Factorizations. John Wiley & Sons, West Sussex, U.K. (2009)
11. Dhillon, I.S., Sra, S.: Generalized nonnegative matrix approximations with Bregman
divergences. In: Advances in Neural Information Processing Systems, pp. 283–290 (2005)
12. 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)
13. Févotte, C., Idier, J.: Algorithms for nonnegative matrix factorization with the βdivergence. Neural Computation 23(9), 2421–2456 (2011)
14. Finesso, L., Spreij, P.: Nonnegative matrix factorization and I-divergence alternating
minimization. Linear Algebra and its Applications 416, 270–287 (2006)
15. Gillis, N., Glineur, F.: Nonnegative factorization and the maximum edge biclique problem. arXiv e-prints (2008)
16. Gonzalez, E.F., Zhang, Y.: Accelerating the Lee-Seung algorithm for non-negative matrix factorization. Dept. Comput. & Appl. Math., Rice Univ., Houston, TX, Tech. Rep.
TR-05-02 (2005)
17. Guan, N., Tao, D., Luo, Z., Yuan, B.: NeNMF: an optimal gradient method for nonnegative matrix factorization. IEEE Transactions on Signal Processing 60(6), 2882–2898(2012)
18. Guillamet, D., Vitria, J.: Non-negative matrix factorization for face recognition. In:
Lecture Notes in Artificial Intelligence, pp. 336–344. Springer (2002)

28

Norikazu Takahashi et al.

19. Hansen, S., Plantenga, T., Kolda, T.G.: Newton-based optimization for Kullback-Leibler
nonnegative tensor factorizations. Optimization Methods and Software 30(5), 1002–
1029 (2015)20. Holzapfel, A., Stylianou, Y.: Musical genre classification using nonnegative matrix
factorization-based features. IEEE Transactions on Audio, Speech, and Language Processing 16(2), 424–434 (2008)
21. Hsieh, C.J., Dhillon, I.S.: Fast coordinate descent methods with variable selection for
non-negative matrix factorization. In: Proceedings of the 17th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 1064–1072. ACM(2011)
22. Katayama, J., Takahashi, N., Takeuchi, J.: Boundedness of modified multiplicative updates for nonnegative matrix factorization. In: Proceedings of the Fifth International
Workshop on Computational Advances in Multi-Sensor Adaptive Processing, pp. 252–
255 (2013)
23. Kim, D., Sra, S., Dhillon, I.S.: Fast newton-type methods for the least squares nonnegative matrix approximation problem. In: Proceedings of the Sixth SIAM International
Conference on Data Mining, pp. 343–354. SIAM (2007)
24. Kim, H., Park, H.: Nonnegative matrix factorization based on alternating nonnegativity
constrained least squares and active set method. SIAM Journal on Matrix Analysis andApplications 30(2), 713–730 (2008)
25. Kim, J., He, Y., Park, H.: Algorithms for nonnegative matrix and tensor factorization: a unified view based on block coordinate descent framework. Journal of Global
Optimization 58(2), 285–319 (2014)
26. Kimura, T., Takahashi, N.: Global convergence of a modified HALS algorithm for nonnegative matrix factorization. In: Proceedings of 2015 IEEE 6th International Workshop
on Computational Advances in Multi-Sensor Adaptive Processing, pp. 21–24 (2015)
27. Kompass, R.: A generalized divergence measure for nonnegative matrix factorization.
Neural Computation 19(3), 780–791 (2007)
28. Lee, D.D., Seung, H.S.: Learning the parts of objects by non-negative matrix factorization. Nature 401, 788–792 (1999)
29. 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)
30. Lin, C.J.: On the convergence of multiplicative update algorithms for nonnegative matrix factorization. IEEE Transactions on Neural Networks 18(6), 1589–1596 (2007)
31. Lin, C.J.: Projected gradient methods for non-negative matrix factorization. NeuralComputation 19(10), 2756–2779 (2007)
32. Paatero, P., Tapper, U.: Positive matrix factorization: a non-negative factor model with
optimal utilization of error estimates of data values. Environmetrics 5(2), 111–126(1994)
33. Panagakis, Y., Kotropoulos, C., Arce, G.R.: Non-negative multilinear principal component analysis of auditory temporal modulations for music genre classification. IEEE
Transactions on Audio, Speech, and Language Processing 18(3), 576–588 (2010)
34. Seki, M., Takahashi, N.: New updates based on Kullback-Leibler, gamma, and Rényi
divergences for nonnegative matrix factorization. In: Proceedings of 2014 International
Symposium on Nonlinear Theory and its Applications, pp. 48–51 (2014)35. 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)
36. Takahashi, N., Hibi, R.: Global convergence of modified multiplicative updates for nonnegative matrix factorization. Computational Optimization and Applications 57, 417–
440 (2014)
37. Takahashi, N., Katayama, J., Takeuchi, J.: A generalized suﬃcient condition for global
convergence of modified multiplicative updates for NMF. In: Proceedings of 2014 International Symposium on Nonlinear Theory and its Applications, pp. 44–47 (2014)
38. Takahashi, N., Nishi, T.: Global convergence of decomposition learning methods forsupport vector machines. IEEE Transactions on Neural Networks 17(6), 1362–1369(2006)

A unified global convergence analysis of multiplicative update rules for NMF

29

39. Vavasis, S.A.: On the complexity of nonnegative matrix factorization. SIAM Journal
on Optimization 20(3), 1364–1377 (2009)
40. Wang, R.S., Zhang, S., Wang, Y., Zhang, X.S., Chen, L.: Clustering complex networks
and biological networks by nonnegative matrix factorization with various similarity
measures. Neurocomputing 72, 134–141 (2008)
41. Wang, Y.X., Zhang, Y.J.: Nonnegative matrix factorization: A comprehensive review.
IEEE Transactions on Knowledge and Data Engineering 25(6), 1336–1353 (2013)
42. Wu, C.F.J.: On the convergence properties of the EM algorithm. The Ananls of Statistics
11(1), 95–103 (1983)
43. Xu, W., Liu, X., Gong, Y.: Document clustering based on non-negative matrix factorization. In: Proceedings of the 26th Annual International ACM SIGIR Conference on
Research and Development in Informaion Retrieval, pp. 267–273. ACM (2003)
44. Yamauchi, S., Kawakita, M., Takeuchi, J.: Botnet detection based on non-negative matrix factorization and the MDL principle. In: Proceedings of 19th International Conference on Neural Information Processing, pp. 400–409. Springer (2012)
45. Yang, Z., Oja, E.: Unified development of multiplicative algorithm for linear andquadratic nonnegative matrix factorization. IEEE Transactions on Neural Networks
22(12), 1878–1891 (2011)
46. Zangwill, W.: Nonlinear Programming: A Unified Approach. Prentice-Hall, Englewood
Cliﬀs, NJ (1969)
47. Zhao, R., Tan, V.Y.: A unified convergence analysis of the multiplicative update algorithm for nonnegative matrix factorization. arXiv preprint arXiv:1609.00951 (2016)

Erratum Sheet
Norikazu Takahashi, Jiro Katayama, Masato Seki and Jun’ichi Takeuchi, “A unified global convergence analysis of multiplicative update rules for nonnegative matrix factorization”, Computational
Optimization and Applications, vol.71, no.1, pp.221–250, 2018.
1. Eq.(10) on Page 226 is not correct. It should be replaced with
(l+1)

Hkj

(l)

= Hkj

((W (l+1) )T X)kj((W (l+1) )T W (l+1) H (l) )kj(10)(Last updated: April 5, 2019)