2011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

October 16-19, 2011, New Paltz, NY

NEW FORMULATIONS AND EFFICIENT ALGORITHMS FOR MULTICHANNEL NMF
Hiroshi Sawada, Hirokazu Kameoka, Shoko Araki, Naonori Ueda
NTT Communication Science Laboratories, NTT Corporation, Japan
ABSTRACTThis paper proposes new formulations and algorithms for a multichannel extension of nonnegative matrix factorization (NMF), intending convolutive sound source separation with multiple microphones. The proposed formulation employs Hermitian positive
semidefinite matrices to represent a multichannel version of nonnegative elements. 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. Efficient optimization algorithms in the form of multiplicative updates are derived byusing properly designed auxiliary functions. Experimental results
show that two instrumental sounds coming from different directions
were successfully separated by the proposed algorithm.
Index Terms— nonnegative matrix factorization, multichannel, positive semidefinite, auxiliary function, source separation
1. INTRODUCTION
Identifying frequent patterns in recorded sounds is an important task
in various kinds of audio signal applications, including sound separation and music transcription. Nonnegative matrix factorization(NMF) [1], which was originally applied to learning parts-based
representations of images, is also a widely-used technique for suchaudio applications (e.g., [2]). When we analyze a recorded audio
signal with NMF, the signal is typically transformed into the representation xij in the time-frequency domain by a short-time Fouriertransform (STFT), where i and j represent frequency and time, respectively. Then the amplitudes |xij | are used to construct a nonnegative matrix for the analysis.
Recently, multichannel extensions of NMF have been receiving attention with a view of realizing sound source localization andseparation with multiple microphones, and have been studied for
instantaneous mixtures [3, 4] and for convolutive mixtures [5, 6].
In a convolutive case, the phase difference information on timefrequency observations between different microphones is crucial for
distinguishing/identifying source directions. Therefore, complexvalued raw STFT observations should be taken into account. In[5, 6], the covariance matrices of multivariate complex Gaussian
distributions are modeled with the NMF scheme, and EM algorithms are derived for a source separation task. However, the algorithms are sensitive to initialization, and thus they used the original
source information for perturbed oracle initializations.
In this paper, we study new types of multichannel NMF formulations for a convolutive mixing case. The new formulations are
based on modeling the means of complex Gaussian distributions,and can be seen as natural extensions of standard Euclidean NMF.
We study two formulations in which Hermitian positive semidefinite matrices H representing spatial properties are considered foreach NMF basis (the first one) and for each source (the secondone). In the second formulation, a source separation task can be per-978-1-4577-0693-6/11/$26.00 c 2011 IEEE

formed by automatically clustering the NMF bases to each source.
This constitutes a considerable advance on our previous study [7].
The derived algorithms are in the form of multiplicative updates [1].
Experimental results show that the algorithms converge favorably.
2. OPTIMIZATION FRAMEWORK FOR NMF
This section reviews the formulation and optimization frameworkof standard NMF [1]. This framework can handle complex values[8] and is important for understanding our multichannel extensionsof NMF explained in Sect. 3. The generative model of STFT amplitudes A, [A]ij = aij = |xij |, in standard NMF can be written
as
Q
P
Q
1
p(A|T, V) = Ii=1 Jj=1 N (aij | K(1)
k=1 tik vkj , 2 )where N represents a Gaussian distribution, and K is the numberof rank-1 basis matrices. Nonnegative matrices T and V, whose
elements are tik ≥ 0 and vkj ≥ 0, have sizes of I × K and K × J,
respectively.
The negative log likelihood of (1)P
P
L(T, V) = − log p(A|T, V) = i,j |aij − k tik vkj |2 , (2)where constant terms are omitted, can be minimized in an iterative
manner by majorization [9] with an auxiliary function
P P
2
1
(3)
L+ (T, V, S) = i,j K
k=1 rijk |sijk − tik vkj |P
where sijk are auxiliary variables
P that satisfy k sijk = aij , andrijk are parameters that satisfy k rijk = 1 and rijk > 0. It canbe verified that the auxiliary function L+ has two properties:
L(T, V) ≤ L+ (T, V, S) ,L(T, V) = minS L+ (T, V, S) .

