STUDENT’S T NONNEGATIVE MATRIX FACTORIZATION AND POSITIVE SEMIDEFINITETENSOR FACTORIZATION FOR SINGLE-CHANNEL AUDIO SOURCE SEPARATION
Kazuyoshi Yoshii1

Katsutoshi Itoyama1

Masataka Goto2

1

2

Graduate School of Informatics, Kyoto University, Japan
National Institute of Advanced Industrial Science and Technology (AIST), Japan{yoshii, itoyama}@kuis.kyoto-u.ac.jp m.goto@aist.go.jp
Approximated byweighted sums ofbasis matrices

ABSTRACTThis paper presents a robust variant of nonnegative matrix factorization (NMF) based on complex Student’s t distributions (t-NMF) forsource separation of single-channel audio signals. The Itakura-Saito
divergence NMF (Gaussian NMF) is justified for this purpose under
an assumption that the complex spectra of source signals and thoseof the mixture signal are complex Gaussian distributed (the additivity of power spectra holds). In fact, however, the source spectra areoften heavy-tailed distributed. When the source spectra are complex
Cauchy distributed, for example, the mixture spectra are also complex Cauchy distributed (the additivity of amplitude spectra holds).
Using the complex t distribution that includes the complex Gaussian
and Cauchy distributions as its special cases, we propose t-NMF asa unified extension of Gaussian NMF and Cauchy NMF. Furthermore, we propose the corresponding variant of positive semidefinite
tensor factorization based on multivariate complex t distributions (tPSDTF). The experimental results showed that while t-NMF and tPSDTF were comparative to Gaussian counterparts in terms of peak
performance, they worked much better on average because they areinsensitive to initialization and tend to avoid local optima.

: Complex spectrum at frame
Focus on only diagonal elements(a set of power spectra)[dB]

A set of PSD matrices

A set of nonnegative vectors

Approximated byweighted sumsof basis vectors

Fig. 1. Comparison between PSDTF and NMF.

One of the most standard approaches to source separation of singlechannel audio signals is to perform nonnegative matrix factorization(NMF) and Wiener filtering in the frequency domain [1–9]. NMF
approximates a nonnegative matrix (a set of nonnegative vectors) as
the product of two nonnegative matrices (a set of basis vectors anda set of activation vectors), i.e., each nonnegative vector is approximated by the weighted sum of nonnegative basis vectors. If NMF
is applied to a nonnegative spectrogram (e.g., an amplitude or power
spectrogram) of piano sounds (mixture), for example, the basis vectors are expected to be average energy spectra of different pitches
used in the piece. The mixture spectrogram is then decomposed into
the sum of source spectrograms via time-frame-and-frequency-binwise Wiener filtering (time-frequency bins are processed independently) according to the source proportions determined by the basis
and activation vectors. However, it is difficult to resynthesize highquality time-domain source signals because the phase information ofsource spectrograms should be recovered by post-processing [10].
To circumvent this problem, positive semidefinite tensor factorization (PSDTF) has recently been proposed for audio source separation [11, 12]. Given a set of positive semidefinite (PSD) matrices

as input data, each PSD matrix is approximated by the weighted sumof PSD basis matrices. As shown in Fig. 1, PSDTF is a mathematically fundamental extension of NMF because positive semidefiniteness of matrices is an extended concept of nonnegativity of scalars
and vectors. Each PSD matrix of input data is obtained by calculating the product of the complex spectrum and its conjugate transpose
in a time frame (window). NMF, on the other hand, focuses on only
the diagonal elements of the PSD matrix, i.e., on nonnegative vectors (power spectra). Since PSDTF can deal with phase information,
the complex spectrograms of source signals are directly obtained via
time-frame-wise Wiener filtering (the frequency bins of each frame
are processed jointly in an interdependent manner).
Among many variants of NMF and PSDTF, Itakura-Saito divergence NMF (IS-NMF or Gaussian NMF) [13] or LogDet divergence
PSDTF (LD-PSDTF or Gaussian PSDTF) [11] is justified for audio source separation under an assumption that the complex spectraof source signals are univariate or multivariate complex Gaussian
distributed, respectively. The underlying probabilistic models havecomplex Gaussian likelihoods for the observed mixture spectra dueto the reproductive property of the Gaussian distribution (the additivity of power spectra). In fact, however, Gaussian NMF often underperforms Kullback-Leibler divergence NMF (KL-NMF or Poisson
NMF) [14] that assumes the additivity of amplitude spectra based on
the physically-meaningless Poisson likelihood [15, 16].
The univariate complex symmetric α-stable (SαS) distribution
with a characteristic-exponent parameter 0 < α ≤ 2 was recently
found to justify the additivity of fractional power spectra1 [17]. It
has the reproductive property for any α and includes the univariate
complex Gaussian (α = 2) and Cauchy (α = 1) distributions as its
special cases. If the source spectra are complex SαS distributed, themixture spectra are also complex SαS distributed. Under this generative process, generalized Wiener filtering was proposed for decomposing the mixture spectra into the sum of source spectra in termsof posterior inference [17]. Maximum-likelihood estimation of SαSdistributions of source signals, however, has been proposed only forThis study was partially supported by JST OngaCREST Project, JSPS
KAKENHI 24220006, 26700020, and 26280089, and Kayamori Foundation.

