This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

1

Dirichlet Process Mixture of Mixtures Modelfor Unsupervised Subword Modeling
Michael Heck, Member, IEEE, Sakriani Sakti, Member, IEEE, and Satoshi Nakamura, Fellow, IEEE

Abstract—We develop a parallelizable Markov chain Monte
Carlo sampler for a Dirichlet process mixture of mixtures model(DPMoMM). Our sampler jointly infers a codebook and clusters.
The codebook is a global collection of components. Clusters are
mixtures, defined over the codebook. We combine a non-ergodic
Gibbs sampler with two layers of split and merge samplers on
codebook and mixture level to form a valid ergodic chain. We
design an additional switch sampler for components that supports
convergence in our experimental results. In the use case of
unsupervised subword modeling, we show that our method infers
complex classes from real speech feature vectors that consistently
show higher quality on several evaluation metrics. At the sametime we infer fewer classes that represent subword units moreconsistently and show longer durations, compared to a standard
DPMM sampler.
Index Terms—acoustic unit discovery, Bayesian nonparametrics, Dirichlet process, Gibbs sampling, mixture of mixtures,
unsupervised subword modeling

I. I NTRODUCTIOND

IRICHLET process mixture models (DPMMs) [1], [2]are firmly established in pattern recognition and machinelearning. Also known as infinite mixture models [3], they elegantly extend finite mixture models by the aspect of automatic
model selection. This property makes them a popular tool forsolving clustering tasks that are challenging with regards toestimating model complexity a priori. Several extensions to
the original concept have been introduced over time, mostnotably hierarchical models [4], [5] and Dirichlet processes(DPs) with dependencies [6], [7].
DPMMs with Gaussian components gained increased interest in the field of low resource automatic speech processing,particularly as method for tackling the task of unsupervised
subword modeling. The task is to infer acoustic units from raw
audio data that are suitable to reliably represent human speech,
i.e., that show low discriminability errors. DPMM samplers
were used for subword model inference in an array of works
related to the zero resource speech challenge [8]–[11]. The
idea is that each Gaussian in a mixture model that was inferredfrom speech data is considered a separate acoustic class. Previous work by the authors [12]–[14] improved unsupervised
subword modeling via DPMM sampling by unsupervisedlytransforming the sampler’s input.
A major impediment for producing better subword models
however is the simplicity of the inferred model. It is a longstanding modeling assumption that speech observations, i.e.,feature vectors that belong to specific sound categories, are
M. Heck, S. Sakti and S. Nakamura are with Nara Institute of Science andTechnology.
Manuscript received xxx xx, 201x; revised xxx xx, 201x.

multimodally distributed [15]. In practice, Gaussian mixture
models (GMMs) are well established to model acoustic units
such as phones [16], [17], and methods such as continuous
hidden Markov models (HMMs) make use of state-dependent
GMMs as multimodal distributions to model emission probabilities of speech observations [15], [18]. It seems therefore
an oversimplification to assume single unimodal distributions
to be a good model representation for individual sounds. This
assumption limits the inferred units to represent generally very
short stationary sound phenomena [19].
DPMMs work very well when the clusters in a data set
are unimodally distributed. But problems arise when clusters
follow more complex, e.g., multimodal, distributions. In suchcases, a model that fits unimodal distributions (e.g., single
Gaussians) to clusters tends to over-fragment the feature space
and to suggest more clusters than actually present. In otherwords, real clusters tend to be represented by multiple components, i.e., “sub-clusters”. However, without dependency
modeling in DPMMs, inferred “sub-clusters” are considered
independent and the relations between them are lost. The
consequence is that the inferred clusters do not reflect theactual structure in the data. To more accurately approximate
multimodally distributed clusters, a model that assigns multiple mixture components to each cluster would be required [20].
We consider unsupervised subword modeling to be sucha problem where the ability to infer multimodal clusters
from a data set can provide models that represent the real
underlying data distribution more accurately. The expectationis that a mixture of clusters, where each cluster is a mixture
itself, should be a better representation of acoustic units withfavorable characteristics. Specifically, one would expect the
number of inferred classes to be lowered, the overall model
size therefore be reduced. At the same time, average durations
when classifying sequences of sounds should be longer. Sound
representations should also be more robust across data ofdifferent speakers and show higher discriminability due tomore natural modeling.
We develop a sampler for a Dirichlet process mixture of
mixtures model (DPMoMM) to overcome the limitations of
DPMMs and to enable the inference of a mixture of multimodal clusters. In our proposed DPMoMM, each cluster is a
mixture of components, and the collection of clusters forms a
global mixture of mixtures. Throughout this paper, we will use
the following terms to describe our model. Each mixture in
the global mixture of mixtures is called a cluster. Each cluster
is a mixture of cluster components. Cluster components canbe shared across clusters, which is why they exist on a global
level. We borrow a term from automatic speech recognitionand call the global collection of components that can be part of

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

Mixture of Mixtures

Parameters:

Clusters

Codebook

Fig. 1. Illustration of the mixture of mixtures model. The codebook is a global
collection of components. Clusters are mixtures, defined over the codebook.

a cluster a codebook. When we speak of codebook components
or simply components, we refer to the global components thatmake up the codebook. A codebook component is a clustercomponent (i.e., a member) in at least one cluster, and eachcluster component is identical with exactly one codebook
component. The difference is that the same component has
different weights in different clusters (if it belongs to morethan one). Besides that, each codebook component has a global
weight. An intuitive illustration of this model is provided inFig. 1. The detailed description of the model and its parameters
is given in Section III.
We build on the idea of Chang et al. [21] and develop a split,
merge and switch sampler with the following characteristics:(1) our sampler jointly infers a global codebook of components, and clusters which are mixtures that are defined overthis codebook; (2) split and merge moves modify codebook
components; (3) split, merge and switch moves modify clustercomponents; (4) all sampling steps can be parallelized acrossclusters and components; (5) the jointly inferred codebook andthe mixture of mixtures provide two alternatives to model thesame underlying data.
We demonstrate in unsupervised subword modeling use caseexperiments on real speech data that our DPMoMM sampler
is superior to a DPMM in terms of inferring models thatare better representations of the underlying data structure.
Specifically, with our method we infer mixtures of Gaussian
mixtures from real speech observations that consistently show
higher quality on several subword model evaluation metrics.
For that, we extract frame-wise speech feature vectors from adata set and use our proposed sampler to cluster these speechobservations into classes. In a standard DPMM, each class isrepresented by a single Gaussian, whereas with our proposed
DPMoMM, each class is represented with its own GMM. Thisway we infer fewer clusters that represent subword units moreconsistently across speakers and that show longer average
durations. We also show that an additional switch sampler
supports the convergence of the algorithm.

2

II. R ELATED W ORK
The hierarchical Dirichlet process (HDP) is a well-known
method for sampling mixture models that employ a hierarchy [4], [5]. The HDP can be used to infer topics that areshared between multiple documents, i.e., groups of data. An
analogy between HDP and DPMoMM can be drawn to thefollowing extent. A document in the HDP is a mixture of
topics, just as a cluster in the DPMoMM is a mixture ofcomponents. Documents in the HDP share topics from a global
set, just as clusters in the DPMoMM can share components
from the codebook. The similarities between the models end
at this point, however. The HDP differs greatly from the
DPMoMM by assuming that a particular grouping of data
into a finite set of documents is known a priori. Topics are
assumed to be shared across documents, and each document isassumed to have its own particular distribution of these topics.
Topic mixtures that describe individual documents can heavily
overlap each other if they have many topics in common. In
contrast, our proposed model does not assume a pre-definedgrouping of data into document-like clusters. Instead, theDPMoMM sampler infers an unknown number of clusters
within ungrouped data, comprised of an unknown number ofcomponents each. Our method infers a mixture of mixtures,i.e., a group of clusters which itself forms a mixture modelwith explicit mixture weights. In the DPMoMM, clusters aregroups of neighboring components, and clusters do not tend
to occupy the same regions in the feature space.
More related to our DPMoMM is the infinite mixture of
infinite Gaussian mixtures (I2 GMM) by Yerebakan et al. [20].
The I2 GMM is a generative model that represents each cluster
within a data set with its own mixture model. Here, a top
layer DP defines meta-clusters, and lower layer DPs model
the cluster data as a mixture of components. The top layer
generates cluster parameters according to a base distribution
H and defines the number and local expansion of clusters.
The cluster parameters in turn define base distributions Hk
for the lower layer which control the number and local
expansion of cluster components. Covariance matrices areshared across components within the same cluster, leaving
components to only differ in their means. In contrast to theI2 GMM, our proposed model does not define a prior overcluster appearances in form of meta-clusters. Instead, DPMoMM clusters can take on any structure that the data might
inform. DPMoMM components also have their own covariancematrix each, which allows more natural approximations to
the data. Further, unlike the I2 GMM, the DPMoMM supports
the sharing of components across clusters, which enables thesharing of statistical strength [4].
Dependent Dirichlet processes are suitable to capture timedependencies between clusters or samples [22]–[24]. Temporaldependencies might in certain cases be reflected by locality
in the feature space. With the DPMoMM, we propose a
new kind of model sampler that explicitly infers a mixtureof multimodal distributions to handle dependencies that arereflected by locality in the feature space.
Split and merge samplers are thoroughly discussed in anarray of publications [25]–[27]. The DPMM sub-cluster al-

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

gorithm of Chang et al. [21] addresses several issues that
previous approaches coped with. Their sampler combines a
non-ergodic restricted Gibbs sampler and split and merge
samplers into a valid Markov chain. The Gibbs sampler isrestricted to non-empty clusters. Splits are proposed from subclusters that are learned jointly by deferred sampling. Moves
are proposed with a Metropolis-Hastings (MH) algorithm. As
instantiated-weight (IW) sampler, cluster weights are explicitly
represented, as opposed to collapsed-weight (CW) samplers.
Like in [28], [29], no finite approximations are used for the
Dirichlet process, contrary to [29], [30]. The authors see theadvantage of IW samplers in the possibility to parallelize
across data points (which they refer to “inter-cluster parallelizable”) and propose to use global split and merge moves tocounter convergence issues. Inspired by these works, Chang etal. [21] propose moves that rely on jointly learned sub-clusters
to reduce computational overhead during the MH steps.
This algorithm was used with success for unsupervised
subword modeling in the scope of the zero resource challenge [8]. Chen et al. [10] inferred a DPMM from raw
speech, with Gaussians as components, and used the Gaussian posteriorgrams extracted after sampling as new speech
representation. In our own work, we extended this approachby developing a framework that unsupervisedly learns featuretransformations from inferred classes [12]–[14]. We showedthat these transformations in turn improve the input to the
DPMM sampler so that even better classes can be inferred,
according to a sound class discriminability measure. We
further demonstrated that the inferred classes can be used tomodel sounds for speech recognition purposes [19]. However,
we found there is need for a method to find more complex
classes that are generalizing across speakers and that coverconsistent sequences with longer durations. We show in thiswork how we enhance the DPMM to be a DPMoMM thatjointly learns components and mixtures of components andhow we successfully use our model to improve unsupervised
subword modeling.
Discriminative (as alternative to generative) non-parametricmodels such as infinite (structured) support vector machines(i(S)SVMs) [31], [32] have also been successfully appliedin the ASR domain to dynamically model speech concepts.
The idea of the iSVM and iSSVM is to divide the feature
space into regions to be handled by a mixture of experts,
i.e., specialized sub-models, where the number of experts isinferred from the data and the mixture underlies a DP prior.
However, the number of actual concepts to be modeled is
known beforehand, and the SVM training is supervised andrelies on labels. In contrast, our DPMoMM infers concepts
from data only without prior knowledge of any sort.
III. DP M IXTURE OF M IXTURES M ODELIn this section, we develop our DP mixture of mixturesmodel (DPMoMM). Definitions and mathematical expressions
are kept general and are not restricted to a specific type of
mixture components. Within the scope of this work, we use
our sampling algorithm to infer mixtures of Gaussian mixturesin the use case of unsupervised subword modeling, for which3

