# Bayesian Clustering via Fusing of Localized Densities

Alexander Dombowsky  Department of Statistical Science, Duke University,
Durham, NC, USA  David B. Dunson  Department of Statistical Science, Duke
University, Durham, NC, USA  Department of Mathematics, Duke University,
Durham, NC, USA

###### Abstract

Bayesian clustering typically relies on mixture models, with each component
interpreted as a different cluster. After defining a prior for the component
parameters and weights, Markov chain Monte Carlo (MCMC) algorithms are
commonly used to produce samples from the posterior distribution of the
component labels. The data are then clustered by minimizing the expectation of
a clustering loss function that favors similarity to the component labels.
Unfortunately, although these approaches are routinely implemented, clustering
results are highly sensitive to kernel misspecification. For example, if
Gaussian kernels are used but the true density of data within a cluster is
even slightly non-Gaussian, then clusters will be broken into multiple
Gaussian components. To address this problem, we develop Fusing of Localized
Densities (FOLD), a novel clustering method that melds components together
using the posterior of the kernels. FOLD has a fully Bayesian decision
theoretic justification, naturally leads to uncertainty quantification, can be
easily implemented as an add-on to MCMC algorithms for mixtures, and favors a
small number of distinct clusters. We provide theoretical support for FOLD
including clustering optimality under kernel misspecification. In simulated
experiments and real data, FOLD outperforms competitors by minimizing the
number of clusters while inferring meaningful group structure.

Keywords: Mixture model; Loss function; Decision theory; Kernel
misspecification; Statistical distance

##  1 Introduction

Clustering data into groups of relatively similar observations is a canonical
task in exploratory data analysis. Algorithmic clustering methods, such as
k-means, k-medoids, and hierarchical clustering, rely on dissimilarity
metrics, an approach which is often heuristic but may perform well in
practice; see Hastie et al., (2009), Jain, (2010), and Kiselev et al., (2019)
for an overview. In comparison, model-based clustering methods utilize
mixtures of probability kernels to cluster data, ordinarily by inferring the
component labels (Fraley and Raftery,, 2002). Choices of kernel depend on the
type of data being considered, with the Gaussian mixture model (GMM)
particularly popular for continuous data. A conceptual advantage of model-
based methods is the ability to express uncertainty in clustering. For
example, the expectation-maximization (EM) algorithm uses maximum likelihood
estimates of the component-specific parameters and weights to calculate each
observations’ posterior component allocation probabilities, which are
interpreted as a measure of clustering uncertainty (Bensmail et al.,, 1997).
As estimation of the weights and component-specific parameters is crucial for
model-based clustering, there is a rich literature on quantifying rates of
convergence for various mixture models, which is usually expressed in terms of
the mixing measure (see Guha et al.,, 2021 and references therein).

Bayesian mixture models have received increased attention as a clustering
method in recent years. The Bayesian framework can account for uncertainty in
the number of components and incorporate prior information on the component-
specific parameters, with Markov chain Monte Carlo (MCMC) algorithms employed
to generate posterior samples for the mixture weights, kernel parameters, and
component labels for each data point. Based on posterior samples of the
component labels, one can obtain Monte Carlo approximations to Bayes
clustering estimators. For $p$-dimensional data $\bm{X}=(X_{1},\dots,X_{n})$,
a Bayes estimator $\bm{c}^{*}$ is the minimizer of an expected loss function
conditional on $\bm{X}$ over all possible clusterings
$\bm{c}=(c_{1},\dots,c_{n})$. There are several popular choices of loss,
including Binder’s (Binder,, 1978) and the Variation of Information (VI)
(Meilă,, 2007), which are distance metrics over the space of partitions and
favor clusterings that are similar to the component labels. Further details on
set partitions and their role in Bayesian clustering can be found in Meilă,
(2007), Wade and Ghahramani, (2018), and Paganin et al., (2021). Several
measures of uncertainty in clustering with Bayesian mixtures exist, including
the posterior similarity matrix and credible balls of partitions (Wade and
Ghahramani,, 2018).

Figure 1: A Bayesian GMM with $10$ components and Dirichlet prior
concentration parameter equal to $1/10$ is fit to data generated from a
mixture of skew Gaussian distributions.

However, model-based clustering approaches, including Bayesian
implementations, can be brittle in applications due to unavoidable impacts of
model misspecification. While the data are assumed to come from a model
$X_{i}\sim f=\sum_{k=1}^{K}a_{k}g(\cdot;\tilde{\theta}_{k})$, in reality, data
are generated via a true process $f^{0}$, which may or may not be in the
support of the prior for $f$. A specific example is kernel misspecification,
which occurs when $f^{0}$ is in fact a mixture model whose components are not
contained in the kernel family mixed over by the fitted model. Mixtures will
often induce over-clustering in this setting by approximating a single
component in $f^{0}$ with many components from $f$. An example of this
phenomena is shown in Figure 1, where a $10$ component Bayesian GMM is fit to
data generated from a mixture of bivariate skew Gaussian distributions.
Despite using a concentration parameter of $1/10$ in the symmetric Dirichlet
prior to induce a small number of clusters (Rousseau and Mengersen,, 2011),
the GMM allocates the data into $10$ poorly defined groups. A group of
observations in the bottom third of the sample space is split into two, while
several observations are placed into their own groups despite being near dense
collections of data.

Several approaches have been proposed to address the issue of over-clustering
due to kernel misspecification. A natural solution is to define a flexible
class of kernels, exemplified by the mixtures in Karlis and Santourian,
(2009), O’Hagan et al., (2016), and Dang et al., (2023). To increase
flexibility further, Rodríguez and Walker, (2014) propose a mixture of
nonparametric unimodal kernels. Similarly, Bartolucci, (2005), Li, (2005), Di
Zio et al., (2007), and Malsiner-Walli et al., (2017) use carefully chosen
mixtures of Gaussians to characterize the data within each cluster. However,
there is an unfortunate pitfall with the general strategy of using flexible
families of kernels. In particular, as the flexibility of the kernel
increases, identifiability and optimal estimation rates for inferring the
mixing measure tend to weaken, especially when the true number of mixture
components is unknown (Nguyen,, 2013; Ho and Nguyen, 2016a, ; Heinrich and
Kahn,, 2018). Even the transition from location Gaussian kernels with known
covariance to location-scale Gaussian kernels can have substantial
consequences on the convergence rate of the component means (Manole and Ho,,
2022). Such problems motivated Ho et al., (2020) to propose an alternative
estimator of the mixing measure that is more robust than maximum likelihood
and Bayesian approaches and can achieve optimal convergence rates.

Alternatively, one can develop generalized Bayesian methods of clustering that
avoid defining a fully generative probabilistic model for the data. For
example, Duan and Dunson, (2021) propose to conduct model-based clustering
based on a pairwise distance matrix instead of the data directly to reduce
sensitivity to kernel misspecification. Rigon et al., (2023) instead define a
Gibbs posterior for clustering, incorporating a clustering loss function in
place of a likelihood function, completely bypassing modeling of the data. An
alternative that maintains a fully generative model while robustifying
inferences to misspecification is to use a coarsened posterior (Miller and
Dunson,, 2019; Gorsky et al.,, 2024). This coarsening approach often has good
practical performance in reducing over-clustering due to kernel
misspecification.

Unfortunately, these approaches for accommodating kernel misspecification have
various drawbacks. Along with slower convergence rates, flexible kernels
typically require a large number of parameters, worsening the already
burdensome computational cost of Bayesian clustering. The generalized Bayes
approaches can perform well in certain settings. However, both Gibbs
posteriors and coarsened posteriors are highly sensitive to key tuning
parameters, which can be difficult to choose objectively in practice.

A key theoretical advance is given in Aragam et al., (2020), which attempts to
solve the problem of clustering based on a mixture model by merging components
in a Gaussian mixture. They rely on a two-stage procedure that lacks partition
uncertainty quantification and assumes that the true number of kernels is
known. However, the approach of viewing clusters as arising from merging
closely overlapping kernels is promising. Related merging ideas have been
implemented in both frequentist and Bayesian approaches in a variety of
settings, and several algorithms exist for deciding how and when to combine
components together (Chan et al.,, 2008; Baudry et al.,, 2010; Hennig,, 2010;
Melnykov,, 2016; Guha et al.,, 2021; Manole and Khalili,, 2021).

In this paper, we propose a novel decision theoretic method for Bayesian
clustering that mitigates the effects of model misspecification. Suppose we
model the data with a Bayesian mixture model, with $K$ components, component
labels $\bm{s}=(s_{1},\dots,s_{n})$, component-specific atoms
$\tilde{\theta}_{1},\dots,\tilde{\theta}_{K}$, and kernels
$g(\cdot;\tilde{\theta}_{1}),\dots,g(\cdot;\tilde{\theta}_{K})$. Rather than
focusing on $\bm{s}$, we compute clusters with the localized densities
$\left\\{g(\cdot;\theta_{i})\right\\}_{i=1}^{n}$, defined by
$\theta_{i}=\tilde{\theta}_{s_{i}}$. We define a loss function for any
clustering $\widehat{\bm{c}}$ that favors allocating $i$ and $j$ to the same
cluster when the statistical distance between $g(\cdot;\theta_{i})$ and
$g(\cdot;\theta_{j})$ is small, encouraging grouping observations with
overlapping component kernels. We cluster the data with a Bayes estimator,
interpreted as a Fusing of Localized Densities (FOLD). Our method has a fully
decision theoretic justification, leads to interpretable uncertainty
quantification, and can be readily implemented using the output of existing
MCMC algorithms for mixtures.

Though previous methods have utilized merging kernels to account for kernel
misspecification, to our knowledge none have a formal Bayesian decision
theoretic justification. Although FOLD requires combinatorial optimization
when obtaining the point estimate, we suggest two reasonable approximations:
one in which we reduce the optimization space to a smaller set of candidate
clusterings, and another in which we use existing greedy algorithms to compute
a locally-optimal solution. We also provide concrete theoretical guarantees
for our procedure, including the result that as the sample size grows, the
clustering point estimate converges to an oracle clustering rule. In addition,
we list intuitive examples of oracle rules and empirically validate our
asymptotic theory using a simple example.

In Section 2, we explain our clustering method from a Bayesian decision
theoretic perspective, provide a framework for uncertainty quantification, and
demonstrate how to implement FOLD in practice. In Section 3, we show
asymptotic concentration of FOLD for misspecified and well-specified regimes.
We apply FOLD to synthetic data and a cell line data set in Section 4, showing
excellent performance in both cases. Finally, we provide concluding remarks
and some extensions in Section 5.

##  2 Clustering with Localized Densities

###  2.1 Notation and Setup

Let $X_{i}^{0}=(X_{i1}^{0},\dots,X_{ip}^{0})\in\mathbb{R}^{p}$ be multivariate
observations, collected into $\bm{X}^{0}=(X_{1}^{0},\dots,X_{n}^{0})$. Assume
that $X_{1}^{0},\dots,X_{n}^{0}$ are random and generated from an unknown
mixture model: $X_{i}^{0}\sim f^{0}$, where
$f^{0}=\sum_{m=1}^{M^{0}}a_{m}^{0}g_{m}^{0}(\cdot)$, $a_{m}^{0}>0$ for
$m=1,\dots,M^{0}$, $\sum_{m=1}^{M^{0}}a_{m}^{0}=1$, and $M^{0}<\infty$. Since
$f^{0}$ is a mixture, the data-generating process can be equivalently stated
with the addition of latent variables
$s_{i}^{0}\in\left\\{1,\dots,M^{0}\right\\}$, so that $\left(X_{i}^{0}\mid
s_{i}^{0}=m\right)\sim g_{m}^{0}(\cdot)$ for all $i=1,\dots,n$.
$\mathbb{P}^{0}$ refers to probability and convergence statements with respect
to the true data generating process, with
$\mathbb{P}^{0}(A)=\int_{A}f^{0}(x)dx$ for any set $A$.

Throughout this article, we will focus on the setting in which $\bm{X}^{0}$ is
modeled with a mixture. Let $X_{i}=(X_{i1},\dots,X_{ip})\in\mathbb{R}^{p}$,
where $\bm{X}=(X_{1},\dots,X_{n})$. A Bayesian mixture model assumes that
$\bm{X}$ are generated via

$\displaystyle\lambda$ $\displaystyle\sim\pi_{\Lambda},$
$\displaystyle\theta_{i}\mid\lambda$ $\displaystyle\sim\lambda,$
$\displaystyle X_{i}\mid\theta_{i}\sim g(X_{i};\theta_{i});$ (1)

for all $i=1,\dots,n$. Here, $\lambda$ is a probability distribution over a
parameter space $\Theta$ and known as the mixing measure, whereas
$\mathcal{G}=\left\\{g(\cdot;\theta):\theta\in\Theta\right\\}$ denotes a
family of probability densities with support on $\mathbb{R}^{p}$. We will
impose more conditions on both $\Theta$ and $\mathcal{G}$ when we develop our
theoretical results in Section 3. $\pi_{\Lambda}$ is the prior distribution
for the mixing measure with support $\Lambda$, and examples of such priors are
implemented in Bayesian finite mixture models, mixtures of finite mixtures
(MFMs) (Miller and Harrison,, 2018), Dirichlet process mixtures (DPMs)
(Ferguson,, 1973; Antoniak,, 1974), and Pitman-Yor process mixtures (Ishwaran
and James,, 2001, 2003). Generally, realizations from $\pi_{\Lambda}$ are
almost-surely discrete (e.g., as in Sethuraman, (1994)) and can be represented
as $\lambda=\sum_{k=1}^{K}a_{k}\delta_{\tilde{\theta}_{k}}$, where $1\leq
K\leq\infty$, $\tilde{\theta}_{k}\in\Theta$, $0<a_{k}<1$, and
$\sum_{k=1}^{K}a_{k}=1$. Hence, sampling $\theta_{i}$ from (1) results in ties
amongst the atoms $\bm{\theta}=(\theta_{1},\dots,\theta_{n})$, which are
represented by unique values
$\tilde{\bm{\theta}}=(\tilde{\theta}_{1},\dots,\tilde{\theta}_{K_{n}})$ for
some $1\leq K_{n}\leq n$. From these ties we can construct a latent partition
$S=(S_{1},\dots,S_{K_{n}})$ of the integers $[n]=\left\\{1,\dots,n\right\\}$,
where $S_{k}=\left\\{i\in[n]:\theta_{i}=\tilde{\theta}_{k}\right\\}$. The
latent partition may also be represented with labels
$\bm{s}=(s_{1},\dots,s_{n})$, where $s_{i}=k\iff i\in S_{k}$. We will assume
that the prior for $\tilde{\bm{\theta}}$ is non-atomic, so the allocation of
two observations to the same component is equivalent to equality in their
atoms, i.e. $s_{i}=s_{j}\iff\theta_{i}=\theta_{j}$.

Ultimately, our goal is to cluster the data $\bm{X}^{0}$ using the mixture
model outlined in (1). The typical approach is to infer the conditional
distribution of $S$ after observing the data. More formally, by assuming that
$\lambda\sim\pi_{\Lambda}$, we have induced a prior distribution for $S$,
$\Pi(S)$. We then cluster the data by inferring $\Pi(S\mid\bm{X}=\bm{X}^{0})$,
which we will refer to for the remainder of this article as
$\Pi(S\mid\bm{X}^{0})$. To construct a point estimate of $S$, we minimize the
posterior expectation of a clustering loss function that favors similarity to
$S$. The decision theoretic approach to Bayesian clustering was first proposed
in Binder, (1978), and was advocated for more recently in Wade and Ghahramani,
(2018) and Dahl et al., (2022). However, Bayesian methods can lead to
homogenous, tiny, and often uninterpretable clusters in practice. These issues
often arise in the presence of kernel misspecification. Even well-separated
true clusters in $\bm{X}^{0}$ will be broken into multiple components from
$\mathcal{G}$ under a minor degree of kernel misspecification. To address this
issue, we propose a loss function focused on simplifying the overfitted
kernels via merging rather than estimating the component labels.

###  2.2 Decision Theory Formulation

We aim to estimate a clustering
$\widehat{C}=\left\\{\hat{C}_{1},\dots,\hat{C}_{\hat{K}_{n}}\right\\}$ of
$\bm{X}^{0}$, or labels $\widehat{\bm{c}}=(\hat{c}_{1},\dots,\hat{c}_{n})$,
based on merging components from the joint posterior that have similar kernels
$g(\cdot;\tilde{\theta}_{k})$. The motivation here is to remedy the cluster
splitting that occurs in using multiple parametric kernels
$g(\cdot;\tilde{\theta}_{k})$ to represent each “true" kernel
$g_{m}^{0}(\cdot)$. We define the localized density of observation $i$ to be
$g(\cdot;\theta_{i})$, a random probability density that takes values
$g(x;\theta_{i})$ for all $x\in\mathbb{R}^{p}$, where $\theta_{i}$ is the atom
of observation $i$ as defined in (1). Observe that the localized density is
the distribution for $X_{i}$ under the fitted mixture model given
$\theta_{i}$, i.e. $(X_{i}\mid\theta_{i}=\tilde{\theta}_{k})\sim
g(\cdot;\tilde{\theta}_{k})$. Prior variation in the localized densities is
modeled by sampling $\lambda\sim\pi_{\lambda}$, and then sampling
$\theta_{i}\sim\lambda$, for $i=1,\ldots,n$.

Figure 2: A Bayesian location Gaussian mixture is fit to a version of the
moons data set. The localized densities are inferred, then merged together to
recover the crescents.

To counteract cluster splitting, we define the loss of assigning two
observations into the same or different clusters as a function of the
statistical distance between their localized densities. Figure 2 shows the
behavior of the posterior distributions of the localized densities when the
model $f(x)=\sum_{k=1}^{30}a_{k}\mathcal{N}_{2}(x;\tilde{\theta}_{k},0.02I)$
is fit to a version of the moons data. The red points plot
$\mathbb{E}_{\Pi}(\theta_{i}\mid\bm{X}^{0})$ for each $i=1,\dots,n$, with gray
circles representing the $95\%$ high density regions of a
$\mathcal{N}_{2}(\mathbb{E}_{\Pi}(\theta_{i}\mid\bm{X}^{0}),0.02I)$
distribution. The crescent clusters are split into multiple overlapping
Gaussian kernels, but are recovered by fusing these kernels together.

We define the loss of any clustering $\widehat{\bm{c}}$ resulting from the
model in (1) to be

$\displaystyle\mathcal{L}(\widehat{\bm{c}},\bm{\theta})=\sum_{i<j}\left\\{\textbf{1}_{\hat{c}_{i}=\hat{c}_{j}}\mathcal{D}_{ij}+\omega\textbf{1}_{\hat{c}_{i}\neq\hat{c}_{j}}\left(1-\mathcal{D}_{ij}\right)\right\\},$
(2)

