EFFICIENT ALGORITHMS FOR MULTICHANNEL EXTENSIONS OF ITAKURA-SAITONONNEGATIVE MATRIX FACTORIZATION
Hiroshi Sawada, Hirokazu Kameoka, Shoko Araki, Naonori Ueda
NTT Communication Science Laboratories, NTT Corporation, Japan
ABSTRACT

2. NMF WITH IS DIVERGENCE

This paper proposes new algorithms for multichannel extensionsof nonnegative matrix factorization (NMF) with the Itakura-Saito(IS) divergence. We employ Hermitian positive definite matrices
for modeling the covariance matrix of a multivariate complex Gaussian distribution. Such matrices are basically estimated for NMF
bases, but a source separation task can be performed by introducing
variables that relate NMF bases and sources. The new algorithms
are derived by using a majorization scheme with properly designed
auxiliary functions. The algorithms are in the form of multiplicative
updates, and exhibit good convergence behavior. We have succeeded in separating a professionally produced music recording into
its vocal and guitar components.

This section reviews the formulation and algorithm of IS-NMF [4,9]. Let us assume that we have a single channel observation, andwe apply a short-time Fourier transform (STFT) to the observation.
Let xij ∈ C be the STFT coefficient at frequency i and time j. The
generative model of all the STFT coefficients X, [X]ij = xij , ofsize I × J can be written as
Q
P
Qp(X|T, V) = Ii=1 Jj=1 N (xij |0, K(1)
k=1 tik vkj )Index Terms— Nonnegative matrix factorization, Itakura-Saito
divergence, Auxiliary function, Multichannel, Source separation
1. INTRODUCTIONk

Identifying frequent spectral patterns in sounds is an important technique for various kinds of audio signal applications, including audio
scene analysis and music transcription. Nonnegative matrix factorization (NMF) [1] is widely used for such purposes (e.g., [2]). There
are many choices for the distance/divergence measure used in the
NMF cost function, such as the Euclidean distance [3], the generalized Kullback-Leibler (KL) divergence [3], and the Itakura-Saito(IS) divergence [4]. It is recognized that the IS divergence is often
preferable for modeling audio signals.
On the other hand, multichannel extensions of NMF have been
receiving attention with a view to realizing sound source separation
with multiple microphones. In [5, 6], the covariance matrices of multivariate complex Gaussian distributions are modeled with the NMF
scheme. This can be seen as a multichannel extension of IS-NMF(NMF with the IS divergence). Although expectation-maximization(EM) algorithms have been derived for source separation tasks, it
was reported that the algorithms were sensitive to the initializationof parameters. In [7], we proposed a multichannel extension of Euclidean NMF, and derived multiplicative update rules [3] based on
majorization [8] with a properly designed auxiliary function. The
algorithm converges favorably.
In this paper, we propose new algorithms for multichannel extensions of IS-NMF. The algorithms are derived from the same majorization scheme that we employed for the multichannel extensionsof Euclidean NMF [7]. We study two formulations in which a spatial property is considered for each NMF basis (the first one) and foreach source (the second one). In the second formulation, a source
separation task is performed by automatically clustering the NMF
bases for each source. The derived algorithms are in the form ofmultiplicative updates [3], and they exhibited efficient convergence
behaviors in the experiments.

978-1-4673-0046-9/12/$26.00 ©2012 IEEE

where N represents a complex Gaussian distribution, and K is thenumber of rank-1 basis matrices. Nonnegative matrices T and V,
whose elements are tik ≥ 0 and vkj ≥ 0, have sizes of I × K andK × J, respectively.
The negative log likelihood − log p(X|T, V) is given by
!
X
X
|x |2
P ij
+ log
tik vkj(2)
L(T, V) =k tik vkj
i,j261

where constant terms are omitted. Minimizing L(T, V) is equivalent to minimizing the IS divergence
P
|x |2
|x |2
P ij
− log P tij v − 1 (3)dIS (|xij |2 , K
k=1 tik vkj ) =t v
k ik kj

k ik kj