1 Element-wise power-of-α of amplitude spectra, e.g., power spectra (α =2) and amplitude spectra (α = 1).

Index Terms— Source separation, nonnegative matrix factorization, positive semidefinite tensor factorization, t distribution.
1. INTRODUCTION

978-1-4799-9988-0/16/$31.00 ©2016 IEEE

51

ICASSP 2016

α = 2 (Gaussian NMF [13]) and α = 1 (Cauchy NMF [18]) because the probability density function (PDF) of the SαS distribution
is analytically expressible only for α = 2 or 1.
In this paper we propose Student’s t PSDTF (t-PSDTF), which
includes Student’s t NMF (t-NMF), Gaussian NMF (Gaussian NMF[13]), Cauchy NMF [18], Gaussian PSDTF (Gaussian PSDTF [11]),and Cauchy PSDTF as its special cases. The complex t distribution
has an analytically expressible PDF for any degree-of-freedom parameter ν and includes the complex Gaussian (ν = ∞) and Cauchy(ν = 1) distributions as its special cases (Fig. 2). Since the t distribution, like the SαS distribution, has heavy tails, it is considered
to be more suitable for modeling real-world audio signals than theGaussian distribution. Although the additivity of fractional power
spectra is not generally justified, i.e., the reproductive property of
the univariate or multivariate complex t distribution holds only forν = 1 or ν = ∞, respectively, other values of ν have the potentialof performing best in practice. The main contributions of this paper
are to propose Cauchy PSDTF as an extension of Cauchy NMF [18]and to formulate a unified probabilistic model that enables us to continuously adjust the value of ν.
To execute t-PSDTF, we derive a convergence-guaranteed multiplicative update algorithm that maximizes the t likelihood for a mixture spectrogram by using an auxiliary function technique [19, 20].
The parameter updating rules of t-PSDTF and t-NMF are found toconverge to those of Gaussian PSDTF [12] and Gaussian NMF [19]when ν → ∞ and to give new updating rules of Cauchy NMF when
ν = 1. The unified view of these updating rules reveals why t-PSDTF
based on the heavy-tailed t distribution is more robust to outliers andavoids over-fitting to the observed mixture spectrogram.
2. MATHEMATICAL FOUNDATIONThis section explains the probabilistic models of Gaussian NMF [13]and Gaussian PSDTF [11] and introduces the characteristics of the
α-stable and Student’s t distributions.

def

hkn wk = yn ,

(3)

k=1
×K
is a set of activation vectors,
where H = [h1 , · · · , hK ] ∈ RN+
K  min(M, N ) is the number of basis vectors, and Yn  0 ∈CM ×M represents a reconstruction matrix. Let Y
kn = hkn Wk bea source reconstruction matrix such that Yn =k Ykn . In LDPSDTF [12] the reconstruction error D(Xn |Yn ) between Xn and
Yn is evaluated using the LogDet divergence [21]:




(4)
DLD (Xn |Yn ) = tr Xn Yn−1 − log Xn Yn−1  − M.