where $\omega>0$,
$\mathcal{D}_{ij}=d\left\\{g(\cdot;\theta_{i}),g(\cdot;\theta_{j})\right\\}$,
and $d$ is any distance between probability distributions assumed to be
bounded in the unit interval, i.e. $0\leq
d\left\\{P(\cdot),Q(\cdot)\right\\}\leq 1$ for all measures $P$ and $Q$, and
with the property that $d\left\\{P(\cdot),Q(\cdot)\right\\}=1$ if and only if
$\textrm{supp}(P)\cap\textrm{supp}(Q)=\emptyset$. A key example is when $d$ is
chosen to be the Hellinger distance, but this need not be the case in general.
Observe that $\mathcal{L}(\widehat{\bm{c}},\bm{\theta})$ is non-negative for
any $\widehat{\bm{c}}$. The loss of assigning $\hat{c}_{i}=\hat{c}_{j}$ is
$\mathcal{D}_{ij}$. If $s_{i}=s_{j}$, then $\theta_{i}=\theta_{j}$ and hence
there is no loss incurred. When $s_{i}\neq s_{j}$, then
$\theta_{i}\neq\theta_{j}$ but the loss will remain small when the kernels
$g(\cdot;\theta_{i}),g(\cdot;\theta_{j})$ are similar under $d$. Conversely,
allocating $i$ and $j$ to different clusters results in a loss of
$\omega(1-\mathcal{D}_{ij})$. If $\mathcal{D}_{ij}=1$, i.e. if the supports of
$g(\cdot;\theta_{i})$ and $g(\cdot;\theta_{j})$ are disjoint, then setting
$\hat{c}_{i}\neq\hat{c_{j}}$ accumulates zero loss. Otherwise, the loss
depends on the separation in $d$ between the localized densities and the loss
parameter $\omega$.

Our loss $\mathcal{L}(\widehat{\bm{c}},\bm{\theta})$ is invariant to
permutations of the data indices and of the labels in either
$\widehat{\bm{c}}$ or $\bm{s}$, which is desirable for clustering losses
(Binder,, 1978). In addition, $\mathcal{L}(\widehat{\bm{c}},\bm{\theta})$ is a
continuous relaxation of Binder’s loss function (Binder,, 1978; Lau and
Green,, 2007),

$\mathcal{L}_{\textrm{B}}(\widehat{\bm{c}},\bm{s})=\sum_{i<j}\left\\{\textbf{1}_{\hat{c}_{i}=\hat{c}_{j}}\textbf{1}_{s_{i}\neq
s_{j}}+\omega\textbf{1}_{\hat{c}_{i}\neq\hat{c}_{j}}\textbf{1}_{s_{i}=s_{j}}\right\\}.$
(3)

A key property of Binder’s loss is that it is a quasimetric over the space of
partitions (Wade and Ghahramani,, 2018; Dahl et al.,, 2022).
$\mathcal{L}(\widehat{\bm{c}},\bm{\theta})$ is not a quasimetric, but instead
aims to mitigate the problem of cluster splitting arising under kernel
misspecification. Note the difference in notation in the arguments of (3) and
(2): we omit $\bm{s}$ from the notion for
$\mathcal{L}(\widehat{\bm{c}},\bm{\theta})$ to reflect that we will
incorporate information from the atoms in clustering $\bm{X}^{0}$, not just
the latent partition. However, we can also state (3) in terms of $\bm{\theta}$
by recalling that $s_{i}=s_{j}\iff\theta_{i}=\theta_{j}$. One can show that
$\mathcal{L}(\widehat{\bm{c}},\bm{\theta})$ can be rewritten as the sum of
$\mathcal{L}_{\textrm{B}}(\widehat{\bm{c}},\bm{s})$ and a remainder term
$\mathcal{B}(\widehat{\bm{c}},\tilde{\bm{\theta}})$ that depends on the unique
values $\tilde{\bm{\theta}}$, and that $\mathcal{L}(\bm{s},\bm{\theta})=0$
only when the components of $f$ are completely separated under $d$. The closed
form of the remainder $\mathcal{B}(\widehat{\bm{c}},\tilde{\bm{\theta}})$ is
given in the Supplement, and is interpreted as an added cost to Binder’s loss
that encourages merging components.

######  Proposition 1.

$\mathcal{L}(\widehat{\bm{c}},\bm{\theta})=\mathcal{L}_{\textrm{B}}(\widehat{\bm{c}},\bm{s})+\mathcal{B}(\widehat{\bm{c}},\tilde{\bm{\theta}}),$
where $\mathcal{B}(\bm{s},\tilde{\bm{\theta}})=0$ if and only if
$d\\{g(\cdot;\tilde{\theta}_{k}),g(\cdot;\tilde{\theta}_{k^{\prime}})\\}=1$
for all component pairs $k\neq k^{\prime}=1,\dots,K$.

The parameter $\omega$ calibrates separation of the clusters. For example,
suppose we compare the clustering $\widehat{\bm{c}}_{1}$, which includes
clusters $\hat{C}_{h}$ and $\hat{C}_{h^{\prime}}$, with the clustering
$\widehat{\bm{c}}_{2}$, which is equivalent to $\widehat{\bm{c}}_{1}$ but now
contains the merged cluster $\hat{C}_{h}\cup\hat{C}_{h^{\prime}}$. The
difference in their losses is
$\mathcal{L}(\widehat{\bm{c}}_{1},\bm{\theta})-\mathcal{L}(\widehat{\bm{c}}_{2},\bm{\theta})=\omega\sum_{i\in\hat{C}_{h},j\in\hat{C}_{h^{\prime}}}\left(1-\mathcal{D}_{ij}\right)-\sum_{i\in\hat{C}_{h},j\in\hat{C}_{h^{\prime}}}\mathcal{D}_{ij}.$
The loss of $\widehat{\bm{c}}_{2}$ is less than that of $\widehat{\bm{c}}_{1}$
when

$\frac{1}{|\hat{C}_{h}||\hat{C}_{h^{\prime}}|}\sum_{i\in\hat{C}_{h},j\in\hat{C}_{h^{\prime}}}\mathcal{D}_{ij}<\gamma:=\frac{\omega}{1+\omega}.$
(4)

This implies that large values of $\omega$ promote fusing clusters, while
smaller values lead to more clusters having less between-cluster
heterogeneity. When $\widehat{\bm{c}}_{1}=\bm{s}$, it is clear that a smaller
loss can be attained by merging components with average pairwise statistical
distance less than $\gamma$. Trivial clusterings are favored when $\omega$ is
taken to its lower and upper limits, similar to both Binder’s loss and the VI
loss (Wade and Ghahramani,, 2018; Dahl et al.,, 2022). As $\omega\to 0$, our
loss is minimized by placing each observation in its own cluster, and, as
$\omega\to\infty$, all observations are placed in a single cluster. Binder’s
loss exhibits a similar interpretation, and one can show that, under the same
settings outlined above,
$\mathcal{L}_{\textrm{B}}(\widehat{\bm{c}}_{2},\bm{s})<\mathcal{L}_{\textrm{B}}(\widehat{\bm{c}}_{1},\bm{s})$
when
$|\hat{C}_{h}|^{-1}|\hat{C}_{h^{\prime}}|^{-1}\sum_{i\in\hat{C}_{h},j\in\hat{C}_{h^{\prime}}}\textbf{1}_{s_{i}\neq
s_{j}}=|\hat{C}_{h}|^{-1}|\hat{C}_{h^{\prime}}|^{-1}m_{hh^{\prime}}<\gamma$,
where $m_{hh^{\prime}}$ counts the number of pairs across clusters $h$ and
$h^{\prime}$ that are allocated to different mixture components. However, the
loss is always minimized when $\widehat{\bm{c}}=\bm{s}$.

The risk of any clustering $\widehat{\bm{c}}$ is the posterior expectation of
its loss, which integrates over the uncertainty in $\bm{\theta}$ after
observing the data;

$\displaystyle\mathcal{R}(\widehat{\bm{c}})=\mathbb{E}_{\Pi}\left\\{\mathcal{L}(\widehat{\bm{c}},\bm{\theta})\mid\bm{X}^{0}\right\\}=\sum_{i<j}\left\\{\textbf{1}_{\hat{c}_{i}=\hat{c}_{j}}\Delta_{ij}+\omega\textbf{1}_{\hat{c}_{i}\neq\hat{c}_{j}}\left(1-\Delta_{ij}\right)\right\\},$
(5)

where $\Delta_{ij}=\mathbb{E}_{\Pi}(\mathcal{D}_{ij}\mid\bm{X}^{0})$ is the
posterior expected distance between the localized densities
$g(\cdot;\theta_{i})$ and $g(\cdot;\theta_{j})$ for observations $i$ and $j$,
respectively. The point estimator of the clustering of $\bm{X}^{0}$ is given
by a Bayes estimator, denoted $\bm{c}_{\textrm{FOLD}}$, which is the minimizer
of (5) over the space of all possible cluster allocations, i.e.
$\bm{c}_{\textrm{FOLD}}=\underset{\widehat{\bm{c}}}{\textrm{argmin
}}\mathcal{R}(\widehat{\bm{c}}).$ The expected statistical distance terms in
the risk (5) exhibit several bounded and metric properties.

######  Proposition 2.

Let $\Delta=(\Delta_{ij})_{i,j}$, where
$\Delta_{ij}=\mathbb{E}_{\Pi}(\mathcal{D}_{ij}\mid\bm{X}^{0})$. Then $\Delta$
satisfies the following with $\mathbb{P}^{0}$-probability equal to $1$ for all
$i,j,l=1,\dots,n$: (a) $0\leq\Delta_{ij}\leq 1$; (b) $\Delta_{ii}=0$,
$\Delta_{ji}=\Delta_{ij}$, and $\Delta_{il}\leq\Delta_{ij}+\Delta_{jl}$; and
(c) $\Delta_{ij}\leq\Pi(s_{i}\neq s_{j}\mid\bm{X}^{0})$.

Properties (a) and (b) follow naturally from the assumptions we have made
regarding the statistical distance $d$. To show (c), recall that, conditional
on $\lambda=\sum_{k=1}^{K}a_{k}\delta_{\tilde{\theta}_{k}}$, variation in
$\theta_{i}$ and $\theta_{j}$ is explained by sampling from $\lambda$ in (1).
This implies that
$\mathbb{E}_{\Pi}(\mathcal{D}_{ij}\mid\bm{X}^{0},\lambda)=\sum_{k\neq
k^{\prime}}d\left\\{g(\cdot;\tilde{\theta}_{k}),g(\cdot;\tilde{\theta}_{k^{\prime}})\right\\}\Pi(s_{i}=k,s_{j}=k^{\prime}\mid\bm{X}^{0},\lambda)$.
Since $0\leq d(\cdot,\cdot)\leq 1$, we have that
$\mathbb{E}_{\Pi}(\mathcal{D}_{ij}\mid\bm{X}^{0},\lambda)\leq\sum_{k\neq
k^{\prime}}\Pi(s_{i}=k,s_{j}=k^{\prime}\mid\bm{X}^{0},\lambda)=\Pi(s_{i}\neq
s_{j}\mid\bm{X}^{0},\lambda)$. Taking the expectation with respect to the
posterior of $\lambda$ on both sides of the inequality results in (c).
Properties (a) and (b) imply that $\mathcal{R}(\widehat{\bm{c}})$ is non-
negative for any clustering and that the sum in (5) may be taken over the
indices $i<j$. Property (c) states that FOLD will generally favor a fewer
number of clusters than Binder’s loss. One can show that Binder’s loss prefers
assigning $i$ and $j$ to the same cluster when $\Pi(s_{i}\neq
s_{j}\mid\bm{X}^{0})<\gamma$. If $\Delta_{ij}<\gamma<\Pi(s_{i}\neq
s_{j}\mid\bm{X}^{0})$, FOLD will disagree with Binder’s loss, with FOLD
preferring to merge $i$ and $j$ into the same cluster, implying that
$\bm{c}_{\textrm{FOLD}}$ is often more coarse then
$\bm{c}_{\textrm{B}}=\underset{\widehat{\bm{c}}}{\textrm{argmin
}}\mathbb{E}_{\Pi}\left\\{\mathcal{L}_{\textrm{B}}(\widehat{\bm{c}},\bm{s})\mid\bm{X}^{0}\right\\}$.

Recently, the Variation of Information (VI) (Meilă,, 2007) has attracted
attention in the literature for use as a clustering loss function (Wade and
Ghahramani,, 2018; De Iorio et al.,, 2023; Denti et al.,, 2023), and we
compare FOLD to the VI from a methodological point of view in the Supplement.
The VI is primarily motivated from an information theoretic perspective, but
it, like Binder’s loss, is a metric on the partition space (Meilă,, 2007; Dahl
et al.,, 2022) and encourages similarity to the component labels when used in
Bayesian clustering. While the VI has been observed to yield fewer clusters
than Binder’s loss empirically on synthetic and real data (Wade and
Ghahramani,, 2018), we show in a simulation study in the Supplementary and an
application in Section 4.2 that $\bm{c}_{\textrm{VI}}$ is still vulnerable to
the over-clustering problem, albeit usually to a lesser degree than
$\bm{c}_{\textrm{B}}$, whereas $\bm{c}_{\textrm{FOLD}}$ is more robust.

###  2.3 Uncertainty Quantification

In applications of Bayesian clustering to real data, we often find that there
is substantial uncertainty in the component labels a posteriori. To avoid
overstating the significance of $\bm{c}_{\textrm{FOLD}}$, it is thus important
to express this uncertainty in a clear and interpretable manner. We focus on
the clustering that minimizes (2), i.e.
$\bm{c}_{\bm{\theta}}=\underset{\widehat{\bm{c}}}{\textrm{argmin
}}\mathcal{L}(\widehat{\bm{c}},\bm{\theta}).$ Observe that
$\bm{c}_{\bm{\theta}}$ depends on $\omega$ and may not be equal to $\bm{s}$.
To see this, we recall (4), which shows that we can acquire a clustering with
loss smaller than $\bm{s}$ by merging any two components with statistical
distance below $\gamma$. Thus, we formulate our measures of uncertainty with
respect to the FOLD posterior $\Pi(\bm{c}_{\bm{\theta}}\mid\bm{X}^{0})$,
instead of $\Pi(\bm{s}\mid\bm{X}^{0})$.

Accordingly, we adapt the notion of a credible ball for Bayesian clustering
estimators (Wade and Ghahramani,, 2018) to our procedure. The 95% credible
ball around $\bm{c}_{\textrm{FOLD}}$ is defined as
$B(\bm{c}_{\textrm{FOLD}})=\left\\{\bm{c}:\mathbb{D}(\bm{c}_{\textrm{FOLD}},\bm{c})\leq\epsilon_{\textrm{FOLD}}\right\\}$,
where $\epsilon_{\textrm{FOLD}}\geq 0$ is the smallest radius such that
$\Pi\left\\{\bm{c}_{\bm{\theta}}\in
B(\bm{c}_{\textrm{FOLD}})\mid\bm{X}^{0}\right\\}\geq 0.95$. We interpret
$B(\bm{c}_{\textrm{FOLD}})$ as a neighborhood of clusterings centered at
$\bm{c}_{\textrm{FOLD}}$ with posterior probability mass of at least $0.95$.
Larger $\epsilon_{\textrm{FOLD}}$ means that the FOLD posterior distributes
its mass across a larger set around $\bm{c}_{\textrm{FOLD}}$, implying that we
are more uncertain about the cluster allocations of the point estimator. In
practice, we characterize the credible ball with bounds that give a sense of
the clusterings contained within it (Wade and Ghahramani,, 2018), much like we
would for an interval on the real line. The horizontal bounds of the credible
ball consist of the clusterings $\bm{c}\in B(\bm{c}_{\textrm{FOLD}})$ for
which $\mathbb{D}(\cdot,\bm{c}_{\textrm{FOLD}})$ attains its maximum value. We
can also impose further restrictions on our bounds, such as requiring that
they contain the minimum or maximum number of clusters over
$B(\bm{c}_{\textrm{FOLD}})$ while also maximizing
$\mathbb{D}(\cdot,\bm{c}_{\textrm{FOLD}})$ amongst groupings with the same
number of clusters, giving rise to the notions of vertical upper and vertical
lower bounds, respectively. Alternatively, one can display the posterior
similarity matrix (PSM) for $\bm{c}_{\bm{\theta}}$,
$\mathcal{P}_{\textrm{FOLD}}=(\Pi\left\\{c_{i\bm{\theta}}=c_{j\bm{\theta}}\mid\bm{X}^{0}\right\\})_{i,j}$,
as a heatmap, with the entries ordered by $\bm{c}_{\textrm{FOLD}}$, as is
common practice for both Binder’s loss and the VI.

###  2.4 Implementation with MCMC Output

Our methodology relies on Markov chain Monte Carlo (MCMC) samples from
$\Pi(\bm{\theta}\mid\bm{X}^{0})$, which are used to estimate pairwise
distances between the localized densities,
$\Delta_{ij}\approx(1/T)\sum_{t=1}^{T}\mathcal{D}_{ij}^{(t)},$ where
$\mathcal{D}_{ij}^{(t)}=d\left\\{g(\cdot;\theta_{i}^{(t)}),g(\cdot;\theta_{j}^{(t)})\right\\}$,
and $\theta_{i}^{(t)},\theta_{j}^{(t)}$ can be computed directly from samples
of the components $\bm{s}^{(t)}$ and unique values $\tilde{\bm{\theta}}^{(t)}$
via $\theta_{i}^{(t)}=\tilde{\theta}^{(t)}_{k}\iff s_{i}^{(t)}=k$. When $d$ is
chosen to be the Hellinger distance, many families of kernels, such as the
Gaussian family, admit a closed form of $\mathcal{D}_{ij}$. The label-
switching problem, which results from the inherent non-identifiability of
mixture models with any permutation of the labels yielding the same likelihood
(Redner and Walker,, 1984), has no impact on FOLD, as $\Delta_{ij}$ is
invariant to labeling. FOLD can be easily applied to virtually any mixture
model for which samples from $\Pi(\bm{\theta}\mid\bm{X}^{0})$ can be obtained
and pairwise statistical distances can be estimated.

Figure 3: Candidate clusterings of the moons data under a location Gaussian
mixture model, generated from average linkage hierarchical clustering on
$\Delta$.

In practice, the minimization of (5) must be approximated due to the enormous
size of the partition space for even modest $n$. We discuss two methods for
obtaining an approximate minimizer. The first is to simply reduce the set of
partitions over which we minimize by only considering a tree of clusterings
produced by hierarchical clustering on $\Delta$. We suggest generating
candidates using average linkage due to the appearance of the average linkage
dissimilarity metric as a criterion for merging clusters in (4). Figure 3
displays candidate clusterings of the moons data resulting from average
linkage clustering on $\Delta$. The true grouping of the data is present
amongst the candidates. Since these groupings are hierarchical, Figure 3 can
be interpreted as an enumeration of clusterings favored by
$\mathcal{R}(\widehat{\bm{c}})$, arranged by increasing values of $\omega$. A
similar method with single-linkage is defined in Aragam et al., (2020), and
Medvedovic and Sivaganesan, (2002) and Fritsch and Ickstadt, (2009) have
previously used hierarchical clustering on
$\textbf{1}_{n}\textbf{1}_{n}^{T}-\mathcal{P}$ to compute clustering point
estimators, where
$\mathcal{P}=(\Pi\left\\{s_{i}=s_{j}\mid\bm{X}^{0}\right\\})_{i,j}$.

