Clustering

Clustering divides a set of unlabelled items into groups, the clusters, so that items in the same group are more alike than items in different groups. It is one of the main tasks of unsupervised machine learning: there are no labels to fit, only a set of feature vectors and a notion of which of them are similar. In computer vision the items are pixels described by colour and position, local descriptors, image patches, bounding-box shapes or whole-image embeddings, and the clusters become palettes, segments, superpixels, visual vocabularies, anchor shapes or pseudo-labels for training a network.

There is no single clustering problem. A method must decide what “alike” means, at what scale two items stop being alike, and what shape a group may take, and different methods answer differently.

Definition

Given $n$ items, a hard clustering assigns each item to exactly one of $k$ disjoint groups whose union is the whole set, a partition. A soft clustering gives each item a degree of membership in every group, for example the posterior probabilities of a mixture model. A hierarchical clustering returns a nested sequence of partitions, from $n$ singletons up to one group holding everything. Some methods also leave items unassigned and call them noise.

What makes a grouping good is not part of the definition. Hastie, Tibshirani and Friedman (2009, §14.3) put the goal as grouping objects so that those in one cluster are more closely related to one another than to objects in other clusters, and they stress that the notion of similarity this requires can only come from subject-matter considerations, much as a loss function does in supervised learning.

Intuition

Look at a scatter of points and the clusters seem obvious, but each judgement uses assumptions. Two blobs separated by a gap look like two clusters at one magnification and like two specks of one cluster when zoomed out. A ring around a disk is two clusters if what matters is connectedness, and one cluster with a hole if what matters is distance to a centre. Two long parallel bars a few pixels apart are two clusters if one follows each bar, and the same bars could be cut into short segments by a method that likes round groups. A sparse cloud next to a dense one is a cluster or noise depending on what density counts.

So a clustering method is a choice of representation (which features describe each item), dissimilarity (how two items are compared), scale (how far apart is “different”) and cluster model (round, elliptical, density-connected, graph-separated). Hastie et al. (§14.3.3) judge the choice of dissimilarity far more important for success than the choice of algorithm. In vision this is the difference between clustering pixel colours in RGB, where shading stretches one surface into a long streak, and in chromaticity, where it collapses to a point, as the k-means article shows.

No method satisfies every reasonable demand

Kleinberg (2002) made the difficulty precise. He defines a clustering function as any rule that maps the pairwise distances between $n$ points to a partition, and asks for three properties:

• Scale invariance: multiplying all distances by the same positive factor does not change the result; the method has no built-in length scale.
• Richness: every partition of the points is the output for some choice of distances.
• Consistency: if distances inside the output’s clusters are shrunk and distances between clusters are stretched, the output does not change.

His Theorem 2.1 states that for every $n \ge 2$ no clustering function has all three. Any two are attainable: single linkage stopped at $k$ clusters is scale-invariant and consistent, stopped at a fixed distance $r$ it is rich and consistent, and stopped at a fixed fraction of the largest distance it is scale-invariant and rich (his Theorem 2.2). He also shows that a family of centroid-based functions, which includes k-median and the k-means objective with centres chosen among the data points, is not consistent once $n$ is large enough relative to $k$ (Theorem 4.1). Kleinberg reads the result as a set of trade-offs inherent in clustering, among them whether the number of clusters is fixed in advance and whether a length scale is built in.

Formal Definition

