Clustering: What Each Algorithm Assumes a Cluster Is

The question nobody asks out loud

Supervised learning has an unambiguous scoreboard: there are labels, and you either predict them or you do not. Clustering has no such thing. You hand an algorithm a pile of unlabelled points, it hands back a partition, and no external fact says whether the partition is right.

That is not because clustering is vague. It is because the word “cluster” is doing an enormous amount of unstated work. Consider three perfectly reasonable definitions:

1. A cluster is a group of points that are all close to some common centre.
2. A cluster is a group of points you can walk between in small steps without leaving the group.
3. A cluster is a region where the data-generating density has a local maximum.

These sound like paraphrases. They are not. On two concentric rings, definition 1 says the two rings are one cluster (they share a centre); definitions 2 and 3 say they are two. On a long thin filament that nearly touches a dense blob, definition 2 fuses them and definition 1 does not.

Every clustering algorithm hard-codes one of these definitions. The algorithm’s failure cases are not bugs; they are the definition being applied faithfully to data that disagrees with it.

k-means: a cluster is a ball around a centroid

The objective

Fix the number of clusters \(k\). A partition assigns each of the \(n\) points \(x_1, \dots, x_n \in \mathbb{R}^d\) to one of \(k\) sets \(C_1, \dots, C_k\), and each set gets a representative \(\mu_j \in \mathbb{R}^d\). k-means minimises the within-cluster sum of squares:

Two things are already decided by this line. First, the distance is Euclidean and squared, so a point twice as far away is penalised four times as much: k-means is dragged around by outliers. Second, each cluster is summarised by a single point. Whatever a cluster is, k-means has committed to it being the things near one location.

Lloyd’s algorithm as alternating minimisation

\(J\) depends on two arguments, and it is hard to optimise jointly — the exact problem is NP-hard for general \(k\) and \(d\). But each argument is trivial given the other, which is exactly the setting for alternating minimisation. Lloyd’s algorithm just alternates the two easy problems:

Assignment step — fix the centroids, optimise the labels. With \(\mu\) held fixed, \(J\) decouples over points, so every point independently picks its nearest centroid:

Update step — fix the labels, optimise the centroids. With \(C\) held fixed, \(J\) decouples over clusters, and \(\sum_{x \in C_j} \lVert x - \mu_j \rVert^2\) is a convex quadratic in \(\mu_j\) minimised at the mean:

That is the whole algorithm. The “means” in k-means is not a design choice; it is what falls out of minimising squared Euclidean distance. Swap the loss for absolute distance and the optimal representative becomes the componentwise median, which is a different algorithm (k-medians) with different failure modes.

A worked iteration, by hand

Take six points in the plane:

Set \(k = 2\) and initialise the centroids at two of the data points, \(\mu_1 = (1,1)\) and \(\mu_2 = (9,2)\). The assignment step needs squared distances, and because we only compare them we never take a square root:

So \(C_1 = \{p_1, p_2, p_3\}\) and \(C_2 = \{p_4, p_5, p_6\}\), and the objective at this moment — still using the old centroids — is \(0 + 5 + 5 + 5 + 2 + 0 = 17\).

Now the update step. Averaging coordinates:

Recomputing the objective with the new centroids, cluster \(C_1\) contributes \((1{+}1) + (0{+}1) + (1{+}0) = 4\) and \(C_2\) contributes the same by symmetry, so \(J = 8\). One iteration took the objective from \(17\) to \(8\). Running the assignment step again returns the identical labels, so the algorithm has converged after a single full sweep.

It converges, but that is a weaker promise than it sounds

Both steps are exact minimisations of \(J\) over their own argument, so \(J\) never increases. There are finitely many partitions of \(n\) points into \(k\) groups, and each distinct partition that the algorithm visits has a strictly lower cost than the last, so no partition can repeat and the algorithm must halt in finitely many steps.

What it halts at is a local optimum, and “local” can be arbitrarily bad. Take four points at the corners of a \(10 \times 1\) rectangle: \((0,0), (0,1), (10,0), (10,1)\), with \(k = 2\). Splitting by the long side — pairing \((0,0)\) with \((0,1)\) and \((10,0)\) with \((10,1)\) — gives centroids \((0, 0.5)\) and \((10, 0.5)\) and cost \(J = 1\). Splitting by the short side gives centroids \((5, 0)\) and \((5, 1)\) and cost \(J = 100\). And that second, hundredfold-worse partition is a genuine fixed point: the point \((0,0)\) sits at squared distance \(25\) from \((5,0)\) and \(26\) from \((5,1)\), so it refuses to move. Lloyd’s algorithm, started there, stops there forever.