Recent developments for minimizing risk functions have focused on greedy
algorithms that take an existing candidate partition and make minor, but
locally-optimal updates, typically by reallocation of objects to clusters.
Examples include a method motivated by the Hasse diagram in Wade and
Ghahramani, (2018), the greedy algorithm in Rastelli and Friel, (2018), and
the SALSO algorithm of Dahl et al., (2022). SALSO provides a point estimate by
first initializing a clustering, then reallocating observations to clusters
and, finally, breaking existing clusters and then reallocating observations
again over several parallel runs. The final point estimate is the partition
that exhibits the smallest risk across all the runs. In our simulations and
application, we implement FOLD with both the hierarchical clustering heuristic
and SALSO, and find that they generally agree on the optimal clustering.

To express uncertainty in $\bm{c}_{\textrm{FOLD}}$, we rely on samples from
$\Pi(\bm{c}_{\bm{\theta}}\mid\bm{X}^{0})$, defined by

$\bm{c}_{\bm{\theta}}^{(t)}=\underset{\widehat{\bm{c}}}{\textrm{argmin
}}\sum_{i<j}\left\\{\textbf{1}_{\hat{c}_{i}=\hat{c}_{j}}\mathcal{D}^{(t)}_{ij}+\omega\textbf{1}_{\hat{c}_{i}\neq\hat{c}_{j}}\left(1-\mathcal{D}_{ij}^{(t)}\right)\right\\}.$
(6)

As in the point estimation case, $\bm{c}_{\bm{\theta}}^{(t)}$ cannot be
computed exactly and so we approximate the minimization in (6) either using a
set of candidates generated from hierarchical clusterings of
$\mathcal{D}^{(t)}=(\mathcal{D}^{(t)}_{ij})_{i,j}$ or via a greedy algorithm
such as SALSO. The posterior probability of a credible ball is approximated by
$\Pi\left\\{\mathbb{D}(\bm{c}_{\bm{\theta}},\bm{c}_{\textrm{FOLD}})\leq\epsilon\mid\bm{X}^{0}\right\\}\approx(1/T)\sum_{t=1}^{T}\textbf{1}_{\mathbb{D}(\bm{c}_{\bm{\theta}}^{(t)},\bm{c}_{\textrm{FOLD}})\leq\epsilon}$
and FOLD computes $\epsilon_{\textrm{FOLD}}$ by incrementally increasing the
value of $\epsilon$ over a grid. The entries of $\mathcal{P}_{\textrm{FOLD}}$
are estimated via
$\Pi(c_{i\bm{\theta}}=c_{j\bm{\theta}}\mid\bm{X}^{0})\approx(1/T)\sum_{t=1}^{T}\textbf{1}_{c^{(t)}_{i\bm{\theta}}=c^{(t)}_{j\bm{\theta}}}$.

###  2.5 Selection of $\omega$

The loss parameter $\omega$ controls the separation among the inferred
clusters in $\bm{c}_{\textrm{FOLD}}$. This parameter is also influential in
$\bm{c}_{\textrm{VI}}$ and $\bm{c}_{\textrm{B}}$ (Dahl et al.,, 2022), though
there has been little discussion in the literature on how to select an
appropriate $\omega$ for these loss functions. The formulation of (2) in terms
of statistical distances leads to two sensible methods for determining the
loss parameter. First, similarly to common practice in choosing key tuning
parameters in algorithmic clustering methods such as k-means and DBSCAN, we
propose an elbow plot diagnostic to select $\omega$. At a grid of possible
$\omega$ values, we calculate
$r_{\omega}=(\sum_{h=1}^{K_{\omega}^{*}}r_{\omega
h})/{\sum_{i<j}\Delta_{ij}}$, where $r_{\omega h}=\sum_{i,j\in C_{\omega
h}^{*}}\Delta_{ij}$, $\bm{c}_{w}^{*}$ is the FOLD point estimator for that
specific $\omega$, and $C^{*}_{\omega}=\left\\{C_{\omega
1}^{*},\dots,C_{\omega K_{w}^{*}}^{*}\right\\}$ is the partition associated
with $\bm{c}_{\omega}^{*}$. The numerator in $r_{w}$ sums over $r_{\omega
h}=\mathbb{E}_{\Pi}\left\\{\sum_{i,j\in C_{\omega
h}^{*}}\mathcal{D}_{ij}\mid\bm{X}^{0}\right\\}$, or the expected statistical
distance between all localized densities associated with the $h$th cluster in
$\bm{c}_{\omega}^{*}$, and the denominator normalizes so that
$r_{\omega}\in[0,1]$ for all $\omega$.

To build some intuition about the utility of the elbow plot, recall that as
$\omega\to\infty$, (5) is minimized by placing all observations in a single
cluster, i.e. $C_{\infty}^{*}=\left\\{[n]\right\\}$ and $r_{\infty}=1$. In
many cases, such as the data in Figure 3, decreasing $\omega$ (and
subsequently adding more clusters in $C_{\omega}^{*}$) will decrease
$r_{\omega}$, e.g. because $\Delta_{ij}\approx 1$ for any $i,j$ pairs in
different crescents. Of course, at the other extreme, if we set $\omega$ to be
very small, we place all objects in their own cluster, making
$C_{0}^{*}=\left\\{\left\\{1\right\\},\dots,\left\\{n\right\\}\right\\}$ and
$r_{0}=0$. We construct the elbow plot to quantify the amount of improvement
each time we increment $\omega$ across this grid. At some point, the
improvement in $r_{\omega}$ becomes negligible as we decrease $\omega$ because
we begin to cut dense clusters into smaller and smaller sub-clusters. On an
elbow plot, this causes the curve of $r_{\omega}$ to bend, creating the so-
called “elbow". It is at this threshold that we fix $\omega$. The threshold is
easier to identify when the elbow plot is monotonic, and we verify this is the
case if one uses hierarchical clustering to generate candidates.

######  Proposition 3.

If candidate clusterings are selected using average linkage hierarchical
clustering on $\Delta$, then there is a positive relationship between $\omega$
and $r_{\omega}$.

A corollary of Proposition 3 is that any jump in the elbow plot corresponds to
decreasing the number of clusters by one, meaning that there is a negative
relationship between $K_{\omega}^{*}$ and $r_{\omega}$. Therefore, one can
relabel the x-axis of the elbow plot in terms of the number of clusters rather
than $\omega$, leading to a more interpretable figure. An example of using an
elbow plot to select $\omega$ is given in Section 4.2.

Alternatively, one can use the default value
$\omega^{\textrm{AVG}}=\gamma^{\textrm{AVG}}/(1-\gamma^{\textrm{AVG}})$, with
$\gamma^{\textrm{AVG}}={n\choose 2}^{-1}\sum_{i<j}\Delta_{ij}$. Under this
choice of $\omega$, candidate clusters are combined when the average of
$\Delta_{ij}$ between the clusters is smaller than the average of
$\Delta_{ij}$ across the entire sample. $\gamma^{\textrm{AVG}}$ estimates
$\bar{\mathcal{D}}={n\choose
2}^{-1}\sum_{i<j}\mathcal{D}_{ij}=\sum_{k<k^{\prime}}(|S_{k}||S_{k^{\prime}}|/{n\choose
2})d\left\\{g(\cdot;\tilde{\theta}_{k}),g(\cdot,\tilde{\theta}_{k^{\prime}})\right\\}$,
a weighted pairwise sum across the components. If we use
$\bar{\mathcal{D}}/(1-\bar{\mathcal{D}})$ in place of $\omega$ in (2), (4)
implies that FOLD will favor merging mixture components when
$d\left\\{g(\cdot;\tilde{\theta}_{k}),g(\cdot;\tilde{\theta}_{k^{\prime}})\right\\}<\bar{\mathcal{D}}$.
Importantly, the decision to merge components will depend on how separated
they are from the others and their sizes. To see this, consider the following
example, in which we fit a mixture with $K=3$ components, where
$d\left\\{g(\cdot;\tilde{\theta}_{1}),g(\cdot;\tilde{\theta}_{2})\right\\}=\epsilon>0$,
$d\left\\{g(\cdot;\tilde{\theta}_{1}),g(\cdot;\tilde{\theta}_{3})\right\\}=d\left\\{g(\cdot;\tilde{\theta}_{2}),g(\cdot;\tilde{\theta}_{3})\right\\}=\delta>0$,
and $|S_{1}|=|S_{2}|=|S_{3}|=n/3$. Then, as $n$ grows large,
$\bar{\mathcal{D}}\to(2/9)\epsilon+(4/9)\delta$, and so FOLD will favor
merging $S_{1}$ and $S_{2}$ into one cluster when $\epsilon<(4/7)\delta$. The
smaller $\delta$ is, the smaller $\epsilon$ must be in order to merge $S_{1}$
and $S_{2}$. Hence, $\omega^{\textrm{AVG}}$ excels at problems in which
$f^{0}$ is composed of well-separated kernels that are approximated by
multiple components in $f$. We provide simulation studies in Section 4.1 and
the Supplement that show that $\omega^{\textrm{AVG}}$ performs very well in
this setting.

##  3 Asymptotic Analysis

Though there is a rich literature on asymptotic properties of Bayesian mixture
models for estimating the density (Ghosal et al.,, 1999; Ghosal and van der
Vaart,, 2007), number of components (Rousseau and Mengersen,, 2011; Miller and
Harrison,, 2014; Cai et al.,, 2021; Ascolani et al.,, 2023), and the mixing
measure (Nguyen,, 2013; Ho and Nguyen, 2016a, ; Ho and Nguyen, 2016b, ; Guha
et al.,, 2021), little attention has been given to the large sample behavior
of clustering estimators. Rajkowski, (2019) focused on the maximum a
posteriori (MAP) estimator for $\bm{s}$ when data are generated from a DPM
with Gaussian components and the kernels are correctly specified. They show
multiple properties of the MAP estimator, including the key result that the
intersection between the convex hulls of two clusters is at most a single
observation. In this section, we show convergence of FOLD towards the oracle
clustering procedure in which the objects are partitioned into groups using
knowledge of the true dating generating process. We take a different approach
than Rajkowski, (2019) by focusing on posterior contraction of the mixing
measure $\lambda$, which allows us to consider misspecifed models, either in
the kernels or number of components.

###  3.1 Assumptions

We first require that the parameter space $\Theta\subset\mathbb{R}^{p}$ is
compact and $f$, defined formally as $f(x)=\int g(x;\theta)\lambda(d\theta)$,
is identifiable with respect to the mixing measure, that is, $\int
g(x;\theta)\lambda_{1}(d\theta)=\int g(x;\theta)\lambda_{2}(d\theta)$ for all
$x\in\chi$ if and only if $\lambda_{1}=\lambda_{2}$ (Teicher,, 1961).
Identifiability of the mixing measure is satisfied by location and location-
scale GMMs, as well as various exponential family mixtures (Barndorff-
Nielsen,, 1965) and location family mixtures (Teicher,, 1961). More
specifically, our results rely on contraction of the mixing measure $\lambda$
to some oracle measure $\lambda^{*}$, which implies convergence of $f$ to an
oracle model $f^{*}=\int g(\cdot;\theta)\lambda^{*}(d\theta)$. An oracle
measure is defined as any Kullback-Liebler (KL) divergence minimizer between
our class of models and the true data generating process, i.e.
$\lambda^{*}\in\underset{\lambda\in\Lambda}{\textrm{argmin
}}\textrm{KL}(f^{0},f)$, where recall $\Lambda=\textrm{supp}(\pi_{\Lambda})$.
We now present the general assumptions required for our results.

######  Assumption 1.

Suppose that the following conditions hold.

  * (A1) There exists an $L>0$ such that $|\textrm{supp}(\lambda)|\leq L$ for all $\lambda\in\Lambda$.

  * (A2) The KL minimizer $\lambda^{*}=\sum_{m=1}^{M^{*}}a_{m}^{*}\delta_{\theta_{m}^{*}}$ exists and is unique.

  * (A3) $f^{*}$ and $f^{0}$ are such that $\textrm{supp}(f^{0})\subseteq\textrm{supp}(f^{*})$.

  * (A4): $\mathcal{G}=\left\\{\tilde{g}(x-\theta):\theta\in\Theta\right\\}$ for some bounded probability density function $\tilde{g}(\cdot)$, $\tilde{g}(\cdot)$ is $\zeta$-Hölder continuous for some $\zeta>0$, and for any fixed $\theta^{\prime}$, $D_{g}(\theta,\theta^{\prime})=d\left\\{g(\cdot;\theta),g(\cdot;\theta^{\prime})\right\\}$ is continuous in $\theta$.

  * (A5) There exists a non-negative sequence $\epsilon_{n}$ so that $\epsilon_{n}\to 0$ and

$\rho_{n}(\bm{X}^{0})=\Pi\left\\{\lambda\in\Lambda:W_{2}(\lambda,\lambda^{*})\underset{\sim}{>}\epsilon_{n}\mid\bm{X}^{0}\right\\}\overset{\mathbb{P}^{0}}{\longrightarrow}0.$
(7)

(A1) is satisfied by any mixture model in which the number of components is
bounded, and examples include a finite mixture model, a truncated MFM, or
truncated DPM. Assumption (A2) ensures that the limiting values of the
clustering procedure are well-defined and exist. A simple example of when (A2)
holds is when the kernels are correctly specified and $\Pi(K=M^{0})>0$, in
which case $\lambda^{*}=\lambda^{0}$ and $f^{*}=f^{0}$. Under any
misspecification, $f^{*}$ is the best possible approximation to $f^{0}$ by the
fitted statistical model. (A3) is a technical assumption to ensure there are
no regions in the sample space under which $f^{0}$ gives positive probability
but $f^{*}$ does not. Assumption (A4) restricts our results to location
families, and the continuity conditions are satisfied by a variety of models
including the location GMM with fixed covariance. Assumption (A5) is a crucial
statement on the estimation of parameters in a mixture model. (A5) will hold
under various Lipschitz and strong identifiability conditions on $\mathcal{G}$
(Nguyen,, 2013; Ho and Nguyen, 2016a, ; Guha et al.,, 2021) that are satisfied
by the location GMM (Manole and Ho,, 2022), and we discuss these conditions in
detail in the Supplementary Material. (A5) will hold with $\epsilon_{n}=(\log
n/n)^{1/4}$ for second-order identifiable families, and first-order
identifiable families with $\epsilon_{n}=(\log n/n)^{1/2}$ if the kernels are
correctly specified and $\Pi(K=M^{0})=1$ (Guha et al.,, 2021); see Ho and
Nguyen, 2016a ; Ho and Nguyen, 2016b  for a characterization of first and
second-order identifiable kernels.

###  3.2 Convergence to The Oracle Rule

Synonymous with the oracle mixing measure is the oracle clustering procedure
or rule. We define the oracle clustering procedure to be the minimizer of the
oracle risk function
$\mathcal{R}^{*}(\widehat{\bm{c}})=\mathbb{E}_{\Pi}[\mathcal{L}(\widehat{\bm{c}},\bm{\theta})\mid\lambda^{*},\bm{X}^{0}]$,

$\bm{c}_{\textrm{FOLD}}^{*}=\underset{\widehat{\bm{c}}}{\textrm{argmin
}}\mathcal{R}^{*}(\widehat{\bm{c}})=\underset{\widehat{\bm{c}}}{\textrm{argmin
}}\sum_{i<j}\left\\{\textbf{1}_{\hat{c}_{i}=\hat{c}_{j}}\Delta_{ij}^{*}+\omega\textbf{1}_{\hat{c}_{i}\neq\hat{c}_{j}}(1-\Delta_{ij}^{*})\right\\},$
(8)

where
$\Delta^{*}_{ij}=\mathbb{E}_{\Pi}(\mathcal{D}_{ij}\mid\lambda^{*},\bm{X}^{0})=\sum_{m<m^{\prime}}d\left\\{g(\cdot;\theta_{m}^{*}),g(\cdot;\theta_{m^{\prime}}^{*})\right\\}q_{ij}^{mm^{\prime}*}$,
with
$q_{ij}^{mm^{\prime}*}=\Pi(s_{i}=m,s_{j}=m^{\prime}\mid\bm{X}^{0},\lambda^{*})+\Pi(s_{i}=m^{\prime},s_{j}=m\mid\bm{X}^{0},\lambda^{*})$,
and

$\Pi(s_{i}=m,s_{j}=m^{\prime}\mid\bm{X}^{0},\lambda^{*})=\frac{a_{m}^{*}a_{m^{\prime}}^{*}g(X_{i}^{0};\theta_{m}^{*})g(X_{j}^{0};\theta_{m^{\prime}}^{*})}{f^{*}(X_{i}^{0})f^{*}(X_{j}^{0})}.$

Observe that $\Delta_{ij}^{*}$ is a weighted sum of the total statistical
distance between the oracle components, where the weight of each component
pair is given by the conditional probability that $s_{i}=m$ and
$s_{j}=m^{\prime}$ (or vice versa) given $\bm{X}=\bm{X}^{0}$ and
$\lambda=\lambda^{*}$. We interpret the oracle clustering procedure as a rule
for grouping the observations in $\bm{X}^{0}$ based on the knowledge of the
optimal parameter values in the mixture model. We can construct two simple
examples of oracle rules which follow from the definition in (8).

######  Proposition 4.

Instances of the oracle clustering procedure include:

  1. (a)

If
$d\left\\{g(\cdot;\theta_{m}^{*}),g(\cdot;\theta_{m^{\prime}}^{*})\right\\}=1$
for all $m,m^{\prime}$, then $\bm{c}_{\textrm{FOLD}}^{*}=\bm{s}^{*}$, where
$s^{*}_{i}=s^{*}_{j}$ if and only if there exists a unique $m\in[M^{*}]$ such
that $g(X_{i}^{0};\theta_{m}^{*})g(X_{j}^{0};\theta_{m}^{*})>0$.

  2. (b)

If $M^{*}=2$ and
$d\left\\{g(\cdot;\theta_{1}^{*}),g(\cdot;\theta_{2}^{*})\right\\}<\gamma$,
then $\bm{c}_{\textrm{FOLD}}^{*}=\bm{c}_{0}$, where $c_{0i}=c_{0j}$ for all
$i,j$.

