CONVERGENCE-GUARANTEED MULTIPLICATIVE ALGORITHMS FORNONNEGATIVE MATRIX FACTORIZATION WITH β-DIVERGENCE
Masahiro Nakano† , Hirokazu Kameoka‡ , Jonathan Le Roux‡ , Yu Kitano† , Nobutaka Ono† , Shigeki Sagayama††

Graduate School of Information Science and Technology, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan‡NTT Communication Science Laboratories, NTT Corporation,
3-1 Morinosato Wakamiya, Atsugi, Kanagawa 243-0198, Japan

ABSTRACTThis paper presents a new multiplicative algorithm for nonnegative matrix factorization with β-divergence. The derived update rules have a similar form to those of the conventional multiplicative algorithm, only differing throughthe presence of an exponent term depending on β. The convergence is theoretically proven for any real-valued β basedon the auxiliary function method. The convergence speed
is experimentally investigated in comparison with previous
works.
1. INTRODUCTION
Nonnegative matrix factorization (NMF) [1] has recently
become a very popular technique in signal processing, andhas been succesfully applied to various problems such assource separation [2, 3, 4], feature extraction, music transcription [5] or dimension reduction. Given a nonnegative
matrix Y, the goal of NMF is to find two nonnegative matrices H and U such that Y ≈ HU. To measure how close
Y and HU are, the Euclidean (EUC) distance, the generalized Kullback-Leibler (KL) divergence and the ItakuraSaito (IS) divergence are often chosen. The three of them
are enclosed in the more general framework of β-divergence[6, 7]. Since the choice of an appropriate divergence depends on the application and the data [2, 8, 9], an algorithm
stable for a wide range of β is desired.
Multiplicative gradient descent [7, 10] is one of the mostpopular approaches for NMF with β-divergence. A proof ofthe convergence of the algorithms for β = 2 (EUC distance)and β = 1 (KL divergence) was given in [10], and it has
been extended to 1 ≤ β ≤ 2 [11]. However, convergencefor β < 1 and β > 2 remains to be proven. A generalized
multiplicative algorithm, which introduces an exponent stepsize, has also recently been proposed in [12], discussing inparticular the local and stable convergence, in the sense of
Lyapunov’s theory, of this algorithm when initialized in agiven neighborhood of a local minimum. However, conver-gence is not guaranteed in general.
Another way to derive optimization algorithms for NMF
is through a statistical approach. The minimization of a βdivergence can indeed be shown to be equivalent, for specific βs, to a Maximum-Likelihood (ML) estimation problem, due to the reproductive properties of some probabilistic distributions. Update equations for NMF with EUC distance (β = 2), KL divergence (β = 1) and IS divergence(β = 0) have been obtained under this framework in [13]based on the expectation maximization (EM) algorithm. Although convergence to a stationary point is then guaranteed,this approach is currently limited to the cases β = 0, 1, and2.
This paper proposes a new multiplicative algorithm for
NMF with β-divergence, for which the monotonic decreaseof the objective function at each iteration is theoreticallyguaranteed for all β. The derivation of this algorithm isbased on the careful design of a so-called auxiliary function [10] for each term of the objective function.
The remainder of this paper is organized as follows. We
will first briefly review the formulation of NMF with βdivergence in Section 2, and give a survey of the previous
methods in Section 3. We will then derive the proposed
multiplicative algorithm in Section 4, and finally present in
Section 5 basic experimental results validating our method
and comparing it to previous works.
2. NONNEGATIVE MATRIX FACTORIZATIONWITH β-DIVERGENCE
Given a matrix Y = (Yω,t )Ω×T ∈ R≥0,Ω×T and an integer
K, NMF is the problem of finding a factorization:
Y ≈ HU,(1)where H = (Hω,k )Ω×K ∈ R≥0,Ω×K and U = (Uk,t )K×T ∈R≥0,K×T are nonnegative matrices of dimensions Ω × K
and K × T , respectively. K is usually chosen such that
ΩK +KT ¿ ΩT , hence reducing the data dimension. Thisproblem can be formulated as the minimization of an objective function´
X ³X
D(Y | HU) =d Yω,t |Hω,k Uk,t ,
(2)
ω,t

k

where d is a scalar divergence.
A common way to measure how close Y and HU are is
to use a so-called β-divergence [14], defined by

yβ
xβ
yxβ−1


