IEICE TRANS. FUNDAMENTALS, VOL.E100–A, NO.12 DECEMBER 2017

2925

PAPER

Gauss-Seidel HALS Algorithm for Nonnegative Matrix
Factorization with Sparseness and Smoothness Constraints
Takumi KIMURA†a) , Nonmember and Norikazu TAKAHASHI†b) , Senior Member

SUMMARY Nonnegative Matrix Factorization (NMF) with sparsenessand smoothness constraints has attracted increasing attention. When theseproperties are considered, NMF is usually formulated as an optimization
problem in which a linear combination of an approximation error term andsome regularization terms must be minimized under the constraint that thefactor matrices are nonnegative. In this paper, we focus our attention on theerror measure based on the Euclidean distance and propose a new iterative
method for solving those optimization problems. The proposed method
is based on the Hierarchical Alternating Least Squares (HALS) algorithm
developed by Cichocki et al. We first present an example to show thatthe original HALS algorithm can increase the objective value. We thenpropose a new algorithm called the Gauss-Seidel HALS algorithm that
decreases the objective value monotonically. We also prove that it has theglobal convergence property in the sense of Zangwill. We finally verify
the eﬀectiveness of the proposed algorithm through numerical experiments
using synthetic and real data.
key words: nonnegative matrix factorization, hierarchical alternating leastsquares algorithm, Euclidean distance, global convergence, sparseness,smoothness

1.

Introduction

Nonnegative Matrix Factorization (NMF) [1]–[3] is an operation that decomposes a given M × N nonnegative matrix
X = [X mn ] into an M × K nonnegative matrix W = [Wmk ]and a K × N nonnegative matrix HT = [Hnk ]T (see Fig. 1).
Because NMF is useful for extraction of nonnegative bases
and dimensionality reduction, it has found many applications
in various fields such as face image processing [2], [4], text
mining [5], recommender systems [6], [7], and so on.
NMF is formulated as a constrained optimization problem in which an error between X and WHT must be minimized under the constraint that all entries of W and H are
nonnegative. The Euclidean distance and various types ofdivergences have been used as the error criterion. The Euclidean distance-based NMF is formulated as the optimization problem:
2
1
X − WHTF
2subject to W ≥ 0 M×K , H ≥ 0 N ×K

minimize

where ∥·∥ F denotes the Frobenius norm, that is,
X − WHT

2
F

=

M ∑
N (
∑

X mn − (WHT )mn)2,

m=1 n=1

Manuscript received July 5, 2017.

† The authors are with Okayama University, Okayama-shi, 700-

8530 Japan.
a) E-mail: kimu@momo.cs.okayama-u.ac.jpb) E-mail: takahashi@cs.okayama-u.ac.jp
DOI: 10.1587/transfun.E100.A.2925

Fig. 1

Nonnegative matrix factorization.

0 M×K (0 N ×K , resp.) is the M × K (N × K, resp.) matrix ofall zeros and the inequality holds componentwise. In general, it is hard to find a global optimal solution of the NMF
optimization problem because the objective function is notconvex. So our goal is to find a local optimal solution.
The most popular algorithm for solving NMF optimization problems is the multiplicative update rules (MURs) developed by Lee and Seung [3]. They considered the Euclidean distance and the I-divergence as the error measureand derived update rules based on the idea of minimizing
a strictly convex function called the auxiliary function instead of the objective function itself. Later on, this idea
was generalized and applied to various types of error measures (see [8] for example). However, because the MURs
are expressed in the form of a fraction, they are not definedfor all pairs of nonnegative matrices W and H. Also, theglobal convergence is not guaranteed because of this problem. By the global convergence, we mean that any sequenceof solutions has at least one convergent subsequence and thelimit of any convergent subsequence is a stationary pointof the corresponding optimization problem [9]. In order toavoid the problem mentioned above, Gillis and Glineur [10]devised a modified version of the MUR for the Euclidean
distance by using the idea of Cichocki et al. [11], which
prevents variables from being less than a small positive constant. Furthermore, it was proved by Takahashi et al. that
this modification guarantees the global convergence of many
MURs [12]–[14].
The MURs have many good properties. They are simple and thus easy to implement. They are applicable tovarious types of error measures. Also, they have the global
convergence property as mentioned above. However, the
MURs are very slow in general. Therefore, during the last
decades, many authors have developed faster algorithms for
NMF [11], [15]–[17] which require less number of iterations
than the MURs. Among them, the hierarchical alternating
least squares (HALS) algorithm proposed by Cichocki etal. [11] is widely known as a simple and fast method for the

Copyright © 2017 The Institute of Electronics, Information and Communication Engineers

IEICE TRANS. FUNDAMENTALS, VOL.E100–A, NO.12 DECEMBER 2017

2926

Euclidean distance based NMF. In this algorithm, W and H
are partitioned into 2K blocks as W = [w1, w2, . . . , wK ] andH = [h1, h2, . . . , hK ], and these 2K blocks are updated one
by one in a cyclic manner. A similar method is rank-one
residue iteration method proposed by Ho [16]. Also, Kim etal. [18] recently showed that some algorithms including the
HALS algorithm can be derived using one common framework of the block coordinate descent method. Because theupdate rule of the HALS algorithm is expressed in the formof a fraction as in the case of the MURs, they are not definedfor all pairs of nonnegative matrices W and H. In order tosolve this problem, Cichocki et al. [11] proposed a modified
HALS algorithm, and recently this algorithm was proved tohave the global convergence property [19].
In many applications of NMF, it is preferable that theobtained factor matrices are sparse or/and smooth. A matrix
is said to be sparse if it has a small number of nonzero entries,
while it is said to be smooth if neighboring entries takesimilar values. A simple way to control the sparseness andthe smoothness is to add regularization terms representing
the L 1 norm and the Frobenius norm of W and H to theobjective function [20]–[24]. In this paper, we focus on
the NMF optimization problem considered by Cichocki etal. [23], [24] which is described asminimize f (W, H)
subject to W ≥ 0 M×K , H ≥ 0 N ×K

(1)

to show that the original HALS algorithm can increase theobjective value. We then propose a new algorithm called theGauss-Seidel HALS (GSHALS) algorithm with which theobjective value decreases monotonically, and prove that it
has the global convergence property in the sense mentioned
above. Finally, we verify the eﬀectiveness of the proposed algorithm through numerical experiments using synthetic andreal data.
Recently, various constraints have been considered forNMF (see [25] and references therein). Although we focus our attention only on the sparseness and smoothness
constraints mentioned above, the idea behind the GSHALS
algorithm may be useful for some other constraints. For example, the algorithm proposed by Liao and Zhang [26] forthe graph regularized NMF [27], which has the same problem as the HALS algorithm, can be easily modified by using
the same idea to guarantee the global convergence.
2. Nonnegative Matrix Factorization with Sparsenessand Smoothness Constraints
2.1 HALS Algorithm
Let us partition matrices W and H into 2K blocks as W =[w1, w2, . . . , wK ] and H = [h1, h2, . . . , hK ]. Then (1) can be
rewritten as follows:

