COMPLEX NMF WITH THE GENERALIZED KULLBACK-LEIBLER DIVERGENCE
Hirokazu Kameoka† , Hideaki Kagami‡ , Masahiro Yukawa‡† NTT Communication Science Laboratories, NTT Corporation, Japan.
‡ Dept. Electronics and Electrical Engineering, Keio University, Japan.
spectra scaled by time-varying amplitudes, and amounts toapproximating the observed spectrogram as the sum of rank-1
spectrograms. In a supervised/semi-supervised setting, NMF
is first employed to train the basis spectra of each soundsource using individually recorded audio samples. At test
time, NMF is applied to the spectrogram of a test mixture
signal, where the subsets of the basis spectra are fixed at the
pretrained spectra.
Although the NMF approach has been shown to be successful, one drawback is that it assumes the additivity ofmagnitude (or power) spectra, which holds only approximately, and does not take account of phase information. To
address this drawback, we have previously proposed a framework called the “Complex NMF (CNMF)” [5], where thecomplex spectrum observed at each time frame is modeled
as the sum of components, each of which is described by
the multiplication of a static basis spectrum, a time-varying
amplitude and a time-varying phase spectrum. With a similar
motivation, Parry and Essa [6] and Fevotte et al. [7] proposeda generative model of the complex spectrogram obtained withthe short-time Fourier transform (STFT) of a mixture signal,
where the power spectrogram of each component is modeled
as a rank-1 matrix whereas the phase spectrogram is treatedas uniformly distributed latent variables. It can be shownthat when we assume that each element of the complex spectrogram independently follows a zero-mean complex normal
distribution, the maximum likelihood estimation of the model
parameters amounts to fitting the NMF model to an observed
power spectrogram using the Itakura-Saito (IS) divergence asa goodness-of-fit criterion. This approach is called IS-NMF.
A similar kind of generative model using a complex Cauchy
distribution instead of a complex normal distribution has
also been proposed [8]. IS-NMF and Cauchy-NMF treat thephase spectrogram of each underlying component as a latentvariable to be marginalized out and the aim is to find the
expectation of the complex spectrogram of each component
taken over all possible phase spectrograms of all the components. Although this estimator is reasonable if the phase
spectrograms are really stochastic, they are in fact deterministic and unique. Indeed, the phase part of this estimator isalways given as that of the observed mixture signal, which
certainly differs from the true value. By contrast, CNMF
allows us to find the jointly optimal estimates of the power
and phase spectrograms of all the components. Thus, we can
expect CNMF to lead to higher source separation accuracy
than the other NMF variants if the parameters can be properly
estimated.
However, one limitation with the CNMF framework isthat the divergence measure used to measure the difference
between an observed spectrogram and the model is limited
to only the Euclidean distance (squared error) due to the factthat the arguments are complex numbers unlike NMF. This is
in contrast to the conventional NMF framework where effi-

ABSTRACT
We previously introduced a phase-aware variant of the
non-negative matrix factorization (NMF) approach for audio source separation, which we call the “Complex NMF(CNMF).” This approach makes it possible to realize NMFlike signal decompositions in the complex time-frequency
domain. One limitation of the CNMF framework is that thedivergence measure is limited to only the Euclidean distance.
Some previous studies have revealed that for source separation tasks with NMF, the generalized Kullback-Leibler (KL)divergence tends to yield higher accuracy than when using
other divergence measures. This motivated us to believe that
CNMF could achieve even greater source separation accuracy
if we could derive an algorithm for a KL divergence counterpart of CNMF. In this paper, we start by defining the notionof the “dual” form of the CNMF formulation, derived fromthe original Euclidean CNMF, and show that a KL divergence
counterpart of CNMF can be developed based on this dual
formulation. We call this “KL-CNMF”. We further derivea convergence-guaranteed iterative algorithm for KL-CNMF
based on a majorization-minimization scheme. The source
separation experiments revealed that the proposed KL-CNMF
yielded higher accuracy than the Euclidean CNMF and NMF
with varying divergences.
Index Terms— Audio source separation, non-negative
matrix factorization (NMF), Complex NMF, generalized
Kullback-Leibler (KL) divergence
1. INTRODUCTION
Audio source separation has long been a challenging task in
the field of audio signal processing. A deep neural networkbased approach has recently proved powerful for supervised
audio source separation tasks where the mixture consists ofspeech and noise [1]. Furthermore, a recently proposed approach called “deep clustering” [2] has made it possible todeal with “cocktail party” scenarios where the interferenceis also speech. Although these methods have been shown
to work well when a large number of training samples areavailable, the non-negative matrix factorization (NMF) approach [3, 4] still remains attractive for audio source separation tasks particularly where only a limited amount of training
data is available or when prior knowledge about the underlying sources is limited. In addition, since NMF is a generative
approach, it can be convenient in semi-supervised scenarios.
With the NMF approach, the magnitude (or power) spectrogram of a mixture signal, interpreted as a non-negative
matrix Y , is modeled as the product of two non-negative matrices H and U . This can be interpreted as approximating theobserved spectrum at each time frame as a linear sum of basis
This work was supported by JSPS KAKENHI Grant Numbers
26730100, 15K06081, 15K13986, 15H02757.

