Interpretable Clustering with the Distinguishability Criterion
Ali Turfah Xiaoquan Wen



Fast and Robust Expectation Propagation MIMO Detection via Preconditioned Conjugated Gradient
    Luca Schmid1, Dominik Sulz2, and Laurent Schmalen1
1Communications Engineering Lab, Karlsruhe Institute of Technology, 76187 Karlsruhe, Germany

2Mathematical Institute, University of Tübingen, 72076 Tübingen, Germany

Email: 1, 2

    April 29, 2024
===========================================================================================================================================================================================================================================





    Cluster analysis is a popular unsupervised learning tool used in many disciplines to identify heterogeneous sub-populations within a sample. 
    However, validating cluster analysis results and determining the number of clusters in a data set remains an outstanding problem. 
    In this work, we present a global criterion called the Distinguishability criterion to quantify the separability of identified clusters and validate inferred cluster configurations. 
    Our computational implementation of the Distinguishability criterion corresponds to the Bayes risk of a randomized classifier under the 0-1 loss. 
    We propose a combined loss function-based computational framework that integrates the Distinguishability criterion with many commonly used clustering procedures, such as hierarchical clustering, k-means, and finite mixture models. 
    We present these new algorithms as well as the results from comprehensive data analysis based on simulation studies and real data applications.








§ INTRODUCTION


Cluster analysis is a ubiquitous unsupervised learning approach to uncover latent structures and patterns in observed data. 
Clustering algorithms have been used in a wide variety of scientific applications, such as animal behavior studies <cit.>, weather anomaly detection <cit.>, disease diagnosis <cit.>, and novel cell type identification <cit.>. 
Often, the identified clusters are interpreted to represent distinct populations from which the corresponding samples originate.

Many challenges with cluster analysis, such as determining the number of clusters, arise from an inability to rigorously quantify desired cluster characteristics. 
While the precise definition of a “meaningful” cluster is usually context-dependent, it is generally accepted that the clusters should display “internal cohesion” (i.e., objects within a cluster are similar to one another) and “external isolation” (i.e., the clusters are well-separated) <cit.>. 
Despite this almost universally agreed-upon principle, quantifying the separability of the clusters, i.e., the level of external isolation with respect to internal cohesion, remains an open problem in cluster analysis.

In this paper, we introduce the Distinguishability criterion to measure the separability of a set of assumed clusters. 
The criterion is motivated by this simple intuition: if all clusters are well separated from each other, then the originating clusters for all data points (whether observed or not) should be easily traceable.
To implement the Distinguishability criterion, we formulate labeling the generating cluster for an arbitrary data point as a probabilistic classification problem.
Naturally, the difficulty (or lack thereof) of this classification problem can be described by an overall misclassification probability averaged over all possible data points.

We employ a statistical viewpoint to define the Distinguishability criterion for cluster analysis.
The partitioned observed data are taken to be realizations from cluster-specific data generative distributions, which are essential for computing the proposed misclassification probability.
Although not all clustering algorithms make explicit distributional assumptions for the presumed clusters, many commonly applied heuristics-based algorithms achieve optimal performance under specific probabilistic generative models 
<cit.>.
Furthermore, the identified cluster structures from cluster analysis are typically expected to be replicated in future datasets—an implicit assumption for consistent data generative distributions. 
As a result, statistical inference procedures based on explicit distributional assumptions have become more popular for their ability to not only enhance clustering performance but also to examine and interpret cluster structures identified by both model and heuristics-based clustering methods <cit.>.

The remainder of this paper is organized as follows. 
We first provide the mathematical definition and properties of the Distinguishability criterion. 
We then discuss its usage with existing clustering algorithms. 
Finally, we illustrate the applications of the Distinguishability criterion using both synthetic and real data from various scientific applications. 
Our implementation of the Distinguishability criterion and the analyses presented in this paper can be found at .




§ RESULTS





 §.§ The Distinguishability Criterion
 

The proposed Distinguishability criterion measures the overall separability of a given cluster configuration and is derived by quantifying the misclassification probability from a multi-class classification problem.

Given a cluster configuration where each of the K clusters corresponds to a distinct class, we denote the class label for an observation  by θ∈{ 1,2,..., K}. We assume a pre-defined classifier, δ(): ℝ^p →{1,...,K}, and evaluate the classification performance using the 0-1 loss function, i.e.,


    L(δ(), θ) = {δ() θ}.
  

The overall misclassification probability under the assumed cluster configuration, denoted by , is defined as the Bayes risk of the classifier δ(), i.e.,


    = E_[  E_θ(   L(δ(), θ) | )  ]


Using the 0-1 loss ensures that the resulting Bayes risk is a valid probability measurement, ranging from 0 to 1. 
It is naturally interpreted as the probability of misclassifying a data point under the given cluster configuration, marginalizing all potential  values and their respective true generating clusters. 
For instance, a  value close to 0 signifies a high degree of cluster separation, indicated by a minimal probability of erroneously assigning an arbitrary data point to an incorrect generating cluster (Supplementary Figure 2).

As we are primarily interested in assessing different cluster configurations, the selection of the required classifier is flexible. 
However, the choice of classifier can impact the computational efficiency of the  evaluation.
Our implementation focuses on the set of classifiers working directly with the probabilities


    π_k () := (θ = k |),   k = 1,...,K.