where f (W, H) is given byf (W, H) =minimize2
αsm
1
X − WHT + αsp ∥H∥ 1 +∥LH∥ 2F ,F
2
2∥·∥ 1 denotes L 1 norm, that is,∥H∥ 1 =

2
K
∑
1
T
X−
wk hk
2
k=1
F
K
K
∑
αsm ∑
+αsp
∥hk ∥ 1 +∥Lhk ∥ 222
k=1
k=1
wk ≥ 0 M×1, hk ≥ 0 N ×1,
k = 1, 2, . . . , K .

N ∑
K
∑

|Hnk |.

n=1 k=1

Also, αsp > 0 and αsm > 0 are regularization parameters
controlling the levels of sparseness and smoothness, respectively. As examples of L, Cichocki et al. [24] presented
1 −1 0 · · · 0 
0 1 −1 ,
L1 =  .
..
..

.
.
 ..

0
1 −1
−1 2 −1 0 · · · 0 
0 −1 2 −1
L2 =  .
..
..
..

.
.
.
.
 .

 0
−1 2 −1
where L1 is a (N−1)×N matrix and L2 is a (N−2)×N matrix.
However, L is not restricted to these specific matrices, but
can be any T × N real matrix with 1 ≤ T ≤ N.
The objective of this paper is to develop a global convergence guaranteed algorithm for solving (1) based on
the HALS algorithm. We first introduce the HALS algorithm [23], [24] to solve (1). We then present an example

subject to

The HALS algorithm [23], [24] tries to minimize theobjective value by updating 2K blocks w1, w2, . . . , wK ,h1, h2, . . . , hK one by one in a cyclic manner. Typical update
orders are
w1 → h1 → w2 → h2 → · · · → wK → hK

(2)

and
w1 → w2 → · · · → wK → h1 → h2 → · · · → hK . (3)When updating wk , other 2K − 1 blocks are considered as
constants and the optimization problem:2
1 T
R − hk wTkF
2 kminimizeϕk (wk ) =subject to

wk ≥ 0 M×1(4)

is solved, where
Rk = X −K∑

wk̃ hTk̃ .

k̃=1, k̃,kSimilarly, when updating hk , other 2K − 1 blocks are considered as constants and the optimization problem:KIMURA and TAKAHASHI: GSHALS ALGORITHM FOR NMF WITH SPARSENESS AND SMOOTHNESS CONSTRAINTS

2927

minimize ψk (hk )subject to hk ≥ 0 N ×1

(5)

is solved, where ψk (hk ) is given byψk (hk ) =2
αsm
1
Rk − wk hTk +αsp ∥hk ∥ 1 +∥Lhk ∥ 22 .
F
2
2

As shown in [18], the problem (4) has a unique optimal
solution given by wk = [Rk hk ]+ /(hTk hk ) if hk , 0 N ×1 ,where [v]+ denotes the vector obtained from v by replacing
all negative entries with zero. As for the problem (5), wecannot obtain the explicit formula for the optimal solution,
but it is easy to find the minimum point of ψk (hk ) by solving∇ψk (hk ) = 0 N ×1 , that is,(wTk wk I N ×N + αsm LT L)hk = RTk wk − αsp 1 N ×1,(6)

where I N ×N is the N × N identity matrix and 1 N ×1 is the
N-dimensional column vector of all ones. Suppose that wk ,0 N ×1 . Then wTk wk I N ×N +αsm LT L is nonsingular and hence,
multiplying both sides of (6) by (wTk wk I N ×N + αsm LT L) −1
from the left, we have
hk = (wTk wk I N ×N +αsm LT L) −1 (RTk wk −αsp 1 N ×1 ).

(7)

The right-hand side may have one or more negative entries,
but we can obtain a nonnegative vector from it by applying
the operator [·]+ .
The HALS algorithm is based on the above idea, anddescribed by the following update rule:
wk ←[Rk hk ]+,hTk hk(8)hk ← [(wTk wk I N ×N + αsm LT L) −1
× (RTk wk − αsp 1 N ×1 )]+ .

(9)Let the solution after l (≥ 0) rounds of updates be denoted by (W(l), H(l) ). When using update rule (8) and (9),it could occur that one or more columns of W(l) and H(l)become zero for some l. If wk = 0 M×1 and L = L2 for
instance, wTk wk I N ×N + αsm LT L in the right-hand side of (9)is not invertible. If hk = 0 N ×1 , the denominator of the righthand side of (8) becomes zero. In these cases, the algorithm
has to be stopped before the solution is obtained. In order toavoid this situation, Cichocki et al. [24] used the update rule
expressed by
Rk hk wk ←  T  , hk hk 
ϵ+

(10)

hk ← [(wTk wk I N ×N + αsm LT L) −1
× (RTk wk − αsp 1 N ×1 )]ϵ +Fig. 2Situation where ψk (h k ) is increased by the HALS algorithm.

This notation is also used for scalars and matrices in later
discussion. When the update rule given by (10) and (11) isused, we have to consider a modified optimization problem:
minimize f (W, H)
subject to W ≥ ϵ1 M×K , H ≥ ϵ1 N ×K ,(12)

where 1 M×K (1 N ×K ) denotes the M × K (N × K) matrix ofall ones. In the following sections, by the HALS algorithm,
we mean the update rule given by (10) and (11).
2.2 Problem of HALS Algorithm
If wk ≥ ϵ1 M×1 then the right-hand side of (7) is the minimum
point of ψk (hk ). However, it is not always true that the righthand side of (11) is the optimal solution of the optimization
problem:minimize ψk (hk )subject to hk ≥ ϵ1 N ×1 .

(13)

To make matters worse, it may occur that the value of ψk isincreased by the update (11). To see this, let us consider thesituation shown in Fig. 2 where hk(0) represents the current
solution and h∗k is the minimum point of ψk (hk ). Because
the second entry of h∗k is less than ϵ, (11) returns hk(1) as the
new solution. However, looking at contours of ψk , we seethat ψk (hk(1) ) is greater than ψk (hk(0) ).
We now give a more concrete example to show that the
update (11) can increase the value of ψk . Let M = 1, N = 3,
K = 1, X = [3, 2, 1], L = [−1, 2, −1], αsp = 3/2, αsm = 1,
ϵ = 1, w1(0) = [1] and h1(0) = [6/5, 1, 1]T . Under this setting,
we have ψk (h1(0) ) = 6.94. If the value of h1 is updated by(11), we have h1(1) = [3/2, 1, 1]T and ψk (h1(1) ) = 7, which isgreater than ψk (h01 ).
3. GSHALS Algorithm(11)

instead of (8) and (9), where [v]ϵ + denotes the vector obtained from v by replacing all entries less than ϵ with ϵ.

3.1 Derivation of Update Rule
In order to solve the problem of the HALS algorithm pointed

IEICE TRANS. FUNDAMENTALS, VOL.E100–A, NO.12 DECEMBER 20172928