978-1-5090-4117-6/17/$31.00 ©2017 IEEE

56

ICASSP 2017

We use F(θ) to denote an objective function that we want
to minimize with respect to θ. F + (θ, α) is defined as an auxiliary function for F(θ) if it satisfies

cient algorithms have been proposed for different divergencemeasures [7–14]. Some previous studies have revealed thatfor source separation tasks with NMF, using the generalized
Kullback-Leibler (KL) divergence tends to yield higher accuracy than using other divergence measures [14] (though consensus has yet to be reached [15]). This motivated us to expectthat we can achieve even higher source separation accuracy if
we can derive an algorithm for a KL divergence counterpartof CNMF. In this paper, we start by defining the notion of
the “dual” form of the CNMF formulation, derived from theauxiliary function of the original Euclidean CNMF, and showthat a KL divergence counterpart of CNMF can be developed
based on this dual formulation. We call this “KL-CNMF”. We
further derive a convergence-guaranteed algorithm for KLCNMF based on a majorization-minimization scheme.

F(θ) = min F + (θ, α).

We call α an auxiliary variable. By using F + (θ, α), F(θ) canbe iteratively decreased according to the following theorem:
Theorem 1. F(θ) is non-increasing under the updates, θ ←argminθ F + (θ, α) and α ← argminα F + (θ, α).
It can be shown [5] that+IEU(θ, α) =2. EUCLIDEAN COMPLEX NMF

L∑

Hk,l Ul,m ,∑ |Xl,k,m − Hk,l Ul,m cl,k,m |2

+λ

∑{

}
2
p|Vl,m |p−2 Ul,m
+ (2 − p)|Vl,m |p , (6)l,mis an auxiliary function of IEU (θ) where ∑α = {X, V }, βl,k,m
can be any positive number satisfying l βl,k,m = 1, andXl,k,m ∈ C and Vl,m ≥ 0 are auxiliary variables satisfying∑
Xl,k,m = Yk,m .
(7)

(1)

l

l=1

+
By using IEU(θ, α), the update rules for α are derived as

CNMF approximates the complex spectrum Yk,m at eachframe m as a linear sum of static basis spectra Hk,1 , . . . , Hk,Lscaled by time-varying amplitudes U1,m , . . . , UL,m and multiplied by time-varying phase spectra ejϕ1,k,m , . . . , ejϕL,k,m∑
Yk,m ≃
Hk,l Ul,m ejϕl,k,m .
(2)

Xl,k,m =Hk,l Ul,m cl,k,m()∑+ βl,k,m Yk,m −Hk,l Ul,m cl,k,m ,(8)lVl,m =Ul,m ,lHere, it is important to emphasize that ϕl,k,m is indexed bym, meaning that this model allows the phase spectrum ofeach component to vary freely over time. For simplicity ofnotation, let us hereafter put cl,k,m = ejϕl,k,m . In [5], we
considered an objective function∑∑
2
IEU (θ) =Yk,m −
Hk,l Ul,m cl,k,m + R(U ), (3)k,mβl,k,ml,k,m

