Clustering

Every model so far, from nearest neighbours to the deep networks of the previous note, learned from pairs $(x, y)$: an input and the answer it should produce. Labels are the expensive part of a dataset, and often there are none: a log of customer sessions, a million unannotated images, the vectors an encoder produces. The data still have structure, and the first question to ask of them is whether they fall into groups. This note is about that question. It shows why “group” needs a definition before it has an answer, develops the most used definition, k-means, as an optimization problem with a simple algorithm, and then turns to definitions that k-means cannot express: dense regions, hierarchies and graphs.

Notation. $x_1, \dots, x_n \in \mathbb R^d$ are the data points, $K$ the number of clusters, $c_i \in \{1, \dots, K\}$ the cluster of point $i$, $C_k = \{i : c_i = k\}$ the points of cluster $k$ and $n_k = \lvert C_k \rvert$ their number. $\mu_1, \dots, \mu_K \in \mathbb R^d$ are cluster centres (centroids), $\lVert \cdot \rVert$ the Euclidean norm.

What is a cluster?#

A clustering is a partition of the points into groups such that points in a group are similar and points in different groups are not. That sentence leaves two choices open: what “similar” means, and how the groups are allowed to look. Three common answers:

• Compact groups. Each cluster is a ball around a centre, and a good clustering has small distances from points to their centre. This is k-means.
• Dense regions. A cluster is a region where points are packed closely, of any shape, separated from other clusters by sparse space; isolated points belong to no cluster. This is DBSCAN.
• Connected structure. A cluster is a set of points linked by chains of near neighbours, or a part of a similarity graph with few edges leaving it. This is single-linkage and spectral clustering.

On well-separated round blobs all three agree. On anything else they give different answers, and none of them is wrong: they answer different questions. Kleinberg (2002) made this precise. He asked for a clustering function of the pairwise distances with three properties: it ignores the unit of distance (scale invariance), every partition is a possible output for some distances (richness), and shrinking distances within clusters while stretching distances between them does not change the answer (consistency). His Theorem 2.1 shows that no function has all three, and Theorem 2.2 that single linkage with a suitable stopping rule achieves any two of them. There is no definition-free clustering.

The choice of distance matters as much as the choice of method. All the methods below see the data only through distances, so the units of the features decide the answer: multiplying one feature by 10 makes it dominate every distance. On the four-blob data of the Code section, k-means recovers the blobs exactly (adjusted Rand index 1.0, defined below), and after the second feature is multiplied by 10 only partially (0.736). Scaling and encoding, the subject of the data and features note, are part of the clustering, and in high dimensions the distances themselves become less informative, for the reason given in the nearest neighbours note and in the shell of the multivariate Gaussian.

k-means: the objective#

k-means makes “compact” precise with a number. Given an assignment $c$ and centres $\mu$, the within-cluster sum of squares is

$$J(c, \mu) = \sum_{i=1}^n \lVert x_i - \mu_{c_i} \rVert^2 = \sum_{k=1}^K \sum_{i \in C_k} \lVert x_i - \mu_k \rVert^2 ,$$

and k-means asks for the $c$ and $\mu$ that minimize it. scikit-learn calls $J$ the inertia. Two facts make the objective concrete.

The best centre of a cluster is its mean. For fixed points $C_k$, the function $\mu \mapsto \sum_{i \in C_k} \lVert x_i - \mu \rVert^2$ is a strictly convex quadratic; its gradient $-2\sum_{i \in C_k} (x_i - \mu)$ vanishes at the mean $\bar x_k = \frac{1}{n_k}\sum_{i \in C_k} x_i$. Expanding around the mean shows how much any other centre loses:

$$\sum_{i \in C_k} \lVert x_i - \mu \rVert^2 = \sum_{i \in C_k} \lVert x_i - \bar x_k \rVert^2 + n_k \lVert \bar x_k - \mu \rVert^2 ,$$

because the cross term $2(\bar x_k - \mu)^\top \sum_{i \in C_k}(x_i - \bar x_k)$ is zero. This is the statement of the conditional expectation note that the mean is the best constant prediction under squared error, applied to each cluster.

Small within means large between. Summing the same identity with $\mu = \bar x$, the mean of all the data, over the clusters gives

$$\underbrace{\sum_{i=1}^n \lVert x_i - \bar x \rVert^2}_{\text{total}} = \underbrace{\sum_{k=1}^K \sum_{i \in C_k} \lVert x_i - \bar x_k \rVert^2}_{\text{within}} + \underbrace{\sum_{k=1}^K n_k \lVert \bar x_k - \bar x \rVert^2}_{\text{between}} ,$$