out in the previous section, we propose to update entries of hk
one by one instead of using (11). In the following discussion,
N entries of hk are denoted by hk1, hk2 . . . , hk N instead ofH1k , H2k , . . . , H N k . In addition, let the n-th columns of X
and Rk be denoted by xn and rkn . Here, rkn is given byK∑

rkn = xn −wk̃ hk̃n .

k̃=1, k̃,kNote that we do not focus on some specific update orderof entries of hk , but only assume that the update order isfixed during the execution of the proposed algorithm. When
updating hkn , other N − 1 entries are considered as constants
and the optimization problem:
minimizesubject toψkn (hkn )hkn ≥ ϵ(14)

is solved, whereψkn (hkn ) =1∥rkn − wk hkn ∥ 222
+ αsp hkn +TN)2
αsm ∑ (∑
L t ñ hk ñ .
2 t=1 ñ=1

∗ be the minimum point of ψ (h ). If h∗ is greater
Let hkn
kn kn
kn
than or equal to ϵ then it is the optimal solution of (14)because ψkn (hkn ) is strictly convex (see Fig. 3(a)). Otherwise, the optimal solution of (14) is ϵ (see Fig. 3(b)). Theseobservations can be summarized in the following lemma.

Lemma 1: If wk is positive then the optimization problem(14) has a unique optimal solution given by

Fig. 3 The minimum point h k∗ n of ψk n (h k n ). (a) The case where h k∗ n ≥ϵ . (b) The case where h k∗ n < ϵ .

T
N


∑
∑
 rTkn wk − αsp − αsm
L tn
L t ñ hk ñ 
t=1
ñ=1, ñ,n .
hkn = 
T
∑


2
T
wk wk + αsm
L tn



ϵ+
t=1
(15)

that minimizes ψkn (hkn ) because ψkn (hkn ) is strictly convex. If this is greater than or equal to ϵ then it is a uniqueoptimal solution of (14). Otherwise, ϵ is a unique optimal
solution of (14) because ψkn (hkn ) is strictly monotone increasing in [ϵ, ∞). Therefore, the optimal solution of (14) is
given by (15).
□
From Lemma 1, we obtain the update rule for entries ofhk , which is expressed as

Proof : The solution of the equation′
ψkn(hkn ) = (wTk wk )hkn − rTkn wk+ αsp + αsm

T
∑

L tn

N
∑

L t ñ hk ñ = 0

t=1

ñ=1

T
∑

N
∑

is given byrTkn wk − αsp − αsm
hkn =

L tn

t=1

wTk wk + αsmñ=1, ñ,n
T
∑
t=1

L 2tn

L t ñ hk ñ

T
N


∑
∑
 rTkn wk − αsp − αsm
L tn
L t ñ hk ñ 
t=1
ñ=1, ñ,n .
hkn ← 
T
∑


2
wTk wk + αsm
L tn



ϵ+
t=1
(16)

It is clear from the above discussion that the objective valueof (12) is not increased by the update (16).
Note that entries of hk are updated one by one by (16)like the Gauss-Seidel method (see [28] for example) for solving linear equations. We thus call this algorithm the GaussSeidel HALS (GSHALS) algorithm.

KIMURA and TAKAHASHI: GSHALS ALGORITHM FOR NMF WITH SPARSENESS AND SMOOTHNESS CONSTRAINTS

2929

3.2

Global Convergence of GSHALS Algorithm

In this section, we prove that the GSHALS algorithm has theglobal convergence property in the sense of Zangwill [9].
Let the feasible region and the set of stationary points of (12)be denoted by Fϵ and Sϵ , respectively. A point (W, H) ∈Fϵ is called a stationary point if it satisfies the following
conditions:∇W f (W, H) ≥ 0 M×K ,∇H f (W, H) ≥ 0 N ×K ,∇W f (W, H) ⊙ (ϵ1 M×K − W) = 0 M×K ,∇H f (W, H) ⊙ (ϵ1 N ×K − H) = 0 N ×K ,

(17)
(18)
(19)
(20)

It is clear from the update rule given by (10) and (16)that the third statement is true. Therefore, we prove that theremaining two statements hold.
First, we prove the second statement by the following
two lemmas. The first lemma shows the relation between Sϵand the optimal solutions of the subproblems. The second
one proves that f is strictly decreased by the mapping A if
the current solution is not a stationary point.
Lemma 2: (W∗, H∗ ) = ((w∗1, . . . , w∗K ), (h∗1, . . . , h∗K )) ∈ Fϵ
is a stationary point of (12) if and only if w∗k is a uniqueoptimal solution of the optimization problem:
minimizesubject toϕk (wk )wk ≥ ϵ1 M×1(21)∗ is a unique optimal solution offor k = 1, 2, . . . , K and hkn
the optimization problem:

where∇W f (W, H) =(WHT − X)H,
∇H f (W, H) =(HWT − XT )W+ αsp 1 N ×K + αsm LT LH
and ⊙ represents componentwise multiplication.
In order to make discussions simple, we define one
round of updates of 2K blocks by using (10) and (16) as a
mapping A : Fϵ → Fϵ . Then the solution after l rounds of
updates is expressed as(W(l), H(l) ) = A(W(l−1), H(l−1) )= A2 (W(l−2), H(l−2) )..
.
= Al (W(0), H(0) ).
Note that we do not focus on some specific update order suchas (2) and (3), but only assume that the update order is fixedduring the execution of the algorithm.
The global convergence property of the GSHALS algorithm is stated as follows.
Theorem 1: For any positive constant ϵ and initial solution∞ generated(W(0), H(0) ) ∈ Fϵ , the sequence {(W(l), H(l) )}l=0
by the GSHALS algorithm expressed by (10) and (16) has
at least one convergent subsequence and the limit of any
convergent subsequence belongs to Sϵ .
In the rest of this section, we prove Theorem 1 by using
Zangwill’s global convergence theorem [9]. Namely, weshow that the following statements hold true.
1. For any initial solution (W(0), H(0) ) ∈ Fϵ , the sequence∞ = { Al (W (0), H (0) )} ∞ is contained in{(W(l), H(l) )}l=0
l=0
a closed bounded subset of Fϵ .
2. The mapping A and the objective function f satisfy thefollowing statements.
a. If (W, H) ∈ Fϵ \Sϵ then f ( A(W, H)) < f (W, H).
b. If (W, H) ∈ Sϵ then f ( A(W, H)) ≤ f (W, H).
3. The mapping A is continuous in Fϵ \ Sϵ .

minimizesubject toψkn (hkn )hkn ≥ ϵ(22)

for k = 1, 2, . . . , K and n = 1, 2, . . . , N.
Proof : The Lagrangian function for the optimization problem (21) is given by
L(wk , λ) = ϕk (wk ) + λ T (ϵ1 M×1 − wk )where λ = (λ 1, . . . , λ M )T is the Lagrange multiplier vector.
Then w∗k ≥ ϵ1 M×1 is a stationary point of (21) if and only ifthere exists a λ that satisfies the following conditions:∂L(w∗ , λ) = 0 M×1,
∂wk k
λ ⊙ (ϵ1 M×1 − w∗k ) = 0 M×1,
λ ≥ 0 M×1 .

(23)
(24)
(25)

Here, the left-hand side of (23) is given by∑∂L(w∗k , λ) = *.w∗k h∗k T +w∗k̃ h∗k̃ T − X+/ h∗k − λ
∂wk
k̃,k
,
∗ ∗T
∗
= (W H − X)hk − λ.
So the conditions (23)–(25) can be rewritten as(W∗ H∗T − X)h∗k ≥ 0 M×1,)(
(W∗ H∗T − X)h∗k ⊙ (ϵ1 M×1 − w∗k ) = 0 M×1 .
Note that the problem (21) has a unique stationary point andit is also a unique optimal solution because ϕk (wk ) is strictly
convex. Therefore, w∗k ≥ ϵ1 M×1 is a unique optimal solutionof (21) for k = 1, 2, . . . , K if and only if(W∗ H∗T − X)H∗ ≥ 0 M×K ,(26))(W∗ H∗T − X)H∗ ⊙ (ϵ1 M×K − W∗ ) = 0 M×K .

