IEICE TRANS. FUNDAMENTALS, VOL.E102–A, NO.2 FEBRUARY 2019458

LETTER

Independent Low-Rank Matrix Analysis Based on Generalized
Kullback–Leibler Divergence∗
Shinichi MOGAMI†a) , Yoshiki MITSUI† , Norihiro TAKAMUNE† , Daichi KITAMURA†† , Nonmembers,
Hiroshi SARUWATARI†b) , Member, Yu TAKAHASHI††† , Nonmember, Kazunobu KONDO††† , Member,
Hiroaki NAKAJIMA††† , Nonmember, and Hirokazu KAMEOKA†††† , Member

SUMMARY In this letter, we propose a new blind source separation method, independent low-rank matrix analysis based on generalized
Kullback–Leibler divergence. This method assumes a time-frequencyvarying complex Poisson distribution as the source generative model, whichyields convex optimization in the spectrogram estimation. The experimental
evaluation confirms the proposed method’s efficacy.
key words: blind source separation, nonnegative matrix factorization,
Poisson distribution, Kullback–Leibler divergence

1.

Introduction

Blind source separation (BSS) [1]–[3] is a technique for extracting specific sources from an observed multichannel mixture signal without knowing a priori information about the
mixing system. Let N and M be the numbers of sources andchannels, respectively. The short-time Fourier transforms(STFTs) of the multichannel source, observed, and estimatedsignals are defined as si j = (si j1, . . . , si j N ) > ∈ C N , xi j =(x i j1, . . . , x i j M ) > ∈ C M , and yi j = ( yi j1, . . . , yi j N ) > ∈ C N ,where i = 1, . . . , I; j = 1, . . . , J; n = 1, . . . , N; andm = 1 . . . , M are the integral indices of the frequency
bins, time frames, sources, and channels, respectively, and> denotes the transpose. We assume the mixing system
xi j = Ai si j , where Ai = (ai1, . . . , ai N ) ∈ C M×N is afrequency-wise mixing matrix and ain is the steering vector
for the nth source. When M = N and Ai is not a singular matrix, the estimated signal yi j can be expressed as yi j = Wi xi j ,Hwhere Wi = A−1
i = (wi1, . . . , wi N ) is the demixing matrix,
win is the demixing filter for the nth source, and H denotes
the Hermitian transpose. Multichannel nonnegative matrix
factorization (MNMF) [4], [5] is a BSS method that simultaneously estimates both the low-rank time-frequency structure
and the spatial covariance matrix for each source, indirectly
Manuscript received August 2, 2018.
Manuscript revised October 10, 2018.
† The authors are with the University of Tokyo, Tokyo, 113-8656
Japan.
†† The author is with National Institute of Technology, Kagawa
College, Takamatsu-shi, 761-8058 Japan.
††† The authors are with Yamaha Corporation, Hamamatsu-shi,
430-8650 Japan.
†††† The author is with Nippon Telegraph and Telephone Corporation, Atsugi-shi, 243-0198 Japan.
∗ This work was supported by SECOM Science and Technology
Foundation and JSPS KAKENHI Grant Number JP17H06572.
a) E-mail: shinichi_mogami@ipc.i.u-tokyo.ac.jpb) E-mail: hiroshi_saruwatari@ipc.i.u-tokyo.ac.jp
DOI: 10.1587/transfun.E102.A.458

identifying the mixing system. Recently, independent lowrank matrix analysis (ILRMA) [6], which is a unification of
direct estimation of the demixing matrix [7] and simple nonnegative matrix factorization (NMF) [8], [9], was proposed
as a state-of-the-art BSS method. In terms of optimization,
ILRMA is faster and more stable than MNMF.
Conventional MNMF and ILRMA assume the timefrequency-varying complex Gaussian distribution in the generative model of each source spectrogram, which corresponds to NMF based on Itakura–Saito (IS) divergence (ISNMF). We refer to the conventional ILRMA based on IS
divergence as IS-ILRMA. Recently, the source generativemodels assumed in ILRMA and MNMF were generalized tothe time-frequency-varying complex Student’s t-distribution(t-MNMF [10] and t-ILRMA [11]). These conventional
ILRMAs always include a non-convex optimization in the
spectrogram modeling that is sensitive to the initializationand diversifies the separation quality. In this letter, we propose ILRMA based on generalized Kullback–Leibler (KL)divergence (KL-ILRMA) by using a time-frequency-varying
complex Poisson distribution as the source generative model.
In KL-ILRMA, the estimation of a low-rank spectrogram
model results in KL-divergence-based NMF (KL-NMF) [8],which is a convex optimization in terms of each decomposed
matrix variable. To our knowledge, the proposed method is
the world’s first attempt to realize convex-optimization-based
ILRMA w.r.t. source modeling. Owing to this property, we
can perform separation robust against the initialization of thesource model.
2.