the sample version of the law of total variance. The total is fixed by the data, so minimizing the within-cluster scatter is the same as maximizing the scatter of the cluster means. The within term can also be written without centres, $\sum_{i \in C_k} \lVert x_i - \bar x_k \rVert^2 = \frac{1}{2 n_k} \sum_{i, j \in C_k} \lVert x_i - x_j \rVert^2$: k-means penalizes the average squared distance between points of the same cluster.

Minimizing $J$ exactly is NP-hard, even for two clusters when the dimension is part of the input (Aloise et al., 2009). What is used in practice is a local search.

Lloyd’s algorithm is coordinate descent#

$J$ depends on two blocks of variables, the discrete assignment $c$ and the continuous centres $\mu$. Each block is easy to optimize when the other is fixed, and Lloyd’s algorithm (Lloyd, 1982; the name k-means comes from MacQueen, 1967) alternates the two exact minimizations:

1. Assignment step. With $\mu$ fixed, $J$ is a sum of separate terms, one per point, so each point independently takes its nearest centre: $c_i \leftarrow \arg\min_k \lVert x_i - \mu_k \rVert^2$ (on a tie, a point keeps its current cluster).
2. Update step. With $c$ fixed, $J$ is a sum of separate terms, one per cluster, and each is minimized by the cluster mean: $\mu_k \leftarrow \bar x_k$. (A cluster that loses all its points keeps its old centre or is restarted at a far point.)

Each step minimizes $J$ over one block, so neither can increase it:

$$J(c^{t+1}, \mu^{t+1}) \le J(c^{t+1}, \mu^{t}) \le J(c^{t}, \mu^{t}) .$$

This is coordinate descent with two coordinates. It also terminates. There are at most $K^n$ assignments; when an assignment step changes $c$, it strictly lowers $J$ (a point moves only to a strictly nearer centre), so no assignment can repeat, and the algorithm stops after finitely many steps at a point that neither step can improve. That point is a local minimum in this coordinate-wise sense, not the global one. Each iteration costs $O(nKd)$ for the distances, and in practice few iterations are needed.

Worked example. Five points on a line, $0, 2, 3, 8, 10$, with $K = 2$ and starting centres $\mu = (0, 3)$.

iteration	centres	clusters	$J$
0	$0,\ 3$	$\{0\}$, $\{2, 3, 8, 10\}$	$0 + 1 + 0 + 25 + 49 = 75$
1	$0,\ 5.75$	$\{0, 2\}$, $\{3, 8, 10\}$	$34.69$
2	$1,\ 7$	$\{0, 2, 3\}$, $\{8, 10\}$	$16$
3	$1.67,\ 9$	$\{0, 2, 3\}$, $\{8, 10\}$	$6.67$

At iteration 3 the assignment no longer changes, the means stay where they are, and the algorithm stops with $J = 6.67$. In two dimensions the same alternation looks like this:

Figure

Lloyd's algorithm: assign, then move the centroids

The assignment step cuts the plane into Voronoi cells: the points nearer to $\mu_k$ than to $\mu_j$ form the half-plane on one side of the perpendicular bisector of $\mu_k$ and $\mu_j$, and a cell is an intersection of such half-planes, a convex polygon. Every k-means clustering is therefore a partition into convex pieces with straight boundaries, whatever shape the data have. This one sentence explains most of the failures below.

Initialization: k-means++#

Lloyd’s algorithm converges to whichever local minimum is downhill from its start. In the figure, two of the three starting centres sit in the same blob, and the algorithm recovers only because one of them can walk to the empty blob. With other starts it does not: on the four-blob data of the Code section, a single start from four random data points reaches the best $J = 1127.4$ in 56% of 200 runs, and the average over runs is $J = 4064.4$, because the failed runs merge two blobs and split another.

k-means++ (Arthur and Vassilvitskii, 2007) chooses the starting centres so that they are spread out. The first centre is a uniformly random data point. Each next centre is a data point drawn with probability proportional to $D(x)^2$, the squared distance from $x$ to the nearest centre chosen so far. Points far from all current centres, which are likely to lie in a cluster without a centre, are the likely picks. Their Theorem 3.1 bounds the result of the seeding alone, before any Lloyd iterations:

$$\mathbb E[J] \le 8(\ln K + 2)\, J_{\text{opt}} .$$

For $K = 4$ the factor is about 27, a loose bound, but a guarantee that no uniform start has; in practice the seeding lands much closer. On the same data every one of 200 single k-means++ starts reaches $J = 1127.4$. (scikit-learn’s k-means++ draws several candidates for each new centre and keeps the one that lowers $J$ most, a greedy variant.) The usual practice adds restarts: run from several seeds and keep the lowest $J$.

Worked example. For the points $0, 2, 3, 8, 10$ and a first centre at $0$, the squared distances are $0, 4, 9, 64, 100$, with sum $177$. The second centre is $10$ with probability $100/177 = 0.565$ and $8$ with probability $0.362$: with probability $0.93$ it lands in the far group, which a uniform draw would hit with probability $0.4$.

How many clusters?#

$J$ cannot choose $K$ by itself. The optimal $J$ never increases with $K$, since one more centre can always be placed on an existing one, and it reaches 0 at $K = n$, one cluster per point. Two common heuristics read $K$ off the data instead.

The elbow. Plot the best $J$ against $K$ and look for the point where adding clusters stops paying. On the four-blob data, $J = 18363,\ 5241,\ 1127,\ 1020,\ 921$ for $K = 2, \dots, 6$: the drop from 3 to 4 is large, the drops after 4 are small, and the elbow is at $K = 4$. On real data the curve often bends smoothly, and the elbow is a matter of taste.

The silhouette (Rousseeuw, 1987) scores each point by how much closer it is to its own cluster than to the next one. With $a(i)$ the mean distance from $x_i$ to the other points of its cluster and $b(i)$ the smallest mean distance from $x_i$ to the points of another cluster,

$$s(i) = \frac{b(i) - a(i)}{\max\big(a(i), b(i)\big)} \in [-1, 1] ,$$

with $s(i) = 0$ for a point alone in its cluster. Values near 1 mean well inside its cluster, near 0 on a border, negative probably in the wrong cluster. The average over all points scores the clustering and can be compared across $K$. For five points $0, 1, 2, 6, 7$ with clusters $\{0, 1, 2\}$ and $\{6, 7\}$, the point $2$ has $a = (2 + 1)/2 = 1.5$, $b = (4 + 5)/2 = 4.5$ and $s = (4.5 - 1.5)/4.5 = 0.67$. On the four-blob data the average silhouette is $0.588, 0.762, 0.811, 0.680, 0.577$ for $K = 2, \dots, 6$, highest at $K = 4$. The silhouette uses all pairwise distances, $O(n^2)$, and like $J$ it favours convex, well-separated clusters: it inherits the k-means notion of a cluster. A probabilistic alternative, a penalized likelihood, comes with the mixture models of the next note.

Where k-means fails#

The objective prefers clusters that are round, of similar spread and of similar size, and the Voronoi geometry cannot express anything else.

• Non-convex clusters. Two interleaving half-moons cannot be separated by a straight line, so k-means cuts each of them in half. On the 400-point moons of the Code section its adjusted Rand index with the true moons is 0.274.
• Unequal spreads. The boundary between two clusters is the bisector of their centres, halfway between them, wherever the spreads are. With a tight cluster (standard deviation 0.5) four units from a wide one (standard deviation 2), the outer points of the wide cluster that lie nearer to the tight centre are assigned to it; the adjusted Rand index is 0.561. A Gaussian mixture, which gives each cluster its own covariance, reaches 0.902 on the same data; it is the subject of the next note.
• Unequal sizes and elongated clusters. Splitting a large cluster can lower $J$ more than respecting a small one, and an elongated cluster is cut across its long axis for the same reason.
• Outliers. A mean moves toward far points, and a single far point can take a centre for itself. Replacing squared distances by absolute distances and means by medians (k-medians), or centres by data points (k-medoids), makes the method more robust at a higher cost.

Figure

Two moons: k-means against DBSCAN

DBSCAN: clusters as dense regions#

DBSCAN (Ester et al., 1996) defines clusters by density and needs two parameters, a radius $\varepsilon$ and a count $\text{minPts}$. Write $N_\varepsilon(x)$ for the data points within distance $\varepsilon$ of $x$, $x$ included.

• A point is a core point if $\lvert N_\varepsilon(x)\rvert \ge \text{minPts}$: it sits in a dense region.
• Two core points within $\varepsilon$ of each other belong to the same cluster, and so, by chaining, do all core points connected by a path of such steps. A cluster is a maximal set of core points connected this way, together with the non-core points within $\varepsilon$ of one of them, its border points.
• Every other point is noise and belongs to no cluster.