(27)

(

IEICE TRANS. FUNDAMENTALS, VOL.E100–A, NO.12 DECEMBER 2017

2930

The Lagrangian function for the optimization problem(22) is given byL(hkn, µ) = ψkn (hkn ) + µ(ϵ − hkn )∗ ≥ ϵ is awhere µ is the Lagrange multiplier. Then hkn
stationary point of (22) if and only if there exists a µ thatsatisfies the following conditions:∂L(h∗ , µ) = 0,∂hkn kn∗µ ⊙ (ϵ − hkn) = 0,µ ≥ 0.

(28)
(29)
(30)

Here, the left-hand side of (28) is given by∂L(h∗ , µ)∂hkn kn

ñ=1

( ∑
K
*.
h∗pn w∗p T − xTn +/ w∗k + αsp
, p=1
)
T
N
∑
∑
∗
∗
+ αsm
L tn
L t ñ hk ñ (ϵ − hkn) = 0.
ñ=1Note that the problem (22) has a unique stationary point
and it is also a unique optimal solution because ψkn (hkn )is strictly convex. Therefore, hkn ≥ ϵ is a unique optimal
solution of the problem (22) for n = 1, 2, . . . , N if and only
if
(H∗ W∗T − XT )w∗k + αsp 1 N ×1(

k = 1, 2, . . . , K,

(H∗ W∗T − XT )w∗k + αsp 1 N ×1 + αsm LT Lh∗k

+ αsm LT LH∗ ≥ 0 N ×K ,((H∗ W∗T − XT )W∗ + αsp 1 N ×K + αsm LT LH∗

(31)

)

⊙ (ϵ1 N ×K − H∗ ) = 0 N ×K .

(32)

Note that the set of conditions (26)–(27) and (31)–(32) isequivalent to set of conditions (17)–(20) with W = W∗ and
H = H∗ , which is the necessary and suﬃcient condition for(W∗, H∗ ) ∈ Fϵ to be a stationary point of (12).
□holds for k = 1, 2, . . . , K andK∑*.
h∗pn w∗p T − xTn +/ w∗k + αsp
, p=1
T
N
∑
∑
+ αsm
L tn
L t ñ hk∗ ñ ≥ 0,

+ αsm LT Lh∗k ≥ 0 N ×1,(H∗ W∗T − XT )W∗ + αsp 1 N ×K
Rk hk wk =  T  hk hk ϵ+So the conditions (28)–(30) can be rewritten as

t=1

These equations can be rewritten as

Proof : Suppose first that (W, H) ∈ Sϵ . Then it follows
from Lemma 2 that wk and hkn are unique optimal solutionsof (21) and (22), respectively, for k = 1, 2, . . . , K and n =1, 2, . . . , N. In other words,

ñ=1

t=1

k = 1, 2, . . . , K .

Lemma 3: If (W, H) ∈ Sϵ then A(W, H) = (W, H). If(W, H) ∈ Fϵ \ Sϵ then f ( A(W, H)) < f (W, H).

∑∗∗= *. hknw∗k T +hk̃nw∗k̃ T − xTn +/ w∗k + αsp
k̃,k
,
T
∑
∑
∗
+ αsm
L tn * L tn hkn+
L t ñ hk∗ ñ + − µ
ñ,n
,
t=1
K
∑
h∗pn w∗p T − xTn +/ w∗k + αsp
= *.
, p=1
T
N
∑
∑
+ αsm
L tn
L t ñ hk∗ ñ − µ.
t=1

⊙ (ϵ1 N ×1 − h∗k ) = 0 N ×1,

)

T
N


∑
∑
 rTkn wk − αsp − αsm
L tn
L t ñ hk ñ 
t=1
ñ=1, ñ,n
hkn = 
T
∑


2
T
wk wk + αsm
L tn


ϵ+

t=1

holds for k = 1, 2, . . . , K and n = 1, 2, . . . , N. We thus
have A(W, H) = (W, H). Suppose next that (W, H) ∈Fϵ \ Sϵ . Let K1 ⊆ {1, 2, . . . , K } be the set of k suchthat wk is not an optimal solution of (21), and K2 ⊆{(1, 1), (1, 2), . . . , (1, N ), (2, 1), (2, 2), . . . , (K, N )} be the setof pairs (k, n) such that hkn is not an optimal solutionof (22). Then it is clear that K1 , ∅ or K2 , ∅holds. If wk̂ is the block first updated among all blocks
in {wk | k ∈ K1 } ∪ {hkn | (k, n) ∈ K2 }, the objective value
does not change before the update of wk̂ , strictly decreases
through the update of wk̂ , and does not increase after the
update of wk̂ . Similarly, if hk̂ n̂ is the block first updated
among all blocks in {wk | k ∈ K1 } ∪ {hkn | (k, n) ∈ K2 }, thevalue of the objective function does not change before theupdate of hk̂ n̂ , strictly decreases through the update of hk̂ n̂ ,and does not increase after the update of hk̂ n̂ . From theseobservations, we conclude that f ( A(W, H)) is strictly lessthan f (W, H).
□
Finally, we prove the first statement.
Lemma 4: For any W(0) ≥ ϵ1 M×K , H(0) ≥ ϵ1 N ×K , the se∞ generated by (10) and (16) is conquence {(W(l), H(l) )}l=0
tained in a closed bounded set

KIMURA and TAKAHASHI: GSHALS ALGORITHM FOR NMF WITH SPARSENESS AND SMOOTHNESS CONSTRAINTS

2931{(W, H) | f (W, H) ≤ f (W(0), H(0) ),W ≥ ϵ1 M×K , H ≥ ϵ1 N ×K }.

with a stopping condition is stated as follows.
(33)

∞ is
Proof : It is clear from Lemma 3 that {(W(l), H(l) )}l=0
contained in (33). So it suﬃces for us to prove that (33) isbounded. Let C = f (W(0), H(0) ). If f (W, H) ≤ C then we
have the following inequality:
2
1
X − WHT ≤ CF2Algorithm 1 GSHALS
Input: X ∈ R+M × N , L ∈ RT × N , K ∈ N, αsp, αsm, ϵ, δ1, δ2 ∈ R++
Output: (W, H) ∈ S̄ϵ1: Choose (W, H) ∈ Fϵ .
2: Update 2K columns of W and H one by one in the fixed order by using(10) and (16).
3: If (W, H) satisfies (37)–(40) then return (W, H) and stop. Otherwise
go to Step 2.

which implies thatK∑√√
X mn − 2C ≤Wmk Hnk ≤ X mn + 2C

(34)

k=1

for all m and n. Moreover, if (W, H) ∈ Fϵ , it follows from(34) that√√
X mn + 2C
X mn + 2C
Wmk ≤
≤
,
(35)
Hnk
ϵ
√
√
X mn + 2C
X mn + 2C
Hnk ≤
≤
(36)
Wmk
ϵ
for all m, n and k. This means that (33) is bounded.

□

In the above analysis, we took the same approach as in
the previous work [19], but the boundedness of solutions wasproved in a diﬀerent way. In [19], a compact subset of Fϵ ,which depends only on X and ϵ, was derived directly from the
update rule, while in this paper we proved the boundednessof solutions by making use of the level set of the objectivefunction which is determined not only from X and ϵ but also
the initial solution (W(0), H(0) ).
3.3

Theorem 2: For any positive constants ϵ, δ1, δ2 and any
initial solution (W(0), H(0) ) ∈ Fϵ , Algorithm 1 with the
stopping condition given by (37)–(40) always stops within afinite number of iterations.
3.4 Fast GSHALS Algorithm
In the HALS algorithm, the computational cost per roundis reduced by updating matrices E = X − WHT and Rk(k = 1, 2, . . . , K ) in a proper way [23]. This technique canbe directly applied to the GSHALS algorithm. An important
point is that Rk can be eﬃciently computed from the current
values of E, wk and hk as
Rk = E + wk hTk
and E can be eﬃciently computed from the current values of
Rk , wk and hk as
E = Rk − wk hTk .

Finite Termination of GSHALS Algorithm

Note that Theorem 1 does not guarantee the convergence of
the whole sequence but only the existence of a subsequence
converging to a stationary point. This is, however, suﬃcientfor us because, by introducing an appropriate stopping condition, we can obtain an algorithm that always stops within
a finite number of iterations after reaching an approximate
stationary point [12].
The necessary and suﬃcient condition described by(17)–(20) for (W, H) ∈ Fϵ to be a stationary point can berelaxed by using any positive constants δ1 , δ2 as follows:∇W f (W, H) ≥ −δ1 1 M×K ,∇H f (W, H) ≥ −δ1 1 N ×K ,Wmk − ϵ ≤ δ2 if (∇W f (W, H)) mk > δ1,
Hnk − ϵ ≤ δ2 if (∇H f (W, H)) nk > δ1 .

For Algorithm 1, we have the following theorem. The
proof is omitted because it is similar to [12].

(37)
(38)
(39)
(40)

Let the set of (W, H) ∈ Fϵ which satisfies (37)–(40) be
denoted by S¯ϵ . Also, let R, R+ and R++ denote the set ofreal numbers, the set of nonnegative real numbers, and theset of positive real numbers, respectively. Furthermore, let
N denote the set of natural numbers. Then the algorithm

Another important point is that ∇W f (W, H) and ∇H f (W, H)
can be eﬃciently computed from the current values of W,
H and E as ∇W f (W, H) = −EH and ∇H f (W, H) = −ET W.
Therefore, we can determine whether the relaxed optimality
condition given by (37)–(40) is satisfied or not, by checking
the following conditions:
EH ≤ δ1 1 M×K ,(41)

E W ≤ δ1 1 N ×K ,Wmk − ϵ ≤ δ2 if (EH)mk < −δ1,

(42)
(43)

Hnk − ϵ ≤ δ2 if (ET W)nk < −δ1 .

(44)

T

Making use of these ideas, we obtain two fast GSHALS
algorithms depending on the update order as follows.

IEICE TRANS. FUNDAMENTALS, VOL.E100–A, NO.12 DECEMBER 2017

2932

Algorithm 2 Fast GSHALS with update order (2)Input: X ∈ R+M ×N , L ∈ RT × N , K ∈ N, αsp, αsm, ϵ, δ1, δ2 ∈ R++
Output: (W, H) ∈ S̄ϵ1: Set M ← LT L.
2: Choose (W, H) ∈ Fϵ and set E ← X − WHT .
3: Set k ← 1.
4: Set R k ← E + w k hTk .
[]5: Set w k ← R k h k /hTk h k.
ϵ+
6: Set n ← 1.
 rT w k − αsp − αsm ∑ N
ñ=1, ñ,n Mn ñ h k ñ kn
7: Set h k n ← 
 .
T

w k w k + αsp Mnn
ϵ+8: If n = N then go to Step 9. Otherwise add 1 to n and go to Step 7.
T
9: Set E ← R k − w k hk .
10: If k = K then go to Step 11. Otherwise add 1 to k and go to Step 4.
11: If (W, H) satisfies (41)–(44) then return (W, H) and stop. Otherwise
go to Step 3.

4. Numerical Experiments

Algorithm 3 Fast GSHALS with update order (3)where ⊘ represents the componentwise division. The computational complexity per round of the MUR is O(M N K +N 2 K ) like the HALS and GSHALS algorithms.
In all experiments, we set L = L2 and use (37)–(40) as astopping condition. The HALS and GSHALS algorithms are
implemented by using the technique described in Sect. 3.4.
In the GSHALS algorithm, hk is updated in the order hk1 →hk2 → · · · → hk N . All methods are implemented in C
language with BLAS and LAPACK libraries, compiled with
gcc 5.3.0, and executed on a PC with Intel Core i7-6700,
16 GB memory and Windows 10.

Input: X ∈ R+M ×N , L ∈ RT × N , K ∈ N, αsp, αsm, ϵ, δ1, δ2 ∈ R++
Output: (W, H) ∈ S̄ϵ1: Set M ← LT L.
2: Choose (W, H) ∈ Fϵ and set E ← X − WHT .
3: Set k ← 1.
4: Set R k ← E + w k hTk .
[]5: Set w k ← R k h k /hTk h k.
ϵ+

6: Set E ← R k − w k hTk .
7: If k = K then go to Step 8. Otherwise add 1 to k and go to Step 4.
8: Set k ← 1.
9: Set R k ← E + w k hTk .
10: Set n ← 1.
 rT w k − αsp − αsm ∑ N
ñ=1, ñ,n Mn ñ h k ñ kn
11: Set h k n ←  .
wTk w k + αsp Mnn
ϵ+12: If n = N then go to Step 13. Otherwise add 1 to n and go to Step 11.
13: Set E ← R k − w k hTk .
14: If k = K then go to Step 15. Otherwise add 1 to k and go to Step 9.
15: If (W, H) satisfies (41)–(44) then return (W, H) and stop. Otherwise
go to Step 3.

It is easy to see that the computational complexity perround of Algorithms 2 and 3 is O(M N K + N 2 K ). This isequal to the computational complexity per round of the Fast
HALS algorithm because, using the eigenvalue decomposition LT L = QΛQT where Q is an orthogonal matrix and Λ
is a nonnegative diagonal matrix, we can rewrite (11) as
hk ← [Q(wTk wk I N ×N + αsm Λ) −1 QT
× (RTk wk − αsp 1 N ×1 )]ϵ +which takes O(M N + N 2 ) time. An important diﬀerence
between Algorithms 2 and 3 is that the former updates Rk
and E K times in each round, while the latter needs to do it
2K times. Therefore, the computation time of Algorithm 2
is shorter than that of Algorithm 3 for the same number ofrounds. However, we cannot conclude from this observation that Algorithm 2 is faster than Algorithm 3, because thenumber of rounds required to reach an approximate stationary point depends on the update order.

In order to examine the eﬀectiveness of the GSHALS algorithm, we compare the performance of the MUR, the HALS
algorithm and the GSHALS algorithm by using synthetic
and real datasets. The MUR can be easily obtained by using
the unified method proposed by Yang and Oja [8], and theresulting update rule is described by[](45)
W ← W ⊙ XH ⊘ WHT H ,ϵ+[TH← H⊙X W()]⊘ HWT W + αsp 1 N ×K + αsm LT LH
,
ϵ+
(46)

4.1 Experiment Using Synthetic Datasets
We first compare the performance of the five methods:
MUR, HALS with (2), HALS with (3), GSHALS with(2), and GSHALS with (3) by using synthetic datasets.
In this experiment, we set M = 100, N = 50, K = 10,
αsm = 0.1, αsp = 0.1, ϵ = 0.001, δ1 ∈ {0.1, 0.01, 0.001} andδ2 ∈ {0.01, 0.001, 0.0001}. For each of the nine pairs of δ1and δ2 values, we applied the five methods to 10 diﬀerent
triples (X, W(0), H(0) ) which were generated in such a waythat each entry of X, W(0) and H(0) was drawn from an independent uniform distribution on the interval [0, 1] and thenall entries of W(0) and H(0) less than ϵ were replaced with ϵso that (W(0), H(0) ) ∈ Fϵ . The performance of the five methods are compared in terms of the number of rounds and thecomputation time. The maximum number of rounds were
set to 60000, that is, each algorithm stops when the numberof rounds reaches 60000 even if the stopping condition is notsatisfied.
Results of the experiment are shown in Tables 1–3. We
first see from these tables that the computational cost is not
so sensitive to the value of δ2 for all algorithms. From
Tables 1 and 2 where the results for δ1 = 0.1 and δ1 = 0.01
are given respectively, we see that i) the HALS and the
GSHALS algorithms are much faster than the MUR, ii) the
update order (2) is faster than (3) for the same algorithm,
and iii) HALS with (2) is faster than GSHALS with (3).
However, some of these properties do not hold for a smaller
value of δ1 . In fact, we see from Table 3, where the results