This set includes the optimal classifier under the 0-1 loss, δ_o, i.e.,


    δ_o() = max_k π_k()


Our default classifier for computing  is a randomized decision function, δ_r, which assigns a label to an observation  by sampling from a categorical distribution based on the probability distribution () = (π_1 (),...,π_K()), i.e.,


    δ_r() ∼ Categorical  (()).


In addition to yielding highly comparable  values to the optimal classifier within the decision-critical ranges of cluster separation (Supplementary Figure 1), the randomized classifier's computational properties enable highly efficient cluster analysis procedures.


The π_k's are the key quantities bridging the observed clustering data and . They are calculated using Bayes rule,


    π_k ()  ∝ α_k(_c)   p(|θ(_c) = k).
 

The notation emphasizes that both the prior, α_k(_c), and the likelihood function, p(|θ(_c)), are directly informed by and estimated from the clustering data.
More specifically, the prior quantifies the relative frequency of the observations arising from each assumed cluster, while the likelihood function encodes the characteristics of the corresponding cluster population, such as its centroid and spread information.
Computing π_k values is straightforward for model-based clustering algorithms such as Gaussian mixture models (GMMs) <cit.>.
For non-model-based clustering algorithms, explicit distributional assumptions specifying the parametric family of likelihood functions are required.
We illustrate these procedures for k-means and hierarchical clustering algorithms in subsequent sections. 

In summary,  is a probability measurement of global separability across inferred clusters. 
It can accommodate a wide range of distributional assumptions, making it compatible with a diverse set of clustering procedures and data modalities. 
Moreover, as a function of clustering data, _c, the estimate of  itself is a valid loss function suitable for selecting optimal cluster configurations in cluster analysis. 



 §.§ Combined Loss Function for Cluster Analysis
 

In cluster analysis, the desired cluster characteristics are often defined by multiple criteria <cit.>. A single criterion on its own, including the proposed Distinguishability criterion, is insufficient to define a practically optimal clustering solution. 
Alternatively, combining multiple loss functions targeting different desired cluster properties can result in more balanced and holistic clustering solutions. 
This observation leads to a principled way to incorporate  with other established clustering criteria and algorithms. 

Specifically, let L_1 denote a loss function associated with existing clustering algorithms. Formally, we consider a compound loss, L, as a weighted linear combination of L_1 and , i.e.,


    L = L_1 + λ,  λ >  0.
 

Because of the scale of , it is often convenient to solve the following equivalent constrained optimization problem,


    ≤τ,


where τ is a pre-specified probability threshold. It is worth noting that the stringency of the τ value may depend on the dimensionality of the clustering data.

The choice of clustering algorithm determines the functional form of L_1. For example, the distortion function or Ward's linkage are natural choices for L_1 when using k-means and hierarchical clustering methods, respectively. 
Alternatively, the negative of the gap statistic <cit.> can also be an excellent choice in these application scenarios.
For model-based clustering algorithms, the L_1 function can be derived from various model selection criteria, e.g., the negative of Bayesian information criterion (BIC). 



 §.§ Connections to Related Approaches


The misclassification probability defined by  falls into the category of internal clustering validity indices <cit.>, which assess the quality of a clustering solution without additional external information beyond the observed data.
This category includes many commonly applied statistical measures, such as the Silhouette index <cit.>, Calinski-Harabaze index <cit.>, Dunn index <cit.>, among others. 
A common behavior of internal clustering validity indices is that they evaluate both the cohesion (or compactness) within a cluster as well as the separation between clusters. 
For , the within-cluster cohesion is quantified through the estimated parameters in the likelihood function, p(|θ(_c) = k).
The separation, relative to the cohesion, is quantified by the overall misclassification probability.

Henning <cit.> and Melnykov <cit.> also compute misclassification probabilities to assess the separation between clusters in the context of mixture models for clustering.
They introduce the metrics—named “directly estimated misclassification probability" (DEMP) and DEMP+—specifically designed to compute misclassification probabilities between pairs of clusters to inform decisions about the local merging of mixture components.
In comparison,  is a global measure of the misclassification probability across all clusters. It, too, can be used to combine mixture components to form interpretable clusters, as illustrated in Section <ref>.
Additionally, the entropy criterion proposed by Celeux and Soromenho <cit.> provides another alternative approach to quantify the separability of clusters using a classification problem set-up.

The Distinguishability criterion also agrees with the principle of stability measures commonly employed in clustering analysis. Specifically, if the underlying cluster distributions are all well-separated, as indicated by low  values, data sampled repeatedly from these distributions are expected to produce consistent clustering outcomes <cit.>.
Empirical evidence has shown that  and various measures of clustering instability are highly correlated, which is demonstrated in the subsequent sections.



 §.§ Finite Mixture Models Incorporating 
 

Finite mixture models (MM) are probabilistic models that can seamlessly incorporate the Distinguishability criterion.
MM-based clustering algorithms primarily infer the distributional characteristics underlying each latent cluster. 
As a result, no additional assumptions are needed to compute  in the MM setting — the required quantities, {α_k,   p( |θ=k),  π_k}, are all direct outputs or by-products from standard MM inference procedures <cit.>, e.g., the EM algorithm.

We make an important distinction between a mixture component and an interpretable cluster, a point previously discussed by <cit.>.
We view mixture models as flexible density estimation devices, where the number of mixture components is chosen to adequately fit the observed data. On the other hand, the distribution of an underlying cluster—characterized by the Distinguishability criterion—may itself be a mixture distribution comprising multiple components. This distinction naturally leads to combining the loss functions represented by - BIC, which evaluates the goodness-of-fit of a mixture density, and , which characterizes the separation between potential clusters.

