IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 21, NO. 5, MAY 2013

971

Multichannel Extensions of Non-Negative Matrix
Factorization With Complex-Valued Data
Hiroshi Sawada, Senior Member, IEEE, Hirokazu Kameoka, Member, IEEE, Shoko Araki, Senior Member, IEEE,
and Naonori Ueda, Member, IEEE

Abstract—This paper presents new formulations and algorithms
for multichannel extensions of non-negative matrix factorization(NMF). The formulations employ Hermitian positive semidefinite
matrices to represent a multichannel version of non-negative
elements. Multichannel Euclidean distance and multichannel
Itakura-Saito (IS) divergence are defined based on appropriate
statistical models utilizing multivariate complex Gaussian distributions. To minimize this distance/divergence, efficient optimization algorithms in the form of multiplicative updates are derived
by using properly designed auxiliary functions. Two methods areproposed for clustering NMF bases according to the estimatedspatial property. Convolutive blind source separation (BSS) is
performed by the multichannel extensions of NMF with the clustering mechanism. Experimental results show that 1) the derived
multiplicative update rules exhibited good convergence behavior,and 2) BSS tasks for several music sources with two microphones
and three instrumental parts were evaluated successfully.
Index Terms—Blind source separation, clustering, convolutive
mixture, multichannel, non-negative matrix factorization.

I. INTRODUCTION
ON-NEGATIVE matrix factorization (NMF) is an unsupervised learning technique with a wide range of applications such as parts-based image representation [3], document
clustering [4], and music transcription [5]. As the top part of
Fig. 1 shows, NMF decomposes a given non-negative matrix
into two smaller non-negative matrices and . When we analyze an audio/music signal with NMF, we typically employ ashort-time Fourier transform (STFT) to obtain complex-valued
representations in the time-frequency domain. Then, we makethem non-negative by calculating the (squared) absolute values(see (1)) in order to apply NMF. The bottom part of Fig. 1 showsthat NMF extracts frequent sound patterns as five NMF bases
from an audio clip containing five different notes.
A typical issue with NMF-based audio signal analysis is how
to cluster the extracted NMF bases for a higher-level interpretation of the audio signal. Various NMF models have been proposed for that purpose. Temporal continuity [6] is consideredN

Manuscript received July 26, 2012; revised November 07, 2012; accepted
December 23, 2012. Date of publication January 14, 2013; date of current version February 13, 2013. Earlier versions of this work were presented at WASPAA
2011 [1] and ICASSP 2012 [2] as workshop/conference papers. The associate
editor coordinating the review of this manuscript and approving it for publication was Dr. Roberto Togneri.
The authors are with the NTT Communication Science Laboratories,
NTT Corporation, Kyoto 619-0237, Japan (e-mail: sawada.hiroshi@lab.
ntt.co.jp;
kameoka.hirokazu@lab.ntt.co.jp;
araki.shoko@lab.ntt.co.jp;
ueda.naonori@lab.ntt.co.jp).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TASL.2013.2239990

Fig. 1. Formulation of NMF (top) and its application to a music signal
along with their(bottom). Frequent sound patterns are identified in matrix
activation periods and strengths shown in matrix .

on activation gains . The sequential modeling of bases
is conducted with convolutive NMF [7] and hidden Markov
models [8], [9]. Shifted NMF [10], [11] identifies the bases that
correspond to the notes played by an identical instrument. All
these methods are related to single-channel (monaural) source
separation.
As humans/animals have two ears, multichannel processing
is a way of realizing a more general source separation capability because the spatial properties (directions or locations) of
source signals can be exploited [12]. Specifically with an NMFbased method, the bases can be clustered according the spatial
property, as Fig. 2 shows. With that in mind, multichannel extensions of NMF have been studied with the aim of realizing
sound source separation and localization. In an instantaneous
mixture case, the authors of [13], [14] were interested in thegain of each source to each microphone. Thus, all the values
are still non-negative, and the notion of non-negativity remains
clear. However, in a convolutive mixture case, the phase difference between different microphones is crucial information
for source localization and separation. Thus, we need to handle
complex-valued multichannel observations for these purposes,
but the notion of non-negativity is not obvious. In [15], [16], thepower spectrums of source signals are modeled with non-negative values, but there is no explicit description regarding nonnegativity for mixing matrices or covariance matrices. In thispaper, we propose that the Hermitian positive semidefinitenessof a matrix is a multichannel counterpart of non-negativity, andextend NMF to a multichannel case in a more generic way.
There are several choices available for the distance/divergence measures used in the NMF cost function including the
Euclidean distance [17], the generalized Kullback-Leibler (KL)divergence [17], and the Itakura-Saito (IS) divergence [18]. In
this paper, we define multichannel extensions of the Euclidean

1558-7916/$31.00 © 2013 IEEE

972

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 21, NO. 5, MAY 2013

we changed the update rules for the multichannel Euclidean
NMF from those shown in [1] to make the link to the standard
single-channel Euclidean NMF clearer.
This paper is organized as follows. Section II explains thebasics of the existing single-channel NMF. The proposed multichannel extensions and the clustering techniques are describedin Sections III and IV, respectively. Section V reports experimental results on the convergence behavior of the algorithms
and source separation performance. Section VI concludes thispaper.
Fig. 2. Multichannel extensions of NMF associate the spatial property witheach NMF basis. This enables us to cluster NMF bases according to the source
location, and thus perform a source separation task.