Proposed Method

2.1

ILRMA Based on Time-Frequency-Varying Complex
Poisson Distribution

KL-ILRMA approximates the source spectrogram with anonnegative low-rank matrix by minimizing their generalized KL divergence. Here we assume that the source model
follows the time-frequency-varying complex Poisson distribution [12], which is the extension of the real-valued Poisson distribution to complex values. The probability density
function of the complex Poisson distribution is defined asp(z) =|z| −1 λ |z | −λe ,2π(|z|)!

(1)

where z is in the set D = {z ∈ C | |z| ∈ N} and λ is a shape

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

LETTER459

parameter characterizing its distribution.
In KL-ILRMA, the following time-frequency-varying
complex Poisson source generative model is assumed:Yp( yi jn ) =i, j,nY yi jnX

λ i jn | yi j n |2π yi jn !

i, j,nλ i jn =−1e−λi j n ,(2)

t ikn vk jn,

2.2

In the cost function Eq.
P (7), the only term related to Tn andH
Vn is DKL (|win
xi j | | k t ikn vk jn ). Therefore, Tn and Vn canbe optimized by minimizing the divergence via the KL-NMF
update rules [8](3)kwhere λ i jn is the shape parameter as well as the mean of
the amplitude of the complex Poisson distribution and isinterpreted as the ith row and jth column value of the source
spectrogram model as in Eq. (3). The variables t ikn andvk jn are the elements of the basis matrix Tn ∈ RI×K
≥0 andthe activation matrix Vn ∈ RK×J,
respectively,
where
R ≥0
≥0
denotes the set of nonnegative real numbers. k = 1, . . . , K
is the index of the bases, and the number of bases K is
usually much less than I and J. Note that the validity of thisgenerative model is discussed in Appendix.
H
Since yi jn = win
xi j in Eq. (2), the negative log−1
likelihood of xi j = Wi yi j is given by
X"
H
H
L=
log(|win
xi j |!) + log|win
xi j |
i, j,nH− |winxi j | log

Xt ikn vk jn +

X

k

− 2J

X

#
t ikn vk jnklog |det Wi | + const.

(4)

i
H
To express the term log(|win
xi j |!) via elementary functions
and extend the domain of the definition to an analog domain,
we apply Stirling’s approximation (its practical validity will
be assessed in Sect. 3.1)

log(|win xi j |!) ≈ |win xi j | log|win xi j | − |win xi j |.
H

H

H

H

(5)

Hence, we obtain the cost function to be minimized as
X"
H
H
H
H
J =
|win
xi j | log|win
xi j | − |win
xi j | + log|win
xi j |
i, j,n− |win xi j | log
HXt ikn vk jn +

k

− 2J

X

X

#
t ikn vk jnklog |det Wi | + const.

XH
DKL (|win
xi j | |P(6)

H
k t ikn vk jn ) + log|win xi j |

i, j,n− 2JXlog |det Wi | + const.,

t ikn ← t ikn

X |w H x |v
 X
 −1
i j k jn
vk jn ,
P in
k 0 t ik 0 n vk 0 jn
j
j

(8)

vk jn ← vk jn

 X
 −1
X |w H x |ti j ikn
t ikn .
P in
k 0 t ik 0 n vk 0 jn
i
i

(9)

Since minimization of the generalized KL divergenceis a convex problem w.r.t. either Tn or Vn [9], KL-ILRMA
is expected to more stably estimate the source spectrograms
than conventional ILRMAs, which do not involve convex
problems w.r.t. either Tn or Vn .
2.3

Update Rule for Demixing Matrix