(4)
(5)

Therefore, the negative log likelihood L is indirectly minimized byrepeating the following updates, each of which corresponds to theminimization of the auxiliary function L+P
1. with respect to S: sijk ← tik vkj + rijk (aij − k tik vkj )
P

â

s

2. with respect to T: tik ← Pj âijij vijk
kjj

3. with respect to S: sijk ← tik vkj + rijk (aij −

4.

P
â s
with respect to V: vkj ← Pi âijij tijk
ik
i

.

P

k tik vkj )These updates are derived from the partial derivative of L+ withrespect to the corresponding variables. If we specify the parameters
rijk asP
t v
rijk = ikâijkj , âij = k tik vkj ,(6)then the update for the auxiliary variables is simplified as sijk ←aij tik vkj /âij . By substituting this in the T and V updates, we
obtain the well-known multiplicative updates for NMF [1]:PP
aij tikj aij vkj
, vkj ← vkj Pi.
(7)tik ← tik P
j âij vkj
i âij tik

1532011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

3. FORMULATIONS AND ALGORITHMS
3.1. Nonnegativity in Multichannel NMF
In a multichannel case, we have a complex-valued vector x =[x1 , . . . , xM ]T ∈ CM for a time and frequency slot, where xm
is the m-th microphone observation. Given that, here we consider
how we guarantee the nonnegativity utilized in the NMF scheme.
A trivial way of guaranteeing nonnegativity is to take the absolute value of each element as [|x1 |, . . . , |xM |]T . However, thisignores the phase information that reflects the directivity of sourcesfrom the microphone array. Another way is to select a reference
channel (e.g., the 1st channel) and to force the corresponding element to be nonnegative as x ← x · sign(x∗1 ), where sign(x) =x/|x| and ·∗ is the complex conjugate. Although this worked wellin our previous study [7], all channels are not treated evenly and theresults may be affected by the reference channel selection.
This paper proposes using the outer product3
2
1. . . |x1 xM | 2 sign(x1 x∗M )
|x1 |
7
6
..
..
..
7 (8)
X=6
.
5
4
.
.
1
∗
|xM |
|xM x1 | 2 sign(xM x1 ) . . .

= x̄x̄H of the amplitude square-rooted observation vector x̄ =1
1[|x1 | 2 sign(x1 ), . . . , |xM | 2 sign(xM )]T . Note that [X]mn contains
information regarding the phase difference between the m-th andn-th microphones. We then introduce an M × M Hermitian positive semidefinite matrix Hk to model the spatial property of the k-th
NMF basis, and describe the outer product asp(X|θ)=

M
M Y
Y

m=1 n=1∝

Nc ([X]mn |exp(−||X −

X
k

K
X

[Hk ]mn tk vk , 1)

(9)

k=1

Hk tk vk ||2F )(10)with θ = {t1 , . . . , tK , v1 , . . . , vK , H1 , . . . , HK }, where Nc represents a complex Gaussianm and n are microphone
P distribution,PM2indexes, and ||B||2F = M
m=1
n=1 |bmn | is the squared Frobenius norm of a matrix B. Formulation (10) allows us to naturally
extend the Euclidean cost function of the standard NMF to multichannel cases, and to derive effective multiplicative updates withfavorable convergence behaviors, as shown below.
The next two subsections present multichannel NMF formulations and their corresponding algorithms. We start with a simple
case where Hermitian positive semidefinite matrices are assigned
to each NMF basis. We then consider a practical case for source
separation by sharing such Hermitian matrices within a source.

October 16-19, 2011, New Paltz, NY

Here we introduce a constraint ||Hik ||F = 1 to fix the scaling ambiguity between Hik and tik , and also to simplify the update rules
to minimize L as shown below.
Following the optimization framework explained in Sect. 2, we
introduce an auxiliary functionP
1
L+ (T, V, H, S) = i,j,k rijk||Sijk − Hik tik vkj ||2F
(13)