Let the data be $X \in \mathbb{R}^{n \times d}$ with one item $\mathbf{x}_i$ per row (scikit-learn’s layout), or more generally a matrix of pairwise dissimilarities $d(\mathbf{x}_i, \mathbf{x}_{i'})$. A hard clustering into $k$ groups is an assignment $c : \{1, \dots, n\} \to \{1, \dots, k\}$, or equivalently a partition $\mathcal{C} = \{C_1, \dots, C_k\}$. Most combinatorial methods choose $c$ to make a loss small, for example the within-cluster point scatter (Hastie et al., §14.3.5)

$$W(c) = \frac{1}{2} \sum_{j=1}^{k} \sum_{c(i) = j} \sum_{c(i') = j} d(\mathbf{x}_i, \mathbf{x}_{i'}) ,$$

where, with squared Euclidean distance, each cluster’s term equals its size times its sum of squared distances to the cluster mean, which ties $W$ to the k-means objective. Because the total scatter is fixed by the data, minimizing $W$ is the same as maximizing the scatter between clusters. Exhaustive search is hopeless: the number of partitions of $n$ items into $k$ groups is already 34,105 for $n = 10$, $k = 4$, and about $10^{10}$ for $n = 19$, $k = 4$ (Hastie et al., §14.3.5), so practical methods descend greedily from a start, build the partition step by step, or optimize a relaxed problem.

The label values themselves carry no meaning. Any permutation of the $k$ cluster indices describes the same partition, which is why clusterings are compared through the contingency table $N = [n_{ij}]$, the number of items in true class $i$ and cluster $j$, rather than label by label.

Families of Methods

Hastie et al. (§14.3.4) divide clustering algorithms into combinatorial methods, mixture modelling and mode seeking; with hierarchical and graph-based methods added, the families below cover the methods vision uses.

Family	Representative methods	Number of clusters	Cluster shape it can represent	Main parameter
Centroid (partitional)	k-means, k-medoids	fixed $k$	convex cells around a centre	$k$, initialization
Mixture models	Gaussian mixture fitted by EM	fixed $k$	ellipsoids	$k$, covariance type
Mode seeking	mean shift	from the data	basins of a density	bandwidth
Density-connected	DBSCAN, OPTICS, HDBSCAN	from the data	any connected dense region; noise	radius and minimum count
Hierarchical	single, complete, average, Ward linkage	a whole tree; cut it	depends on the linkage	where to cut
Graph-based	spectral clustering, normalized cuts, Felzenszwalb–Huttenlocher	fixed $k$ or from a threshold	any set weakly connected to the rest	affinity graph

Centroid methods

k-means represents each cluster by a centre and assigns every item to the nearest one, minimizing the sum of squared distances to the centres; Lloyd’s alternation of assignment and mean updates descends to a local minimum, and k-means++ seeding and restarts lower the risk of a poor one. Its clusters are the Voronoi cells of the centres, so they are always convex, and the equal-distance boundary favours clusters of similar spread. k-medoids restricts each centre to be one of the cluster’s own items, which allows any dissimilarity and reduces the pull of outliers at the cost of a slower update (Hastie et al., §14.3.10).

Mixture models

A Gaussian mixture model treats the data as samples from a weighted sum of Gaussians, one per cluster, and assigns each item the posterior probability of each component, a soft clustering. Fitting by maximum likelihood with expectation–maximization estimates a covariance per component, so clusters can be elongated, tilted and of different spread; k-means is its limit with equal, vanishing spherical covariances.

Mode seeking and density-based methods

Mean shift estimates the density of the data with a kernel, moves every item uphill towards a mode, and puts items that end at the same mode in one cluster. The number of clusters follows from the bandwidth, and a cluster can have any shape that is the basin of one density peak.

DBSCAN (Ester, Kriegel, Sander & Xu, 1996) defines clusters as regions where points are dense. With a radius $\varepsilon$ (Eps) and a count MinPts, a point is a core point if its $\varepsilon$-neighbourhood, which includes the point itself, holds at least MinPts points; points within $\varepsilon$ of a core point are directly density-reachable from it (Definitions 1–2); and a cluster is a maximal set of points connected through chains of core points, with border points attached to the cluster of a nearby core point. Points in no cluster are noise (Definitions 5–6). The authors argue that partitioning methods are restricted to convex clusters, one per Voronoi cell, while density-connected clusters can have any shape. They suggest setting MinPts to 4 for two-dimensional data and reading $\varepsilon$ off the sorted distances to each point’s fourth nearest neighbour, at the first “valley”, and they note the cost of global parameters: they are meant to describe the least dense cluster, anything sparser is noise, and two clusters of different density that lie close together may be merged (§4). With a spatial index, they estimate the average cost as $O(n \log n)$.

scikit-learn 1.7.2 provides two descendants. OPTICS computes an ordering of the points and a reachability distance for each, so DBSCAN-like clusterings for a range of $\varepsilon$ can be extracted from one run; HDBSCAN builds a minimum spanning tree of a “mutual reachability” graph, in which the distance between two points is at least each one’s distance to its min_samples-th nearest neighbour, and extracts clusters across all density levels, so it needs no $\varepsilon$; components smaller than min_cluster_size become noise (scikit-learn User Guide §2.3.8–§2.3.9).

Hierarchical clustering

Agglomerative clustering starts with every item as its own cluster and repeatedly merges the two closest clusters, giving $n - 1$ levels; divisive clustering starts with one cluster and splits recursively (Hastie et al., §14.3.12). The result is a binary tree, usually drawn as a dendrogram with each merge at a height equal to the dissimilarity of the two merged groups, and cutting it at a height gives a partition. What “closest clusters” means is set by the linkage:

• Single linkage: the distance between the closest pair of items, one from each cluster.
• Complete linkage: the distance between the farthest such pair.
• Average (group-average) linkage: the mean distance over all such pairs.
• Ward linkage: merges so as to keep the within-cluster sum of squared differences small, a variance-minimizing, bottom-up counterpart of the k-means objective (scikit-learn User Guide §2.3.6).

The linkage decides which clusters can form. Single linkage needs only one short link to join two groups, so it follows thin, curved or elongated structures but tends to chain groups together through a few intermediate points, and its clusters can have very large diameters; complete linkage produces compact clusters but can place an item closer to another cluster’s members than to some of its own; group average lies between them (Hastie et al., §14.3.12). Single and complete linkage depend only on the order of the dissimilarities, so any monotone transformation of them leaves the result unchanged, while group average does not have this invariance. Single linkage is a graph method in disguise. Kleinberg (2002) describes it as adding the edges of the complete distance graph in order of increasing length and reading the clusters off the connected components; in our reading, this is Kruskal’s algorithm for the minimum spanning tree, so the single-linkage clusters at $k$ are the components left after deleting the $k - 1$ longest edges of that tree, and scikit-learn 1.7.2 computes single linkage from a minimum spanning tree (our reading of its source).

Graph-based and spectral methods

Graph methods treat items as the nodes of a weighted graph whose edges say how similar two items are, and look for sets of nodes that are strongly connected inside and weakly connected to the rest. Normalized cuts measure a cut relative to the total connection of each side and approximate the best partition with generalized eigenvectors of the graph Laplacian. Spectral clustering generalizes this: compute the first $k$ eigenvectors of a graph Laplacian, use their entries as $k$ new coordinates for each item, and run k-means on those coordinates (von Luxburg, 2007, §4, who gives the unnormalized algorithm and the normalized ones of Shi and Malik and of Ng, Jordan and Weiss). In the ideal case of $k$ groups with no edges between them, eigenvalue 0 has multiplicity $k$ and, for the unnormalized Laplacian and the random-walk one, its eigenspace is spanned by the groups’ indicator vectors (von Luxburg, Propositions 2 and 4), so all items of a group map to one point, whatever the group’s shape in the original space. In practice the result depends on the similarity graph, its connectivity parameter and the number of clusters, the practical choices von Luxburg discusses in §8.

The graph-based segmentation of Felzenszwalb and Huttenlocher (2004) works directly on the pixel grid. It processes edges in order of increasing weight and merges two regions unless the edge between them is heavier than the largest edge inside either region’s minimum spanning tree plus a size-dependent tolerance. The authors call the algorithm closely related to Kruskal’s minimum spanning tree algorithm; in our reading it is single linkage with an adaptive stopping rule. The superpixels and image segmentation articles describe its behaviour.

Evaluating a Clustering

Evaluation is harder than for classification, since cluster indices are arbitrary and there may be no reference answer. External indices compare a clustering with a reference partition, such as the true regions of a synthetic image or human annotations. Internal indices score a clustering from the data alone.

External indices

Every external index is computed from the contingency table between the reference classes $U$ and the clusters $V$, so it ignores how clusters are numbered.

Purity counts, for each cluster, the items of its most common class, and divides by $n$. It is easy to read but rewards fragmentation: $n$ singleton clusters have purity 1. Vinh, Epps and Bailey (2010, §2) note two long-known problems of such set-matching scores: the two partitions may have different numbers of clusters, and the unmatched part of each matched pair is ignored.

Rand index and adjusted Rand index. Over all $\binom{n}{2}$ pairs of items, let $N_{11}$ count the pairs placed together in both partitions and $N_{00}$ the pairs separated in both. The Rand index $\mathrm{RI} = (N_{00} + N_{11}) / \binom{n}{2}$ is the fraction of pairs on which the partitions agree; in practice it often lies between 0.5 and 1, and its value for unrelated partitions is high and not constant (Vinh et al., 2010, §2). Hubert and Arabie (1985) corrected it for chance with the general form (as Vinh et al., 2010, §4.1, give it)

$$\text{adjusted index} = \frac{\text{index} - \mathbb{E}[\text{index}]}{\max(\text{index}) - \mathbb{E}[\text{index}]},$$

where the expectation is taken over random partitions with the same cluster sizes. The adjusted Rand index (ARI) is 1 for identical partitions and close to 0 for random ones, and the scikit-learn User Guide gives its range as $[-0.5, 1]$.

Mutual information. With $P(i) = |U_i|/n$, $P'(j) = |V_j|/n$ and $P(i, j) = n_{ij}/n$, the mutual information $I(U, V) = \sum_{i,j} P(i, j) \log \frac{P(i, j)}{P(i) P'(j)}$ measures how much knowing one partition tells about the other. Dividing by a mean of the entropies $H(U)$ and $H(V)$ gives a normalized mutual information (NMI) in $[0, 1]$; Vinh et al. (2010, Table 2) list variants normalized by the minimum, geometric mean, arithmetic mean, maximum or joint entropy. NMI is not corrected for chance: for random partitions it is near zero when there are many items per cluster but grows as that number falls (Vinh et al., 2010, §4). Using the expected mutual information of two random partitions with the given cluster sizes, Vinh et al. define the adjusted mutual information

$$\mathrm{AMI}(U, V) = \frac{I(U, V) - \mathbb{E}[I(U, V)]}{\operatorname{mean}\big(H(U), H(V)\big) - \mathbb{E}[I(U, V)]} .$$

In their experiments the expected mutual information was close to zero once the data held at least about 100 items per cluster, which they suggest as the regime where the correction matters little (§4.2). In scikit-learn 1.7.2, normalized_mutual_info_score and adjusted_mutual_info_score both default to average_method="arithmetic", the variant Vinh et al. call “sum”.

ARI, NMI and AMI are symmetric in their two arguments, and ARI and AMI are close to zero for random labelings whatever the number of clusters (scikit-learn User Guide §2.3.11). Because they count pairs or information over all items, large clusters dominate them.

Internal indices

Without a reference, a clustering can only be scored by how compact and separated its clusters are under the chosen distance. The scikit-learn User Guide (§2.3.11.5–§2.3.11.7) defines three common indices.

• Silhouette (Rousseeuw, 1987, as the guide cites it). For item $i$, let $a(i)$ be its mean distance to the other members of its cluster and $b(i)$ its mean distance to the members of the nearest other cluster; then $s(i) = \big(b(i) - a(i)\big) / \max\big(a(i), b(i)\big) \in [-1, 1]$, and the score is the mean over items; higher is better, and values near 0 indicate overlapping clusters.
• Calinski–Harabasz index: the ratio of the between-cluster to the within-cluster dispersion, $\frac{\operatorname{tr}(B_k)}{\operatorname{tr}(W_k)} \cdot \frac{n - k}{k - 1}$, with $W_k$ and $B_k$ the within- and between-cluster scatter matrices; higher is better.
• Davies–Bouldin index: for each cluster, the largest over other clusters of $(s_i + s_j)/d_{ij}$, where $s_i$ is the mean distance of a cluster’s items to its centroid and $d_{ij}$ the distance between centroids, averaged over clusters; lower is better, and the centroids tie it to Euclidean space.

All three reward compact, well-separated clusters: the guide notes that the silhouette and the Calinski–Harabasz index are generally higher for convex clusters than for density-based ones such as DBSCAN’s. They are therefore best used to compare clusterings of the same kind, for example k-means with different $k$, and not to decide whether a convex or a density-based clustering is right.

Choosing the Number of Clusters

Every method needs an answer to “how many?”, either directly or through a scale parameter, and no criterion settles it in general. The methods’ own articles cover the standard tools:

• k-means: the elbow of the within-cluster sum of squares, the mean silhouette and the gap statistic (k-means).
• Gaussian mixtures: penalized likelihood, BIC and AIC, or held-out likelihood (Gaussian mixture model).
• Mean shift: the bandwidth, chosen by a plug-in rule, by stability over a range of bandwidths or by the task (mean shift).
• Hierarchical clustering: the height at which to cut the dendrogram; groups that merge at heights much greater than their subgroups are candidates for natural clusters, possibly at several levels, and the gap statistic can choose the level (Hastie et al., §14.3.12).
• Spectral clustering: the eigengap heuristic, choosing $k$ so that the first $k$ Laplacian eigenvalues are small and the next one clearly larger; von Luxburg (2007, §8.3) notes that it works when clusters are well pronounced and gives ambiguous answers when they are not.
• DBSCAN: the knee of the sorted $k$-nearest-neighbour distances fixes $\varepsilon$, and the number of clusters follows.

In vision the number is often fixed by the use rather than discovered: a 16-colour palette, a 1,000-word vocabulary, a budget of superpixels, five anchor shapes.

Scalability

k-means costs $O(nkd)$ per iteration and has mini-batch variants for very large $n$; mean shift with every point as a seed is quadratic in $n$ per iteration without a spatial index; a full pairwise dissimilarity matrix has $n(n-1)/2$ entries, which is what makes naive agglomerative clustering and dense spectral clustering expensive, and spectral methods must also compute eigenvectors of an $n \times n$ matrix. The scikit-learn overview table (User Guide §2.3.1) classes k-means (with its mini-batch code) and DBSCAN as suited to very large numbers of samples, Ward and agglomerative clustering to large ones (agglomerative clustering becomes expensive without connectivity constraints, since each step considers all possible merges), spectral clustering to medium numbers of samples and few clusters, and mean shift and Gaussian mixtures as not scalable. Sparse graphs help: connectivity constraints that allow only neighbouring pixels to merge make agglomerative clustering faster (User Guide §2.3.6), and Felzenszwalb and Huttenlocher’s algorithm runs in $O(m \log m)$ time for $m$ edges, which on a pixel grid grows in proportion to the number of pixels.

A second distinction is whether a method can label new data. k-means, mixtures and mean shift define a rule (nearest centre, largest posterior, nearest mode) that applies to new points; scikit-learn calls them inductive. Agglomerative, spectral, DBSCAN and HDBSCAN clustering label only the data they were run on (transductive).

Examples

Shapes that defeat different families

The example renders a $160 \times 160$ stippled image with six shapes: a disk inside a ring, two long parallel bars 4 pixels apart, a small dense blob and a large sparse blob. Each pixel inside a shape is lit with probability 0.35, except in the sparse blob, where the probability is 0.06. The items to cluster are the $(x, y)$ positions of the 1,145 lit pixels, the kind of grouping a bottom-up method faces when it must assemble detected points into objects. Because the shape of every pixel is known, each clustering can be scored against the truth. The methods with a fixed number of clusters are given the true $k = 6$; mean shift (bandwidth from scikit-learn’s estimate_bandwidth), DBSCAN at three radii and HDBSCAN choose their own. The code uses scikit-learn 1.7.2, single-threaded (OMP_NUM_THREADS=1), and also scores the true partition with the internal indices.

import warnings
import numpy as np
import sklearn
from sklearn.cluster import (KMeans, MeanShift, estimate_bandwidth, AgglomerativeClustering,
                             DBSCAN, HDBSCAN, SpectralClustering)
from sklearn.mixture import GaussianMixture
from sklearn.metrics import (adjusted_rand_score, adjusted_mutual_info_score,
                             silhouette_score, davies_bouldin_score)
from sklearn.metrics.cluster import contingency_matrix
from sklearn.neighbors import kneighbors_graph
from scipy.sparse.csgraph import connected_components
warnings.filterwarnings("ignore", message="Graph is not fully connected")  # see the text

def render(seed, H=160, W=160):
    """A 160 x 160 stippled image: each pixel inside a shape is lit with probability
    0.35 (0.06 inside the large blob). Returns the (x, y) positions of the lit pixels,
    one per row, and the shape each one belongs to."""
    rng = np.random.default_rng(seed)
    yy, xx = np.mgrid[0:H, 0:W]
    shape = np.full((H, W), -1)
    r = np.hypot(yy - 42, xx - 42)
    shape[r < 10] = 0                                             # disk
    shape[(r >= 20) & (r < 26)] = 1                               # ring around the disk
    shape[(yy >= 100) & (yy < 106) & (xx >= 12) & (xx < 148)] = 2  # long bar
    shape[(yy >= 110) & (yy < 116) & (xx >= 12) & (xx < 148)] = 3  # parallel bar, 4 px gap
    shape[np.hypot(yy - 30, xx - 88) < 10] = 4                    # small dense blob
    shape[np.hypot(yy - 52, xx - 126) < 22] = 5                   # large sparse blob
    lit = (shape >= 0) & (rng.random((H, W)) < np.where(shape == 5, 0.06, 0.35))
    return np.stack([xx[lit], yy[lit]], 1).astype(float), shape[lit]

def methods(seed):
    return {
        "k-means (k=6)":       KMeans(6, random_state=seed),
        "GMM, full cov. (6)":  GaussianMixture(6, n_init=5, random_state=seed),
        "mean shift":          "meanshift",
        "single linkage (6)":  AgglomerativeClustering(6, linkage="single"),
        "Ward linkage (6)":    AgglomerativeClustering(6, linkage="ward"),
        "DBSCAN eps=2.5":      DBSCAN(eps=2.5, min_samples=5),
        "DBSCAN eps=4":        DBSCAN(eps=4, min_samples=5),
        "DBSCAN eps=6":        DBSCAN(eps=6, min_samples=5),
        "HDBSCAN":             HDBSCAN(min_cluster_size=20),
        "spectral, 10-NN (6)": SpectralClustering(6, affinity="nearest_neighbors",
                                                  n_neighbors=10, random_state=seed),
    }

def fit(m, X):
    if m == "meanshift":                       # bandwidth from scikit-learn's estimator
        m = MeanShift(bandwidth=estimate_bandwidth(X, random_state=0))
    return m.fit_predict(X)

def purity(y, labels):
    """Fraction of points that belong to the majority shape of their cluster
    (DBSCAN's noise label -1 counts as one cluster, as in the scikit-learn metrics)."""
    return contingency_matrix(y, labels).max(0).sum() / len(y)

X, y = render(0)
print(f"scikit-learn {sklearn.__version__}; {len(X)} lit pixels, per shape {np.bincount(y)}")
print(f"{'':20s} clusters noise   ARI    AMI  purity  silhouette  Davies-Bouldin")
rows = [("true shapes", y)] + [(n, fit(m, X)) for n, m in methods(0).items()]
for name, lab in rows:
    k = len(set(lab) - {-1})
    print(f"{name:20s} {k:5d} {np.sum(lab == -1):6d}  {adjusted_rand_score(y, lab):5.3f}  "
          f"{adjusted_mutual_info_score(y, lab):5.3f}  {purity(y, lab):5.3f}  "
          f"{silhouette_score(X, lab):8.3f}  {davies_bouldin_score(X, lab):10.3f}")

# Which shapes share their main cluster, and what fraction of each shape that cluster holds.
print("\nshare of each shape in its main cluster (noise counts as no cluster)")
print("shape:                 0     1     2     3     4     5    shapes sharing a main cluster")
for name, lab in rows[1:]:
    M = contingency_matrix(y, lab)              # rows: shapes 0-5, columns: np.unique(lab)
    M = M[:, np.unique(lab) != -1]              # noise points belong to no cluster
    main = M.argmax(1)
    shared = [g for g in ("+".join(str(s) for s in np.flatnonzero(main == c)) for c in set(main))
              if "+" in g]
    print(f"{name:20s}" + "".join(f"{M[s, main[s]] / np.sum(y == s):6.2f}" for s in range(6))
          + "    " + (", ".join(shared) or "none"))

# The 10-nearest-neighbour graph that spectral clustering builds: its connected components.
A = kneighbors_graph(X, 10, include_self=True)
n_comp, comp = connected_components(A + A.T, directed=False)
print(f"\n10-NN graph: {n_comp} connected components with sizes {np.bincount(comp)}")

# Cluster ids are arbitrary: compare k-means ids with shape ids directly, then renamed.
km = rows[1][1]
renamed = np.random.default_rng(1).permutation(6)[km]
print(f"k-means: label == shape id for {np.mean(km == y):.3f} of pixels, "
      f"{np.mean(renamed == y):.3f} after renaming the clusters; "
      f"ARI {adjusted_rand_score(y, km):.3f} and {adjusted_rand_score(y, renamed):.3f}")

Output:

scikit-learn 1.7.2; 1145 lit pixels, per shape [109 300 264 266 117  89]
                     clusters noise   ARI    AMI  purity  silhouette  Davies-Bouldin
true shapes              6      0  1.000  1.000  1.000     0.141      11.026
k-means (k=6)            6      0  0.403  0.620  0.602     0.507       0.638
GMM, full cov. (6)       6      0  0.523  0.689  0.678     0.517       0.548
mean shift               3      0  0.346  0.511  0.494     0.513       0.618
single linkage (6)       6      0  0.721  0.893  0.769     0.217       8.987
Ward linkage (6)         6      0  0.392  0.618  0.591     0.468       0.678
DBSCAN eps=2.5          24    119  0.470  0.702  0.969     0.193       0.974
DBSCAN eps=4            12     42  0.986  0.962  1.000     0.064       5.937
DBSCAN eps=6             7      3  0.715  0.878  0.769     0.193       7.039
HDBSCAN                  5     15  0.716  0.884  0.769     0.214       9.306
spectral, 10-NN (6)      6      0  0.828  0.906  0.905     0.146       3.868

share of each shape in its main cluster (noise counts as no cluster)
shape:                 0     1     2     3     4     5    shapes sharing a main cluster
k-means (k=6)         1.00  0.70  0.36  0.37  1.00  1.00    0+1, 2+3
GMM, full cov. (6)    1.00  1.00  0.36  0.37  1.00  1.00    0+1, 2+3
mean shift            1.00  1.00  0.54  0.54  1.00  0.88    0+1+4, 2+3+5
single linkage (6)    1.00  1.00  1.00  1.00  1.00  0.99    2+3
Ward linkage (6)      1.00  0.60  0.51  0.42  1.00  1.00    0+1, 2+3
DBSCAN eps=2.5        0.99  0.58  0.33  0.41  1.00  0.06    none
DBSCAN eps=4          1.00  1.00  1.00  1.00  1.00  0.11    none
DBSCAN eps=6          1.00  1.00  1.00  1.00  1.00  0.81    2+3
HDBSCAN               1.00  1.00  1.00  1.00  1.00  0.83    2+3
spectral, 10-NN (6)   1.00  0.54  1.00  1.00  1.00  1.00    0+4

10-NN graph: 5 connected components with sizes [300 117  89 109 530]
k-means: label == shape id for 0.086 of pixels, 0.000 after renaming the clusters; ARI 0.403 and 0.403

The methods that build convex or elliptical clusters fail in the same places. k-means and Ward linkage put the disk and part of the ring into one cluster, and they cut the bars crosswise into segments, with the main segment of each bar shared with the other bar; by our reasoning, a compact group straddling the 4-pixel gap has a smaller within-cluster sum of squares than a bar 136 pixels long. The full-covariance mixture could fit each bar with one elongated ellipse, but having spent one component on the disk and ring together, it also slices the bars crosswise. Mean shift’s estimated bandwidth leaves only three modes: one for the disk, the ring and the small blob, and two that split the bars between them, one of which also takes most of the sparse blob. Single linkage keeps the disk, the ring and the small blob whole, and the sparse blob nearly so, but chains the two bars together. By the minimum-spanning-tree reading above, it cut the tree’s five longest edges; the shortest link between the bars, about 5 pixels, was not among them, while a link to a stray point of the sparse blob was, so the sixth cluster is a single stray point.

DBSCAN shows the single-scale problem that Ester et al. describe. With $\varepsilon = 2.5$ the ring and the bars break into pieces and no cluster holds more than 6% of the sparse blob; with $\varepsilon = 4$ every dense shape is recovered whole, but the sparse blob falls apart: 11% of it is in its largest cluster, and a purity of 1.000, which counts noise as one cluster, shows that all 42 noise points come from it; with $\varepsilon = 6$ the sparse blob becomes one cluster holding 81% of its points, and the bars merge. None of the three radii both separates the bars and keeps the sparse blob. HDBSCAN, which searches across densities, keeps the sparse blob but also merges the bars at this minimum cluster size. Spectral clustering separates the bars and the ring from the disk, but it splits the ring and puts the disk and the small blob into one cluster. Its 10-nearest-neighbour graph has five connected components (the two bars form one, of 530 points), so the Laplacian has five zero eigenvalues, one per component (von Luxburg, 2007, Proposition 4); von Luxburg (§8.1) advises making the similarity graph connected, or nearly so, before clustering with it. scikit-learn warns that the graph is not fully connected, and the code silences the warning.

The internal and external indices disagree completely. The true partition has the second-lowest silhouette (0.141) and the worst Davies–Bouldin index, while k-means, the mixture, mean shift and Ward linkage, which have four of the five lowest ARIs, have the four best scores on both internal indices. These indices measure compactness and separation in Euclidean distance, and a ring around a disk, or two bars 4 pixels apart, is neither compact nor well separated in that sense. Purity misleads in the other direction: DBSCAN at $\varepsilon = 2.5$ scores 0.969 by fragmenting the shapes into 24 clusters plus noise, while its ARI is 0.470. ARI and AMI rank the methods much as the per-shape shares do, but both barely penalize the loss of the sparse blob at $\varepsilon = 4$ (ARI 0.986), since its 89 of 1,145 points carry little weight. Finally, the k-means labels equal the shape numbers for 8.6% of pixels, and for none after the clusters are renamed, while ARI is unchanged: agreement of raw labels says nothing about a clustering.

The same comparison over ten images

The second block draws ten stipple patterns with seeds 0 to 9 (seed 0 is the image above), using the same seed for the algorithms, and records each method’s ARI and silhouette, and how often the method’s silhouette exceeds that of the true shapes.

import numpy as np
from sklearn.metrics import adjusted_rand_score, silhouette_score
# Uses render, methods and fit from the first block.

ari, sil, beats = {}, {}, {}
for seed in range(10):                          # new stipple pattern and algorithm seed
    X, y = render(seed)
    s_true = silhouette_score(X, y)
    sil.setdefault("true shapes", []).append(s_true)
    for name, m in methods(seed).items():
        lab = fit(m, X)
        ari.setdefault(name, []).append(adjusted_rand_score(y, lab))
        sil.setdefault(name, []).append(silhouette_score(X, lab))
        beats[name] = beats.get(name, 0) + (sil[name][-1] > s_true)

print("seeds 0-9           ARI min  median  max   silhouette min-max   silhouette above truth")
print(f"{'true shapes':20s} {1:5.3f}  {1:6.3f}  {1:5.3f}     "
      f"{min(sil['true shapes']):6.3f} to {max(sil['true shapes']):5.3f}")
for name in ari:
    a = np.array(ari[name])
    print(f"{name:20s} {a.min():5.3f}  {np.median(a):6.3f}  {a.max():5.3f}     "
          f"{min(sil[name]):6.3f} to {max(sil[name]):5.3f}      {beats[name]:2d} of 10")

Output:

seeds 0-9           ARI min  median  max   silhouette min-max   silhouette above truth
true shapes          1.000   1.000  1.000      0.124 to 0.143
k-means (k=6)        0.356   0.403  0.444      0.448 to 0.513      10 of 10
GMM, full cov. (6)   0.446   0.502  0.523      0.470 to 0.524      10 of 10
mean shift           0.330   0.346  0.408      0.467 to 0.516      10 of 10
single linkage (6)   0.681   0.702  0.727      0.124 to 0.218      10 of 10
Ward linkage (6)     0.392   0.507  0.542      0.468 to 0.509      10 of 10
DBSCAN eps=2.5       0.368   0.448  0.498      0.121 to 0.288       9 of 10
DBSCAN eps=4         0.929   0.986  0.991      0.046 to 0.082       0 of 10
DBSCAN eps=6         0.677   0.698  0.720      0.156 to 0.218      10 of 10
HDBSCAN              0.676   0.698  0.724      0.192 to 0.224      10 of 10
spectral, 10-NN (6)  0.334   0.614  0.868     -0.057 to 0.275       3 of 10

The ranking is stable across images. DBSCAN with $\varepsilon = 4$ has the best ARI on every image (its lowest, 0.929, is above every other method’s highest) and is the only method whose silhouette never exceeds the truth’s; that radius, though, was picked by trying several against the known answer, which an unlabelled problem does not allow. The convex and elliptical methods stay between 0.33 and 0.54 in ARI, while their silhouettes, from 0.45 to 0.52, are more than three times the truth’s and above every other method’s highest, so choosing by silhouette would have picked one of them on every image; only a reference partition, or a clear idea of what shape a cluster should have, can make that choice. Spectral clustering is the least stable, with ARI from 0.334 to 0.868.

Common Misconceptions

• “The data have a true clustering that a good algorithm finds.” The data and a similarity measure together admit many defensible groupings; Kleinberg’s theorem shows that any method must give up scale invariance, richness or consistency. Decide what a cluster should look like for the task, then choose the method.
• “A higher silhouette means a better clustering.” It means more compact, better separated clusters in the given distance; in the example the true shapes scored below every convex method.
• “k-means is the default clustering method.” It is the default for compact, similarly sized, convex groups and for quantization, where any partition into cells is acceptable. For connected, curved or nested structures, density-based, single-linkage or graph methods fit better, and for groups of very different density none of these may be enough.
• “Cluster label 3 means the same thing in two runs.” Labels are arbitrary; compare partitions through the contingency table, with ARI or AMI, or match clusters explicitly (for example with the Hungarian algorithm) before comparing them label by label.
• “One external score is enough.” Purity rises as clusters fragment, and ARI and AMI weigh clusters by size, so a lost small cluster barely registers; look at the contingency table as well.
• “The dendrogram reveals the structure of the data.” It shows the structure imposed by the chosen linkage, a hierarchy whether or not the data have one, and different linkages or small changes in the data can change it considerably (Hastie et al., §14.3.12).

Where It Is Used

• Colour quantization and colour segmentation. k-means on pixel colours gives a $k$-colour palette, and clustering colour, or colour and position, is the simplest image segmentation; mean shift in the joint spatial–colour domain gives edge-preserving filtering and segmentation (k-means, mean shift).
• Superpixels. SLIC runs a localized k-means in colour and position, quick shift and mean shift seek modes, and Felzenszwalb–Huttenlocher merges regions on the pixel graph (superpixels).
• Graph-based segmentation. Normalized cuts and spectral clustering partition a pixel-affinity graph; graph cuts find a minimum-energy labelling of the pixel graph, for example from user seeds in interactive segmentation.
• Background and colour models. Per-pixel Gaussian mixtures model the background of a video, and colour mixtures for foreground and background drive interactive segmentation (Gaussian mixture model).
• Visual vocabularies and descriptor indexing. k-means codebooks quantize local descriptors into visual words (bag of visual words), and inverted-file indexes and product quantization use k-means centroids to search large collections (approximate nearest neighbor search).
• Anchor boxes. YOLOv2 chose its anchor boxes by k-means on box widths and heights with the distance $1 - \mathrm{IoU}$ instead of the Euclidean one.
• Deep clustering and self-supervised learning. DeepCluster alternates k-means on a network’s features with training the network to predict the cluster indices, and SwAV predicts balanced cluster assignments of one view from another (self-supervised learning).
• Clustering embeddings. Bottom-up instance segmentation maps pixels to embeddings and groups them with a variant of mean shift. FaceNet’s paper shows a cluster of one person’s faces from a personal photo collection, produced by agglomerative clustering of its face embeddings, without a quantitative evaluation (FaceNet); face clustering is one of the tasks of face recognition benchmarks such as IJB-B and IJB-C. Clustering image features also finds near-duplicate videos when datasets are built, as for Kinetics.