This divergence is never less than zero and is zero only when Xn =Yn . To estimate W and H such that the cost function D(X |Y) =n D(Xn |Yn ) is minimized, a convergence-guaranteed multiplicative update (MU) algorithm was derived [12].
The similarity between Eq. (1) and Eq. (3) indicates that Gaussian PSDTF reduces to Gaussian NMF when all PSD matrices are
restricted to diagonal matrices (diag(Xn ) = xn , diag(Yn ) = yn ,diag(Ykn ) = ykn , and diag(Wk ) = wk ). Note that the diagonal
elements of PSD matrices always take nonnegative values.
2.3. Application to audio source separation
To formulate a probabilistic model for audio source separation, it isnecessary to represent the generative process of the complex spectrogram of a target mixture signal. Let Sk = [sk1 , · · · , skN ] ∈CM ×N be the complex spectrogram of the k’th source signal, where
M is the number of frequency bins (window size) and N is the number of frames. Let S = [s1 , · · · , sN ] ∈ CM ×N be the complex
spectrogram of the mixture signal such that S = K
k=1 Sk .
In Gaussian PSDTF, skn is assumed to be multivariate complex
Gaussian distributed with a covariance matrix Ykn as follows:(6)

More specifically, the log-likelihood function is given by−1
log p(sn |Yn ) = −M log(π) − log |Yn | − sHn Yn sn= − log |Yn | − tr(Xn Yn−1 ) = L(Xn |Yn )
c

(1)

def

(7)

where Xn = sn sH
n is an observed PSD matrix obtained by calcu-

×K
is a set of activation vectors,
where H = [h1 , · · · , hK ] ∈ RN+
K  min(M, N ) is the number of basis vectors, and yn ∈ RM+
represents a reconstruction vector. Let y
kn = hkn wk be a sourcereconstruction vector such that yn =k ykn . In IS-NMF [13]the reconstruction error D(xn |yn ) between xn and yn is evaluatedusing the Itakura-Saito divergence:
M



−1
−1
− log xnm ynm−1
xnm ynmdef

hkn Wk = Yn ,sn ∼ Nc (0, Yn ).

k=1

DIS (xn |yn ) =

K


(5)
skn ∼ Nc (0, Ykn ).
KSince S = k=1 Sk and Yn = k=1 Ykn , the reproductive property of the Gaussian distribution leads to×N
Given a set of nonnegative vectors X = [x1 , · · · , xN ] ∈ RM+
as input data, NMF approximates each nonnegative vector xn ∈RM
+ by the weighted sum of a limited number of nonnegative basis
×K
vectors W = [w1 , · · · , wK ] ∈ RM
as follows:
+
K


Xn ≈

K

2.1. Nonnegative matrix factorizationxn ≈Xn  0 ∈ CM ×M by the weighted sum of a limited number of PSD
basis matrices W = [W1 , · · · , WK ] ∈ CM ×M ×K as follows:

lating the product of the complex spectrum and its conjugate transpose. The maximum-likelihood estimates of W and
H are obtainedby maximizing the total log-likelihood L(X |Y) = n L(Xn |Yn ).
To do this, Gaussian PSDTF can be used because maximization of
Eq. (7) is equivalent to minimization of Eq. (4).
Given the mixture spectrogram S, each source spectrogram S k
can be estimated via Wiener filtering, i.e., calculation of posterior
Gaussians of source spectra based on estimated W and H.

(2)

E[skn |sn ] = Ykn Yn−1 sn .

m=1

(8)

This divergence is never less than zero and is zero only when xn =yn .To estimate W and H such that the total cost function D(X|Y )=n D(xn |yn ) is minimized, a convergence-guaranteed multiplicative update (MU) algorithm was derived [19].

In Gaussian NMF, in contrast, the elements of skn are assumed
to be independently distributed, i.e., sknm ∼ Nc (0, yknm ) and snm∼ Nc (0, ynm ). This means that the correlations between frequency
bins are ignored, resulting in inferior quality of source separation.

2.2. Positive semidefinite tensor factorization

2.4. Elliptically contoured α-stable distribution

Given a set of complex-valued PSD matrices X = [X1 , · · · , XN ] ∈CM ×M ×N as input data, PSDTF approximates each PSD matrix

The elliptically contoured (sub-Gaussian) version of the multivariate(d-variate) complex α-stable distribution, Sαc (μ, Σ), is specified by52