To optimize the combined loss function, we first find P(), the optimal mixture distribution with κ components, by maximizing the BIC. Subsequently, we merge the mixture components into clusters until  falls below a pre-specified threshold τ. 
Since merging mixture components into clusters does not alter the mixture component distributions in any way, the BIC is unchanged by the merging process. 
It can be shown that merging existing clusters in this manner always decreases  (Appendix <ref>). 
Specifically, under the default randomized classifier δ_r, the reduction of  by combining clusters i and j is given by


    Δ^(i, j) = 2 ∫π_i() π_j()   P(d ).
 

For clusters with little to no overlap, i.e., π_i() π_j() → 0 for all  values, merging (i,j) results in minimal changes in . On the other hand, for clusters with significant overlap, Δ^(i,j) can be substantial. 


By further utilizing the cluster merging property of  (Proposition <ref>, Methods Section), i.e.,  


    = ∑_i<jΔ^(i, j),


we propose an efficient  Hierarchical Merging (PHM) algorithm to sequentially amalgamate mixture components into clusters (Algorithm <ref>). 
Briefly, starting by assigning each of the κ mixture components to individual clusters, the PHM algorithm pre-computes Δ^(i, j) for all (i, j) cluster pairs. 
It then sequentially combines the pairs of clusters with the largest Δ^(i, j) into a single cluster and updates the Δ values for the remaining clusters. 
The process is repeated until the updated  falls at or below a pre-defined τ value. 
Intuitively, this procedure prioritizes merging the most similar or closely related clusters at each step, quantified by their Δ^(i, j) value.



By setting τ = 0, the algorithm runs until all mixture components have been merged into a single cluster. 
The complete merging process can be visualized using a dendrogram (Appendix <ref>), characterizing the hierarchical merging orders between individual mixture components and merged clusters. 
In many scientific applications, e.g., genetics and single-cell data analysis, such a dendrogram can provide a snapshot of the underlying continuous differentiation process that forms the identified clusters. We provide two examples in our real data applications (Section <ref>).


To illustrate the PHM algorithm with Gaussian mixture models (GMMs), we use the synthetic data from Section 4.1 in Baudry et al. <cit.>.
Specifically, 600 observations are drawn from six Gaussian distributions arranged along the corners of a square in the following manner. 
Two overlapping Gaussian distributions are placed in the top left corner of the square, each with 1/5 of the samples. 
The bottom left and top right corners each have a single Gaussian distribution, each contributing 1/5 of the samples. 
Finally, two overlapping Gaussian distributions are placed in the bottom right corner, each with 1/10 of the samples.


A GMM with six components is selected by BIC using the  package  <cit.>. 
The observed clustering data is shown in the left panel of Figure <ref>, with colors corresponding to an observation's assignment to one of the κ GMM components. 
The initial cluster configuration labeling each mixture component as a single cluster has = 0.139.
The heatmap in the right panel of Figure <ref> visualizes the Δ^(i, j) contributions from each pair of mixture components, showing that the main sources of  come from the overlapping components in the top left and bottom right corners.



With a threshold of τ = 0.01, the PHM algorithm sequentially combines the mixture components in the top left and bottom right corners, reducing the values of  to 0.049 and 0.004, respectively. 
In the end, the algorithm returns four clusters: one corresponding to the components in each of the four corners.



 §.§ Applications in k-means Clustering
 

Applying the Distinguishability criterion to heuristics-based clustering algorithms requires additional distributional assumptions to compute .  
Although popular algorithms of this kind—such as k-means and hierarchical clustering—rely on intuitive heuristics, their implicit connections to probability models have been well studied.
These results provide insights into the data types for which certain non-model-based clustering algorithms are expected to be optimal. 
The k-means algorithm, in particular, has been shown as equivalent to approximately optimizing a multivariate Gaussian classification likelihood function <cit.>.
Hence, when applying the Distinguishability criterion with the usual k-means distortion function, it seems natural to assume that data within each inferred partition are normally distributed.
It then becomes straightforward to estimate the necessary parameters and compute  given the partitioned data from the k-means output.


We illustrate a procedure to determine the number of clusters (K) in k-means clustering by optimizing the combined loss of  and the gap statistic <cit.>.
For a given data partitioning from the k-means algorithm, the gap statistic compares the observed within-cluster dispersion to the expected dispersion under a null reference distribution.
It subsequently estimates the optimal number of clusters corresponding to the largest gap statistic value from a range of potential K values. 
We consider observations drawn from three Gaussian clusters in ℝ^2 centered at (0, 0), (1.75, 1.75), and (-4, 4) with identity covariance matrices (n_k = 150). The observed data are shown in Figure <ref>, with substantial overlap between the data generated from two of the three distributions.
The k-means algorithm is performed for K = 1 to 7 using the R package  <cit.> with  initialization <cit.>. In addition to the  values, we also present the averaged Silhouette index <cit.> for each value of K=2, …, 7.