KIMURA and TAKAHASHI: GSHALS ALGORITHM FOR NMF WITH SPARSENESS AND SMOOTHNESS CONSTRAINTS

2933

for δ1 = 0.001 are given, that the HALS algorithm did notsatisfy the stopping condition before the number of rounds
reached 60000 in all cases. This is because the objective
value eventually stops decreasing as explained in Sect. 2.2.
In contrast, the MUR and the GSHALS algorithm stopped
within 60000 rounds in all cases. The GSHALS algorithm ismuch faster than the MUR. In particular, the GSHALS with
Results for synthetic data (δ1 = 0.1).

Table 1
δ2

Method0.01

MUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)MUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)MUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)0.0010.0001

Method0.01

MUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)MUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)MUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)0.0001

Table 3
δ2

Number of Rounds
Ave.
Min.
Max.
15793.2 10982
24186
1229.1
694
2296
2664.8
2271
3702
962.8
342
1467
2788.22564
2927
16336.0 10403
33477
1051.0571
2233
2564.42295
2870
1163.6487
2512
2790.72632
2919
20129.8 10458
38751
1092.5571
1601
2631.52470
3015
1150.8641
1583
2945.8
2683
3848

0.001

0.0001MUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)MUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)MUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)Computation Time [s]
Ave.
Min.
Max.
6.018
4.104
9.168
0.374
0.204
0.687
1.083
0.922
1.548
0.313
0.103
0.473
1.211
1.118
1.307
6.141
3.642
12.894
0.314
0.172
0.672
1.030
0.922
1.156
0.377
0.152
0.779
1.207
1.110
1.312
7.268
3.657
14.1400.335 0.172
0.485
1.077
1.000
1.235
0.355 0.192
0.481
1.229
1.107
1.584