a location vector μ ∈ Cd and a PSD scale matrix Σ  0 ∈ Cd×d ,and Sαc (0, Σ) is called the symmetric α-stable (SαS) distribution[22]. When α = 2 and α = 1, it reduces to the complex Gaussian
and Cauchy distributions given by Nc (μ, 2Σ) and Cc (μ, Σ), where1exp −z H Σ−1 z ,(9)Nc (x|μ, Σ) = d
π |Σ|− 2d+12d Γ( 2d+1 ) 2
1 + 2z H Σ−1 z.
(10)
Cc (x|μ, Σ) = 2d+1 2
π 2 |Σ|where z = x − μ. While the PDF of a sub-Gaussian α-stable random vector x ∈ Cd cannot be given in closed form, its characteristic
function (CF; Fourier transform of PDF) is analytically given by
 

α
H
ϕ(t) = E eiRe(t X) = exp − 12 tH Σt 2 + iRe(tHμ) . (11)

0.2
0.1
-4

-3

-2

-1

0

1

2

3

4

5

-5

-4

-3

-2

-1

0

1

5

Instead of the multivariate complex Gaussian log-likelihood given by
Eqs. (6) and (7), we use the multivariate complex t log-likelihood:c
L(Xn |Yn ) = − log |Yn | − 2M2+ν log 1 + ν2 tr Xn Yn−1 . (15)cwhere = denotes equality up to a constant and skn ∼ Tνc (0, Ykn )⇒ sn ∼ Tνc (0, Yn ) does not hold for finite ν (cf. Eqs. (5) and (6)).
To derive convergence-guaranteed updating rules
 of W and H
that maximize the total log-likelihood L(X |Y) = n L(Xn |Yn ),we use an auxiliary function technique [19] that maximizes the lower
bound of L(X |Y), F (X |Y). First, for a convex function f (Z) =− log |Z| (Z  0 ∈ CM ×M ), we calculate a tangent plane at arbitrary Ω  0 by using a first-order Taylor expansion as follows:

(12)

(16)

where the equality holds when Ω = Z. Second, for a concave function g(Z) = −tr(Z −1 A) with any PSD matrix A  0, we use thefollowing inequality [24]:
K
−1




K
,
(17)
A ≥−
tr Zk−1 Λk AΛH−trk
k=1 Zkified by a degree of freedom ν, a location vector μ ∈ C , and a PSD
scale matrix Σ  0 ∈ Cd×d . Its PDF is explicitly given by− 2d+ν
2

4

3.1. Maximum-likelihood estimation− log |Z| ≥ − log |Ω| − tr(Ω−1 Z) + M,

2 H −1z Σ zν3Fig. 2. PDFs of univariate real SαS and Student’s t distributions; thescale parameters of the SαS distribution are adjusted for comparison,considering S2 (0, 0.5) = T∞ (0, 1) and S1 (0, 1) = T1 (0, 1).

The multivariate complex student’s t distribution Tνc (μ, Σ) is specd

1+

2Standard distribution:

S S distribution:

2.5. Student’s t distribution

2d Γ( 2d+ν)2νd(πν) Γ( 2 )|Σ|

0.3

0.2

-5

x1 ∼ Nc (μ1 , Σ1 )⇒ x1 +x2 ∼ Nc (μ1 +μ2 , Σ1 +Σ2 ) . (13)x2 ∼ Nc (μ2 , Σ2 )Tνc (x|μ, Σ) =

0.4

0.3

0.1

This indicates that while the reproductive property holds for any α
in the univariate case (d = 1), it holds only for α = 2 in the multivariate case (d > 1)2 , i.e.,x1 ∼ Sc (μ1 , σ1α )⇒ x1 + x2 ∼ Sc (μ1 + μ2 , σ1α + σ2α ) ,x2 ∼ Sc (μ2 , σ2α )

0.4. (14)

k=1

Like the α-stable distribution, it converges to Nc (μ, Σ) as ν → ∞and reduces to Cc (μ, Σ) when ν = 1. The reproductive property
holds only when ν = 1 and d = 1.
For comparison, Fig. 2 shows the PDFs of univariate real SαSand Student’s t distributions, Sα (0, σ α ) and Tν (0, σ 2 ), with different values of α and ν, where S2 (0, σ 2 ) = T∞ (0, 2σ 2 ) = N (0, 2σ 2 )and S1 (0, σ) = T1 (0, σ 2 ) = C(0, σ). Since both distributions haveheavy tails, in this paper we focus on the t distribution as a substitutefor the α-stable distribution because of its tractability.