The values of the gap statistic, , and the Silhouette index for different K are shown in Figure <ref>. When used alone, the gap statistic selects K = 3, coinciding with the number of generating distributions.
There is a noticeable difference in the inferred clusters' separability between K=2 (= 0.002) and K=3 (= 0.062).
For a  threshold τ = 0.01, optimizing the combined loss function (<ref>) leads to selecting K = 2, which seems most reasonable judging by the data visualization.
This decision is also consistent with the Silhouette index, which takes a maximal value of 0.632 at K= 2 compared to 0.522 at K = 3.
Additionally we find that  is highly negatively correlated with cluster stability measures. We present the values of the stability measure proposed by Lange et al. <cit.> using the adjusted Rand index to compute cluster similarities as well as the prediction strength <cit.> in Supplementary Table 1.


As noticed in the original paper <cit.>, the gap statistic can struggle to determine K when the underlying generative distributions have substantial overlapping support.   
The above numerical example illustrates that the Distinguishability criterion can alleviate this challenge in k-means clustering. 





 §.§ Hypothesis Testing in Hierarchical Clustering


It has become increasingly common to perform formal statistical testing in post-clustering analysis to reduce cluster over-identification. 
Recent works by Gao et al. <cit.>, Chen et al. <cit.>, and Grabski et al. <cit.> highlight the necessity and importance of such analysis in scientific applications, where false positive findings of clusters are considered more costly than false negative findings. 
In many application contexts, it is often of interest to assess whether the partitions of the observed data output by heuristics-based clustering algorithms (e.g., hierarchical clustering) could arise from a single homogeneous distribution (most commonly, a single Gaussian distribution) by chance.  
 emerges as a natural test statistic in this parametric hypothesis testing framework because its values are expected to be quantitatively different under the null and alternative scenarios.


Formally, we consider testing the null hypothesis,


    H_0:


in the hierarchical clustering setting <cit.>.
Following hierarchical clustering of observed data, we compute  for the pair of inferred clusters resulting from the first split of the dendrogram (i.e., K=2). 
The null distribution of  can be estimated by a simple Monte Carlo procedure that repeatedly samples data from H_0, performs hierarchical clustering, and computes  for the partitions defined by the first split.
Alternatively, the p-value of the  statistic can be derived from the standard parametric bootstrap procedure <cit.>. 

We present the results of simulation studies examining the performance of  in this hypothesis testing setting.



First, we illustrate 's ability to control for spurious cluster detection. For 5,000 simulation replicates, we draw 150 observations from a N(0, 1) distribution and perform hierarchical clustering based on the squared Euclidean distance with the Ward linkage. 
We compute  and derive the corresponding p-values based on the estimated null distribution from 5,000 additional Monte Carlo simulations, each with a sample size of 150. 
We also compute the  p-values using a bootstrap procedure, using 500 bootstrap replicates to estimate the p-value for each simulation replicate. 
For comparison, we also compute the p-values from the selective inference procedure proposed by Gao et al. <cit.> as well as the standard two-sample t-test.
The results are shown in the top panel of Figure <ref>. 
As expected, the p-values from the two-sample t-test are all in the range of 10^-20 to 10^-40, providing no control for false positive findings.
This is because the naive t-test fails to take into account that the hierarchical clustering procedure always partitions data according to their observed values even under H_0, as noted by <cit.>.
In contrast, p-values derived from  and Gao et al.'s methods are both roughly uniformly distributed under H_0, suggesting well-controlled type I error rates.
Specifically, we find that = 0.094 corresponds to the cutoff at the 5% α level, and the realized type I error rate is 4.8%.
The type I error rate from the bootstrap procedure using  at the same control level is 1%, which is more conservative. 


Second, we examine the power of the -based hypothesis test to identify truly separated clusters. 
Specifically, we draw samples from two distinct Gaussian distributions, N(μ_1, σ^2) and N(μ_2, σ^2), where a range of |μ_1 - μ_2|/σ values are selected from the set {1, 1.5, …, 7}. 
In each alternative scenario represented by a unique |μ_1 - μ_2|/σ value, we draw 75 observations from each cluster distribution and compute the p-values using  as well as Gao et al.'s procedure. 
We calculate the power based on 500 replicates for each alternative scenario. 
As can be seen in the bottom panel of Figure <ref>,  exhibits higher power than Gao et al.'s procedure for moderate-to-large degrees of cluster separation (defined by |μ_1 - μ_2|/σ∈ [2, 6]).



 §.§ Real Data Applications
 



  §.§.§ Cluster Analysis of Palmer Penguin Data


We analyze the  data from the  package <cit.>, which consists of bill, flipper, and mass measurements from three species of penguins in the Palmer Archipelago between 2007 and 2009. 
Following the analysis by <cit.>, we consider the subset of female penguins with complete data for bill and flipper length (both measured in millimeters), leaving us with 165 observations. Prior to performing the clustering, we center and scale the observations so that the measurements have zero mean and unit variance.
Following the procedure laid out in Section <ref>, we perform k-means clustering of the scaled observations into K = 1, …, 8 clusters and compute the gap statistic as well as  for each grouping of the observations. 
Figure <ref> shows the data used for clustering and the values of the gap statistic and  for different values of K.

For the  threshold of 0.05, the combined loss defined by the gap statistic and  is optimized at K=3.  
Both the visualization of the clustering data and the external species information suggest that the result is reasonable. We observe that the values of  in this example are strongly negatively correlated with other cluster validity indices <cit.> to evaluate the k-means clustering partitions (Supplementary Table 2).
We repeat this analysis using hierarchical clustering based on the squared Euclidean distance with the Ward linkage and come to the same conclusions (Supplementary Figure 3 and Supplementary Table 3).






  §.§.§ Inferring Population Structure from HGDP data
 

