14th International Society for Music Information Retrieval Conference (ISMIR 2013)

BEYOND NMF: TIME-DOMAIN AUDIO SOURCE SEPARATIONWITHOUT PHASE RECONSTRUCTION
Kazuyoshi Yoshii1 Ryota Tomioka2 Daichi Mochihashi3 Masataka Goto1
1
National Institute of Advanced Industrial Science and Technology (AIST)2
The University of Tokyo 3 The Institute of Statistical Mathematics (ISM){k.yoshii, m.goto}@aist.go.jp tomioka@mist.i.u-tokyo.ac.jp daichi@ism.ac.jp

PSDTF

ABSTRACT

Factorize

This paper presents a new fundamental technique for source
separation of single-channel audio signals. Although nonnegative matrix factorization (NMF) has recently become
very popular for music source separation, it deals only withthe amplitude or power of the spectrogram of a given mixture signal and completely discards the phase. The component spectrograms are typically estimated using a Wiener
ﬁlter that reuses the phase of the mixture spectrogram, butsuch rough phase reconstruction makes it hard to recoverhigh-quality source signals because the estimated spectrograms are inconsistent, i.e., they do not correspond to any
real time-domain signals. To avoid the frequency-domain
phase reconstruction, we use positive semideﬁnite tensor
factorization (PSDTF) for directly estimating source signals from the mixture signal in the time domain. Since PSDTF is a natural extension of NMF, an efﬁcient multiplicative update algorithm for PSDTF can be derived. Experimental results show that PSDTF outperforms conventional
NMF variants in terms of source separation quality.

: Discrete Fourier transform (DFT) matrix: Local signal of frame: Complex spectrum of frame
NMF focuses only on diagonal elements(noticeable correlation structures are ignored)[dB]

Tensor data: a set of covariance matrices(a set of positive semideﬁnite matrices)

Factorize
NMF

Matrix data: a set of power spectra(a set of nonnegative vectors)

Figure 1. PSDTF is a natural extension of NMF.

1. INTRODUCTION
Source separation of music audio signals is a fundamental
task for music information retrieval (MIR). High-quality
source separation could help users ﬁnd their favorite songs
according to the content (such as vocals or instruments) [1].
It would also let them enjoy active music listening [2] basedon the remixing of existing instrumental parts [1–4].
Nonnegative matrix factorization (NMF) [5] has recently
played a key role in the source separation of single-channel
audio signals. It can approximate a nonnegative matrix (theamplitude or power spectrogram of a given mixture signal)as the product of two nonnegative matrices— a set of basis
spectra and a set of the corresponding activations. Then thecomplex spectrogram of the mixture signal is decomposed
into a sum of source spectrograms by using a Wiener ﬁlter
that simply reuses the original phase. However, we cannot
recover high-quality source signals from the decomposed
spectrograms having the unreal phase.
Permission to make digital or hard copies of all or part of this work forpersonal or classroom use is granted without fee provided that copies are
not made or distributed for proﬁt or commercial advantage and that copies
bear this notice and the full citation on the ﬁrst page.
c 2013 International Society for Music Information Retrieval.


369

Considerable effort has been devoted to estimating consistent complex spectrograms that correspond to real timedomain signals. To reconstruct the phase of a given amplitude spectrogram, Grifﬁn and Lim [6] proposed an iterative
short-time Fourier transform (STFT) method that estimatesa time-domain signal such that its amplitude spectrogramis closest to the given spectrogram. Le Roux et al. [7] proposed a cost function that evaluates the inconsistency of
a complex spectrogram and derived an efﬁcient algorithm
for minimizing the cost function [8]. Kameoka et al. [9],on the other hand, formulated complex NMF for directly
factorizing a complex spectrogram. The cost function evaluating the inconsistency could be integrated into complex
NMF as suggested in [10]. Note that improved consistency
does not always result in improved sound quality.
To circumvent the phase reconstruction, we use positive
semideﬁnite tensor factorization (PSDTF) [11] for timedomain separation of single-channel audio signals. Insteadof explicitly considering the wave shapes and phases of basis signals, we focus on the statistical characteristics (e.g.,
periodicity and whiteness) of those signals. More speciﬁcally, we assume that each basis signal follows a Gaussian
process (GP) having a stationary kernel. A given mixture
signal consisting of multiple basis signals is thus locallymodeled by a GP with a convex combination of the corresponding kernels. These kernels can be estimated froma set of local covariances of the mixture signal by using a
multiplicative update algorithm.
We can show that the time-domain PSDTF has an equivalent frequency-domain representation used for factorizing
a mixture spectrogram like NMF. As shown in Figure 1,
PSDTF deals with a set of Hermitian positive semideﬁnite
matrices (outer products of complex spectra) for considering the correlations between frequency components. Thisis reasonable because the discrete Fourier transform (DFT)14th International Society for Music Information Retrieval Conference (ISMIR 2013)

