2017 IEEE INTERNATIONAL WORKSHOP ON MACHINE LEARNING FOR SIGNAL PROCESSING, SEPT. 25–28, 2017, TOKYO, JAPANINFINITE PROBABILISTIC LATENT COMPONENT ANALYSIS
FOR AUDIO SOURCE SEPARATION
Kazuyoshi Yoshii1,2
1Eita Nakamura1

Katsutoshi Itoyama1

Masataka Goto3

Kyoto University 2 RIKEN 3 National Institute of Advanced Industrial Science and Technology (AIST){yoshii, enakamura, itoyama}@sap.ist.i.kyoto-u.ac.jp m.goto@aist.go.jp
ABSTRACT

1. INTRODUCTION
Statistical matrix factorization forms the basis of modern signal processing. Given a matrix X ∈ RM ×N , the typical goal of factorization is to estimate two matrices A ∈ RM ×K and B ∈ RK×N
such that X ≈ AB, where K ≪ M, N . Such low-rank representation with a limited degree of freedom enables us to extract essential
information from the original redundant data X. To perform matrix factorization, it is important to carefully examine the statisticalcharacteristics of X, A, and B and define an appropriate likelihood
function of A and B for X to be maximized. In principal component analysis (PCA), for example, X is a set of observed variables,
B is a set of the corresponding latent variables that are Gaussian distributed, A is an orthogonal transformation matrix, and A and B are
obtained by maximizing the Gaussian likelihood for X. In independent component analysis (ICA) [1], X is a set of mixture signals,
B is a set of source signals that are independently super-Gaussian(e.g., Laplacian) distributed, A is a mixing matrix, and A and B are
obtained by maximizing the Gaussian likelihood for X.
In the field of music signal processing, nonnegative matrix factorization (NMF) [2] and probabilistic latent component analysis
This work was partly supported by JST ACCEL No. JPMJAC1602 and
JSPS KAKENHI Numbers 26700020 and 16H01744, Japan.

c
978-1-5090-6341-3/17/$31.00 ⃝2017
IEEE)*+",-.-/0/12/'*3,4#0* ,5*89:)

7.;&/2%4#
1+#'2"%3</12,;".3!"#$%

Index Terms— Source separation, nonparametric Bayes, probabilistic latent component analysis, Dirichlet process, Gibbs sampling, variational Bayes.

)*+",-.-/0/12/'*3,4#0* ,5*67!

&'($)%#*+

This paper presents a statistical method of audio source separation
based on a nonparametric Bayesian extension of probabilistic latentcomponent analysis (PLCA). A major approach to audio source separation is to use nonnegative matrix factorization (NMF) that approximates the magnitude spectrum of a mixture signal at each frame
as the weighted sum of fewer source spectra. Another approachis to use PLCA that regards the magnitude spectrogram as a twodimensional histogram of “sound quanta” and classifies each quantum into one of sources. While NMF has a physically-natural interpretation, PLCA has been used successfully for music signal analysis. To enable PLCA to estimate the number of sources, we propose
Dirichlet process PLCA (DP-PLCA) and derive two kinds of learning methods based on variational Bayes and collapsed Gibbs sampling. Unlike existing learning methods for nonparametric Bayesian
NMF based on the beta or gamma processes (BP-NMF and GaPNMF), our sampling method can efficiently search for the optimal
number of sources without truncating the number of sources to beconsidered. Experimental results showed that DP-PLCA is superior
to GaP-NMF in terms of source number estimation.

!"#$%#&'(!"#$%#&'(