we originally developed this method (see Section VIII). We
begin by reviewing the standard DPMM and the augmented
DPMM of Chang et al. [21].
A. Graphical Representation
Fig. 2 is a representation of the general DPMM in plate
notation. xi is an observed data point i out of N data points,
and zi is the corresponding discrete label for that data point. πdenotes the theoretically infinite dimensional vector of mixture
weights. α is commonly referred to as the concentration parameter for the Dirichlet process, which governs the likelihoodfor new classes to be generated during sampling, and λ is thehyper-parameter for the Dirichlet process base measure. θkdenotes the parameters of cluster k, e.g., mean and covariancein the case of Gaussians. The generative process of the DPMM
is expressed as follows:xi ∼ p(xi |θzi ),zi ∼ Discrete(π),θk ∼ H(λ),π ∼ GEM(α),

(1)
(2)

where GEM(·) denotes the stick-breaking process, and H isthe DP base measure. The generative story of a data point
xi is this. A discrete cluster label zi is sampled from the setof all possible clusters, which are distributed according to theweights in π. Given the cluster label, xi is drawn from thecluster with parameters θzi .
Fig. 3 is Chang et al.’s [21] augmented DPMM using
auxiliary variables. Each regular cluster is augmented withtwo explicit sub-clusters, denoted as l for “left” and r for“right”. The goal is to design a model that is tailored toward
splitting clusters. By picking suitable distributions for these
sub-clusters, then they can provide good split proposals fortheir regular parent cluster. Each data point is assigned to either
the “left” or “right” sub-cluster with a sub-cluster label z̄i ∈{l, r}. The naming convention implies that the sub-clusters
are designed towards separating the data points into distinctgroups within the parent cluster. Sub-clusters have their own
weights π̄k = {π̄kl , π̄kr } and parameters θ̄k = {θ̄kl , θ̄kr }. It is
important to note that in this auxiliary space the data points xi
generate the labels z̄i , in contrast to the regular space where
zi generate xi .
Our proposed DPMoMM is depicted in Fig. 4. As before,
xi is an observed data point i out of N data points that belong
to a cluster k. β governs the global mixture proportions, andπ is the vector of weights for clusters, sampled according
to a stick-breaking process. zi denotes the cluster assignment
label for the corresponding data point. In our model, clusters
are composed of components, which are represented by thefollowing variables. c̃ki is the label of the cluster component
conditioned on cluster k that the corresponding data point is
assigned to. π̃k are the cluster component weights, governedby β and conditioned on cluster k. θc are the parameters for thecomponent that generated xi . All cluster components are defined globally in the codebook. The codebook is the weighted
collection of all existing components in the DPMoMM, i.e.,
it is a global mixture of components. Given the codebook,sampling operations can be performed on component level,
independent of their respective cluster memberships. ċi assigns

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

4

Fig. 2. Standard DPMM.

Fig. 3. Augmented DPMM by Chang et al. [21]. Auxiliary variables aredenoted by dotted circles.

Fig. 4. Proposed DP mixture of mixtures model (DPMoMM). Auxiliary
variables are denoted by dotted circles.

a data point i to a codebook component and is derived
from the cluster component label c̃ki . Codebook components
have global weights π̇, governed by a separate concentration
parameter α. λ is again the hyper-parameter for the DP basemeasure. The generative process of the DPMoMM is formally
expressed as follows:

Algorithm 1 DPMoMM sampling algorithm
Randomly initialize Kinit clusters with 1 component eachwhile stop criterion not met do
Propose cluster merges and splits. Sec. VI
Propose cluster component switches. Sec. VI
Propose component merges and splits. Sec. V
for all clusters with split components do
Duplicate and update. Sec. V
end forSample parameters and labels. Sec. IV
end whilexi ∼ p(xi |θc̃zi i ),

θc ∼ H(λ),(3)

zi ∼ Discrete(π),
c̃ki ∼ Discrete(π̃k ),ċi = c̃zi i ,π ∼ GEM(β),π̃k ∼ GEM(β),π̇ ∼ GEM(α).

(4)
(5)
(6)