k-means++: spend a little to initialise well

Because the destination depends entirely on the start, initialisation matters more than any other detail. k-means++ picks the first centre uniformly at random, then picks each subsequent centre from the data with probability proportional to the squared distance to the nearest centre already chosen:

The effect is to spread the initial centres out — far-away points are exponentially more likely to be picked than near ones, because of the squaring — while remaining random enough not to latch onto a lone outlier every time. Arthur and Vassilvitskii proved that this seeding alone, before a single Lloyd iteration, gives an expected cost within a factor of \(O(\log k)\) of the optimum. In practice you still run it several times with different seeds and keep the lowest \(J\); that is not a workaround, it is the standard recipe.

What k-means assumes, and what breaks

The assignment rule compares \(\lVert x - \mu_j \rVert\) across \(j\), so the boundary between two clusters is the perpendicular bisector of their centroids. Every k-means cluster is therefore an intersection of half-spaces: a convex Voronoi cell. That single fact generates the whole list of failure cases.

• Clusters must be roughly spherical. An elongated, diagonal cloud has points at its two ends that are further from each other than they are from a neighbouring cluster’s centroid, so k-means slices it crosswise.
• Clusters must be roughly the same size and spread. The objective sums squared distances over all points, so a large, diffuse cluster contributes far more cost than a small tight one. k-means buys cost reductions by carving up the big cluster and donating its edges to the small one.
• Clusters must be convex. Two moons, two concentric rings, an S-shape: none of these is a Voronoi cell, so none of them is reachable.
• \(k\) must be supplied. k-means will happily return seventeen clusters from data with three.
• Every point gets a label. There is no “none of the above”. A single distant outlier either drags a centroid towards it or forms a singleton cluster, and either way it costs you.

The concentric-rings case is the cleanest demonstration, because it fails for a reason you can state in one line rather than observe empirically.

Why is k-means guaranteed to lose here? Because both rings are centred on the origin, their means are essentially the same point — in the data above they are \(0.056\) apart, against ring radii of \(2\) and \(5\). The correct partition is not merely a bad local optimum of \(J\); it is not a fixed point of Lloyd’s algorithm at all, because with two coincident centroids the assignment step cannot reproduce it. No amount of restarting, no cleverer initialisation, and no better optimiser will find it. The right answer is outside the model’s vocabulary.

Hierarchical clustering: a cluster is whatever the linkage rule fuses

Agglomerative merging

Agglomerative clustering starts with \(n\) singleton clusters and repeatedly merges the two closest, until one cluster remains. The only thing you have to specify is what “closest” means for two sets of points, and that choice — the linkage — is the entire algorithm.

Let \(A\) and \(B\) be clusters and \(d(a,b)\) the distance between points.

These are four different beliefs about what makes a cluster coherent. Single linkage believes in connectivity: two groups belong together if any two of their members are close, so it can trace out long, thin, curved shapes. Complete linkage believes in compactness: a merge is only allowed if it keeps the worst-case diameter small, so it produces tight, round groups and refuses to string things out. Average linkage splits the difference. Ward is the odd one out — its criterion is exactly the increase in within-cluster sum of squares caused by the merge, which makes it a greedy hierarchical version of the k-means objective. That is why Ward’s dendrograms look so much like k-means solutions, and why they inherit the same spherical bias.

The output is a dendrogram, not a partition

This is the genuinely different thing about hierarchical clustering, and it is easy to miss because the library call returns labels. The algorithm’s actual product is the dendrogram: a binary tree recording every merge and the distance at which it happened. Cutting the tree at a height gives a partition; cutting at a different height gives a different, strictly nested one. You get the family of all of them at once.

That is a real advantage when the structure genuinely is nested — taxonomies, document topics, gene expression — because you can choose the level that matters for your question after seeing the hierarchy rather than before. You never have to commit to \(k\) up front. You do still have to commit eventually, and the dendrogram will not tell you where to cut; it only lays out the options with their merge heights, so a big vertical gap between successive merges hints that the structure below it is more stable than the structure above.

What hierarchical clustering does not buy you