Results for synthetic data (δ1 = 0.001).

MethodAve.

0.01

Computation Time [s]
Ave.
Min.
Max.
2.191
1.032
3.891
0.113
0.063
0.157
0.433
0.375
0.469
0.121
0.063
0.235
0.512
0.484
0.672
2.704
1.719
3.985
0.155
0.093
0.235
0.464
0.422
0.515
0.107
0.063
0.172
0.478
0.438
0.516
2.977
1.708
5.678
0.131
0.078
0.235
0.469
0.375
0.609
0.096
0.037
0.160
0.507
0.451
0.670

Results for synthetic data (δ1 = 0.01).

Table 2δ20.001

Number of Rounds
Ave.
Min.
Max.
6173.3 2969
10985
375.9
210
539
1049.9
919
1178
379.5
204
739
1237.8 1126
1728
7618.8 4819
11179
495.8
319
752
1146.8
1030
1268
354.7
220
585
1165.51060
1214
8203.5 4888
15777
433.8
267
772
1155.5
925
1521
306.5
125
515
1212.4
1096
1577

Number of Rounds
Min.
Max.

27251.0
60000.0
60000.0
1633.04554.124850.6
60000.0
60000.02152.14577.627367.8
60000.0
60000.01195.54554.4

14184
60000
60000
852
4281
11136
60000
60000
1232
4426
14684
60000
60000
985
4444

46474
60000
60000
4352
4795
53797
60000
60000
3723
4798
41168
60000
60000
1995
4921

Table 4Computation Time [s]
Ave.
Min.
Max.
9.474
18.52224.867
0.526
1.930
8.676
18.49624.705
0.655
1.894
9.453
18.63624.939
0.370
1.930

4.647
18.17424.424
0.249
1.754
3.672
18.18924.471
0.371
1.838
5.047
18.36124.5340.296
1.812

18.34018.79925.2601.377
2.162
19.90818.69025.003
1.141
1.968
13.37718.95525.5190.609
2.172