K
where {Zk  0}K
k=1 is a set of arbitrary PSD matrices, {Λk }k=1 isa set of auxiliary matrices that sum to the identity
 matrix ( k Λk =I), and the equality holds when Λk = Zk ( k Zk )−1 .
Using Eqs. (16) and (17), F (X |Y) can be derived as follows: 
L(X |Y) ≥ n − log |Ωn | − tr(Ω−1n Yn ) + M





− 2M2+ν ψn + ψn−1 1 + ν2 tr Xn Yn−1 − 1 (18)c    −1 −1
H
2M +ν −1
≥ n −tr(Ω−1ψnn Yn ) −k tr hkn Wk Λkn Xn Λkn ,ν

3. PROPOSED METHOD

where the equality holds (F (X |Y) is maximized) when2 
Ωn = Yn , ψn = 1 + tr Xn Yn−1 , Λkn = Ykn Yn−1 . (19)ν
Letting the partial derivatives of F (X |Y) with respect to Wk
and hkn be zero, we get the following updating rules:

This section proposes a novel variant of PSDTF based on the multivariate complex t likelihood (t-PSDTF). It has a degree-of-freedom
parameter ν for decomposing the complex spectrogram of a mixture
signal into the sum of the spectrograms of source signals via generalized Wiener filtering [17] as in Gaussian PSDTF [12] (Eq. (8)). If thetime-frequency bins of source spectrograms are independently univariate complex Gaussian or Cauchy distributed (ν = ∞ or 1), theadditivity of power or amplitude spectra is satisfied [23] (Eq. (12)).
This theoretically justifies the generative process of a mixture spectrogram behind Gaussian NMF [13] and Cauchy NMF [18]. If eachframe of source spectrograms is multivariate complex Gaussian distributed (ν = ∞), the additivity of power spectra forms the basisof Gaussian PSDTF [12] (Eq. (13)). Although the additivity of fractional power spectra does not hold in the other cases, the unified
probabilistic formulation of t-PSDTF enables us to flexibly tune νor perform a kind of annealing by gradually increasing ν as an inverse temperature to avoid local optima of Gaussian PSDTF.

1

1

1

1

1

(20)
Wk ← Wk Qk2 (Qk2 Wk Pk Wk Qk2 )− 2 Qk2 Wk ,
 
1

2
tr πn Xn Yn−1 Wk Yn−1hkn ← hkn,
(21)−1
tr Wk Yn

N
−1
−1
−1
where Pk = N
n=1 hkn Yn , Qk =n=1 hkn Yn (πn Xn )Yn ,2M +ν
. If tr(Wk ) = s, we disambiguate theand πn =2tr(Xn Yn−1 )+ν
scales of W and H by Wk ← 1s Wk and hkn ← shkn .
3.2. Time-domain formulation
t-PSDTF can be defined in the time domain. Let Ŝ = [ŝ1 , · · · , ŝN ] ∈RM ×N be a set of windowed signals of N frames. We assume thatŝn ∼ Tν (0, Ŷn ), where Ŷn ∈ RM ×M is a PSD matrix and Tν isa multivariate real t distribution. Let F ∈ CM ×M be the DFT ma-

2 The convolution of two PDFs (the PDF of the sum of two random variables) is equivalent to the product of the corresponding CFs.

53

30

SDR

[dB]

SIR

Observed mixure signal

SAR
0

20

0.5 1 2 5 10 20 ∞ 1 ∞t -NMF0.5 1 2 5 10 20 ∞ 1 ∞t -PSDTF

t -NMF

4

t -NMF

Basis matrix

t -PSDTF

14[s]

384

Period of C4

512

Activation vector




c
L(X̂n |Yˆn ) = − log |Ŷn | − M2+ν log 1 + ν1 tr X̂n Ŷn−1 . (22)0The updating rules that iteratively maximize L(X̂n |Yˆn ) are the same
M +νas Eqs. (20) and (21) except that πn =. This is nottr(Xn Yn−1 )+ν
equivalent to the frequency-domain formulation given in Section 3.1.