distance and the IS divergence, and extend the NMF modelwith these two definitions. We show that minimizing these distance/divergence is equivalent to maximizing the log-likelihoodof the observations with appropriate statistical models utilizing
multivariate complex Gaussian distributions.
The wide popularity of standard single-channel NMF comesfrom the fact that the algorithm is very easy to implement andworks efficiently. In particular, multiplicative update rules [17]provide rapid convergence and have the attractive property ofguaranteeing the non-negativity of all the matrix elements once
they are initialized with non-negative values. In previous workon multichannel NMF [15], [16], however, expectation-maximization (EM) algorithms have been derived. It was reported
that the algorithms were sensitive to parameter initialization,and thus they used the original source information for perturbed
oracle initializations. This paper presents multiplicative update
rules for multichannel extensions of NMF, and shows experimentally that their convergence behavior is similar to that ofsingle-channel NMFs.
With the multichannel extensions of NMF presented in thispaper, we have the estimations of the spatial properties for eachbasis. To perform a source separation task, we need to cluster
the NMF bases for a source according to the similarity of thespatial properties. This paper proposes an automatic clustering
mechanism that is built into the NMF model with cluster-indicator latent variables. The update rules are slightly changed but
still multiplicative in form.
The main contributions of this paper are summarized asfollows.
1) The notion of non-negativity is defined for a complexvalued vector (Section III-A).
2) Multiplicative update rules are derived for multichannel
extensions of NMF (Section III-C). These updates provide faster convergence than the previous EM algorithms(Fig. 9).
3) Multichannel extensions are found for the Euclidean distance and IS divergence (Section III).
4) Methods for clustering NMF bases are proposed for asource separation task (Section IV).
In our previous work [1], [2], we succeeded merely in separating
two sources and found it difficult to separate more sources.
This paper newly proposes a bottom-up clustering method anda source separation procedure, in which redundant spatial properties are allowed. Consequently, we have succeeded in separating three sources for a variety of music sources. In addition,

II. NON-NEGATIVE MATRIX FACTORIZATIONThis section reviews the formulation and algorithm of standard single-channel NMF [17]–[19]. Let us assume that we havea single-channel audio observation, to which we apply a shorttime Fourier transform (STFT).
A. FormulationLet
be the STFT coefficient at frequency bin andtime frame . To apply NMF, we need to convert
to a nonnegative value
via preprocessing. Typically, we takethe absolute value or its squared value as(1)
represents a complex conjugate. Then, a matrix ,, is constructed with all the preprocessed values
for
and
.
NMF factorizes thematrixinto the product of an
matrixand amatrix . The parameter
specifies the number of NMF bases, and is generally determined
empirically by the user. All the elements of the two matrices,, should be non-negative, i.e.,and.
NMF algorithms are designed to minimize the distance/divergence between the given matrix
and its factored form
.
Let

where

(2)
be the factorization approximation of. Then the distance/
divergence can be defined in a general form
(3)
specifies an element-wise distance/divergence. The
where
following three types of distance/divergence are widely used:
Squared Euclidean distance(4)
Generalized Kullback-Leibler (KL) divergence(5)

SAWADA et al.: MULTICHANNEL EXTENSIONS OF NMF WITH COMPLEX-VALUED DATA

973

C. Probability Distributions Related to Distance/Divergence
There are relations between the three distance/divergences(4)–(6) and specific probability distributions [18], [20], namelya Gaussian distribution , a complex Gaussian distribution
and a Poisson distribution
as shown below. Studying these
relationships helps us to consider multichannel extensions of
NMF in the Section III.
Minimizing the distance/divergenceisequivalent to maximizing the log-likelihoodor, where ,, is a matrix of STFT
coefficients.
Squared Euclidean distanceFig. 3. Three types of distance/divergence withtance, generalized KL divergence and IS divergence.

: squared Euclid dis-

Itakura-Saito (IS) divergence

(10)
(6)

Fig. 3 shows examples of these distance/divergences with. We observe that KL and IS divergences are less sensitive toover-approximation than under-approximation. And from (6),
we observe that IS divergence depends only on the ratio.
Thus, for instance,. This property
is favorable when analyzing most audio signals such as music
and speech, where low frequency components have much higher
energy than high frequency components. This is because low
and high frequency components are treated equally with similar
importance according to the property.

KL divergence

(11)
where
is the Gamma function.
IS divergence

B. Algorithm: Multiplicative Update Rules(12)

We can minimize the distance/divergence according to (3)together with (4), (5), or (6) in the following manner. First, theelements of and are randomly initialized with non-negativevalues. Then, the following update rules [17], [19] are iteratively
applied until convergence.
Squared Euclidean distance(7)

KL divergence

Regarding IS divergence, the likelihood
is calculated not for the matrix of preprocessed non-negative values
but for the matrix of complex-valued STFT coefficients, and
it is thus necessary to specify
in a preprocessing
step for the connection to (6).
When
, the distance/divergence (4), (5) or (6) becomes 0 and each -term of the log-likelihood defined above ismaximized. Therefore, the distance/divergence can be derived
as the difference between the log-likelihoods ofand. We
show the IS divergence case as an example:

(8)
(13)
IS divergence
III. MULTICHANNEL EXTENSIONS OF NMF(9)These update rules are called multiplicative, since each element
is updated by multiplying a scalar value, which is guaranteed to
be non-negative.

This section presents our multichannel extensions of NMF.
Fig. 4 shows an overview of the multichannel extensions (inred), in contrast with standard single-channel NMF (in blue).
We begin with the multichannel extension of IS divergence(6), since this extension is the most natural. We then extend
Euclidean distance (4) to a multichannel case. Unfortunately,

974IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 21, NO. 5, MAY 2013

Fig. 5. Illustrative example of multichannel NMF:,,,
. Non-negative values are shown in gray and complex values are shownin red.

Fig. 4. Variations of NMF presented in this paper. Items that correspond to
standard single-channel NMF are shown in blue. Three distance/divergences
are discussed, and their corresponding probability distributions are presented.
Items that correspond to multichannel extensions of NMF are shown in red.
Two distance/divergence are extended to multichannel.

we have not found a multichannel counterpart for generalized
KL divergence (5).
A. Formulation (IS Divergence)be the number of microphones, andbe a complex-valued vector for atime-frequency slot, withbeing the STFT coefficient at the-th microphone. Letbe such a vector at frequency bin
and time frame . Now, let us introduce a multivariate complex
Gaussian distributionthat extends (12)
Let

(14)
where
is an
covariance matrix that should be Hermitian positive definite. Let
or
..
.

..

.

..
.

(15)

where
and
are non-negative scalars as in the singlechannel case. To solve the scaling ambiguity between
and
, let
have a unit trace.
In a matrix-wise notation, let and beandhierarchical matrices whose elements arematrices, i.e.,and. Fig. 5 provides an illustrative
example in which multichannel NMF factorizes a hierarchicallystructured matrix into the product ofand , where
represents the Hadamard product, i.e.,. The
multichannel NMF is formulated to minimize the total multichannel divergence similar to (3)(18)
is the element-wise multichannel divergence such as
where(16) in the IS divergence case.
B. Formulation (Squared Euclidean Distance)
In this subsection, we consider a multichannel extension of
Euclidean NMF. Thanks to the versatility of a univariate complex Gaussian distribution, we can model the preprocessed
observations ,, as

be the outer product of a complex-valued vector. We then define
the multichannel IS divergence similarly to (13)

(19)
(16)
where
is the trace of a square matrix .
We assume that the source locations are time-invariant in asource separation task (see Fig. 2). Therefore, we introduce amatrixthat models the spatial property of the -th NMF
basis at frequency bin . The matrix is of size
to bematched with the size of. Also, the matrix
is Hermitian
positive semidefinite to possess the non-negativity in a multiwith a sum-of-product formchannel sense. Then, we model

(17)

where
is the squared Frobenius
norm of matrix . Maximizing the log of the likelihood (19)is equivalent to minimizing the distance (18) with the elementwise multichannel distance(20)
Therefore, multichannel Euclidean NMF has been formulated
as minimizing (18) with (20).
When applying standard single-channel Euclidean NMF, it
is typical that the absolute valuein (1) is employedrather than the squared value to prevent some observations frombeing unnecessarily enhanced. In the same sense, an amplitude

SAWADA et al.: MULTICHANNEL EXTENSIONS OF NMF WITH COMPLEX-VALUED DATA

square-rooted version of the outer product (15) would be usefulin the multichannel Euclidean NMF:

975

For both the IS and Euclidean cases, unit-trace normalization
should follow.
D. Derivation of Algorithm (IS-NMF)

..
.

..

.

..
.
(21)

where

.

This subsection explains the derivation of the multiplicative
update rules (22)–(24) for IS divergence. For a given observation , the total distance (18) together with (16) can be writtenas(28)

C. Algorithm: (Multiplicative Update Rules)
As shown in the next two subsections, the following multiplicative update rules are derived to minimize the total distance/divergence (18) with (16) or (20). These update rules reduce to their single channel counterparts (9) and (7) if,
and. Therefore, sets of these updates constitute multichannel extensions of NMF.
IS-NMF (IS divergence)

where constant terms are omitted. To minimize this function, we follow the optimization scheme of majorization [21], [22], in which an auxiliary (majorization) function isused. Let us define an auxiliary function

(22)
(29)

(23)

To updateAppendix I), we solve an algebraic Riccati equation (see
(24)

with

where
is the target matrix before the update.
EU-NMF (Squared Euclidean distance)
(25)

and
where
definiteness,size , andauxiliary function
1)
2)and the equality

are auxiliary variables that satisfy positive
with being an identity matrix of(Hermitian). It can be verified that thehas two properties:
is satisfied when(30)

(see Appendix II for the proof).
The function is indirectly minimized by repeating the following two steps:
1) Minimizing
with respect toandby (30), which
makes.
2) Minimizing
with respect to ,or , which also
minimizes .
For the second step, we calculate the partial derivatives ofw.r.t. the variables
and . Setting these derivatives atzero, we have the following equations.

(26)
(31)
(27)
Post-processing is needed to make
Hermitian and positive semidefinite. This can be accomplished byand then by performing eigenvalue decomposition as, setting all the negative elementsofat zero, and updatingwith the new .
We confirmed empirically that the update (27) followed by the
post-processing always decreases the squared Euclidean distance. However, we have not yet found a theoretical guarantee.