P
where Sijk are auxiliary variables
P that satisfy k Sijk = Xij , andrijk are parameters that satisfy k rijk = 1 and rijk > 0. As shownin Appendix 6.1, the auxiliary function L+ satisfies two conditions
likewise (4) and (5), and the equality L+ = L is satisfied whenP
Sijk = Hik tik vkj + rijk Eij , Eij = Xij − k Hik tik vkj . (14)The minimization updates with respect to the main variables
T, V and H are derived from the partial derivatives of L+ (seeAppendix 6.2), and given with the definition (6) of âij byP
j âij tr(Sijk Hik )
P
tik ←
,
(15)
j âij vkjPij tr(Sijk Hik )
i âP
,
(16)
vkj ←
i âij tik
P(17)
Hik ← j âij Sijk , Hik ← Hik /||Hik || ,PM
where tr(B) =m=1 bmm is the matrix trace. It is not ensured
that the update for Hik keeps it positive semidefinite, and it sometimes has nonnegative eigenvalues after the update. If so, we apply
eigenvalue decomposition as Hik = UDUH , set all the negative
elements of D to zero, and update Hik ← UDUH with the new D.
The above T and V updates (15) and (16) are not in the form ofmultiplicative updates [1] and are not very efficient. If we utilize arelation âij tr(Sijk Hik ) = tik vkj [âij + tr(Eij Hik )], these updates
become multiplicative as shown below. In summary, the negative
log-likelihood (12) is iteratively minimized by repeating (18) and
one of the three updates (19)-(21).
PP
Eij = Xij − k Hik tik vkj , âij = k tik vkj
(18)
"
#
P
j vkj tr(Eij Hik )P(19)
tik ← tik · 1 +j vkj âijP»–
tik tr(Eij Hik )vkj ← vkj · 1 + i P(20)
i tik âij
PP
Hik ← Hik j âij vkj + j Eij vkj , Hik ← Hik /||Hik || (21)3.3. Source-wise spatial property Hio

3.2. NMF basis-wise spatial property Hik
This subsection studies a multichannel NMF formulation in which
spatial property Hik depends on frequency i and NMF basis k, anddoes not depend on time j.
Let xij be M -channel observations at frequency i and time j,and let Xij be the corresponding outer product defined by (8). Following (10) and the NMF scheme, the likelihood of the total observations X = {X11 , . . . , X1J , . . . , XIJ } is given as
Q
P
p(X|T, V, H) ∝ i,j exp(−||Xij − k Hik tik vkj ||2F ) , (11)and the negative log-likelihood is given byP
P
L(T, V, H) = i,j ||Xij − k Hik tik vkj ||2F .

(12)

When multichannel NMF is applied to a source separation task, it isdesired that NMF bases originating from the same source are clustered together. To realize such functionality, this subsection studies another formulation in which spatial property Hio depends onfrequency i and source o, where o = 1, . . . , N with N being thenumber of sources.
The generative model is a modification of (11) by introducing
variables zko , which relates NMF basis k and source o:PQp(X|T, V, H, Z) ∝ i,j exp(−||Xij − k,o zko Hio tik vkj ||2F )
(22)
P
z
=
1.
It
is
interpreted
that
the
k-th
with zko ≥ 0 and N
koo=1
NMF basis belongs to the o-th source if zko is close to 1.

154

2011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

Along with discussions similar to those in the previous subsection, we obtain the following results.
The negative log likelihood to be minimized is given byP
P(23)
L(T, V, H, Z) = i,j ||Xij − k,o zko Hio tik vkj ||2F

Sijko = zko Hio tik vkj + rijko Eij ,P
Eij = Xij − k,o zko Hio tik vkj .

(25)
(26)

If we specify the parameters rijko asPPz t vkjrijko = koâik, âij = k,o zko tik vkj = k tik vkj , (27)ijmultiplicative updates for T, V and Z are obtained as
"
#
P
P
j vkj
o zko tr(Eij Hio )P
tik ← tik · 1 +,j vkj âij