The algorithm grows clusters by breadth-first search from unvisited core points. Clusters can have any shape, because a chain of $\varepsilon$-steps can follow a curve; the number of clusters is an output, not an input; and outliers are labelled as such instead of being forced into a cluster. A border point within $\varepsilon$ of two clusters goes to whichever reaches it first, the only order dependence. Each step needs a neighbourhood query: with all pairwise distances the cost is $O(n^2)$, and a spatial index makes it much cheaper on low-dimensional data.

Worked example. On the line, take the points $0, 0.5, 1, 1.5, 5, 9, 9.4, 9.8$ with $\varepsilon = 0.6$ and $\text{minPts} = 3$. The neighbourhood of $0.5$ is $\{0, 0.5, 1\}$ and that of $1$ is $\{0.5, 1, 1.5\}$, so both are core points; $0$ and $1.5$ have only two points within $0.6$ and are border points of the same cluster. Likewise $9.4$ is a core point with border points $9$ and $9.8$. The point $5$ has no neighbour: it is noise. Two clusters and one outlier, with no $K$ given.

Choosing the parameters. The original paper sets $\text{minPts} = 4$ for two-dimensional data and reads $\varepsilon$ off the sorted $k$-distance graph: the distance of every point to its $k$-th nearest neighbour, $k = \text{minPts}$, sorted in decreasing order. Noise points have large $k$-distances, cluster points small ones, and $\varepsilon$ goes at the knee between them (Section 4.2). The method is sensitive to $\varepsilon$: on the moons of the Code section, $\varepsilon = 0.1$ gives four clusters and 10 noise points, $\varepsilon = 0.2$ recovers the two moons exactly, and $\varepsilon = 0.3$ merges everything into one cluster.

One density for all clusters. A single $\varepsilon$ fixes one density threshold. When one cluster is much denser than another, a radius small enough to separate the dense cluster from its neighbours declares the sparse cluster noise. HDBSCAN (Campello, Moulavi and Sander, 2013) considers all radii at once: the clusterings for decreasing $\varepsilon$ form a hierarchy, and it keeps the clusters that are most stable, those that persist over a wide range of radii.

Hierarchical clustering and linkages#

Optional depth: a whole hierarchy of clusterings and the linkages that build it. The rest of the note does not depend on it.

Often there is no single right $K$: documents form topics, topics split into subtopics. Agglomerative clustering builds the whole hierarchy. It starts with $n$ clusters of one point each and repeatedly merges the two closest clusters until one remains. The $n - 1$ merges form a tree, the dendrogram, whose merge heights record how far apart the merged clusters were; cutting it at a height, or at the level where $K$ clusters remain, gives a flat clustering.

“Closest clusters” needs a distance between sets of points, the linkage:

Linkage	Distance between clusters $A$ and $B$	Behaviour
single	$\min_{a \in A, b \in B} \lVert a - b \rVert$	follows chains; finds curved clusters, but one bridge of noise merges two clusters
complete	$\max_{a \in A, b \in B} \lVert a - b \rVert$	compact clusters of similar diameter
average	mean of $\lVert a - b \rVert$ over all pairs	between single and complete
Ward	increase of $J$ caused by the merge	the k-means objective, greedily

Ward’s linkage (Ward, 1963) merges the pair that increases the within-cluster sum of squares least. By the decomposition of the k-means section, merging $A$ and $B$ increases it by

$$\Delta(A, B) = \frac{n_A n_B}{n_A + n_B}\, \lVert \bar x_A - \bar x_B \rVert^2 .$$

For $A = \{0, 2, 3\}$ and $B = \{8, 10\}$ from the k-means example, $\Delta = \frac{3 \cdot 2}{5}(9 - 1.67)^2 = 64.53$, exactly the difference between the scatter of all five points around their mean, $71.2$, and the within-cluster scatter $6.67$.

Single linkage is the minimum spanning tree in disguise: its clustering into $K$ groups is obtained by deleting the $K - 1$ longest edges of a minimum spanning tree of the points (Gower and Ross, 1969). That is why it follows chains of near neighbours along the moons (adjusted Rand index 1.0 in the Code section), while complete, average and Ward linkage cut them like k-means does (0.285, 0.329 and 0.239). On round blobs all four recover the clusters. The price of the hierarchy is the matrix of pairwise distances, $O(n^2)$ memory, and between $O(n^2)$ and $O(n^3)$ time depending on the linkage, which limits it to tens of thousands of points.