(32)
(33)
By substituting (30) into these equations, we obtain the multiplicative update rules (22)–(24).

976IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 21, NO. 5, MAY 2013

E. Derivation of Algorithm (EU-NMF)
The EU-NMF updates (25)–(27) can be derived in a similarmanner. For a given observation , the total distance (18) together with (20) and (17) can be written as(34)
where constant terms are omitted. To minimize this function, we again follow the optimization scheme of majorization [21], [22]. Let us define an auxiliary function

Fig. 6. An example of learned spatial properties. They are represented as, the phase difference between the first and second microphones,for each frequency bin and NMF basis .

To understand and interpret, we detail it by using therank-1 convolutive model [23] as(35)
that satisfy Hermitian positive defwith auxiliary variables
initeness and. It can be verified that the auxiliary
functionhas two properties:
1)2)
and the equality
is satisfied when(36)
(see Appendix III for the proof).
The function is indirectly minimized by repeating the following two steps:
1) Minimizing
with respect toby (36), which makes.
2) Minimizing
with respect to ,or , which also
minimizes .
For the second step, we calculate the partial derivatives ofw.r.t. the variables
and . Setting these derivatives atzero, we have the following equations.
(37)
(38)

(39)
By substituting (36) into these equations, we obtain the multiplicative update rules (25)–(27).
F. Interpretation of Learned Matrices
The multichannel NMF algorithm, (22)–(24) or (25)–(27),learns matrices , and . The interpretation of the matrices
and is the same as with standard single-channel NMF. This
subsection provides an interpretation of the matrices , which
are particular to the multichannel NMF.

where
is a mixing vector whose-th element
is the Fourier transform of a windowed impulse response from the source to the -th microphone,
is asmall positive scalar, and is an
identity matrix. Then,
we see that the diagonal elements ofrepresent the power
gain of the -th basis at the -th frequency bin to each microphone. And the off-diagonal elements represent the phase differences between microphones.
With a small microphone array, phase differences among
microphones are typically more visible than gain differences.
Fig. 6 shows an example of phase differences that appearedon learned matrices . We interpret that the bases of number
have similar spatial properties, and thus these
are coming from the same direction. There is another groupof NMF bases regardingthat constitute anothersource coming from a different direction. How to cluster NMF
bases according to such spatial properties will be discussed in
Section IV.
IV. CLUSTERING NMF BASES
This section presents two techniques for clustering NMF
bases for a source separation task. The first is a top-down
approach whose clustering mechanism is built in the NMF
model. The second is a bottom-up approach that performs
sequential pair-wise merges. Later, in Section V-C, we describe
a robust source separation procedure that combines these twotechniques.
A. Top-Down Approach With Modified NMF Modelinto
Let us consider clustering matrices
classes, and letindicate whether the -th
matrix belongs to the -th cluster
or not.
in the NMF model can be replaced with,
Then,
and the sum-of-product form (17) is changed to(40)

SAWADA et al.: MULTICHANNEL EXTENSIONS OF NMF WITH COMPLEX-VALUED DATA

977

with
(46)

(47)
is the target matrix before the update.
where
EU-NMF (Squared Euclidean distance)(48)

Fig. 7. The matrices shown in Fig. 6 are clustered into two sets of matrices.
. Latent variables show that -th basis
They are represented as
belongs to the first or second cluster in a soft sense.

The multichannel Euclidean distance (20) and IS divergence(16) can still be employed in the NMF model.
Now, we want to optimize the cluster-indicator latent vari, in the same manner as and . For thatables ,
purpose, let us allow
to have a continuous value such thatand
. We consider that this relaxation
corresponds to estimating the expectation ofaccording to
the posterior probability
instead of the value
itself. But for simplicity, we do not change the notation ofin (40) even with this relaxation. Fig. 7 shows an example,
where ten sets of matrices shown in Fig. 6 are clustered into two
sets of matrices.
The algorithms to minimize(41)
with the element-wise distance/divergence (20) or (16) and withthe model (40) can be derived in a similar manner as explainedin Section III-E. The update rules are in the following forms.
IS-NMF (IS divergence)

(49)

(50)

(51)
For both the IS and Euclidean cases, unit-trace normalization
and unit-sum normalization
should follow.
The clustering result obtained with the top-down approachheavily depends on the initial values of the cluster-indicator latent variables . To prevent an important cluster from disappearing by chance, it is a good idea to have some redundant
clusters by setting the cluster number at a larger than expected.
Then later, the redundant clusters can be merged by employing
the bottom-up clustering shown in Section IV-B.
B. Bottom-Up Clustering by Sequential Merge Operation

(42)

(43)

This subsection explains another way to cluster NMF bases.
It is based on a pair-wise merge operation, in which a pair withthe minimum distance is identified and merged. The pair-wise
distance between the -th set and the -th set is defined byusing the Frobenius norm as
(52)

(44)
For

, we solve an algebraic Riccati equation(45)

Algorithm 1 shows the whole procedure. Inside the basic updates for NMF, the pair-wise merge operation is interleaved at a
. In the pair-wise merge oprate specified by a variablewith the minimum distance is found, and aeration, a pair
new setis calculated as element-wise weighted
means of the two sets. Then, the number of clusters is decreased

978IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 21, NO. 5, MAY 2013by 1 as the new set is added and the -th and -th sets are removed. The main loop is repeated until the number of clusters.
becomes a specified number
Algorithm 1 Multichannel NMF with bottom-up clustering
1: Procedure MchNMF_BottomUpClustering2:
3: Whiledo4:
5:
update by (42) or (48)6:
update by (43) or (49)
7:
If
then
8:
9:
10:
end if11:
update by (45) or (51)12:
update by (44) or (50)13:end while
14: end procedure

Fig. 8. Experimental setup for room impulse responses.
TABLE IMUSIC SOURCES

15: Procedure
16:
17:
18:
19:
20:
21:
22: end procedure
C. Source Separation
If NMF bases are appropriately clustered, source separation
can be performed by using Wiener filters. Remember thatis an -dimensional complex vector representing STFT coefficients at frequency bin and time frame , and letbethe STFT coefficient vector for the -th separated signal. Then,the separated signals are obtained by the single-channel Wiener
filter for the -th channel

parts. Sets of stereo mixtures were generated by convolving themusic parts and the impulse responses measured in a real room
whose conditions are shown in Fig. 8. The impulse responses
were measured by using a maximum length sequence generated
by a 17-th order polynomial over GF (2). We made four sets of
mixtures using the music sources listed in Table I, which canbe found at the professionally produced music recordings pageof the Signal Separation Evaluation Campaign (SiSEC 2011)[24]. The mixtures were down-sampled to 16 kHz. The STFT
frame size was 64 ms and the frame shift was 16 ms. The algorithms were coded with Matlab and run on an Intel Xeon W3690(3.46 GHz) processor.
B. Convergence Behavior(53)

or by the multichannel Wiener filter
(54)
where

is the sum-of-product form defined in (40).
V. EXPERIMENTS

A. Experimental Setups
We examined the proposed multichannel extensions of NMF
with stereo
music mixtures that contained three music

Let us first show the convergence behavior of the multichannel NMF algorithms proposed in Section III. For comparison, we run the algorithms of single-channel NMF (7) and(9), multichannel NMF (25)–(27) and (22)–(24), and the EM
algorithm for multichannel IS-NMF shown in Appendix IV.
We set the number of NMF bases, and used a music
mixture with.
Fig. 9 shows the convergence behavior for 1000 iterations,and Table II shows the computational time. NMF algorithms
with IS divergence generally take more time than EU-NMF.
We observe that the convergence behavior of the single-channel
NMF algorithms is similar to that of the proposed multichannel
NMF algorithms. With respect to multichannel IS-NMF, theproposed algorithms show faster convergence and computation
than the EM algorithm.

SAWADA et al.: MULTICHANNEL EXTENSIONS OF NMF WITH COMPLEX-VALUED DATA

Fig. 9. Convergence behavior shown in log-log plots: EU-NMF (top) and
IS-NMF (bottom), single-channel (left) and multichannel (right).

TABLE IICOMPUTATIONAL TIME (IN SECONDS) FOR 1000
ITERATIONS WITH 14-SECOND SIGNALS

C. Source Separation Procedure
This subsection explains the source separation procedure
with the multichannel NMF proposed in Section IV. We
adopted the following procedure so that the spatial properties
of sources were well extracted.
1) Preprocessing: For EU-NMF: we normalized the powerof the observation vectors at each frequency bin such that
. This was to prevent the low frequency components from being dominant in the total distance (18). Then, we
generated matrices
as the amplitude square-rooted outer
products (21). For IS-NMF: we generated matrices
by theouter products (15). To prevent the determinant of multichannel
IS divergence (16) from being zero, we then added a regularization term to each matrix aswithand being an identity matrix.
2) Initialization: Matrices andwere randomly initialized with non-negative entries. The diagonal elements ofwere initially all set at, and the off-diagonal elements were
initially all set at zero. The elements of matrix were initializedwith random values around.
3) Multichannel NMF: We set the number of NMF bases
at 30 (ten times the number of sources). As for clustering thebases for each source, we employ both top-down and bottom-up
approaches as follows.
i) 20 iterations to update and .
ii) 200 iterations to update , , and by the top-down
approach with.
iii) Bottom-up clustering with
until.
iv) 200 iterations to update , , and by the top-down
approach with.

979

Fig. 10. Convergence behavior of two different clustering strategies. The blue
corresponds to the procedure described in Section V-C. It
lineused both top-down and bottom-up approaches. The zigzag pattern (iterations
220–280) shows that the IS divergence was increased by a pair-wise merge opcorresponds to using only the top-down
eration. The red line
approach.

Having redundantspatial properties (step ii) followedby the bottom-up clustering (step iii) contributes to robust estimations of the spatial properties. Some results will be shown in
the next subsection with Fig. 12.
4) Separation: For EU-NMF, a single-channel Wiener filter(53) was used for each channel. For IS-NMF, a multichannel
Wiener filter (54) was used. We selected these configurations
because each of these produced better results empirically.
Fig. 10 shows an example in which the IS divergence wasminimized by the procedure. The blue line shows how stepsi)–iv) work, especially when parameteris decreased from
9 to 3. The red line shows the case where only the top-down
approach was employed. In this example, the IS divergence
was better minimized by having redundantspatial properties.
D. Source Separation Results
The separation performance was numerically evaluated interms of the signal-to-distortion ratio (SDR) [25]. We need toknow all the source images
for all microphones
and sources. To calculate
for the-th source, we first decompose the time-domain multichannel
signals
as
(55)
,
, and
are unwanted error comwhere
ponents that correspond to spatial (filtering) distortion, interferences, and artifacts, respectively. These are calculated by using
a least-squares projection [25]. Then,is calculated as thepower ratio between the wanted and unwanted components