We refer to the procedure on Proposition 4(a) as the “match rule" in that it
allocates each observation to a mixture component with compatible support, and
call the procedure in Proposition 4(b) the “merge rule", i.e., as the
statistical distance between oracle components becomes smaller, we are more
likely to combine them and subsequently place all observations in one cluster.
A similar quantity exists for Binder’s loss function, denoted
$\bm{c}_{\textrm{B}}^{*}=\underset{\widehat{\bm{c}}}{\textrm{argmin
}}\mathbb{E}_{\Pi}[\mathcal{L}_{\textrm{B}}(\widehat{\bm{c}},\bm{s})\mid\lambda^{*},\bm{X}^{0}]$.
Under the settings of Proposition 4(a),
$\bm{c}_{\textrm{B}}^{*}=\bm{c}_{\textrm{FOLD}}^{*}=\bm{s}^{*}$, meaning that
the match rule is the optimal clustering procedure under both FOLD and
Binder’s loss when the oracle components are perfectly separated. Intuitively,
this is because these settings reflect an ideal scenario in which all mixture
components are perfectly separated and no merging is required. Otherwise, say
in the setting of Proposition 4(b), the two methods can differ. A sufficient
condition for $\bm{c}_{\textrm{B}}^{*}=\bm{c}_{0}$ is
$q_{ij}^{12}<\omega(1-q_{ij}^{12})$ for all pairs $i,j$. Any observations in
the tails of the two oracle components will blow up the value of
$q_{ij}^{12}$, hence, it is possible that
$\bm{c}_{\textrm{B}}^{*}\neq\bm{c}_{0}$ even for relatively large $\omega$.
This implies that a concrete advantage of our approach is robustness to
outliers when the oracle components are not well separated. We now state our
primary theoretical result on the relationship between
$\bm{c}_{\textrm{FOLD}}$ and $\bm{c}_{\textrm{FOLD}}^{*}$.

######  Theorem 1.

Fix $0<\delta<1$. Then under assumptions (A1)-(A5), with
$\mathbb{P}^{0}$-probability tending to $1$,

${n\choose
2}^{-1}|\mathcal{R}(\bm{c}_{\textrm{FOLD}})-\mathcal{R}^{*}(\bm{c}_{\textrm{FOLD}}^{*})|\underset{\sim}{<}\max(\omega,1)\left\\{\rho_{n}(\bm{X}^{0})\bar{\Delta}^{*}+\mathcal{Q}(\bm{X}^{0})\right\\},$
(9)

where $\bar{\Delta}^{*}={n\choose 2}^{-1}\sum_{i<j}\Delta^{*}_{ij}$,
$\mathcal{Q}(\bm{X}^{0})={n\choose
2}^{-1}\sum_{i<j}\tau^{*}_{n}(X_{i}^{0},X_{j}^{0})$, and
$\tau_{n}^{*}(X_{i}^{0},X_{j}^{0})$ is a function of
$(X_{i}^{0},X_{j}^{0},\epsilon_{n},\rho_{n}(\bm{X}^{0}),\lambda^{*},\zeta,\delta)$
so that $\tau_{n}^{*}(X_{i},X_{j})\overset{\mathbb{P}^{0}}{\longrightarrow}0$.

To prove Theorem 1, we first note that we can write $\Delta_{ij}$ as the
posterior expectation of a quantity that only depends on
$\lambda=\sum_{k=1}^{K}a_{k}\delta_{\tilde{\theta}_{k}}$,

$\Delta_{ij}=\int_{\lambda\in\Lambda}\Pi(d\lambda\mid\bm{X}^{0})\sum_{k<k^{\prime}}D_{g}(\tilde{\theta}_{k},\tilde{\theta}_{k^{\prime}})a_{k}a_{k^{\prime}}\frac{g(X_{i}^{0};\tilde{\theta}_{k})g(X_{j}^{0};\tilde{\theta}_{k^{\prime}})+g(X_{j}^{0};\tilde{\theta}_{k})g(X_{i}^{0};\tilde{\theta}_{k^{\prime}})}{f(X_{i}^{0})f(X_{j}^{0})}.$
(10)

Next, we decompose $\Lambda=B_{\epsilon_{n}}(\lambda^{*})\cup(\Lambda\setminus
B_{\epsilon_{n}}(\lambda^{*}))$, where
$B_{\epsilon_{n}}(\lambda^{*})=\left\\{\lambda\in\Lambda:W_{2}(\lambda,\lambda^{*})\underset{\sim}{<}\epsilon_{n}\right\\}$,
and split (10) into a sum of integrals over these domains. Since $0\leq
D_{g}(\theta,\theta^{\prime})\leq 1$, the integral over $\Lambda\setminus
B_{\epsilon_{n}}(\lambda)^{*}$ is bounded above by $\rho_{n}(\bm{X}^{0})$,
which goes to $0$ in probability by (A5). Then, using a similar proof
technique given in Nguyen, (2013), we show that there exists a finite $N$ so
that for all $n\geq N$ and $\lambda\in B_{\epsilon_{n}}(\lambda^{*})$,
$\lambda$ has at least $M^{*}$ support points and for any $m\in[M^{*}]$ there
exists a set of indices $\mathcal{I}_{m}^{(\lambda)}\subset[K]$ so that
$\max_{k\in\mathcal{I}_{m}^{(\lambda)}}\left\lVert\tilde{\theta}_{k}-\theta_{m}^{*}\right\rVert\underset{\sim}{<}\epsilon_{n}$
and
$\bigg{|}\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}-a_{m}^{*}\bigg{|}\underset{\sim}{<}\max(\epsilon_{n},\epsilon_{n}^{2\delta})$.
Furthermore, if we denote
$\mathcal{I}^{(\lambda)}=\bigcup_{m=1}^{M^{*}}\mathcal{I}^{(\lambda)}_{m}$,
(A5) also implies that
$\sum_{k\not\in\mathcal{I}^{(\lambda)}}a_{k}\underset{\sim}{<}\epsilon_{n}^{2\delta}$.
By (A4), these bounds translate into similar upper and lower bounds for
$\Delta_{ij}$. We then sum over all $i,j$ pairs and show that ${n\choose
2}^{-1}|\mathcal{R}(\widehat{\bm{c}})-\mathcal{R}^{*}(\widehat{\bm{c}})|$ is
uniformly bounded by the right hand side of (9) for all $\widehat{\bm{c}}$.
That is, the rate of convergence of $\mathcal{R}(\widehat{\bm{c}})$ to
$\mathcal{R}^{*}(\widehat{\bm{c}})$ does not depend on the clustering
estimator itself. Finally, we show that this implies convergence of the
minimizers at the same rate.

The bound in (9) is split into two remainder terms. It follows from (A5) that
$\rho_{n}(\bm{X}^{0})\bar{\Delta}^{*}\overset{\mathbb{P}^{0}}{\longrightarrow}0$
since $0\leq\bar{\Delta}^{*}\leq 1$. Therefore, smaller values of
$\bar{\Delta}^{*}$ will cause this term to diminish faster. This can arise
when $f^{*}$ is comprised of many closely overlapping components. Hence,
convergence to the oracle rule will actually benefit from overfitting.
$\mathcal{Q}(\bm{X}^{0})$ is a more general remainder term that results from
translating the posterior contraction of $\tilde{\theta}_{k}$ and $a_{k}$ to
the sum in (10). In particular, this remainder term will depend on $\delta$,
which controls the rate at which the irrelevant mixture components are
emptied, the Hölder constant $\zeta$, which controls convergence of
$g(X_{i}^{0};\tilde{\theta}_{k})$ to $g(X_{i}^{0};\theta_{m}^{*})$, and
$L-M^{*}$, or the degree to which the number of components is misspecified,
with larger $L$ implying a slower rate of convergence.

In the Supplementary Material we empirically validate Theorem 1 on a simple
simulated example. We simulate data from $M^{0}=4$ Gaussian components with
covariance equal to $0.25I$, then fit a $4$-component Bayesian mixture model,
meaning that
$f^{0}(x)=f^{*}(x)=\sum_{m=1}^{4}a_{m}^{0}\mathcal{N}(x;\mu_{m}^{0},0.25I)$.
We compute $\bm{c}_{\textrm{FOLD}}$, $\bm{c}_{\textrm{FOLD}}^{*}$,
$\bm{c}_{\textrm{B}}$, and $\bm{c}_{\textrm{B}}^{*}$, the latter of which
denotes the oracle clusters corresponding to Binder’s loss, for increasing
sample size. We show that both
$\mathcal{R}(\widehat{\bm{c}})\to\mathcal{R}^{*}(\widehat{\bm{c}})$ and
$\bm{c}_{\textrm{FOLD}}\to\bm{c}_{\textrm{FOLD}}^{*}$ as $n\to\infty$.
Additionally, there are substantial differences between
$\bm{c}_{\textrm{FOLD}}^{*}$ and $\bm{c}_{\textrm{B}}^{*}$ because oracle FOLD
corresponds to the merge rule for some of the Gaussian components, while
oracle Binder prefers estimation of each individual component.

##  4 Empirical Performance

###  4.1 Synthetic Data

| $n$ | FOLD | VI | Binder | Mclust
---|---|---|---|---|---
ARI | 100 | 0.982 (0.029) | 0.974 (0.026) | 0.800 (0.094) | 0.884 (0.177)
| 500 | 0.995 (0.032) | 0.921 (0.136) | 0.799 (0.147) | 0.715 (0.159)
| 1000 | 0.975 (0.085) | 0.691 (0.057) | 0.665 (0.037) | 0.574 (0.054)
| 2500 | 0.967 (0.098) | 0.679 (0.012) | 0.639 (0.027) | 0.512 (0.056)
$\hat{K}_{n}$ | 100 | 3.25 (0.500) | 3.89 (1.014) | 12.62 (4.292) | 3.56 (0.925)
| 500 | 3.23 (0.601) | 3.64 (0.894) | 13.94 (4.552) | 4.17 (0.779)
| 1000 | 3.10 (0.302) | 4.83 (1.025) | 15.73 (3.795) | 5.08 (0.442)
| 2500 | 3.10 (0.302) | 4.58 (0.843) | 18.99 (2.634) | 6.20 (0.876)
Table 1: Averages (and standard deviations) for the adjusted Rand index (ARI)
and number of clusters ($\hat{K}_{n}$) with $\bm{s}^{0}$ on 100 replications
from the skew-symmetric mixture.

To gauge overall performance, we evaluate FOLD in comparison to clusterings
produced by other model-based methods in the presence of kernel
misspecification. We focus on a particular choice for $f^{0}$ we call the
“skew-symmetric mixture". The skew-symmetric mixture consists of $M^{0}=3$
well separated clusters, where one cluster is a Gaussian kernel, another is a
skew Gaussian kernel, and the third cluster is itself a mixture over two
Gaussian components and one skew-Gaussian component. We vary the sample size
$n\in\left\\{100,500,1000,2500\right\\}$ and, for each sample size, generate
$100$ independent replications. In each replication, we fit a Bayesian finite
Gaussian mixture model with $K=30$ components, then use the MCMC samples to
estimate clusterings obtained from FOLD with $\omega=\omega^{\textrm{AVG}}$,
Binder’s loss, and the VI (where both have $\omega=1$, the default value in
the mcclust package (Fritsch,, 2022)). In addition, we use mclust (Scrucca et
al.,, 2016) to implement the EM algorithm as another point of comparison, with
the number of clusters chosen by minimizing the Bayesian information criterion
(BIC). Extensive details on the simulations are given in the Supplementary
Material, including two additional scenarios where we sample from a Gaussian
mixture and a skew Gaussian mixture.

The adjusted Rand index (ARI) (Hubert and Arabie,, 1985) to the true cluster
labels and the number of clusters ($\hat{K}_{n}$) for each method are given in
Table 1. Over all sample sizes, FOLD has the highest average ARI and lowest
number of clusters, consistently scoring above $0.95$ in the former and near
$3$ in the latter. While the VI is comparable in these metrics for $n\leq
500$, its ARI drops dramatically for larger $n$ as it prefers generally
$4$-$5$ clusters, causing it to achieve an ARI typically less than $0.7$.
Binder’s loss tends to substantially overestimate the number of clusters and
its ARI with the true cluster labels steadily decreases with increasing sample
size. While mclust tends to do better than Binder’s loss at controlling the
number of clusters, it attains the lowest average ARI across all methods for
larger sample size. Hence, FOLD excels over widely used model-based methods in
the misspecified kernel regime.

###  4.2 GSE81861 Cell Line Dataset

We apply FOLD to the GSE81861 cell line dataset (Li et al.,, 2017), which
measures single cell counts in $57,241$ genes from $630$ single-cell
transcriptomes. There are $7$ distinct cell lines present in the dataset, so
we compare FOLD and other model-based clustering methods to the true cell line
labels as a performance benchmark. First, we apply routine pre-processing
steps for RNA sequence analysis (as in Chandra et al., (2023)). We discard
cells with low read counts, giving $n=519$ total cells. We then normalize the
data using SCRAN (Lun et al.,, 2016) and select informative genes with M3Drop
(Andrews and Hemberg,, 2019). We use principal component analysis (PCA) for
dimension reduction by taking $\bm{X}^{0}$ to be the projection of the
normalized cell counts onto the first $p=5$ principal components, then scale
the projections to have zero mean and unit variance.

We fit a $p$-dimensional Bayesian Gaussian finite mixture to $\bm{X}^{0}$ with
$K=50$ components, and simulate $25,000$ posterior samples after a burn-in of
$1,000$. Every fourth iteration is discarded, leaving $6,000$ samples
remaining. Along with $\bm{c}_{\textrm{FOLD}}$, we calculate model-based
clusterings with the VI loss, Binder’s loss, and mclust. We also compare to
algorithmic clustering methods, including average linkage hierarchical
clustering (HCA), k-means, DBSCAN, and spectral clustering. Details on
hyperparameter choices for these methods are given in the Supplement. For
FOLD, we create candidate clusterings by average-linkage clustering on
$\Delta$, then choose $\omega$ by consulting the elbow plot in Figure 4. The
plot suggests $6$ clusters, which corresponds to $\omega=25$. We also
implement SALSO over $150$ independent runs to minimize (5). However, over
repeated replications we find that the SALSO estimate is either equal to or
has a higher risk than the point estimate given by hierarchical clustering on
$\Delta$. Binder’s loss and the VI loss are implemented with the default value
$\omega=1$, which results in $\bm{c}_{\textrm{VI}}=\bm{c}_{\textrm{B}}$.

Figure 4: Elbow plot for choosing the number of clusters in
$\bm{c}_{\textrm{FOLD}}$ with the cell line dataset.  Figure 5: UMAP plots of
the cell line dataset with colors corresponding to the true cell types,
$\bm{c}_{\textrm{FOLD}}$, VI/Binder, and mclust, respectively.

Figure 5 shows the UMAP plots (McInnes et al.,, 2018) of the original
normalized count data along with colors indicating the true cell types and the
clusterings from the three model-based methods. The ARI with the true cell
types and number of clusters for $\bm{c}_{\textrm{FOLD}}$ and all other
clusterings are given in Table 2. First, observe that the Bayesian clustering
methods outperform the competitors in the ARI, indicating that the overfitted
Gaussian components are accurately approximating the distribution within each
type. However, FOLD notably excels over the other Bayesian methods with an ARI
of $0.995$, only misclassifying six observations. Note that the VI, Binder’s
loss, and mclust split the GM12878 and H1 cell types into two clusters each,
while FOLD correctly identifies both types, keeping each as single clusters.
The splitting of types is not unexpected in this application since the GM12878
and H1 types were each sequenced in two separate batches (Li et al.,, 2017),
indicating that FOLD is more robust to the batch effect than other Bayesian
methods. However, $\bm{c}_{\textrm{FOLD}}$ underestimates the number of groups
by combining the H1437 and IMR90 cell types. The other model-based methods
merge these types as well, with the VI and Binder’s loss producing $9$
clusters and mclust giving $7$ clusters. Of the algorithmic methods, k-means
performs the best, estimating the correct number of types while achieving an
ARI of $0.904$. However, k-means splits the GM12878 and A549 types into two
clusters each, the latter of which is always kept as one cluster by the model-
based methods.

| FOLD | VI | Binder | Mclust | HCA | K-Means | DBSCAN | Spectral
---|---|---|---|---|---|---|---|---
ARI | 0.995 | 0.915 | 0.915 | 0.854 | 0.622 | 0.904 | 0.679 | 0.620
$\hat{K}_{n}$ | 6 | 9 | 9 | 7 | 7 | 7 | 5 | 5
Table 2: Adjusted Rand index (ARI) with the true cell lines and the number of
clusters ($\hat{K}_{n}$) for FOLD and competitors on the cell line dataset.
Figure 6: The clustering $\bm{c}_{\textrm{FOLD}}$ along with the horizontal,
vertical lower, and vertical upper bounds for the cell line dataset. Here,
$\mathbb{D}(\cdot,\cdot)$ is the VI and the bounds are unique.

The 95% credible ball around $\bm{c}_{\textrm{FOLD}}$ is displayed in Figure
6. The credible ball communicates that there is substantial uncertainty in
cell types where batch effects occur. The horizontal and vertical upper bounds
effectively merge the GM12878 type with the H1437 cell type, which is likely
due to the proximity of these two types in the sample space. Conversely, the
vertical lower bound splits the GM12878 cell type into two clusters. We are
also uncertain in our classification of the H1 cell type, which is similarly
split into multiple clusters by the vertical lower bound because of separate
batching. Interestingly, in all bounds, the H1437 and IMR90 cell types are
allocated to the same cluster. This type merging could be explained by the
fact that IMR90 consists of a small number of cells or that both types are
isolated from the lung (Li et al.,, 2017).

##  5 Discussion

Fusing of Localized Densities (FOLD) is a Bayesian method for cluster analysis
that characterizes clusters as possibly containing multiple mixture
components. We first detailed the decision theoretic justification behind our
approach, in which we obtain a point estimate of the clustering by minimizing
a novel loss function. Our loss function has several appealing properties,
such as favoring the merging of overlapping kernels, simplification to
Binder’s loss when all the mixture components are well separated, and
invariance to permutations of the labels. Uncertainty in cluster allocations
is expressed with a credible ball and posterior similarity matrix. We have
given concrete guidance on tuning the loss parameter, including an elbow plot
method and default value that performs excellently on simulated examples.

Throughout the article, we have primarily focused on the Gaussian mixture
model because of its ubiquity in the literature and useful theoretical
properties. However, FOLD can be applied to any parametric mixture in which a
bounded statistical distance between kernels is simple to compute. For the
Hellinger distance, this includes the beta, exponential, gamma, and Weibull
families. A near identical approach can be applied to discrete families, where
localized mass functions would replace the role of the localized densities.

We have shown that under regularity conditions, our clustering point estimator
converges to the oracle clustering procedure as $n\to\infty$ for both
misspecified and well-specified kernel regimes. The oracle rule categorizes
observations into clusters using knowledge about the KL-minimizer. We gave two
concrete examples of oracle rules, including the match rule and the merge
rule. Our findings imply similar behavior for Binder’s loss, meaning that
analogous rules could be potentially derived for other methods such as the VI.
These results would provide valuable insight into the overall asymptotic
performance of Bayesian clustering. In addition, our results could be extended
by focusing on contraction of the mixing measure in terms of the Voronoi loss
function (Manole and Ho,, 2022), a loss specifically intended for strongly
identifiable families that has been shown to contract at a faster rate.

