Cluster Validity, Stability, and Choosing the Number of Clusters

The question “how many clusters are in this data” has no answer that is independent of the definition of a cluster. For a centroid objective, more groups always reduce the objective; for a density criterion, groups appear and disappear as the radius crosses a gap; for a hierarchical method, every level of the tree is a valid partition. An impossibility result makes the difficulty precise: no clustering function can satisfy scale invariance, richness (every partition is recoverable for some distance matrix), and consistency (a transformation that shrinks within-cluster distances and expands between-cluster distances does not change the output) simultaneously. Validation therefore always measures agreement between a partition and some particular expectation about what clusters look like.

Internal indices and the geometry they assume

The silhouette of a point compares its mean distance to the members of its own cluster, a(i), with the smallest mean distance to the members of any other cluster, b(i), producing

$$s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))}$$

which lies between -1 and 1, is close to 1 for a well-placed point, near 0 for a point on a boundary, and negative for a point that is closer to another group on average. Averaging over points gives a score for a partition, and the value for each k can be plotted. The cost is O(n^2) distances, since each point compares itself with every other point, unless cluster statistics are used to summarize the means.

Several other indices follow the same pattern with different aggregations. The Calinski-Harabasz criterion is a variance ratio: between-cluster scatter divided by within-cluster scatter, scaled by (n - k) / (k - 1). The Davies-Bouldin index averages, over clusters, the worst ratio of a within-cluster spread to the distance between that cluster’s centroid and another’s. Every one of these definitions is written in terms of distances to centroids or mean distances within a group, so they reward compact, roughly spherical, similarly sized groups.

That bias is the central trap. A density-based partition such as the one produced by DBSCAN describes crescents, rings, and elongated regions, and a silhouette computed on it scores those shapes poorly against a competing partition of two round blobs, because the nearest competing group on average may be a neighbouring part of the same crescent. The same applies to hierarchical partitions produced with single linkage, where chaining is a feature rather than a defect. The index must match the model class: geometric indices for centroid methods, and either external criteria or graph and density statistics for the others.

Indices are also not monotone in k, so a maximum is a genuine maximum rather than a truncation artefact. The within-cluster sum of squares used by k-means decreases monotonically, which is why its elbow is a heuristic reading of a curve and not a selection criterion. The silhouette and the variance ratio typically peak at some k, and the peak moves with the metric, so the same data standardized differently can produce a different preferred number of groups.

A reference distribution gives the gap statistic a null model

The gap statistic replaces the elbow reading with a comparison against a null. The within-cluster scatter W_k is computed for the observed data, and the same quantity is computed for several reference datasets drawn uniformly over the data’s bounding box, or over a principal-component-aligned box for the rotated case. The statistic is

$$Gap(k) = E^*[\log W_k] - \log W_k$$

where the expectation is over the reference samples, and the reported number of clusters is the smallest k satisfying Gap(k) >= Gap(k+1) - s_{k+1}, with s the standard error of the Monte Carlo estimate. The standard error term stops the rule from chasing small, noisy improvements after the true structure stops being added. The test is a comparison against a specific null hypothesis, namely uniform data with no clusters, and it loses its meaning if the data are far from that reference geometry; heavy tails and strong dependence between coordinates both shift the reference scatter.

A model-based alternative replaces the geometric index with a likelihood. Fitting mixtures with a range of component counts and comparing penalized likelihoods, as described in Gaussian Mixture Models and the EM Lower Bound, gives a criterion that is defined for the model that generates the partition, at the price of relying on asymptotic approximations that the degeneracy of the likelihood complicates.

Resampling, agreement, and what stability does not prove

A resampling approach asks whether a partition is reproducible instead of whether it is geometrically good. Split or bootstrap the data, cluster each resample, and compare the partitions on the observations that both runs share. The comparison needs a chance correction, because two arbitrary partitions of the same points agree substantially by accident when k is small. The adjusted Rand index subtracts the expected agreement under a random model, so a value of 0 means the partitions are no better than random and 1 means they coincide. Consensus clustering goes further and aggregates many runs: the entry of a co-association matrix is the fraction of runs in which two points fall into the same group, and the matrix itself is clustered to get a final partition and a picture of which points are consistently grouped.

The interpretation of stability is where the method is usually overstated. Stability measures agreement among clusterings produced from similar data under one algorithm and one assumed k, so it cannot distinguish a genuine cluster from a region the algorithm always splits the same way. For several objective functions, stability can be maximized by a solution that separates a single cluster into two parts, and the two-part solution may be more stable than the correct one. The theory that makes stability a consistent selector needs assumptions about the algorithm and the data that are not under the analyst’s control, and without them the criterion is a diagnostic rather than a proof. Splitting the data and checking whether the labels transferred from one half still hold in the other, which is the prediction-strength idea, narrows the question: it measures whether the partition’s own boundaries generalize, not whether the groups exist.

Practical ordering of the evidence

The order that avoids wasted work is to fix the model class first, then the geometry, then the count. If the groups are expected to be convex and similar, a centroid method, an internal index, and a stability check all speak the same language, and disagreement between them is informative. If the groups are expected to be non-convex or to differ in density, no centroid-based index should be consulted, and the parameters of the density criterion should be swept instead. If the constraint is a downstream task, the honest criterion is performance on that task, since no internal index can rank partitions by usefulness for a purpose it does not know about.

Whichever route is taken, three checks cost little. Report more than one index and record where they disagree, because a single number hides the geometry it rewarded. Verify that the preferred k does not move when the parameters of the representation move, such as the number of retained components in a linear projection or the neighbourhood size in an embedding. And check the sensitivity of the answer to the random seed of any iterative method, since the difference between two local minima can exceed the difference between two candidate values of k.

import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, adjusted_rand_score

def summary(X, ks):
    out = []
    for k in ks:
        km = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
        out.append((k, silhouette_score(X, km.labels_), km.inertia_))
    return out                                   # silhouette has a maximum; inertia never does

def stability(X, k, rng, rounds=20, frac=0.8):
    n = len(X)
    scores = []
    for _ in range(rounds):
        a = rng.choice(n, int(frac * n), replace=False)
        b = rng.choice(n, int(frac * n), replace=False)
        common = np.intersect1d(a, b)
        la = KMeans(n_clusters=k, n_init=5, random_state=0).fit(X[a]).predict(X[common])
        lb = KMeans(n_clusters=k, n_init=5, random_state=0).fit(X[b]).predict(X[common])
        scores.append(adjusted_rand_score(la, lb))   # chance-corrected agreement
    return np.mean(scores), np.std(scores)

Two further limitations apply to every index discussed here. None of them detects structure that the chosen algorithm cannot represent, so a value of 1 for the silhouette at k = 2 on strongly non-convex data means only that two round groups fit well. And none of them distinguishes a real but small group from a group created by a few outliers, because nothing in the arithmetic knows how many points a group must contain to be interesting; that judgement belongs to the model of the data, not to the index.