In this illustration, we apply the PHM algorithm to perform cluster analysis on the genetic data from the Human Genome Diversity Project (HGDP) <cit.>, aiming to identify population structures. 
The dataset comprises 927 unrelated individuals sampled worldwide and genotyped at 2,543 autosomal SNPs. 
The geographic sampling locations are broadly divided into 7 continental groups: Europe, Central/South Asia (C/S Asia), Africa, Middle East, the Americas, East Asia, and Oceania.

Following the standard procedures for genetic data analysis, we pre-process the genotype matrix using principal component analysis (PCA) and select the first 5 PCs for our analysis based on the elbow point of the scree plot (Supplementary Figure 4). 
To apply the PHM algorithm, we fit a GMM to the dimension-reduced PC score matrix and select the model with 9 components based on the BIC. 
Each individual's posterior component assignment probability π_k is shown in the distruct plot <cit.> in Figure <ref>. 
Except for the Europe, Central/South Asia, and Middle East groups, the remaining continental groups tend to correspond to unique mixture components.



Starting with assigning each mixture component to its own cluster and an initial = 0.1437, we proceed with the steps of the PHM algorithm with τ = 0, i.e., merging until all observations belong to a single cluster and  is decreased to 0.
Figure <ref> visualizes the merging process as a dendrogram. The numbers in each node indicate the reduction of , i.e., Δ, by combining the corresponding branches into a single cluster. 
A smaller value of Δ indicates that the clusters being combined are more distinct. 
Because the algorithm always prioritizes merging the most overlapping clusters, the merging sequence reflects the relative genetic dissimilarities between different clusters. This relationship can be straightforwardly interpreted from the dendrogram.
For example, the first two merges — reducing  to 0.063 and 0.021, respectively — form a cluster representing samples from Europe, the Middle East, and Central/South Asia, reflecting a close genetic relationship and noticeable genetic admixture among these population groups.
The general structure of the dendrogram can be roughly explained by the likely path of historical human migrations: from Africa into the Middle East, from the Middle East to Europe and Central/South Asia, from Central/South Asia to East Asia, and from East Asia to Oceania and the Americas. 
The overall pattern of human genetic diversity among the continental groups identified from our analysis is also corroborated by more sophisticated genetic analysis using additional information (e.g., haplotype analysis) <cit.> and high-coverage genome sequencing data <cit.>.



  §.§.§ Cluster Analysis of Single-cell RNA Sequence Data


In this illustration, we apply the PHM algorithm to single-cell RNA sequencing (scRNA-seq) data from a sample of peripheral blood mononuclear cells to identify the different cell types. 
The data (sequenced on the Illumina NextSeq 500 and freely available from 10x Genomics) consist of gene expression counts for 2,700 single cells at 13,714 genes.

The raw sequence data are quality-controlled and pre-processed using standard procedures implemented in the  package <cit.>, leaving us with 2,638 cells. 
Principal component analysis is subsequently performed on the normalized and scaled expression counts for dimension reduction. The first 10 components are selected based on the elbow point of the scree plot (Supplementary Figure 5).
The resulting 2,638 × 10 data matrix is used for our cluster analysis. 
Additionally, each cell is annotated as one of the following cell types: Naïve CD4+ T cells, Memory CD4+ T cells, CD8+ T cells, Natural Killer (NK) cells, B cells, CD14+ monocytes, FCGR3A+ monocytes, Dendritic cells (DC), and Platelets, using known biomarkers. 
The annotated cell type information is not used in our cluster analysis procedure.

We fit a GMM on the dimension-reduced PC matrix and select the model with 9 components based on the optimal BIC. The posterior component assignment probability π_k for each cell is visualized in Figure <ref>. With the exception of the Naïve and Memory CD4+ T cells, each cell type corresponds to a unique mixture component.

Starting with each mixture component as its own cluster and an initial = 0.0703, we perform the merging procedure until all components have been combined into a single cluster. The merging process is represented by the dendrogram in Figure <ref>. 
We note that the pattern shown in the dendrogram has a striking similarity to the known immune cell differentiation trajectories. 
For example, the first merge combines the two types of CD4+ T cells, reducing  from 0.070 to 0.042.
Furthermore, CD8+ T cells, CD4+ T cells, NK cells, and B cells—all derived from Lymphoid progenitor cells—are grouped together in the dendrogram and separated from the branch consisting of FCGR4A+, CD14+ monocytes, and Dendritic cells—all of which are derived from Myeloid progenitor cells.





§ CONCLUSION AND DISCUSSION


In this work, we introduce the Distinguishability criterion, , to quantify the separability of clusters inferred from cluster analysis procedures.
We discuss the intuition behind the criterion, as well as the derivation and properties of .
We propose a combined loss function-based computational framework that integrates the Distinguishability criterion with available model and heuristics-based clustering algorithms and demonstrate its use with synthetic and real data applications.


The proposed  is an internal clustering validity index to assess the separability of the clustering results, with unique advantages over alternative validity indices. 
Since  is measured on the probability scale, the threshold for our proposed constrained optimization problem is interpretable. 
Additionally,  measurements are directly comparable across different datasets and various types of clustering applications,  enabling future work on assessing the replicability of clustering analysis.