Three more limitations worth stating plainly. Merges are greedy and irreversible: a bad early merge, made when the algorithm has seen almost no structure, is never revisited. The naive implementation needs the full pairwise distance matrix, so memory is \(O(n^2)\) and time is \(O(n^3)\) (or \(O(n^2 \log n)\) with a heap; single linkage can be done in \(O(n^2)\)), which puts hierarchical clustering out of reach for large \(n\) without approximation. And there is no notion of noise — every point ends up somewhere in the tree, so outliers appear as singleton branches that merge in very late.

DBSCAN: a cluster is a connected dense region

Density instead of distance-to-centre

DBSCAN throws away the centroid entirely. It asks a local question of each point: are there enough other points nearby? Two parameters set the meaning of “nearby” and “enough”:

• \(\varepsilon\) — the radius of the neighbourhood.
• \(\text{minPts}\) — how many points must lie inside that radius (counting the point itself) for the neighbourhood to count as dense.

A point with at least \(\text{minPts}\) neighbours within \(\varepsilon\) is a core point. A point that is not core but lies within \(\varepsilon\) of a core point is a border point. Everything else is noise.

The clusters are then the connected components of the core points: two core points are in the same cluster if you can hop from one to the other in steps of length at most \(\varepsilon\), always landing on core points. Border points join whichever cluster reached them.

What DBSCAN buys, and what it costs

Two capabilities here are genuinely unavailable from k-means. First, arbitrary shapes: rings, crescents, spirals, anything density-connected. Second, explicit noise: points in sparse regions are labelled as noise rather than being forced into a cluster, which means one distant outlier no longer distorts anything. The number of clusters is also an output rather than an input — you specify a density scale and DBSCAN reports how many components exist at that scale.

The bill comes due in three places.

• Varying density is the fundamental weakness. A single global \(\varepsilon\) sets one density threshold for the whole dataset. If one cluster is tight and another is diffuse, no \(\varepsilon\) works: small values shatter the diffuse cluster into fragments and label most of it noise, large values merge the tight cluster with its neighbours. This is the mirror image of the rings problem — a case where the definition of “cluster” is wrong for the data. HDBSCAN exists precisely to relax this by considering all density scales at once.
• The parameters are not intuitive in more than two or three dimensions, because \(\varepsilon\) is an absolute distance and you rarely have a feel for absolute distances in high-dimensional feature space. The usual diagnostic is a \(k\)-distance plot: sort every point by its distance to its \(\text{minPts}\)-th nearest neighbour, plot the sorted curve, and look for the knee. It is a heuristic, and it is the same kind of eyeballing as the elbow method below.
• High dimensions degrade it, as they degrade every distance-based method. When distances concentrate, the ratio between the nearest and farthest neighbour of a point approaches one, and the notion of an \(\varepsilon\)-ball that is “dense but not everything” evaporates.

Border points also make DBSCAN mildly order-dependent: a border point reachable from two clusters is assigned to whichever is processed first. It is a small effect, but it means the output is not a deterministic function of the data alone.

Gaussian mixtures: a cluster is a bump in a density

The generative story

The methods so far are procedures. A Gaussian mixture model (GMM) is instead a model of how the data were produced: pick a component \(j\) with probability \(\pi_j\), then draw a point from \(\mathcal{N}(\mu_j, \Sigma_j)\). The density is

Because it is a probability model, it answers questions the others cannot: how likely is this new point under the fitted structure, how confident is this assignment, what happens if I generate fresh samples. Instead of a hard label, each point gets a vector of responsibilities — the posterior probability of each component given the point:

A point squarely inside one component gets responsibilities near \((1, 0, \dots)\); a point in the overlap gets something like \((0.55, 0.45, \dots)\) — honest information about ambiguity that a hard partition throws away.

EM is the same alternating trick

Maximising the log-likelihood directly is awkward because of the sum inside the logarithm. Expectation–Maximisation alternates, exactly as Lloyd’s algorithm does:

• E step — compute the responsibilities \(\gamma_{ij}\) from the current parameters, using the formula above.
• M step — refit each component using the responsibilities as soft counts. Writing \(N_j = \sum_i \gamma_{ij}\) for the effective number of points in component \(j\):

Each round is guaranteed not to decrease the log-likelihood, which — as with k-means — buys convergence to a stationary point and nothing about its quality. Multiple restarts are again standard, and k-means is the usual way to initialise a GMM.

How it generalises k-means

Put the three restrictions on a GMM at once: force every covariance to be \(\Sigma_j = \sigma^2 I\) (spherical, equal across components), fix all mixing weights equal, and let \(\sigma^2 \to 0\). As the components become infinitely tight, the responsibility of the nearest component tends to \(1\) and all others to \(0\) — the E step hardens into “assign to nearest centroid”. The M step, with responsibilities that are \(0\) or \(1\), reduces to “average the assigned points”. That is Lloyd’s algorithm exactly. k-means is the hard-assignment, equal-spherical-covariance limit of EM for a Gaussian mixture.