plex spectrogram of the mixture signal and S k ∈ CM ×N
be that of the k-th source signal. If the mixture signal is an
instantaneous mixture of K source signals, we can say

cannot perfectly decorrelate the frequency components ofthe mixture signal. On the other hand, NMF focus only ona set of the nonnegative diagonal elements of those matrices (power spectra) by discarding the correlations between
frequency components. This indicates that PSDTF is a natural and elegant extension of NMF.

S=

k
E[Smn
|Smn ] =The goal of NMF is to approximate a nonnegative matrix
X ∈ RM ×N as the product of two nonnegative matrices
W ∈ RM ×K and H ∈ RK×N as follows:(1)

2.3 Source Separation based on KL-NMF

where W and H respectively represent a set of basis vectors and a set of the corresponding weight vectors, K min(M, N ) is the number of basis vectors, and Y ∈ RM ×N
represents a reconstruction matrix. Eq. (1) can be rewritten
in an element-wise manner as follows:
K


def

Wmk Hkn = Ymn .

kYmn
Wmk Hkn
Smn = Smn . (5)Ymn
k Wmk Hkn

Eq. (5) is known as Wiener ﬁltering in which the original
phase of S is directly attached to each S k . The real-valued
source signal can then be recovered from E[S k |S] by using the overlap-add synthesis method [16]. Note that thecomplex spectrogram of the resulting source signal is unlikely to be equal to E[S k |S] because in general E[S k |S]is an inconsistent spectrogram that does not correspond toany actual time-domain signals.

2.1 Nonnegative Matrix Factorization

Xmn ≈

(4)

Given the mixture spectrogram S (observed variable), eachcomponent spectrogram S k (latent variable) can be estimated in a probabilistic manner as follows:This section aims to reveal the theoretical basis underlying
nonnegative matrix factorization (NMF) in the context ofsource separation. We review two popular variants of NMF
called KL-NMF [12] and IS-NMF [13] and how these variants can be used for source separation.

def

Sk.

k=1

2. FREQUENCY-DOMAIN SOURCE SEPARATION

X ≈ WH = Y ,K

KL-NMF is used for factorizing the amplitude spectrogram[12], i.e., Xmn = |Smn |. The cost function is given bymn
CKL (Xmn |Ymn ) = Xmn log X
Ymn −Xmn +Ymn . Note that
KL-NMF is not theoretically justiﬁed for source separation
because the cost function is scale-variant, i.e., CKL (X|Y ) =
CKL (αX|αY ) for a positive number α.
The probabilistic model of KL-NMF can be formulatedk| is Poissonby assuming that each latent component |Smnkdistributed with a mean parameter Ymn as follows:

k  k
k
|Smn
| Ymn ∼ Poisson(Ymn
).
(6)

(2)

k=1
k
We here let Ymn= Wmk Hkn be a component reconstrucktion such that Ymn = k Ymn. A popular way to evaluate
the reconstruction error C(Xmn |Ymn ) between Xmn and
Ymn is to use the Bregman divergence [14] deﬁned as

We here assume that the condition for amplitude additivity
is satisﬁed, i.e., that the phase of each S k isequal to that ofk
S. Eq. (4) can then be written as|Smn | = k |Smn|. Usk, the reproducing
ing Xmn = |Smn | and Ymn = k Ymn
property of the Poisson distribution gives
Xmn Ymn ∼ Poisson(Ymn ).
(7)