While our numerical illustrations primarily use Gaussian or mixtures of Gaussian distributions to evaluate , it is important to highlight the proposed computational framework's flexibility and compatibility with any valid parametric likelihood function. 
As a result, the applications of  can be extended to a more diverse class of latent variable models, e.g., latent Dirichlet allocation (LDA) and generalized factor analysis models, in order to help address the similar model selection problems that arise with their use.
We will explore these extensions in our future work.


The properties of  are best utilized in our proposed PHM algorithm, which is motivated by Baudry et al.'s entropy criterion-based merging procedure <cit.>.
By combining mixture components into clusters, both algorithms enjoy excellent model fit and interpretability for the inferred clusters. 
However, the constrained optimization formulation described above leads to a more interpretable stopping rule for the PHM algorithm compared to Baudry et al.'s procedure.
Furthermore, the PHM algorithm shows improved computational efficiency—due to the cluster merging property of , the complexity of the PHM algorithm is invariant to the sample size of the observed clustering data, making it more suitable for analyzing large-scale data.
Finally, we note that the dendrogram visualizing the complete merging procedure can be used to represent a coalescent process that has many applications across scientific disciplines, such as developmental biology, human genetics, and evolutionary biology. 
One possible future research direction is shifting the focus from cluster analysis to uncovering the underlying coalescent trees by incorporating more context-specific information, with the PHM algorithm as a natural starting point.

It is also possible to extend the PHM algorithm to work with a general class of hard clustering algorithms that output optimal partitions of observed data.
One strategy is to treat each output partition as a distinct population sample and estimate its distribution using a finite mixture model. We can then re-normalize the mixture proportions globally and initilize the clusters at the level of the output partitions, at which point applying the PHM algorithm is straightforward.
This simple strategy extends the applications of the Distinguishability criterion to a more diverse class of clustering algorithms, including density-based and graph clustering algorithms.



§ METHODS




 §.§ Overview of Clustering Methods


We briefly review existing clustering methods, focusing on the algorithms used in this paper. 
The available clustering methods can be roughly classified into two categories: heuristics-based clustering methods, represented by the k-means and hierarchical clustering algorithms, and model-based clustering methods, represented by finite-mixture models. 
The heuristics-based clustering methods typically do not make explicit distributional assumptions and instead perform hard clustering by outputting optimal partitions of the observed data based on their corresponding objective functions. 
The model-based clustering methods perform formal statistical inference of the latent cluster structures underlying the observed data.



  §.§.§ k-means clustering


k-means clustering is a widely used heuristics-based clustering algorithm for partitioning a sample into clusters. 
As the name suggests, this procedure identifies clusters by their centroids (defined as the mean of all points in the cluster) and assigns observations to the cluster with the nearest centroid. More formally, for a specified number of clusters K, the algorithm groups the observations into disjoint sets 𝒞_1, …, 𝒞_K to minimize the distortion function, i.e., 


    min_𝒞_1, …, 𝒞_K∑_k=1^K ∑_i ∈𝒞_k |  | _i - ∑_j ∈𝒞_k_j / |𝒞_k|  |  |_2^2,


which corresponds to minimizing the within-cluster variances. The clustering is usually performed using Lloyd's algorithm <cit.>, which alternates between assigning observations to the cluster with the nearest centroid and updating the cluster centroid to eventually arrive at a locally optimal solution. 

It has been shown that the k-means algorithm is equivalent to an approximate maximum likelihood procedure, where the underlying probability model assumes a multivariate Gaussian distribution for each latent cluster <cit.>.
Hence, the observation frequency of each cluster and the corresponding parameters for the Gaussian distribution can be estimated straightforwardly using the partitioned output from the k-means algorithm. 



  §.§.§ Hierarchical clustering


Hierarchical clustering produces a sequence of partitions for observations. In the commonly used agglomerative hierarchical clustering procedure, each observation comprises its own cluster in the initial partition. Subsequent clusterings in the sequence are formed by repeatedly combining the most similar pair of groups of observations into a single group until all observations have been combined into a single cluster.
The similarity between groups of observations in hierarchical clustering is typically defined by a distance measure between pairs of observations, such as the squared Euclidean distance, in addition to a function that generalizes this similarity to groups of observations (referred to as a linkage function). Commonly used linkage functions between groups of observations are single linkage, complete linkage, average linkage, and the Ward linkage.

The hierarchical clustering algorithm is also often connected to probability models assuming a Gaussian data distribution underlying each cluster <cit.>, and explicit Gaussian assumptions are commonly used for model-selection or post-selection inference in hierarchical clustering <cit.>.  



  §.§.§ Model-based Clustering


Finite mixture models are the most representative approach for model-based clustering, where a mixture distribution with finite components models the observed clustering data.
Traditionally, each mixture component is taken to correspond to a homogeneous subpopulation or cluster. 
The goal of inference in the mixture model setting is to estimate the characteristics of each mixture component, i.e., p(|θ = k), and the corresponding mixture proportion, i.e., α_k.
The expectation-maximization (EM) algorithm is commonly used to find maximum likelihood estimates of the parameters of interest. 
The Gaussian mixture model (GMM), which models each mixture component using a unique Gaussian distribution, is probably the most commonly used mixture model in practice.
This is because, among other reasons, the GMM is considered to be a universal approximator <cit.> and is flexible enough to fit diverse data types.