Fig. 1. A comparison of NMF and PLCA. NMF is a factor model
based on the sum of random variables (values of magnitude) generated from probability distributions (e.g., Poisson distributions) and
PLCA is a mixture model based on the sum of probability distributions (i.e., categorical distribution) used for generating random variables (sound quanta).
(PLCA) [3], which restrict X, A, and B to nonnegative matrices,
have been widely used for source separation of music audio signals.
Note that PLCA has often been mistakenly understood as a probabilistic version of NMF. More precisely and technically, the probabilistic models underlying NMF and PLCA are a factor model anda mixture model, respectively. This makes a clear difference in how
to model the observed data X (Fig. 1). A factor model represents
each “sample” in X as a weighted sum of all sources and can beused for solving a decomposition problem. A mixture model, on theother hand, assumes each “sample” to exclusively belong to one of
the sources and can be used for solving a clustering problem.
When NMF is used for source separation, the magnitude spectrum of each frame in an observed mixture spectrogram is regarded
as a “sample” and is approximated as the sum of source spectra.
Various probabilistic models of NMF can be formulated by specifying a probability distribution that generates the sample. Amongothers, Euclidean NMF (EU-NMF) based on the Gaussian distribution, Kullback-Leibler NMF (KL-NMF) based on the Poisson distribution [4], and Itakura-Saito NMF (IS-NMF) based on the complex
Gaussian distribution [5] are popular variants of NMF. Although ISNMF can be theoretically justified, KL-NMF has been empiricallyknown to work best for source separation. In KL-NMF, the magnitude is assumed to be Poisson distributed and must take an integer
value (can take real value in practice). This implies a deep connection between KL-NMF and PLCA.
When PLCA is used for source separation, the magnitude spectrogram in the time-frequency plane is regarded as a two-dimensional

This section reviews nonparametric Bayesian infinite mixture andfactor models of matrix factorization based on the Dirichlet, gamma,
and beta processes (DP, GaP, and BP).
2.1. Infinite mixture models
An infinite mixture model (e.g., Gaussian mixture model) for N
observations X = {xn }N
n=1 is defined by an infinite number ofcomponent distributions (e.g., Gaussian distribution) {p(x|φk )}∞k=1with mixing ratios π = {πk }∞k=1 as follows:
∞
!

k=1

πk p(xn |φk ),(1)

and covariance matrix) and Θ = {π, φ} is a set of all parameters.
Using a Dirichlet process (DP) prior DP(G0 , α) with a base measure
G0 (e.g., Gaussian-Wishart distribution) and a concentration parameter α, the generative process of xn is represented as follows:(2)

φ̂n ∼ G, xn ∼ p(xn |φ̂n ),2.1.1. Stick-breaking process∞!

πk δφk ,Ω

Fig. 2. A probability measure G drawn from a DP with a base measure G0 and a concentration parameter α.
%!

%"

!!

%#

⋯

1 # !!
!"

1 # !"!#

1 # !#

Fig. 3. A stick-breaking process for recursively generating infinitely
many weights π that sum to unity.
fragments, π, sum to unity. This process is governed by the concentration parameter α as follows (Fig. 3):
ηk ∼ Beta(1, α), πk = ηk

k−1
!

k′ =1

(1 − ηk′ ).

π ∼ SBP(α).

(5)

(6)

Using the SBP, Blei and Jordan [11] proposed a deterministic
learning method based on variational Bayes (VB) that approximatesthe complicated true posterior of the parameters Θ as a tractablefactorized distribution that can be iteratively optimized. Anotherapproximation required in practice is that the number of components considered should be truncated at a sufficiently large level because infinitely many parameters cannot be dealt with in reality. The
method is thus initialized with sufficiently many components andunnecessary components are gradually removed in each iteration toestimate the appropriate number of components, K +Another constructive interpretation of the DP is known as a Chinese
restaurant process (CRP) [12] that sequentially generates component
parameters φ̂ = {φ̂n }N
n=1 (φ̂n ∈ φ) used for generating X ={xn }Nn=1 as follows:nφ̂n+1 |φ̂1 , · · · , φ̂n ∼(3)where φ̂n ∈ φ is the parameter of a component distribution used forgenerating xn and G can be explicitly written as the sum of infinitely
many Dirac measures as follows (Fig. 2):
φk ∼ G0 , G =#"2.1.2. Chinese restaurant process

where φ = {φk }∞k=1 is a set of component parameters (e.g., meanG ∼ DP(G0 , α),)"&*+$,%&$"#'($%"$"

More simply, the SBP is often written as follows:

2. PRIOR ART

p(xn |Θ) =!"#$%&$"#'($%"!

⋯