for all the STFT coefficients X.
The negative log likelihood L(T, V) can be minimized in an
iterative manner by majorization [8] with an auxiliary function
L+ (T, V, R, U) =!
P
2
X rijkk tik vkj − uij+ log uij +|xij |(4)tik vkjuiji,jkPwhere rijk and uij are auxiliary variables that satisfy k rijk = 1,
rijk > 0 and uij > 0. It can be verified [10, 9] that the auxiliary
function L+ cannot be smaller than the original function L(T, V) ≤L+ (T, V, R, U) and the equality L = L+ is satisfied whenX
tik vkj, uij =
tik vkj .
(5)rijk = Pk tik vkj
X

2kTherefore, the negative log likelihood L is indirectly minimized byrepeating the following two steps:
1. Minimizing L+ with respect to R and U by (5), which makes
L(T, V) = L+ (T, V, R, U).
2. Minimizing L+ with respect to T or V, which also minimizes L. The update rules are derived from the partial derivative of L+ with respect to the corresponding variables asv
v
u P r2 |x |2
u P r2 |x |2uijk iju
ijk ij
u j
u i
vkj
t
, vkj ← t P tik. (6)tik ← t P vkj
ik
j uij

i uij

ICASSP 2012

By substituting (5) into (6), we obtain the update rules for IS-NMF:v
v
uP
u P v |x |2
tik |xij |2
kj
ij
u
u
u j x̂ij x̂ij
u i x̂ij x̂ijtt
, vkj ← vkj(7)
tik ← tik
P vkj
P tik
j x̂ijwith x̂ij =i x̂ijP

k tik vkj .

3. PROPOSED MULTICHANNEL EXTENSIONSThis section presents our proposed extensions of IS-NMF fora multichannel case, where we have a complex-valued vector
x = [x1 , . . . , xM ]T ∈ CM for a time-frequency slot, with xm
being the m-th microphone observation. Let xij be such a vector at
frequency i and time j. Then, the generative model of all the STFT
coefficients X, [X]ij = xij ∈ CM , can be written asp(X|θ) =

QI

i=1

QJ

j=1 N (xij |0, X̂ij )(8)where N represents a multivariate complex Gaussian distribution,and X̂ij is an M × M Hermitian positive definite covariance matrix
that models the multichannel observation at frequency i and time j.
We will detail the parameter set θ and the covariance matrix X̂ij in
the next two subsections.

The minimization updates with respect to the main variables T
and V are derived from the partial derivatives of L+ (see Appendix
6.2), and given byv
uP−1
u j vkj tr(X̂−1
ij Xij X̂ij Hik )tik ← tik tP−1
j vkj tr(X̂ij Hik )

(13)

v
uP
−1
u
tik tr(X̂−1
ij Xij X̂ij Hik )vkj ← vkj t i P−1i tik tr(X̂ij Hik )(14)

We observe that these update rules reduce to (7) if M = 1, Xij =|xij |2 and Hik = 1. Thus, we understand that the set of these updates constitutes a multichannel extension of IS-NMF.
The minimization update with respect to Hik is performed bysolving an algebraic Riccati equation as shown in Appendix 6.2. The
procedure is as follows. First, two matrices
A=
B = HikP−1
j vkj X̂ij ,

“P

−1
−1
j vkj X̂ij Xij X̂ij

(15)
”

Hik

(16)

3.1. NMF basis-wise spatial property Hik

are calculated. Then, we perform an eigenvalue decomposition of a2M × 2M matrix
»
–
0
−A
(17)
−B
0

We start with a simple case where each NMF basis has its own spatial property in this multichannel scenario. Let Hik be an M × M
Hermitian positive definite matrix that models the spatial property of
the k-th NMF basis at frequency i. Then, we detail the covariancematrix asP(9)X̂ij = K
k=1 Hik tik vkj ,and let e1 , . . . , eM be eigenvectors with negative eigenvalues. It is
guaranteed that there are exactly M negative eigenvalues. Then, let
us decompose the 2M -dimensional eigenvectors as»–fem = mfor m = 1, . . . , M(18)gm

where tik and vkj are nonnegative scalars that constitute NMF bases
T and V. To fix the scaling ambiguity between Hik and tik , a constraint ||Hik ||F = 1 is introduced.
The negative log likelihood − log p(X|θ) with θ = {T, V, H}
is given by (omitting constant terms)iP hL(T, V, H) = i,j tr(Xij X̂−1
(10)ij ) + log det X̂ij