2

4

6

Period of E4

512

128 256 384 512

8 10 12 14

Basis matrix
256

256

384

k ĥkn Ŵk . The log-likelihood function is given by

128 256 384 512

384

Period of G4

512

128 256 384 512

Activation vector

Activation vector

2

2

0

4

6

8 10 12 14 0

4

6

8 10 12 14

Fig. 4. Time-domain Cauchy PSDTF for a piano mixture signal.
N
The local signals {ŝn }Nn=1 and complex spectra {sn }n=1 were extracted using a Gaussian window with a width of 512 samples (M =512) and a shifting interval of 160 samples (N = 1400).
We tested Poisson NMF [14], t-NMF with ν = 0.5, 1 (Cauchy
NMF [18]), 2, 5, 10, 20, ∞ (Gaussian NMF [13]), and t-PSDTF
with ν = 1 (Cauchy PSDTF), ∞ (Gaussian PSDTF [12]). Eachvariant of NMF was executed 100 times with random initializationand t-PSDTF was initialized by using the average results of Cauchy
NMF for fast convergence. The separation quality was evaluatedin terms of source-to-distortion ratio (SDR), source-to-interferences
ratio (SIR), and sources-to-artifacts ratio (SAR) [26].

3.3. Connection to Gaussian PSDTF and Gaussian NMF
The updating rules of t-PSDTF given by Eqs. (20) and (21) converge
to those of Gaussian PSDTF [12] as ν → ∞, where πn Xn plays arole in controlling the impact of the observed data Xn in each iteration. Since πn Xn → Xn as ν → ∞, Gaussian PSDTF totally relieson the observed data Xn . Interestingly, we found that t-PSDTF withfinite ν can be interpreted as Gaussian PSDTF that virtually regards
πn Xn as observed data varying over updating iterations, where πndepends on the current reconstruction Yn .
As shown in Section 2.2, the updating rules of t-NMF are obtained by restricting PSD matrices to diagonal matrices (diag(Xn ) =xn , diag(Yn ) = yn , and diag(Wk ) = wk ) as follows:12
2
nm xnm )hkn /ynm
n (π

,
(23)
wkm ← wkm
n hkn /ynm
1
2
2
nm xnm )wkm /ynmm (πhkn ← hkn,
(24)m wkm /ynm
and πnm xnm is given by
where πnm = 2xnm2+ν
/ynm +ν
−1

−1
2
ν
ynm
+ 2+ν
x−1
.
πnm xnm = 2+ν
nm

12

PSD matrix

128

128

256

trix. Given X̂n = ŝn ŝTn = F HXn F as observed data, we can estimate Ŵk = F H Wk F and ĥkn = hkn by using t-PSDTF such that

10

PSD matrix

Basis matrix

128

Fig. 3. Source separation performance.

X̂n ≈ Ŷn =8

PSD matrix

Discover basis PSD matrices such that the multivariate Studentʼs t likelihood is maximized

0.5 1 2 5 10 20 ∞ 1 ∞t -PSDTF

6

PSD matrix

Gaussian

Gaussian
CauchyCauchyPoissonCauchyGaussian

Gaussian

CauchyPoissonCauchyGaussian

Gaussian

CauchyPoisson

10

0

2

PSD matrix

4.2. Experimental results
The experimental results showed the clear superiority of t-NMF with
ν = 1, 2 for source separation (Fig. 3). Although the peak performance of t-NMF with any ν was comparable to Gaussian NMF [13],t-NMF with small ν was robust to initialization and stably attained
good performance. In addition, t-PSDTF significantly outperformed
t-NMF. If t-PSDTF was initialized appropriately, almost the same
performance was achieved regardless of ν. Interestingly, t-NMF
with ν = 2 slightly worked better than Cauchy NMF [18] assuming
the additivity of amplitude spectra in a theoretically-founded way.
This indicates the practical effectiveness of our unified formulation
with a tunable parameter ν. To reduce the prohibitive computational
cost of PSDTF, O(KN M 3 ), we plan to use low-rank approximation of basis matrices. Maximum-likelihood estimation of ν would
be feasible by starting with ν = 1 and gradually increasing ν.

(25)

Eqs. (23) and (24) converge to the updating rules of Gaussian NMF[19] as ν → ∞ and give new updating rules of Cauchy NMF asalternatives to those in [18]. Eq. (25) is the harmonic mean of observation xnm and reconstruction ynm with a ratio of ν to 2. Unlike
Gaussian NMF, t-NMF makes reconstruction ynm close to virtual
observation πnm xnm depending on ynm , not real obeservation xnm .
This prevents t-NMF from over-fitting to xnm .

5. CONCLUSION
This paper presents t-PSDTF, a robust version of positive semidefinite tensor factorization for single-channel audio source separation.
Based on the multivariate complex t likelihood, it includes t-NMF,
Gaussian NMF (IS-NMF [13]), Cauchy NMF [18], Gaussian PSDTF(LD-PSDTF [11]), and Cauchy PSDTF (newly proposed) as its special cases. We found that while the peak performances of t-NMF and
t-PSDTF were on a par with those of Gaussian NMF and PSDTF, tNMF and t-PSDTF were less likely to get stuck in bad local optima.
We will try to find multiplicative updating rules of PSDTF basedon the multivariate complex symmetric α-stable likelihood (SαSPSDTF). Its NMF counterpart called SαS-NMF [27] is theoreticallyan ideal approach to source separation because the additivity of fractional power spectra holds true for any α [17]. Although W , H, αcan be optimized or marginalized out in a Bayesian manner, only a
MCMC-based method can be used at the momemt.

4. EVALUATIONThis section reports a comparative experiment evaluating the source
separation performance of t-PSDTF and t-NMF.
4.1. Experimental conditions
We used three audio recordings each of which was synthesized using
piano sounds (011PFNOM), electric guitar sounds (131EGLPM), or
clarinet sounds (311CLNOM) from the RWC Music Database: Musical Instrument Sound [25]. Each recording (16 kHz, mono) of 14 s
was made by concatenating seven single tones and chords of 2 s (C4,E4, G4, C4+E4, C4+G4, E4+G4, and C4+E4+G4). Each recording
was separated into three source signals of C4, E4, and G4 (K = 3).

54

[14] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization for polyphonic music transcription,” in IEEE Workshop
on Applications of Signal Processing to Audio and Acoustics(WASPAA), 2003, pp. 177–180.

6. REFERENCES[1] P. Smaragdis, C. Févotte, G. Mysore, N. Mohammadiha, andM. Hoffman, “Dynamic source separation using nonnegative
factorizations: A unified view,” IEEE Signal Processing Magazine, vol. 31, no. 3, pp. 66–75, 2014.

[15] D. Fitzgerald, M. Cranitch, and E. Coyle, “On the use of thebeta divergence for musical source separation,” in Irish Signals
and Systems Conference (ISSC), 2009, pp. 1–6.

[2] A. T. Cemgil, “Bayesian inference for nonnegative matrix
factorisation models,” Computational Intelligence and Neuroscience, vol. 2009, pp. 1–17, 2009.

[16] B. King, C. Févotte, and P. Smaragdis, “Optimal cost function
and magnitude power for NMF-based speech separation andmusic interpolation,” in International Workshop on Machine
Learning for Signal Processing (MLSP), 2012, pp. 1–6.

[3] A. Ozerov, A. Liutkus, R. Badeau, and G. Richard, “Codingbased informed source separation: Nonnegative tensor factorization approach,” IEEE Transactions on Audio, Speech, andLanguage Processing, vol. 21, no. 8, pp. 1699–1712, 2013.

[17] A. Liutkus and R. Badeau, “Generalized Wiener filtering withfractional power spectrograms,” in International Conferenceon Acoustics, Speech, and Signal Processing (ICASSP), 2015,
pp. 266–270.

[4] D. Liang, M. Hoffman, and D. Ellis, “Beta process sparse nonnegative matrix factorization for music,” in International Society for Music Information Retrieval Conference (ISMIR), 2013,
pp. 375–380.

[18] A. Liutkus, D. Fitzgerald, and R. Badeau, “Cauchy nonnegative matrix factorization,” in IEEE Workshop on Applicationsof Signal Processing to Audio and Acoustics (WASPAA), 2015.

[5] F. J. Canadas-Quesada, P. Vera-Candeas, N. Ruiz-Reyes, J. J.
Carabias-Orti, and P. Cabanas-Molero, “Percussive/harmonicsound separation by non-negative matrix factorization withsmoothness/sparseness constraints,” EURASIP Journal on Audio, Speech, and Music Processing, vol. 26, no. 1, pp. 1–17,
2014.

[19] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono, and
S. Sagayama, “Convergence-guaranteed multiplicative algorithms for non-negative matrix factorization with beta divergence,” in International Workshop on Machine Learning for
Signal Processing (MLSP), 2010, pp. 283–288.

[6] N. Souviraà-Labastie, E. Vincent, and F. Bimbot, “Music
separation guided by cover tracks: Designing the joint NMF
model,” in International Conference on Acoustics, Speech, andSignal Processing (ICASSP), 2015, pp. 484–488.

[20] M. Hoffman, D. Blei, and P. Cook, “Bayesian nonparametric matrix factorization for recorded music,” in International
Conference on Machine Learning (ICML), 2010, pp. 439–446.
[21] B. Kulis, M. Sustik, and I. Dhillon, “Low-rank kernel learning with Bregman matrix divergences,” Journal of Machine
Learning Research (JMLR), vol. 10, pp. 341–376, 2009.

[7] N. Q. K. Duong D. El Badawy, A. Ozerov, “Relative group
sparsity for non-negative matrix factorization with application
to on-the-fly audio source separation,” in International Conference on Acoustics, Speech, and Signal Processing (ICASSP),
2015, pp. 256–260.

[22] G. Samorodnitsky and M. Taqqu, Stable Non-Gaussian Random Processes: Stochastic Models with Infinite Variance,
Chapman & Hall/CRC, 1994.

[8] K. O’Hanlon, M. Sandler, and M. D. Plumbley, “Matrix factorisation incorporating greedy Hellinger sparse coding applied to polyphonic music transcription,” in International Conference on Acoustics, Speech, and Signal Processing(ICASSP), 2015, pp. 3112–3116.

[23] G.A. Tsihrintzis, P. Tsakalides, and C.L. Nikias, “Spectral
methods for stationary harmonizable alpha-stable processes,”in European Signal Processing Conference (EUSIPCO), 1998,
pp. 1833–1836.
[24] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Efficientalgorithms for multichannel extensions of Itakura-Saito nonnegative matrix factorization,” in International Conference on
Acoustics, Speech, and Signal Processing (ICASSP), 2012, pp.
261–264.

[9] C.-W. Wu and A. Lerch, “Drum transcription using partially
fixed non-negative matrix factorization with template adaptation,” in International Society for Music Information Retrieval
Conference (ISMIR), 2015.
[10] J. Le Roux, E. Vincent, Y. Mizuno, H. Kameoka, N. Ono, and
S. Sagayama, “Consistent Wiener filtering: Generalized timefrequency masking respecting spectrogram consistency,” in International Conference on Latent Variable Analysis and Signal
Separation (LVA/ICA), 2010, pp. 89–96.

[25] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka, “RWC
music database: Music genre database and musical instrument
sound database,” in International Conference on Music Information Retrieval (ISMIR), 2003, pp. 229–230.

[11] K. Yoshii, R. Tomioka, D. Mochihashi, and M. Goto, “Beyond
NMF: Time-domain audio source separation without phase reconstruction,” in International Society for Music Information
Retrieval Conference (ISMIR), 2013, pp. 369–374.

[26] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Transactions
on Acoustics, Speech, and Signal Processing, vol. 14, no. 4, pp.
1462–1469, 2006.

[12] K. Yoshii, R. Tomioka, D. Mochihashi, and M. Goto, “Infinite
positive semidefinite tensor factorization for source separationof mixture signals,” in International Conference on Machine
Learning (ICML), 2013, pp. 576–584.

[27] U. Şimşekli, A. Liutkus, and A. T. Cemgil, “Alpha-stable matrix factorization,” IEEE Signal Processing Letters, vol. 22, no.
12, pp. 2289–2293, 2015.

[13] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix
factorization with the Itakura-Saito divergence: With application to music analysis,” Neural Computation, vol. 21, no. 3,
pp. 793–830, 2009.

55