FORMULATIONS AND ALGORITHMS FOR MULTICHANNEL COMPLEX NMF
Hiroshi Sawada, Hirokazu Kameoka, Shoko Araki, Naonori Ueda
NTT Communication Science Laboratories, NTT Corporation, Japan
ABSTRACTThis paper studies some formulations and algorithms for the multichannel extension of nonnegative matrix factorization (NMF). We
model the inter-channel characteristics of each NMF basis, including
both the amplitude ratios and the phase differences on a channel pair.
The learned inter-channel characteristics provide useful information
for binding each NMF basis to each source component in such asituation that multiple sources are mixed in a convolutive mannerand observed at multiple microphones. Effective optimization algorithms based on majorization are derived by using properly designed
auxiliary functions. Experimental results show that the algorithms
converged favorably regardless of the initialization.
Index Terms— Nonnegative matrix factorization, Complex
phase, Inter-channel characteristic, Majorization, Auxiliary function
1. INTRODUCTION
Extracting or identifying frequent sound 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 and to the analysis of documents [2], is also a widely-used technique for such audio applications (e.g., [3]). When we analyze a time-domain single-channel
audio signal x(t) with NMF, the signal is typically transformed into
the representation xij in the time-frequency domain by a short-time
Fourier transform (STFT), where i and j represent time and frequency, respectively. Then the amplitude |xij | is used to construct anonnegative matrix for the NMF analysis.
Some extensions of NMF have been proposed to improve thecapabilities for audio applications. Complex NMF [4] takes accountof the phase information xij /|xij |, which is discarded with standard NMF. It improves the analysis accuracy by faithfully obeying
the mixing process in the complex domain. Another avenue for extension is to take advantage of multichannel recordings (typicallystereo). Research effort has been devoted to considering an instantaneous mixing case [5, 6] and a convolutive mixing case [7]. In
this paper, we study a convolutive mixing case, with particular focus
on extracting the inter-channel characteristics, including the phase
information, of the NMF decomposed components.
2. PRELIMINARIES
2.1. NMF and Complex NMF

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.
If we preserve the phase information of X, complex NMF [4]can be applied. The generative model can be written as
Q
P
p(X|T, V, G) = i,j Nc (xij | K(2)
k=1 gijk tik vkj , 1)
where Nc is a complex Gaussian distribution, and G is a third-ordertensor whose (i, j, k)-element is a complex gijk ∈ C with a unit
amplitude constraint |gijk | = 1.
2.2. Iterative optimization
The negative log likelihood of (1) or (2) can be minimized in an iterative manner by majorization [8] with a properly designed auxiliary
function. This section briefly explains the procedure for standard
NMF. A similar procedure can be derived for complex NMF as presented in [4].
The negative log-likelihood of (1) is given byPP
2
L(T, V) = i,j |yij − K(3)
k=1 tik vkj | .
Now, we introduce an auxiliary function
»
–
2
P
P
P t2 vkj
2
L+ (T, V, R) = i,j yij(4)− 2yij ( k tik vkj ) + k ik
rijkthat satisfies
1. L(T, V) ≤ L+ (T, V, R) ,2. L(T, V) = minR L+ (T, V, R) .
These are guaranteed according to the following equation, which isderived from Jensen’s inequality“P”2 P(tik vkj )2
K
≤ K
(5)
k=1 tik vkj
k=1
r
ijk

PK

where rijk ≥ 0 and k=1 rijk = 1.
The negative log likelihood L is indirectly minimized by repeating the following updates, each of which corresponds to the minimization of the auxiliary function L+t

Let us start with a single channel case before moving on to a multichannel case in Sect. 3. Let x(t) be a microphone observation for a
specific time duration, say t = 1, . . . , T . By applying an STFT to
the observation, we have a complex matrix X whose (i, j)-element
is denoted by xij ∈ C. Again, i and j represent time and frequency,
respectively. Let the matrix size be I × J.