Spectral clustering#

Optional depth: clusters as weakly connected parts of a similarity graph, found with the eigenvectors of its Laplacian; it assumes eigenvalues and eigenvectors of a symmetric matrix. The rest of the note does not depend on it.

Single linkage looks at chains; spectral clustering looks at the whole graph of near neighbours and finds the parts with few edges between them. Build a similarity graph on the points, with weights $w_{ij} \ge 0$: for example $w_{ij} = 1$ if $x_j$ is among the $m$ nearest neighbours of $x_i$ or the other way round, or $w_{ij} = \exp(-\lVert x_i - x_j \rVert^2 / 2\sigma^2)$. With the degrees $d_i = \sum_j w_{ij}$ and $D = \operatorname{diag}(d_1, \dots, d_n)$, the graph Laplacian is $L = D - W$. For any vector $f \in \mathbb R^n$, one value per point,

$$f^\top L f = \frac12 \sum_{i,j=1}^n w_{ij} (f_i - f_j)^2 ,$$

which is small when $f$ is nearly constant on strongly connected points. So $L$ is symmetric and positive semidefinite, the constant vector has $Lf = 0$, and a vector $f$ with $f^\top L f = 0$ must be constant along every edge. Hence the key fact (von Luxburg, 2007, Proposition 2): the multiplicity of the eigenvalue 0 of $L$ equals the number of connected components of the graph, and the indicator vectors of the components span its eigenspace.

If the graph falls into exactly $K$ components, the $K$ eigenvectors with eigenvalue 0 are, up to rotation, the indicators of the components; in the coordinates they give each point, all points of a component coincide. If the components are joined by a few weak edges, the eigenvalues are small rather than zero and the eigenvectors are nearly piecewise constant. The algorithm uses exactly that:

1. Build the similarity graph and its Laplacian.
2. Compute the eigenvectors $u_1, \dots, u_K$ for the $K$ smallest eigenvalues.
3. Give point $i$ the new coordinates $(u_{1i}, \dots, u_{Ki}) \in \mathbb R^K$ and run k-means on them.

In the new coordinates the moons become two tight blobs, and k-means separates them. On the moons of the Code section, the 10-nearest-neighbour graph has exactly two components, the Laplacian has two zero eigenvalues followed by $0.02$, and k-means on the two eigenvectors recovers the moons (adjusted Rand index 1.0). The method can also be derived as a relaxation of a graph cut: finding $K$ groups with little edge weight between them relative to their size is NP-hard (RatioCut, normalized cut), and relaxing the discrete indicator vectors to real vectors turns it into the eigenvalue problem above; normalized versions of $L$ correspond to the normalized cut of Shi and Malik (2000) and the algorithm of Ng, Jordan and Weiss (2001) (von Luxburg, Section 5). The cost is the eigendecomposition, $O(n^3)$ for a dense matrix and much less for a sparse nearest-neighbour graph, and the answer depends on how the graph is built.

Evaluating a clustering#

Without labels, a clustering can only be judged against the definition it was built from. $J$ compares solutions with the same $K$; the silhouette compares across $K$ but inherits the k-means notion of a cluster: on the moons of the Code section it scores the true split $0.33$ and the k-means split, which cuts both moons, $0.49$. A practical check is stability: cluster several bootstrap samples or random subsets and measure how much the results agree. Stability is evidence, not proof. Structure that changes with every resample is suspect, but real clusters that overlap or are weak can be unstable too, and a stable answer can be imposed by the method: with $K = 2$, k-means cuts a single elongated Gaussian blob across its long axis into nearly the same two halves on every resample. Whether a clustering is useful also depends on the task and on the representation the distances were computed in.

With labels for some data, withheld from the clustering (the points themselves are clustered, only their labels are hidden), the clustering can be compared with them, but not by accuracy: cluster numbers are arbitrary, and the clustering $\{1, 1, 2, 2\}$ is identical to $\{2, 2, 1, 1\}$. Pair-counting measures avoid the matching. Over all $\binom n2$ pairs of points, count the pairs on which the two partitions agree (both put the pair together, or both apart). The share of agreements is the Rand index. It has a flaw: two random partitions agree on many pairs by chance, most of them “apart”, so its values crowd near 1. The adjusted Rand index (ARI) of Hubert and Arabie (1985) subtracts the agreement expected by chance. With $n_{ij}$ the number of points in true class $i$ and cluster $j$, $a_i$ and $b_j$ the row and column sums,