Fig. 11 shows the source separation results obtained withthe procedure described in the last subsection (EU-NMF and980

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 21, NO. 5, MAY 2013

Fig. 11. Source separation performance evaluated in terms of SDRs averand
aged over all the three sources. With NMF-based methods (), ten trials were conducted for each mixture ID. The error bars
represent one standard deviation. For comparison, the results obtained with anareexisting underdetermined blind source separation method [26]also shown.

TABLE IIICOMPUTATIONAL TIME (IN SECONDS) FOR THE SEPARATION PROCEDUREDESCRIBED IN SECTION V-C AND AN EXISTING METHOD [26]IS-NMF). For comparison, results obtained with the Underdetermined Blind Source Separation (UBSS) method [26] arealso shown. Four sets of mixtures whose sources are listed inTable III were examined. The source separation result obtainedwith the NMF-based method depends on the initial values of, and . Therefore, we conducted ten trials with different
initializations for each set of mixtures. Table III shows thecomputational time. Sound examples can be found at our web
page [27].
From these results, we observe the following. For music
recordings with frequent sound patterns, the new NMF-based
methods generally performed better than the existing method[26] that relies on the spatial property and simple time-wise
activity information of each source. However, the computational burden of the NMF-based methods was heavy. Thiswas because many operations related to matrix inversions
and eigenvalue decompositions were involved in the NMF
updates. Among the NMF-based methods, IS-NMF produced
clearly better separation results than EU-NMF with increased
computational effort. This result supports the superiority of
IS divergence for audio signal modeling [18] also in a multichannel scenario.
Fig. 12 show the effect of having the redundant spatial properties mentioned in the previous subsection. We observe thateven a small number of redundant spatial properties contributed
to a better separation result than not having the redundancy. Having too much redundancy also worked well,
but the computational demands also increased. Therefore, we
employed
for the experiments whose results areshown in Fig. 11.

Fig. 12. Source separation performance with various numbers of redundant. Ten trials were conducted for each mixture ID, andspatial properties
thus the error bars represent one standard deviation over 40 averaged SDRs.

VI. CONCLUSION
We have presented new formulations and algorithms for multichannel NMF. We started with a simple model (17) where each
NMF basis ,, has its own spatial properties.
Then, to cluster the NMF bases according to their spatial properties, we introduced the second model (40), which includes
cluster-indicator variables . Multiplicative update rules were
derived to minimize the multichannel IS divergence (16) or multichannel Euclidean distance (20). Experimental results showthat the derived multichannel algorithms were as efficient as thestandard single-channel algorithms in terms of the number of
iterations to converge. Multichannel NMF with IS divergenceproduced better source separation results than with Euclidean
distance for several stereo music mixtures. Future work will include the automatic determination of such model complexity
parameters as the number of NMF bases and sources (spatial
properties) , for example by employing Bayesian nonparametrics [28], [29]. Computationally efficient implementations of themultichannel NMF algorithms, such as one using a general-purpose graphics processing unit (GPGPU), will also constitute future work.
APPENDIX ISOLVING AN ALGEBRAIC RICCATI EQUATION
To solve (24), we perform an eigenvalue decomposition of amatrix(56)
and letbe eigenvectors with negative eigenvalues.
It is theoretically guaranteed that there are exactly
negative
eigenvalues. But in reality there may be computer arithmetic
errors. Thus, we actually sort the
eigenvectors according to
the corresponding eigenvalues in ascending order and employ
the first
eigenvectors.
Then, let us decompose the-dimensional eigenvectors
as(57)

SAWADA et al.: MULTICHANNEL EXTENSIONS OF NMF WITH COMPLEX-VALUED DATA

with
and
calculated by

being-dimensional vectors. The new

is

981

A solution for this equation is given by
(66)

(58)

and adding this for

gives

with
and
. Again to compensate for computer arithmetic errors, we ensure
is Hermitian
by
.

(67)
with the fact that
auxiliary function

APPENDIX IIPROOF FOR THE AUXILIARY FUNCTION (29)defined in (29) with
Let us consider the minimization ofrespect to and subject to the constraint. By
introducing Lagrange multipliersof size, we have(68)
and the minimum value is equal towith respect toat zero(60). Adding this for
giveswith the fact that.
Therefore, the minimum of the auxiliary function
is obtained
when
we have(61)
and the minimum value is equal to defined in (28).
The partial derivative ofwith respect to Hermitian matrix
is given by [30]defined in (34).

APPENDIX IV
EM ALGORITHM FOR COMPARISON(59)
By setting the partial derivative of. Therefore, the minimum of theis obtained whenThis appendix shows an EM algorithm designed to minimize
the total multichannel IS divergence, (18) with (16), according
to the NMF model (17). The EM algorithm shown here is a
simplification of the EM algorithm shown in [16]. For the STFT
coefficient vectors, let
be latent vectorsthat satisfy.
E-step: calculate the expectation of the outer product of
by
(69)
with
(70)
M-step: update the NMF model parameters by(62)
Setting this zero gives