Reading that in the other direction tells you what each relaxation buys:

• Letting \(\Sigma_j\) be full lets a component be an ellipse at any orientation, so elongated diagonal clusters — one of k-means’ clean failures — are handled natively.
• Letting \(\Sigma_j\) differ between components handles clusters of different sizes and spreads.
• Keeping responsibilities soft stops points in the overlap from being arbitrarily forced one way, and gives you a calibrated measure of that ambiguity.
• Having a likelihood gives you a principled model-selection criterion for \(k\), which no other method here has.

What a GMM does not buy you

Beyond that: a GMM still assumes each cluster is Gaussian, so it fails on the two rings just as k-means does — a ring is not an ellipse-shaped bump, and a mixture of two Gaussians centred at the origin cannot represent one. A full covariance matrix costs \(d(d+1)/2\) parameters per component, so in high dimensions you run out of data and are pushed back to diagonal or spherical covariances, which is a return to k-means-like assumptions. And \(k\) must still be chosen, though at least here there is a criterion rather than a squint.

The comparison, in one table

Choosing \(k\), honestly

“Number of clusters” is not a property of the data. It is a property of the data plus a definition of cluster plus a scale of interest. A satellite image has two clusters (land, sea), or nine (terrain types), or thousands (individual fields), and all three are correct answers to different questions. The tools below narrow the range; they do not settle it.

The elbow method. Plot the k-means objective \(J\) against \(k\). It decreases monotonically — more centroids can never fit worse — so you are looking for the kink where the marginal improvement flattens. If the data really has \(k^{\star}\) well-separated spherical clusters, adding the \(k^{\star}\)-th centroid removes a large chunk of cost and the next one removes little, so the kink is real. On data that is merely a continuum of density, the curve is a smooth arc with no kink at all, and you will find yourself squinting at a plot in order to justify a decision you have already made.

Silhouette. For each point \(i\), let \(a(i)\) be its mean distance to the other points in its own cluster and \(b(i)\) the mean distance to the points of the nearest other cluster. Then

and the score for a clustering is the mean of \(s(i)\) over all points. Near \(1\) means points sit much closer to their own cluster than to any other; near \(0\) means they sit on a boundary; negative means they would be better off elsewhere. Unlike the elbow, silhouette does not decrease automatically with \(k\), so its maximum is a defensible choice. But note what \(a\) and \(b\) measure: mean distances, which is a compactness-and-separation criterion. Silhouette is therefore biased towards exactly the round, well-separated clusters that k-means likes, and will score a correct DBSCAN partition of two rings poorly. Per-point silhouettes beat the average — they show which points are badly placed.

BIC, for models with a likelihood. A GMM has a likelihood, so you can penalise complexity honestly. With \(\hat{L}\) the maximised likelihood, \(p\) the number of free parameters and \(n\) the sample size:

and you pick the \(k\) minimising it. The \(p \ln n\) term makes extra components pay for themselves, which is why BIC has an interior minimum where raw likelihood does not. This is the most principled of the three — and it still assumes the model class is right. Fit a GMM to two rings and BIC will confidently select some number of Gaussians that tile the rings; it is answering “how many Gaussians?”, which was never your question.

Which one to reach for

• Large \(n\), features already scaled, you expect compact blobs, and you need speed: k-means with k-means++ and several restarts. It is \(O(nkd)\) per iteration and hard to beat on throughput.
• You suspect nested structure, or \(n\) is small enough to hold a distance matrix, or you want to see the whole hierarchy before committing: agglomerative clustering. Ward if you expect compact groups, single linkage if you expect filaments and are willing to risk chaining.
• Clusters may be oddly shaped, and there is real junk in the data you want excluded: DBSCAN, provided density is broadly uniform within each cluster. If it is not, reach for HDBSCAN.
• You want probabilities, calibrated uncertainty, elliptical clusters, or a criterion for \(k\): a Gaussian mixture. It is the right default when clusters overlap and you need to say how much.

One practical note that outranks all of the above: every method here except the hierarchical ones with a custom metric depends on the scale of your features. A feature measured in metres and one measured in millimetres are not comparable under Euclidean distance, and the second will dominate every distance in the dataset. Standardise first, or decide deliberately not to and be able to say why.