Though we implemented FOLD with the Hellinger distance in our simulations and
application, other distribution metrics could be used instead. For a general
statistical distance $d$, one could set
$\mathcal{D}_{ij}=1-\exp\left\\{-d(g(\cdot;\theta_{i}),g(\cdot;\theta_{j}))\right\\}$.
The loss of co-clustering $i$ and $j$ is
$1-\exp\left\\{-d(g(\cdot;\theta_{i}),g(\cdot;\theta_{j}))\right\\}$. This
implies that we would favor co-clustering $i$ and $j$ when
$d(g(\cdot;\theta_{i}),g(\cdot;\theta_{j}))<-\log(1-\gamma)$. For
$\theta_{i}\in\mathbb{R}^{p}$, an even simpler variant is to choose
$\mathcal{D}_{ij}=1-\exp\left\\{-\rho\left\lVert\theta_{i}-\theta_{j}\right\rVert^{2}\right\\}$,
where $\rho>0$ is some fixed bandwidth, which would promote merging components
with similar atoms. Alternatively, one can view $\Delta$ as a stochastic
distance matrix and employ k-medoids, hierarchical clustering, or spectral
clustering to cluster the data. Finally, the general notion of specifying loss
functions for Bayesian clustering using $\bm{\theta}$ and not just $\bm{s}$
provides a promising direction for future research.

## Acknowledgements

This work was supported by the National Institutes of Health (NIH) under
Grants 1R01AI155733, 1R01ES035625, and 5R01ES027498; the United States Office
of Naval Research (ONR) Grant N00014-21-1-2510; and the European Research
Council (ERC) under Grant 856506.