Let Yk,m ∈ C denote the complex spectrogram of an observed
mixture signal, where k and m denote the frequency and time(frame) indices, respectively. While the NMF approach approximates the magnitude (or power) spectrum |Yk,m | at eachframe m as a linear sum of static basis spectra Hk,1 , . . . , Hk,Lscaled by time-varying amplitudes U1,m , . . . , UL,m
|Yk,m | ≃(5)

αand the update rules for θ are derived as∑
m Ul,m |Xl,k,m |/βl,k,m,Hk,l = √∑ ∑2(km Ul,m |Xl,k,m |/βl,k,m )∑Hk,l |Xl,k,m |/βl,k,m
Ul,m = ∑ k 2
p−2 ,k Hk,l /βl,k,m + λpVl,mcl,k,m =Xl,k,m /|Xl,k,m |.

lwhere θ = {H, U , C} and R(U ) is a regularization termfor U . It is important to note that the complex NMF modelallows the components to cancel each other out, and so someconstraint is needed to induce the sparsity of U . For this purpose, we define R(U ) using the ℓp norm∑
R(U ) = 2λ|Ul,m |p ,
(4)

(9)

(10)
(11)
(12)

The CNMF algorithm can thus be summarized as follows:
1. Initialize H, U and C.
2. Update X and V using (8) and (9).
3. Update H, U and C using (10)–(12) and return to 2.
We call this “Euclidean CNMF (EU-CNMF)”.
Here, the auxiliary variable Xl,k,m can be viewed as an
estimate of the complex spectrogram of the l-th signal component. At step 2, Xl,k,m is updated by adding the portionof the error between the observed spectrogram and the modelto the current estimate of Hk,l Ul,m cl,k,m . cl,k,m is then updated at its argument Xl,k,m /|Xl,k,m |, and H and U are updated using its magnitude |Xl,k,m |. Although the details areomitted owing to space limitations, it can be shown that (8) isan MMSE estimator of Xl,k,m , i.e., E[Xl,k,m |Y , H, U , C],
when seen from a generative model perspective. This is incontrast to the Wiener filter E[Xl,k,m |Y , H, U ], which onlyuses the estimate of the power spectrogram.

l,mwhere λ > 0 weighs the importance of the sparsity cost relative to the fitting cost. When 0 < p < 2, R(U ) promotes
sparsity
if the norm∑
∑ of U is bounded. To bound U , we assume2
k Hk,l = 1 ork Hk,l = 1.
Although it is difficult to solve the above optimization
problem analytically, a convergence-guaranteed algorithm forfinding a stationary point can be developed based on the auxiliary function concept. To derive the algorithm, we first introduce the general principle of the auxiliary function approach(majorization-minimization approach) [9, 16, 17].

57

3. KULLBACK-LEIBLER COMPLEX NMF
Interestingly, it can be shown that the algorithm presented
above also converges to a stationary point of the following
optimization problem:
minimize JEU (θ̄)∑
Xl,k,m = Yk,m ,subject to

(13)

l

Fig. 1. KL-CNMF derivation process.

where θ̄ = {H, U , X} and∑ (|Xl,k,m | − Hk,l Ul,m )2

JEU (θ̄) =l,k,m

βl,k,m+ R(U ).

subject to(14)=|Xl,k,m |

2

(18)

which we call “KL-CNMF.”
The objective JKL (θ̄) is non-differentiable with respect
to Xl,k,m . In the following, we construct an easy-to-optimize
auxiliary function to obtain a closed form update rule for
Xl,k,m . First, we can show that

(15)

2
2
− 2Hk,l Ul,m Re[c∗l,k,m Xl,k,m ] + Hk,lUl,m.

|Xl,k,m | log |Xl,k,m | − |Xl,k,m | log Hk,l Ul,m − |Xl,k,m |{}|Xl,k,m |−Zl,k,m≤|Xl,k,m |
+
log
Z
l,k,m
Zl,k,m∗Now, by using the fact that a function g(z) = −Re(c z) withcomplex arguments z and |c| = 1 is a tangent plane to a conef (z) = −|z|, where c indicates the direction of the tangent
line, we obtain an inequality−Re(c∗ z) ≥ −|z|.