(71)

.

(72)

APPENDIX III
PROOF FOR THE AUXILIARY FUNCTION (35)Let us consider the minimization ofrespect tosubject to the constraint
ducing Lagrange multipliersof size(73)

defined in (35) with
. By intro, we have(63)

The partial derivative ofis given by [30]with respect to Hermitian matrix(64)
Setting this zero and introducing a matrix

we have(65)

REFERENCES[1] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “New formulations
and efficient algorithms for multichannel NMF,” in Proc. WASPAA ’11,
Oct. 2011, pp. 153–156.
[2] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Efficient algorithms
for multichannel extensions of Itakura-Saito nonnegative matrix factorization,” in Proc. ICASSP ’12, Mar. 2012, pp. 261–264.
[3] D. D. Lee and H. S. Seung, “Learning the parts of objects with nonnegative matrix factorization,” Nature, vol. 401, pp. 788–791, 1999.
[4] W. Xu, X. Liu, and Y. Gong, “Document clustering based on non-negative matrix factorization,” in Proc. ACM SIGIR, 2003, pp. 267–273.
[5] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization for
polyphonic music transcription,” in Proc. WASPAA ’03, Oct. 2003, pp.
177–180.
[6] T. Virtanen, “Monaural sound source separation by nonnegative matrix
factorization with temporal continuity and sparseness criteria,” IEEE
Trans. Audio, Speech, Lang. Process., vol. 15, no. 3, pp. 1066–1074,
Mar. 2007.

982

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 21, NO. 5, MAY 2013[7] P. Smaragdis, “Convolutive speech bases and their application tosupervised speech separation,” IEEE Trans. Audio, Speech, Lang.
Process., vol. 15, no. 1, pp. 1–12, Jan. 2007.
[8] M. Nakano, J. L. Roux, H. Kameoka, Y. Kitano, N. Ono, and S.
Sagayama, “Nonnegative matrix factorization with Markov-chainedbases for modeling time-varying patterns in music spectrograms,” in
Latent Variable Analysis and Signal Separation. New York, NY,
USA: Springer, 2010, pp. 149–156.
[9] G. Mysore, P. Smaragdis, and B. Raj, “Non-negative hidden Markov
modeling of audio with application to source separation,” in Latent
Variable Analysis and Signal Separation. New York, NY, USA:
Springer, 2010, pp. 140–148.
[10] D. Fitzgerald, M. Cranitch, and E. Coyle, “Shifted non-negative matrix factorisation for sound source separation,” in Proc. IEEE/SP 13thWorkshop Statist. Signal Process., 2005, pp. 1132–1137.
[11] R. Jaiswal, D. FitzGerald, D. Barry, E. Coyle, and S. Rickard, “Clustering NMF basis functions using shifted NMF for monaural soundsource separation,” in Proc. ICASSP ’11, May 2011, pp. 245–248.
[12] J. Blauert, Spatial Hearing: The Psychophysics of Human Sound Localization. Cambridge, MA, USA: MIT Press, 1997.
[13] D. FitzGerald, M. Cranitch, and E. Coyle, “Non-negative tensor factorisation for sound source separation,” in Proc. Irish Signals Syst.
Conf., Sep. 2005, pp. 8–12.
[14] R. M. Parry and I. A. Essa, “Estimating the spatial position of spectral components in audio,” in Proc. ICA ’06, Mar. 2006, pp. 666–673,
Springer.
[15] A. Ozerov and C. Févotte, “Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation,” IEEE Trans.
Audio, Speech, Lang. Process., vol. 18, no. 3, pp. 550–563, Mar. 2010.
[16] 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 ’10, May 2010, pp. 1–4.
[17] D. Lee and H. Seung, “Algorithms for non-negative matrix factorization,” Adv. Neural Inf. Process. Syst., vol. 13, pp. 556–562, 2001.
[18] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music
analysis,” Neural Comput., vol. 21, no. 3, pp. 793–830, 2009.
[19] M. Nakano, H. Kameoka, J. L. Roux, Y. Kitano, N. Ono, and S.
Sagayama, “Convergence-guaranteed multiplicative algorithms for
non-negative matrix factorization with beta-divergence,” in Proc.
MLSP ’10, Aug. 2010, pp. 283–288.
[20] A. Cemgil, “Bayesian inference for nonnegative matrix factorisation
models,” Comput. Intell. Neurosci., vol. 2009, 2009.
[21] J. de Leeuw, “Block-relaxation methods in statistics,” in Information
Systems and Data Analysis, H. H. Bock, W. Lenski, and M. M. Richter,
Eds. Berlin, Germany: Springer-Verlag, 1994, pp. 308–324.
[22] A. Marshall, I. Olkin, and B. Arnold, Inequalities: Theory of majorization and its applications. Berlin, Germany: Springer-Verlag, 2010.
[23] N. Duong, E. Vincent, and R. Gribonval, “Under-determined reverberant audio source separation using a full-rank spatial covariancemodel,” IEEE Trans. Audio, Speech, Lang. Process., vol. 18, no. 7,
pp. 1830–1840, Sep. 2010.
[24] S. Araki, F. Nesta, E. Vincent, Z. Koldovsky, G. Nolte, A. Ziehe,and A. Benichoux, “The 2011 signal separation evaluation campaign(SiSEC2011):-audio source separation,” Latent Variable Analysis andSignal Separation, pp. 414–422, 2012.
[25] E. Vincent, H. Sawada, P. Bofill, S. Makino, and J. Rosca, “First stereo
audio source separation evaluation campaign: Data, algorithms and results,” in Proc. ICA, 2007, pp. 552–559 [Online]. Available: http://
www.irisa.fr/metiss/SASSEC07/[26] H. Sawada, S. Araki, and S. Makino, “Underdetermined convolutive
blind source separation via frequency bin-wise clustering and permutation alignment,” IEEE Trans. Audio, Speech, Lang. Process., vol. 19,
no. 3, pp. 516–527, Mar. 2011.
[27] [Online]. Available: http://www.kecl.ntt.co.jp/icl/signal/sawada/
demo/mchnmf/[28] M. Hoffman, D. Blei, and P. Cook, “Bayesian nonparametric matrix
factorization for recorded music,” in Proc. ICML, 2010, pp. 641–648.
[29] M. Nakano, J. Roux, H. Kameoka, T. Nakamura, N. Ono, and S.
Sagayama, “Bayesian nonparametric spectrogram modeling based on
infinite factorial infinite hidden Markov model,” in Proc. WASPAA’11, Oct. 2011, pp. 325–328.
[30] A. Hjorungnes, Complex-Valued Matrix Derivatives. Cambridge,
U.K.: Cambridge Univ. Press, 2011.