Cφ (Xmn |Ymn )= φ(Xmn ) − φ(Ymn ) − φ (Ymn )(Xmn − Ymn ), (3)where φ is a strictly convex function. This divergence isno less than zero and is zero only when Xmn = Ymn . The
Kullback-Leibler (KL) divergence (φ(x) = x log x − x)and the Itakura-Saito (IS) divergence (φ(x) = − log x) are
well-known special cases of the Bregman divergence. To
estimate
W and H such that the cost function Cφ (X|Y ) =
mn Cφ (Xmn |Ymn ) is minimized, we can use an efﬁcient
multiplicative update (MU) algorithm [15].

This probabilistic model based on superimposed Poiskk|}K
son variables {|Smnk=1 satisfying |Smn | =k |Smn |enables us to calculate the expectation of each latent varik| in a principled manner as follows:able |Smn

k 
k
−1
E[|Smn| |Smn |] = Ymn
Ymn|Smn |.
(8)

2.2 Application to Source Separation

Since the phase is assumed to be preserved, we get Eq. (5).

The goal of source separation is to decompose a given mixture signal into the sum of K source signals. NMF enables
us to perform source separation in the frequency domain.
We regard the nonnegative spectrogram of the mixture signal as an X for which M is the number of frequency bins
and N is the number of frames. We then factorize the given
spectrogram X as X ≈ W H, where W and H respectively represent a set of basis nonnegative spectra and a setof the corresponding temporal activations.
A probabilistic formulation of NMF enables us to estimate latent source signals. Let S ∈ CM ×N be the com-

2.4 Source Separation based on IS-NMF
IS-NMF is used for factorizing the power spectrogram [13],i.e., Xmn = |Smn |2 . The cost function is CIS (Xmn |Ymn ) =Xmn
Xmn
Ymn − log Ymn − 1. IS-NMF is suitable for source separation because the cost function is scale-invariant.
The probabilistic model of IS-NMF can be formulatedkis complex
by assuming that each latent component SmnkGaussian distributed with a variance Ymn as follows:
k
k
k
Smn
|Ymn
∼ Nc (0, Ymn
).

370

(9)

14th International Society for Music Information Retrieval Conference (ISMIR 2013)

3.2 Probabilistic Formulation

Local signals extracted
by using a short window

We formulate a probabilistic model of Eq. (13). A key feature is to focus on the stationary characteristics of the basis
signal. Since the stationarity means that {φkn }N
n=1 are expected to have the same covariance, we put a multivariate
Gaussian prior shared over all frames as follows:Time
Sine wave with time-varying scale

Figure 2. Local signals xkn and xkn have different wave
shapes and phases, but share the same periodicity.

φkn ∼ N (0, V k ), k
 kUsing Smn = k Smn
and Ymn = k Ymn, the reproducing property of the Gaussian distribution gives
Smn |Ymn ∼ Nc (0, Ymn ).

where V k ∈ RM ×M is a full covariance matrix. The meanis set to a zero vector because an audio signal is recorded
as real numbers distributed on both sides of zero.
We can say that the k-th basis signal is Gaussian process(GP) distributed with a stationary (shift-invariant) kernel
V k . Since in reality the signal exists over continuous time,
it is essential to consider a probability distribution of thecontinuous signal. Such a distribution is a GP by deﬁnition
because its marginal over any M discrete time points is a
Gaussian distribution, as indicated by Eq. (14). If V k is aperiodic kernel, {φkn }N
n=1 are expected to have a certain
period but their phases and wave shapes can be different.
We will derive a likelihood of the observed signal xn .
The linear relationship xkn = πkn φkn and Eq. (14) leadto a likelihood of xkn as follows:(10)

Using Xmn = |Smn |2 , we get an exponential distribution:
Xmn |Ymn ∼ Exponential(Ymn ).

(11)

The probabilistic model based on superimposed kGausk}K
satisfying
S=sian variables {Smnmn
k=1k Smn enables us to represent a posterior distribution of latent varik
able Smn
as a Gaussian distribution whose mean and variance are given by Eq. (5) and
k
k
k
−1 k
V[Smn
|Smn ] = Ymn
− Ymn
Ymn
Ymn .

(12)

3. TIME-DOMAIN SOURCE SEPARATION2
xkn |π, V ∼ N (0, πkn
V k ).

This section recasts the problem of source separation in thetime domain. We propose a probabilistic model of sourcesignal superimposition and show how latent source signals
can be estimated in a probabilistic manner.

The goal of source separation is to decompose a given mixture signal into the sum of K source signals. This decomposition is performed on a frame-by-frame basis. Suppose
we have a set of N samples O = [x1 , · · · , xN ] ∈ RM ×N ,where xn ∈ RM is a local signal in the n-th frame extracted from the mixture signal by using a window of size
M . Each xn can be decomposed as follows:
K


xkn ,

(15)

Then, using the reproducing property of the Gaussian distribution and the linear relationship given by Eq. (13), we
get the likelihood of xn as follows:

 K

2
xn |π, V ∼ N 0,
(16)
πkn
Vk .

3.1 Problem Speciﬁcation

xn =

(14)

k=1

Note that Eq. (16) does not include φkn , i.e., all possibilities of φkn are taken into account. This formulation frees
us from explicitly considering the phase of φkn . We here2
≥ 0, X n = xn xTn deﬁne some symbols
as Hkn = πkn
0, and Y n = k Hkn V k  0, where Ψ  0 means that
Ψ is a positive semideﬁnite (PSD) matrix. Then Eq. (16)gives the log-likelihood of X n as follows:

(13)

1
1
c
log p(X n |Y n ) = − log |Y n | − tr(X n Y −1n ), (17)22k=1where xkn is a local signal extracted from the k-th source
signal. This time-domain formulation is equivalent to thefrequency-domain formulation given by Eq. (4). AlthoughN{xkn }Nn=1 are different, we assume that {xkn }n=1 sharesome characteristics. For example, suppose the k-th source
signal is a sine wave whose scale varies over time as shownin Figure 2. Note that xkn and xkn (n = n ) have different scales but have the same period. We factorize xkn into
nonstationary and stationary factors as xkn = πkn φkn ,where πkn is a coefﬁcient (scale) of a local signal φkn extracted from the k-th stationary basis signal. Note that thebasis signal is assumed to vary over time according to stationary characteristics (e.g., periodicity and whiteness).
Given O as observed data, we aim to estimate a set oflatent signals {xkn }N
n=1 for each k. The k-th source signal
can be obtained by the overlap-add synthesis method [16].
We do not need any frequency analysis such as short-term
Fourier transform (STFT) or inverse STFT.

c

where = represents equality except for the constant terms.
Given a tensor X = [X 1 , · · · , X N ] ∈ RM ×M ×N , weaim to estimate H ∈ RK×N and V =[V 1 , · · · , V K ] ∈RM ×M ×K such that the log-likelihood n log p(X n |Y n )is maximized. As shown in Section 4, this is a special caseof PSDTF in which each X n is restricted to a rank-1 PSD
matrix (X n = xn xTn ). We can therefore use a multiplicative update algorithm described in Section 4.3.
3.3 Probabilistic Decomposition
After H and V are obtained, we can estimate a local signal
xkn = πkn φkn in a probabilistic manner. Instead of estimating φkn , we can directly calculate a Gaussian posteriorof xkn whose mean and covariance are given by
E[xkn |xn , H, V ] = Y nk Y −1
n xn ,

V[xkn |xn , H, V ] = Y nk − Y nk Y −1n Y nk ,

371

(18)
(19)

14th International Society for Music Information Retrieval Conference (ISMIR 2013)
where Y nk = Hkn V k  0 such that Y n = k Y nk  0.
To evaluate the reconstruction error between PSD maEq. (18) works as a Wiener ﬁlter that passes only a comtrices X n and Y n , we propose to use a Bregman matrix
divergence [14] deﬁned as follows:ponent signal of xn matching the characteristics of kernel V k without explicitly considering the phase and wave
Cφ (X n |Y n )shape. Eqs. (18) and (19) formulated in the time domain= φ(X n ) − φ(Y n ) − tr ∇φ(Y n )T (X n − Y n ) , (22)look similar in form to Eqs. (5) and (12) formulated in thefrequency domain. The key difference is that we consider
the covariance structure over frequency components (e.g.,harmonic partials) when decomposing xn .

where φ is a strictly convex matrix function. In this paper
we focus on the log-determinant (LD) divergence (φ(Z) =− log |Z|) [17] given by + tr X n Y −1
CLD (X n |Y n ) = − log X n Y −1
− M.(23)nn

3.4 Frequency-Domain Representation

This divergence is always nonnegative and is zero if andonly if X n = Y n . The Itakura-Saito (IS) divergence overnonnegative numbers given by CIS (x|y) = − log(x/y) +x/y − 1 is a well-known special case when M = 1, and it
is often used for audio source separation.
Our goal is to estimate H = [h1 , · · · , hK ] ∈ RN ×K
and V = [V 1 , · · · , V K ]∈ RM ×M ×K such that the costfunction CLD (X|Y ) =n CLD (X n |Y n ) is minimized.
Note that our model imposes the nonnegativity constraint
on H and the positive semideﬁniteness constraint on V .
This speciﬁc model is called LD-PSDTF.

We discuss a frequency-domain representation (Figure 1).
Let F ∈ CM ×M be the DFT matrix. Eq. (16) means thatthe complex spectrum F xn (linear transformation of xn )is complex-Gaussian distributed as follows:
 K


H
F xn |H, V ∼ Nc 0,, (20)
Hkn F V k F
k=1

where the full covariance structure between frequency bins
is considered. Note that F V k F H becomes a diagonal matrix if V k is a circulant matrix. A trivial example is a casethat V k is an identity matrix, i.e., φkn is stationary white
Gaussian noise. If V k is a periodic kernel and its size M
is much larger than its period, V k can be roughly viewed
as a circulant matrix. If F V k F H is diagonal, Eq. (20) reduces to a probabilistic model of IS-NMF discarding thecovariance structure between frequency bins [13]. In reality, however, V k is considerably different from a circulant
matrix as shown in Section 5 (Figure 4 and Figure 5).

4.2 Auxiliary Function Approach
We use the auxiliary function approach [15] for indirectly
maximizing the cost function CLD (X|Y ) with respect to Y
(i.e., H and V ) because of its analytical tractability. Let
F(θ) is an objective function to be minimized with respect
to a parameter θ. A function F + (θ, φ) satisfying
F(θ) ≤ F + (θ, φ)

4. POSITIVE SEMIDEFINITE TENSOR
FACTORIZATION

is an auxiliary function for F(θ), where φ is an auxiliary
parameter. It can be proved that F(θ) is non-increasing
through the following iterative update rules:

This section explains a new tensor factorization techniquecalled positive semideﬁnite tensor factorization (PSDTF),
in a general-purpose way. As NMF approximates N nonnegative vectors (a matrix) as the convex combinations of
K nonnegative vectors, PSDTF approximates N PSD matrices (a tensor) as the convex combinations of K PSD matrices. PSDTF is therefore a natural extension of NMF.

φnew ← argminφ F + (θ old , φ),θWe formalize the problem of PSDTF. Suppose we have as
observed data a three-mode tensor X = [X 1 , · · · , X N ] ∈RM ×M ×N , where each slice X n ∈ RM ×M is a real symmetric positive semideﬁnite (PSD) matrix. Note that a similar discussion can be applied even if X n ∈ CM ×M is acomplex Hermitian PSD matrix such that X n = X Hn.
The goal of PSDTF is to approximate each PSD matrix
X n by a convex combination of PSD matrices {V k }K
k=1
(K basis matrices) as follows:
Xn ≈

def

Hkn V k = Y n ,new+← argminθ F (θ, φ

new

).

(25)
(26)

This algorithm is theoretically guaranteed to converge andin fact a similar idea was used for IS-NMF [18].
To derive an auxiliary function U(X|Y ) for CLD (X|Y ),we need to use matrix-variate inequalities based on function concavity and convexity. First, since f (Z) = log |Z|
is concave, we calculate a tangent plane of f (Z) by using
a ﬁrst-order Taylor expansion as follows:

4.1 Problem Speciﬁcation

K


(24)

log |Z| ≤ log |Ω| + tr(Ω−1 Z) − M,(27)

where Ω is an auxiliary PSD matrix (tangent point), M isthe size of Z, and the equality holds when Ω = Z. For
a convex function g(Z) = tr(Z −1 A) for any PSD matrix
A, we can use a sophisticated inequality [19] as follows:
tr

(21)

k=1

K

k=1 Z k

−1

A

≤

K

k=1

T
tr Z −1k Φk AΦk , (28)K
where {Z k }K
k=1 is a set of arbitrary PSD matrices, {Φk }k=1where Hkn ≥ 0 is a weight at the n-th slice. Eq. (21) can
defalso be represented as X ≈ k hk ⊗ V k = Y , where ⊗indicates the Kronecker product.

is a set ofauxiliary matrices that sum to the identity matrix (i.e.,k Φk = I), and the equality holds when Φk =
Z k ( k Z k )−1 .

372

14th International Society for Music Information Retrieval Conference (ISMIR 2013)30 [dB]

Using Inequalities (27) and (28), we can derive an auxiliary function U(X n |Y n ) for Eq. (23) as follows:c
CLD (X n |Y n ) = log |Y n | + tr X n Y −1n−M
≤ log |Ωn | + tr Y n Ω−1
n

T
(29)
+ k tr Y −1nk Φnk X n Φnk−1≤ log |Ωn | + k tr Hkn V k Ωn − M

def
−1 −1
V k Φnk X n ΦTnk = U(X n |Y n ),+ k tr Hkn

3.7 dB25
3.9 dB
20

10

(30)

n=1

N


We used three mixture audio signals each of which wassynthesized using piano sounds (011PFNOM), electric guitar sounds (131EGLPM), or clarinet sounds (311CLNOM)
recorded in the RWC Music Database: Musical Instrument
Sound [20]. Each mixture signal was made by concatenating seven 2.0-s isolated or mixture sounds (C4, E4, G4,
C4+E4, C4+G4, E4+G4, and C4+E4+G4). The resulting
14.0-s signals were sampled at 16kHz.
The task was to separate each mixture signal into threesource signals corresponding to C4, E4, and G4. The local
signals {xn }N
n=1 were extracted by using a Gaussian window with a width of 512 samples (M = 512) and a shifting
interval of 160 samples (N = 1400). The PSD matrices V
and their activations H were estimated by using the MU
algorithm with K = 3. For comparison, we used KLNMF for amplitude-spectrogram decomposition and ISNMF for power-spectrogram decomposition (Section 2.3
and Section 2.4). The number of iterations was 100 in eachmethod. We also tested the iterative STFT method [6] asa phase reconstructor for NMF. We evaluated the qualityof separated signals in terms of source-to-distortion ratio(SDR), source-to-interferences ratio (SIR), and sources-toartifacts ratio (SAR) using the BSS Eval toolbox [21].

(32)

n=1

Eq. (32) can be solved analytically by using the Cholesky
decomposition Qk = Lk LTk , where Lk is a lower triangular matrix. Finally, we get the following update rule:
1

SAR

5.1 Experimental Conditions−1
Hkn Y −1n X n Y n . (33)V k ← V k Lk (LTk V k P k V k Lk )− 2 LTk V k ,SIRThis section reports a comparative experiment evaluating
the source separation performance of LD-PSDTF.

where P k and Qk are PSD matrices given by
Hkn Y −1
n , Qk =SDR

5. EVALUATION

Then, letting the partial derivative with respect to V k beequal to be zero and using Eq. (30), we getN

KL-NMF
KL-NMF with iterSTFT
IS-NMF
IS-NMF with iterSTFT
LD-PSDTF (proposed)

To use LD-PSDTF for source separation, we set X n =xn xTn (rank-1 matrix) as shown in Section 3.2. Since thenegative of the log-likelihood given by Eq. (17) is equal
to the cost function given by Eq. (23) except for constantterms, the maximum-likelihood estimates of H and V canbe obtained by the MU algorithm for LD-PSDTF. Note thatgeneral LD-PSDTF is formulated for any-rank matrix X n .

We can derive multiplicative update (MU) rules that monotonically decrease the total auxiliary function U(X|Y ) =n U(X n |Y n ). We here assume tr(V k ) = 1 (unit trace)
to remove the scale arbitrariness. If tr(V k ) = s, the scale
adjustments V k ← 1s V k and Hkn ← sHkn do not change
CLD (X n |Y n ) and U(X n |Y n ). Letting the partial derivative of Eq. (29) with respect to Hkn be equal to be zero andusing Eq. (30), we get the following update rule:
−1
tr Y −1n V kY n Xn


.
(31)
Hkn ← Hkn
tr Y −1n Vk

Pk =1.4 dB

Figure 3. Source separation performance.

4.3 Multiplicative Updateold
V k P k V k = V old
k Qk V k ,The iterative STFT method
degraded the quality.

1.3 dB

15

where Ωn is a PSD matrix and {Φnk }K
k=1 is a set of auxiliary matrices that sum to the identity matrix. The equality
holds, i.e., U(X n |Y n ) is minimized, whenΩn = Y n , Φnk = Y nk Y −1n .

4.1 dB

1.8 dB

LD-PSDTF achieved thesignificant improvements
over KL-NMF and IS-NMF.

(34)

where the positive semideﬁniteness of V k is ensured. Notethat a real matrix A is said to be positive semideﬁnite if andonly if A = ZZ T is satisﬁed for some real matrix Z.
4.4 Connection to IS-NMF and Source Separation
LD-PSDTF reduces to IS-NMF if PSD matrices X n andV k are restricted to diagonal matrices. Since the diagonal
elements of an arbitrary PSD matrix are always nonnegative, the cost function
 given by Eq. (23) is decomposed as
CLD (X n |Y n ) = m CIS (Xnmm |Ynmm ) and the MU rules
given by Eq. (31) and Eq. (34) thus reduce to the MU rulesof IS-NMF [15]. As described in [15], several algorithms
such as an expectation-maximization (EM) algorithm anda convergence-guaranteed MU algorithm can be used for
IS-NMF. The same is true for LD-PSDTF, and for fasterconvergence we derived the MU algorithm.

5.2 Experimental Results
The experimental results showed the clear superiority of
LD-PSDTF for source separation (Figure 3). The average
SDR, SIR, and SAR were 17.7 dB, 22.2 dB, and 19.7 dB in
KL-NMF, 19.1 dB, 24.0 dB, and 21.0 dB in IS-NMF, and23.0 dB, 27.7 dB, and 25.1 dB in LD-PSDTF.1 We foundthat the iterative STFT method degraded the quality of separated signals. This implies that the spectrogram consistency does not always lead to the perceived quality of au1 Audio ﬁles and a sample code are at the ﬁrst author’s website.

373

14th International Society for Music Information Retrieval Conference (ISMIR 2013)
Basis matrix

Observed mixure signal

0

2

4

Observed matrix

68

Observed matrix

1012

Observed matrix

14[s]Observed matrix

128

128

256

256

256

384

384

384

512

512

0

Basis matrix

Basis matrix

128

128

128

256

256

256

384

Period of C4

512

128

384

512

Activation vector

0.10

0

256

0

2

4 6

8 10 12 14384

Period of E4

512

128

384

512

Activation vector

0.14

0

256

0

2

4 6

8 10 12 14

384

128

384

512

Activation vector

0.35

0

256

0

2

4 6

384

512

128

256

384

512

Activation vector
0.06

0

2

4 6

8 10 12 14

0

0

2

4 6

8 10 12 14

512

128

0

256

384

512

Activation vector

0.08

0

2

4 6

2

8 10 12 14 x10[2] M. Goto. Active music listening interfaces based on signal
processing. ICASSP, volume 4, pp. 1441–1444, 2007.
[3] K. Yoshii et al. Drumix: An audio player with real-timedrum-part rearrangement functions for active music listening. IPSJ Digital Courier, 3:134–144, 2007.
[4] N. Sturmel et al. Linear mixing models for active listeningof music productions in realistic studio conditions. AES Convention, 2012.
[5] D. Lee and H. Seung. Algorithms for non-negative matrix
factorization. NIPS, pp. 556–562, 2000.
[6] D. W. Grifﬁn and J. S. Lim. Signal estimation from modiﬁed short-time Fourier transform. IEEE Trans. on ASLP,
32(2):236–243, 1984.
[7] J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama. Explicit
consistency constraints for STFT spectrograms and their application to phase reconstruction. SAPA, pp. 23–28, 2008.
[8] J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama. Fast signal reconstruction from magnitude STFT spectrogram basedon spectrogram consistency. DAFx, pp. 397–403, 2010.
[9] H. Kameoka et al. Complex NMF: A new sparse representation for acoustic signals. ICASSP, pp. 45–48, 2009.
[10] J. Le Roux et al. Consistent Wiener ﬁltering: Generalized
time-frequency masking respecting spectrogram consistency.
LVA/ICA, pp. 89–96, 2010.
[11] K. Yoshii, R. Romioka, D. Mochihashi, and M. Goto. Inﬁnite
positive semideﬁnite tensor factorization for source separation of mixture signals. ICML, pp. 576–584, 2013.
[12] P. Smaragdis and J. Brown. Nonnegative matrix factorization
for polyphonic music transcription. WASPAA, 2003.
[13] C. Févotte, N. Bertin, and J.-L. Durrieu. Nonnegative matrix
factorization with the Itakura-Saito divergence: With application to music analysis. Neu. Comp., 21(3):793–830, 2009.
[14] L. M. Bregman. The relaxation method of ﬁnding the common points of convex sets and its application to the solutionof problems in convex programming. USSR CMMP, 1967.
[15] M. Nakano et al. Convergence-guaranteed multiplicative algorithms for non-negative matrix factorization with beta divergence. MLSP, pp. 283–288, 2010.
[16] J. Allen and L. Rabiner. A uniﬁed approach to short-time
Fourier analysis and synthesis. IEEE, 1977.
[17] B. Kulis, M. Sustik, and I. Dhillon. Low-rank kernel learningwith Bregman matrix divergences. JMLR, 10:341–376, 2009.
[18] M. Hoffman et al. Bayesian nonparametric matrix factorization for recorded music. ICML, pp. 439–446, 2010.
[19] H. Sawada, H. Kameoka, S. Araki, and N. Ueda. Efﬁcientalgorithms for multichannel extensions of Itakura-Saito nonnegative matrix factorization. ICASSP, pp. 261–264, 2012.
[20] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka. RWC
music database: Music genre database and musical instrument sound database. ISMIR, pp. 229–230, 2003.
[21] E. Vincent, R. Gribonval, and C. Févotte. Performance measurement in blind audio source separation. IEEE Trans. on
ASSP, 14(4):1462–1469, 2006.
[22] H. Kameoka and K. Kashino. Composite autoregressive system for sparse source-ﬁlter representation of speech. ISCAS,
pp. 2477–2480, 2009.

Period of G4

512

256

Activation vector

Figure 5. Factorization of a clarinet mixture signal.

Basis matrix

128

Basis matrix

128

0.06

Discover basis matrices such that the log-determinant divergence is minimized

Basis matrix

2

8 10 12 14 x10

Figure 4. Factorization of a piano mixture signal.
dio signals, as suggested in [10]. We conﬁrmed that LDPSDTF can appropriately estimate V and H from both decaying and sustained sounds (Figure 4 and Figure 5). The
reason that each X n does not appear to be well approximated by Y n (a convex combination of {V k }Kk=1 ) is thatthe cost function based on the LD divergence allows Y n tooverestimate X n with a smaller penalty. A main limitationof LD-PSDTF is that its computational cost is O(KN M 3 )while the computational cost of NMF is O(KN M ). In thisexperiment, LD-PSDTF spent several hours for analyzingeach mixture signal on Xeon X5492 (3.4 GHz). Therefore,
we think that initializing LD-PSDTF by using basis vectors and their activations obtained by IS-NMF can reducethe computational cost and help avoid local minima.
6. CONCLUSION
This paper presented log-determinant positive semideﬁnite
tensor factorization (LD-PSDTF) as a natural extension of
Itakura-Saito NMF (IS-NMF). We derived a convergenceguaranteed multiplicative update algorithm and showed theclear superiority of LD-PSDTF over NMF variants in termsof source separation quality.
There are several interesting directions. To separate music signals into instrument parts, we plan to fuse the sourceﬁlter model into the framework of LD-PSDTF as in thecomposite autoregressive system [22]. We also plan to investigate another variant of PSDTF based on the von Neumann divergence (φ(Z) = tr(Z log Z − Z) in Eq. (22))
that can be viewed as an extension of KL-NMF.
Acknowledgment: This study was supported in part by JSPS
KAKENHI 23700184, MEXT KAKENHI 25870192, and JST
CREST OngaCREST.

7. REFERENCES[1] K. Itoyama et al. Query-by-example music information retrieval by score-informed source separation and remixing
technologies. EURASIP Journal, 2010. Article ID 172961.

374