Unlike most heuristics-based clustering approaches, model-based clustering algorithms perform soft (or fuzzy) clustering, as they do not directly partition the observed data. 
Instead, every data point has an associated (posterior) probability distribution over all possible cluster assignments.
A post-hoc classification procedure, with a pre-specified decision rule using the cluster assignment probabilities, can be applied to partition the observed data.

For a more thorough review of model-based clustering, we refer the reader to <cit.>



 §.§ Derivation and Estimation of 
  

The Bayes risk for a general classifier δ(): ℝ^p ↦{1, …, K} assigning a point  to one of the K clusters can be derived as follows,


    = 𝔼_[ 
            𝔼_θ[ 
                L(δ(), θ)
                |] ]  
        =  ∫( ∑_j=1^K ∑_i ≠ jπ_i()  ( δ() = j |) )  P(d) 
        =  ∫( ∑_j=1^K  (1 - π_j())  ( δ() = j |) )  P(d)


For a more detailed derivation, see Appendix <ref>. The marginal data distribution P(), with loose notation, is given by


    P() = ∑_k=1^K (θ=k)   p(|θ = k) = ∑_k=1^K α_k   p(|θ = k).



For the default randomized decision rule, δ_r() ∼  ( ()),


    (δ() = j |) = π_j().


Thus,


    P_ mc,  δ_r   = ∫( ∑_j=1^K  π_j() (1-π_j()) ) P(d) 
        = 2 ∑_i<j∫π_i() π_j()  P(d)


To compute  using the optimal decision rule, δ_o() = max_k π_k(), we define a partition of the sample space, ∪_k=1^K ℛ_k, such that ℛ_k := {: δ_o() = k}. 
It follows that, 


    (δ() =j |) = 1{∈ℛ_j},


and  


    P_ mc,  δ_o    = ∫( ∑_j=1^K ∑_i ≠ jπ_i()  1{∈ℛ_j}) P(d) 
        = ∫(1 - max_k π_k())   P(d )


Note that, 


    ∑_j=1^K  π_j() (1-π_j()) ≥∑_j=1^k π_j() (1-max_k π_k())  = 1 - max_k π_k(), ∀


Hence, 


    P_ mc,  δ_r≥ P_ mc,  δ_o


Unless otherwise specificied, we use the notation  to refer to P_ mc,  δ_r by default.
In this paper, we estimate  by plugging in the point estimates of the α_k's as well as the key distributional parameters in the corresponding likelihood functions obtained from the observed data.



 §.§ Lower and Upper Bounds of 


When all clusters are well-separated,  approaches its lower bound at 0. 
More specifically, assuming well-separated clusters, both of the following conditions should hold:

 
    π_i() π_j() →  0,  ∀  (i ,j) , 
    and   max_k π_k() → 1,  ∀ .


Hence, both decision rules (δ_r and δ_o) approach perfect classification accuracy.

 is maximized when all clusters are completely overlapping, i.e., 


    p(|θ = i) = p(|θ = j)   ∀  (i ,j) .


Thus, π_k(x) = α_k,  ∀.
It follows that 


    max P_ mc,  δ_r = ∑_k=1^K α_k   (1 - α_k),


and 


    max P_ mc,  δ_o = 1 - max_k α_k.


In the special case that α_k = 1/K for all k values,


    max P_ mc,  δ_r = max P_ mc,  δ_o = K-1/K.
 



 §.§ The Cluster Merging Property of 


The merging property is specific to the default , evaluated using the randomized decision rule δ_r.
For a given cluster configuration with K ≥ 2, consider merging two arbitrary clusters to form a new combined cluster.
Let  and ^† denote the misclassification probabilities before and after the merging, respectively. 
The following proposition summarizes the merging property.

 
   Merging two existing clusters indexed by i and j leads to  
   
    Δ^(i,j) :=  - ^†  = 2 ∫π_i () π_j ()   P(d) ≥ 0.
 
   Furthermore, 
   
    = ∑_i<jΔ^(i,j)




    Appendix <ref>.


The cluster merging property forms the basis of the PHM algorithm. As Δ^(i,j) can be pre-computed for all pairs of clusters from the initial configuration, the subsequent updates for —merging one pair of clusters at a time—become straightforward to compute. 




 §.§ Numerical Evaluation of 


Evaluating  numerically can be challenging, especially when clustering data are high-dimensional.
For low-dimensional data, it is possible to evaluate Eqn (<ref>) by numerical integration using various quadrature methods. 
However, they generally do not scale well when the clustering data dimensionality becomes larger than 5. 
When the marginal data distribution, P(), can be directly sampled from (as in the case in all examples presented in this paper), Monte Carlo (MC) integration becomes an efficient solution. 
Specifically, we sample M data points from P() and approximate  by


    P̂_ mc = 1/M∑_i=1^M( ∑_j=1^K  (1 - π_j(_i) )  ( δ(_i) = j |) )


The unique advantage of the Monte Carlo integration method is that its error bound is always O(1/√(M)) regardless of the dimensionality of . 
We provide comparisons between the Monte Carlo method and the numerical integration for evaluating  in some low-dimensional settings (Supplementary Method 1 and Supplementary Table 4), indicating that the MC integration method is accurate and efficient. 



unsrt
 




    
    