## References

  * Andrews and Hemberg, (2019) Andrews, T. S. and Hemberg, M. (2019).  M3Drop: dropout-based feature selection for scRNAseq.  Bioinformatics, 35(16):2865–2867.
  * Antoniak, (1974) Antoniak, C. E. (1974).  Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems.  The Annals of Statistics, 2(6):1152–1174.
  * Aragam et al., (2020) Aragam, B., Dan, C., Xing, E. P., and Ravikumar, P. (2020).  Identifiability of nonparametric mixture models and Bayes optimal clustering.  The Annals of Statistics, 48(4):2277–2302.
  * Ascolani et al., (2023) Ascolani, F., Lijoi, A., Rebaudo, G., and Zanella, G. (2023).  Clustering consistency with Dirichlet process mixtures.  Biometrika, 110(2):551–558.
  * Azzalini and Capitanio, (1999) Azzalini, A. and Capitanio, A. (1999).  Statistical applications of the multivariate skew normal distribution.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3):579–602.
  * Azzalini and Valle, (1996) Azzalini, A. and Valle, A. D. (1996).  The multivariate skew-normal distribution.  Biometrika, 83(4):715–726.
  * Barndorff-Nielsen, (1965) Barndorff-Nielsen, O. (1965).  Identifiability of mixtures of exponential families.  Journal of Mathematical Analysis and Applications, 12(1):115–121.
  * Bartolucci, (2005) Bartolucci, F. (2005).  Clustering univariate observations via mixtures of unimodal normal mixtures.  Journal of Classification, 22(2):203–219.
  * Baudry et al., (2010) Baudry, J.-P., Raftery, A. E., Celeux, G., Lo, K., and Gottardo, R. (2010).  Combining mixture components for clustering.  Journal of Computational and Graphical Statistics, 19(2):332–353.
  * Bensmail et al., (1997) Bensmail, H., Celeux, G., Raftery, A. E., and Robert, C. P. (1997).  Inference in model-based cluster analysis.  Statistics and Computing, 7(1):1–10.
  * Berenfeld et al., (2022) Berenfeld, C., Rosa, P., and Rousseau, J. (2022).  Estimating a density near an unknown manifold: a Bayesian nonparametric approach.  arXiv preprint arXiv:2205.15717.
  * Binder, (1978) Binder, D. A. (1978).  Bayesian cluster analysis.  Biometrika, 65(1):31–38.
  * Cacciatore and Tenori, (2022) Cacciatore, S. and Tenori, L. (2022).  KODAMA: Knowledge discovery by accuracy maximization.  R package version 1.9.
  * Cai et al., (2021) Cai, D., Campbell, T., and Broderick, T. (2021).  Finite mixture models do not reliably learn the number of components.  In Meila, M. and Zhang, T., editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 1158–1169. PMLR.
  * Chan et al., (2008) Chan, C., Feng, F., Ottinger, J., Foster, D., West, M., and Kepler, T. B. (2008).  Statistical mixture modeling for cell subtype identification in flow cytometry.  Cytometry Part A: The Journal of the International Society for Analytical Cytology, 73(8):693–701.
  * Chandra et al., (2023) Chandra, N. K., Canale, A., and Dunson, D. B. (2023).  Escaping the curse of dimensionality in Bayesian model-based clustering.  Journal of Machine Learning Research, 24(144):1–42.
  * Chen, (1995) Chen, J. (1995).  Optimal rate of convergence for finite mixture models.  The Annals of Statistics, 23(1):221–233.
  * Dahl et al., (2022) Dahl, D. B., Johnson, D. J., and Müller, P. (2022).  Search algorithms and loss functions for Bayesian clustering.  Journal of Computational and Graphical Statistics, 31(4):1189–1201.
  * Dang et al., (2023) Dang, U. J., Gallaugher, M. P., Browne, R. P., and McNicholas, P. D. (2023).  Model-based clustering and classification using mixtures of multivariate skewed power exponential distributions.  Journal of Classification, pages 1–23.
  * De Iorio et al., (2023) De Iorio, M., Favaro, S., Guglielmi, A., and Ye, L. (2023).  Bayesian nonparametric mixture modeling for temporal dynamics of gender stereotypes.  The Annals of Applied Statistics, 17(3):2256–2278.
  * Denti et al., (2023) Denti, F., Camerlenghi, F., Guindani, M., and Mira, A. (2023).  A common atoms model for the Bayesian nonparametric analysis of nested data.  Journal of the American Statistical Association, 118(541):405–416.
  * Di Zio et al., (2007) Di Zio, M., Guarnera, U., and Rocci, R. (2007).  A mixture of mixture models for a classification problem: The unity measure error.  Computational Statistics & Data Analysis, 51(5):2573–2585.
  * Duan and Dunson, (2021) Duan, L. L. and Dunson, D. B. (2021).  Bayesian distance clustering.  Journal of Machine Learning Research, 22(224):1–27.
  * Ferguson, (1973) Ferguson, T. S. (1973).  A Bayesian analysis of some nonparametric problems.  The Annals of Statistics, 1(2):209–230.
  * Fraley and Raftery, (2002) Fraley, C. and Raftery, A. E. (2002).  Model-based clustering, discriminant analysis, and density estimation.  Journal of the American Statistical Association, 97(458):611–631.
  * Fritsch, (2022) Fritsch, A. (2022).  Mcclust: Process an MCMC sample of clusterings.  R package version 1.0.1.
  * Fritsch and Ickstadt, (2009) Fritsch, A. and Ickstadt, K. (2009).  Improved criteria for clustering based on the posterior similarity matrix.  Bayesian Analysis, 4(2):367–391.
  * Ghosal et al., (1999) Ghosal, S., Ghosh, J. K., and Ramamoorthi, R. (1999).  Posterior consistency of Dirichlet mixtures in density estimation.  The Annals of Statistics, 27(1):143–158.
  * Ghosal and van der Vaart, (2007) Ghosal, S. and van der Vaart, A. (2007).  Posterior convergence rates of Dirichlet mixtures at smooth densities.  The Annals of Statistics, 35(2):697–723.
  * Gorsky et al., (2024) Gorsky, S., Chan, C., and Ma, L. (2024).  Coarsened mixtures of hierarchical skew normal kernels for flow cytometry analyses.  Bayesian Analysis, Forthcoming.
  * Guha et al., (2021) Guha, A., Ho, N., and Nguyen, X. (2021).  On posterior contraction of parameters and interpretability in Bayesian mixture modeling.  Bernoulli, 27(4):2159–2188.
  * Hahsler et al., (2019) Hahsler, M., Piekenbrock, M., and Doran, D. (2019).  dbscan: Fast density-based clustering with R.  Journal of Statistical Software, 91:1–30.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., and Friedman, J. H. (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction, volume 2.  Springer.
  * Heinrich and Kahn, (2018) Heinrich, P. and Kahn, J. (2018).  Strong identifiability and optimal minimax rates for finite mixture estimation.  The Annals of Statistics, 46(6A):2844–2870.
  * Hennig, (2010) Hennig, C. (2010).  Methods for merging Gaussian mixture components.  Advances in Data Analysis and Classification, 4(1):3–34.
  * (36) Ho, N. and Nguyen, X. (2016a).  Convergence rates of parameter estimation for some weakly identifiable finite mixtures.  The Annals of Statistics, 44(6):2726–2755.
  * (37) Ho, N. and Nguyen, X. (2016b).  On strong identifiability and convergence rates of parameter estimation in finite mixtures.  Electronic Journal of Statistics, 10(1):271–307.
  * Ho et al., (2020) Ho, N., Nguyen, X., and Ritov, Y. (2020).  Robust estimation of mixing measures in finite mixture models.  Bernoulli, 26(2):828–857.
  * Hubert and Arabie, (1985) Hubert, L. and Arabie, P. (1985).  Comparing partitions.  Journal of Classification, 2(1):193–218.
  * Ishwaran and James, (2001) Ishwaran, H. and James, L. F. (2001).  Gibbs sampling methods for stick-breaking priors.  Journal of the American Statistical Association, 96(453):161–173.
  * Ishwaran and James, (2003) Ishwaran, H. and James, L. F. (2003).  Some further developments for stick-breaking priors: Finite and infinite clustering and classification.  Sankhyā: The Indian Journal of Statistics (2003-2007), 65(3):577–592.
  * Jain, (2010) Jain, A. K. (2010).  Data clustering: 50 years beyond K-means.  Pattern Recognition Letters, 31(8):651–666.
  * Karlis and Santourian, (2009) Karlis, D. and Santourian, A. (2009).  Model-based clustering with non-elliptically contoured distributions.  Statistics and Computing, 19(1):73–83.
  * Kiselev et al., (2019) Kiselev, V. Y., Andrews, T. S., and Hemberg, M. (2019).  Challenges in unsupervised clustering of single-cell RNA-seq data.  Nature Reviews Genetics, 20(5):273–282.
  * Krijthe, (2016) Krijthe, J. H. (2016).  RSSL: R package for semi-supervised learning.  In Reproducible Research in Pattern Recognition. RRPR 2016. Lecture Notes in Computer Science, vol 10214, pages 104–115.
  * Krijthe and Loog, (2015) Krijthe, J. H. and Loog, M. (2015).  Implicitly constrained semi-supervised least squares classification.  In 14th International Symposium on Advances in Intelligent Data Analysis XIV (Lecture Notes in Computer Science Volume 9385), pages 158–169.
  * Lau and Green, (2007) Lau, J. W. and Green, P. J. (2007).  Bayesian model-based clustering procedures.  Journal of Computational and Graphical Statistics, 16(3):526–558.
  * Li et al., (2017) Li, H., Courtois, E. T., Sengupta, D., Tan, Y., Chen, K. H., Goh, J. J. L., Kong, S. L., Chua, C., Hon, L. K., Tan, W. S., et al. (2017).  Reference component analysis of single-cell transcriptomes elucidates cellular heterogeneity in human colorectal tumors.  Nature Genetics, 49(5):708–718.  PMID: 28319088.
  * Li, (2005) Li, J. (2005).  Clustering based on a multilayer mixture model.  Journal of Computational and Graphical Statistics, 14(3):547–568.
  * Lun et al., (2016) Lun, A. T., Bach, K., and Marioni, J. C. (2016).  Pooling across cells to normalize single-cell RNA sequencing data with many zero counts.  Genome Biology, 17(1):1–14.
  * Malsiner-Walli et al., (2017) Malsiner-Walli, G., Frühwirth-Schnatter, S., and Grün, B. (2017).  Identifying mixtures of mixtures using Bayesian estimation.  Journal of Computational and Graphical Statistics, 26(2):285–295.
  * Manole and Ho, (2022) Manole, T. and Ho, N. (2022).  Refined convergence rates for maximum likelihood estimation under finite mixture models.  In Chaudhuri, K., Jegelka, S., Song, L., Szepesvari, C., Niu, G., and Sabato, S., editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 14979–15006. PMLR.
  * Manole and Khalili, (2021) Manole, T. and Khalili, A. (2021).  Estimating the number of components in finite mixture models via the group-sort-fuse procedure.  The Annals of Statistics, 49(6):3043–3069.
  * McInnes et al., (2018) McInnes, L., Healy, J., Saul, N., and Großberger, L. (2018).  UMAP: Uniform manifold approximation and projection.  Journal of Open Source Software, 3(29):861.
  * Medvedovic and Sivaganesan, (2002) Medvedovic, M. and Sivaganesan, S. (2002).  Bayesian infinite mixture model based clustering of gene expression profiles.  Bioinformatics, 18(9):1194–1206.
  * Meilă, (2007) Meilă, M. (2007).  Comparing clusterings—an information based distance.  Journal of Multivariate Analysis, 98(5):873–895.
  * Melnykov, (2016) Melnykov, V. (2016).  Merging mixture components for clustering through pairwise overlap.  Journal of Computational and Graphical Statistics, 25(1):66–90.
  * Miller and Dunson, (2019) Miller, J. W. and Dunson, D. B. (2019).  Robust Bayesian inference via coarsening.  Journal of the American Statistical Association, 114(527):1113–1125.
  * Miller and Harrison, (2014) Miller, J. W. and Harrison, M. T. (2014).  Inconsistency of Pitman-Yor process mixtures for the number of components.  The Journal of Machine Learning Research, 15(1):3333–3370.
  * Miller and Harrison, (2018) Miller, J. W. and Harrison, M. T. (2018).  Mixture models with a prior on the number of components.  Journal of the American Statistical Association, 113(521):340–356.
  * Nguyen, (2011) Nguyen, X. (2011).  Wasserstein distances for discrete measures and convergence in nonparametric mixture models.  Technical report, Univ. of Michigan Dept. of Statistics.
  * Nguyen, (2013) Nguyen, X. (2013).  Convergence of latent mixing measures in finite and infinite mixture models.  The Annals of Statistics, 41(1):370–400.
  * O’Hagan et al., (2016) O’Hagan, A., Murphy, T. B., Gormley, I. C., McNicholas, P. D., and Karlis, D. (2016).  Clustering with the multivariate normal inverse Gaussian distribution.  Computational Statistics & Data Analysis, 93:18–30.
  * Paganin et al., (2021) Paganin, S., Herring, A. H., Olshan, A. F., and Dunson, D. B. (2021).  Centered partition processes: Informative priors for clustering (with discussion).  Bayesian Analysis, 16(1):301 – 670.
  * Rajkowski, (2019) Rajkowski, Ł. (2019).  Analysis of the maximal a posteriori partition in the Gaussian Dirichlet process mixture model.  Bayesian Analysis, 14(2):477–494.
  * Rand, (1971) Rand, W. M. (1971).  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical Association, 66(336):846–850.
  * Rastelli and Friel, (2018) Rastelli, R. and Friel, N. (2018).  Optimal Bayesian estimators for latent variable cluster models.  Statistics and Computing, 28(6):1169–1186.
  * Redner and Walker, (1984) Redner, R. A. and Walker, H. F. (1984).  Mixture densities, maximum likelihood and the EM algorithm.  SIAM Review, 26(2):195–239.
  * Rigon et al., (2023) Rigon, T., Herring, A. H., and Dunson, D. B. (2023).  A generalized Bayes framework for probabilistic clustering.  Biometrika, 110(3):559–578.
  * Rodríguez and Walker, (2014) Rodríguez, C. E. and Walker, S. G. (2014).  Univariate Bayesian nonparametric mixture modeling with unimodal kernels.  Statistics and Computing, 24(1):35–49.
  * Rousseau and Mengersen, (2011) Rousseau, J. and Mengersen, K. (2011).  Asymptotic behaviour of the posterior distribution in overfitted mixture models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(5):689–710.
  * Scrucca et al., (2016) Scrucca, L., Fop, M., Murphy, T. B., and Raftery, A. E. (2016).  Mclust 5: Clustering, classification and density estimation using Gaussian finite mixture models.  The R Journal, 8(1):289–317.
  * Sethuraman, (1994) Sethuraman, J. (1994).  A constructive definition of Dirichlet priors.  Statistica Sinica, 4(2):639–650.
  * Teicher, (1961) Teicher, H. (1961).  Identifiability of mixtures.  The Annals of Mathematical Statistics, 32(1):244–248.
  * Wade, (2015) Wade, S. (2015).  Mcclust.ext: Point estimation and credible balls for Bayesian cluster analysis.  R package version 1.0.
  * Wade, (2023) Wade, S. (2023).  Bayesian cluster analysis.  Philosophical Transactions of the Royal Society A, 381(2247):20220149.
  * Wade and Ghahramani, (2018) Wade, S. and Ghahramani, Z. (2018).  Bayesian cluster analysis: Point estimation and credible balls (with discussion).  Bayesian Analysis, 13(2):559–626.

## Supplementary Material

The Supplementary Material contains additional simulation studies, proofs for
all results, and technical details. Code to implement FOLD is available at
https://github.com/adombowsky/FOLD.

##  Appendix A Point Estimation Simulation Studies

In light of the asymptotic results shown in Section 3, we examine the effect
of increasing $n$ on FOLD and other model-based clustering methods. For a
given $f^{0}$, we repeatedly simulate observations with varying
$n\in\left\\{100,500,1000,2500\right\\}$ for $100$ replications each. We
compare FOLD with $\omega=\omega^{\textrm{AVG}}$ to clusterings returned by
minimizing Binder’s loss and the VI loss using the packages mcclust (Fritsch,,
2022) and mcclust.ext (Wade,, 2015), respectively. Along with these Bayesian
approaches, we also cluster the data using mclust (Scrucca et al.,, 2016),
which groups observations based on the EM algorithm for a Gaussian mixture.

For each replication and $n$, we run the EM algorithm and fit a Bayesian
location-scale GMM,

$\displaystyle X_{i}^{0}\mid\bm{a},\bm{\mu},\bm{\Sigma}$
$\displaystyle\sim\sum_{k=1}^{K}a_{l}\mathcal{N}_{2}(\tilde{\mu}_{k},\tilde{\Sigma}_{k});$
$\displaystyle(\tilde{\mu}_{k},\tilde{\Sigma}_{k})$
$\displaystyle\sim\mathcal{NIW}(\mu,\kappa,\nu,\Psi);$ $\displaystyle\bm{a}$
$\displaystyle\sim\textrm{Dir}\left(\alpha 1_{K}\right);$

where $1_{K}$ is the $K$-dimensional vector consisting of ones. We then
compute clusterings with FOLD, VI, and Binder’s loss using the posterior
samples from the Bayesian model. For each clustering, we save the number of
clusters and the adjusted Rand index (Rand,, 1971; Hubert and Arabie,, 1985)
with $\bm{s}^{0}$. The hyperparameters are set at $K=30$, $\alpha=1/2$, and
$\mu=0_{p}$, $\kappa=1$, $\nu=p+2$, and $\Psi=I$. We run a Gibbs sampler for
$9,000$ iterations with a burn-in of $1,000$, then use every third iteration
for computing the Bayesian clusterings. For mclust, the number of clusters is
automatically selected each replication by minimizing the Bayesian information
criterion (BIC).

Figure 7: Contour plots of the (a) Gaussian mixture, (b) skew Gaussian mixture and (c) skew-symmetric mixture. | $n$ | FOLD | Oracle FOLD | VI | Binder’s | Mclust
---|---|---|---|---|---|---
ARI | 100 | 0.903 (0.087) | 0.985 (0.021) | 0.852 (0.088) | 0.822 (0.061) | 0.962 (0.052)
| 500 | 0.979 (0.013) | 0.998 (0.008) | 0.974 (0.014) | 0.947 (0.012) | 0.987 (0.008)
| 1000 | 0.985 (0.006) | 0.987 (0.006) | 0.982 (0.007) | 0.968 (0.007) | 0.985 (0.024)
| 2500 | 0.986 (0.004) | 0.987 (0.004) | 0.984 (0.004) | 0.979 (0.005) | 0.987 (0.004)
$\hat{K}_{n}$ | 100 | 2.93 (0.383) | 3.00 (0.00) | 8.02 (5.714) | 19.76 (0.818) | 3.10 (0.414)
| 500 | 3.06 (0.278) | 3.00 (0.00) | 3.53 (0.926) | 19.89 (0.852) | 3.00 (0.00)
| 1000 | 3.05 (0.219) | 3.00 (0.00) | 3.30 (0.541) | 19.96 (0.400) | 3.01 (0.100)
| 2500 | 3.05 (0.219) | 3.00 (0.00) | 3.08 (0.273) | 20.00 (0.00) | 3.00 (0.00)
Table 3: Averages and standard deviations (in parentheses) for the adjusted Rand index and number of clusters on 100 replications from a mixture of bivariate Gaussian kernels. | $n$ | FOLD | VI | Binder’s | Mclust
---|---|---|---|---|---
ARI | 100 | 0.992 (0.016) | 0.987 (0.021) | 0.915 (0.043) | 0.989 (0.030)
| 500 | 0.999 (0.003) | 0.998 (0.003) | 0.992 (0.006) | 0.980 (0.049)
| 1000 | 0.999 (0.002) | 0.997 (0.007) | 0.990 (0.009) | 0.866 (0.149)
| 2500 | 0.999 (0.001) | 0.990 (0.020) | 0.952 (0.014) | 0.576 (0.071)
$\hat{K}_{n}$ | 100 | 3.02 (0.141) | 3.16 (0.395) | 13.84 (6.080) | 3.14 (0.377)
| 500 | 3.03 (0.171) | 3.42 (0.755) | 6.15 (2.258) | 3.23 (0.529)
| 1000 | 3.20 (0.449) | 3.53 (0.881) | 9.32 (3.816) | 4.33 (1.429)
| 2500 | 3.19 (0.443) | 3.40 (0.682) | 8.93 (2.409) | 7.99 (0.959)
Table 4: Averages and standard deviations (in parentheses) for the adjusted Rand index and number of clusters on 100 replications from a mixture of bivariate skew Gaussian kernels. | $n$ | FOLD | VI | Binder’s | Mclust
---|---|---|---|---|---
ARI | 100 | 0.982 (0.029) | 0.974 (0.026) | 0.800 (0.094) | 0.884 (0.177)
| 500 | 0.995 (0.032) | 0.921 (0.136) | 0.799 (0.147) | 0.715 (0.159)
| 1000 | 0.975 (0.085) | 0.691 (0.057) | 0.665 (0.037) | 0.574 (0.054)
| 2500 | 0.967 (0.098) | 0.679 (0.012) | 0.639 (0.027) | 0.512 (0.056)
$\hat{K}_{n}$ | 100 | 3.25 (0.500) | 3.89 (1.014) | 12.62 (4.292) | 3.56 (0.925)
| 500 | 3.23 (0.601) | 3.64 (0.894) | 13.94 (4.552) | 4.17 (0.779)
| 1000 | 3.10 (0.302) | 4.83 (1.025) | 15.73 (3.795) | 5.08 (0.442)
| 2500 | 3.10 (0.302) | 4.58 (0.843) | 18.99 (2.634) | 6.20 (0.876)
Table 5: Averages and standard deviations (in parentheses) for the adjusted
Rand index and number of clusters on 100 replications from the skew-symmetric
mixture. Figure 8: Comparison of the number of clusters and the adjusted Rand
index with $\bm{s}^{0}$ on $100$ replications from a mixture of well-separated
multivariate Gaussian kernels. Figure 9: Comparison of the number of clusters
and the adjusted Rand index with $\bm{s}^{0}$ on $100$ replications from a
mixture of multivariate skew Gaussian distributions. Figure 10: Comparison of
the number of clusters and the adjusted Rand index with $\bm{s}^{0}$ on $100$
replications from the skew-symmetric mixture.

In the first scenario,
$f^{0}(x)=\sum_{m=1}^{3}a_{m}^{0}\mathcal{N}_{2}(x;\mu_{m}^{0},\Sigma_{m}^{0})$,
with weights $\bm{a}^{0}=(0.45,0.25,0.3)$ and atoms $\mu_{1}^{0}=(6.5,5)$,
$\Sigma_{1}^{0}=I$, $\mu_{2}^{0}=(0,0)$, $\Sigma_{2}^{0}=\textrm{diag}(5,2)$,
and $\mu_{3}^{0}=(-5,-5)$, $\Sigma_{3}^{0}=\textrm{diag}(3,1)$. As shown in
Figure 7(a), the three kernels are well separated. The averages and standard
deviations for the number of clusters and adjusted Rand index are given in
Table 3, with boxplots of these benchmarks displayed in Figure 8. FOLD, oracle
FOLD, and mclust achieve high adjusted Rand index with increasing sample size.
Oracle FOLD is perfect at computing the true number of clusters, and mclust
and FOLD are also excellent. Furthermore, the number of clusters and adjusted
Rand index for $\bm{c}_{\textrm{FOLD}}$ gradually approach the values produced
by $\bm{c}_{\textrm{FOLD}}^{*}$. The VI loss improves in the number of
clusters as $n$ increases, but Binder’s loss falls short, consistently
producing $20$ clusters with diminishing variation across replications. VI and
Binder’s loss produce high adjusted Rand indices that improve with $n$.

Next, we let
$f^{0}(x)=\sum_{m=1}^{3}a_{m}^{0}\mathcal{SN}_{2}(x;\mu_{m}^{0},\Sigma_{m}^{0},\psi_{m}^{0})$,
where $\mathcal{SN}_{2}(x;\cdot,\cdot,\cdot)$ denotes the PDF of a bivariate
skew Gaussian distribution (Azzalini and Valle,, 1996; Azzalini and
Capitanio,, 1999) and $\bm{a}^{0}=(0.45,0.25,0.3).$ These kernels have the
same location and scale parameters as the bivariate Gaussian mixture in
simulation case 1, and skewness parameters $\psi_{1}^{0}=(1,1)$,
$\psi_{2}^{0}=(-10,15)$, and $\psi_{3}^{0}=(4,-17)$. Figure 7(b) shows a
contour plot of $f^{0}$. Results across the replications are reported in
Figure 9 and Table 4. For all values of $n$, FOLD attains high accuracy with a
small number of clusters. As $n$ increases, all four clustering methods
achieve high levels of the adjusted Rand index, though FOLD results in values
close to $1$ while the VI loss and Binder’s loss decline slightly. FOLD tends
to produce between $3$ and $4$ clusters for each $n$ and reports less than or
equal to the number of clusters induced by the VI loss in $94.5\%$ of all
instances.

In the final simulation case, we take $f^{0}$ to be a bivariate mixture with
weights $\bm{a}^{0}=(0.55,0.3,0.15)$, and kernels

$\displaystyle
g_{1}^{0}=0.364\cdot\mathcal{SN}_{p}\left(\begin{pmatrix}2.50\\\
3.50\end{pmatrix},I,\begin{pmatrix}-10\\\ 15\end{pmatrix}\right)$
$\displaystyle+0.212\cdot\mathcal{N}_{p}\left(\begin{pmatrix}2.325\\\
4.381\end{pmatrix},\textrm{diag}(0.20,0.80)\right)$
$\displaystyle+0.424\cdot\mathcal{N}_{p}\left(\begin{pmatrix}1.085\\\
2.009\end{pmatrix},\textrm{diag}(0.70,0.60)\right);$ $\displaystyle
g_{2}^{0}=\mathcal{SN}_{p}\left(\begin{pmatrix}0\\\
-3.50\end{pmatrix},\textrm{diag}(5,2),\begin{pmatrix}4\\\
-17\end{pmatrix}\right);$ $\displaystyle
g_{3}^{0}=\mathcal{N}_{p}\left(\begin{pmatrix}-4\\\
-2.50\end{pmatrix},\begin{pmatrix}0.50&0.50\\\ 0.50&2.50\end{pmatrix}\right).$

$g_{1}^{0}$ is a mixture of one skew Gaussian kernel and two Gaussian kernels,
$g_{2}^{0}$ is a skew Gaussian kernel, and $g_{3}^{0}$ is a Gaussian kernel.
We refer to $f^{0}$ in this case as a skew-symmetric mixture. As can be seen
in Figure 7(c), $g_{1}^{0}$ is a non-Gaussian, multimodal density. Similarly,
data generated from $g_{2}^{0}$, the oblong kernel in the lower third of the
sample space, will most likely be approximated by multiple components, despite
only constituting one cluster. The results of the simulations are summarized
in Table 5 and Figure 10. FOLD generally allocates observations to $3$
clusters, but the VI loss favors between $3$ and $6$ clusters. In $96.5\%$ of
the replicates, the number of clusters for FOLD is less than or equal to the
number for VI. As in simulation case 2, Binder’s loss and mclust frequently
return a larger number of clusters than the truth for larger $n$. The adjusted
Rand index between $\bm{c}_{\textrm{FOLD}}$ and $\bm{s}^{0}$ is close to $1$
across all sample sizes. In contrast, the three other methods achieve high
adjusted Rand index for small $n$, but these values sharply drop for $n\geq
1000$. This is usually the result of splitting $g_{1}^{0}$ or $g_{2}^{0}$ into
multiple clusters.

##  Appendix B Illustration of Theorem 1: Location GMMs

We now empirically validate the consistency result in Theorem 1 with a simple
example. We set $M^{0}=M^{*}=4$ and
$f^{0}=f^{*}=\sum_{m=1}^{4}(1/4)\mathcal{N}_{2}(\mu_{m}^{0},(1/4)I)$, where
$\mu_{1}^{0}=(1,1)$, $\mu_{2}^{0}=(1.75,1.75)$, $\mu_{3}^{0}=(-1.75,-1.75)$,
and $\mu_{4}^{0}=(-1,-1)$, then fit the following location Gaussian model with
$K=4$ components,

$\displaystyle X_{i}$ $\displaystyle\sim\mathcal{N}_{2}(\theta_{i},(1/4)I);$
$\displaystyle\tilde{\theta}_{k}$ $\displaystyle\sim\mathcal{N}_{2}(0,I);$
$\displaystyle\bm{a}$ $\displaystyle\sim\textrm{Dir}(1/4,\dots,1/4).$ (11)

We vary $n\in\left\\{50,100,500,1000\right\\}$, then fit a Gibbs sampler to
extract samples from $\pi(\bm{\theta},\bm{s}\mid\bm{X})$. We implement FOLD,
oracle FOLD, Binder’s loss, and Binder’s oracle loss using SALSO to minimize
the risk functions, and we set $d(\cdot,\cdot)$ to be the Hellinger distance.
For FOLD, we set $\omega=\omega^{\textrm{AVG}}$ and for the FOLD oracle, we
take $\omega=\omega^{*\textrm{AVG}}$, which is constructed by setting
$\gamma=\bar{\Delta}^{*}$ in the oracle risk. For Binder’s loss, we take
$\omega=1$, matching the choice used in packages such as mcclust. Figure 11
shows $\mathcal{R}(\widehat{\bm{c}})$ and $\mathcal{R}^{*}(\widehat{\bm{c}})$
for increasing $n$, where $\widehat{\bm{c}}$ are taken to be clusterings
generated by average linkage hierarchical clustering on $\Delta$. We can see
that as $n$ grows,
$\mathcal{R}(\widehat{\bm{c}})\to\mathcal{R}^{*}(\widehat{\bm{c}})$. In all
plots, the minimum risk over the candidates is attained at $2$ clusters,
though for larger $n$ the risk of $3$ clusters becomes similar to that of $2$
clusters. In addition, the point estimates $\bm{c}_{\textrm{FOLD}}$ and
$\bm{c}_{\textrm{FOLD}}^{*}$ are very similar as well even for smaller sample
size, with adjusted Rand indicies between the two partitions being $0.98$ for
$n=50$, then $1.0$ for $n=100,500,1000$. At $n=50$, $\bm{c}_{\textrm{FOLD}}$
has $3$ clusters, but for $n>50$ both $\bm{c}_{\textrm{FOLD}}$ and
$\bm{c}_{\textrm{FOLD}}^{*}$ have exactly $2$ clusters. Figure 12 shows point
estimates for $\bm{c}_{\textrm{FOLD}}$ and $\bm{c}_{\textrm{FOLD}}^{*}$ along
with $\bm{c}_{\textrm{B}}$ and $\bm{c}_{\textrm{B}}^{*}$ at the largest sample
size. Binder’s loss and its oracle prefer clusterings that approximately
correspond to the original $4$ Gaussian components, while FOLD and its oracle
prefer to merge these. Additionally, several small clusters are present in
Binder’s clusterings which seem to exist in the tails of the Gaussian
components, leading to $5$ clusters in $\bm{c}_{\textrm{B}}$ and $7$ clusters
in $\bm{c}_{\textrm{B}}^{*}$.

Figure 11: $\mathcal{R}(\widehat{\bm{c}})$ in red and
$\mathcal{R}^{*}(\widehat{\bm{c}})$ in green over a set of partitions produced
by average linkage hierarchical clustering on $\Delta$ and ordered by the
number of clusters for varying $n\in\left\\{100,250,500,1000\right\\}$. As $n$
grows, the two risk functions and their minimizers begin to closesly resemble
each other. Figure 12: $\bm{c}_{\textrm{FOLD}}$, $\bm{c}_{\textrm{FOLD}}^{*}$,
$\bm{c}_{\textrm{B}}$, and $\bm{c}_{\textrm{B}}^{*}$ for $n=1000$ simulated
data points from a mixture of $4$ Gaussian kernels. The number of clusters for
both $\bm{c}_{\textrm{FOLD}}$ and $\bm{c}_{\textrm{FOLD}}^{*}$ is $2$, whereas
for $\bm{c}_{\textrm{B}}$ there are $5$ clusters and for
$\bm{c}_{\textrm{B}}^{*}$ there are $7$ clusters.

To better understand the results in Figure 12, note that the Hellinger
distance between the kernels in $f^{*}$ is directly related to the Euclidean
distance between their centers,

$\displaystyle
1-d^{2}\left\\{\mathcal{N}_{2}(\cdot;\theta_{m}^{0},1/4),\mathcal{N}_{2}(\cdot;\theta_{m^{\prime}}^{0},1/4)\right\\}\propto\exp\left\\{-\frac{1}{2}\left\lVert\theta_{m}^{0}-\theta_{m^{\prime}}^{0}\right\rVert^{2}\right\\}.$

We have deliberately chosen the true component means so that $\theta_{1}^{0}$
and $\theta_{2}^{0}$ are close, and $\theta_{3}^{0}$ and $\theta_{4}^{0}$ are
close. This means that, despite the fact that we originally sampled from $4$
mixture components, the scatterplot of the data in Figure 12 resembles $2$
ovular clusters rather than $4$ spherical ones. $\bm{c}_{\textrm{FOLD}}^{*}$
creates ovular clusters because it is less reliant on the fitted model in (11)
than Binder’s loss. For example, if $X_{i}^{0}$ is in a high density region of
$\mathcal{N}_{2}(\cdot;\theta_{1}^{0},1/4)$ and $X_{j}^{0}$ is in a high
density region of $\mathcal{N}_{2}(\cdot;\theta_{2}^{0},1/4)$, then the
Hellinger distance terms in (10) shrink the value of $\Delta_{ij}^{*}$ close
to $0$, ultimately promoting co-clustering of objects $i$ and $j$. In
contrast, Binder’s oracle is more reliant on assigning the observations to
mixture components and would most likely place $X_{i}^{0}$ and $X_{j}^{0}$ in
different clusters.

##  Appendix C Basic Properties of FOLD

###  C.1 Proof of Proposition 1

###### Proof.

We proceed as in Dahl et al., (2022) and decompose the loss function into
separate terms that relate to the contingency table between $\widehat{\bm{c}}$
and $\bm{s}$. For any $k,k^{\prime}\in[K]$, let
$\eta_{kk^{\prime}}=d(g(\cdot;\tilde{\theta}_{k}),g(\cdot;\tilde{\theta}_{k^{\prime}}))$.
For index sets $\hat{C}_{h}$ and $S_{k}$, set $n_{h\cdot}=|\hat{C}_{h}|$,
$n_{\cdot k}=|S_{k}|$, and $n_{hk}=|\hat{C}_{h}\cap S_{k}|$.

$\displaystyle\mathcal{L}(\widehat{\bm{c}},\bm{\theta})=\sum_{i<j}\left\\{\textbf{1}_{\hat{c}_{i}=\hat{c}_{j}}\mathcal{D}_{ij}+\omega\textbf{1}_{\hat{c}_{i}\neq\hat{c}_{j}}\left(1-\mathcal{D}_{ij}\right)\right\\}$
$\displaystyle=\sum_{i<j}\left\\{\textbf{1}_{\hat{c}_{i}=\hat{c}_{j}}\left(1-\left(1-\mathcal{D}_{ij}\right)\right)+\omega\left(1-\textbf{1}_{\hat{c}_{i}=\hat{c}_{j}}\right)\left(1-\mathcal{D}_{ij}\right)\right\\}$
$\displaystyle=\sum_{i<j}\textbf{1}(\hat{c}_{i}=\hat{c}_{j})+\omega\sum_{i<j}\left(1-\mathcal{D}_{ij}\right)-(1+\omega)\sum_{i<j}\textbf{1}_{\hat{c}_{i}=\hat{c}_{j}}\left(1-\mathcal{D}_{ij}\right)$
$\displaystyle=\sum_{h=1}^{\hat{K}_{n}}{n_{h\cdot}\choose
2}+\omega\sum_{k=1}^{K}{n_{\cdot k}\choose
2}+\omega\sum_{k<k^{\prime}}n_{\cdot k}n_{\cdot
k^{\prime}}\left(1-\eta_{kk^{\prime}}\right)$
$\displaystyle-(1+\omega)\sum_{h=1}^{\hat{K}_{n}}\sum_{k=1}^{K}{n_{hk}\choose
2}-(1+\omega)\sum_{h=1}^{\hat{K}_{n}}\sum_{k<k^{\prime}}n_{hk}n_{hk^{\prime}}\left(1-\eta_{kk^{\prime}}\right).$

Consequently, this shows that the loss is related to the Binder’s loss with
unit costs $a=1$ and $b=\omega$,

$\displaystyle\mathcal{L}(\widehat{\bm{c}},\bm{\theta})=\mathcal{L}_{\textrm{B}}(\widehat{\bm{c}},\bm{s})+\mathcal{B}(\widehat{\bm{c}},\tilde{\bm{\theta}})$

where

$\displaystyle\mathcal{B}(\widehat{\bm{c}},\tilde{\bm{\theta}})$
$\displaystyle=\omega\sum_{k<k^{\prime}}n_{\cdot k}n_{\cdot
k^{\prime}}\left(1-\eta_{kk^{\prime}}\right)$
$\displaystyle-(1+\omega)\sum_{h=1}^{\hat{K}_{n}}\sum_{k<k^{\prime}}n_{hk}n_{hk^{\prime}}\left(1-\eta_{kk^{\prime}}\right).$

It is clear that

$\mathcal{B}(\bm{s},\tilde{\bm{\theta}})=\omega\sum_{k<k^{\prime}}n_{\cdot
k}n_{\cdot k^{\prime}}\left(1-\eta_{kk^{\prime}}\right)$

which is equal to zero if and only if $\eta_{kk^{\prime}}=1$ for all $k\neq
k^{\prime}$ pairs. ∎

###  C.2 Proof of Proposition 2

###### Proof.

We proceed by showing that the bounded metric properties of the Hellinger
distance are preserved by taking the posterior expectation.

  1. 1.

For each $i\leq j$, $0\leq d(g(\cdot;\theta_{i}),g(\cdot;\theta_{j}))\leq 1$
almost surely, and so $0\leq\Delta_{ij}\leq 1.$

  2. 2.

We next verify the metric axioms on the entries of $\Delta$.

    1. (a)

For each $i=1,\dots,n$, $d(g(\cdot;\theta_{i}),g(\cdot;\theta_{i}))=0$ almost
surely, and so $\Delta_{ii}=0$.

    2. (b)

$\Delta_{ij}=\Delta_{ji}$ since the Hellinger distance is symmetric and
$(\theta_{1},\dots,\theta_{n})$ are exchangeable conditional on $\bm{X}$.

    3. (c)

$\Delta_{il}\leq\Delta_{ij}+\Delta_{jl}$ since the triangle inequality in the
Hellinger distance holds almost surely for $\theta_{i},\theta_{j},\theta_{l}$
conditional on $\bm{X}$.

  3. 3.

Observe that

$\displaystyle\Delta_{ij}=\mathbb{E}_{\Pi}[\mathbb{E}_{\Pi}[d(g(\cdot;\theta_{i}),g(\cdot;\theta_{j}))\mid\lambda,\bm{X}]\mid\bm{X}]$
$\displaystyle=\int_{\Lambda}\sum_{k<k^{\prime}}\eta_{kk^{\prime}}q_{ij}^{kk^{\prime}}\>\Pi(d\lambda\mid\bm{X}),$

where

$\displaystyle
q^{kk^{\prime}}_{ij}=\Pi(\left\\{s_{i}=k,s_{j}=k^{\prime}\right\\}\cup\left\\{s_{i}=k^{\prime},s_{j}=k\right\\}\mid\lambda,\bm{X}),$
$\displaystyle\eta_{kk^{\prime}}=d(g(\cdot;\tilde{\theta}_{k}),g(\cdot;\tilde{\theta}_{k^{\prime}})),$

for all $1\leq k<k^{\prime}\leq K$. Since $\eta_{kk^{\prime}}\leq 1$,

$\displaystyle\Delta_{ij}\leq\int_{\Lambda}\sum_{k<k^{\prime}}q_{ij}^{kk^{\prime}}\Pi(d\lambda\mid\bm{X})=\Pi(s_{i}\neq
s_{j}\mid\bm{X}).$

∎

###  C.3 Proof of Proposition 3

###### Proof.

To see why there is a positive relationship between $\omega$ and $r_{\omega}$,
recall that the candidate clusterings over which we optimize
$\mathcal{R}(\widehat{\bm{c}})$ are nested. Given a value of $\omega$ and
minimizer $\bm{c}_{\omega}^{*}$, suppose we now increment the tuning parameter
from $\omega$ to $\omega+\epsilon$ for some $\epsilon>0$. Denote
$\bm{c}^{\prime}$ to be the child clustering of $\bm{c}^{*}_{\omega}$ (i.e.,
it is next in the agglomerative hierarchy), and assume that it results from
merging, say, $C^{*}_{\omega 1}$ and $C^{*}_{\omega 2}$ in
$\bm{c}_{\omega}^{*}$. Then
$\mathcal{R}(\bm{c}^{\prime})<\mathcal{R}(\bm{c}^{*}_{\omega})$ if $\sum_{i\in
C^{*}_{\omega 1},j\in C^{*}_{\omega
2}}\Delta_{ij}<(\omega+\epsilon)/(1+\omega+\epsilon)$. For large enough
$\epsilon$, this will hold, and so $\bm{c}_{\omega+\epsilon}^{*}$ will be
equal to $\bm{c}^{\prime}$ (or one of its child clusterings in the hierarchy).
Suppose then that $\epsilon$ is such that
$\bm{c}^{*}_{\omega+\epsilon}=\bm{c}^{\prime}$. Then clearly,
$r_{\omega+\epsilon}>r_{\omega}$ because any pairing $i\in C^{*}_{\omega 1}$
and $j\in C^{*}_{\omega 2}$ will now contribute a non-negative term to
$r_{\omega+\epsilon}$. Put differently, we have that
$r_{\omega+\epsilon}=r_{\omega}+(\sum_{i\in C^{*}_{\omega 1},j\in
C^{*}_{\omega 2}}\Delta_{ij})/(\sum_{i<j}\Delta_{ij})$. ∎

###  C.4 Proof of Proposition 4

###### Proof.

  1. (a)

This oracle rule follows by noting that if
$d\left\\{g(\cdot;\theta_{m}^{*}),g(\cdot;\theta_{m^{\prime}}^{*})\right\\}=1$
for all component pairs, then all components have disjoint support. Hence,
$\mathcal{R}^{*}(\widehat{\bm{c}})$ is just Binder’s loss of
$\widehat{\bm{c}}$ to $\bm{s}^{*}$, where $s_{i}^{*}=m$ if and only
$X_{i}^{0}\in\textrm{supp}(g(\cdot;\theta_{m}^{*}))$, which is uniquely
minimized at $\widehat{\bm{c}}=\bm{s}^{*}$.

  2. (b)

We can verify this by noting that, under the described conditions,
$\bm{c}_{0}$ minimizes $\mathcal{R}^{*}(\widehat{\bm{c}})$ if

$d\left\\{g(\cdot;\theta_{1}^{*}),g(\cdot;\theta_{2}^{*})\right\\}\sum_{i<j}q_{ij}^{12*}\leq\mathcal{R}^{*}(\widehat{\bm{c}}).$

A sufficient condition for this to hold is when
$d\left\\{g(\cdot;\theta_{1}^{*}),g(\cdot;\theta_{2}^{*})\right\\}<\gamma/q_{ij}^{12*}$
for all $i,j$ pairs. Since $0\leq q_{ij}^{12*}\leq 1$, this criteria can be
satisfied when
$d\left\\{g(\cdot;\theta_{1}^{*}),g(\cdot;\theta_{2}^{*})\right\\}<\gamma$.

∎

##  Appendix D Preliminaries for Asymptotic Analysis

Recall the assumptions for Theorem 1 that we make in the main article.

######  Assumption 2.

Suppose that the following conditions hold.

  * (A1) $\Lambda$ consists of measures with uniformly finite support, i.e. there exists an $L>0$ so that $|\textrm{supp}(\lambda)|\leq L$ for all $\lambda\in\Lambda$.

  * (A2) The KL minimizer $\lambda^{*}=\sum_{m=1}^{M^{*}}a_{m}^{*}\delta_{\theta_{m}^{*}}$ exists and is unique.

  * (A3) $f^{*}$ and $f^{0}$ are such that $\textrm{supp}(f^{0})\subseteq\textrm{supp}(f^{*})$.

  * (A4): $\mathcal{G}=\left\\{\tilde{g}(x-\theta):\theta\in\Theta\right\\}$ for some bounded probability density function $\tilde{g}(\cdot)$, there exists a $\zeta>0$ so that $\tilde{g}(z)$ is $\zeta$-Hölder continuous, and for any fixed $\theta^{\prime}$, $D_{g}(\theta,\theta^{\prime})=d\left\\{g(\cdot;\theta),g(\cdot;\theta^{\prime})\right\\}$ is continuous in $\theta$.

  * (A5) The mixing measure of the fitted model contracts to the oracle in the Wasserstein 2-distance, i.e. there exists a non-negative sequence $\epsilon_{n}$ so that $\epsilon_{n}\to 0$ and

$\rho_{n}(\bm{X})=\Pi\left\\{\lambda\in\Lambda:W_{2}(\lambda,\lambda^{*})\underset{\sim}{>}\epsilon_{n}\mid\bm{X}\right\\}\overset{\mathbb{P}^{0}}{\longrightarrow}0.$
(12)

In this section, we will be primarily concerned with introducing concepts
necessary for understanding the proof of Theorem 1. We will also detail
conditions under which (A5) holds, which is a statement on the posterior
contraction of the mixing measure. The main contribution that (A5) has to the
proof of Theorem 1 is that it establishes consistency of the cluster-specific
parameters $\bm{a}$ and $\tilde{\bm{\theta}}$. Conditions under which (A5)
hold have been extensively explored in Nguyen, (2013), Ho and Nguyen, 2016b ,
and Guha et al., (2021).

###  D.1 Additional Notation and Definitions

Given $x=(x_{1},\dots,x_{p})^{T}\in\mathbb{R}^{p}$, let $\left\lVert
x\right\rVert=\sqrt{\sum_{s=1}^{p}x_{s}^{2}}$. For $q\geq 1$, the $q$-th order
Wasserstein distance between the finite measures
$\lambda=\sum_{k=1}^{K}a_{k}\delta_{\theta_{k}}$ and
$\lambda^{\prime}=\sum_{k^{\prime}=1}^{K^{\prime}}a^{\prime}_{k^{\prime}}\delta_{\theta_{k^{\prime}}^{\prime}}$
is

$W_{q}(\lambda,\tilde{\lambda})=\inf_{\bm{b}\in\mathcal{C}(\bm{a},\bm{a}^{\prime})}\left(\sum_{k,k^{\prime}}b_{kk^{\prime}}\left\lVert\theta_{k}-\theta^{\prime}_{k^{\prime}}\right\rVert^{q}\right)^{1/q},$

where $\mathcal{C}(\bm{a},\bm{a}^{\prime})$ is the set of all couplings with
marginal probability vectors $\bm{a}$ and $\bm{a}^{\prime}$. That is, for all
$k=1,\dots,K$, $\sum_{k^{\prime}=1}^{K^{\prime}}b_{kk^{\prime}}=a_{k}$ and for
all $k^{\prime}=1,\dots,K^{\prime}$,
$\sum_{k=1}^{K}b_{kk^{\prime}}=a_{kk^{\prime}}$. The Wasserstein q-distance is
a metric over the space of probability measures on $\Theta$, and therefore
$W_{q}(\lambda,\lambda)=0$, if $\lambda\neq\lambda^{\prime}$ then
$W_{q}(\lambda,\lambda^{\prime})>0$,
$W_{q}(\lambda,\lambda^{\prime})=W_{q}(\lambda^{\prime},\lambda)$, and $W_{q}$
obeys the triangle inequality. We will only consider mixtures having
identifiable mixing measures. That is, if $\int
g(x;\theta)\lambda(d\theta)=\int g(x;\theta)\lambda^{\prime}(d\theta)$ for all
$x$, then $\lambda=\lambda^{\prime}$. This identifiability condition is
satisfied by Gaussian mixtures, along with several exponential family mixtures
(Barndorff-Nielsen,, 1965) and location family mixtures (Teicher,, 1961).

Recall that $\mathbb{P}^{0}$ refers to probability statements with respect to
the true data generating process, with $\mathbb{P}^{0}(A)=\int_{A}f^{0}(x)dx$.
For any sequences $y_{n}$ and $z_{n}$, $y_{n}\underset{\sim}{<}z_{n}$ if there
exists a $N\in\mathbb{N}$ and constant $B>0$ so that $|y_{n}|\leq B|z_{n}|$
for all $n\geq N$; $y_{n}\asymp z_{n}$ if $y_{n}\underset{\sim}{<}z_{n}$ and
$z_{n}\underset{\sim}{<}y_{n}$; and $y_{n}\sim z_{n}$ if $y_{n}/z_{n}\to 1$. A
random variable $Z_{n}=o_{\mathbb{P}^{0}}(1)$ if $Z_{n}\to 0$ in
$\mathbb{P}^{0}$-probability. Finally, we say that an event $A_{n}$ holds with
high probability if $\mathbb{P}^{0}(A_{n})\to 1$ as $n\to\infty$.

###  D.2 Strong Identifiability and Lipschitz Conditions

As stated in Section 3, Wasserstein contraction of $\lambda$ will generally
require that $\mathcal{G}$ is strongly identifiable, referring to a linear
independence condition on the derivatives of $g(x;\theta)$. We formally define
these concepts below. Recall that $\Theta\subset\mathbb{R}^{p}$ is assumed to
be compact and $f$ to be identifiable in $\lambda$. In these definitions,
$\nabla_{g}(\theta)$ and $\textbf{H}_{g}(\theta)$ denote the gradient and
Hessian of $g(x;\theta)$ with respect to $\theta$.

######  Definition D.1.

The family $\mathcal{G}=\left\\{g(\theta):\theta\in\Theta\right\\}$ is first
order identifiable if $g(x;\theta)$ is differentiable in $\theta$ and, given
$r$ atoms $\theta_{1},\dots,\theta_{r}\in\Theta$, if there exists
$\beta_{1},\dots,\beta_{r}\in\mathbb{R}$ and
$\rho_{1},\dots,\rho_{r}\in\mathbb{R}^{p}$ so that for all $x$

$\sum_{s=1}^{r}\left\\{\beta_{s}g(x;\theta_{s})+\rho_{s}^{T}\nabla_{g}(\theta_{s})\right\\}=0,$

then $\beta_{s}=0$ and $\rho_{s}=\mathbf{0}$ for all $s=1,\dots,r$.

######  Definition D.2.

The family $\mathcal{G}=\left\\{g(\theta):\theta\in\Theta\right\\}$ is second
order identifiable if $g(x;\theta)$ is twice differentiable in $\theta$ and,
given $r$ atoms $\theta_{1},\dots,\theta_{r}\in\Theta$, if there exists
$\beta_{1},\dots,\beta_{r}\in\mathbb{R}$,
$\rho_{1},\dots,\rho_{r}\in\mathbb{R}^{p}$, and
$\nu_{1},\dots,\nu_{r}\in\mathbb{R}^{p}$ so that for all $x$

$\sum_{s=1}^{r}\left\\{\beta_{s}g(x;\theta_{s})+\rho_{s}^{T}\nabla_{g}(\theta_{s})+\nu_{s}^{T}\textbf{H}_{g}(\theta_{s})\nu_{s}\right\\}=0,$

then $\beta_{s}=0$ and $\rho_{s}=\nu_{s}=\mathbf{0}$ for all $s=1,\dots,r$.

It is clear from the two definitions that second-order identifiability implies
first-order identifiability. These notions are commonly used in the literature
because they bridge the connection between convergence from $f$ to $f^{0}$ (or
$f^{*}$) in the Hellinger distance and convergence from $\lambda$ to
$\lambda^{0}$ (or $\lambda^{*}$) in Wasserstein distance. This will be useful
for showing concentration of the weights and atoms, as this is straightforward
to extract from convergence of the mixing measure. A canonical example of a
second-order identifable family is Gaussian location mixtures for a known
covariance matrix $\Sigma$, i.e.
$\mathcal{G}=\left\\{\mathcal{N}_{p}(\theta,\Sigma):\theta\in\Theta\right\\}$,
whereas location-scale Gaussian families
$\mathcal{G}=\left\\{\mathcal{N}_{p}(\mu,\Sigma):\mu\in\mathbb{R}^{p},\Sigma\in\mathbb{S}^{p}\right\\}$,
where $\mathbb{S}^{p}$ is the set of $p\times p$ positive semi-definite
matrices, are only first-order identifiable.

In addition to strong identifiability, various Lipschitz conditions on the
derivatives of $g(x;\theta)$ are required for contraction of the mixing
measure. The strongest condition is the second order integral Lipschitz
property, as introduced in Guha et al., (2021).

######  Definition D.3.

$\mathcal{G}$ satisfies the integral Lipschitz property of the second order
with respect to the mixing measures $\lambda^{0}$ and $\lambda^{*}$ if
$g(x;\theta)$ is twice differentiable in $\theta$ and for all $x\in\chi$,

$|\zeta^{T}(\textbf{H}_{g}(\theta_{1})-\textbf{H}_{g}(\theta_{2}))\zeta|\leq
M(x)\left\lVert\theta_{1}-\theta_{2}\right\rVert^{\xi}\left\lVert\zeta\right\rVert^{2}$
(13)

for any $\zeta\in\mathbb{R}^{p}$ and for some $\xi>0$ that does not depend on
$x$ and $\theta_{1},\theta_{2}$, where $M(x)$ is a function of $x$ so that
$\mathbb{E}_{0}[M(X)/f^{*}(X)]<\infty$.

Note that this is not only a condition on $\mathcal{G}$, but also on
$\lambda^{0}$ and $f^{0}$. The above inequality is effectively a smoothness
condition on the Hessian of $g(x;\theta)$, where the Lipschitz constant is a
function of $x$. Guha et al., (2021) invokes the integral Lipschitz property
when studying the misspecificed regime. Related to the integral Lipschitz
property is the uniform Lipschitz property, which was introduced in Ho and
Nguyen, 2016b .

######  Definition D.4.

$\mathcal{G}$ satisfies the uniform Lipschitz property of the second order if
$g(x;\theta)$ is twice differentiable in $\theta$ and for all $x\in\chi$,

$|\zeta^{T}(\textbf{H}_{g}(\theta_{1})-\textbf{H}_{g}(\theta_{2}))\zeta|\leq
M\left\lVert\theta_{1}-\theta_{2}\right\rVert^{\xi}\left\lVert\zeta\right\rVert^{2}$

for any $\zeta\in\mathbb{R}^{p}$ and for some $\xi,M>0$ that do not depend on
$x$ and $\theta_{1},\theta_{2}$.

The uniform Lipschitz property is generally invoked in well-specified regimes.
As in (13), it is also a Lipschitz condition on the Hessian of $g(x;\theta)$.
The second order Lipschitz property is satisfied by the location-scale
Gaussian family and location-scale Student’s t family (Ho et al.,, 2020). A
similar notion exists for the gradient, known as the first-order Lipschitz
property.

Lipschitz conditions on the derivatives, particularly on the Hessian, of a
probability density function are frequently assumed in the literature on
mixing measure concentration. See, for example, Lemma 2 in Chen, (1995) and
Theorem 1 in Nguyen, (2013). While it is difficult to characterize what
choices of $f^{0}$ and $\mathcal{G}$ can admit the integral Lipschitz property
in the misspecified regime, the property will hold if $\mathcal{G}$ is
uniformly Lipschitz of the second order and
$\int_{\mathbb{R}^{p}}\frac{f^{0}(x)}{f^{*}(x)}dx<\infty$. Hence, for
uniformly Lipschitz families, such as location Gaussians, the integral
Lipschitz property can be reformulated as a tail condition on $f^{*}$ and
$f^{0}$.

###  D.3 Properties of Wasserstein Concentration

The proof of Theorem 1 relies on characterizing the mixing measures in the set

$B_{\epsilon_{n}}(\lambda^{*})=\left\\{\lambda\in\Lambda:W_{2}(\lambda,\lambda^{*})\underset{\sim}{<}\epsilon_{n}\right\\},$

or the $\epsilon_{n}$ ball around $\lambda^{*}$. (A1) becomes a crucial
assumption for doing so, as it tells us that each $\lambda\in\Lambda$ has at
most $L$ support points, where $L<\infty$. In the exact fitted case where
$\Pi(K=M^{*})=1$, Nguyen, (2013) shows that for all $m\in[M^{*}]$, there
exists a unique $k\in[K]$ so that
$\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert\underset{\sim}{<}\epsilon_{n}$
and $|a_{k}-a_{m}^{*}|\underset{\sim}{<}\epsilon_{n}$. In this next
proposition, we will show that a similar result holds for the overfitted case,
i.e., when $\Pi(K>M^{*})>0$.

######  Proposition 5.

Let $\lambda=\sum_{k=1}^{K}a_{k}\delta_{\theta_{k}}\in
B_{\epsilon_{n}}(\lambda^{*})$, where $M^{*}\leq K\leq L<\infty$ and
$\epsilon_{n}\to 0$, and fix $0<\delta<1$.

  1. 1.

For large enough $n$, for all $m\in[M^{*}]$ there exists an index set
$\mathcal{I}^{(\lambda)}_{m}\subset[K]$ so that

$\max_{k\in\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert\underset{\sim}{<}\epsilon_{n},\qquad\bigg{|}\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}-a_{m}^{*}\bigg{|}\underset{\sim}{<}\max(\epsilon_{n},\epsilon_{n}^{2\delta}),$
(14)

and
$\mathcal{I}^{(\lambda)}_{m}\cap\mathcal{I}^{(\lambda)}_{m^{\prime}}=\emptyset$
for all $m\neq m^{\prime}=1,\dots,M^{*}$.

  2. 2.

Let
$\mathcal{I}^{(\lambda)}=\bigcup_{m=1}^{M^{*}}\mathcal{I}^{(\lambda)}_{m}$.
Then,

$\sum_{k\not\in\mathcal{I}^{(\lambda)}}a_{k}\underset{\sim}{<}\epsilon_{n}^{2\delta}.$
(15)

Proposition 5 gives us the necessary tools to establish convergence to the
oracle. The statement in (14) tells us that for large enough $n$, we can
allocate a subset of the components in $\lambda$ to $M^{*}$ groups, given by
$\mathcal{I}^{(\lambda)}_{m}$ for $m\in[M^{*}]$. In the $m$th group, all of
the atoms contract to $\theta_{m}^{*}$, and the total weight of the $m$th
group contracts to $a_{m}^{*}$. (15) tells us that any components left over,
i.e. any components with atoms that do not contract to an atom in
$\lambda^{*}$, will have total weight diminishing at an
$\epsilon_{n}^{2\delta}$ rate. This behavior is comparable to Theorem 1 in
Rousseau and Mengersen, (2011), which established that under certain settings
of overfitted finite mixture models, the total weight of the overfitted
components empty at an approximately $n^{-1/2}$ rate. Under the settings of
Theorem 4.3 from Guha et al., (2021), (15) implies that the redundant
components empty at an approximately $\left(\log n/n\right)^{1/2}$ rate for
misspecified models. We now present the proof of the proposition.

###  D.4 Proof of Proposition 5

###### Proof.

Let $\lambda\in B_{\epsilon}(\lambda^{*})$, where $K\geq M^{*}$. We now show
concentration of the weights and atoms using techniques from Nguyen, (2011)
and Nguyen, (2013). Let $0<\delta<1$ be a small positive constant. For large
enough $n$, $a_{\min}^{*}\underset{\sim}{>}\epsilon_{n}^{2\delta}$, implying

$W_{2}^{2(1-\delta)}(\lambda,\lambda^{*})\geq\min_{k}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}.$