with fm and gm being M -dimensional vectors. Finally the new Hik
is calculated by
Hik ← GF−1
(19)with Xij = xij xH
ij and tr(·) being the trace of a matrix. This negative log-likelihood can be seen as a multichannel extension of (2).
Therefore, we consider the minimization of L(T, V, H) as a multichannel extension of IS-NMF.
Following the optimization framework explained in Sect. 2, we
introduce an auxiliary functionX(11)
L+ (T, V, H, R, U) =3.2. Source-wise spatial property Hio"i,j

X tr(Xij Rijk H−1 Rijk )ikktik vkj

det X̂ij − det Uij+ log det Uij +det Uij#

where Rijk and Uij are auxiliary
variables that satisfy HermitianP
positive definiteness and k Rijk = I, with I being the identity
matrix of size M . As shown in Appendix 6.1, the auxiliary function L+ cannot be smaller than the original function L(T, V, H) ≤L+ (T, V, H, R, U) and the equality L = L+ is satisfied when
Rijk = tik vkj Hik X̂−1
ij ,Uij = X̂ij .

(12)

262

with F = [f1 , . . . , fM ] and G = [g1 , . . . , gM ].
In summary, the negative log-likelihood (10) with (9) is iteratively minimized by the updates (13), (14) and (19).

When multichannel NMF is applied to a source separation task,
NMF bases originating from the same source should be clustered
together. To realize such functionality, this subsection proposes
another model for the covariance matrix
X̂ij =

PK

k=1

PN

o=1 zko Hio tik vkj(20)

by extending (9). In this model, spatial property matrices Hio depend on frequency i and source o, where o = 1, . . . , N with N
being the number of sources. To relate
P NMF bases and sources, we
introduce variables zko that satisfy N
o=1 zko = 1 and zko ≥ 0. It
is interpreted that the k-th NMF basis belongs to the o-th source if
zko is close to 1.
The negative log-likelihood L(T, V, H, Z), − log p(X|θ) with
θ = {T, V, H, Z}, in this model is literally the same as (10) but thedefinition of X̂ij is different and should follow (20). L(T, V, H, Z)
can be minimized in a similar manner to that shown in the previous
subsection.

We introduce an auxiliary functionX
L+ (T, V, H, Z, R, U) =

io

k,o

zko tik vkj+ log det Uij +3

10

det X̂ij − det Uij 5
det Uij

variables that satisfy Hermiwhere Rijko and Uij are auxiliaryPtian positive definiteness and k,o Rijko = I. It is verified that
L(T, V, H, Z) ≤ L+ (T, V, H, Z, R, U) and the equality is satisfied when
Rijko = zko tik vkj Hio X̂−1
ij ,Uij = X̂ij .

6.8

6.7

10

6

10

6.6

10

6.5

10
0

50
# iterations(22)

The minimization updates with respect to the main variables T,
V and Z are given by
v
P
uP
−1
u j vkj o zko tr(X̂−1
ij Xij X̂ij Hio )t,(23)
tik ← tik
P
P
−1
j vkj
o zko tr(X̂ij Hio )
v
P
uP
−1
u
tik
zko tr(X̂−1
ij Xij X̂ij Hio )vkj ← vkj t i P o P
,
(24)
−1
i tik
o zko tr(X̂ij Hio )v
uP−1
u i,j tik vkj tr(X̂−1
ij Xij X̂ij Hio )zko ← zko t
.
(25)
P
−1
i,j tik vkj tr(X̂ij Hio )P
Regarding zko , a normalization zko ← zko / ( o zko ) should beP
conducted after the updates to satisfy the constraint N
o=1 zko = 1.
Again, we observe that these updates reduce to (7) if N = 1, M =1, Xij = |xij |2 and Hik = 1. Thus, we understand that the set ofthese updates constitutes another multichannel extension of IS-NMF.
The minimization procedure for Hio is basically the same as thatshown in the previous subsection, except that the calculation of theA, B matrices should be performed as follows
P
A = j,k zko tik vkj X̂−1
(26)
ij ,
“P
”
−1
−1
Hio .
(27)
B = Hioj,k zko tik vkj X̂ij Xij X̂ij
4. EXPERIMENTS
We used a sound mixture available at the Signal Separation Evaluation Campaign (SiSEC 2010) [11]. More specifically, we used
dev1 tamy-que pena tanto faz snip 6 19 mix.wav
found in the development data of the Professionally produced music
recording dataset. The sound file has 2 channels (microphones) andis a mixture of 2 sources (vocal and guitar). We down-sampled the
mixture from 44.1 to 16 kHz, and applied STFT (with a 64 ms frame
size and a 16 ms frame shift) to the mixture.
Figure 1 shows the convergence behavior when we applied NMF
algorithms to the mixture. For comparison, we ran four different algorithms: 1) single-channel NMF described in Sect. 2 applied to thefirst microphone observation, 2) multichannel NMF with a basiswise spatial property proposed in Sect. 3.1 (Basis-wise H), 3) thesame multichannel NMF with a basis-wise spatial property but optimized by the EM algorithm proposed in [6] (EM algorithm), and4) multichannel NMF with a source-wise spatial property proposed
in Sect. 3.2 (Source-wise H). The number K of rank-1 basis matrices was set at 10 in all cases. For each algorithm, we performed