In conventional IS-ILRMA, the demixing matrix Wi can be
updated by applying iterative projection (IP) [3], which is afast and stable optimization algorithm that can be applied toH
the sum of |win
xi j | 2 and − log|det Wi |. In KL-ILRMA, however, IP cannot be applied, i.e., the cost function Eq. (6) does
not satisfy the necessary condition for the use of IP. Instead,
we apply a majorization-minimization (MM) algorithm [13]to derive the update rule of win . First, we apply the tangent
line inequality
H
log|winxi j | ≤1H(|win
xi j | − αi jn ) + log αi jn
αi jn

(10)

H
to the term log|win
xi j | in Eq. (6), where αi jn > 0 is an
auxiliary variable. Thus, the majorization function can be
designed as
#
X" 1
H
H
|win
xi j | 2 + d i jn |win
xi j |J ≤ J1 =αi jni, j,nX− 2J
log |det Wi | + const.,
(11)

i

i

=

Update Rule for Source Model(7)iwhere DKL ( y | x) = y log y − y log x − y + x is the generalized KL divergence. Therefore, the minimization of the costfunction J simultaneously achieves high independence between the sources and the low-rank modeling of each source
spectrogram based on the generalized KL divergence.

P1
d i jn =+ log(αi jn / k t ikn vk jn ) − 2,
αi jn

(12)

H
where J and J1 become equal only when αi jn = |win
xi j |.
Second, we design the further majorization function ofH
d i jn |win
xi j | to make it differentiable w.r.t. win . As shown in
Fig. 1, we can compose a majorization function of d | y | thatcoincides with d | y | at an arbitrary point y0 (black dotted line
in Fig. 1) by branching into “paraboloid type” (d ≥ 0) and“plane type” (d < 0) cases. Thus, the majorization functionsHof d i jn |win
xi j | are obtained as

IEICE TRANS. FUNDAMENTALS, VOL.E102–A, NO.2 FEBRUARY 2019460

of win , Eq. (17) can be regarded as a function of win by
fixing the other row vectors of Wi . The partial derivative ofH
Eq. (17) w.r.t. win
is
bin
∂J2
= Uin win + rin − H.
H∂winwin binFig. 1 Shapes of majorization function (blue) of original function d | y |(green) with contact point y0 (black dotted line). When d > 0, d | y | canbe majorized by paraboloid of revolution (left). When d < 0, d | y | can be
majorized by plane (right).

d i jn H
1|win xi j | 2 + d i jn βi jn



2
β
2
i
jn




(d i jn ≥ 0)H
d i jn |win xi j | ≤ 
,


1

H
H


d
ω
w
w
x
x
+
ω

i
jn
i
jn
i
j
i
j
i
jn
in
in


2



(d
i jn < 0)
(13)

H
From ∂J2 /∂win= 0, the stationary point is given in [14] as
"
#
1

−1


U
b
−
r
(if ubr = 0)
in
in

in √


ubb

s




4ubb +
win = −1  ubr *1 −
 ,

b
−
r
U
1
+

in
in 

in  2u2|ubr |  bb ,
(otherwise)(19)
H−1 b and u
H−1
where ubb = bin
Uinin
br = bin Uin rin . Substituting
the equality condition of the auxiliary variables for Eq. (19),the update rule of win can be obtained as follows:P
1H+ log(|win
xi j |/ k t ikn vk jn ) − 2, (20)H|win
xi j |
max(0, d i jn ) +1 X* 1Uin ←+xi j xiHj , (21)H
H
J j , |win xi j |2|win
xi j | d i jn ←where ∗ is the complex conjugate of ∗, βi jn > 0 is a real
auxiliary variable, and ωi jn is a complex auxiliary variable that satisfies ωi jn = 1. The equality of Eq. (13)H
holds only when βi jn = |win
xi j | (d i jn ≥ 0) and ωi jn =HHwin xi j /|win xi j | (d i jn < 0). Applying Eq. (13) to Eq. (11),
we obtain the majorization function J2 as follows:
Xf
g
H
H
H
J1 ≤ J2 = Jwin
Uin win + win
rin + rinwin

r̃in ←− 2J

log |det Wi | + const.,(14)

max(0, d i jn )1X 1+
xi j xiHj ,J j αi jn
2 βi jn(15)!

1 X ωi jnmin(d i jn, 0) xi j .
rin =
J j 2