$$\text{ARI} = \frac{\sum_{ij} \binom{n_{ij}}{2} - E}{\frac12\big[\sum_i \binom{a_i}{2} + \sum_j \binom{b_j}{2}\big] - E}, \qquad E = \frac{\sum_i \binom{a_i}{2} \sum_j \binom{b_j}{2}}{\binom n2} .$$

$E$ is the expected number of pairs together in both partitions if the cluster labels were shuffled with the cluster sizes kept, a hypergeometric model of chance. The ARI is 1 for identical partitions, about 0 for unrelated ones, and can be negative. Points that DBSCAN calls noise need a convention. The Code section passes the noise label $-1$ to the ARI as if it were one more cluster: for $\varepsilon = 0.1$ that gives $0.731$, while scoring only the 390 clustered points gives $0.769$, a higher number that ignores the 10 discarded points. A reported ARI should say which one it is.

Worked example. Six points with true classes $AAABBB$ and clusters $112222$. The table $n_{ij}$ has rows $(2, 1)$ and $(0, 3)$, so $\sum \binom{n_{ij}}{2} = 1 + 0 + 0 + 3 = 4$; the classes give $\binom32 + \binom32 = 6$, the clusters $\binom22 + \binom42 = 7$, and $\binom62 = 15$. Then $E = 6 \cdot 7 / 15 = 2.8$ and $\text{ARI} = (4 - 2.8)/(6.5 - 2.8) = 0.32$. The Rand index of the same pair is $10/15 = 0.67$, which sounds far better than the clustering is.

Labels are one structure among many. Images of digits can be grouped by digit, by writer or by stroke width; a low ARI against the digit labels may mean a different and still useful grouping. An external index measures agreement with a chosen structure, not quality in general.

Where k-means reappears#

k-means is the standard way to turn continuous vectors into a small set of discrete symbols, and the generative models course uses it that way several times. The VQ-VAE replaces each encoder vector by the nearest entry of a learned codebook, and its moving-average codebook update is k-means computed online, with the same failure: a centre that no point selects never moves. Neural audio codecs stack such codebooks in residual vector quantization, each stage clustering what the previous ones left. ImageGPT turns pixels into tokens with a palette of 512 colours found by k-means, and the precision–recall curve for generative models compares histograms over k-means clusters of real and generated embeddings. In vector search, the approximate indexes mentioned in the semantic search note include product quantization (Jégou, Douze and Schmid, 2011), which splits each vector into pieces and replaces each piece by the nearest of a set of k-means centroids.

What changed#

	Supervised learning	k-means	DBSCAN	Hierarchical	Spectral
Input	pairs $(x, y)$	points, $K$	points, $\varepsilon$, minPts	points, linkage	similarity graph, $K$
A cluster is		a Voronoi cell around a mean	a dense connected region	a subtree of the dendrogram	a weakly connected part of the graph
Shapes		convex	any	any (single) or compact	any
Outliers		pull centres	labelled as noise	singleton branches (single)	depends on the graph
Cost		$O(nKd)$ per iteration	$O(n^2)$ without an index	$O(n^2)$ memory	eigenvectors of an $n \times n$ matrix
Judged by	test error	$J$, silhouette	ARI with labels, stability	the dendrogram, ARI	ARI with labels, stability

Limitations#

• The answer depends on the definition. Every method encodes its own notion of a cluster, and no method is right on all data; the choice is part of the modelling.
• Distances must be meaningful. Feature scaling and encoding change every result, and in high dimensions distances concentrate; clustering an embedding learned for the purpose, or a few principal components, often works better than clustering raw features.
• Hyperparameters without a target. $K$, $\varepsilon$, minPts and the graph construction cannot be tuned by validation error, because there is none; the heuristics for choosing them encode the same assumptions as the methods.
• Hard assignments. Each method above puts a point in exactly one cluster (or none), with no measure of how sure it is. A point on the border of two clusters is treated like a point at a centre.
• Local optima and cost. k-means finds a local minimum of an NP-hard problem; hierarchical and spectral clustering scale quadratically or worse in $n$.

What you should remember#