The generative story for any data point xi in a DPMoMM
is as follows. A cluster label zi is sampled from the set of
all possible clusters, which are distributed according to theweights in π. Then, a component label c̃zi i is sampled fromthe set of components that belong to the cluster with label zi ,which are distributed according to the weights in π̃zi . Given
the cluster and cluster component labels, xi is drawn fromthe cluster component with parameters θc̃zi i . The codebook
is a by-product of this generative process. The codebook
component labels are copies of the cluster component labels
according to (6), and the weights of the codebook components
are conditioned on these labels. Fig. 1 is an illustration of thehierarchy within a DPMoMM. The sampling of the full model
is described in detail in Section IV.
An intuition of what the DPMoMM represents can be given
as follows. Assuming the tackled task is topic clustering, one
can view DPMoMM clusters as groups of closely related
topics (modeled by cluster components). If for instance there
are three components modeling the topics “cars”, “trucks” and“motorbikes”, then a cluster that contains these components
would model the meta-topic “personal vehicles”.
Similar to the augmented DPMM, we define an auxiliary
space to enable a split and merge sampling approach. Eachcomponent in the DPMoMM is augmented with two subcomponents, parametrized by θ̄c = {θ̄cl , θ̄cr }, to provide good
split candidates for a component split move proposal. c̄ki ∈{l, r} is the label that assigns the corresponding data point to alr}
sub-component of c̃ki within cluster k, and π̄kc = {π̄kc, π̄kcdenotes the weights for sub-components of component cwithin cluster k. On codebook level, sub-component labelsc̄˙i are derived from the cluster sub-component labels c̄ki , i.e.,c̄˙i = c̄zi i . π̄˙ c = {π̄˙ cl , π̄˙ cr } are the sub-component weights foreach codebook component, governed by α. The choice of theauxiliary parameter distributions follows Chang et al. [21],which is reflected in the way we sample these variables in
Section IV.
B. Sampling Algorithm
Chang et al.’s DPMM sampler [21] is an instantiated-weightsampler that combines non-ergodic Markov chains into an
ergodic chain and proposes splits from learned sub-clusters
and merges of clusters. Their algorithm runs a Gibbs sampler,which samples the parameters and weights of each cluster
and its sub-clusters, followed by sampling the cluster andsub-cluster labels for each data point. The Gibbs sampler isrestricted to non-empty clusters. After Gibbs sampling, a split
and merge sampler proposes with an MH algorithm to either
split or merge clusters into new clusters. Gibbs sampling andsplit and merge sampling iterate until convergence or until astop criterion is fulfilled.
Our proposed sampler uses a similar structure. Algorithm 1
is an outline of our algorithm in pseudo-code. We combine
a restricted Gibbs sampler with a split, merge and switch
sampler for clusters and a split and merge sampler for

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

5components. The Gibbs sampler samples the parameters andweights of each codebook component and its sub-components,
the weights of each cluster, and the weights of each clustercomponent. This is followed by sampling the cluster, component and sub-component labels for each data point. A split,
merge and switch sampler proposes to either split or merge
clusters or to move a cluster component from one cluster toanother. A split and merge sampler for components proposesto either split or merge codebook components. Illustrationsof the possible component and cluster moves are given in
Fig. 6. Gibbs sampling, split, merge and switch sampling forclusters, and split and merge sampling for components iterate
until convergence or until a stop criterion is fulfilled.
The Gibbs sampling steps and the split and merge moves
for the codebook components are equivalent to the original
DPMM sampler of Chang et al. [21], and the codebook
together with the global component weights is exactly themodel that the original sampler would infer.
The following sections explain in detail the individual
non-ergodic samplers that make up our proposed algorithm.
The non-restricted Gibbs sampler is explained in Section IV.
Section V explains the component split and merge sampler
which are technically identical with the sampler in [21], but
applied to the codebook components. Because changes of thecodebook components lead to changes of the clusters, we
introduce novel update steps for clusters. These are explained
accordingly in the respective subsections. Lastly, we propose
cluster split, merge and switch moves in Section VI. Fig. 5
is an illustration of how our proposed DPMoMM algorithm
behaves, compared to a DPMM.

IV. R ESTRICTED G IBBS S AMPLING
In this section we lay out the details of the restricted Gibbs
sampler that we employ. The Dirichlet process uses an infinite
length prior on the cluster labels zi , cluster component labelsc̃ki and codebook component labels ċi . However, any label canonly point to a finite number of entities, i.e., the clusters andthe components that exist in any current state of the model.
Because the restricted Gibbs sampler does not create newclusters and components itself, the dimensions of the infinite
vectors π, π̃k , π̇ and θ are technically finite during Gibbs
sampling. Posterior distributions of weights are conditionedon the assignments of data points. The restricted conditional
distributions of the DPMoMM are

with
Nk =
Ṅc =

N
X
i=1
K
X

Nkc =

1zi =k ,NX

1 zi =k ,c̃ki =c
i=1
( c
Nk +β
if Nkc > 0,
Ck
Bkc =

Nkc ,

0

k=1

else,

(7)

p(π̃k |z, c̃, α, β) = Dir(Bk1 , . . . , BkC , α),(8)p(π̇|ċ, α) = Dir(Ṅ1 , . . . , ṄC , α),p(θc |x, ċ, λ) ∝ f ({x}c , θc )f (θc , λ),

(9)
(10)

p(zi = k|x, π, θ) ∝ πkCXπ̃kc f (xi , θc ),(11)(π1 , . . . , πK , πK+1 ) ∼ Dir(N1 , . . . , NK , β),(16)(π̃k1 , . . . , π̃kC , π̃kC+1 ) ∼ Dir(Bk1 , . . . , BkC , α),(17)(π̇1 , . . . , π̇C , π̇C+1 ) ∼ Dir(Ṅ1 , . . . , ṄC , α),

(18)

∝

θc ∼ f ({x}c , θc )f (θc , λ),
!
C
X
∝
c
zi ∼
πk
πk f (xi , θc ) 1zi =k ,
K
X

c=1
C
X
∝
c̃ki ∼
πkc f (xi , θc )1c̃ki =c ,
c=1

(12)
(13)

(19)
(20)

k=1

(21)
(22)

with
πK+1 = 1 −π̃kC+1 = 1 −c=1p(c̃ki = c|x, z, π̃, θ) ∝ π̃zci f (xi , θc ),p(ċi = c|x, π̇, θ) ∝ π̇c f (xi , θc ),(15)

where K is the current number of non-empty clusters, C is thecurrent number of non-empty codebook components and Ck isthe current number of non-empty cluster components in cluster
k. x is the vector of data points, z is the vector of cluster
labels, c̃ is the vector of cluster component labels, ċ is thevector of codebook component labels, and π̃ = {π̃1 , . . . , π̃K }.
{x}c denotes all data points assigned to the global codebook
component c. f (·) is a particular parametrized probability
density function. E.g., f ({x}c , θc ) is the likelihood of the datasubset {x}c given the cluster parameters θc , and f (xi , θc ) is
the likelihood of the single data point xi given θc . 1(·) is anindicator function that equals to 1 if the condition (·) holds
true and is 0 otherwise.
Given these probabilities, the sampling is as follows. Conditioned on the labels in the current state, sample the parametersof each codebook component, and all cluster and component
weights. Conditioned on all cluster and component parametersin the current state, for each data point, sample a label for acluster, then sample a label for a component within the cluster.
The conditional distribution in (11) shows that the generative
process described further above prefers cluster components tobe neighbors in the feature space. During the model sampling
described in this and the following sections, the algorithm willcause the clusters to be groups of nearby components so as tomaximize the likelihood of the data.
For our proposed algorithm, we lay out the Gibbs sampler asfollows. Given (7)-(13), the posterior distributions of weights,
labels and component parameters are expressed asċi = c̃zi i ,p(π|z, β) = Dir(N1 , . . . , NK , β),(14)π̇C+1 = 1 −

K
X
k=1
C
X
c=1
C
X

πk ,

(23)

π̃kc ,

(24)

π̇c .

(25)

c=1

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

×× × ××
××
××××
× ×
×
×
×
×××
×
×
×××
× ×××
×××
×
×××
× ×××
×
×
×××××
×
× ××
××
×
×
×
×
×
×
××
××
×××××
×
×
×××
××××
××
×××
××
×× ×
×
×
×
×
×
×
×
×
×
×
××××
× ××××
××××
×
× ×
×
×
× × ××
×
××
×
××
×
× ×
×
×
×
×
×
××

× ×
×× ×
××
××× ×
× ×× × × ×××××
× ××× ×
× ×
××
×
×× ×
×
× ×
×
×
×
×
×
×
××
××
×× ×
× ××
×
×
×
××
×
×
×
×
××××
×
××
×
× ××
× ××
×
×
××
××
×
××
×
×
×××
××× ×

×
××××
×
×
× ××
×
××
× ××
×
×
×
××
××
××
××
×
×
×× ××
×
××× ×
××
×
×
×
××
×
×× ××××××
× ×× × ×
××
×
××
× ××

× ×
× ××
×
××
×
×
××
×
××××
×× × × ×
× ×××
×
×× ××
×
×××
×
×
×
×× × ×× ×××××× ××××× ×××× ×××× ×
×
××
××
×
×
×
× ××
×× ×
× ××
×
×××
×
×
×××
××
×
×
××
× × ×× ××
× ××
××
×××
×
×
× ××
×××× ×
× ××
×
×
×
×
××
×
×
×
×
×
××
×××
×
×××
×
× ××××
×
× ×
××
××
××
× ××

6

×
× × ××
××
××××
× ×
×
×
×
×××
×
×
× ××
× ×××
×××
×
×××
× ×××
×
×
×× ×××
×
× ××
××
×
×
×
×
×
×
××
××
×××××
×
×
× ××
××××
××
×××
××
×× ×
×
×
×
×
×
×
×
×
×
×
× ×××
× ×× ××
××× ×
×
× ×
×
×
× × ××
×
××
×
××
×
× ×
×
×
×
×
×
××

× ×
×× ×
××
××× ×
× ××× × × ×××××
× ××× ×
× ×
××
×× ×
×
× ×
×
×
×
×
×
×
××
××
×× ×
× ××
×
×
×
××
×
×
×
×
××××
×
××
×
× ××
× ××
×
×
××
××
×
××
×
×
×××
××× ×

×
××××
×
×
× ××
×
××
× ××
×
×
×
××
××
××
××
×
×
×× ××
×
××× ×
××
×
×
×
××
×
×× ×××× ××× ×× × ×
××
×
××
× ××

× ×
× ××
×
××
×
×
××
×
××××
×× × × ×
× ×××
×
×× ××
×
×× ×
×
×
×
×× × ×× ×××××× ××××××××
× ×××× ×
×
××
××
×
×
×
× ××
×× ×
× ××
×
×××
×
×
×××
××
×
×
××
× × ×× ××
× ××
××
×××
×
×
× ××
×××× ×
× ××
×
×
×
×
××
×
×
×
×
×
××
×××
×
× ××
×
× ××××
×
× ×
××
××
× ××
××

(a)

×
× × ××
××
××××
× ×
×
×
×
×××
×
×
×××
× ×××
×××
×
×××
× ×××
×
×
×××××
×
× ××
××
×
×
×
×
×
×
××
××
×××××
×
×
×××
××××
××
×××
××
×× ×
×
×
×
×
×
×
×
×
×
×
××××
× ××××
××××
×
× ×
×
×
× × ××
×
××
×
××
×
× ×
×
×
×
×
×
××

× ×
×× ×
××
××× ×
× ××× × × ×××××
× ××× ×
× ×
××
×× ×
×
× ×
×
×
×
×
×
×
××
××
×× ×
× ××
×
×
×
××
×
×
×
×
××××
×
××
×
× ××
× ××
×
×
××
××
×
××
×
×
×××
××× ×

×
××××
×
×
× ××
×
××
× ××
×
×
×
××
××
××
××
×
×
×× ××
×
××× ×
××
×
×
×
××
×
×× ××××××
× ×× × ×
××
×
××
× ××

(b)

× ×
× ××
×
××
×
×
××
×
××××
×× × × ×
× ×××
×
×× ××
×
×××
×
×
×
×× × ×× ×××××× ××××× ×××× ×××× ×
×
××
××
×
×
×
× ××
×× ×
× ××
×
×××
×
×
×××
××
×
×
××
× × ×× ××
× ××
××
×××
×
×
× ××
×××× ×
× ××
×
×
×
×
××
×
×
×
×
×
××
×××
×
×××
×
× ××××
×
× ×
××
××
× ××
××

(c)

Fig. 5. Illustration of the algorithm. Components are drawn with solid lines, sub-components with dotted lines. Only some exemplary sub-components
are illustrated. (a): Single component clusters inferred by a DPMM. This also corresponds to the codebook of the DPMoMM; (b). Clusters inferred by a
DPMoMM, where components of the same color belong to the same cluster; (c): The same DPMoMM after a component split in the upper right cluster anda component merge in the lower right cluster. The original clusters are duplicated, and the bi-colored components are now shared across clusters.
∝

component splitXc̄ki ∼

Sec. V-Aπ̄c̃ski f (xi , θ̄c̃ski )1c̄ki =s ,c̄˙i = c̄zi i ,component merge
Sec. V-B

Sec. VI-A

cluster merge
Sec. VI-B

(31)

with
Nkc,s =

cluster split(30)s∈{l,r}NX
i=1

1 zi =k ,c̃ki =cc̄ki =sṄcs =K
X

Nkc,s ,

(32)

k=1

where {x}lc and {x}rc denote the subsets of data points thatare assigned to the left and right sub-components of c. Note
that, analogous to (22), the labels c̄˙i are not sampled explicitly
to not break the assignment of data points to sub-components
on the global level.

cluster switch
Sec. VI-C

Fig. 6. Overview of the component moves (top) and cluster moves (bottom).

Note that the codebook component labels ċi are not sampled
explicitly but derived from the cluster component labels.
This is done to not break the assignment of data points tocomponents on the global level; a data point that is assigned
to component c within a cluster must also belong to component
c within the codebook for the split and merge sampling logic
to work properly.
For the split and merge sampling steps described in Section V, we make use of auxiliary variables that are jointly
sampled with the regular variables. These auxiliary variables
describe the sub-components which augment all the regularcomponents. The sampled sub-components serve as good
split candidates for eventual component splits. The auxiliary
variables are sampled as follows:
Nkc,l + α Nkc,r + α
,
),
2
2
l
r
Ṅ + α Ṅc + α(π̄˙ cl , π̄˙ cr ) ∼ Dir( c
,
),
2
2
∝
θ̄cl ∼ f ({x}lc , θ̄cl )f (θ̄cl , λ),

l
r
(π̄kc
, π̄kc) ∼ Dir(∝θ̄cr ∼ f ({x}rc , θ̄cr )f (θ̄cr , λ),

(26)
(27)
(28)
(29)

V. C OMPONENT S PLIT AND M ERGE S AMPLER
The split and merge moves for components are performedon the global codebook level and therefore rely on thecodebook level variables that are jointly sampled with theother model variables. The moves are designed for efficiency
by reducing the overhead of computational costs during the
MH step and enabling parallelization across components.
Components are equipped with auxiliary variables for subcomponents. The parameters of the sub-components are sampled in the same fashion as the parameters for the regular“parent” components. Conveniently, samples for the variablesof the regular components can be obtained by sampling theauxiliary variables, since we draw from a joint parameter
space.
For performing split and merge moves for components withan MH-MCMC method, candidate moves, or proposals are
required. Let O = {π̇, θ, ċ} be the set of component variables˙ θ̄, c̄}˙ be the set of sub-component variables. We
and Ō = {π̄,ˆ } for components
propose a new set of random variables {Ô, Ō
and sub-components and compute the Hastings ratio of theformˆ)ˆ |x, ċ)ˆ q(O, Ō|Ô, Ōp(Ô, x)p(Ō
HR =
,
(33)
ˆ |O, Ō)p(O, x)p(Ō|x, ċ) q(Ô, Ō
where q is called the proposal distribution and x denotes thecollection of all observations. The Hastings ratio weights thestate of the model before and after actually performing a move,

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 20157with the numerator standing for the post-move state and thedenominator the pre-move state. As can be seen requires theHastings ratio a reverse proposal to the proposed move. In thecase of proposing a split move, that would be a merge, andvice versa. A proposed move is accepted with the probabilitymin(1, HR).

(34)

A. Split Moves

sampler is also called deferred MH sampler [21] and conveˆ |x, Ô) = p(Ō
ˆ |x, Ô). Components are marked
niently sets q(Ōas “splittable” after the variables show signs of a burn-in,which is the case when the likelihood f ({x}lc , θcl )f ({x}rc , θcr )for all data points assigned to component c begins to oscillatewith the iterations. It can be shown [21] that the Hastings ratio
for a component split can be expressed asˆ )Γ(Ṅˆ )ˆ αΓ(Ṅˆ p(ċ)p(x|ċ)mnp(x|ċ) p(ċ)Γ(Ṅc )ˆ )Γ(Ṅˆ )ˆ θ̂m )f (x|ċ,ˆ θ̂n ) αΓ(Ṅf (x|ċ,mn=.
f (x|ċ, θc )Γ(Nc )

HRc-split =The sub-components are utilized as good split move candidates for the MH algorithm. Proposing a split move is typicallya non-trivial task. The construction of a prospective move isnecessary for an MH framework, where a proposal is weighedagainst the status quo and accepted or rejected with a certain
probability. Theoretically, any kind of split proposal can lead toan ergodic chain [21]. However, proposals with low probabilityof being accepted unnecessarily increase the computationalload, since in case of a rejection, all previous computational
efforts are in vain. Iterative fitting of sub-components with thehelp of the auxiliary variables introduced above circumvents
the risk of wasted computational time. The sub-components
are sampled jointly with the normal components. During
the MH-step, sub-components pose good proposals for split
moves. Moreover, split move computations can be parallelized
across components.
The proposal distribution for proposing a component split
move with the help of the auxiliary variables is defined asfollows. First, a split or merge move Q ∈ {Qcc-split , Qm,nc-merge }is selected randomly. Qcc-split denotes a move for splitting
component c into m, n, and Qm,nc-merge is a merge of components
m, n into c. New sets of model variables are sampled asfollows, each conditioned on Q.
If Q = Qcc-split :ˆ m , {ċ}ˆ n ) = split (ċ, c̄),
˙
({ċ}
c
l
ˆ
ˆ
(π̇m , π̇n ) = π̇c · (π̄˙ c , π̄˙ cr ),ˆ˙ˆ c̄),(θ̂m , θ̂n ) ∼ q(θ̂m , θ̂n |x, ċ,ˆ , Ōˆ ) ∼ p(Ō
ˆ , Ō
ˆ |x, ċ),
ˆ
(Ō

(35)
(36)
(37)

(43)

Splits of codebook components affect the clusters that
contain them as cluster components. All clusters that containa split component require an update of their variables. The
update step for clusters is subject to design, as there arevarious ways to execute it. The naive approach is to split
the affected component within each cluster and update thecorresponding variables. We opted for an update scheme thatkeeps the cluster sizes unchanged so as to separate codebook
growth and mixture of mixtures growth entirely.
After performing a component split, we proceed as follows.
All clusters that contain the respective component are duplicated with a function({c̃ˆk }c , {c̃ˆk0 }c ) = duplicateck (c̃, c̄),(44)

where c̄ is the vector of cluster sub-component labels.
The original cluster k keeps all unchanged components
and the “left” split result. The duplicate k 0 also keeps all
unchanged components and the “right” split result. Data points
with zi = k that are assigned to the split component canbe reassigned so that the ones labeled with c̄ki = l belong
to cluster k, and data points with and c̄ki = r now belong
to cluster k 0 . Data points that are assigned to unchanged
components within the duplicated cluster however can not bereassigned unambiguously. Because a data point can neverbelong to two clusters at the same time, the duplication
automatically invalidates the labels zi and weights π. To reestablish a valid state for the sampler, we re-sample all variables. Therefore, the split sampler for components is the last
step before the next iteration of the restricted Gibbs samplingin our implementation, as can be seen in Algorithm 1. Fig. 5band 5c illustrate this step.

n

(38)

ˆ c = merge (ċ),{ċ}m,nˆc = π̇m + π̇n ,
π̇

(39)

ˆ˙
ˆ c̄),
θ̂c ∼ q(θ̂c |x, ċ,ˆ ∼ p(Ōˆ |x, ċ).
ˆ
Ō

(41)

B. Merge Moves(42)

The Hastings ratio for a proposed component merge also
requires a reverse proposal. Where there is only one way tomerge two sets of label assignments, there are 2Ṅc −1 −1 ways
to split a set of labels into non-empty partitions. However,since split proposals in this algorithm are determined by the
sub-components, the Hastings ratio will be zero if the labels
after a proposed reverse split do not match the pre-merge
labels. Therefore, the probability for accepting a component
merge rapidly diminishes with increasing number of assigned
data points. This behavior is approximated by automaticallyrejecting all component merges.

m

n

m

ˆ , Ṅ
ˆ ).
with (π̄˙ cl , π̄˙ cr ) ∼ Dir(Ṅmnm,nIf Q = Qc-merge :

c

c

(40)

The function splitc (·) splits the label assignments of component c so that labels are assigned to the two new components
m and n, whereas the function mergem,n (·) does the reverse
move and merges the label assignments of components m andn so that the respective data points are assigned to a newcomponent c. Sampling new parameters given the new label
assignments is done with the restricted Gibbs sampler. The
sub-component auxiliary variables are sampled jointly withthe variables for the regular components. This specific joint

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 20158

In order to mitigate for slow convergence in certain situations, a random merge sampler is introduced instead to propose
component merge moves whose reverse move is a random
split, in contrast to the sub-component based deterministic split
proposals. Two random components are sampled and a merge
proposal is computed. The reverse split proposal is generated
by a random partitioning of the data points assigned to therespective component. The split proposal will generally have adiminishing acceptance probability, whereas the corresponding
merge move is much more likely to be sensible. The Hastings
ratio for a random merge proposal is as follows:ˆ p(ċ)ˆ Γ(α)Γ(Ṅm )Γ(Ṅn )p(x|ċ)HRc-merge =ˆ )Γ( α )2p(x|ċ) p(ċ) Γ(α + Ṅ
c

2

ˆ )
ˆ
p(x|ċ)Γ(ṄΓ(α)Γ(Ṅm )Γ(Ṅn )c=ˆ )Γ( α )2p(x|ċ) αΓ(Ṅm )Γ(Ṅn ) Γ(α + Ṅc2

ˆ )ˆ θ̂c )
Γ(α)Γ(Ṅf (x|ċ,c=. (45)ˆ )Γ( α )2f (x|ċ, θm )f (x|ċ, θn ) αΓ(α + Ṅc2

A random split of component ċˆ is sampled for the reverse
split proposal, therefore the weights for the split results are
Dirichlet distributed.
Analogous to the split moves, all clusters that contain amerged component need to be updated. Fortunately, the cluster
updates after component merges are much simpler. All clusters
that contained any of the two merged components replaces therespective component with the merge result. Cluster labels zi
and weights π remain unchanged. Cluster component labelsc̃ki and weights π̃k are updated if cluster k contains bothcomponents involved in the merge:{c̃ˆk }c = mergem,n (c̃k ),ˆkc = π̃
ˆkm + π̃
ˆkn .
π̃

(46)
(47)

New values for cluster sub-component auxiliary variablesŌk = {π̄k , c̄k } are sampled for all clusters to be updated
according toˆ c ∼ p(Ō
ˆ c |x, c̃ˆ ).
Ō
(48)
k
k
k
Again, Fig. 5b and 5c illustrate the outcome of this step.

is proposed by any of the possible moves. The Hastings ratio
for a move is of the formˆ)
ˆ |x, c̃)ˆ q(M, M̄ |M̂ , M̄p(M̂ , x)p(M̄
HR =
.
(49)
ˆ |M, M̄ )p(M, x)p(M̄ |x, c̃) q(M̂ , M̄As before, a proposed move is accepted with the probability
defined in (34). Note that the component parameters θc are
not subject to updates during cluster moves. That is because
neither assignments of data points to codebook components
nor codebook component parameters are modified.
A. Split Moves
Analogous to component splits, where good split candidates
are provided by auxiliary sub-components, we resort to thecluster components themselves as support for producing good
proposals. To generate a good split proposal for cluster k,
we consider all 2Ck −1 − 1 possible non-empty partitions into
two separate mixtures and propose the most promising split,
according to the Hastings ratio.
The computational overhead grows with the number ofcomponents in a cluster. For instance, a cluster with 16 components can be partitioned into two clusters in 32.767 different
ways, which would require a same amount of computationsfor Hastings ratios. In order to control the expected maximum
computational load for splits, we introduce a parameter cmax
into the algorithm which caps the maximum size of clusters.
Setting this parameter accordingly prevents clusters to growtoo large. By setting the value arbitrarily high, mixtures may
grow to any size. To maintain computation feasible for verylarge clusters as well, we could also limit the consideredpartitions to a random subset of possibilities smaller than the
Stirling number above.
The cluster split sampler’s design essentially follows theconsiderations of the component split sampler in Section V-A.
As is the case for components, we can also parallelize thecluster split move computations. The proposal distribution fora cluster split move is defined as follows. First, we randomlya,b
select a split move or a merge move Q ∈ {Qmk-split , Qk-merge },where Qm
k-split denotes a move for splitting cluster m into a, b,
and Qa,b
k-merge is a merge of clusters a, b into m. New sets ofmodel variables are sampled as follows, conditioned on Q.
If Q = Qm
k-split :VI. C LUSTER S PLIT, M ERGE AND S WITCH S AMPLER
Cluster split, merge and switch moves modify the assignments of components to clusters. A split move splits a cluster –which is a mixture of components – into two smaller mixtures.
A merge move merges two clusters into one larger mixture.
A switch move moves one component from one cluster toanother cluster. In all these cases the components themselves,i.e., their parameters, are not modified. Cluster moves rely on
the local, cluster dependent component parameters to produce
good move proposals for a MH step. Cluster moves are
efficient because they can be easily computed based on theexisting data partitions.
Let M = {π, π̃, z, c̃} be the set of cluster and clustercomponent variables and M̄ = {π̄, c̄} be the set of cluster subˆ}
component variables. A new set of random variables {M̂ , M̄({ẑ}a , {ẑ}b ) = splitm (z, c̃),

(50)

a
b
(π̂a , π̂b ) = πm · (πm
, πm
),
1
C
1
ˆa , . . . , π̃ˆa ) ∼ Dir(B̂a , . . . , B̂aC ),(π̃ˆb1 , . . . , π̃ˆbC ) ∼ Dir(B̂b1 , . . . , B̂bC ),
(π̃

(51)
(52)
(53)

a
a
b
with (πm, πnb ) ∼ Dir(N̂m
, N̂m
).
a,b
If Q = Qk-merge :{ẑ}m = mergea,b (z, c̃),π̂m = πa + πb ,1C1Cˆm , . . . , π̃ˆm(π̃
) ∼ Dir(B̂m, . . . , B̂m
),
with

(
B̂kc =

N̂kc +βCk0

if N̂kc > 0,else.

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

(54)
(55)
(56)

(57)

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

9

The function splitm (·) splits the label assignments of clusterm so that labels are assigned to the two new clusters aand b, whereas the function mergea,b (·) does the reverse
move and merges the label assignments of clusters a andb so that the respective data points are assigned to a newcluster m. The component labels c̃ki do not require an update
since component IDs are valid globally. To promote a morestable splitting behavior, clusters are marked “splittable” if
all components within the respective clusters are also marked“splittable” (see Section V-A). We express the Hastings ratio
for a cluster split proposal as follows:N̂a)Γ(β NN̂mb )p(x|ẑ) p(ẑ) Γ(β + Nm )Γ(β Nm
HRk-split =p(x|c̃) p(z)Γ(β)Γ(N̂a )Γ(N̂b )

=
=

N̂a
)Γ(β NN̂mb )p(x|ẑ) βΓ(N̂a )Γ(N̂b ) Γ(β + Nm )Γ(β Nmp(x|c̃)Γ(Nm )
Γ(β)Γ(N̂a )Γ(N̂b )N̂a
ˆ Θ̂a )f (x|c̃,ˆ Θ̂b ) βΓ(β + Nm )Γ(β N)Γ(β NN̂mb )f (x|c̃,m,f (x|c̃, Θm )Γ(Nm )Γ(β)(58)

many component splits. Another factor is the maximal cluster
size cmax . The larger clusters can get, the more likely they
duplicate due to eventual component splits. In cases where alarge amount of clusters overlap (that is, many clusters sharethe same components) and therefore cover the same region in
the feature space, the algorithm can suffer from convergenceissues due to high ambiguity during label sampling.
To mitigate this issue, we developed a switch move samplerthat supports algorithm convergence. A switch move is anoperation, where all data points that are assigned to component
a in cluster m are re-assigned to the same component a incluster n. Given a pair of clusters, we consider a switch movefor a in both directions and propose the direction which ismost promising, according to the Hastings ratio. Switch moveproposals are sampled by a random switch sampler, where
for a prospective split, two random clusters are sampled andthe proposal is computed. New sets of model variables are
sampled as follows:ˆ m , {c̃}ˆ n ) = switcham,n (z, c̃),({ẑ}m , {ẑ}n , {c̃}aa(π̂m , π̂n ) = (πm −πm π̃m, πn + πm π̃m
),
1
C
1
C
ˆ
ˆ
(π̃m , . . . , π̃m ) ∼ Dir(B̂m , . . . , B̂m ),ˆn1 , . . . , π̃ˆnC ) ∼ Dir(B̂n1 , . . . , B̂nC ).
(π̃with Θ̂k = {π̂k , {θ}k }, and {θ}k being the parameters of allcodebook components that are also cluster components in k.
B. Merge Moves
For a prospective merge, two random clusters are sampled
and a merge proposal is computed. A cluster merge is only
permitted if the component size of the merge result is notexceeding cmax . The reverse proposal is a random partition of
the cluster components into two separate mixtures. Analogous
to the cluster split proposal above, the Hastings ratio for acluster merge proposal is expressed as follows:f (x|ẑ, Θ̂m )HRk-merge =f (x|z, Θa )f (x|z, Θb )Γ(Na + Nb )Γ(β)×.
Na)Γ(β N̂Nb )βΓ(β + Na + Nb )Γ(β N̂
m

(60)
(61)
(62)
(63)

The function switcham,n (·) re-assigns all data points incomponent a of cluster m to component a of cluster n by
updating the component labels c̃ki and the cluster labels zi .
A switch move can be interpreted as splitting one component a off of a cluster m and merging a single-component
cluster with a cluster n that already contains a as component.
The Hastings ratio for a switch proposal is therefore expressed
as follows:
HRswitch
=

(59)

m

As mentioned earlier, (11) suggests that the likelihood of
a cluster is higher if the components of the respective cluster
are located closer to each other in the feature space. During
inference, β controls the importance of proximity for grouping
nearby components into clusters. The Hastings ratios (58) and(59) suggest that with larger β the probability of accepting
a split proposal becomes higher, and the probability of thereverse merge move becomes smaller.

a)Γ(N̂m )Γ(Nn ) Γ(β + Nm )Γ(β + Nn )p(x|ẑ) βΓ(Nmp(x|z)
Γ(Nm )Γ(Nn )
Γ(β)
Na

×
×
×
=

m
n
m
)Γ(β NmN̂+N)Γ(β NmN+N)
Γ(β Nm +Nnnna )Γ(N̂ )Γ(N )Γ(Nmmn
Γ(N̂m )Γ(N̂n )Γ(β)
a )Γ(N̂ )Γ(N ) Γ(β + N̂ )Γ(β + N̂ )
βΓ(Nm
m
n
m
n
a
Γ(Nm
)Γ(N̂m )Γ(Nn )nmΓ(β NmN̂+N)Γ(β NmN̂+N)nnp(x|ẑ) Γ(N̂m )Γ(N̂n ) Γ(β + Nm )Γ(β + Nn )p(x|z) Γ(Nm )Γ(Nn ) Γ(β + N̂m )Γ(β + N̂n )Na

C. Switch Moves
The component split and merge moves (see Section V)
produce clusters with shared components due to the specificsof the cluster update steps. After a component split, allclusters that contain the split component are duplicated. Aftera component merge, all clusters who contained the previous
components now share the new component (unless a cluster
already contained both original components).
The number of clusters that share the same component
can grow very quickly if the algorithm decides to perform×mnΓ(β Nm +N)Γ(β NmN+N)nnnΓ(β NmN̂+N
)
n

.

(64)

The reverse move requires to split the data points thatare assigned to component a in cluster n and assign thesedata points to cluster m. For the same reasons than in thecase of normal component merges (see Section V-B), theprobability for the reverse proposal quickly approaches zerowith increasing data size and the reverse move will be rejected.
We approximate this behavior by automatically rejecting all
reverse switch moves.

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

VII. U SE C ASE : U NSUPERVISED S UBWORD M ODELINGIn the following section, we will demonstrate the benefitsof our DPMoMM sampler on real speech data in the use caseof unsupervised subword modeling.
The objective of unsupervised subword modeling is toconstruct a representation of speech sounds that is robust
to variation within and across speakers and that maximizes
class discrimination [8]. Previous work of Chen et al. [10],[11] and ourselves [12]–[14] already achieved good results
using a DPMM sampler to tackle this task in the context
of the zero resource speech challenges [8], [9]. The generalprocedure is to cluster speech feature vectors into classes
by sampling a Dirichlet process Gaussian mixture model(DPGMM) [5], which is an infinite Gaussian mixture model(IGMM) [3]. Speech would then be represented by framelevel posteriorgrams, for instance, or simply by a sequence of
textual class labels, where the classes are the components inthe sampled mixture.
The official evaluation metric for the zero resource speechchallenge is the minimal pair ABX phone discriminability
between phonemic minimal pairs [33]. Another popular metricfor evaluating speech representations especially with respect
to information and sequentiality is the normalized mutual
information (NMI). We use both metrics to evaluate the output
of our sampler in the following section. The details of theevaluation measures are explained in this section.

10and gives a measure of how good X can be predicted, given
the knowledge about Y . In our use case, X corresponds to
the random variable over the “true” distribution of sounds,
approximated by the sequence of phones for any target data
set. Y is the random variable over the estimated distributionof sounds, given for instance by the label output of a DPMM
sampler for the same data. H(X) is the entropy of thetrue transcription and is used as normalizing factor in the
denominator.
NMI(·) is not symmetric, which makes the comparison of
random variables that are defined over very different inventory
sizes difficult. To guarantee a fair comparison, we suggest to
use a symmetric NMI [34] of the form
NMImax (X; Y ) =H(X) − H(X|Y )max(H(X), H(Y ))(67)

With this, we have a fair measure of correlation at hand forthe frame-level phone transcriptions of any target data and theframe-level label output of a DPMM or DPMoMM.
VIII. E XPERIMENTAL E VALUATION
A. Data1(1d(A,X)<d(B,X) + 1d(A,X)=d(B,X) )(65)
2
where d(·, ·) is the dynamic time warping (DTW) distancedefined over sequences of frame based speech representations(posteriorgrams, textual labels, etc.). Any proper distance measure can be used for computing the DTW distance. We evaluate
two types of representations, textual labels and posteriorgrams.
For the comparison of label sequences, we use the Levenshtein
distance (ABXls ). For evaluating posteriorgram sequences, we
use the Kullback-Leibler divergence (ABXkl ).
We collect discrimination errors for all possible pairings ofphone triplets and average them over all contexts for a given
pair of central phones, over all pairs of central phones andover all speakers.

We use two separate data sets known from the zero resource
speech challenge [8]. Because the sets vary in size, language
and speech quality, comparable experimental results should
be a good indicator for the robustness of our sampler. One set
contains spontaneous, conversational interview-style American
English (4h 59min), extracted from the Buckeye corpus [35].
The other set contains carefully uttered, read speech in Xitsonga (2h 29min), a southern African Bantu language. The
latter is an excerpt of the NCHLT corpus [36]. The referencesfor English and Xitsonga contain 165k and 72k phone-triplet
annotations for 39 and 53 unique center phones, respectively.
From each of the two data sets, we extract two sets of speechobservations. Specifically, from each data set, we extract two
types of frame-wise feature vectors using the Kaldi speech
recognition toolkit [37]. We can use about 1.7M frames for English and 0.8M frames for Xitsonga as input to the DPMM andDPMoMM samplers. The frame width is 25 milliseconds andthe frame shift is 10 milliseconds. The first type of features isperceptual linear predictive (PLP) speech feature vectors [38]with first and second order derivatives (PLP+∆+∆∆). The
second type is stacked PLP vectors that were transformed
unsupervisedly by linear discriminant analysis (PLP+LDA)
with the method described in [13]. LDA is a commonly usedin speech recognition to optimize speech features towards
discriminability [39]. We conduct experiments for each of thetwo input feature vector types. Overall, we conduct all our
experiments four times, once for each data set and feature
vector type.

B. Symmetric Normalized Mutual Information

B. Procedure

The mutual information of two random variables is a measure for mutual dependence. Normalized mutual information
is defined as
H(X) − H(X|Y )NMI(X; Y ) =
(66)
H(X)

We use our DPMoMM sampler to cluster a set of speechfeature vectors into classes. In our use case, our algorithm
jointly samples Gaussian mixtures as well as a global codebook of Gaussians. This codebook is precisely what theoriginal DPMM sampler of Chang et al. [21] would infer. Due

A. ABX Phone Discriminability
Let A and X be two speech representations of soundcategories a, and B a speech representation of sound categoryb. The ABX phone discrimination error for categories a and
b is
XX X
1
c(a, b) = 1 −|a| · |b| · (|a| − 1)A∈a B∈b X∈a\{A}

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500

200
β

300

50

100

200
β

300

-0,01-0

0,02-0,03
3300

0-0,010-0,015

50

200
β

50

100

300

200
β

300

200
β

300

2
400

8

0,9-1

0,15-0,18

0,7-0,8

4
0,12-0,15

0,09-0,12

3

100

200
β

300

0-0,030,5-0,6
30,3-0,40,2-0,30,03-0,06
2400

4

0,6-0,70,4-0,5

3

0,06-0,09

0,01-0,02
50

100

(e)

4

0,04-0,05

0-0,010-0,1 50

0,8-0,90,03-0,04

2400

0,1-0,2

8

0,02-0,03

100

0,2-0,3
2
400

0,18-0,21

0,06-0,075
3

3

(d)

0,05-0,06
cmax

0,075-0,09

0,015-0,03

2
400

0-0,030,06-0,07

0,03-0,045

300

200
β

cmax

cmax

100

0,5-0,60,3-0,4

0,03-0,06

2
400

4

0,6-0,70,4-0,5

3

8

4

0,045-0,06

0,01-0,02

0,12-0,15
0,09-0,12

0,07-0,08

0,09-0,105
cmax

0,03-0,04

0,15-0,18(c)

0,105-0,12

200
β

50

8

4

0,7-0,8

4

0,06-0,09(b)

0,04-0,05

100

30-0,012
400

0,12-0,135

50

0,04-0,05

0,01-0,02

8-0,01-0

0,05-0,06

0,02-0,03(a)
0,05-0,06

0,18-0,21

4

0,03-0,04

3

0,8-0,90,21-0,24

0,06-0,07

cmax

100

2
400

0,07-0,08
4

8

0,9-1

0,24-0,27

cmax

3

0,12-0,135
0,105-0,12
0,09-0,105
0,075-0,09
0,06-0,075
0,045-0,06
0,03-0,045
0,015-0,03
0-0,015

cmax

cmax

4

8

8

0,08-0,09

cmax

8

cmax

8

50

100

200
β

300

0,1-0,22
400

0-0,1 50

100

200
β

300

2
400

JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

0,09-0,1
0,08-0,09
0,07-0,08
0,06-0,07
0,05-0,06
0,04-0,05
0,03-0,04
0,02-0,03
0,01-0,02
0-0,01-0,01-0 50

(f)
(g)
(h)
(i)
(j)
Fig. 7. Relative performance improvements when interpreting clusters instead of codebook components as acoustic classes. The inputs are 39 dimensional PLP+∆+∆∆ speech features. (a)-(e): Results on Xitsonga;(f)-(j): Results on English; (a),(f) Relative improvement of NMImax , (b),(g) Relative improvement of ABXls , (c),(h) Relative improvement of ABXkl , (d),(i) Relative unit length increase, (e),(j) Cluster to component
amount ratio. Paired t-tests on all ABX task outputs yield p  0.0001.

300

50

0,04-0,05

4
cmax

0,03-0,04

0,02-0,03
3

0,01-0,02
0-0,01
-0,01
0
0,02
0,01
0,03
0,04
0,06
0,05

200
β

cmax

0-0,01
100

200
β

300

2
400

-0,01
0
0,01
0,03
0,02
0,04
0,06
0,05

-0,01-0 50

100

300

2
400

0,09-0,1
0,08-0,09
0,07-0,08
0,06-0,07
0,05-0,06
0,04-0,05
0,03-0,04
0,02-0,03
0,01-0,02
0-0,01
-0,01
0,01
0
0,03
0,02
0,05
0,04
0,08
0,07
0,06
0,09
0,1
-0,01-0 50

200
β

300

2
400

0,7-0,8

40,2-0,3
100

200
β

300

0,1-0,2
0,2
0,1
0,4
0,3
0,6
0,5
0,9
0,8
0,7
10

2
400

0-0,1 50

0-0,01
3-0,01-0

-0,01
0
0,01
0,02

50

100

200
β

100

300

2
400

0,24-0,27

0,9-1

0,21-0,24

0,8-0,90,18-0,21

300

2
400

0,12-0,15
3

0,09-0,128

0,7-0,8

4

0,15-0,18

4

0,6-0,70,5-0,60,4-0,5

30,3-0,40,2-0,30,03-0,06
0
0,06
0,03
0,12
0,09
0,18
0,15
0,24
0,21
0,27

0-0,03200
β(e)0,06-0,09
2400

3

0,3-0,4

8

4

3

300

0,5-0,6(d)

0,01-0,02

200
β

4

0,6-0,70,4-0,5

3

8

4

100

8

0,9-1
0,8-0,9

(c)
8

cmax

0,05-0,06

100

30,01-0,02

(b)
8

50

0,02-0,03

8

cmax

200
β

3

cmax

100

2
400

4
0,03-0,04

0,3-0,33
0,27-0,3
0,24-0,27
0,21-0,24
0,18-0,21
0,15-0,18
0,12-0,15
0,09-0,12
0,06-0,09
0,03-0,06
0
0,09
0,06
0,03
0,15
0,12
0,21
0,18
0,27
0,24
0,3
0,33
0-0,03 50cmax

30,04-0,05
4
cmax

cmax

4

(a)

-0,01-0

8
0,05-0,06

cmax

8

0,07-0,08
0,06-0,07
0,05-0,06
0,04-0,05
0,03-0,04
0,02-0,03
0,01-0,02
0-0,01-0,01-0
-0,01
0
0,02
0,01
0,04
0,03
0,06
0,05
0,08
0,07

cmax

8
0,08-0,09
0,07-0,08
0,06-0,07
0,05-0,06
0,04-0,05
0,03-0,04
0,02-0,03
0,01-0,02
0-0,01
-0,01
0
0,03
0,02
0,01
0,05
0,04
0,07
0,06
0,09
0,08
-0,01-0
50

50

100

200
β

300

2
400

0,1-0,2
0,2
0,1
0,4
0,3
0,6
0,5
0,9
0,8
0,7
10
0-0,1 50

100

200
β

300

2
400

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

11

(f)
(g)
(h)
(i)
(j)
Fig. 8. Relative performance improvements when interpreting clusters instead of codebook components as acoustic classes. The inputs are 20 dimensional PLP+LDA speech features. (a)-(e): Results on Xitsonga;(f)-(j): Results on English; (a),(f) Relative improvement of NMImax , (b),(g) Relative improvement of ABXls , (c),(h) Relative improvement of ABXkl , (d),(i) Relative unit length increase, (e),(j) Cluster to component
amount ratio. Paired t-tests on all ABX task outputs yield p  0.0001.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

β=300

β=400

#clusters#components

250

53

18
3

19

15

16

8

12

9

150
100
50

0

1-4

5-8
Components in cluster

250 -2,64-2,66
200 -2,69-2,71
150 -2,74
-2,76
100
-2,79
50 -2,81-2,84
0 -2,86500 0#Components

200

0

100

200
300
Iterations

400

likelihood
-2,64
-2,66
-2,69
-2,71
-2,74
-2,76
-2,79
-2,81
-2,84
-2,86
100

200
300
Iterations(a)Fig. 9. Distribution of cluster sizes by the example of sampling Xitsonga
data for 500 iterations. Higher β values lead to more clusters with fewer
components. The same behavior was observed on the English data set.

to the joint sampling we can directly compare the modeling
quality of inferred Gaussian mixtures and the single Gaussians,
or in other words the DPMoMM and the codebook, i.e., a
DPGMM/IGMM.
After sampling the DPMoMM, the data can be represented
by frame-wise labels. We use the symmetric NMI to compare
the quality of label sequences. For that, we compute thesymmetric NMI once using the labels for clusters and once
using the labels for codebook components and calculate therelative improvement. In the same way, we also compare the
ABX phone discriminability using the frame-wise labels as
representation and calculate the relative improvement fromusing cluster labels instead of codebook component labels.
The ABX phone discriminability can also be computedfor posteriorgrams as representation for the data. In thatcase, either a posteriorgram over clusters or over codebook
components is computed for each speech frame. The two kindsof posteriorgrams are scored and compared to get a value
for the relative improvement by using cluster posteriorgrams
instead of component posteriorgrams.
We run every sampling for 1000 iterations. Each sampling
step is parallelized across 30 threads. For all conducted experiments we sample each model 5 times, score each output
and average the results. It is known that the influence of αdiminishes in very high data regimes [21]. Chen et al. [10]conducted an experience study and confirmed that the value of
α does not impact the outcome of sampling a DPMM givenhigh dimensional speech feature vectors. Their samples are
extracted from the same data sets that we use in our workand are similar in nature. In several informal experiments inconnection with earlier work we also observed this behavior
and therefore set α = 1 for all our experiments.
C. The Impact of β
We compare the use of clusters versus using the codebook
components as model for the underlying data. The lattercorresponds to output that the original sampler of Chang etal. [21] produces. We test on both data sets, English andXitsonga, and use either PLP+LDA features or PLP+∆+∆∆
as input. Fig. 7 and 8 show the relative improvements that our
proposed method achieved as contour plots.
We observed that using a very small value for the mixture
concentration parameter β tends to result in few sampled

400

Log likelihood (×107)

β=200
67

Log likelihood (×107)

β=100#Components

Number of clusters

β=50
70
60
50
40
30
20
10
0

12

500

(b)

Fig. 10. Sampling behavior on Xitsonga data when (a) not using switch
moves for cluster components, and (b) using switch moves. Without switch
moves, the number of clusters might grow fast, and convergence is slow. Withswitch moves, this issue does not occur. The same behavior was also observedon the English data set.

mixtures that each contain a maximum amount of components.
In other words, the sampler is over-confident in grouping components together based on minimal proximity. In contrast tosynthetic data, real data tends to be comprised of overlappingclasses. The aggressive grouping is one of the consequencesof this fact. With larger β, we observed that fewer components
are grouped together to form mixtures, which results in alarger number of clusters that contain fewer components.
Fig. 9 exemplarily plots the distribution of cluster sizes for the
Xitsonga data with different values for β. Fig. 7e, 7j and 8e, 8jshow that the number of clusters approximates the number ofcodebook components as β increases.
The behavior of the cluster inference dependent on β isbest explained by analyzing the sampling of weights during
restricted Gibbs sampling and the Hastings ratio for cluster
split and merge moves. According to (4)-(5), (7)-(8) and (16)(17), the distributions of cluster and the cluster component
weights are governed by β, whose impact is twofold. Itsvalue determines the probability mass that is reserved forgenerating a new cluster by the split sampler, and it regulates
the weights of the cluster components. A large β will motivate
the generation of more clusters and cause cluster component
weights to take on more similar values, therefore keeping
more cluster components alive for a longer time. This inturn encourages more cluster splits, which is also reflectedin the Hastings ratios for cluster moves. With larger β, (58)takes on a larger value, i.e., the probability of accepting a
split proposal becomes higher, and (59) takes on a smaller
value, i.e., the probability of the reverse merge move becomes
smaller. Intuitively, the effect is that only closely related
components remain grouped in form of a cluster, and lessdense clusters are likely to be split to form new clusters withless components. The resulting DPMoMM tends to be madeup by many clusters with mostly low amounts of components,and only few clusters with higher amounts of components, if
the data suggests so. A small value for β has the exact opposite
effect and a sampled DPMoMM will have few clusters withmostly large amounts of cluster components.
D. Convergence and the Switch Sampler
During our experiments, we found that under certain conditions, the number of clusters can grow rapidly and stay

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

13

TABLE IP ERFORMANCE COMPARISON OF STANDARD DPMM S AND PROPOSED DPM O MM S APPLIED TO BOTH DATA SETS .
PAIRED T- TESTS ON ALL ABX TASK OUTPUTS YIELD p  0.0001.
Features

Sampler

PLP+∆+∆∆DPMoMM

PLP+LDADPMoMM

PLP+fMLLR

DPMM [13]DPMoMM
DPMoMM

PLP+∆+∆∆DPMoMM

PLP+LDADPMoMM

PLP+fMLLR

DPMM [13]DPMoMM
DPMoMM

Clusters modeled bySingle Gaussians
Gaussian Mixtures
Single Gaussians
Gaussian Mixtures
Single Gaussians
Single Gaussians
Gaussian Mixtures
Single Gaussians
Gaussian Mixtures
Single Gaussians
Gaussian Mixtures
Single Gaussians
Single Gaussians
Gaussian Mixtures

β
cmax
K
C
K/C
Results for the Xitsonga data set
154
1.0
200
8
114
154
0.74
142
1.0
400
8
120
142
0.848
139
144
1.0
400
8
126
144
0.876
Results for the English data set
232
1.0
100
4
110
232
0.473
152
1.0
300
4
98
152
0.649
156
157
1.0
300
4
93
157
0.591

large for a long stretch of iterations. This is the case when
the allowed cluster size cmax is large and β is set to facilitate
many clusters with a large number of components. Under suchcircumstances, the duplication of clusters during the cluster
update step in the component split sampler (see Section V-A)
becomes more likely and more frequent, especially if multiple
components in the same cluster are subject to splitting.
We developed the switch move sampler to mitigate thisissue and support convergence. Fig. 10 exemplarily shows the
sampling behavior of our DPMoMM sampler on Xitsonga datawith and without using switch moves for cluster components.
As can be seen, the number of clusters might grow fast withoutswitch moves, and convergence is slow. With switch moves,
however, this issue does not occur. The same behavior wasobservable for both of our data sets. Without switch moves, aconsiderably larger amount of time would be required by the
sampler to converge. With switch moves, convergence is fast.
E. Use Case PerformanceFig. 7a, 7f and 8a, 8f show the relative improvementsof cluster label sequences over codebook component label
sequences on the symmetric NMI metric. The best results
are always achieved by allowing larger clusters. The optimalvalues for β seems to lie within a certain range. This range
is the same for our two data sets, English and Xitsonga,
but it seems to be input feature dependent. Sampling from
PLP+∆+∆∆ features benefits from a β value between 100
and 300, where sampling from PLP+LDA features might evenbenefit from a β larger than 400, which during our experiments
was the highest value that we tested. These observationstransfer to the evaluation of cluster and component labelswith the ABX phone discriminability test, as can be seen in
Fig. 7b, 7g and 8b, 8g.
Interestingly, we observed considerable performance improvements when we extracted cluster posteriorgrams andcompared them to component posteriorgrams with the ABX
discriminability test, especially with PLP+∆+∆∆ as input(see Fig. 7c, 7h and 8c, 8h). This is noteworthy because despite
the much lower dimensionality of cluster posteriorgrams,

NMImax

ABXls

ABXkl

avg. seg. len.

0.275
0.302
0.283
0.31
0.293
0.315

21.31%
18.44%
20.29%
18.74%
19.33%
18.06%

14.03%
12%
13.15%
12.47%
12.2%
12.44%
11.93%

2.16 frames
2.76 frames
2.22 frames
2.68 frames
2.29 frames
2.7 frames

0.232
0.241
0.248
0.262
0.253
0.266

28.85%
25.65%
25.88%
23.74%
25.28%
23.28%

18.99%
17.61%
16.18%
15.99%
15.7%
15.84%
15.47%

2.06 frames
2.37 frames
2.33 frames
2.73 frames
2.31 frames
2.64 frames

TABLE IIP ERFORMANCE COMPARISON OF ACOUSTIC UNIT RECOGNIZERS TRAINEDON DPMM AND DPM O MM LABELS . PAIRED T- TESTS ON ALL ABX TASKOUTPUTS YIELD p  0.0001.
Labels ofGaussians
Mixtures
Gaussians
MixturesUnits

NMImax
ABXls
ABXkl
Results for the Xitsonga data set
144
0.319
18.06%
12.09%
126
0.327
16.45%
11.17%
Results for the English data set
157
0.263
24.48%
16.5%
93
0.278
22.21%
15.76%

avg. seg. len.
2.59 frames
2.95 frames
2.77 frames
3.12 frames

performance not just equals, but even increases, comparedto the component posteriorgrams. This is an indicator thatthe clusters inferred by our proposed DPMoMM are notjust accumulations of related classes, but in fact a better
approximation to the real underlying multimodal distributionsof more complex classes in the data.
The number of inferred clusters does not necessarily correlate with the average length of sample sequences that use
the clusters as classes. Fig. 7d, 7i and 8d, 8i show that it is
the maximum number of allowed components per cluster cmax
that governs the average unit length, rather than β. The average
unit lengths is higher if the clusters are allowed to grow larger.
A too large β somewhat slows down this trend. The increased
unit lengths for larger clusters is an indicator that the groupingof components into mixtures allows the mixture of mixturesmodel to capture wider acoustic phenomena, an ability that is
highly valued in unsupervised subword modeling.
Overall, our results as depicted in Fig. 7 and 8 reveal thatwith some tuning of β and cmax , relative improvements on
all evaluation metrics can be as high as 13.5%. The contour
plots also suggest that larger mixtures might result in evenlarger improvements. It is a strong support for our initial
motivation for this work that good results can be achieved
with settings that lead to the inference of as little as 27% of
the amount of clusters than there are codebook components,
which considerably reduces model complexity. At the sametime, the found clusters can serve as models with average
durations of sound instances that are up to 31% longer. Table I

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

compares the performance of standard DPMMs and goodDPMoMMs that we sampled in absolute numbers. During
our experiments we sampled models that are best on one
metric, but sub-optimal on others. The exemplary DPMoMMs
in Table I perform reasonably well on all evaluation metrics
and give a good impression of what can be expected if theparameters are tuned reasonably. Note that we did not yet
exhaust the exploration of the hyper-parameter space. We
expect that even better results are possible with larger valuesfor β and especially cmax .
F. State-of-the-art Zero Resource Challenge PerformanceIn [13], we transformed PLP features with LDA, maximum
likelihood linear transforms (MLLT) [40], [41] and featurespace maximum likelihood linear regression (fMLLR) [42],[43] (denoted as “PLP+fMLLR”) – a method commonly used
for speaker adaptation – to improve the input to a DPMM
sampler. The extracted posteriorgrams achieved the to-datebest results on the zero resource speech challenge 2015
Xitsonga and English data sets.
We repeated these experiments with our novel sampler,using the hyper-parameters that we tuned for the PLP+LDA
features. Table I compares the previously published performance of [13] with using a DPMoMM for sampling instead.
We could infer more compact models having fewer clusters
and at the same time reduce the discriminability errors – theofficial evaluation method of the challenge – even further, thus
establishing a new state-of-the-art with our proposed method.
G. Acoustic Unit Recognition PerformanceIn [19], we conducted experiments where the output ofa DPMM sampler served as labels for training a contextdependent “triphone” acoustic unit recognizer. In the proposed
framework, the DPMM defines the number and distributionsof units dynamically and without any prior supervision, given
only extracted speech observations, i.e., frame-based feature
vectors. We trained an HMM acoustic model and n-gram language model using collapsed DPMM component labels. The
resulting DPMM-HMM acoustic unit recognizer was evaluatedby solving the ABX sound class discriminability task. Our
results showed that it is possible to build a DPMM-HMM
acoustic unit recognizer that is competitive with supervisedly
trained phone recognizers.
For this work, we conducted experiments with a comparable setup, but using the DPMoMM for subword modelinginstead. The labels that we used come from the models
that were inferred from PLP+fMLLR-transformed features
as listed in Table I. We trained language-dependent acoustic
unit recognizers given the DPMoMM cluster labels. We thencompared the performance to recognizers that were trained on
the DPMoMM component labels instead, which correspondto standard DPMM labels. Since we sampled each DPMoMM
five times, we also trained each recognizer five times, scoredfive times and averaged the results. In addition to the ABX
test, we also evaluated the NMI scores. Training and testing
was done with 12-fold cross-validation. Table II compares theresults of the trained acoustic unit recognizers using cluster

14

labels vs. using single component labels as training references.
The recognizers trained on cluster labels perform significantly
better across data sets and across all evaluation metrics.
IX. C ONCLUSION
We developed a Dirichlet process mixture of mixturesmodel (DPMoMM) sampler that jointly infers a codebook ofglobal components and a mixture of clusters. The clusters
are themselves mixtures that are defined over the codebook
components. Split and merge samplers for modifying codebook components were complemented with split and merge
samplers for modifying clusters. We introduced an additionalswitch sampler for cluster components to support and accelerate convergence, which has shown to be effective inexperiments on real data. We demonstrated in the use case of
unsupervised subword modeling on two separate data sets, thatclasses represented by a mixture of mixtures model reliably
outperform the output of a standard DPMM. For both data
sets and on all our evaluation metrics, the models inferredwith our proposed DPMoMM sampler consistently achieved
significant performance improvements of up to 13.5% relative.
At the same time, considerably fewer classes that show longer
average durations were sampled, a behavior that is desired
for unsupervised subword modeling. Our experiments suggestthat the inferred mixture of mixtures approximates the true
underlying distributions of our experimental data much better
than a standard DPMM. Lastly, parallelization of sampling
steps allows the algorithm to work well even on larger data
sets in higher data regimes. Although we originally developed
the DPMoMM for improving unsupervised subword modeling,
we think that our sampler will be useful not only for handlingspeech data, but for many tasks where classes are defined bymultimodal distributions.
ACKNOWLEDGMENT
Part of this work was supported by JSPS KAKENHI Grant
Numbers JP17H06101 and JP17K00237.
R EFERENCES[1] T. S. Ferguson, “A Bayesian analysis of some nonparametric problems,”
Annals of Statistics, pp. 209–230, 1973.
[2] C. E. Antoniak, “Mixtures of Dirichlet processes with applications toBayesian nonparametric problems,” Annals of Statistics, pp. 1152–1174,
1974.
[3] C. E. Rasmussen, “The infinite Gaussian mixture model,” in Advances in
Neural Information Processing Systems. Neural Information Processing
Systems Foundation, 2000, pp. 554–560.
[4] Y. W. Teh, M. I. Jordan, M. J. Beal, and D. M. Blei, “Sharing clusters
among related groups: Hierarchical Dirichlet processes,” in Advances in
Neural Information Processing Systems. Neural Information Processing
Systems Foundation, 2005, pp. 1385–1392.
[5] J. Chang and J. W. Fisher III, “Parallel sampling of HDPs using subcluster splits,” in Advances in Neural Information Processing Systems.
Neural Information Processing Systems Foundation, 2014, pp. 235–243.
[6] S. N. MacEachern, “Dependent nonparametric processes,” in Proceedings of the Section on Bayesian Statistical Science. American Statistical
Association, 1999, pp. 50–55.
[7] ——, “Dependent Dirichlet processes,” Unpublished manuscript, Department of Statistics, The Ohio State University, pp. 1–40, 2000.
[8] M. Versteegh, R. Thiolliere, T. Schatz, X. N. Cao, X. Anguera,
A. Jansen, and E. Dupoux, “The zero resource speech challenge 2015,”
in Proceedings of Interspeech. International Speech Communication
Association, 2015, pp. 3169–3173.

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015[9] E. Dunbar, X. N. Cao, J. Benjumea, J. Karadyi, M. Bernard, L. Besacier,X. Anguerra, and E. Dupoux, “The zero resource speech challenge
2017,” in Proceedings of the Automatic Speech Recognition and Understanding Workshop. Institute of Electrical and Electronics Engineers,
2017, pp. 323–330.
[10] H. Chen, C.-C. Leung, L. Xie, B. Ma, and H. Li, “Parallel inferenceof Dirichlet process Gaussian mixture models for unsupervised acoustic
modeling: A feasibility study,” in Proceedings of Interspeech. International Speech Communication Association, 2015, pp. 3189–3193.
[11] ——, “Multilingual bottle-neck feature learning from untranscribedspeech,” in Proceedings of the Automatic Speech Recognition and Understanding Workshop. Institute of Electrical and Electronics Engineers,
2017.
[12] M. Heck, S. Sakti, and S. Nakamura, “Unsupervised linear discriminant analysis for supporting DPGMM clustering in the zero resource
scenario,” in Proceedings of the International Workshop on Spoken
Language Technologies for Under-resourced Languages, 2016, pp. 73–
79.
[13] ——, “Supervised learning of acoustic models in a zero resource
setting to improve DPGMM clustering.” in Proceedings of Interspeech.
International Speech Communication Association, 2016, pp. 1310–1314.
[14] ——, “Feature optimized DPGMM clustering for unsupervised subword
modeling: A contribution to Zerospeech 2017,” in Proceedings of theAutomatic Speech Recognition and Understanding Workshop. Instituteof Electrical and Electronics Engineers, 2017, pp. 740–746.
[15] L. R. Rabiner and B.-H. Juang, Fundamentals of speech recognition.
Prentice Hall, 1993.
[16] R. Schwartz, Y. Chow, O. Kimball, S. Roucos, M. Krasner, andJ. Makhoul, “Context-dependent modeling for acoustic-phonetic recognition of continuous speech,” in Proceedings of the International Conference on Acoustics, Speech and Signal Processing, vol. 10. Instituteof Electrical and Electronics Engineers, 1985, pp. 1205–1208.
[17] K.-F. Lee, “On large-vocabulary speaker-independent continuous speech
recognition,” Speech Communication, vol. 7, no. 4, pp. 375–379, 1988.
[18] A. Waibel and K.-F. Lee, Eds., Readings in speech recognition. Morgan
Kaufmann, 1990.
[19] M. Heck, S. Sakti, and S. Nakamura, “Iterative training of a DPGMMHMM acoustic unit recognizer in a zero resource scenario,” in Proceedings of the Spoken Language Technology Workshop. Institute ofElectrical and Electronics Engineers, 2016, pp. 57–63.
[20] H. Z. Yerebakan, B. Rajwa, and M. Dundar, “The infinite mixture of
infinite Gaussian mixtures,” in Advances in Neural Information Processing Systems. Neural Information Processing Systems Foundation, 2014,
pp. 28–36.
[21] J. Chang and J. W. Fisher III, “Parallel sampling of DP mixture models
using sub-cluster splits,” in Advances in Neural Information Processing
Systems. Neural Information Processing Systems Foundation, 2013,
pp. 620–628.
[22] X. Zhu, Z. Ghahramani, and J. Lafferty, “Time-sensitive Dirichlet
process mixture models,” School of Computer Science, Carnegie Mellon
University, Tech. Rep., 2005.
[23] J. E. Griffin and M. J. Steel, “Order-based dependent Dirichlet processes,” Journal of the American Statistical Association, vol. 101, no.
473, pp. 179–194, 2006.
[24] D. M. Blei and P. I. Frazier, “Distance dependent Chinese restaurant
processes,” Journal of Machine Learning Research, vol. 12, pp. 2461–
2488, 2011.
[25] D. B. Dahl, “An improved merge-split sampler for conjugate Dirichlet
process mixture models,” Department of Statistics, University of Wisconsin Madison, Tech. Rep., 2003.
[26] S. Jain and R. M. Neal, “A split-merge Markov chain Monte Carlo procedure for the Dirichlet process mixture model,” Journal of Computational
and Graphical Statistics, vol. 13, no. 1, pp. 158–182, 2004.
[27] ——, “Splitting and merging components of a nonconjugate Dirichlet
process mixture model,” Bayesian Analysis, vol. 2, no. 3, pp. 445–472,
2007.
[28] O. Papaspiliopoulos and G. O. Roberts, “Retrospective Markov chain
Monte Carlo methods for Dirichlet process hierarchical models,”
Biometrika, vol. 95, no. 1, pp. 169–186, 2008.
[29] S. Favaro and Y. W. Teh, “MCMC for normalized random measuremixture models,” Statistical Science, pp. 335–359, 2013.
[30] H. Ishwaran and L. F. James, “Gibbs sampling methods for stickbreaking priors,” Journal of the American Statistical Association, vol. 96,
no. 453, pp. 161–173, 2001.
[31] J. Yang, R. C. Van Dalen, and M. J. Gales, “Infinite support vector
machines in speech recognition,” in Proceedings of Interspeech. International Speech Communication Association, 2013.

15[32] J. Yang, R. C. Van Dalen, S.-X. Zhang, and M. J. Gales, “Infinite
structured support vector machines for speech recognition,” in Proceedings of the International Conference on Acoustics, Speech and Signal
Processing. Institute of Electrical and Electronics Engineers, 2014, pp.
3320–3324.
[33] T. Schatz, V. Peddinti, F. Bach, A. Jansen, H. Hermansky, and
E. Dupoux, “Evaluating speech features with the minimal-pair ABX
task: Analysis of the classical MFC/PLP pipeline,” in Proceedings ofInterspeech. International Speech Communication Association, 2013,
pp. 1781–1785.
[34] A. F. McDaid, D. Greene, and N. Hurley, “Normalized mutual information to evaluate overlapping community finding algorithms,” arXiv
preprint arXiv:1110.2515, 2011.
[35] M. A. Pitt, K. Johnson, E. Hume, S. Kiesling, and W. Raymond, “The
Buckeye corpus of conversational speech: Labeling conventions and atest of transcriber reliability,” Speech Communication, vol. 45, no. 1, pp.
89–95, 2005.
[36] N. J. De Vries, M. H. Davel, J. Badenhorst, W. D. Basson, F. De Wet,
E. Barnard, and A. De Waal, “A smartphone-based ASR data collection
tool for under-resourced languages,” Speech Communication, vol. 56,
pp. 119–131, 2014.
[37] D. Povey, A. Ghoshal, G. Boulianne, L. Burget, O. Glembek, N. Goel,
M. Hannemann, P. Motlicek, Y. Qian, P. Schwarz et al., “The Kaldi
speech recognition toolkit,” in Proceedings of the Automatic Speech
Recognition and Understanding Workshop. Institute of Electrical andElectronics Engineers, 2011.
[38] H. Hermansky, “Perceptual linear predictive (PLP) analysis of speech,”
The Journal of the Acoustical Society of America, vol. 87, no. 4, pp.
1738–1752, 1990.
[39] R. Haeb-Umbach and H. Ney, “Linear discriminant analysis for improved large vocabulary continuous speech recognition,” in Proceedingsof the International Conference on Acoustics, Speech and Signal Processing, vol. 1. Institute of Electrical and Electronics Engineers, 1992,
pp. 13–16.
[40] R. A. Gopinath, “Maximum likelihood modeling with Gaussian distributions for classification,” in Proceedings of the International Conferenceon Acoustics, Speech and Signal Processing. Institute of Electrical andElectronics Engineers, 1998, pp. 661–664.
[41] M. J. Gales, “Semi-tied covariance matrices for hidden Markov models,”IEEE Transactions on Speech and Audio Processing, vol. 7, no. 3, pp.
272–281, 1999.
[42] T. Anastasakos, J. McDonough, R. Schwartz, and J. Makhoul, “A
compact model for speaker-adaptive training,” in Proceedings of theInternational Conference on Spoken Language. Institute of Electricaland Electronics Engineers, 1996, pp. 1137–1140.
[43] M. J. Gales, “Maximum likelihood linear transformations for HMMbased speech recognition,” Computer Speech and Language, vol. 12,
no. 2, pp. 75–98, 1998.

Michael Heck received his diploma degree in computer science from Karlsruhe Institute of Technology (KIT), Germany, in 2012. He was granted
the “Baden-Württemberg” scholarship and was supported by the Global Initiatives Program for a stay at
Nara Institute of Science and Technology (NAIST),
Japan, during his thesis work. Between 2013-2015,
he worked as research assistant at the Interactive
Systems Labs at KIT on unsupervised training andadaptation of acoustic models for automatic speech
recognition. He was involved in several speech
recognition evaluations and international programs including Quaero, Babeland EU-BRIDGE. While at KIT, he was a regular visiting researcher at theAugmented Human Communication Laboratory (AHC Lab) at NAIST. Since
2015, he is a research assistant and doctoral course student at AHC Lab asrecipient of the NAIST international scholarship. In 2017, he was a researchintern at IBM Research – Tokyo. His research interests include automatic
speech recognition, unsupervised learning and the zero resource scenario.

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at
http://dx.doi.org/10.1109/TASLP.2018.2852500
JOURNAL OF LATEX CLASS FILES, VOL. 14, NO. 8, AUGUST 2015

Sakriani Sakti received her B.E. degree in Informatics (cum laude) from Bandung Institute of Technology, Indonesia, in 1999. In 2000, she received
DAAD-Siemens Program Asia 21st Century Award
to study in Communication Technology, University
of Ulm, Germany, and received her MSc degreein 2002. During her thesis work, she worked with
Speech Understanding Department, DaimlerChrysler
Research Center, Ulm, Germany. Between 20032009, she worked as a researcher at ATR SLC Labs,
Japan, and during 2006-2011, she worked as anexpert researcher at NICT SLC Groups, Japan. While working with ATRNICT, Japan, she continued her study (2005-2008) with Dialog Systems Group
University of Ulm, Germany, and received her PhD degree in 2008. She wasactively involved in collaboration activities such as Asian Pacific Telecommunity Project (2003-2007), A-STAR and U-STAR (2006-2011). In 2009-2011,
she served as a visiting professor of Computer Science Department, Universityof Indonesia (UI), Indonesia. From 2011, she has been an assistant professor
at the Augmented Human Communication Laboratory, NAIST, Japan. She
served also as a visiting scientific researcher of INRIA Paris-Rocquencourt,
France, in 2015-2016, under ”JSPS Strategic Young Researcher Overseas
Visits Program for Accelerating Brain Circulation”. She is a member of JNS,
SFN, ASJ, ISCA, IEICE and IEEE. Her research interests include statistical
pattern recognition, speech recognition, spoken language translation, cognitive
communication, and graphical modeling framework.

Satoshi Nakamura is Professor at the Graduate
School of Information Science, Nara Institute of
Science and Technology, Japan, Honorarprofessorof Karlsruhe Institute of Technology, Germany, andATR Fellow. He received his B.S. from Kyoto
Institute of Technology in 1981 and Ph.D. from
Kyoto University in 1992. He was Associate Professor of Graduate School of Information Science at
Nara Institute of Science and Technology in 19942000. He was Director of ATR Spoken Language
Communication Research Laboratories in 2000-2008
and Vice president of ATR in 2007-2008. He was Director General of
Keihanna Research Laboratories and the Executive Director of Knowledge
Creating Communication Research Center, National Institute of Information
and Communications Technology, Japan in 2009-2010. He is currently Director of Augmented Human Communication laboratory and a full professorof Graduate School of Information Science at Nara Institute of Science andTechnology. He is interested in modeling and systems of speech-to-speechtranslation and speech recognition. He is one of the leaders of speech-tospeech translation research and has been serving for various speech-to-speechtranslation research projects in the world including C-STAR, IWSLT and
A-STAR. He received Yamashita Research Award, Kiyasu Award from theInformation Processing Society of Japan, Telecom System Award, AAMT
Nagao Award, Docomo Mobile Science Award in 2007, ASJ Award forDistinguished Achievements in Acoustics. He received the Commendationfor Science and Technology by the Minister of Education, Science andTechnology, and the Commendation for Science and Technology by the
Minister of Internal Affairs and Communications. He also received LREC
Antonio Zampoli Award 2012. He has been Elected Board Member ofInternational Speech Communication Association, ISCA, since June 2011,
IEEE Signal Processing Magazine Editorial Board Member since April 2012,
IEEE SPS Speech and Language Technical Committee Member since 2013,
and IEEE Fellow since 2016.

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

16