PP–»
tik o zko tr(Eij Hio )P,
vkj ← vkj · 1 + i
i tik âijPzko ← zko i,j tik vkj [âij + tr(Eij Hio )] .

(28)

(29)
P

(30)

Regarding zko , a normalization procedure zko ← zko / o zko
Pshould be conducted to satisfy the constraint N
o=1 zko = 1. The
update rule for H is obtained asPPPP
Hio ← Hio k zko tik j âij vkj + k zko tik j Eij vkj . (31)After this update, the procedure for ensuring positive semidefiniteness, described just after (17), and Hio ← Hio /||Hio || must be performed. In summary, the negative log-likelihood(23) is iterativelyPminimized by repeating (26) with âij = k tik vkj and one of thefour updates (28)-(31).
4. EXPERIMENTS
We considered a situation where two 6-second sounds of instruments coming from different directions (70◦ and 150◦ ) were mixed
and observed at two omnidirectional microphones 4 cm apart. The
mixtures were made by convolving measured impulse responses,
whose reverberation time was 130 ms, and the sound sources. The
sampling rate was 16 kHz. Then, we applied STFT (with a 64 ms
frame size and a 16 ms frame shift) to the mixtures and obtained
I × J = 513 × 372 outer products (8).
Figure 1 shows the convergence behavior when minimizing thenegative log-likelihoods of NMF basis-wise (12) and source-wise(23). In both cases, the number K of rank-1 basis matrices was set
at 10, T, V, and Z matrices were initialized randomly (but with thesame random numbers for both runs), and H elements were all set
to 1/M , with M being the number of microphones. We see that thenegative log-likelihoods were minimized efficiently. The program
was coded in Matlab and run on an Intel Core i7 965 (3.2-GHz)processor. The computation times for 500 iterations were 83 and 49
seconds, respectively.

4.9

Negative log−likelihoodand its auxiliary function can be defined byP
1||Sijko − zko Hio tik vkj ||2F
L+ (T, V, H, Z, S) = i,j,k,o rijko(24)
P
where Sijko are auxiliary variables that satisfy k,o Sijko = Xij ,P
and rijko are parameters that satisfy k,o rijko = 1 and rijko > 0.
The equality L+ = L is satisfied when

October 16-19, 2011, New Paltz, NY

10

NMF basis−wise H
Source−wise H

4.6

10

4.3

10

0

100

200
300
# iterations

400

500

Figure 1: Convergence behavior
In the case (23) of modeling source-wise spatial properties Hio ,a source separation task was performed. Figure 2 shows spectrograms of two separated signals. They were constructed with learned
matrices T, V, and Z asP
k zko tik vkj
yijo = xij P(32)
k tik vkj

for each source o and every time-frequency j, i slot. Figure 3 shows
the relationship zko between NMF bases k and sources o. From thefigure, it is interpreted that the 2nd basis contributes to the blue output and the remaining bases contribute mainly to the red output. It isan advantage of the proposed method that the source-wise number
of NMF bases does not need to be specified but only the total number needs to be specified. As a result of the source separation, thesignal-to-distortion ratios [10] were improved from (−7.37, 7.37)
dB at the mixture to (4.83, 12.21) dB.
Figure 4 shows the learned spatial properties, or more specifically the arguments of [Hio ]12 , which represent the phase difference between two microphones of each source component o at eachfrequency i. We observe that the blue and red dots correspond to
the sources coming from 70◦ and 150◦ , respectively (the anechoicmodel lines of these two directions are drawn for reference). Thisway, the multichannel NMF also contributes to source localization.
5. CONCLUSION
We proposed new formulations and algorithms for multichannel
NMF. We start with a simple case where each NMF basis k has
its own spatial properties Hik . We then introduced a constraint thatsuch properties Hio should belong to each source o. The constraint
forces the learned NMF bases to be clustered into the sources, andcontributes to solving a source separation task with the multichannel NMF framework. Experimental results show the effectivenessof the proposed algorithms. Future work will include intensive verification of the algorithms with a variety of sound sources.
6. APPENDIX
6.1. Proof for the auxiliary function condition of (13)+Let us consider the minimization
P of L defined in (13) with respect
to S subject to the constraint k Sijk = Xij . Introducing Lagrange
multipliers λij , we havePP
F = L+ + ij λij ( k S∗ijk − X∗ij ) .
(33)