Xl,k,m = Yk,m ,lThis can be confirmed as follows. The first term of the auxiliary function I + (θ, α) of EU-CNMF can be written as|Xl,k,m − Hk,l Ul,m cl,k,m |2∑− |Xl,k,m | log Hk,l Ul,m − |Xl,k,m |()
Zl,k,m
|X
|2
+
|X
|
log
−
2
,
= Zl,k,ml,k,mHk,l Ul,m
l,k,m

(16)

(19)

by using the fact that a logarithmic function is a concave function and that for any concave function f (x), f (x) is below orequal to its tangent at any point, namely f (x) ≤ f ′ (z)(x −z)+f (z). The equality of (19) holds when Zl,k,m = |Xl,k,m |.
The right-hand side of (19) still involves a non-differentiable
term |Xl,k,m |. Here, the coefficient

We therefore obtain
2
2|Xl,k,m |2 − 2Hk,l Ul,m Re[c∗l,k,m Xl,k,m ] + Hk,lUl,m2
2≥|Xl,k,m |2 − 2Hk,l Ul,m |Xl,k,m | + Hk,lUl,m=(|Xl,k,m | − Hk,l Ul,m )2 .

Z

dl,k,m = log Hk,ll,k,m
Ul,m − 2,Hence, JEU (θ̄) ≤ I + (θ, α). This implies that (6) is also an
auxiliary function of JEU (θ̄). Here, the difference is that theroles of X = {Xl,k,m } and C = {cl,k,m } are reversed: X
is a model parameter and C is an auxiliary variable, while
X is an auxiliary variable and C is a model parameter inthe original CNMF. Thus, there is a duality between the optimization problem of the original CNMF and (13). (13)can be interpreted as a problem of decomposing an observed
complex spectrogram Yk,m into the sum of L components
X1,k,m , . . . , XL,k,m such that the magnitude spectrogram ofeach component is as close as possible to a rank-1 structure.
This gives a different explanation to the objective of CNMF.
We call this formulation technique “dual formulation.”
Since the KL divergence only allows non-negative arguments, it cannot be straightforwardly
used to measure the∑difference between Yk,m and l Hk,l Ul,m cl,k,m . However,since the dual CNMF formulation uses the Euclidean distance
between non-negative values, |Xl,k,m | and Hk,l Ul,m , as thecost function, we can also use the KL divergence to measure
their difference, which leads us to an objective function∑
JKL (θ̄) =DKL (|Xl,k,m |∥Hk,l Ul,m ) + R(U ), (17)(20)

can be either non-negative or negative. According to the signof dl,k,m , we can use the following inequalities|X|2

Wl,k,m,2(21)−|Xl,k,m | ≤ −Re[c∗l,k,m Xl,k,m ],(22)|Xl,k,m | ≤ 2Wl,k,m+
l,k,m
to obtain|Xl,k,m |2
Zl,k,m + dl,k,m |Xl,k,m |∗≤Al,k,m |Xl,k,m |2 − 2Re[Bl,k,m
Xl,k,m ] + Dl,k,m ,

where