• A clustering needs a definition of a cluster and a distance; different definitions give different, equally valid answers (Kleinberg’s impossibility theorem).
• k-means minimizes $J = \sum_i \lVert x_i - \mu_{c_i} \rVert^2$. Lloyd’s algorithm alternates nearest-centre assignment and mean updates; it is coordinate descent, $J$ never increases, and it stops at a local minimum.
• k-means++ seeds with probability proportional to $D(x)^2$ and guarantees $\mathbb E[J] \le 8(\ln K + 2) J_{\text{opt}}$; use it with a few restarts.
• Choose $K$ with the elbow of $J$ or the average silhouette $s = (b - a)/\max(a, b)$; $J$ alone always prefers more clusters.
• k-means cells are convex and split the space halfway between centres, so it fails on curved clusters, unequal spreads and outliers.
• DBSCAN: core points have at least minPts neighbours within $\varepsilon$; clusters are chains of core points plus their border points; the rest is noise.
• Agglomerative clustering builds a dendrogram; single linkage follows chains (a minimum spanning tree), Ward’s linkage greedily minimizes $J$.
• Spectral clustering embeds the points with the bottom eigenvectors of $L = D - W$, whose zero eigenvalues count the connected components, and runs k-means there.
• Evaluate against labels with the adjusted Rand index, which corrects pair agreement for chance; without labels use the silhouette and stability, knowing what they assume.

Code#

Lloyd’s algorithm, seeding, the choice of $K$, the failure cases, and the other methods on the moons, with scikit-learn:

import numpy as np
from scipy.sparse.csgraph import laplacian
from sklearn.datasets import make_blobs, make_moons
from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering
from sklearn.mixture import GaussianMixture
from sklearn.neighbors import kneighbors_graph
from sklearn.metrics import adjusted_rand_score as ari, silhouette_score

X, y = make_blobs(n_samples=600, centers=4, cluster_std=1.0, random_state=7)

# Lloyd's algorithm by hand from four random data points: J never increases
rng = np.random.default_rng(0)
mu = X[rng.choice(len(X), 4, replace=False)]
for t in range(4):
    d = ((X[:, None, :] - mu[None]) ** 2).sum(axis=2)          # squared distances, n x K
    c = d.argmin(axis=1)                                         # assignment step
    print(t, round(d[np.arange(len(X)), c].sum(), 1))
    mu = np.array([X[c == k].mean(axis=0) for k in range(4)])   # update step: cluster means
# 0 22942.3 | 1 2406.5 | 2 1127.4 | 3 1127.4

# one random start against one k-means++ start, 200 seeds each: mean J and share of runs at the best J
for init in ('random', 'k-means++'):
    J = np.array([KMeans(4, init=init, n_init=1, random_state=s).fit(X).inertia_ for s in range(200)])
    print(init, round(J.mean(), 1), (J < 1.001 * J.min()).mean())
# random 4064.4 0.56 | k-means++ 1127.4 1.0

# the elbow and the silhouette for K = 2..7
for K in range(2, 8):
    km = KMeans(K, n_init=10, random_state=0).fit(X)
    print(K, round(km.inertia_, 1), round(silhouette_score(X, km.labels_), 3))
# 2 18363.1 0.588 | 3 5241.2 0.762 | 4 1127.4 0.811 | 5 1019.9 0.68 | 6 921.1 0.577 | 7 813.9 0.446

# the answer depends on the units: the second feature multiplied by 10
km = lambda Z, K: KMeans(K, n_init=10, random_state=0).fit_predict(Z)
print(ari(y, km(X, 4)), round(ari(y, km(X * [1, 10], 4)), 3))       # 1.0 0.736

# unequal spreads: k-means against a Gaussian mixture (the next note)
Xu, yu = make_blobs(n_samples=[200, 200], centers=[[0, 0], [4, 0]], cluster_std=[0.5, 2.0], random_state=2)
print(round(ari(yu, km(Xu, 2)), 3), round(ari(yu, GaussianMixture(2, random_state=0).fit(Xu).predict(Xu)), 3))   # 0.561 0.902

# two moons: k-means, DBSCAN, four linkages
Xm, ym = make_moons(n_samples=400, noise=0.06, random_state=0)
print(round(ari(ym, km(Xm, 2)), 3))                                  # 0.274
print(round(silhouette_score(Xm, ym), 2), round(silhouette_score(Xm, km(Xm, 2)), 2))   # 0.33 0.49: the silhouette prefers k-means
for eps in (0.1, 0.2, 0.3):
    lab = DBSCAN(eps=eps, min_samples=5).fit_predict(Xm)
    print(eps, round(ari(ym, lab), 3), lab.max() + 1, (lab == -1).sum())   # ARI, clusters, noise points