For any $m\in[M^{*}]$ let

$\mathcal{I}^{(\lambda)}_{m}=\left\\{k\in[K]:\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\leq
W_{2}^{2(1-\delta)}(\lambda,\lambda^{*})\right\\}$

and
$\mathcal{I}^{(\lambda)}=\bigcup_{m=1}^{M^{*}}\mathcal{I}^{(\lambda)}_{m}$. It
follows that

$\displaystyle
W_{2}^{2}(\lambda,\lambda^{*})\geq\sum_{k\not\in\mathcal{I}^{(\lambda)}}a_{k}\min_{m}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}>W_{2}^{2(1-\delta)}(\lambda,\lambda^{*})\sum_{k\not\in\mathcal{I}^{(\lambda)}}a_{k}$
$\displaystyle\implies\sum_{k\not\in\mathcal{I}^{(\lambda)}}a_{k}<W_{2}^{2\delta}(\lambda,\lambda^{*})\underset{\sim}{<}\epsilon_{n}^{2\delta}.$
(16)

Now, we will show concentration of the weights towards their oracle values.
For fixed $m$,

$\displaystyle\epsilon_{n}^{2}\underset{\sim}{>}W_{2}(\lambda,\lambda^{*})\geq\sum_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}_{m}^{(\lambda)}}b_{km}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\geq\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\sum_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}_{m}^{(\lambda)}}b_{km}$
$\displaystyle=\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\left(a_{m}^{*}-\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}b_{km}-\sum_{k\not\in\mathcal{I}^{(\lambda)}}b_{km}\right)$
$\displaystyle=\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{h}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\left(a_{m}^{*}-\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}+\sum_{m\in\mathcal{I}^{(\lambda)}_{h}}\sum_{m^{\prime}\neq
m}b_{km^{\prime}}-\sum_{k\not\in\mathcal{I}^{(\lambda)}}b_{km}\right)$
$\displaystyle\geq\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\left(a_{m}^{*}-\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}-\sum_{k\not\in\mathcal{I}^{(\lambda)}}b_{km}\right)$
$\displaystyle\geq\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\left(a_{m}^{*}-\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}-\sum_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}b_{km}-\sum_{k\not\in\mathcal{I}^{(\lambda)}}b_{km}\right)$
$\displaystyle=\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\left(\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}b_{km}-\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}\right)$
$\displaystyle=-\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}\sum_{m^{\prime}\neq
m}b_{km^{\prime}}.$