978-1-4577-0539-7/11/$26.00 ©2011 IEEE

When we employ the standard NMF [1] for the analysis, we only
consider the amplitudes Y = |X|, or equivalently yij = |xij |, ∀ i, j.
The generative model can be written as
QP
1p(Y|T, V) = i,j N (yij | K(1)
k=1 tik vkj , 2 )229v

1. with respect to R: rijk ← P ikt kjv
P

k ik kjyv

2. with respect to T: tik ← P jv2ij/rkj
j

3.
4.

kj

ijk

t v
with respect to R: rijk ← P iktikkjvkj
kP
y twith respect to V: vkj ← P it2 ij/rik
ijk
i ik

.

ICASSP 2011

The R update is derived from the equality condition of the Jensen’s
inequality, and the T and V updates are derived from the partial
derivative of L+ with respect to tik or vkj . By substituting rijk
in the T and V updates and rearranging them, we obtain the wellknown multiplicative updates for NMF [1].
3. MULTICHANNEL COMPLEX NMF

3.1.3. Update rules

This section proposes two types of multichannel extensions of NMF.
One is from complex NMF and the other is from standard NMF. For
the sake of simplicity, we consider a stereo case here. Generalization
to more channels is straightforward. Let X(1) and X(2) be complex
matrices representing time-frequency domain observations of the 1st
and 2nd microphone. Both extensions basically model multichannel
observations with common nonnegative matrices T and V, together
with a newly introduced complex matrix H that represents the interchannel characteristics corresponding to the matrix T.
3.1. Extension from Complex NMF

Let us start with the extension from complex NMF. Without lossof generality, let the 1st channel be the reference channel, where all
the inter-channel characteristics are normalized to one. Then, the 1st
channel is modeled exactly by complex NMF (2):Q
(1) Pp(X(1) |T, V, G) = i,j Nc (xij | k gijk tik vkj , 1) . (6)The 2nd channel is modeled by introducing an I × K matrix H
whose elements are complex hik ∈ C and by using the same
T, V, G matrices commonly:Q(2) Pp(X(2) |T, V, G, H) = i,j Nc (xij | k gijk hik tik vkj , 1) .
(7)
It should be emphasized that the inter-channel characteristic hik is acomplex number and thus contains information on the phase difference as well as the amplitude ratio between the 1st and 2nd channels.
When we use a microphone array whose microphone inter-distances
are not very large, e.g., 4 cm, the phase information is more important than the amplitude information for distinguishing sources coming from different directions. The introduction of hik in this manneris the main contribution of this paper over the previous work on multichannel NMF [5, 6, 7].
3.1.2. Negative log-likelihood and its auxiliary function
What we want to maximize is the multichannel likelihoodp(X(1) , X(2) |T, V, G, H) =p(X|T, V, G) · p(X|T, V, G, H)

(8)

(9)
Lmc (T, V, G, H) =
!
X
X
X
(1)
(2)
gijk tik vkj |2 + |xij −gijk hik tik vkj |2 .
|xij −kk