# 0.1 0.731 4 10 | 0.2 1.0 2 0 | 0.3 0.0 1 0
lab = DBSCAN(eps=0.1, min_samples=5).fit_predict(Xm); kept = lab != -1
print(round(ari(ym[kept], lab[kept]), 3), kept.sum())               # 0.769 390: noise points left out
for link in ('single', 'complete', 'average', 'ward'):
    agg = lambda Z, K: AgglomerativeClustering(K, linkage=link).fit_predict(Z)
    print(link, round(ari(ym, agg(Xm, 2)), 3), round(ari(y, agg(X, 4)), 3))
# single 1.0 1.0 | complete 0.285 1.0 | average 0.329 1.0 | ward 0.239 1.0

# spectral clustering by hand: a 10-nearest-neighbour graph, its Laplacian, k-means on two eigenvectors
W = kneighbors_graph(Xm, 10).toarray()
W = np.maximum(W, W.T)                                               # an edge if either point is a neighbour of the other
lam, V = np.linalg.eigh(laplacian(W))                                # L = D - W
print(lam[:3].round(4), round(ari(ym, km(V[:, :2], 2)), 3))          # [0. 0. 0.02] 1.0: two components, two zeros

Next#

Every method in this note gives hard answers: a point belongs to one cluster, and nothing says how sure that is or how spread out each cluster may be. k-means, in particular, behaves as if all clusters were round and equally wide. The next note makes clustering probabilistic. Each cluster becomes a Gaussian with its own mean and covariance, each point gets a probability of belonging to each cluster, and the parameters are fitted by maximum likelihood with the EM algorithm, of which Lloyd’s algorithm will turn out to be a limiting case.

References#

• J. Kleinberg. An Impossibility Theorem for Clustering. NeurIPS, 2002. Section 2 (the three properties, Theorem 2.1, Theorem 2.2).
• S. P. Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2), 129–137, 1982.
• J. MacQueen. Some methods for classification and analysis of multivariate observations. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, 1, 281–297, 1967.
• D. Aloise, A. Deshpande, P. Hansen, P. Popat. NP-hardness of Euclidean sum-of-squares clustering. Machine Learning, 75(2), 245–248, 2009.
• D. Arthur, S. Vassilvitskii. k-means++: The Advantages of Careful Seeding. SODA, 2007. Section 2.2 (the seeding), Theorem 3.1.
• P. J. Rousseeuw. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics, 20, 53–65, 1987.
• M. Ester, H.-P. Kriegel, J. Sander, X. Xu. A density-based algorithm for discovering clusters in large spatial databases with noise. KDD, 226–231, 1996. Section 3 (definitions), Section 4.2 (choosing Eps and MinPts).
• R. J. G. B. Campello, D. Moulavi, J. Sander. Density-Based Clustering Based on Hierarchical Density Estimates. PAKDD, LNCS 7819, 160–172, 2013.
• J. H. Ward. Hierarchical Grouping to Optimize an Objective Function. Journal of the American Statistical Association, 58(301), 236–244, 1963.
• J. C. Gower, G. J. S. Ross. Minimum Spanning Trees and Single Linkage Cluster Analysis. Applied Statistics, 18(1), 54–64, 1969.
• U. von Luxburg. A Tutorial on Spectral Clustering. Statistics and Computing, 17(4), 2007. Section 3 (Propositions 1 and 2), Section 4 (the algorithms), Section 5 (the graph-cut view).
• J. Shi, J. Malik. Normalized cuts and image segmentation. IEEE TPAMI, 22(8), 888–905, 2000.
• A. Y. Ng, M. I. Jordan, Y. Weiss. On Spectral Clustering: Analysis and an algorithm. NeurIPS, 2001.
• L. Hubert, P. Arabie. Comparing partitions. Journal of Classification, 2(1), 193–218, 1985.
• H. Jégou, M. Douze, C. Schmid. Product Quantization for Nearest Neighbor Search. IEEE TPAMI, 33(1), 117–128, 2011.
• T. Hastie, R. Tibshirani, J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. Section 14.3 (cluster analysis).
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Section 9.1 (k-means clustering).
• E. Sokolov et al. Машинное обучение (HSE), lecture notes: lecture 12, unsupervised learning (2021); lecture 17, clustering (2025).