§ DERIVATION OF 
  
            
    To compute the Bayes risk for a general classifier δ() under the 0-1 loss, we first evaluate its posterior expected loss 𝔼_θ[ L(δ(), θ) |] as follows: 
    
    𝔼_θ[ L(δ(), θ) |]     = ( θδ() |) 
        = ∑_j=1^K (δ() =j,  θ j |) 
        = ∑_j=1^K ( θ j  | )  (δ() = j |) 
        =  ∑_j=1^K ∑_i  j(θ = i |  )  (δ() = j |)  
        = ∑_j=1^k ∑_i  jπ_i() (δ() = j |) 
        =  ∑_j=1^k  (1-π_j())  (δ() = j |)

    
    Note that (θ j |, δ()) = (θ j |).
    Subsequently,
    
    = 𝔼_[ 
                        𝔼_θ[ 
                            L(δ(), θ)
                            |] ] 
    
                            =  ∫( ∑_j=1^K ∑_i ≠ jπ_i()  ( δ() = j |) )  P(d) 
        =  ∫( ∑_j=1^K  (1 - π_j())  ( δ() = j |) )  P(d)

            
    

§ PROOF OF PROPOSITION 1
 
    
    
    Consider merging two existing clusters i,j to a new combined cluster k'. It follows that

    
    α_k' = α_i + α_j
    and    p( |θ=k) = α_i   p( |θ=i) + α_j   p( |θ=j)/α_k'

    
    Consequently, by applying Bayes rule,
    
    
    π_k'() = π_i() + π_j()

            
    Let S denote the set of indices of the existing clusters not impacted by the merge, where |S| = K-2. By Eqn (<ref>),  can be written as 

    
    = 2 ∑_m,n ∈ S,   m < n∫π_m () π_n () P(d)  
        + 2 ∑_l ∈ S∫π_l () π_i () P(d) +  2 ∑_l ∈ S∫π_l () π_j () P(d) 
        + 2 ∫π_i () π_j () P(d)
 

    By Eqn (<ref>),
    
    2 ∑_l ∈ S∫π_l () π_i () P(d) +  2 ∑_l ∈ S∫π_l () π_j () P(d) = 2 ∑_l ∈ S∫π_l() π_k' ()  P(d)


    and note that, 

    
    ^† = 2 ∑_m,n ∈ S,   m < n∫π_m () π_n () P(d)  +  2 ∑_l ∈ S∫π_l() π_k' ()  P(d)


    It becomes evident that 
        
    Δ^(i,j) =  - ^† = 2 ∫π_i () π_j () P(d) ≥ 0


        Plugging in the expression of Δ^(i,j) into Eqn  (<ref>) yields
        
    = ∑_i<jΔ^(i,j)
            
    
            
    

  
Remark The merging property is specific to the randomized decision rule δ_r under the 0-1 loss. For the optimal decision rule, δ_o, it can be shown that ^†≤ after merging a pair of existing clusters. 
    However, the quantitive expression for Δ^(i,j) is analytically intractable. 
            
    

§  DENDROGRAM CONSTRUCTION

            
    We use a dendrogram to visualize the PHM procedure described in Algorithm <ref>. The leaf nodes in the tree correspond to the individual components of the mixture model (MM) used to fit the data.
    The edges between parent and child nodes correspond to a cluster merge step; we present the Δ reduction from the merge on the parent node for the merged clusters.
    In this way, it is possible to determine the value of  after each successive merge by subtracting the cumulative Δ values from the leaf nodes of the tree up to the height of a specific merge from the  of the initial cluster configuration.
            
    The height of a merge in the tree is determined as follows. 
    Let ^0 denote the  value at the initial cluster configuration and, for a given merge, let ^ be the value of the criterion prior to that merge taking place.
    We place the merge in the dendrogram at a height corresponding to the log_10 scaled ratio of these values, i.e., log_10^0 / ^.
    The log_10 transformation prevents the merges from early on in the procedure from being “squashed" to the bottom of the tree.
    We opt to use the  value before the merge occurs rather than after to avoid dividing by zero when calculating the height of the final merge (which would result in = 0). 

    
















Supplementary Figures
























Supplementary Tables




















Supplementary Methods





§ MONTE CARLO  COMPARISON




Here we compare  computed using standard cubature methods to P̂_ mc estimated using a Monte Carlo (MC) integral (as described in Section 4.5). We use 50 replicates to obtain the timing measurements and quantify the uncertainty in the MC integral. All values presented are means over these replicates unless otherwise specified. For the MC integration procedure we use M = 10^5 sample points, and parallelize the computation over 8 cores.

The distribution in question consists of three Gaussian clusters with π_k = 1/3 in ℝ^p with I_p variance. The cluster are centered at 0^p and ± d^p, where d^p is the p-dimensional vector whose elements are all d. d is set so that the Euclidean distance between d^p and the origin is fixed to be 3; i.e. d = √(3^2/p). This is so that the value of  remains fixed across dimensions and we do not need to worry about the dimensionality affecting the true value of . The results for dimension p = 1, …, 5 are presented in Supplementary Table <ref>.

Both approaches produce highly similar values of  across dimensions. The standard deviation of the MC estimates is quite low, indicating stability in the estimation procedure. However, while the cubature method evaluation time rapidly increases in p, the time for the MC procedure is roughly constant. These together highlight the MC estimation procedure as a viable and accurate approach to estimate , especially in moderate dimensional data (i.e., p ≥ 3) where cubature methods may struggle to reach a solution in a reasonable time.