+
−
β ∈ R\{0, 1} β(β − 1)ββ−1.
dβ (y|x) = y(log y − log x) + (x − y) β = 1


y
y

 − log − 1β=0xx
It can be shown to be continuous in terms of β through thefollowing identities:
µ β−1
¶
y
− xβ−1xβ − y βlim dβ (y|x) = lim y
+
β→1
β→1
β−1
β
= y(log y − log x) + (x − y),µ β−1¶y β − xβyβ
x−+lim dβ (y|x) = lim y
β→0
β→0
1−β
β
β−1
y
y
= − log − 1.
xx
The choice of β should be driven by the type of data being analyzed and the application considered. In the NMFrelated literature, β = 1 is for example often used for soundsource separation [2], while β = 0.5 is used for the estimation of time-frequency activations [8] and β = 0 for
multipitch estimation [9]. How to choose β for multipitch
estimation and musical source separation is discussed in [5]and [3], respectively.
3. CONVENTIONAL ALGORITHMS
3.1. Multiplicative algorithms
The multiplicative gradient descent approach [10, 7] consists in updating each parameter by multiplying its value atthe previous iteration by a certain coefficient. Here, let θ
denote the set of all parameters {(Hω,k )Ω×K , (Uk,t )K×T }.
The derivative
P with respect
Pto Hω,k of the objective function
Jβ (θ) = ω,t dβ (Yω,t | k Hω,k Uk,t ) isÃ
!β−1
∂Jβ (θ) X X
Hω,k Uk,t
Uk,t
=
∂Hω,k
t
k
Ã
!β−2
X
X
Yω,t
−
Hω,k Uk,t
Uk,t .
t(5)

we obtain the following update rule for Hω,k :PHω,k ← Hω,kt Yω,t (

P

β−2

Uk,t
k Hω,k Uk,t )
.
P P
β−1
Uk,t
t(
k Hω,k Uk,t )(6)

A similar update rule can be obtained for Uk,t . Altogether,
the algorithm can be summarized as follows:H←H·(Y · (HU)β−2 )UT,(HU)β−1 UT

(7)

U←U·

HT (Y · (HU)β−2 ),HT (HU)β−1(8)where the symbol · and the fraction bar denote entrywise
matrix product and division respectively, and the exponentiations are also performed entrywise. Nonnegativity of theparameters is preserved through the updates, provided they
are initialized with nonnegative values.
The convergence of the conventional multiplicative updates have been proven only for 1 ≤ β ≤ 2, i.e, whendβ (y|x) is convex w.r.t. x [10, 11].
3.2. EM-based algorithms
In [13], NMF with EUC distance (β = 2), KL divergence(β = 1) and IS divergence (β = 0) is shown to be implicitin the following generative model of superimposed components,
X
Yω,t =
Cω,t,k .
(9)kThe components Cω,t,k act as latent variables and may be
used as complete data in the EM algorithm. We briefly
review here the update rules obtained through this method
successively for β = 2, β = 1 and β = 0.

NMF with EUC distance (β = 2) is equivalent to constrained ML estimation for the generative model (9) with(3)

Considering the following simple additive update for Hω,k ,∂Jβ (θ)
,
∂Hω,k

Hω,k
ηω,k = P P
,
β−1
Uk,t
t(
k Hω,k Uk,t )3.2.1. EUC-NMF

k

Hω,k ← Hω,k − ηω,kthe objective function will be decreasing if the coefficientsηω,k are all set equal to a sufficiently small positive number,as this corresponds to the conventional gradient descent. If
we now set

(4)

³ 1
K´
exp − (Cω,t,k − Hω,k Uk,t )2 2 .
K
2
σ
(10)
“Constrained” means here that the parameters H and U areestimated under the assumption that they are non-negative.
Cω,t,k ∼³ 2πσ 2 ´− 21Thereby, the following update rules are obtained based on
the EM algorithm:
"P
¡1
¢#
U
(Y
−
W
)
+
H
U
k,t
ω,t
ω,t
ω,k
k,t
t
K
P 2
, (11)
Hω,k =t Uk,t+"P
Uk,t =