histogram of “sound quanta,” each of which is assumed to be asample generated from one of the sources. If each time-frequency
bin is regarded as a single sample as in NMF, it cannot belong tomultiple sources. Instead, each bin is considered to include multiple sources by collecting the sound quanta generated from multiple
sources. This idea was inspired by topic models proposed for natural
language processing. While each word in a document is generated
from one of the topics, the document can be considered to include
multiple topics. PLCA was named after a basic topic model called
probabilistic latent semantic analysis (PLSA) [6].
Considering the difference between NMF and PLCA, we propose nonparametric Bayesian PLCA based on the Dirichlet process(DP) called DP-PLCA that can estimate the number of sources. The
beta or gamma process (BP or GaP), on the other hand, has beenused for formulating nonparametric Bayesian factor models such as
BP-NMF [7, 8] and GaP-NMF [9]. While in theory infinitely manysources are supposed to exist, only a limited number of them are effectively used for representing the finite observed data. We derive
two kinds of learning methods for DP-PLCA. One is a deterministic method based on variational Bayes and the other is a stochastic
method based on collapsed Gibbs sampling. Unlike existing learning
methods for BP-NMF and GaP-NMF, our sampling method can efficiently search for the optimal number of sources without truncating
the number of sources to be considered. In this paper, we examine
the capability of DP-PLCA for estimating the number of sources.

(4)

k=1

A popular constructive interpretations of the DP is a stick-breaking
process (SBP) [10] that generates the mixing ratios π by recursively
breaking off a stick with a unit length such that the lengths of the

1 !
α
δφ̂ ′ +
G0
n
α+n ′
α+n
n =1
Kn

=

1 !
α
G0 ,
nk δφk +
α+n
α+n

(7)

k=1

where Kn is the number of different components (classes) used forgenerating {xn′ }nn′ =1 and nk is the number of samples generated
from class k. More simply, we often sayφ̂ ∼ CRP(G0 , α) or Z ∼ CRP(α),(8)

where Z = {zn }N
n=1 is a set of latent variables (class indicators)
and each zn is represented as a one-hot vector. If xn is generated
from component k, i.e., φ̂n = φk , the k-th dimension of zn takes 1
and the other dimentions take 0.

A key advantage of this representation involving only φ̂ is that
it is unnecessary to deal with the infinite-dimensional vector π bymarginalizing out G from Eq. (2) and Eq. (3). Instead, we consider
at most a finite number of components φ̂ actually used for generating
a finite amount of observed data X.
Using the CRP, Neal [13] proposed a stochastic learning method
based on Gibbs sampling (GS) that is used for generating samples(values of Θ) from the complicated posterior of Θ without calculating its intractable normalizing factor. Unlike the VB method, theeffective number of component K + used for representing X can be
stochastically estimated in each iteration and finally the posterior of
K is obtained. Although the convergence is often hard to judge, in
general the GS method is more efficient (more iterations are needed,
but each iteration can be performed much faster) and more robust to
local maxima than the VB method.
2.1.3. Topic models
Topic models are an important family of mixture models originallyused for natural language processing. The most basic model is probabilistic latent semantic analysis (PLSA) [6]. Let M be the numberof different words in a dictionary. Given a set of N documents asobserved data X, PLSA aims to estimate K topics (i.e., unigramK
probabilities of M words), {{p(m|k)}Mm=1 }k=1 , and the mixing raKNtios of those topics {{p(k|n)}k=1 }n=1 in an unsupervised manner.
The probabilistic model of PLSA is often written as follows:p(n, m) =K!