4.2 Experiment Using Real Datasets
We next compare the performance of the three methods:
MUR, GSHALS with (2), and GSHALS with (3) by using
three kinds of real datasets. We do not consider the HALS
algorithm in this experiment because we have observed in theprevious experiment that the HALS algorithm always showsa similar or worse performance than the GSHALS algorithm
and, more importantly, the HALS algorithm does not always
reach an approximate stationary point.
The first dataset is the ORL Database of Faces† which isa facial image dataset oﬀered by AT&T Laboratories Cambridge. This dataset contains 400 grayscale facial images
and the size of each image is 92 × 112. Reducing the sizeof all images to 46 × 56, transforming them into column
vectors, we obtain a 2576 × 400 nonnegative matrix X. We
set αsm = 0.1, αsp = 0.1, ϵ = 1.0, δ1 = 10.0, δ2 = 1.0
and K = 5, and run the three algorithms for 10 diﬀerent
initial solutions (W(0), H(0) ) which are generated in such a
way that each entry is drawn from an independent uniform
distribution on the interval [0, 10] and then all entries lessthan ϵ are replaced with ϵ so that (W(0), H(0) ) ∈ Fϵ .
The second dataset is the Wisconsin Diagnostic Breast
Cancer (WDBC) dataset which can be obtained from the UCI
Machine Learning Repository†† . This dataset contains 569
instances each of which consists of 32 features. Excluding
the first two features (the first is ID number and the second
is diagnosis), normalizing each of the remaining 30 features
so that the values belong to the interval [0, 1], we obtain a
30×569 nonnegative matrix X. We set αsm = 0.1, αsp = 0.1,
ϵ = 0.001, δ1 = 0.005, δ2 = 0.001 and K = 2, and run thethree algorithms for 10 diﬀerent initial solutions (W(0), H(0) )which are generated in the same way as in the case of the† http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.

html†† http://archive.ics.uci.edu/ml/

Final objective values for synthetic data (δ1 = 0.1, δ2 = 0.0001).

MethodMUR
HALS (2)
HALS (3)
GSHALS (2)GSHALS (3)(2) is the fastest.
It should be noted that the objective value when the
stopping condition is satisfied diﬀers depending on the algorithm, even if the same parameter values and the same initialcondition are used. It often occurs that the final objective
value reached by an algorithm is less than that by a faster algorithm. In fact, we see from Table 4 that the final objective
value obtained by the GSHALS algorithm with the update
order (2), which is the fastest among five methods as shownin Table 1, is not the smallest for all settings. This meansthat diﬀerent algorithms may reach diﬀerent approximate
stationary points.

Setting 1

Setting 2

Setting 3

Setting 4

Final Objective Value
Setting 5
Setting 6

Setting 7

Setting 8

Setting 9

Setting 10127.268
127.285
127.107127.607127.080

124.773124.595
124.482
124.732124.599125.814126.137125.736
126.211125.854130.235129.923130.007130.167130.128

127.546
127.276
127.714127.889127.625127.147
127.428
127.449
128.165
127.182

126.615
126.627126.471126.610126.397124.520125.085
124.609125.574124.541130.190130.319130.028130.470130.049

126.461
126.836127.004
126.309126.216

IEICE TRANS. FUNDAMENTALS, VOL.E100–A, NO.12 DECEMBER 2017

2934
Table 5
MethodMUR
GSHALS (2)GSHALS (3)Ave.

7197.31509.03662.3466.31964.733222.477

Table 6
MethodMUR
GSHALS (2)GSHALS (3)MUR
GSHALS (2)GSHALS (3)

6366
1121
3264

8909
1587
3834

Computation Time [s]
Min.
Max.
420.269
46.488207.611560.45475.769237.239

Results for WDBC dataset.

Number of Rounds
Ave.
Min.
Max.
27032.4 24760
27801
14780.810457
17660
11109.1
8289
15332

Table 7
MethodResults for ORL dataset.

Number of Rounds
Ave.
Min.
Max.

Computation Time [s]
Ave.
Min.
Max.
20.92419.03323.36213.60810.37615.752
11.1238.500
15.987

Results for CLUTO (tr23) dataset.

Number of Rounds
Ave.
Min.
Max.

Ave.

3928.9830.0
784.5

306.357
40.00956.870

2063
742
764

6152
887
804

Computation Time [s]
Min.
Max.
162.23435.86554.368480.680
43.44261.082

ORL dataset except that the interval [0, 1] is used instead of[0, 10].
The third dataset is “tr23” in the CLUTO datasets. The
CLUTO datasets have been used in evaluating the performance of document clustering algorithms and can be obtained from the CLUTO web site† . The dataset “tr23” contains 204 instances corresponding to 204 documents from six
categories, and each instance contains 5832 features. Hencea 5832 × 204 nonnegative matrix X can be obtained fromthis dataset. We set αsm = 0.1, αsp = 0.1, ϵ = 0.1, δ1 = 1.0,
δ2 = 0.1 and K = 6, and run the three algorithms for 10
diﬀerent initial solutions (W(0), H(0) ) which are generated
in the same way as in the case of the ORL dataset.
Results of the experiment are shown in Tables 5–7. It isseen from these tables that the GSHALS algorithm is fasterthan the MUR for all datasets, as in the case of synthetic data.
As for the update order of the GSHALS algorithm, which isfaster depends on the dataset. The computation time of the
GSHALS algorithm with (2) is shorter than that with (3) for
the ORL and the CLUTO datasets, but the latter is faster for
the WDBC dataset.
5.

Conclusion

The HALS algorithm for NMF with sparseness and smoothness constraints has been studied in this paper. We haveproposed the GSHALS algorithm based on the HALS algorithm, and proved that it has the global convergence property.
We have also showed experimentally that the GSHALS algorithm outperforms the HALS algorithm and the MUR interms of computation time. The GSHALS algorithm can beeasily applied to the case where regularization terms for W
added to the objective function, though it is not consideredin this paper. A future problem is the global convergence
analysis of the GSHALS algorithm when the update order is
not fixed.
† http://glaros.dtc.umn.edu/gkhome/cluto/cluto/overview/