263

Basis−wise H
EM algorithm
Source−wise H

6.9

10

100

0

50
# iterations

100

Fig. 1. Convergence behavior
12
Average SDR (dB)4i,jX tr(Xij Rijko H−1 Rijko )Multichannel NMF(21)
IS divergence2

Single−channel NMF

10
8
6
4
2
0

100

200
300
# iterations

400

500

Fig. 2. Source separation performance

five runs starting from different random initial matrices T,√ V, and
Z. The diagonal elements of H were initially all set to 1/ M , with
M being the number of microphones, and the off-diagonal elements
were initially all set to zero.
From the result shown in Fig. 1, we observe the followings. The
convergence behavior of Basis-wise H multichannel NMF was favorable, and appears similar to that of the single-channel IS-NMF.
With Source-wise H, the convergence was slowed slightly by theconstraints caused by sharing the spatial property Hio among some
NMF bases, but its behavior was still satisfactory. If we compare theconvergence of Basis-wise H and EM algorithm, the effectivenessof the proposed algorithms over the conventional EM algorithms becomes clear. The algorithms were coded with Matlab and run on
an Intel Core i7 965 (3.2 GHz) processor. The average computation time for 100 iterations by Basis-wise H, EM algorithm andSource-wise H were about 138, 267 and 89 seconds, respectively.
Figure 2 shows the source separation performance in terms of
the average signal-to-distortion ratio (SDR) [11]. The separation was
performed by the multichannel NMF (source-wise Hio ) proposed in
Sect. 3.2. We obtained the separated signals by Wiener filteringyijo =“PK

k=1 zko tik vkj”

Hio X̂−1
ij xij(28)

for each source o and every time-frequency j, i slot. A fairly highSDR of around 11 dB was attained with the proposed multichannel
NMF, starting from randomly initialized parameters.
Figure 3 shows the behavior of variables zko in the source separation task. It can be seen that randomly initialized variables zko
approached 1 or 0 as the algorithm converged. From the final statusof zko , we understood that six NMF bases constituted the first output (vocal) and the remaining four NMF bases mainly constituted
the second output (guitar).

Setting these derivatives to zero and solving them with respect to thecorresponding variables yields update rules that are similar to (6).
Then, by substituting Rijk and Uij with (12), we have the updates(13) and (14).
The partial derivative with respect to Hik is»–
P det X̂ij −1
H−1
R
X R
H−1
∂L+
ik ijk ij ijk ik
=
t
v
−
X̂
. (34)
ik
kj
ij
j
∂Hik
det Uij
tik vkj

1

z value

0.8

zk1

0.6

zk2

0.4
0.2
0
0

100

200
300
# iterations

400

Setting the derivative to zero, solving it with respect to Hik , andsubstituting Rijk and Uij with (12) yields

500

Fig. 3. Behavior of variables zko , k = 1, . . . , 10, o = 1, 2

Hik AHik = B

5. CONCLUSION

where matrices A and B are defined in (15) and (16). Equation (35)
is a special case of an algebraic Riccati equation, and it can be solved
by using the procedure described in the main text.