ω Hω,k¡1− Wω,t ) + Hω,k Uk,tK (Yω,tP 2
t Hω,k¢#, (12)+

where [x]+ = max{x, 0} andWω,t =

X

Hω,k Uk,t .

(13)

4. NEW MULTIPLICATIVE ALGORITHMS
We consider the following optimization problem:P
P
Minimize Jβ (θ) = ω,t dβ (Yω,t | k Hω,k Uk,t )subject to∀ω, k, Hω,k ≥ 0, ∀k, t, Uk,t ≥ 0,

where θ denotes the set {(Hω,k )Ω×K , (Uk,t )K×T } of all parameters. The main result of this paper can be summarized
in the following
Theorem 1. The objective function Jβ (θ) is non-increasing
under the following updates:kµ

3.2.2. KL-NMF

H←H·

NMF with the generalized KL divergence (β = 1) is equivalent to ML estimation for the model (9) with

U←U·

Cω,t,k ∼ exp(−Hω,k Uk,t )(Hω,k Uk,t )Cω,t,k,Γ(Cω,t,k + 1)

(14)

where Γ denotes the Gamma function. Thereby, the following update rules are obtained based on the EM algorithm:
P
P
Uk,t (Yω,t / k Hω,k Uk,t )P,(15)
Hω,k = Hω,k tt Uk,tPPω,t /k Hω,k Uk,t )ω Hω,k (YP
Uk,t = Uk,t,(16)t Hω,kwhich coincide with the classical multiplicative updates.
3.2.3. IS-NMF
NMF with IS divergence (β = 0) is equivalent to constrained ML estimation for the model (9) with¡¢−1
Cω,t,k ∼ |πHω,k Uk,t | exp −|Cω,t,k |2 |Hω,k Uk,t |−1 .
(17)Thereby, the following update rules are obtained based on
the EM algorithm:
Hω,k =

Uk,t =
with1 X µ2ω,t,k + νω,t,k
,
T t
Uk,t
1 X µ2ω,t,k + νω,t,k
Ω

ω

Hω,k

,

Hω,k Uk,tµω,t,k = P
Yω,t ,k Hω,k Uk,tHω,k Uk,t X
νω,t,k = P
Hω,l Ul,t .
k Hω,k Uk,t
l6=k

(18)

µ

(Y · (HU)β−2 )UT(HU)β−1 UTHT (Y · (HU)β−2 )HT (HU)β−1

¶ϕ(β)
,

(22)

,

(23)

¶ϕ(β)

with


1/(2 − β) (β < 1)
ϕ(β) = 1(1 ≤ β ≤ 2) .


1/(β − 1) (β > 2)

The proof of the above theorem will be based on the auxiliary function approach, similarly to [10]. In the followings, we first explain the principle of the auxiliary function
method. Next, we prepare Lemma 1 and Lemma 2 for theconstruction of an auxiliary function. Lemma 3 gives anauxiliary function of the objective function Jβ (θ). Finally,
Theorem 1 is proven based on the principle of the auxiliary
function method and Lemma 3.
Let us first briefly review the concept of this approach.
Let G(θ) denote an objective function to be minimized w.r.t.
a parameter θ. A function G+ (θ, θ̂) which satisfies
G(θ) = min G+ (θ, θ̂)

(25)

θ̂

is then called an auxiliary function for G(θ), and θ̂ an auxiliary variable. The function G(θ) is non-increasing throughthe following iterative update rules:θ̂(s+1) ←argminθ̂ G+ (θ(s) , θ̂),

(19)

(24)

θ

(s+1)

+

← argminθ G (θ, θ̂

(s+1)

(26)
),

(27)

where θ̂(s+1) and θ(s+1) denote the updated values of θ̂ andθ after the s-th step. Then by construction, we have
(20)
(21)

G(θ(s) ) = G+ (θ(s) , θ̂(s) )≥ G+ (θ(s) , θ̂(s+1) )≥ G+ (θ(s+1) , θ̂(s+1) ) = G(θ(s+1) ).

(28)

This proves that G(θ) is non-increasing. By iteratively updating θ and θ̂, G(θ) will thus converge to a stationary point.
The auxiliary function will be suitably designed depending on the value of β. For β such that 1 ≤ β ≤ 2 [10, 11],such an auxiliary function can be constructed thanks to thefollowing lemma, refered to as Jensen’s inequality:
Lemma 1. Let f : R 7→ R be a convexPfunction. If λk (k =1, 2, · · ·, K) satisfies ∀k, λk > 0 and k λk = 1, then forxk (k = 1, 2, · · ·, K) ∈ R,
¡X ¢ X¡ xk ¢
.
(29)
f
xk ≤
λk f
λk
k
k
P
Equality holds when λk = xk / k xk .
Note that minimizing the auxiliary function obtained through the above lemma then leads to update equations which
are none other than the classical multiplicative update equations. However, Jensen’s inequality cannot be applied toβ < 1 and β > 2, because dβ (y|x) is then not convex w.r.t.
x. We are going to alleviate this problem by decomposing
the objective function into several terms which are going
to be either convex or concave depending on the value ofβ, and then use adequate inequalities to build an auxiliary
function. Indeed, if we write the objective function as
!β
Ã
X
1
1X X
Jβ (θ) =Hω,k Uk,tYω,t +β(β − 1) ω,t
β ω,t
k
Ã
!β−1
X
1 X
−
Hω,k Uk,tYω,t, (30)
β − 1 ω,tIf β now satisfies β ≤ 1, we apply Lemma 2 and obtainthe following inequality for the second term:´β1³XHω,k Uk,tβkβ−1≤ Zω,tk

Lemma 2. Let f : R 7→ R be a continuously differentiable
and concave function. Then, for any point z,
f (x) ≤ f 0 (z)(x − z) + f (z).

(31)

If β satisfies β ≥ 1, Lemma 1 leads to the following
inequality for the second term:
Ã
!β
¶β
µ
1 X
Hω,k Uk,t
1X
Hω,k Uk,tλω,t,k, (32)≤ββ
λω,t,kkkP
where ∀k, λω,t,k ≥ 0 and k λω,t,k = 1. Let θ̂ denotethe set of auxiliary variables {(λω,t,k )Ω×T ×K , (Zω,t )Ω×T },(β)
where Zω,t ∈ R will be used later on. We define Qω,t (θ, θ̂)as the right-hand side of Eq. (32). The equality holds whenHω,k Uk,tλω,t,k = P.
k Hω,k Uk,t

(33)

´ Zβ
ω,t
Hω,k Uk,t − Zω,t +
.
β

(34)

(β)

We define Rω,t (θ, θ̂) as the right-hand side of Eq. (34). The
equality holds whenX
Zω,t =
Hω,k Uk,t .
(35)
k
(β)

(β)

Note that Qω,t (θ, θ̂) = Rω,t (θ, θ̂) when β = 1.
The following inequalities for the third term can be derived similarly:
−

´β−1
1 ³X
Hω,k Uk,tβ−1k−Q(β−1)(θ, θ̂) (β ≤ 2)
ω,t≤
. (36)
−R(β−1) (θ, θ̂) (β ≥ 2)
ω,tThe equality holds when λω,t,k and Zω,t satisfy Eq. (33)and Eq. (35).
We can deduce the following lemma from the above.
Lemma 3. The functionkwe can see that, with respect to each parameter, the second
term is convex for β ≥ 1 and concave for β < 1, while
the third term is convex for β ≤ 2 and concave for β > 2.
To cope with concave terms, as in [15, 4], we shall use thefollowing lemma:³X

Jβ+ (θ, θ̂) =

X

X (β)
Yω,t
+
Sω,t (θ, θ̂),β(β − 1)
ω,t
ω,t

(37)

where
(β)

Sω,t (θ, θ̂)
 (β)
(β−1)

Rω,t (θ, θ̂) − Yω,t Qω,t (θ, θ̂)=(β < 1)(β)(β−1)
Qω,t (θ, θ̂) − Yω,t Qω,t (θ, θ̂) (1 ≤ β ≤ 2) , (38)



 (β)
(β−1)
Qω,t (θ, θ̂) − Yω,t Rω,t (θ, θ̂) (β > 2)

is an auxiliary function for Jβ (θ). Jβ+ (θ, θ̂) is minimized
w.r.t. θ̂ when θ̂ satisfies Eq. (33) and Eq. (35).
Proof of Lemma 3. Eq. (32) and Eq. (34) show that
Jβ (θ) ≤ Jβ+ (θ, θ̂).

(39)

The equality holds when θ̂ satisfies Eq. (33) and Eq. (35).
Thus, Eq. (33) and Eq. (35) minimizes Jβ+ (θ, θ̂) w.r.t. θ̂. ¤We are now ready to prove Theorem 1.
Proof of Theorem 1. Lemma 3 gives us an auxiliary function of Jβ (θ). According to the principle of the auxiliary∂Jβ+ (θ, θ̂)∂Hω,k= Vβ (θ) − Wβ (θ),

6

10

AUX
EM
MU
Objectivefunction method, we need to prove that minimizing Jβ+ (θ, θ̂)
w.r.t. θ and θ̂ iteratively lead to the update rules, Eq. (22)and Eq. (23).
First, we focus on minimizing Jβ+ (θ, θ̂) w.r.t. θ. The
derivative of Jβ+ (θ, θ̂) w.r.t. Hω,k is

5

10

(40)
0

10

1

10

2

10

3

10

Iteration #where
Vβ (θ) =(P

β−1t Zω,t Uk,t
β−1 P
1−β
β
Hω,k
t λω,t,k Uk,t(β < 1),(β ≥ 1)

(
β−2 P
β−1
Hω,k
λ2−β Yω,t Uk,tWβ (θ) = P β−2t ω,t,kt Zω,t Yω,t Uk,t(41)

(β ≤ 2). (42)(β > 2)The second derivative is∂ 2 Jβ+ (θ, θ̂)∂Hω,k ∂Hω0 ,k0= {Vβ0 (θ) − Wβ0 (θ)}δω,ω0 δk,k0 ,(43)where δi,j is 1 if i = j, otherwise 0 and(
0(β < 1)Vβ0 (θ) =
,
β−2 P
1−β
β
(β − 1)Hω,k(β ≥ 1)t λω,t,k Uk,t
(
2−β
β−1
β−3 P
(β ≤ 2)(β − 2)Hω,k0t λω,t,k Yω,t Uk,t.
Wβ (θ) =0(β > 2)Thus, Jβ+ (θ, θ̂) is a convex function in H. Setting the first
derivative to zero, we obtain the update rule for Hω,k :
Ã
1
! 2−β
P 2−β
β−1



t λω,t,k Yω,t Uk,t(β < 1)P β−1Zω,t Uk,t

t

P

2−β
β−1

t λω,t,k Yω,t Uk,tHω,k =(1 ≤ β ≤ 2) .
 P λ1−β U β


t
ω,t,k
k,t


Ã P β−2! 1Zω,t Yω,t Uk,t β−1t(β > 2)P 1−β β


λ
U
t

ω,t,k

k,t

(44)
U can be discussed similarly.
Next, we consider the auxiliary variables θ̂. Eq. (33)and Eq. (35) minimize Jβ+ (θ, θ̂) w.r.t. θ̂. Thus, minimizing Eq. (33) and Eq. (35) into Eq. (44) gives the following
update rule:
!ϕ(β)
ÃP
P
β−2
Y
(
H
U
)
U
ω,t
ω,k
k,t
k,t
t
k
.
Hω,k ← Hω,k
P P
β−1
(
H
Uk,t
ω,k Uk,t )tk
The update rule for Uk,t can be obtained similarly. The update rules for H and U can be simply rewritten as Eq. (22)and Eq. (23).
¤Fig. 1. Evolution in log-log scale of the objective function
with β = 0.
5. EXPERIMENTS
The convergence speed is compared with existing algorithms.
The classical multiplicative algorithm for NMF will be denoted as “MU”, the EM-based algorithm as “EM”, and the
proposed multiplicative algorithm based on the auxiliary
function method as “AUX”. NMF is often applied to analysis of audio signals. Here, we use as input data matrix
the magnitude spectrogram of 8 second length music signal (generated from RWC-MDB-P-2001 No.25 [16]) downmixed to monaural and downsampled to 16kHz. It was computed using the short time Fourier transform with a 32 ms
long Hanning window and with 16 ms overlap.
We compared the performances of all the algorithms forthree different values of β, namely β = 0, 0.5, 2. Fig. 1
shows the results for β = 0. In this case, EM is the slowest and MU the fastest, while AUX is slightly slower than
MU. As shown in Fig. 2, for β = 2, AUX (which is thenequivalent to MU) is again faster than EM. Finally, the results for β = 0.5 are shown in Fig. 3. In all cases, MU isslightly faster than AUX, however, the convergence of ouralgorithm is theoretically proven.
6. CONCLUSIONSIn this paper, we proposed a convergence-guaranteed multiplicative algorithm for NMF with β-divergence. The formof the updates is similar to that of the conventional multiplicative algorithm but with a different exponent term. We
confirmed through basic experiments that the proposed algorithms converge faster than EM algorithms. Future work
will include the extension of our auxiliary function approachto the stable algorithms for constrained NMF methods whichhave objective function and penalty terms, such as sparsity
or smoothness.
ACKNOWLEDGMENTS
We thank the reviewers for very helpful suggestions andcomments.

10

10

AUX
EM

9

Objective

10[6] S. Eguchi and Y. Kano, “Robustifying maximum likelihood estimation,” Tokyo Institute of Statistical Mathematics, Tech. Rep., 2001.

8

10

7

10

6

10

5

10

0

10

1

10

2

10

3

10

[7] A. Cichocki, R. Zdunek, and S. Amari, “Csiszárs divergences for non-negative matrix factorization : Family of new algorithms,” in Proc. Int. Conf. Independent Component Analysis and Blind Signal Separation, Mar. 2006, pp. 32–39.

Iteration #Fig. 2. Evolution in log-log scale of the objective function
with β = 2. MU is equivalent to AUX.
7

10

AUX
MU

[9] N. Bertin, R. Badeau, and E. Vincent, “Enforcing harmonicity and smoothness in Bayesian non-negative
matrix factorization applied to polyphonic music transcription,” IEEE Trans. on Audio, Speech, and Language Processing, pp. 538–549, 2010.

6

Objective[8] R. Hennequin, R. Badeau, and B. David, “NMF with
time-frequency activations to model non stationary
audio events,” in Proc. International Conference on
Acoustics, Speech and Signal Processing, Mar. 2010,
pp. 445–448.

10

5

10

0

10

1

10

2

10

3

10

Iteration #Fig. 3. Evolution in log-log scale of the objective function
with β = 0.5.
7. REFERENCES[1] D. D. Lee and H. S. Seung, “Learning the parts ofobjects by non-negative matrix factorization,” Nature,
vol. 401, pp. 788–791, Oct. 1999.
[2] T. Virtanen, “Monaural sound source separation by
nonnegative matrix factorization with temporal continuity and sparseness criteria,” IEEE Transactions on
Audio, Speech, and Language Processing, vol. 15, pp.
1066–1074, Mar. 2007.
[3] D. FitzGerald, M. Cranitch, and E. Coyle, “On theuse of the beta divergence for musical source separation,” in Proc. of Irish Signals and Systems Conference, 2009.

[10] D. D. Lee and H. S. Seung, “Algorithms for nonnegative matrix factorization,” in Proc. of the Conference on Advances in Neural Information Processing
Systems, vol. 13, Dec. 2001, pp. 556–562.
[11] R. Kompass, “A generalized divergence measure fornonnegative matrix factorization,” Neural Computation, vol. 19, no. 3, pp. 780–791, Mar. 2007.
[12] R. Badeau, N. Bertin, and E. Vincent, “On the stability of multiplicative update algorithms. application to
non-negative matrix factorization,” in Telecom ParisTech, Technical report, 2009.
[13] C. Févotte and A. T. Cemgil, “Nonnegative matrix
factorizations as probabilistic inference in composite
models,” in Proc. European Signal Processing Conference, vol. 47, 2009, pp. 1913–1917.
[14] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative
matrix factorization with the Itakura-Saito divergence.
with application to music analysis,” Neural Computation, vol. 21, pp. 793–830, Mar. 2009.

[4] H. Kameoka, N. Ono, K. Kashino, and S. Sagayama,“Complex NMF: A new sparse representation for
acoustic signals,” in Proc. of International Conferenceon Acoustics, Speech and Signal Processing, 2009, pp.
3437–3440.

[15] H. Kameoka, N. Ono, and S. Sagayama, “Auxiliary
function approach to parameter estimation of constrained sinusoidal model,” in Proc. of International
Conference on Acoustics, Speech and Signal Processing, Apr. 2008, pp. 29–32.

[5] S. A. Raczynski, N. Ono, and S. Sagayama, “Extending nonnegative matrix factorization – a discussion in
the context of multipitch frequency estimation of musical signals,” in Proc. of European Signal Processing
Conference, Aug. 2009, pp. 934–938.

[16] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka,“RWC music database: Popular, classical, and jazzmusic database,” in Proc. International Conference onMusic Information Retrieval, 2002, pp. 287–288.