Selecting the number of clusters remains a fundamental challenge in unsupervised learning. Existing approaches typically focus on identifying a single "optimal" partition, often overlooking statistically meaningful structure present across multiple resolutions. We introduce ElbowSig, a general inferential framework for assessing clustering structure over a range of resolutions. The method formalizes the elbow heuristic by defining a normalized discrete curvature statistic based on the sequence of within-cluster heterogeneity values, and evaluates its significance relative to a null distribution of unstructured data. This yields hypothesis tests across resolutions, enabling simultaneous inference at multiple clustering scales. We derive the asymptotic behavior of the null statistic in both large-sample and high-dimensional regimes, characterizing its limiting form and variability. Because it depends only on the heterogeneity sequence, ElbowSig is compatible with a wide range of clustering algorithms, including hard, fuzzy, and model-based methods. Experiments on synthetic and real datasets show that the procedure controls Type-I error under unstructured data while providing power to detect multiscale organization, revealing structure that is often missed by single-resolution selection criteria.
Table 1: Point-wise heterogeneity measures hk(xi;θ) for the clustering models considered in this paper. For hard clustering (agglomerative and k-means) and FCM, μj denotes the centroid of cluster Cj . For FCM, wij∈[0,1] are the membership weights satisfying ∑j=1kwij=1 , m>1 is the fuzzifier exponent. For GMM, each component has mixing weight πj ( ∑j=1kπj=1 ), mean μj , and covariance Σj ; ϕ(xi∣μj,Σj) denotes the corresponding Gaussian density.
Figure 1: Mean value and variance of the elbow statistic δk(r) for unstructured data as a function of the dimension D for agglomerative, k-means, GMM and FCM clustering methods. The expectation and variance are computed from NR=500 reference datasets of size N=30 uniformly distributed on [0,1]D , assuming k=3 clusters. The dashed lines indicate the predicted asymptotic trends for large D . The dotted line in panel (a) indicates the asymptotic value of E[δk(r)] for FCM in the limit of large fuzzifier parameter m .
Figure 2: Clustering of two-dimensional data. (a) Dataset with N=200 observations consisting of M=3 Gaussian components centred at the locations indicated by crosses, each with spread σc=1 . Solid circles show data points coloured according to their generating component. Grey crosses show one realization of a bounding-box uniform reference dataset. Partitions with k∈{1,2,…,10} were evaluated using five standard methods for estimating the number of clusters: (b) Calinski–Harabasz, (c) Davies–Bouldin, and (d) silhouette scores as functions of k . Vertical dashed lines indicate the value k^ at which each index attains its maximum. (e) Gap statistic as a function of k , with vertical lines marking the optimal number of clusters according to the Gap (I) and Gap (II) criteria described in the text. Error bars indicate the standard error estimate sk ( Tibshirani et al., 2001 ) .
Figure 3: Application of ElbowSig to the two-dimensional dataset shown in Fig. 2 . (a) Total within-cluster heterogeneity (inertia) as a function of k . (b) Observed elbow statistic (circles) together with the 100(1−psig) percentile of the δk distribution obtained from unstructured reference datasets (squares) at level q1=0.05 . (c) Empirical p-values pk comparing the observed elbow statistic at each k with its reference distribution. The horizontal dashed line indicates the per-scale significance threshold psig .
D
M
σc
Method
Number of clusters, k
1
2
3
4
5
6
7
8
9
10
2
3
1
ElbowSig (per-k)
–
28
82
0
3
13
0
1
5
2
ElbowSig (FDR)
–
25
80
0
1
11
0
1
5
2
Gap (I)
6
12
81
1
0
0
0
0
0
0
Gap (II)
0
8
76
1
2
5
0
1
3
4
CH
–
16
70
0
1
3
0
1
4
5
Table 2: Results for synthetic clustered datasets with N=200 observations across dimensions D∈{2,5,20} , with M Gaussian components randomly centred in [−10,10]D and within-cluster standard deviation σc as indicated. Each entry reports how many of the 100 simulated datasets led a method to select the corresponding number of clusters k . Results for the ElbowSig and Gap methods are based on 200 PCA-aligned uniform reference datasets. ElbowSig results are shown for both the individual scale (per- k ) and global (FDR) significance criteria. Gap (I), Gap (II), Calinski–Harabasz (CH), Davies–Bouldin (DB), and the Silhouette index are included for comparison. Dashes at k=1 for all methods except the gap-based ones indicate that those methods cannot return k=1 by definition.
Clust. method
Reference
Control
Number of clusters, k
2
3
4
5
6
7
8
9
10
Agglomerative
BBU
per- k
49
100
2
4
43
9
1
10
6
FDR
47
100
2
3
41
7
1
8
5
PCA
per- k
2
99
0
0
14
2
1
9
4
FDR
2
99
0
0
12
1
0
7
3
K-means
BBU
per- k
38
99
14
9
25
5
7
1
3
Table 3: Effect of clustering algorithm, reference dataset construction, and multiple-testing correction on the ElbowSig outcomes. Shown are results from 100 datasets generated from M=3 Gaussian components with σc=1 in D=5 dimensions. Each realization consists of N=200 observations using randomly located centres of the Gaussian components in [−10,10]5 . “BBU” indicates bounding-box uniform reference datasets, and “PCA” denotes PCA-aligned uniform reference datasets. The numbers in each cell give how many datasets led a method to select the corresponding number of clusters k .
D
Data
Method
Reference
Number of clusters, k
1
2
3
4
5
6
7
8
9
10
2
Uniform
ElbowSig (per-k)
BBU
(78)
2
5
4
4
0
1
1
1
4
ElbowSig (FDR)
BBU
(97)
1
0
0
0
0
0
0
0
2
Gap (I)
BBU
99
1
0
0
0
0
0
0
0
0
Gap (II)
BBU
11
3
14
19
14
6
9
8
3
13
SigClust
Gaussian
27
-
-
-
-
-
-
-
-
-
Table 4: Results for unstructured data in dimensions D∈{2,20} generated from either a Uniform distribution on [0,1]D or a Gaussian distribution with variance σc=1 in each dimension. For each configuration, 100 independent datasets were simulated. Each entry reports how many of the 100 datasets led a method to select the indicated number of clusters. For ElbowSig, the value displayed in parentheses in the k=1 column denotes the number of datasets for which no statistically significant clustering was detected at any tested value k>1 . Results are shown for both the individual scale (per- k ) and global (FDR) significance criteria. Gap (I) and Gap (II) correspond to the two gap-statistic selection rules described in the main text. For SigClust, only the number of cases in which the null hypothesis of a single Gaussian component was not rejected is shown.
Data
N
D
Classes
Clustering
Reference
k (per- k )
k (FDR)
Iris
150
4
3
Agglom.
BBU
2,3,5
2,3,5
PCA
3
–
k-means
BBU
2,3,5
2,3,5
PCA
3
3
GMM
BBU
2,3,5,6,8
2,3,5,6,8
PCA
2
2
Table 5: ElbowSig results for real datasets using several clustering algorithms and reference constructions.
Determining the number of clusters is a central challenge in unsupervised learning, where ground-truth labels are unavailable. The Silhouette coefficient is a widely used internal validation metric for this task, yet its standard micro-averaged form tends to favor larger clusters under size imbalance. Macro-averaging mitigates this bias by weighting clusters equally, but may overemphasize noise from under-represented groups. We introduce Composite Silhouette, an internal criterion for cluster-count selection that aggregates evidence across repeated subsampled clusterings rather than relying on a single partition. For each subsample, micro- and macro-averaged Silhouette scores are combined through an adaptive convex weight determined by their normalized discrepancy and smoothed by a bounded nonlinearity; the final score is then obtained by averaging these subsample-level composites. We establish key properties of the criterion and derive finite-sample concentration guarantees for its subsampling estimate. Experiments on synthetic and real-world datasets show that Composite Silhouette effectively reconciles the strengths of micro- and macro-averaging, yielding more accurate recovery of the ground-truth number of clusters.
Aggelos Semoglou, Aristidis Likas, John Pavlopoulos
Athens University of Economics and Business, Greece. · Archimedes, Athena Research Center, Greece. · University of Ioannina, Greece.
Unsupervised clustering is widely used to discover underlying group structure, after which researchers often seek to validate these groups via post-hoc hypothesis testing on the population means. However, using the same data both to form clusters and to test for a difference in means inflates the Type I error rate. Recent work in selective inference has proposed solutions when the underlying data are drawn from a multivariate Gaussian distribution with a known covariance matrix or an unknown isotropic covariance matrix. We propose a new selective test for the difference in means between two clusters that provably controls the Type I error rate under an arbitrary unknown covariance matrix. Our test statistic does not require estimating the feature covariance matrix; instead, it leverages a Monte Carlo sampling procedure over the set of orthogonal matrices. We show that this new selective test can be computed efficiently, and we derive closed-form solutions when clusters are formed using k-means or hierarchical clustering. In simulations, we demonstrate that our proposed p-value controls the selective Type I error rate and achieves substantial power, and we illustrate its application on gene expression and drug review data sets.
Yiqun T. Chen, Lucy L. Gao
Departments of Biostatistics and Computer Science, Johns Hopkins University · Department of Statistics and Actuarial Science, University of Waterloo
Deep clustering is a powerful approach for discovering meaningful structures in high-dimensional data by learning a low-dimensional latent representation prior to clustering. Despite its empirical success, assessing the statistical reliability of the resulting clusters remains challenging. Testing discovered clusters on the same data induces selection bias and invalidates classical p-values. Selective inference (SI) provides a principled framework for correcting this bias, but existing methods focus on clustering performed directly on the observed features. In this work, we develop an SI framework for deep clustering with a fixed pretrained encoder. The key challenge is that cluster assignments are determined through a nonlinear transformation from the original data space to the latent space, resulting in a substantially more complex selection process than in conventional clustering. Our method provides a computationally tractable way to account for this process and enables valid statistical testing of differences between clusters identified in the latent space. Synthetic experiments demonstrate that the proposed method controls the Type I error rate while achieving higher power than valid but conservative baselines, and genomic applications show that it can identify significant cluster differences while appropriately accounting for selection bias. Our framework provides a principled approach to quantifying the statistical reliability of structures discovered by deep clustering.
Eina Mizui, Tomohiro Shiraishi, Shunichi Nishino +1