Acknowledgments
The authors would like to thank anonymous reviewers fortheir valuable comments to improve the quality of this paper.
This work was supported by JSPS KAKENHI Grant Number
JP15K00035.
References[1] P. Paatero and U. Tapper, “Positive matrix factorization: A nonnegative factor model with optimal utilization of error estimates ofdata values,” Environmetrics, vol.5, no.2, pp.111–126, June 1994.
DOI: 10.1002/env.3170050203[2] D.D. Lee and H.S. Seung, “Learning the parts of objects by nonnegative matrix factorization,” Nature, vol.401, pp.788–792, 1999.
DOI: 10.1038/44565[3] D.D. Lee and H.S. Seung, “Algorithms for non-negative matrix factorization,” in Advances in Neural Information Processing Systems,
T.K. Leen, T.G. Dietterich, and V. Tresp, Eds., vol.13, pp.556–562,
2001.
[4] P.O. Hoyer, “Non-negative matrix factorization with sparseness constraints,” J. Machine Learning Research, vol.5, pp.1457–1469, 2004.
[5] V.P. Pauca, F. Shahnaz, M.W. Berry, and R.J. Plemmons, “Textmining using non-negative matrix factorizations,” Proc. 4th SIAM
International Conference on Data Mining, pp.452–456, 2004. DOI:
10.1137/1.9781611972740.45[6] S. Zhang, W. Wang, J. Ford, and F. Makedon, “Learning from incomplete ratings using non-negative matrix factorization,” Proc. 6th
SIAM International Conference on Data Mining, pp.548–552, 2006.
DOI: 10.1137/1.9781611972764.58[7] X. Luo, M. Zhou, Y. Xia, and Q. Zhu, “An eﬃcient non-negative
matrix-factorization-based approach to collaborative filtering forrecommender systems,” IEEE Trans. Ind. Informat., vol.10, no.2,
pp.1273–1284, 2014. DOI: 10.1109/TII.2014.2308433[8] Z. Yang and E. Oja, “Unified development of multiplicative algorithms for linear and quadratic nonnegative matrix factorization,”
IEEE Trans. Neural Netw., vol.22, no.12, pp.1878–1891, Dec. 2011.
DOI: 10.1109/TNN.2011.2170094[9] W.I. Zangwill, Nonlinear Programming: A Unified Approach,Prentice-Hall, 1969.
[10] N. Gillis and F. Glineur, “Nonnegative matrix factorization and the
maximum edge biclique problem,” arXiv0810.4225, 2008.
[11] A. Cichocki, R. Zdunek, and S. Amari, “Hierarchical ALS algorithms for nonnegative matrix and 3D tensor factorization,” Springer
Lecture Notes in Computer Science, vol.4666, pp.169–176, 2007.
DOI: 10.1007/978-3-540-74494-8_22[12] N. Takahashi and R. Hibi, “Global convergence of modified multiplicative updates for nonnegative matrix factorization,” Computational Optimization and Applications, vol.57, no.2, pp.417–440,
2014. DOI: 10.1007/s10589-013-9593-0[13] N. Takahashi, J. Katayama, and J. Takeuchi, “A generalized suﬃcientcondition for global convergence of modified multiplicative updates
for NMF,” Proc. 2014 International Symposium on Nonlinear Theory
and its Applications, pp.44–47, 2014.
[14] N. Takahashi and M. Seki, “Multiplicative update for a class ofconstrained optimization problems related to NMF and its global
convergence,” Proc. 2016 European Signal Processing Conference,
pp.438–442, 2016. DOI: 10.1109/EUSIPCO.2016.7760286[15] C.-J. Lin, “Projected gradient methods for nonnegative matrix factorization,” Neural Computation, vol.19, no.10, pp.2756–2779, 2007.
DOI: 10.1162/neco.2007.19.10.2756[16] N.D. Ho, Nonnegative Matrix Factorization Algorithms and Applications, Ph.D thesis, Univ. Catholique de Louvain, 2008.
[17] N. Guan, D. Tao, Z. Luo, and B. Yuan, “NeNMF: An optimal gradient method for nonnegative matrix factorization,” IEEEKIMURA and TAKAHASHI: GSHALS ALGORITHM FOR NMF WITH SPARSENESS AND SMOOTHNESS CONSTRAINTS

2935

Trans. Signal Process., vol.60, no.6, pp.2882–2898, 2012. DOI:
10.1109/TSP.2012.2190406[18] J. Kim, Y. He, and H. Park, “Algorithms for nonnegative matrix
and tensor factorization: A unified view based on block coordinate
descent framework,” J. Global Optimization, vol.58, no.2, pp.285–
319, 2014. DOI: 0.1007/s10898-013-0035-4[19] T. Kimura and N. Takahashi, “Global convergence of a modified HALS algorithm for nonnegative matrix factorization,” Proc.
2015 IEEE 6th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, pp.21–24, 2015. DOI:
10.1109/CAMSAP.2015.7383726[20] P.O. Hoyer, “Non-negative sparse coding,” Proc. 2002 IEEE 12th
International Workshop on Neural Networks for Signal Processing,
pp.557–565, 2002. DOI: 10.1109/NNSP.2002.1030067[21] M.W. Berry, M. Browne, A.N. Langville, V.P. Pauca, and R.J. Plemmons, “Algorithms and applications for approximate nonnegative
matrix factorization,” Computational Statistics & Data Analysis,
vol.52, no.1, pp.155–173, 2007. DOI: 10.1016/j.csda.2006.11.006[22] A. Cichocki, S. Amari, R. Zdunek, R. Kompass, G. Hori, and Z.
He, “Extended SMART algorithms for non-negative matrix factorization,” Lecture Notes in Computer Science, Springer, vol.4029,
pp.548–562, 2006. DOI:10.1007/11785231_58[23] A. Cichocki and A.H. Phan, “Fast local algorithms for large scalenonnegative matrix and tensor factorization,” IEICE Trans. Fundamentals, vol.E92-A, no.3, pp.708–721, March 2009.
[24] A. Cichocki, R. Zdunek, A.H. Phan, and S. Amari, Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory MultiWay Data Analysis and Blind Source Separation, John Wiley & Sons,
2009.
[25] Y. Wang and Y. Zhang, “Nonnegative matrix factorization: A comprehensive review,” IEEE Trans. Knowl. Data Eng., vol.25, no.6,
pp.1336–1353, June 2013. DOI: 10.1109/TKDE.2012.51[26] Q. Liao and Q. Zhang, “Eﬃcient rank-one residue approximation
method for graph regularized non-negative matrix factorization,”
Proc. European Conference on Machine Learning and Principles andPractice of Knowledge Discovery in Databases, Part II, pp.242–255,
2013. DOI: 10.1007/978-3-642-40991-2_16[27] D. Cai, X. He, J. Han, and T.S. Huang, “Graph regularized nonnegative matrix factorization for data representation,” IEEE Trans.
Pattern Analysis and Machine Intelligence, vol.33, no.8, pp.1548–
1560, Aug. 2011. DOI: 10.1109/TPAMI.2010.231[28] O. Axelsson, Iterative Solution Methods, Cambridge University
Press, 1996.

Takumi Kimura
received the B.E. degree
in information technology and M.E. degree inelectronic and information systems engineeringfrom Okayama University, Japan in 2015 and2017, respectively. He is now with ASTEC Co.,
Ltd., Osaka, Japan. His research interests include nonnegative matrix factorization and its
applications.

Norikazu Takahashi
received the B.E.,
M.E. and D.E. degrees from Kyushu University, Japan, in 1991, 1993 and 1996, respectively. He is currently a Professor in the Department of Computer Science, Okayama University,
Japan. His research interests include optimization theory, nonlinear systems, multiagent systems, graph theory, and neural networks. He is amember of IEEE and Japanese Neural Network
Society.