To extend IS-NMF to a multichannel case, we have designed an auxiliary function (11) and derived efficient multiplicative update rules(13), (14) and (19). We then proposed a second covariance model(20) by introducing variables zko that relate NMF bases tik vkj anda source specific spatial property Hio . The model forces the learned
NMF bases to cluster into the sources, and contributes to solving asource separation task. Experimental results show that the proposed
algorithms converged favorably, and source separation was accomplished with a high SDR value. Future work will include intensive
experiments with a variety of sound sources.
6. APPENDIX
6.1. Proof for the auxiliary function condition of (11)+Let us consider the minimization
Pof L defined in (11) with respect
to R subject to the constraint k Rijk = I, and also with respect
to U. By introducing Lagrange multipliers Λij of size M × M , we
have
ˆ P
˜
P
(29)
F = L+ + ij tr ( k Rijk − I)H Λij .

The partial derivative of F with respect to R∗ijk is given as
∂F
∂R∗
ijk

= 2(Hik tik vkj )−1 Rijk Xij + Λij .

(30)

Setting this to zero gives 2Rijk = −tik vkj Hik Λij X−1
ij and addingX.
Byeliminatingthese for k = 1, . . . , K gives Λij = −2X̂−1ijij−1Λij , we have Rijk = tik vkj Hik X̂ij .
The partial derivative of L+ with respect to U∗ij is given as
det X̂ij −1
∂L+
= U−1
ij − det Uij Uij .
∂U∗
ij

(31)

Setting this to zero gives Uij = X̂ij .
In summary, the minimum L+ is obtained by (12) and the minimum value is equal to L defined in (10).
6.2. Derivation of update rules
The partial derivatives of L+ defined in (11) with respect to tik andvkj are»–P det X̂ij
tr(Rijk H−1
R
X )
−1
∂L+
ik ijk ij
=
tr(
X̂
H
)v
−
, (32)
ik
kj
2
ij
j
∂t
det Uij
t v
ik

∂L+
∂vkj

ik kj

=

P
i

»

tr(Rijk H−1RX )det X̂ij
ik ijk ijtr(X̂−1
ij Hik )tik −det Uij
tik v 2kj–
. (33)

264(35)

7. REFERENCES[1] D. D. Lee and H. S. Seung, “Learning the parts of objects
with nonnegative matrix factorization,” Nature, vol. 401, pp.
788–791, 1999.
[2] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization for polyphonic music transcription,” in Proc. WASPAA
2003, Oct. 2003, pp. 177–180.
[3] D.D. Lee and H.S. Seung, “Algorithms for non-negative matrix
factorization,” in Advances in Neural Information Processing
Systems, 2001, vol. 13, pp. 556–562.
[4] 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.
[5] A. Ozerov and C. Févotte, “Multichannel nonnegative matrix
factorization in convolutive mixtures for audio source separation,” IEEE Trans. Audio, Speech and Language Processing,vol. 18, no. 3, pp. 550–563, Mar. 2010.
[6] S. Arberet, A. Ozerov, N.Q.K. Duong, E. Vincent, R. Gribonval, F Bimbot, and P. Vandergheynst, “Nonnegative matrix factorization and spatial covariance model for under-determined
reverberant audio source separation,” in Proc. ISSPA 2010,
May 2010, pp. 1–4.
[7] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “New formulations and efficient algorithms for multichannel NMF,” inProc. WASPAA 2011, Oct. 2011, pp. 153–156.
[8] J. de Leeuw, “Block-relaxation methods in statistics,” in Information Systems and Data Analysis, H. H. Bock, W. Lenski,
and M. M. Richter, Eds., pp. 308–324. Springer Verlag, 1994.
[9] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono,and S. Sagayama, “Convergence-guaranteed multiplicativealgorithms for non-negative matrix factorization with betadivergence,” in Proc. MLSP 2010, Aug. 2010, pp. 283–288.
[10] H. Kameoka, M. Goto, and S. Sagayama, “Selective amplifierof periodic and non-periodic components in concurrent audio
signals with spectral control envelopes,” in IPSJ SIG Technical
Reports, 2006-MUS-66-13, Aug. 2006, pp. 77–84, in Japanese.
[11] S. Araki, A. Ozerov, V. Gowreesunker, H. Sawada, F. Theis,
G. Nolte, D. Lutter, and N. Duong, “The 2010 signal separation
evaluation campaign (SiSEC2010): Audio source separation,”in Latent Variable Analysis and Signal Separation, vol. 6365
of Lecture Notes in Computer Science, pp. 114–122. 2010.