For any $m^{\prime}\neq m$,
$k^{\prime}\in\mathcal{I}^{(\lambda)}_{m^{\prime}}$, and
$k\in\mathcal{I}^{(\lambda)}_{m}$,

$\displaystyle\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert\leq\left\lVert\theta_{k^{\prime}}-\theta_{m}^{*}\right\rVert$
$\displaystyle\leq\left\lVert\theta_{k^{\prime}}-\theta_{m^{\prime}}^{*}\right\rVert+\left\lVert\theta_{k}-\theta_{m^{\prime}}^{*}\right\rVert+\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert\underset{\sim}{<}\epsilon_{n}+\left\lVert\theta_{k}-\theta_{m^{\prime}}^{*}\right\rVert.$

Since $\Theta$ is compact, this shows that

$\displaystyle-\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}\sum_{m^{\prime}\neq
m}b_{km^{\prime}}$
$\displaystyle\underset{\sim}{>}-\epsilon_{n}-\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}\sum_{m^{\prime}\neq
m}b_{km^{\prime}}\left\lVert\theta_{k}-\theta_{m^{\prime}}^{*}\right\rVert^{2}\underset{\sim}{>}-\epsilon_{n}-W_{2}(\lambda,\lambda^{*})^{2}.$

Therefore,

$\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\bigg{|}a_{m}^{*}-\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}-\sum_{k\not\in\mathcal{I}^{(\lambda)}}b_{km}\bigg{|}\underset{\sim}{<}\epsilon_{n}.$

Observe that there exists $m^{\prime}\neq m$ so that

$\min_{k\in\mathcal{I}^{(\lambda)}\setminus\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\underset{\sim}{>}\left\lVert\theta_{m}^{*}-\theta_{m^{\prime}}^{*}\right\rVert^{2}-\epsilon_{n}.$

The above derivation implies that there exists constants $M,M^{\prime},$ and
$M^{\prime\prime}>0$ so that

$\displaystyle\bigg{|}a_{m}^{*}-\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}-\sum_{k\not\in\mathcal{I}^{(\lambda)}}b_{km}\bigg{|}\leq\frac{M\epsilon_{n}}{\left\lVert\theta_{m}^{*}-\theta_{m^{\prime}}^{*}\right\rVert^{2}-M^{\prime}\epsilon_{n}}$
$\displaystyle\implies\bigg{|}a_{m}^{*}-\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}\bigg{|}\leq\nu_{n}:=\frac{M\epsilon_{n}}{\left\lVert\theta_{m}^{*}-\theta_{m^{\prime}}^{*}\right\rVert^{2}-M^{\prime}\epsilon_{n}}+(K-M^{*})M^{\prime\prime}\epsilon_{n}^{2\delta}\underset{\sim}{<}\max(\epsilon_{n},\epsilon_{n}^{2\delta}).$

∎

##  Appendix E Proof of Theorem 1

###### Proof.

Recall that $\lambda^{*}=\sum_{m=1}^{M^{*}}a_{m}^{*}\delta_{\theta_{m}^{*}}$
and that by (A1), if $\lambda\in\Lambda$ then
$\lambda=\sum_{k=1}^{K}a_{k}\delta_{\theta_{k}}$ for some $K\leq L$. Let
$B_{\epsilon_{n}}(\lambda^{*})=\left\\{\lambda\in\Lambda:W_{2}(\lambda,\lambda^{*})\underset{\sim}{<}\epsilon_{n}\right\\}$
and $\Lambda_{M^{*},L}=\left\\{\lambda\in\Lambda:|\textrm{supp}(\lambda)|\geq
M^{*}\right\\}$. We decompose the posterior expectation of $\mathcal{D}_{ij}$
into

$\displaystyle\Delta_{ij}=\mathbb{E}[\mathcal{D}_{ij}\mid\bm{X}]$
$\displaystyle=\int_{\Lambda_{M^{*},L}\cap
B_{\epsilon_{n}}(\lambda^{*})}\mathcal{D}_{ij}\Pi(d\phi_{ij}\mid\bm{X})+\int_{\Lambda_{M^{*},L}^{c}\cap
B_{\epsilon_{n}}(\lambda^{*})}\mathcal{D}_{ij}\Pi(d\phi_{ij}\mid\bm{X})+\int_{B_{\epsilon_{n}}^{c}(\lambda^{*})}\mathcal{D}_{ij}\Pi(d\phi_{ij}\mid\bm{X}),$

where $\phi_{ij}=(\theta_{i},\theta_{j},\lambda)$. This implies

$\int_{\Lambda_{M^{*},L}\cap
B_{\epsilon_{n}}(\lambda^{*})}\mathcal{D}_{ij}\Pi(d\phi_{ij}\mid\bm{X})\leq\Delta_{ij}\leq\int_{\Lambda_{M^{*},L}\cap
B_{\epsilon_{n}}(\lambda^{*})}\mathcal{D}_{ij}\Pi(d\phi_{ij}\mid\bm{X})+\nu_{n}(\bm{X})+\rho_{n}(\bm{X}),$

where
$\nu_{n}(\bm{X})=\Pi(|\textrm{supp}(\lambda)|<M^{*},W_{2}(\lambda,\lambda^{*})\underset{\sim}{<}\epsilon_{n}\mid\bm{X})$
and
$\rho_{n}(\bm{X})=\Pi(W_{2}(\lambda,\lambda^{*})\underset{\sim}{>}\epsilon_{n}\mid\bm{X})$.
By assumption (A5),
$\rho_{n}(\bm{X})\overset{\mathbb{P}^{0}}{\longrightarrow}0$. In addition,
$\nu_{n}(\bm{X})\overset{\mathbb{P}^{0}}{\longrightarrow}0$, and we can show
that $\mathbb{P}^{0}\left\\{\nu_{n}(\bm{X})=0\right\\}=1$ for large enough
$n$. To see this, note that $\lambda\in W_{2}(\lambda,\lambda^{*})$ implies
that for all $\theta_{m}^{*}\in\textrm{supp}(\lambda)$, there exists a
$\theta_{k}\in\textrm{supp}(\lambda)$ so that
$\left\lVert\theta_{m}^{*}-\theta_{k}\right\rVert\underset{\sim}{<}\epsilon_{n}$.
More formally,

$\displaystyle\epsilon_{n}^{2}\underset{\sim}{>}W_{2}(\lambda,\lambda^{*})=\sum_{k,m}b_{km}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}\geq\sum_{k}b_{km}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}$
$\displaystyle\geq\min_{k}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}a_{m}^{*}\geq\min_{k}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert^{2}a_{\min}^{*},$

where $a_{\min}^{*}=\min_{m}a_{m}^{*}$. If we set $n$ large enough so that
$\min_{m\neq
m^{\prime}}\left\lVert\theta_{m}^{*}-\theta_{m^{\prime}}^{*}\right\rVert\underset{\sim}{>}2\epsilon_{n}$,
then we have that there must be at least $M^{*}$ support points in $\lambda$.
Hence, for large enough $n$, $\Lambda_{M^{*},L}^{c}\cap
B_{\epsilon_{n}}(\lambda^{*})=\emptyset$. Therefore, we will omit
$\nu_{n}(\bm{X})$ from the remainder of the proof by assuming that $n$ is
large enough so that $\min_{m\neq
m^{\prime}}\left\lVert\theta_{m}^{*}-\theta_{m^{\prime}}^{*}\right\rVert\underset{\sim}{>}2\epsilon_{n}$.
We then formulate the bounds of $\Delta_{ij}$ via

$\int_{\Lambda_{M^{*},L}\cap
B_{\epsilon_{n}}(\lambda^{*})}\mathcal{D}_{ij}\Pi(d\phi_{ij}\mid\bm{X})\leq\Delta_{ij}\leq\int_{\Lambda_{M^{*},L}\cap
B_{\epsilon_{n}}(\lambda^{*})}\mathcal{D}_{ij}\Pi(d\phi_{ij}\mid\bm{X})+\rho_{n}(\bm{X}).$
(17)

Let $0<\delta<1$ and $\lambda\in\Lambda_{M^{*},L}\cap
B_{\epsilon_{n}}(\lambda^{*})$. Proposition 5 states that for any
$m\in[M^{*}]$ there exists a set of indices $\mathcal{I}^{(\lambda)}_{m}$ so
that
$\max_{k\in\mathcal{I}^{(\lambda)}_{m}}\left\lVert\theta_{k}-\theta_{m}^{*}\right\rVert\underset{\sim}{<}\epsilon_{n}$
and
$\bigg{|}\sum_{k\in\mathcal{I}^{(\lambda)}_{m}}a_{k}-a_{m}^{*}\bigg{|}\underset{\sim}{<}\max(\epsilon_{n},\epsilon_{n}^{2\delta})$,
and in addition that
$\sum_{k\not\in\mathcal{I}^{(\lambda)}}a_{k}\underset{\sim}{<}\epsilon_{n}^{2\delta}$,
where
$\mathcal{I}^{(\lambda)}=\bigcup_{m=1}^{M^{*}}\mathcal{I}^{(\lambda)}_{m}$. We
will now translate contraction of the atoms and weights into contraction of
the $\Delta_{ij}$ terms. Recall that we can write

$\int_{\Lambda_{M^{*},L}\cap
B_{\epsilon_{n}}(\lambda^{*})}\mathcal{D}_{ij}\Pi(d\phi_{ij}\mid\bm{X})=\int_{\Lambda_{M^{*},L}\cap
B_{\epsilon_{n}}(\lambda^{*})}\sum_{k<k^{\prime}}\eta_{kk^{\prime}}\>q_{ij}^{kk^{\prime}}\Pi(d\lambda\mid\bm{X}),$
(18)

We then decompose the integrand in (18):

$\displaystyle\sum_{k<k^{\prime}}\eta_{kk^{\prime}}q_{ij}^{kk^{\prime}}$
$\displaystyle=\sum_{k<k^{\prime}\in\mathcal{I}^{(\lambda)}}\eta_{kk^{\prime}}q_{ij}^{kk^{\prime}}+\sum_{(k\not\in\mathcal{I}^{(\lambda)})\textrm{
or
}(k^{\prime}\not\in\mathcal{I}^{(\lambda)})}\eta_{kk^{\prime}}q_{ij}^{kk^{\prime}}$
$\displaystyle=\sum_{m=1}^{M^{*}}\sum_{k<k^{\prime}\in\mathcal{I}^{(\lambda)}_{m}}\eta_{kk^{\prime}}q_{ij}^{kk^{\prime}}+\sum_{m<m^{\prime}}\sum_{\begin{subarray}{c}k\in\mathcal{I}^{(\lambda)}_{m}\\\
k^{\prime}\in\mathcal{I}^{(\lambda)}_{m^{\prime}}\end{subarray}}\eta_{kk^{\prime}}q_{ij}^{kk^{\prime}}+\sum_{(k\not\in\mathcal{I}^{(\lambda)})\textrm{
or
}(k^{\prime}\not\in\mathcal{I}^{(\lambda)})}\eta_{kk^{\prime}}q_{ij}^{kk^{\prime}}.$
(19)

Let $\psi_{m}^{*}(\theta)=D(\theta,\theta_{m}^{*})$. By (A4), Lemma 1, and
Remark 1, for any $1\leq m\leq M^{*}$,

$\displaystyle\max_{k\in\mathcal{I}^{(\lambda)}_{m}}D(\theta_{k},\theta_{m}^{*})\leq\kappa_{nm}=\max_{\theta:\left\lVert\theta-\theta_{h}^{*}\right\rVert\underset{\sim}{<}\epsilon_{n}}\varphi^{*}_{m}(\theta)$

and $\kappa_{nm}\to 0$. Let $\kappa_{n}:=\max_{m=1,\dots,M^{*}}\kappa_{nm}$.
Then, for any $k,k^{\prime}\in\mathcal{I}_{h}^{(\lambda)}$,