(16)

Since Eq. (14) contains a linear term in win , we stillcannot apply IP to Eq. (14). Instead, we apply another typeof vectorwise coordinate descent to minimize functions suchas Eq. (14). In this algorithm, we focus on win , namely,
the Hermitian transpose of the particular row vector of Wi .
Equation (14) can be transformed as follows by cofactor
expansion:
Xf
g
H
H
H
J2 = Jwin
Uin win + win
rin + rin
win
i,n

−J

X

H
log| binwin | 2 + const.,

(17)

i

where bin is the column vector of Bi = (bi1, . . . , bi N ), whichis the adjugate matrix of Wi . bin can also be written as
bin = (det Wi )Wi−1 e n , where e n is an N-dimensional vector
whose nth element is one and whose other elements are zero.
Since bin only depends on win0 (n 0 , n) and is independent

(22)

−1 −1
w̃in ← Uin
Wi e n,(23)

uww ← w̃in Uin w̃in,(24)

uwr ← w̃in Uin r̃in,w̃in− r̃in(if uwr det Wi = 0)

√


uww


s


uwr *
4uww +.
win ← 1− 1+w̃in − r̃in


2

2u

|uwr | ww ,
(otherwise)


(25)

H

i

Uin =

H
xi j
1 −1 X 1 win
Uin
min(d i jn, 0) xi j ,
H
J
2
|w
x
in i j |
j

H

i,n

X

(18)

(26)

In KL-ILRMA, the cost function Eq. (6) is minimizedby alternately repeating the update of the source spectrograms Tn and Vn using Eqs. (8) and (9), and the update ofthe demixing matrix Wi using Eqs. (20)–(26). Since all update rules in KL-ILRMA are derived by the MM algorithm,
a monotonic decrease in the cost is guaranteed.
3.

Numerical Simulation

3.1

Evaluation of KL-ILRMA with Toy Model

We confirmed that KL-ILRMA is valid for separating
sources that follow a time-frequency-varying complex Poisson distribution by using an artificial sound source model,i.e., a toy model. We created the toy model via the following
×K andthree steps. First, we generated each entry of T̃n ∈ RI≥0
K×J
Ṽn ∈ R ≥0 with independent gamma distributions. Second, we calculated the normalized low-rank matrix R n by

LETTER

461
Table 1 Music sources obtained from SiSEC2011.
Index
Source (1/2)
Impulse response
No. 1
A. guitar/vocal
IR1
No. 2
A. guitar/vocal
IR2
No. 3
A. guitar/piano
IR1
No. 4
A. guitar/piano
IR2

Fig. 2 Recording conditions of impulse responses E2A (T60 = 300 ms)
obtained from RWCP database [16]: (a) IR1 and (b) IR2.

Fig. 3 Average SDR improvements of IS-ILRMA and KL-ILRMA for
toy model separation.

dividing T̃n Ṽn by the maximum value of T̃n Ṽn . Third, we obtained the source spectrogram matrix Sn whose ith row andjth column entry is generated from the independent complex
Poisson distribution with the shape (mean) parameter λ i jn
given by each entry of cR n , where c is the scale parameter of
the toy model. In this experiment, the toy model was createdwith a size of I = 257 and J = 514. The number of bases
was K = 10, and the kurtosis of each gamma distribution
was set to 12 in T̃n and 2 in Ṽn . Note that in the timefrequency-varying complex Poisson distribution, λ i jn controls both the shape of the distribution and the spectrogram
strength, and consequently, various settings for the scale parameter c should be tested; we set c = 1, 2, 3, 4, 5, 10, 15, and20.
We produced the two-channel observed signals by convoluting the IR1 impulse response (shown in Fig. 2(a))
with each source. We used the source-to-distortion ratio(SDR) [15] as the total separation performance. The number of bases of the source model in each ILRMA was 10,
which is the same as the number of bases used to generate
the toy model, and the number of iterations was 200. The
initial demixing matrices Wi were set to the ideal value with
5% noise, and the initial source model matrices Tn and Vn
were set to uniformly distributed random values.
Figure 3 shows the SDR improvement for each method
plotted against the scale parameter c, where the plotted values are the average of 20 trials with different initial values ofTn and Vn , and the error bar represents the standard deviation.
KL-ILRMA outperforms IS-ILRMA for all values of c, andthe standard deviation of KL-ILRMA is within 0.03 dB andmuch less than that of IS-ILRMA. This result clearly implies
the robustness of KL-ILRMA against the initialization of thesource model, which is due to the convex property of KL-