To follow the iterative optimization procedure explained in Sect. 2.2,
we define an auxiliary function(1), S(2) ) =L+
mc (T, V, G, H, S

X
i,j

2
X |s(1)
ijk − gijk tik vkj |
k

rijk

+

sition of xij mainly with the k-th component gijk tik vkj plus aP(1)
fraction of the error term xij − k gijk tik vkj .
Regardingthe parameter rijk , it can be specified arbitrarily asPlong as k rijk = 1 and rijk ≥ 0. However, we specify themusing a specific formtik vkjtik vkjrijk = P
=
t
v
x̂ij
ik
kj
k

(14)

that was suggested in [4].
The minimization updates for L+
mc with respect to T, V, H, G
are derived from the partial derivatives of L+
mc (see appendix), andgiven by
P
(1)
(2) ∗
∗
j x̂ij Re[gijk (sijk + sijk hik )]
P
tik =
(15)
(1 + |hik |2 ) j x̂ij vkj
P
(1)
(2) ∗
∗
i x̂ij Re[gijk (sijk + sijk hik )]
P
(16)
vkj =
2
i x̂ij tik (1 + |hik | )P(2) ∗j x̂ij sijk gijk
P
(17)
hik =
tik j x̂ij vkj
gijk =

(1)

(2)

(1)

(2)

sijk + sijk h∗ik
|sijk + sijk h∗ik |(18)

P
where x̂ij = k tik vkj , and Re[·] extracts the real part.
In summary, the negative log-likelihood Lmc is iteratively minimized by repeating (12)-(13) and one of the four updates (15)-(18).
3.1.4. Discussion(2)

whose negative log-likelihood is given byi,j

The negative log-likelihood Lmc is minimized in the same manneras that explained in Sect. 2.2. The minimization updates for the aux(1)iliary function L+
and S(2) are given by
mc with respect to S
P
(1)
(1)
(12)
sijk = gijk tik vkj + rijk (xij − k gijk tik vkj )
P
(2)
(2)
sijk = gijk hik tik vkj + rijk (xij − k gijk hik tik vkj ) (13)
It is easily verified that with the above assignment the equality(1), S(2) ) and condition
Lmc (T, V, G, H) = L+
mc (T, V, G, H, S
(1)(11) holds. We can interpret this to mean that sijk is a decompo(1)3.1.1. Model(1)

where S(1) and S(2) are newly introduced tensors of size I × J × K(1)(2)
whose elements are complex sijk , sijk ∈ C satisfying
P (1)
P (2)
(1)
(2)
(11)
k sijk = xij ,k sijk = xij
Pand rijk , ∀ i, j, k are parameters satisfying k rijk = 1, rijk ≥ 0.

(10)
!

2
X |s(2)
ijk − gijk hik tik vkj |
k

rijk

230

In our preliminary experiments, the multichannel extension of complex NMF proposed here did not work very well, in the sense thathik did not successfully reflect the inter-channel characteristics. The
reason is that the complex phases have too much freedom to model
the observations precisely, as both gijk and hik can have any phase.
There are two possible approaches to overcoming this problem. The
first is to impose some prior or structural constraints to lower the
freedom, for instance, a sparsity constraint on the excitation vkj as
we did in [4], or a clustering constraint on hik . However, this line ofstudy will constitute our future work. The second approach is to ignore the phase gijk of each component tik vkj , and consider only theinter-channel characteristic hik . This leads to another type of multichannel extension of NMF, which is an extension from standard
NMF that we detail in the next subsection.

3.2. Extension from standard NMF
3.2.1. Preprocessing and model
When we apply standard NMF, we only consider the amplitudes of
the time-frequency domain observations. This can be performed on
the 1st channel (reference channel) by„ (1) «∗
xij
(1)
(1)
, ∀ i, j ,(19)
xij ← xij ·(1)

Table 1. Experimental conditions
Sampling rate
16 kHzSTFT frame size and shift 1024 (64 ms) and 256 (16 ms)
Reverberation time
RT60 = 130 ms
6 s of simple instrumental music
Source signals
Direction of sources
70◦ and 150◦ (2 sources)
4 cm
Microphone distance|xij ||xij |

Please note that with the above preprocessing X(1) becomes a nonnegative matrix but X(2) remains a complex matrix.
Thus, the 1st channel is modeled exactly by (1):
Q
(1) P
(21)
p(X(1) |T, V) = i,j N (xij | k tik vkj , 12 )and the 2nd channel is modeled by using an I × K matrix H, whose
elements are complex hik ∈ C, and the same T, V matrices as the
1st channel:
Q(2) Pp(X(2) |T, V, H) = i,j Nc (xij | k hik tik vkj , 1) .
(22)
The difference between the model shown in Sect. 3.1.1 and the above
model lies only in the existence of the complex coefficients gijk .
3.2.2. Negative log-likelihood and update rules
Along with the derivations similar to those described in Sect. 3.1.2
and 3.1.3, we have the following results.
The negative log-likelihood to be minimized is given by
Lm (T, V, H) =
!
X
X
X
(1)
(2)
2
2
tik vkj | + |xij −hik tik vkj | , (23)|xij −i,jkk

and its auxiliary function can be defined by(1), S(2) ) =L+m (T, V, H, S
2
X X |s(1)ijk − tik vkj |
i,j

k

rijk

+

2
X |s(2)
ijk − hik tik vkj |
k

rijk

!
(24)

where S(1) and S(2) are tensors of size I × J × K whose elements(1)(2)
are real sijk and complex sijk ∈ C and satisfying (11).
The updates with respect to S(1) and S(2) are given by
P
(1)
(1)
sijk = tik vkj + rijk (xij − k tik vkj )
(25)
P
(2)
(2)
(26)
sijk = hik tik vkj + rijk (xij − k hik tik vkj )and those with respect to T, V, H are given by
P
(1)
(2) ∗
j x̂ij Re(sijk + sijk hik )P
tik =(1 + |hik |2 ) j x̂ij vkj
P
(1)
(2) ∗
i x̂ij Re(sijk + sijk hik )vkj = P2i x̂ij tik (1 + |hik | )P(2)
j x̂ij sijkPhik =tik j x̂ij vkj

(27)

(28)
(29)

231

3

10
Negative log−likelihood

where ∗ denotes the complex conjugate. To preserve the interchannel characteristics, we apply the same operation to the 2nd„ (1) «∗
channel as
xij
(2)
(2)
, ∀ i, j .
(20)
xij ← xij ·(1)

Random initialization
NMF initialized (5 ite.)
NMF initialized (100 ite.)

2

10

0

20

40
60
# iterations

80

100

Fig. 1. Convergence behavior
P
where x̂ij = k tik vkj , and Re(·) extracts the real part.
In summary, the negative log-likelihood Lm is iteratively minimized by repeating (25)-(26) and one of the three updates (27)-(29).
The set of these updates (25)-(29) can be considered a special caseof the set of updates (12)-(18) shown in Sect. 3.1.3 where gijk isalways fixed at 1 for all i, j, k.
4. EXPERIMENTAL RESULTS
We performed experiments to observe the behavior of the proposed
multichannel NMF. As mentioned in Sect. 3.1.4, the multichannel
extension from complex NMF did not work well so far. Therefore,
we report here the experimental results obtained by the extension
from standard NMF shown in Sect. 3.2.
The experimental conditions are summarized in Table 1. We
considered a situation where two sounds of instruments comingfrom different directions were mixed and observed at two microphones. The mixtures were made by convolving measured impulse
responses, whose reverberation time was 130 ms, and the soundsources. With the application of STFT, we had two matrices X(1)and X(2) both 513×404 in size. Then, we applied the preprocessing(19) and (20), followed by the multichannel NMF updates (25)-(29).
The number K of rank-1 basis matrices was set at 10.
Figure 1 shows the convergence behavior in minimizing the negative log-likelihood (23). We have examined three ways to initialize
T, V matrices. The first was to initialize them with random matrices. The other two involved initializing them with the result ofstandard NMF applied to X(1) . The second and third were different
as regards the number of standard NMF iterations. The H matrix
was initialized randomly in all cases. We see that the proposed updates (25)-(29) successfully minimized the negative log-likelihood
regardless of the initial matrices.
Figure 2 shows the spectrogram reconstructed by the learned T,
V matrices, and Fig. 3 details the T matrix. Although these two matrices were learned appropriately, there was no explicit information
on which basis corresponds to which source. The learned H matrix provided useful information for solving the basis-source binding
problem. Figure 4 shows the phase hik /|hik | of each element. We

500

400

400

Frequency bin index

Frequency bin index

500

300
200
100
100

200
300
Time frame index

400

ti7

20

ti9

0
0
40
0
0
40
20
0
0

400

200

400

200

400

ti6

ti4

200

200

0
0
40
20
0
0
40
20

400

200
400
Frequency bin index

0
0
40
20
0
0
40
20

ti8

0
0
40
20
0
0
40
20

ti2

20

ti10

ti1
ti3
ti5

40

20

0
0

200

400

200

400

200

400

200

400

2
1

300
0
200

−1

100

−2

2

Fig. 2. Spectrogram reconstructed by learned T, V matrices
40

3

4
6
Basis index (k)

8

10

−3

Fig. 4. Phase of each element of H matrix
X 2 2
∂L+
(1) ∗
mc
=
[tik vkj − Re(sijk gijk
)tik
∂vkj
r
ijk
i
(2)

∗
+|hik |2 t2ik vkj − Re(sijk gijkh∗ik )tik ] (31)

X tik vkj
∂L+
(2) ∗
mc
=
(hik tik vkj − sijk gijk
)
∂h∗ik
rijk
j

(32)

∂L+
tik vkj (1)(2)
mc=(sijk − gijk tik vkj + sijk h∗ik − gijk |hik |2 tik vkj )
∗
∂gijk
rijk
(33)
Setting these derivatives to zero and substituting rijk by (14) yields
the updates (15)-(17). For the G update, we use the unit amplitude
constraint |gijk | = 1 to derive (18).

200
400Frequency bin index

7. REFERENCES

Fig. 3. Each column of T matrix

can see that the 1st and 2nd column bases have a tendency to negative phases whereas the others tend towards positive phases. Byclustering those column bases, it was automatically decided that the
1st and 2nd bases belonged to the 1st source and the others belonged
to the 2nd source. As a result of the source separation, the signal-todistortion ratios [9] were improved from (−12.04, 12.04) dB at the
mixture to (4.22, 16.21) dB at the separated sounds.
5. CONCLUSION
We proposed multichannel extensions of complex NMF. We started
with a general formulation where the phase gijk of each component
tik vkj is considered. We then specialized the formulation so thatgijk was always fixed at 1. Experimental results for the specialized version of multichannel NMF were provided to demonstrate theeffectiveness. Future work will include the design of a more sophisticated model advanced from the general formulation.
6. APPENDIX
The partial derivatives of L+
mc (10) with respect to T, V, H, G areshown below.
X 2
∂L+
(1) ∗
2
mc
=
[tik vkj− Re(sijk gijk
)vkj
∂tik
r
ijk
j
(2)

2
∗
+|hik |2 tik vkj− Re(sijk gijkh∗ik )vkj ] (30)232

[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] W. Xu, X. Liu, and Y. Gong, “Document clustering based on
non-negative matrix factorization,” in Proc. ACM SIGIR, 2003,
pp. 267–273.
[3] 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.
[4] H. Kameoka, N. Ono, K. Kashino, and S. Sagayama, “Complex NMF: A new sparse representation for acoustic signals,” inProc. ICASSP 2009, Apr. 2009, pp. 3437–3440.
[5] 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.
[6] R. M. Parry and I. A. Essa, “Estimating the spatial position ofspectral components in audio,” in Proc. ICA 2006. Mar. 2006,
pp. 666–673, Springer.
[7] 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.
[8] J. de Leeuw, “Block-relaxation methods in statistics,” in Information Systems and Data Analysis, H. H. Bock, W. Lenski, andM. M. Richter, Eds., pp. 308–324. Springer Verlag, 1994.
[9] E. Vincent, H. Sawada, P. Bofill, S. Makino, and J. P. Rosca,“First stereo audio source separation evaluation campaign: data,
algorithms and results,” in Proc. ICA 2007, 2007, pp. 552–559.