p(m|k)p(k|n)p(n)!"#$%&'()"*+ /(0#'(1.#(&20 6)&,7&2$2#+/(0#'(1.#(&208!"#, %|' ( 1*!"#, %|' ( 2*=p(m|k)p(n|k)p(k),=")#&'(>$?$2$'"#$
3(0#&%'",+ &4+0&.2/+5."2#"+6,"%2(#./$+ 07$)#'&%'",89"):+1(2+:"0"2+(2#$%$'+ )&.2#
90#(,"#$+

Fig. 4. The overview of PLCA. The magnitude spectrogram is assumed to follow a two-dimensional categorical mixture distribution.
model of PLCA is given by

(9)
(10)

=

where p(n, m) is a joint probability distribution used for generating
X, i.e., a two-dimensional histogram of word counts over N documents and M words and p(n) is a probability distribution over N
documents, which is set to the empirical ratios of word counts overN documents. Eq. (9) is usually converted to an equivalent representation given by Eq. (10) for mathematical convenience. The parameters, Θ, of all the three categorical distributions p(m|k), p(n|k), andp(k) can be estimated by using the expectation-maximization (EM)algorithm such that the likelihood of Θ for X given by Eq. (10) is
maximized (maximum likelihood estimation).
To deal with new documents that are not included in X, i.e.,
to formulate a complete generative model of documents, a Bayesian
extension of PLSA called latent Dirichlet allocation (LDA) was proposed [14]. LDA is based on Eq. (9) and assumes that p(n) is a uniform distribution. Putting Dirichlet priors on p(m|k) and p(k|n),respectively, a topic distribution p(k|n′ ) for a new document n′ canbe generated.
To estimate the effective number of topics K + , Teh et al. [15]proposed nonparametric Bayesian LDA that puts a hierarchical DP(HDP) prior on a set of infinitely many topic distributions (mixing raNtios) {{p(k|n)}∞k=1 }n=1 and word distributions (component distri∞butions) {{p(m|k)}Mm=1 }k=1 . If an independent DP prior is put foreach document n, a document-specific topic distribution {p(k|n)}∞k=1∞
and word distributions {{p(m|k, n)}Mm=1 }k=1 are generated. To
avoid this, a higher-level DP prior is put on document-wise DP priors∞
for sharing {{p(m|k)}Mm=1 }k=1 over all documents.
2.1.4. Probabilistic latent component analysis
Probabilistic latent component analysis (PLCA) [3] is a generalization of PLSA for D-dimensional data (D ≥ 1). The probabilistic! "# ! $#9"):+1(2+:"0"+7'&1"1(*(#;p(x1 , · · · , xD ) =k=1<(-(2% !&,7&2$2#'"#(& /(0#'(1.#(&2!"#$%&'()"*+ ,(-#.'$+ /(0#'(1.#(&2

k=1
K
!

!"#, %|' ( 3*

K
!

k=1
K
!

k=1

p(k)p(x1 , · · · , xD |k)
p(k)

D
"

d=1

p(xd |k),(11)

where xd is a random variable of any type in dimension d (1 ≤ d ≤D) and a sample generated from Eq. (11) is represented by a tuple(x1 , · · · , xD ). When d = 2 and x1 and x2 are discrete randomvariables, PLCA reduces to PLSA (x1 and x2 correspond to n andm, respectively, as shown in Fig. 4).
PLCA have successfully been used for audio source separation
by regarding n and m as time and frequency indices, respectively[3, 16, 17]. A pair (n, m) denotes the time-frequency position of a“sound quantum” generated from Eq. (11) or Eq. (10). Given themagnitude spectrogram of a music signal, i.e., a two-dimensional
histogram of sound quanta, in the time-frequency plane with N frames
and M frequency bins, the frequency distributions of K sources,
K
K
N{{p(m|k)}Mm=1 }k=1 , and their mixing ratios, {{p(k|n)}k=1 }n=1 ,can be estimated by the EM algorithm as in PLSA. The frequency
distributions can be trained from isolated musical instrument sounds
in advance [16,17]. Although HDP-PLCA could be formulated as in
HDP-LDA, in this paper we show that a simpler and more easy-toimplement extension called DP-PLCA is feasible by using Eq. (11)(Eq. (10) instead of Eq. (9)) and a DP prior for efficient Bayesianinference based on VB or GS.
2.2. Infinite factor models
A finite factor model (e.g., KL-NMF) for N observations X ={xn }Nn=1 is defined by a particular probability distribution (e.g.,
Poisson distribution) that is specified by the sum of basis parametersK(factors) φ = {φk }K
k=1 with their local weights ωn = {ωnk }k=1
specific to sample n as follows:
# $ K
%
$!
$
p(xn |Θ) = p xn $
ωnk φk .
(12)
$
k=1

The gamma process (GaP) [9] or beta process (BP) [7, 8] canbe used for formulating an infinite factor model by taking the limitof Eq. (12) when K goes to infinity. A key feature of infinite factor
models is that although in theory infinitely many factors are assumed
to exist, at most a finite number of factors are effectively used for
representing X. The BP was also used for estimating the number ofbasis patterns (atoms) in dictionary learning from images [18, 19].
2.2.1. Gamma process
Introducing a K-dimensional nonnegative vector π = {πk }K
k=1 into
Eq. (12) and taking the infinite limit when K → ∞, an infinite factor
model based on the GaP is given by
# $ ∞
%
$!
$
p(xn |Θ) = p xn $πk ωnk φk ,
(13)
$
k=1

where πk is a global weight of factor k and π does not need to sum
to unity. To make π sparse, i.e., to make fewer factors effective,a GaP prior GaP(G0 , α) with a base measure G0 over a space Ω(G0 (Ω) = γ is the total mass of G0 over Ω) and a concentration
parameter α can be used as follows:G ∼ GaP(G0 , α),(14)where G consists of infinitely many atom as follows:ωk , φk ∼∞!
1
G0 , G =πk δωk ,φk ,γ(15)where G and H over Ω can be explicitly written as
ωk , φk ∼
H=

∞
!

1
G0 ,
γ

πk δωk ,φk , Gn =

k=1

(19)
∞
!

znk δωk ,φk ,

(20)

k=1

where πk ∈ [0, 1] indicates the probability of activating factor k
and znk ∼ Ber(πk ). Since the infinite-dimensional vector π ={πk }∞
k=1 is extremely sparse, only a limited number of different factors are activated in Z.
Several kinds of SBPs have been proposed for explicitly representing π. Teh et al. [21], for example, proposed a SBP similar
to Eq. (5) and derived a GS method that can estimate Θ without
truncated approximation by using slice sampling. Gupta et al. [22]proposed BP-NMF based on the same SBP representation. Paisley
et al. [23, 24] proposed another SBP as the infinite limit of a finite
model and derived a VB method. Another simpler representationof the BP is weak-limit approximation that independently puts an
extremely sparse beta prior on each πk as follows:& αγ &
γ '',α 1 −,(21)
πk ∼ Beta
K
K
where K is set to a sufficiently large number as in the GaP. Liang etal. [7,8] proposed BP-NMF based on this representation and derived
VB and GS methods with truncated approximation. A CRP-like representation of the BP is known as an Indian buffet process (IBP) [25].

k=1

While a SBP for generating π was recently derived [20], a CRP
representation has not been proposed yet. Although a VB method
based on the SBP representation was proposed for estimating theeffective number of factors K + in matrix factorization [20], it has
scarcely been used in practice due to its complexity. Another simpler representation of the GaP is weak-limit approximation that independently puts an extremely sparse gamma prior on each πk as
follows:& αγ ',α .
(16)
πk ∼ Gamma
K
Since the GaP is obtained as the limit of Eq. (16) when K → ∞,K is set to a sufficiently large number in practice in exchange for
unnecessarily increase of computational cost.
A VB method based on Eq. (16) was proposed for GaP-NMF
and applied to music signal analysis [9]. It is, however, often difficult
to determine a threshold for gradually removing ineffective factors
to estimate K + . In addition, the results obtained by the weak-limit
approximation tend to be sensitive to the truncation level.
2.2.2. Beta process
Introducing a binary matrix Z = {znk }N,K
n=1,k=1 into Eq. (12) andtaking the infinite limit when K → ∞, an infinite factor model
based on the BP is given by
# $ ∞
%
$!
$
p(xn |Θ) = p xn $znk ωnk φk ,
(17)
$
k=1

where znk is a latent binary variable that indicates the presence orabsence of factor k in sample n. To make Z sparse, a BP prior
BP(G0 , α) with a base measure G0 over a space Ω (G0 (Ω) = γ)and a concentration parameter α is used with a Bernoulli process
BeP(H) with a base measure H as follows:H ∼ BP(G0 , α), Gn ∼ BeP(H),(18)

3. DIRICHLET PROCESS PLCA
This section describes a nonparametric Bayesian extension of PLCA
based on the Dirichlet process (DP) for audio source separation andthen derives two kinds of learning methods based on VB and collapsed GS.
3.1. Data preparation
To use PLCA for audio source separation, the magnitude spectrogram of a mixture signal is regarded as a two-dimensional histogramof “sound quanta” in the time-frequency plane with N frames andM frequecy bins, as explained in Section 2.1.4. We redefine X ={xi }Ii=1 as a set of observed variables and let Z = {zi }Ii=1 be a setof the corresponding latent variables, where I is the total number ofsound quanta, i.e., the sum of the values of magnitude over the timefrequency plane. Since each quantum i is located at a time-frequency
bin, xi is represented as an N M -dimensional one-hot vector thattakes 1 in a dimension corresponding to the time-frequency location and 0 in the other dimensions. In addition, each quantum i isassumed to be generated from one of K sources (K → ∞), zi isrepresented as a K-dimensional one-hot vector.
To make the observed data X, we need to quantize the valueof magnitude at each time-frequency bin according to an appropriate
resolution. In practice, PLCA is found to work by scaling and rounding the magnitude spectrogram such that the average of magnitude
is around 1. If the average of magnitude is set to a larger value, thenumber of observations I takes a larger value, which affects the posterior uncertanty of the parameters of PLCA in Bayesian estimation.
This is more problematic for nonparametric Bayesian PLCA because
the effective number of sources, K + , is considered to increase logarithmically according to the increase of observed data. To solve
this problem, we could use a method of optimizing the resolution ofmagnitude quantization that was originally proposed for KL-NMF
based on the discrete Poisson distribution [26].

3.2. Model formulation

Finally, q(φ) and q(θ) are given by

The probabilistic model of PLCA for two-dimensional data is given
by Eq. (10), where p(k) represents the mixing ratios of K sources,
and p(n|k) and p(m|k) represent the time and frequency distributions of source k, respectively. The parameters of these categoricalN
distributions are given by π = {πk }K
k=1 , φk = {φkn }n=1 , andθk = {θkm }M.
The
probabilistic
generative
model
of
Z
and X
m=1
(the complete likelihood function of π, φ, and θ for Z and X) is
thus defined as follows:p(Z|π) =

K
I "
"

z

(22)

πkik ,q(φk ) = Dir(λk ),q(θk ) = Dir(ωk ),(30)

where λk and ωk are given by!
!
λkn = β +xinm E[zik ], ωkm = γ +xinm E[zik ]. (31)imin

As in GaP-NMF, K is initialized as a sufficiently large number anduneffective sources are gradually removed in each iteration.
3.4. Collapsed Gibbs sampling

To let K go to infinity, we use a DP prior. There are two major
constructions of the DP (see Sections 2.1.1 and 2.1.2). Using a SBP
and conjugate Dirichlet priors, a nonparametric Baysian model canbe formulated as follows:

Formulating a Bayesian model given by Eqs. (22), (23), and (25)and marginalizing out φ and θ by leveraging the conjugacy between
the Dirichlet and categorical distributions, we aim to draw samples
from the posterior distribution p(Z|X) because at most N different
classes appear in Z (the number of different classes K + is usuallymuch fewer than N ). This enables us to avoid dealing with infinitely
many parameters on computers. We use collapsed Gibbs sampling
for updating each zi in a random order according to the following
conditional distribution:π ∼ SBP(α), φk ∼ Dir(β), θk ∼ Dir(γ),

p(zi |Z¬i , X) ∝ p(zi , xi |Z¬i , X¬i ),i=1 k=1

p(X|Z, φ, θ) =

I "
N
"

M
"

K
"

(φkn θkm )xinm zik .

(23)

i=1 n=1 m=1 k=1(24)

where α, β, and γ are hyperparameters. The SBP representation isconvenient for deriving the VB method.
Using a CRP, i.e., marginalizing π out analytically, we obtain
another representation as follows:
Z ∼ CRP(α), φk ∼ Dir(β), θk ∼ Dir(γ).

(25)

Note that φ and θ can also be marginalized out for efficient Bayesianinference (only Z is considered). The CRP representation is convenient when the (collapsed) GS method is used.
3.3. Variational Bayes
Using a Bayesian model specified by Eqs. (22), (23), and (24), we
aim to calculate the posterior distribution p(Z, η, φ, θ|X), where ηis considered instead of π based on Eq. (5). Since the true posterior
is analytically intractable, it is approximated as a factorizable variational posterior distribution q(Z, η, φ, θ) = q(Z)q(η)q(φ)q(θ)
such that the KL divergence from q(Z, η, φ, θ) to p(Z, η, φ, θ|X)
is minimized, i.e., the lower bound of the log-evidence p(X) is maximized. In each iteration, q(Z) is updated as follows:q(Z) ∝ exp(Eq(η,φ,θ) [log p(X, Z, η, φ, θ)]).

(26)

Similarly, q(η), q(φ), and q(θ) can also be updated alternately.
Each posterior is found to be the same type of the prior distribution
because of the conjugacy.
First, q(Z) is given by a categorical distribution as follows:q(zi ) = Categorical(ζi ),where ζik = !K

ρik

k′ =1

and ρik is given byρik′log ρik = E[log ηk ] +
+

!

(27)

k−1
!

k′ =1

xinm E[log φkn ] +nmE[log(1 − ηk′ )]!

xinm E[log θkm ].

(28)

nm

Then, q(η) is given by
⎛
⎞
K
!
! !
E[zik ], α +E[zik′ ]⎠ ,q(ηk ) = Beta ⎝1 +iik′ =k+1

(29)

(32)

where X¬x indicates a set of parameters X except for x. When zi
is updated, K + might be decremented, incremented, or unchanged.
More specifically, the probability that zi is generated from “existing” source k (zik = 1) and xi is then located at time n and frequency bin m (xinm = 1) is given by,′i′ ̸=i zi kp(zik = 1, xinm = 1|Z¬i , X¬i ) ∝I −1+α, ,, ,′′xi′ nm zi′ k + γ′
mi′ ̸=i xi nm zi k + βn,, i ̸=i. (33)′ k + βN
′
z
z
′
′
i
i ̸=i
i ̸=i i k + γM
On the other hand, the probability that zi is generated from a newsource knew (ziknew = 1) and xi is located at time n and frequency
bin m (xinm = 1) is given byp(ziknew = 1, xinm = 1|Z¬i , X¬i ) ∝α1.
I − 1 + α MN

(34)

4. EVALUATIONThis section reports a comparative experiment evaluatining the capability of DP-PLCA in estimation of the number of sources K + .
4.1. Experimental conditions
We used three mixture signals each of which was synthesized using
the sounds of piano (011PFNOM), electric guitar (131EGLPM), or
clarinet (311CLNOM) recorded in the RWC Music Database: Musical Instrument Sound [27]. Each signal (14 s) was made by concatenating seven 2-s isolated or mixture sounds (C4, E4, G4, C4+E4,
C4+G4, E4+G4, and C4+E4+G4). The sampling rate was 16 kHz.
Each mixture signal was expected to be separated into three sources
corresponding to C4, E4, and G4 (K + = 3). The short-time Fouriertransform (STFT) with a Gaussian window was performed with awindow length of 512 pts and a shifting interval of 160 pts. The
magnitude spectrogram of each signal was represented as a nonnegative matrix with N = 1400 and M = 257 and scaled such that the
average value of magnitude over the time-frequency plane, µ, wasequal to 1 or 10. We tested the VB and GS methods for learning the
DP-PLCA model. For comparison, we tested GaP-NMF based on
the Poisson distribution (GaP extension of KL-NMF, Section 2.2.1)

! ! %"%"(5'6789! " 1%":;9<=5>'?89! " 1

2@'?56A>89! " 1

%!

%!

%!

$"

$"

$"

$!

$!

$!

#"

#"

#"

#!

#!

#!

"

"

"

!

"!

!

#!!

%"

"!

(5'6789! " 10

!

#!!

%"
:;9<=5>'?89! " 10
%!

%!

$"

$"

$"

$!

$!

#"

#"
#!

"

"
!

"!

#!!

!

#!!

2@'?56A>89! " 10

%!

#!

"!

%"

$!
#"
0()(123-2&4&'()*+,-./
0()(123-./
"!

#!
"
#!!

!

"!

#!!

!"#$%"&'(

Fig. 5. Experimental results of K + estimation.
based on the VB method. The truncation level of the VB methods
was set to 30 and the GS method started from 30 sources. The performance of each method was evaluated in terms of accuracy of K +estimation. When K + = 3, these methods achieved comaprable
source separation performance.
4.2. Experimental results
Fig. 5 shows the experimental results. When µ = 1 (i.e., the total
number of sound quanta was I = N M ), all methods converged toa reasonable solution around K + = 4, where an extra source represents a noise spectrum (atack sound). When µ = 10, all methods
tended to overestimate K + . Although this is a natural behavior of
nonparametric Bayesian models, the essential complexity of X remains the same. DP-PLCA based on the VB method was found to beless sensitive to the apparent data size. The computational cost of the
VB method for DP-PLCA was much smaller than that for GaP-NMF.
Note that the computational costs of the VB methods are independent from the data size I, that of the GS method is linearly increased.
One solution is to use non-collapsed GS that can adaptively truncate
π without any approximation by using slice sampling. This enables
parall computation as in the VB methods.
5. CONCLUSION
This paper presented a nonparametric Bayesian extension of PLCA
called DP-PLCA to estimate the number of sources in audio source
separation. One of the major contributions of this paper is to clarify the essential difference between two major matrix factorization
techniques, NMF (factor model) and PLCA (mixture model), which
are often mistakenly considered to be always identical, in terms ofprobabilistic modeling. This is a reason why the DP is used with
PLCA while the GaP or BP is used with NMF. We derived the learning methods based on VB and GS and confirmed that these methods
outperformed GaP-NMF in terms of the capability of autonomousmodel complexity control.
6. REFERENCES[1] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component Analysis, John Wiley & Sons, 2004.
[2] P. Smaragdis et al., “Static and dynamic source separation using nonnegative factorizations: A unified view,” IEEE Sig.
Proc. Mag., vol. 31, no. 3, pp. 66–75, 2014.

[3] M. Shashanka, B. Raj, and P. Smaragdis, “Probabilistic latentvariable models as nonnegative factorizations,” Computational
Intelligence and Neuroscience, vol. 2008, pp. 1–8, 2008.
[4] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization for polyphonic music transcription,” WASPAA, 2003.
[5] C. Févotte et al., “Nonnegative matrix factorization with the
Itakura-Saito divergence: With application to music analysis,”
Neural Computation, vol. 21, no. 3, pp. 793–830, 2009.
[6] T. Hofmann, “Learning the similarity of documents: An
information-geometric approach to document retrieval and categorization,” NIPS, 2000, pp. 914–920.
[7] D. Liang, M. Hoffman, and D. Ellis, “Beta process sparse nonnegative matrix factorization for music,” ISMIR, 2013.
[8] D. Liang and M. D. Hoffman, “Beta process non-negative matrix factorization with stochastic structured mean-field variational inference,” NIPS Workshop, 2010.
[9] M. Hoffman, D. Blei, and P. Cook, “Bayesian nonparametricmatrix factorization for recorded music,” ICML, 2010.
[10] J. Sethuraman, “A constructive definition of Dirichlet priors,”
Statistica Sinica, vol. 4, pp. 639–650, 1994.
[11] D. M. Blei and M. I. Jordan, “Variational inference for Dirichlet process mixtures,” Bayesian Analysis, vol. 1, no. 1, pp.
121–144, 2006.
[12] T. Ferguson, “Bayesian analysis of some nonparametric problems,” Annals of Statistics, vol. 1, no. 2, pp. 209–230, 1973.
[13] R. M. Neal, “Markov chain sampling methods for Dirichlet process mixture models,” Journal of Computational andGraphical Statistics, vol. 9, no. 2, pp. 249–265, 2000.
[14] D. M. Blei, A. Y. Ng, and M. I. Jordan, “Latent Dirichlet allocation,” JMLR, vol. 3, pp. 993–1022, 2003.
[15] Y. W. Teh et al., “Hierarchical Dirichlet processes,” JASA, vol.
101, pp. 1566–1581, 2006.
[16] B. Fuentes, R. Badeau, and G. Richard, “Harmonic adaptive latent component analysis of audio and application to music transcription,” IEEE TASLP, vol. 21, no. 9, pp. 1854–1866, 2013.
[17] E. Benetos and S. Dixon, “Multiple-instrument polyphonicmusic transcription using a temporally constrained shiftinvariant model,” JASA, vol. 133, no. 3, pp. 1727–1741, 2013.
[18] M. Zhou et al., “Nonparametric Bayesian dictionary learning
for analysis of noisy and incomplete images,” IEEE Trans. on
Image Proc., vol. 21, no. 1, pp. 130–144, 2012.
[19] H. P. Dang and P. Chainais, “A Bayesian non parametric approach to learn dictionaries with adapted numbers of atoms,”
MLSP, 2015, pp. 1–6.
[20] A. Roychowdhury and B. Kulis, “Gamma processes, stickbreaking, and variational inference,” AISTATS, 2015, pp. 800–808.
[21] Y. W. Teh et al., “Stick-breaking construction for the Indian
buffet process,” AISTATS, 2007, pp. 556–563.
[22] S. K. Gupta et al., “A nonparametric Bayesian Poisson gamma
model for count data,” ICPR, 2012, pp. 1815–1818.
[23] J. Paisley et al., “A stick-breaking construction of the beta
process,” ICML, 2010, pp. 847–854.
[24] J. Paisley, L. Carin, and D. Blei, “Variational inference for
stick-breaking beta process priors,” ICML, 2011, pp. 889–896.
[25] T. L. Griffiths and Z. Ghahramani, “Infinite latent feature models and the Indian buffet process,” NIPS, 2006, pp. 475–482.
[26] M. D. Hoffman, “Poisson-uniform nonnegative matrix factorization,” ICASSP, 2012, pp. 5361–5364.
[27] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka, “RWC
music database: Popular, classical, and jazz music database,”
ISMIR, 2002, pp. 287–288.