Fig. 4 Original and separated spectrograms for No. 2 in KL-ILRMA. We
only show truncated spectrograms with three-second duration.

NMF. Although KL-ILRMA causes a mismatch between thecost function Eq. (6) and the log-likelihood function Eq. (4)of the generative model due to the Stirling approximation
Eq. (5), the approximation has no practical effect because theSDR improvement of KL-ILRMA is greater than or around
20 dB.

IEICE TRANS. FUNDAMENTALS, VOL.E102–A, NO.2 FEBRUARY 2019462

outperforms the other methods for No. 2–No. 4, and thestandard deviation of KL-ILRMA is the smallest among themethods for all mixed signals. This result confirms the ability of KL-ILRMA for the stable estimation of sources.
4.

Conclusion

We proposed a new BSS method, KL-ILRMA, which estimates the low-rank source model via a convex optimization
based on generalized KL divergence. We derived the update
rule in KL-ILRMA using the MM algorithm to guarantee
a monotonic decrease in the cost function. From the experimental evaluation, we confirmed the robustness of KLILRMA against the parameter initialization and its ability
for the stable estimation of sources.
References
Fig. 5 Average SDR improvements for MNMF, t-MNMF, IS-ILRMA,
t-ILRMA, and KL-ILRMA.

3.2

Evaluation of KL-ILRMA with Audio Source Separation

We compared the separation performance of KL-ILRMA
with those of other conventional methods: MNMF [5], tMNMF [10], IS-ILRMA [6], and t-ILRMA [11]. We used
the music signal bearlin-roads in SiSEC2011 [17] as thedry sources, where acoustic guitar, vocal, and piano wereused. We produced four different two-channel observed signals (No. 1–No. 4) by convoluting the IR1 and IR2 impulse
responses (respectively shown in Figs. 2(a) and (b)) with eachsource. Table 1 shows the pair of instruments and impulse
response for each source. In IS-ILRMA, t-ILRMA, and KLILRMA, the initial demixing matrices were set to identity
matrices and the entries of the initial source model matrices
were set to uniformly distributed random values. In MNMF
and t-MNMF, the initial values of the parameters were set
as in [5]. The degree of freedom parameter in t-MNMF was
set to ν = 1 and the degree of freedom parameter and thedomain parameter in t-ILRMA were set to ν = 3 and p = 2,
respectively, which were the best settings for this experiment. In t-ILRMA, we applied a tempering technique basedon [11]. The sampling frequency was 16 kHz. An STFT
was performed using a 512-ms-long Hamming window witha 128-ms-long shift. The numbers of iterations, bases, andtrials were set to 200, 10, and 20, respectively.
Figure 4 shows an example of separated spectrogramsfor No. 2, where we depict the original source spectrogramsas reference and their low-rank models estimated in KLILRMA. It is confirmed that the proposed method can approximate the source spectrograms by low-rank matrices appropriately.
Figure 5 shows the average SDR improvements for eachmethod, where the error bar represents the standard deviation. In terms of the average improvements, KL-ILRMA[1] P. Smaragdis, “Blind separation of convolved mixtures in the frequency domain,” Neurocomputing, vol.22, no.1, pp.21–34, 1998.
[2] H. Saruwatari, T. Kawamura, T. Nishikawa, A. Lee, and K. Shikano,“Blind source separation based on a fast-convergence algorithm combining ICA and beamforming,” IEEE Trans. Audio, Speech, Language Process., vol.14, no.2, pp.666–678, 2006.
[3] N. Ono, “Stable and fast update rules for independent vector analysis
based on auxiliary function technique,” Proc. WASPAA, pp.189–
192, 2011.
[4] A. Ozerov and C. Févotte, “Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation,” IEEE
Trans. Audio, Speech, Language Process., vol.18, no.3, pp.550–563,
2010.
[5] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Multichannel
extensions of non-negative matrix factorization with complex-valued
data,” IEEE Trans. Audio, Speech, Language Process., vol.21, no.5,
pp.971–982, 2013.
[6] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and H. Saruwatari,“Determined blind source separation unifying independent vector
analysis and nonnegative matrix factorization,” IEEE/ACM Trans.
Audio, Speech, Language Process., vol.24, no.9, pp.1626–1641,
2016.
[7] H. Kameoka, T. Yoshioka, M. Hamamura, J.L. Roux, and K. Kashino,“Statistical model of speech signals based on composite autoregressive system with application to blind source separation,” Proc.
LVA/ICA, pp.245–253, 2010.
[8] D.D. Lee and H.S. Seung, “Algorithms for non-negative matrix factorization,” Proc. NIPS, pp.556–562, 2000.
[9] 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,
UK, 2009.
[10] K. Kitamura, Y. Bando, K. Itoyama, and K. Yoshii, “Student’s tmultichannel nonnegative matrix factorization for blind source separation,” Proc. IWAENC, 2016.
[11] S. Mogami, D. Kitamura, Y. Mitsui, N. Takamune, H. Saruwatari,
and N. Ono, “Independent low-rank matrix analysis based on complex Student’s t-distribution for blind audio source separation,” Proc.
MLSP, 2017.
[12] H. Kameoka, “Towards a statistical audio signal processing framework based on the I-divergence,” Proc. 2011 Spring Meeting of
Acoustical Society of Japan, pp.813–814, 2011.
[13] D.R. Hunter and K. Lange, “Quantile regression via an MM algorithm,” J. Comput. Graph. Stat., vol.9, no.1, pp.60–77, 2000.
[14] Y. Mitsui, N. Takamune, D. Kitamura, H. Saruwatari, Y. Takahashi, and K. Kondo, “Vectorwise coordinate descent algorithm for