{ d

Al,k,m :={Bl,k,m :=l,k,m1
2Wl,k,m + Zl,k,m
1
Zl,k,m0−dl,k,m cl,k,m /2(23)(dl,k,m ≥ 0),(dl,k,m < 0)(24)

(dl,k,m ≥ 0).
(dl,k,m < 0)(25)

Dl,k,m is given by dl,k,m Wl,k,m /2 when dl,k,m ≥ 0 and
0 otherwise. The equalities of (21) and (22) hold whenWl,k,m = |Xl,k,m | and cl,k,m = Xl,k,m /|Xl,k,m |, respectively. Thus, we obtain an auxiliary function of JKL (θ̄) as

l,k,m

where DKL (x∥y) = x log xy − x + y. Thus, we can consider
the following optimization problem+∗
JKL(θ̄, ᾱ) =Al,k,m |Xl,k,m |2 − 2Re[Bl,k,m
Xl,k,m ]minimize JKL (θ̄)

58

+ Dl,k,m + Hk,l Ul,m + R+ (U , V ),(26)

where θ̄ = {H, U , X} is a set of parameters, ᾱ = {Z, W , C, V }is a set of auxiliary variables, and R+ (U , V ) is an auxiliary
function of R(U ), given by
R+ (U , V )
{ ∑
p−1
p
2λ l,m {pVl,m(Ul,m − Vl,m ) + Vl,m
}
∑
=
p−2 2
p
λ l,m {pVl,m Ul,m + (2 − p)Vl,m }(27)(0 < p ≤ 1).
(1 < p ≤ 2)+
As can be seen from (26), JKL(θ̄, ᾱ) is a quadratic function of
Xl,k,m . Thus, by using the method of Lagrange multipliers,
we obtain the update rule for Xl,k,m analytically as

Xl,k,m =(1
Al,k,mYk,m −∑
Bl,k,m +∑ Bl,k,m )l Al,k,m
1

.

(28)

l Al,k,mThis provides yet another form of filter that optimally decomposes Yk,m into L complex spectrograms X1,k,m , . . . , XL,k,musing the estimates of H, U and c. When 0 < p ≤ 1, theupdate rules for H and U can be derived as∑
m |Xl,k,m |,(29)
Hk,l = ∑ ∑km |Xl,k,m |∑|Xl,k,m |Ul,m = k(30)
p−1 .
1 + λpVl,m∑Note that here we have used the ℓ1 norm constraint k Hk,l =1 for H. The update rule for U when 1 < p ≤ 2 can also
be derived in closed form. The update rules for the auxiliary
variables ᾱ are derived as
Zl,k,m = |Xl,k,m |,Wl,k,m = |Xl,k,m |,cl,k,m = Xl,k,m /|Xl,k,m |,Vl,m = Ul,m .

Fig. 2. Average SNR improvements obtained with KL-CNMF (proposed), EU-CNMF, KL-NMF, EU-NMF and IS-NMF. We used threemusic recordings from the SiSEC 2013 database [18], “Ultimate NZ
Tour”(top), “Bearlin - Roads”(middle) and “Fort Minor - Remember
the Name”(bottom), for the test data.

We performed 3-fold cross validation. We partitioned eachrecording into three segments, used one segment as the test
data and the other two segments as the training data, repeatedsignal-to-noise (SNR) evaluations three times with different
test segments, and took the average of the SNR improvements
obtained with the three repeated rounds. With all these methods, the basis spectra were pretrained using the individual
tracks of the training data, and then source separation was
performed on the test data. 6 basis spectra were assigned
to each track. Thus, a total of 30 basis spectra were used
for the separation. All the audio samples were monauraland sampled at 22.05kHz. An STFT was computed using
a square-root Hanning window that was 32ms long with a
16ms overlap. Fig. 2 shows the SNR improvements after
the separations with all the methods. From these results, we
confirmed that KL-CNMF outperformed the other methods.

(31)
(32)
(33)
(34)

Overall, the optimization algorithm for KL-CNMF can besummarized as follows:
1. Initialize H, U and X.
2. Update Z, W , C and V using (31)–(34).
3. Update H, U and X using (28)–(30) and return to 2.
For the initialization (step 1), we can use conventional NMF
algorithms followed by Wiener filtering or the EU-CNMF algorithm to obtain the estimates of H, U and X.
The derivation process of KL-CNMF and its relationshipto EU-CNMF are illustrated in Fig. 1.

5. CONCLUSIONSCNMF is a phase-aware variant of the NMF approach foraudio source separation, which makes it possible to realize
NMF-like signal separations in the complex time-frequency
domain. One limitation of the conventional CNMF is that thedivergence measure is limited to the Euclidean distance. This
paper proposed a KL divergence counterpart of CNMF, which
we call “KL-CNMF,” and derived an algorithm for finding a
locally optimal solution. We confirmed through supervised
source separation experiments that KL-CNMF outperformed
other NMF variants.

4. EXPERIMENTS
We conducted supervised source separation experiments tocompare the source separation accuracy of KL-CNMF (proposed), EU-CNMF, KL-NMF, EU-NMF and IS-NMF. We
used three music recordings from the SiSEC 2013 database,available at [18], as the experimental data. Each recording isa mixture of 5 tracks, each of which is produced by a single
instrument or singer. The separated tracks are also available.

59[13] C. Févotte and J. Idier, “Algorithms for nonnegative matrix factorization with the β-divergence,” Neural Computation, vol. 23, no. 9, 2011.

6. REFERENCES[1] Y. Wang and D. Wang,“Towards scaling up
classification-based speech separation,” IEEE Trans.
Audio, Speech, Language Process., vol. 21, no. 7, pp.
1381–1390, 2013.

[14] D FitzGerald, M Cranitch, and E Coyle, “On the useof the beta divergence for musical source separation,”in Proceedings of the Irish Signals and Systems Conference (ISSC), 2009.

[2] J. R. Hershey, Z. Chen, J. Le Roux, and S. Watanabe,“Deep clustering: Discriminative embeddings for segmentation and separation,” in Proceedings of the International Conference on Acoustics, Speech and Signal
Processing (ICASSP), 2016, pp. 31–35.

[15] B. King, C. Févotte, and P. Smaragdis, “Optimal costand magnitude power for nmf-based speech separationand music interpolation,” in Proceedings of the IEEE
International Workshop on Machine Learning for Signal
Processing (MLSP), 2012, pp. 23–26.

[3] P. Smaragdis, B. Raj, and M Shashanka, “Supervisedand semi-supervised separation of sounds from singlechannel mixtures,” in Proceedings of the International
Conference on Independent Component Analysis andSignal Separation (ICA), 2007, pp. 414–421.

[16] J. M. Ortega and W. C. Rheinboldt, Iterative solutionsof nonlinear equations in several variables, Academic
Press, New York, 1970.
[17] D. R. Hunter and K. Lange, “Quantile regression via an
MM algorithm,” Journal of Computational and Graphical Statistics, vol. 9, pp. 60–77, 2000.

[4] F. Weninger, J. Le Roux, J. R. Hershey, and S. Watanabe, “Discriminative NMF and its application to singlechannel source separation,” in Proceedings of the Annual Conference of the International Speech Communication Association (Interspeech), 2014, pp. 536–543.

[18] “The main page for the fourth community-based signal separation evaluation campaign (SiSEC 2013),”https://sisec.wiki.irisa.fr/.

[5] H. Kameoka, N. Ono, K. Kashino, and S. Sagayama,“Complex NMF: A new sparse representation for acoustic signals,” in Proceedings of the International Conference on Acoustics, Speech and Signal Processing(ICASSP), 2009, pp. 3437–3440.
[6] R. M. Parry and I. Essa, “Phase-aware non-negative
spectrogram factorization,” in Proceedings of the International Conference on Independent Component Analysis and Signal Separation (ICA), 2007, pp. 536–543.
[7] 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, no. 3, pp. 793–830, 2009.
[8] A. Liutkus, D. Fitzgerald, and R. Badeau, “Cauchy nonnegative matrix factorization,” in Proceedings of theIEEE Workshop on Applications of Signal Processing to
Audio and Acoustics (WASPAA), 2015.
[9] D. D. Lee and H. S. Seung, “Algorithms for nonnegative matrix factorization,” in Advances in Neural
Information Processing Systems (NIPS). 2000, pp. 556–
562, MIT Press.
[10] S. Dhillon and S. Sra, “Generalized nonnegative matrix approximations with Bregman divergences,” inAdvances in Neural Information Processing Systems(NIPS), 2005, pp. 283–290.
[11] A. Cichocki, R. Zdunek, and S. Amari, “Csiszar’s divergences for non-negative matrix factorization: Family of new algorithms,” in Proceedings of the International Conference on Independent Component Analysis
and Signal Separation (ICA), 2006, pp. 32–39.
[12] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano,
N. Ono, and S. Sagayama, “Convergence-guaranteed
multiplicative algorithms for non-negative matrix factorization with beta-divergence,” in Proceedings of theIEEE International Workshop on Machine Learning for
Signal Processing (MLSP), 2010, pp. 283–288.

60