Hiroshi Sawada (M’02–SM’04) received the B.E.,
M.E,. and Ph.D. degrees in information science from
Kyoto University, Kyoto, Japan, in 1991, 1993 and2001, respectively.
He joined NTT Corporation in 1993. He is now
the group leader of Learning and Intelligent Systems
Research Group at the NTT Communication Science
Laboratories, Kyoto, Japan. His research interests include statistical signal processing, audio source separation, array signal processing, machine learning, latent variable model, graph-based data structure, andcomputer architecture.
From 2006 to 2009, he served as an associate editor of the IEEE Transactions
on Audio, Speech & Language Processing. He is a member of the Audio andAcoustic Signal Processing Technical Committee of the IEEE SP Society. He
received the 9th TELECOM System Technology Award for Student from the
Telecommunications Advancement Foundation in 1994, the Best Paper Awardof the IEEE Circuit and System Society in 2000, and the MLSP Data Analysis
Competition Award in 2007. Dr. Sawada is a member of the IEICE and the ASJ.
Hirokazu Kameoka (M’07) received B.E., M.E.
and Ph.D. degrees all from the University of Tokyo,
Japan, in 2002, 2004 and 2007, respectively. He
is currently a research scientist at the NTT Communication Science Laboratories and an Adjunct
Associate Professor at the University of Tokyo. His
research interests include computational auditory
scene analysis, statistical signal processing, speech
and music processing, and machine learning. He isa member of the IEICE, the IPSJ and the ASJ. He
received 13 awards over the past 9 years, including
the IEEE Signal Processing Society 2008 SPS Young Author Best Paper Award.
Shoko Araki (M’01–SM’12) is with NTT Communication Science Laboratories, NTT Corporation,
Japan. She received the B.E. and the M.E. degrees
from the University of Tokyo, Japan, in 1998 and
2000, respectively, and the Ph.D. degree from
Hokkaido University, Japan in 2007.
Since she joined NTT in 2000, she has been engaged in research on acoustic signal processing, arraysignal processing, blind source separation (BSS) applied to speech signals, meeting diarization and auditory scene analysis.
She was a member of the organizing committee of the ICA 2003, the finance chair of IWAENC 2003, the co-chair of a special session on undetermined sparse audio source separation in EUSIPCO 2006, the registration chairof WASPAA 2007, and the evaluation co-chair of SiSEC2008, 2010 and 2011.
She received the 19th Awaya Prize from Acoustical Society of Japan (ASJ) in
2001, the Best Paper Award of the IWAENC in 2003, the TELECOM System
Technology Award from the Telecommunications Advancement Foundation in
2004, the Academic Encouraging Prize from the Institute of Electronics, Information and Communication Engineers (IEICE) in 2006, and the Itakura Prize
Innovative Young Researcher Award from (ASJ) in 2008. She is a member of
the IEICE and the ASJ.
Naonori Ueda (M’02) received B.S., M.S., andPh.D. degrees in communication engineering from
Osaka University, Osaka, Japan, in 1982, 1984,
and 1992, respectively. In 1984, he joined the NTT
Electrical Communication Laboratories, Kanagawa,Japan. In 1991, he joined the NTT Communication
Science Laboratories, Kyoto, Japan, as a seniorresearch scientist. His current research interests
are statistical machine learning and its applications
to pattern recognition, signal processing and data
mining. He is the director of the NTT Communication Science Laboratories. He is a Guest Professor of Nara Advanced Instituteof Science and Technology (NAIST) and National Institute of Informatics(NII). He is also the sub-project leader of Funding Program for World-Leading
Innovative R&D on Science and Technology (First Program), cabinet office,government of Japan, March 2010–February 2014. From 1993 to 1994, he
was a visiting scholar at Purdue University, West Lafayette, Indiana. He is
an associate editor of Neurocomputing, and is a member of the IPSJ. He is afellow of the IEICE.