LETTER

463

spatially regularized independent low-rank matrix analysis,” Proc.
ICASSP, pp.746–750, 2018.
[15] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Trans. Audio, Speech,
Language Process., vol.14, no.4, pp.1462–1469, 2006.
[16] S. Nakamura, K. Hiyane, F. Asano, T. Nishiura, and T. Yamada,“Acoustical sound database in real environments for sound scene
understanding and hands-free speech recognition,” Proc. LREC,
pp.965–968, 2000.
[17] S. Araki, F. Nesta, E. Vincent, Z. Koldovský, G. Nolte, A. Ziehe,and A. Benichoux, “The 2011 signal separation evaluation campaign(SiSEC2011): - audio source separation -,” Proc. LVA/ICA, pp.414–
422, 2012.
[18] Y. Takahashi, T. Takatani, K. Osako, H. Saruwatari, and K. Shikano,“Blind spatial subtraction array for speech enhancement in noisy environment,” IEEE Trans. Audio, Speech, Language Process., vol.17,
no.4, pp.650–664, 2009.
[19] H. Saruwatari, Y. Ishikawa, Y. Takahashi, T. Inoue, K. Shikano, andK. Kondo, “Musical noise controllable algorithm of channelwise
spectral subtraction and adaptive beamforming based on higher order
statistics,” IEEE Trans. Audio, Speech, Language Process., vol.19,
no.6, pp.1457–1466, 2011.

Appendix:

Validity of Time-Frequency-Varying Complex Poisson Distribution ModelAs shown in Sect. 2.2, the low-rank approximation based on
the time-frequency-varying complex Poisson distribution isequivalent to KL-NMF. In this appendix, we clarify whatkind of real-world signal follows the statistical model.
We compared the accuracy for modeling various realworld signals with the time-frequency-varying complex Poisson distribution (KL-NMF) and the time-frequency-varying
complex Gaussian distribution (IS-NMF). We used threekinds of signals for evaluation: 18 music signals, 18 speech
signals [17], and 18 environmental sounds [18], [19]. We
plotted the source-to-artifact ratio (SAR) of the signals approximated via both NMFs in Fig. A· 1. We can confirm that
SAR of the signals approximated via KL-NMF is higher thanthat approximated via IS-NMF. This result indicates that thegenerative model based on the time-frequency-varying complex Poisson distribution can be valid for real-world signals,
especially for music signals.

Fig. A· 1 Box-and-whisker plot on SAR of signals approximated via ISNMF and KL-NMF: (a) music signals, (b) speech signals, and (c) environmental sounds.