The partial derivatives of F with respect to S are given as
1∂F=(Sijk − Hik tik vkj ) + λij .
∂S∗ijk
rijk

155

(34)

2011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

October 16-19, 2011, New Paltz, NY

3
2

400
300

1
Argument

Frequency bin index (i)

y1
500

200
100

0
−1

50

100

150

200

250

300

350
−2

Frequency bin index (i)

y2
500

−3
100

400

200
300
Frequency bin index (i)

400

500

300

Figure 4: Arguments of learned [Hio ]12 .

200
100
50

100

150
200
250
Time frame index (j)

300

350

Figure 2: Two separated signals reconstructed by learned T, V matrices weighted with Z (see Eq. 32 for details).
1

0.5

0

2

4

6
Basis index (k)810Figure 3: Learned values of zko .
rijk = 0 and
Setting this to zero gives Sijk − Hik tik vkj + λijP
adding these for k = 1, . . . , K gives −λij = Xij − k Hik tik vkj .
Therefore the minimum L+ is obtained whenP
Sijk = Hik tik vkj + rijk (Xij − k Hik tik vkj ) ,(35)
and the minimum value is equal to L defined in (12).
6.2. Derivation of the update rules
Utilizing ||B||2F = tr(BH B), the auxiliary function L+ defined in(13) can be simplified as
˜
ˆ
P
2
1
− 2tik vkj tr(Sijk Hik )||Sijk ||2F + t2ik vkj
L+ = i,j,k rijk

where the constraint ||Hik ||F = 1 and Hermitian identities Hik =H
HH
ik and Sijk = Sijk are utilized. Then, the partial derivatives withrespect to T and V are
˜
P 2 ˆ
2
∂L+
= j rijk− vkj tr(Sijk Hik ) ,(36)
tik vkj
∂tik
∂L+
∂vkj

=

P

2
i rijk˜
ˆ2tik vkj − tik tr(Sijk Hik ) .

(37)

Setting these derivatives to zero and substituting rijk with (6) yields
the updates (15) and (16). The partial derivatives w.r.t. H areP∂L+= j r 1 (Sijk − Hik tik vkj )(−tik vkj ) .
(38)
∂H∗
ik

Setting the derivative to zero and substituting rijk with (6), together
with the constraint ||Hik ||F = 1, yields the update (17).
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. IEEE
Workshop Applicat. Signal Process. Audio Acoust. (WASPAA),
Oct. 2003, pp. 177–180.
[3] D. FitzGerald, M. Cranitch, and E. Coyle, “Non-negative tensor factorisation for sound source separation,” in Proc. Irish
Signals Syst. Conf., Sept. 2005, pp. 8–12.
[4] R. M. Parry and I. A. Essa, “Estimating the spatial position ofspectral components in audio,” in Proc. ICA 2006. Springer,
Mar. 2006, pp. 666–673.
[5] A. Ozerov and C. Fevotte, “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. 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, “Formulations and algorithms for multichannel complex NMF,” inProc. ICASSP 2011, May 2011, pp. 229–232.
[8] H. Kameoka, N. Ono, K. Kashino, and S. Sagayama, “Complex NMF: A new sparse representation for acoustic signals,”in Proc. ICASSP 2009, Apr. 2009, pp. 3437–3440.
[9] J. de Leeuw, “Block-relaxation methods in statistics,” in Information Systems and Data Analysis, H. H. Bock, W. Lenski,
and M. M. Richter, Eds. Springer Verlag, 1994, pp. 308–324.
[10] E. Vincent, H. Sawada, P. Bofill, S. Makino, andJ. Rosca, “First stereo audio source separation evaluation
campaign: Data, algorithms and results,” in Proc. ICA
2007, 2007, pp. 552–559. [Online]. Available: http://www.irisa.fr/metiss/SASSEC07/

ijk

156