45.08.08 · mathematical-statistics / 08-learning-methods

k-Means and Model-Based Clustering

shipped3 tiersLean: none

Anchor (Master): Hastie, Tibshirani & Friedman 2009 The Elements of Statistical Learning 2e (Springer) §14.3.6-14.3.12 (model-based clustering with Gaussian mixtures and BIC, hierarchical/agglomerative clustering and the dendrogram, spectral clustering via the graph Laplacian); Arthur & Vassilvitskii 2007 Proc. SODA pp. 1027-1035 (k-means++: the $D^{2}$ seeding and its $O (lo g k)$ -competitive expected approximation guarantee); Aloise, Deshpande, Hansen & Popat 2009 Machine Learning 75:245-248 (NP-hardness of Euclidean sum-of-squares clustering for general dimension)

Intuition Beginner

Imagine you have a scatter of dots on a page and you suspect they fall into a few natural clumps, but nobody has labelled which clump each dot belongs to. Clustering is the task of discovering those clumps from the dots alone. The most common recipe, called k-means, fixes the number of clumps in advance — say three — and then hunts for three good "centre" points and an assignment of every dot to its nearest centre.

The recipe runs as a back-and-forth loop. Drop three centres down somewhere to start. First move: hand each dot to whichever centre is closest to it, splitting the page into three territories. Second move: now that each centre owns a set of dots, slide that centre to the exact middle (the average position) of the dots it owns. Sliding the centres changes who is nearest to what, so repeat the first move, then the second, and keep going. Each pass tidies the clumps a little more.

The loop stops when nothing moves: every dot is already closest to its own centre, and every centre already sits at the middle of its dots. At that point you have a clustering. The quantity the loop quietly shrinks at every step is the total spread — the sum, over all dots, of the squared distance from each dot to its centre. Tight clumps mean a small total spread, so driving that number down is the same as finding compact groups.

One catch: where you drop the starting centres changes where you end up. A bad start can settle into a poor clumping, so people run the loop several times from different starts and keep the best.

Visual Beginner

Picture the page after a few passes. There are three small crosses for the three centres, and every dot is coloured to match the cross it was handed to. The colour boundaries look like straight fences halfway between neighbouring crosses, so the page is carved into three straight-edged territories, and each cross sits right in the middle of its own colour.

step	what you do	what comes out
assign	with the centres fixed, colour each dot by its nearest centre	a grouping of the dots into territories
update	with the colours fixed, move each centre to the average of its dots	new centre positions

The takeaway: k-means is two alternating moves — colour the dots by nearest centre, then move each centre to the middle of its dots — and repeating the pair shrinks the total spread until the picture stops changing.

Worked example Beginner

Six numbers on a line: $1, 2, 3, 10, 11, 12$ . We look for two clusters ( $k = 2$ ). Start with centres $c_{1} = 2$ and $c_{2} = 9$ .

Step 1, assign each number to its nearer centre. Numbers $1, 2, 3$ are closer to $c_{1} = 2$ (distances $1, 0, 1$ ) than to $c_{2} = 9$ . Numbers $10, 11, 12$ are closer to $c_{2} = 9$ . So cluster one is ${1, 2, 3}$ and cluster two is ${10, 11, 12}$ .

Step 2, recompute each centre as the average of its cluster. Cluster one averages $(1 + 2 + 3) /3 = 6/3 = 2$ . Cluster two averages $(10 + 11 + 12) /3 = 33/3 = 11$ . The centres move to $c_{1} = 2$ and $c_{2} = 11$ .

Step 3, reassign with the new centres. Numbers $1, 2, 3$ are still nearest $c_{1} = 2$ ; numbers $10, 11, 12$ are still nearest $c_{2} = 11$ . Nobody switches teams, so the next update leaves the centres unchanged. The loop has settled.

Now total the spread to see what was minimised. With $c_{1} = 2$ , cluster one contributes $(1 - 2)^{2} + (2 - 2)^{2} + (3 - 2)^{2} = 1 + 0 + 1 = 2$ . With $c_{2} = 11$ , cluster two contributes $(10 - 11)^{2} + (11 - 11)^{2} + (12 - 11)^{2} = 1 + 0 + 1 = 2$ . The total spread is $2 + 2 = 4$ — small, because the two clumps are tight and far apart.

What this tells us: from a rough start, one assign-then-update pass already found the two obvious groups, and a second pass confirmed nothing should change. The number the loop drove down, the total squared distance to centres, ended at $4$ , the smallest it can be for this split.

Check your understanding Beginner

Exercise (easy, multiple-choice).

What happens to the total spread (the sum of squared distances from points to their centres) over the course of the k-means loop?

It can rise or fall unpredictably from pass to pass.
It never increases, and it stops changing once the loop settles.
It always reaches the single best possible value, no matter the start.
It always increases until the loop stops.

Hint

Each of the two moves is chosen precisely to lower, or at worst keep, the spread.

Answer

B. Both the assign move and the update move can only lower or hold the total spread, so it never increases and levels off when the loop settles. It is not guaranteed to reach the single best value — different starts can settle at different clusterings. Feedback-correct: monotone non-increase is the defining behaviour of the loop. Feedback-wrong: A and D misstate the direction; C ignores the dependence on the starting centres.

Formal definition Intermediate+

Let $x_{1}, \dots, x_{n} \in R^{p}$ be observations and fix the number of clusters $K$ . A clustering is an assignment map $C : {1, \dots, n} \to {1, \dots, K}$ partitioning the indices into clusters $C_{k} = {i : C (i) = k}$ of sizes $N_{k} = ∣ C_{k} ∣$ .

Definition (within-cluster scatter; k-means objective). With squared Euclidean dissimilarity, the pooled within-cluster sum of squares is

W (C) = k = 1 \sum K i \in C_{k} \sum x_{i} - \overset{x}{ˉ}_{k}^{2}, \overset{x}{ˉ}_{k} = \frac{1}{N _{k}} i \in C_{k} \sum x_{i},

where $\overset{x}{ˉ}_{k}$ is the cluster centroid. Introducing free representatives $m_{1}, \dots, m_{K} \in R^{p}$ , this is the minimised value of the joint objective

J (C, m) = k = 1 \sum K i \in C_{k} \sum x_{i} - m_{k}^{2} = i = 1 \sum n x_{i} - m_{C (i)}^{2},

so $W (C) = min_{m} J (C, m)$ . The k-means problem is $min_{C, m} J (C, m)$ , equivalently $min_{C} W (C)$ .

Definition (Lloyd's algorithm). Starting from initial representatives $m^{(0)} = (m_{1}^{(0)}, \dots, m_{K}^{(0)})$ , iterate for $t = 0, 1, 2, \dots$ :

Assignment step. Holding $m^{(t)}$ fixed, assign each point to its nearest representative,

C^{(t + 1)} (i) \in ar g k min x_{i} - m_{k}^{(t)}^{2}, i = 1, \dots, n .

Update step. Holding $C^{(t + 1)}$ fixed, set each representative to its cluster mean,

m_{k}^{(t + 1)} = \frac{1}{N _{k}^{(t + 1)}} i \in C_{k}^{(t + 1)} \sum x_{i} .

Iterate until the assignment is unchanged. (Ties in the assignment step are broken by a fixed rule; an emptied cluster is reseeded, e.g. at the point farthest from its centre.)

The identity $\sum_{k} \sum_{i \in C_{k}} ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2} = \frac{1}{2} \sum_{k} \frac{1}{N _{k}} \sum_{i, i^{'} \in C_{k}} ∥ x_{i} - x_{i^{'}} ∥^{2}$ shows $W (C)$ is equivalently half the average within-cluster pairwise squared distance, weighted by cluster size; this is the form used when only a dissimilarity matrix is available.

Definition (k-medoids / PAM). Replacing the squared Euclidean distance by a general dissimilarity $d (x_{i}, x_{i^{'}})$ and constraining each representative to be one of the data points (a medoid $m_{k} \in {x_{i}}$ ) gives k-medoids; the update step becomes a search over candidate medoids minimising within-cluster total dissimilarity, which is robust to outliers and applies to non-Euclidean data.

Definition (Gaussian mixture model; model-based clustering). A $K$ -component Gaussian mixture has density $p (x ∣ θ) = \sum_{k = 1}^{K} π_{k} ϕ (x; μ_{k}, Σ_{k})$ with mixing weights $π_{k} \geq 0$ , $\sum_{k} π_{k} = 1$ , means $μ_{k}$ , covariances $Σ_{k}$ , and $ϕ (\cdot; μ, Σ)$ the multivariate normal density 45.08.07. Model-based clustering fits this mixture by maximum likelihood (via EM) and assigns each point to the component of highest posterior responsibility $γ_{ik} = π_{k} ϕ (x_{i}; μ_{k}, Σ_{k}) / \sum_{j} π_{j} ϕ (x_{i}; μ_{j}, Σ_{j})$ . Constraining the covariances yields a family of models — spherical $Σ_{k} = σ^{2} I$ , diagonal, full, and tied ( $Σ_{k} = Σ$ ) — of increasing flexibility and parameter count.

The new symbols are recorded in _meta/NOTATION.md: $W (C)$ (within-cluster scatter), $J (C, m)$ (joint k-means objective), $C (i)$ (assignment map), $\overset{x}{ˉ}_{k}$ (cluster centroid), $N_{k}$ (cluster size), $π_{k}, μ_{k}, Σ_{k}$ (mixture parameters), $γ_{ik}$ (responsibility).

Counterexamples to common slips Intermediate+

"Lloyd's algorithm finds the clustering with the smallest $W (C)$ ." It finds a local optimum: a partition no single assignment-or-update move can improve. The global k-means problem is NP-hard, and the local optimum reached depends on initialisation. Different starts give different, generally suboptimal, clusterings; multiple restarts (or k-means++ seeding) mitigate but do not remove this.
"k-means works with any distance." The centroid-update step is the minimiser of within-cluster scatter only for squared Euclidean distance; for an arbitrary dissimilarity the mean need not minimise total dissimilarity, and the monotone-descent guarantee fails. Use k-medoids, where the representative is an actual data point chosen to minimise total dissimilarity, for general $d$ .
"k-means recovers the true clusters if they exist." The objective rewards spherical, comparably sized, comparably scaled clusters. Elongated, unequal-variance, or unequal-population clusters are systematically mis-split because the nearest-centroid rule draws straight (Voronoi) boundaries; a Gaussian mixture with full covariances, which fits orientation and scale, is the corresponding fix.
"More clusters is always a better fit." $W (C)$ is non-increasing in $K$ and reaches zero at $K = n$ (every point its own cluster), so $W$ alone cannot choose $K$ . A complexity-aware criterion is required — the gap statistic for k-means, BIC for the mixture-model likelihood 45.06.09.

Key theorem with proof Intermediate+

The signature result is that Lloyd's algorithm is a coordinate-descent on $J (C, m)$ that decreases the objective at every step and terminates after finitely many iterations at a partition-local optimum. The two steps are the two block-minimisations of $J$ , and finiteness follows because $J$ takes finitely many values along the run ^{[Lloyd, S. P. — Least Squares Quantization in PCM]}.

Theorem (monotone descent and finite termination of Lloyd's algorithm). Let $J (C, m) = \sum_{i} ∥ x_{i} - m_{C (i)} ∥^{2}$ . The Lloyd iteration satisfies

J (C^{(t + 1)}, m^{(t + 1)}) \leq J (C^{(t)}, m^{(t)}) for every t,

and the algorithm terminates after finitely many iterations at a pair $(C^{⋆}, m^{⋆})$ that is simultaneously optimal in each block: $C^{⋆} (i) \in ar g min_{k} ∥ x_{i} - m_{k}^{⋆} ∥^{2}$ and $m_{k}^{⋆} = \overset{x}{ˉ}_{k}^{⋆}$ . Such a fixed point is a local minimum of $W (\cdot)$ in the sense that no reassignment of a single point and no recentring lowers $J$ .

Proof. Both steps minimise $J$ over one block with the other held fixed. The assignment step minimises $C \mapsto J (C, m^{(t)}) = \sum_{i} ∥ x_{i} - m_{C (i)}^{(t)} ∥^{2}$ : the summands are independent across $i$ , and each is minimised by sending $i$ to the nearest representative, which is exactly $C^{(t + 1)}$ . Hence

J (C^{(t + 1)}, m^{(t)}) \leq J (C^{(t)}, m^{(t)}) .

The update step minimises $m \mapsto J (C^{(t + 1)}, m) = \sum_{k} \sum_{i \in C_{k}^{(t + 1)}} ∥ x_{i} - m_{k} ∥^{2}$ , which separates over $k$ . For each $k$ the function $m_{k} \mapsto \sum_{i \in C_{k}} ∥ x_{i} - m_{k} ∥^{2}$ is a sum of squared distances; its gradient $- 2 \sum_{i \in C_{k}} (x_{i} - m_{k})$ vanishes at $m_{k} = \overset{x}{ˉ}_{k}$ , and the Hessian $2 N_{k} I$ is positive definite, so the centroid is the unique minimiser. Hence

J (C^{(t + 1)}, m^{(t + 1)}) \leq J (C^{(t + 1)}, m^{(t)}) \leq J (C^{(t)}, m^{(t)}),

establishing the monotone decrease. For termination: the value $J (C^{(t + 1)}, m^{(t + 1)}) = W (C^{(t + 1)})$ depends only on the assignment $C^{(t + 1)}$ , and there are at most $K^{n}$ assignments, so $J$ takes finitely many values along the run. A strictly decreasing finite sequence cannot repeat, so after finitely many steps no assignment change occurs; at that point $C^{⋆}$ already assigns every point to its nearest current centroid and $m_{k}^{⋆} = \overset{x}{ˉ}_{k}^{⋆}$ , the two fixed-point conditions. Because each summand is independently assigned-optimal and each centroid is recentring-optimal, no single reassignment and no recentring lowers $J$ . $□$

Bridge. This monotone-descent guarantee is the foundational reason Lloyd's algorithm is usable despite the global k-means problem being NP-hard, and it builds toward the model-based generalisation: k-means is exactly the hard-assignment, zero-variance limit of the EM algorithm for a spherical Gaussian mixture 45.08.07, where the assignment step is the degenerate E-step (responsibilities collapsing to nearest-centroid indicators) and the centroid update is the M-step. This is dual to the soft picture of the mixture model — finite component variance gives fractional responsibilities and a smooth descent on the log-likelihood rather than the combinatorial $W (C)$ . The central insight is that alternating block-minimisation of a single objective forces monotone descent whether the blocks are (assignment, centroids) here or (latent posterior, parameters) in EM, and this is exactly the coordinate-ascent structure proved for EM in 45.08.07. The two-condition fixed point — nearest-neighbour partition and centroid representative — generalises Lloyd's original quantizer-optimality conditions, and putting these together, the descent terminates only at a partition-local optimum, which is why initialisation (the k-means++ seeding of the Advanced results) and the choice of $K$ (the gap statistic and BIC) appear again as the two questions the descent itself cannot answer.

Exercises Intermediate+

Exercise 3 (medium, symbolic).

Prove the identity $\sum_{i \in C_{k}} ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2} = \frac{1}{2 N _{k}} \sum_{i, i^{'} \in C_{k}} ∥ x_{i} - x_{i^{'}} ∥^{2}$ , so that the within-cluster scatter equals (half, per point) the average within-cluster pairwise squared distance.

Hint

Expand $∥ x_{i} - x_{i^{'}} ∥^{2} = ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2} + ∥ x_{i^{'}} - \overset{x}{ˉ}_{k} ∥^{2} - 2 ⟨ x_{i} - \overset{x}{ˉ}_{k}, x_{i^{'}} - \overset{x}{ˉ}_{k} ⟩$ and sum, using $\sum_{i \in C_{k}} (x_{i} - \overset{x}{ˉ}_{k}) = 0$ .

Answer

Write $u_{i} = x_{i} - \overset{x}{ˉ}_{k}$ , so $\sum_{i \in C_{k}} u_{i} = 0$ . Then $\sum_{i, i^{'}} ∥ x_{i} - x_{i^{'}} ∥^{2} = \sum_{i, i^{'}} ∥ u_{i} - u_{i^{'}} ∥^{2} = \sum_{i, i^{'}} (∥ u_{i} ∥^{2} + ∥ u_{i^{'}} ∥^{2} - 2 ⟨ u_{i}, u_{i^{'}} ⟩)$ . The first two terms each give $N_{k} \sum_{i} ∥ u_{i} ∥^{2}$ , summing to $2 N_{k} \sum_{i} ∥ u_{i} ∥^{2}$ . The cross term is $- 2 \sum_{i, i^{'}} ⟨ u_{i}, u_{i^{'}} ⟩ = - 2 ⟨ \sum_{i} u_{i}, \sum_{i^{'}} u_{i^{'}} ⟩ = - 2∥ \sum_{i} u_{i} ∥^{2} = 0$ . Hence $\sum_{i, i^{'}} ∥ x_{i} - x_{i^{'}} ∥^{2} = 2 N_{k} \sum_{i} ∥ u_{i} ∥^{2} = 2 N_{k} \sum_{i} ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2}$ , which rearranges to the claim. Full credit for the expansion, the vanishing cross term via $\sum u_{i} = 0$ , and the rearrangement.

Exercise 4 (medium, symbolic).

Show that for a single cluster ${x_{i}}_{i \in C_{k}}$ the representative minimising $\sum_{i \in C_{k}} ∥ x_{i} - m ∥^{2}$ is the centroid $\overset{x}{ˉ}_{k}$ , and that the minimised value is $\sum_{i} ∥ x_{i} ∥^{2} - N_{k} ∥ \overset{x}{ˉ}_{k} ∥^{2}$ .

Hint

Set the gradient in $m$ to zero; then expand the squared norms at $m = \overset{x}{ˉ}_{k}$ .

Answer

The gradient $\nabla_{m} \sum_{i} ∥ x_{i} - m ∥^{2} = - 2 \sum_{i} (x_{i} - m) = - 2 (\sum_{i} x_{i} - N_{k} m)$ vanishes at $m = \frac{1}{N _{k}} \sum_{i} x_{i} = \overset{x}{ˉ}_{k}$ , and the Hessian $2 N_{k} I ≻ 0$ confirms the unique minimum. At $m = \overset{x}{ˉ}_{k}$ , $\sum_{i} ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2} = \sum_{i} ∥ x_{i} ∥^{2} - 2 ⟨ \sum_{i} x_{i}, \overset{x}{ˉ}_{k} ⟩ + N_{k} ∥ \overset{x}{ˉ}_{k} ∥^{2} = \sum_{i} ∥ x_{i} ∥^{2} - 2 N_{k} ∥ \overset{x}{ˉ}_{k} ∥^{2} + N_{k} ∥ \overset{x}{ˉ}_{k} ∥^{2} = \sum_{i} ∥ x_{i} ∥^{2} - N_{k} ∥ \overset{x}{ˉ}_{k} ∥^{2}$ , using $\sum_{i} x_{i} = N_{k} \overset{x}{ˉ}_{k}$ . Full credit for the first-order condition and the simplified minimum.

Exercise 6 (hard, symbolic).

Show that hard k-means (with $Σ_{k} = σ^{2} I$ , equal weights) is the $σ^{2} \to 0$ limit of EM for a spherical Gaussian mixture: the responsibilities $γ_{ik}$ converge to nearest-centroid indicators and the M-step mean reduces to the cluster average.

Hint

Write $γ_{ik} = exp (- ∥ x_{i} - μ_{k} ∥^{2} /2 σ^{2}) / \sum_{j} exp (- ∥ x_{i} - μ_{j} ∥^{2} /2 σ^{2})$ and divide through by the numerator of the nearest centre.

Answer

With equal weights and $Σ_{k} = σ^{2} I$ , $γ_{ik} = exp (- d_{ik} /2 σ^{2}) / \sum_{j} exp (- d_{ij} /2 σ^{2})$ where $d_{ij} = ∥ x_{i} - μ_{j} ∥^{2}$ . Let $k_{i}^{⋆} = ar g min_{j} d_{ij}$ (assumed unique) and factor out $exp (- d_{i k_{i}^{⋆}} /2 σ^{2})$ : $γ_{ik} = exp (- (d_{ik} - d_{i k_{i}^{⋆}}) /2 σ^{2}) / \sum_{j} exp (- (d_{ij} - d_{i k_{i}^{⋆}}) /2 σ^{2})$ . For $k = k_{i}^{⋆}$ the numerator is $1$ and every other denominator term has $d_{ij} - d_{i k_{i}^{⋆}} > 0$ , so its exponent $\to - \infty$ as $σ^{2} \to 0$ , leaving the denominator $\to 1$ ; hence $γ_{i k_{i}^{⋆}} \to 1$ . For $k \neq = k_{i}^{⋆}$ the numerator $\to 0$ , so $γ_{ik} \to 0$ . The responsibilities become $1 {k = k_{i}^{⋆}}$ , exactly the nearest-centroid assignment; the M-step mean $μ_{k} = \sum_{i} γ_{ik} x_{i} / \sum_{i} γ_{ik}$ then becomes the unweighted average of the points assigned to $k$ , the Lloyd update. Full credit for the limit of the responsibilities and the reduction of the M-step to the cluster mean.

Exercise 7 (hard, short-answer).

Lloyd's algorithm decreases $W (C)$ monotonically, yet the global k-means problem is NP-hard. Explain why these facts are not in tension, and what the monotone-descent guarantee does and does not promise about the solution found.

Hint

Distinguish convergence to a fixed point of the iteration from convergence to the global minimiser of a non-convex combinatorial objective.

Answer

Monotone descent guarantees only that each step does not increase $W$ and that the iteration reaches, in finitely many steps, a partition-local optimum — a fixed point where no single reassignment or recentring helps. It says nothing about how that local value compares to the global minimum. The objective $W (C)$ over partitions is non-convex and combinatorial; it has many local optima, and which one the descent reaches is determined by the initialisation. NP-hardness is a statement about the worst-case difficulty of finding the global minimiser over all partitions (proved for sum-of-squares clustering even with $k = 2$ in general dimension), and a polynomial-time descent to a local optimum does not contradict it — the descent simply offers no approximation guarantee on its own. This gap is exactly what k-means++ addresses, providing an $O (lo g k)$ -competitive seeding before the descent so that the local optimum reached is, in expectation, within a logarithmic factor of optimal. Full credit for separating local-fixed-point convergence from global optimisation and for noting that monotone descent carries no global guarantee absent careful seeding.

Advanced results Master

The results below organise clustering around the descent theorem of the Key theorem: the seeding that controls which local optimum the descent reaches, the model-based mixture view and its complexity-penalised selection of the number of components, and the two structurally different families — hierarchical and spectral — that abandon a fixed- $K$ centroid objective entirely. The hard/soft EM connection of 45.08.07 and the model-selection machinery of 45.06.09 run throughout.

Theorem 1 (NP-hardness of the global k-means problem). Minimising $W (C) = min_{C} \sum_{k} \sum_{i \in C_{k}} ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2}$ over all partitions into $K$ clusters is NP-hard: for $K = 2$ in dimension $p$ growing with $n$ (Aloise-Deshpande-Hansen-Popat), and for general $K$ in the plane (Mahajan-Nimbhorkar-Varadarajan). Consequently no polynomial-time algorithm solves k-means exactly unless P $=$ NP, and Lloyd's polynomial-per-iteration descent is necessarily a heuristic for the global problem. Lloyd's iteration count can itself be superpolynomial in the worst case (a $2^{Ω (n)}$ lower bound) although it is fast in practice and has polynomial smoothed complexity ^{[Arthur, D. & Vassilvitskii, S. — k-means++: The Advantages of Careful Seeding]}.

Theorem 2 (k-means++ seeding and its $O (lo g k)$ guarantee). Seed the centres by $D^{2}$ -weighting: pick $m_{1}$ uniformly at random from the data; for $j = 2, \dots, K$ pick $m_{j} = x_{i}$ with probability proportional to $D (x_{i})^{2} = min_{ℓ < j} ∥ x_{i} - m_{ℓ} ∥^{2}$ , the squared distance to the nearest already-chosen centre. Let $ϕ$ be the k-means cost of the resulting seeds (before any Lloyd step) and $ϕ_{OPT}$ the optimal cost. Then

E [ϕ] \leq 8 (ln K + 2) ϕ_{OPT},

so k-means++ is $O (lo g K)$ -competitive in expectation, a guarantee uniform seeding lacks. Running Lloyd's monotone descent from these seeds only decreases the cost further, so the bound is preserved; this seeding-then-descent is the standard modern k-means ^{[Arthur, D. & Vassilvitskii, S. — k-means++: The Advantages of Careful Seeding]}.

Theorem 3 (model-based clustering and BIC component selection). Fitting a $K$ -component Gaussian mixture by EM yields, for each covariance parameterisation (spherical, diagonal, full, tied) and each $K$ , a maximised log-likelihood $\hat{ℓ}_{K}$ and a parameter count $d_{K}$ . The Bayesian information criterion

BIC (K) = - 2 \hat{ℓ}_{K} + (lo g n) d_{K}

is minimised over $K$ (and over the covariance family) to select the number and shape of components 45.06.09. Unlike $W (C)$ , which decreases monotonically in $K$ , the likelihood-plus-penalty trades fit against parameters and admits an interior optimum; under regularity BIC is a consistent order selector. Because the free-covariance mixture likelihood is unbounded (a component mean on a single point with vanishing variance sends $\hat{ℓ} \to \infty$ ), the EM fit targets a sensible interior local maximum, the same degeneracy flagged for mixtures in 45.08.07 ^{[Hastie, T., Tibshirani, R. & Friedman, J. — The Elements of Statistical Learning (2nd ed.)]}.

Theorem 4 (agglomerative clustering and the dendrogram). Agglomerative hierarchical clustering starts with $n$ singleton clusters and repeatedly merges the two closest clusters under an inter-cluster dissimilarity $d (A, B)$ — single linkage $min_{a \in A, b \in B} d (a, b)$ , complete linkage $max$ , group-average linkage $\frac{1}{∣ A ∣∣ B ∣} \sum_{a, b} d (a, b)$ , or Ward's minimum-variance criterion (merge the pair giving the least increase in $W$ ). The merge sequence is recorded as a dendrogram, a rooted binary tree whose merge heights equal the dissimilarities at which merges occur; for the linkages above (single, complete, average, Ward) the heights are monotone non-decreasing up the tree, so the dendrogram has no inversions and a clustering at any $K$ is read off by cutting at the corresponding height. The Lance-Williams recurrence updates $d (A \cup B, \cdot)$ from $d (A, \cdot), d (B, \cdot)$ in $O (1)$ per merged pair, giving an $O (n^{2} lo g n)$ algorithm without committing to a fixed $K$ in advance ^{[Hastie, T., Tibshirani, R. & Friedman, J. — The Elements of Statistical Learning (2nd ed.)]}.

Theorem 5 (spectral clustering via the graph Laplacian). Build a similarity graph with weight matrix $W ⪰ 0$ , degree matrix $D = diag (\sum_{j} W_{ij})$ , and unnormalised graph Laplacian $L = D - W$ . For any partition indicator $f \in R^{n}$ , $f^{⊤} L f = \frac{1}{2} \sum_{i, j} W_{ij} (f_{i} - f_{j})^{2} \geq 0$ , so $L$ is positive semidefinite with constant null vector $1$ on a connected graph; the multiplicity of the eigenvalue $0$ equals the number of connected components. The balanced graph-cut objectives RatioCut and the normalised cut Ncut are NP-hard, and their continuous relaxations are solved by the bottom nonzero-eigenvalue eigenvectors of $L$ (respectively the normalised Laplacian $L_{sym} = D^{- 1/2} L D^{- 1/2}$ ); spectral clustering embeds the points in the span of those eigenvectors and runs k-means in that embedding. The eigenvector embedding linearises clusters that are not linearly separable in the original space, which is why spectral methods recover non-convex (e.g. nested-ring) structure that centroid k-means cannot. This connects to the eigen-decomposition machinery of principal components 45.06.08 ^{[Hastie, T., Tibshirani, R. & Friedman, J. — The Elements of Statistical Learning (2nd ed.)]}.

Synthesis. Clustering is one descent idea worked in several geometries, and the foundational reason the centroid method coheres is the Key theorem: alternating block-minimisation of $\sum_{i} ∥ x_{i} - m_{C (i)} ∥^{2}$ forces monotone descent to a partition-local optimum, exactly the coordinate-ascent structure that drives EM in 45.08.07, of which hard k-means is the zero-variance limit. The central insight is that the descent answers neither of the two questions a clustering poses — which local optimum (controlled by k-means++ seeding, whose $D^{2}$ -weighting is $O (lo g K)$ -competitive against the NP-hard global problem) and how many clusters (controlled by the gap statistic for k-means and by BIC for the mixture likelihood 45.06.09, which unlike $W (C)$ penalises complexity and so admits an interior optimum). This is dual to the model-based view: where k-means imposes spherical, equal-scale Voronoi cells, the Gaussian mixture with full covariances fits orientation and scale and reports soft responsibilities, so the same data can be clustered by a combinatorial objective or by a likelihood, and the two coincide in the spherical zero-variance limit.

Putting these together, hierarchical clustering abandons fixed $K$ for a whole dendrogram of nested partitions read at any cut height, and spectral clustering abandons centroid geometry for the graph Laplacian, whose bottom eigenvectors relax the NP-hard balanced-cut objectives and whose embedding — built from the same spectral theorem behind principal components 45.06.08 — linearises clusters that centroid methods cannot separate. The bridge across all four is that each replaces an intractable global combinatorial optimum with a tractable surrogate — a descent, a seeding, a greedy merge, or an eigenvector relaxation — and that none of them can be validated against a ground truth that, by the unsupervised premise, does not exist; cluster validation (stability under resampling, gap, silhouette, the dimensionless co-membership indices) is therefore an open, criterion-relative judgement rather than a test against a known answer.

Full proof set Master

Proposition 1 (centroid minimises within-cluster scatter; the bias-variance-style decomposition). For points ${x_{i}}_{i \in C_{k}}$ and any $m \in R^{p}$ , $\sum_{i \in C_{k}} ∥ x_{i} - m ∥^{2} = \sum_{i \in C_{k}} ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2} + N_{k} ∥ \overset{x}{ˉ}_{k} - m ∥^{2}$ ; hence $m = \overset{x}{ˉ}_{k}$ is the unique minimiser, with minimum the within-cluster scatter.

Proof. Insert $\pm \overset{x}{ˉ}_{k}$ : $∥ x_{i} - m ∥^{2} = ∥ (x_{i} - \overset{x}{ˉ}_{k}) + (\overset{x}{ˉ}_{k} - m) ∥^{2} = ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2} + 2 ⟨ x_{i} - \overset{x}{ˉ}_{k}, \overset{x}{ˉ}_{k} - m ⟩ + ∥ \overset{x}{ˉ}_{k} - m ∥^{2}$ . Summing over $i \in C_{k}$ , the cross term is $2 ⟨ \sum_{i} (x_{i} - \overset{x}{ˉ}_{k}), \overset{x}{ˉ}_{k} - m ⟩ = 0$ because $\sum_{i \in C_{k}} (x_{i} - \overset{x}{ˉ}_{k}) = 0$ by the definition of the centroid. The remaining terms give $\sum_{i} ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2} + N_{k} ∥ \overset{x}{ˉ}_{k} - m ∥^{2}$ . The first summand is independent of $m$ and the second is non-negative, vanishing iff $m = \overset{x}{ˉ}_{k}$ ; hence the centroid is the unique minimiser and the minimum is the within-cluster scatter. $□$

Proposition 2 (monotone descent and finite termination of Lloyd's algorithm). The Lloyd iteration satisfies $J (C^{(t + 1)}, m^{(t + 1)}) \leq J (C^{(t)}, m^{(t)})$ for all $t$ and terminates in finitely many steps at a pair satisfying the nearest-neighbour and centroid fixed-point conditions.

Proof. The assignment step minimises $J (\cdot, m^{(t)})$ over assignments: $J (C, m^{(t)}) = \sum_{i} ∥ x_{i} - m_{C (i)}^{(t)} ∥^{2}$ is a sum of independent per-point terms, each minimised by assigning $i$ to its nearest representative, so $J (C^{(t + 1)}, m^{(t)}) \leq J (C^{(t)}, m^{(t)})$ . The update step minimises $J (C^{(t + 1)}, \cdot)$ over representatives: by Proposition 1 the per-cluster minimiser is the centroid $m_{k}^{(t + 1)} = \overset{x}{ˉ}_{k}^{(t + 1)}$ , so $J (C^{(t + 1)}, m^{(t + 1)}) \leq J (C^{(t + 1)}, m^{(t)})$ . Chaining gives the monotone decrease. The realised values $J (C^{(t + 1)}, m^{(t + 1)}) = W (C^{(t + 1)})$ are determined by the assignment, of which there are at most $K^{n}$ ; a non-increasing sequence taking finitely many values and strictly decreasing whenever the assignment changes must, after finitely many steps, leave the assignment fixed. At that fixed assignment $C^{⋆}$ , every point is assigned to its nearest current centroid (else the assignment step would change it) and every $m_{k}^{⋆} = \overset{x}{ˉ}_{k}^{⋆}$ (the update fixed point), the two stated conditions. $□$

Proposition 3 (hard k-means as the $σ^{2} \to 0$ limit of spherical-mixture EM). For the equal-weight spherical Gaussian mixture $Σ_{k} = σ^{2} I$ , as $σ^{2} \to 0^{+}$ the EM responsibilities $γ_{ik}$ converge to the nearest-centroid indicators $1 {k = ar g min_{j} ∥ x_{i} - μ_{j} ∥^{2}}$ (for a unique nearest centre), and the EM M-step coincides with the Lloyd update.

Proof. With $d_{ij} = ∥ x_{i} - μ_{j} ∥^{2}$ and equal weights, $γ_{ik} = exp (- d_{ik} /2 σ^{2}) / \sum_{j} exp (- d_{ij} /2 σ^{2})$ . Let $k_{i}^{⋆} = ar g min_{j} d_{ij}$ , unique by assumption, and write $γ_{ik} = 1/ \sum_{j} exp ((d_{ik} - d_{ij}) /2 σ^{2})$ . For $k = k_{i}^{⋆}$ , every exponent $(d_{i k_{i}^{⋆}} - d_{ij}) /2 σ^{2} \leq 0$ with the $j = k_{i}^{⋆}$ term equal to $1$ and all others $\to 0$ as $σ^{2} \to 0^{+}$ ; hence $γ_{i k_{i}^{⋆}} \to 1$ . For $k \neq = k_{i}^{⋆}$ , the $j = k_{i}^{⋆}$ denominator term is $exp ((d_{ik} - d_{i k_{i}^{⋆}}) /2 σ^{2}) \to + \infty$ since $d_{ik} > d_{i k_{i}^{⋆}}$ , so $γ_{ik} \to 0$ . The responsibilities become hard indicators, and the M-step mean $μ_{k} = \sum_{i} γ_{ik} x_{i} / \sum_{i} γ_{ik}$ becomes the average of the points whose nearest centre is $k$ — the Lloyd update. The EM monotonicity of 45.08.07 specialises in this limit to the Lloyd descent of Proposition 2. $□$

Proposition 4 (positive semidefiniteness of the graph Laplacian and the quadratic-form identity). For a symmetric non-negative weight matrix $W$ with degree matrix $D = diag (d_{i})$ , $d_{i} = \sum_{j} W_{ij}$ , the Laplacian $L = D - W$ satisfies $f^{⊤} L f = \frac{1}{2} \sum_{i, j} W_{ij} (f_{i} - f_{j})^{2}$ for all $f \in R^{n}$ ; consequently $L ⪰ 0$ , $L 1 = 0$ , and on a connected graph $dim ker L = 1$ .

Proof. Expand $\frac{1}{2} \sum_{i, j} W_{ij} (f_{i} - f_{j})^{2} = \frac{1}{2} \sum_{i, j} W_{ij} (f_{i}^{2} - 2 f_{i} f_{j} + f_{j}^{2})$ . The $f_{i}^{2}$ and $f_{j}^{2}$ terms each give $\frac{1}{2} \sum_{i} d_{i} f_{i}^{2}$ , summing to $\sum_{i} d_{i} f_{i}^{2} = f^{⊤} D f$ ; the cross term gives $- \sum_{i, j} W_{ij} f_{i} f_{j} = - f^{⊤} W f$ . Hence $\frac{1}{2} \sum_{i, j} W_{ij} (f_{i} - f_{j})^{2} = f^{⊤} D f - f^{⊤} W f = f^{⊤} L f$ . Since $W_{ij} \geq 0$ , the right side is non-negative, so $L ⪰ 0$ . Taking $f = 1$ makes every $(f_{i} - f_{j}) = 0$ , so $1 \in ker L$ . If $f^{⊤} L f = 0$ then $W_{ij} (f_{i} - f_{j})^{2} = 0$ for all $i, j$ , so $f$ is constant across every edge; on a connected graph this forces $f$ constant, so $ker L = span {1}$ is one-dimensional. (On a graph with $c$ components the same argument gives $dim ker L = c$ , the indicator vectors of the components spanning the kernel.) $□$

Proposition 5 (within/between decomposition: k-means maximises between-cluster scatter). Let $T = \sum_{i} ∥ x_{i} - \overset{x}{ˉ} ∥^{2}$ be the total scatter about the grand mean $\overset{x}{ˉ}$ . For any clustering, $T = W (C) + B (C)$ with $B (C) = \sum_{k} N_{k} ∥ \overset{x}{ˉ}_{k} - \overset{x}{ˉ} ∥^{2}$ the between-cluster scatter; since $T$ is fixed, minimising $W (C)$ is equivalent to maximising $B (C)$ .

Proof. Insert $\pm \overset{x}{ˉ}_{k}$ into each total-scatter term: for $i \in C_{k}$ , $∥ x_{i} - \overset{x}{ˉ} ∥^{2} = ∥ (x_{i} - \overset{x}{ˉ}_{k}) + (\overset{x}{ˉ}_{k} - \overset{x}{ˉ}) ∥^{2} = ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2} + 2 ⟨ x_{i} - \overset{x}{ˉ}_{k}, \overset{x}{ˉ}_{k} - \overset{x}{ˉ} ⟩ + ∥ \overset{x}{ˉ}_{k} - \overset{x}{ˉ} ∥^{2}$ . Summing over $i \in C_{k}$ , the cross term vanishes because $\sum_{i \in C_{k}} (x_{i} - \overset{x}{ˉ}_{k}) = 0$ , leaving $\sum_{i \in C_{k}} ∥ x_{i} - \overset{x}{ˉ}_{k} ∥^{2} + N_{k} ∥ \overset{x}{ˉ}_{k} - \overset{x}{ˉ} ∥^{2}$ . Summing over $k$ gives $T = W (C) + B (C)$ . The total scatter $T$ depends only on the data, not on the clustering, so over partitions $min_{C} W (C) = min_{C} (T - B (C)) = T - max_{C} B (C)$ ; the minimiser of within-cluster scatter is the maximiser of between-cluster scatter. $□$

Connections Master

k-means is the hard-assignment, zero-variance limit of the EM algorithm for a spherical Gaussian mixture developed in 45.08.07: the assignment step is the degenerate E-step whose responsibilities collapse to nearest-centroid indicators (Proposition 3), and the centroid update is the M-step. The monotone-descent guarantee of the Key theorem is the specialisation of EM's monotone ascent, and the foundational reason both algorithms converge is the same — alternating block-minimisation of a single objective — so the soft mixture model and hard k-means are two ends of one temperature scale.
Choosing the number of clusters reuses the model-assessment and selection machinery of 45.06.09: BIC $= - 2 \hat{ℓ}_{K} + (lo g n) d_{K}$ selects the number of Gaussian-mixture components by trading likelihood against parameter count, exactly the complexity-penalised selection studied there, and the gap statistic plays the analogous role for the centroid objective $W (C)$ , which on its own decreases monotonically in $K$ and so cannot select it.
Spectral clustering embeds points using the bottom eigenvectors of the graph Laplacian $L = D - W$ , built from the same symmetric-matrix spectral theorem that underlies principal components in 45.06.08: where PCA takes the top eigenvectors of a covariance (directions of maximal variance), spectral clustering takes the bottom nonzero-eigenvalue eigenvectors of the Laplacian (the smoothest functions on the graph), and both are eigen-embeddings that linearise structure before a downstream method (regression there, k-means here) is applied.
The soft-assignment / hard-assignment duality that relates k-means to its mixture model recurs as the relationship between the discriminant boundaries of 45.08.03 and the generative Gaussian model behind them: model-based clustering is the unsupervised counterpart of Gaussian discriminant analysis, fitting the same per-class Gaussians but inferring the class labels rather than observing them, so the central insight that a generative Gaussian model induces a geometric (quadratic or linear) decision rule appears again with the labels treated as latent.

Historical & philosophical context Master

The centroid algorithm has several independent origins. Stuart Lloyd described the alternating nearest-neighbour / centroid procedure for least-squares scalar quantization in a 1957 Bell Laboratories technical note, published only in 1982 in the IEEE Transactions on Information Theory ^{[Lloyd, S. P. — Least Squares Quantization in PCM]}; the name k-means is due to James MacQueen (1967), and Hugo Steinhaus (1956) and Edward Forgy (1965) gave closely related formulations, so the method is genuinely a case of multiple discovery across quantization, statistics, and pattern recognition. The vector-quantization generalisation of Linde, Buzo, and Gray (1980) is identical in form and is the version used in signal coding.

The computational status of the global problem was settled much later: Aloise, Deshpande, Hansen, and Popat (2009) proved that Euclidean sum-of-squares clustering is NP-hard already for two clusters in general dimension, and Mahajan, Nimbhorkar, and Varadarajan (2009) proved NP-hardness for general $K$ in the plane, so Lloyd's polynomial-per-iteration descent is necessarily a heuristic for an intractable optimum. The practical response was the $D^{2}$ -seeding of David Arthur and Sergei Vassilvitskii (2007), whose k-means++ gives the $O (lo g K)$ -competitive expected guarantee that uniform seeding lacks ^{[Arthur, D. & Vassilvitskii, S. — k-means++: The Advantages of Careful Seeding]}. The model-based programme — clustering as maximum-likelihood fitting of a finite mixture, with covariance constraints and BIC selection — was developed by Banfield and Raftery (1993) and Fraley and Raftery (2002) in the mclust framework, and the spectral approach via the graph Laplacian was synthesised by Shi and Malik (2000) (normalised cuts), Ng, Jordan, and Weiss (2002), and the survey of von Luxburg (2007). The recurring methodological difficulty, emphasised throughout Hastie, Tibshirani, and Friedman's treatment ^{[Hastie, T., Tibshirani, R. & Friedman, J. — The Elements of Statistical Learning (2nd ed.)]}, is that unsupervised clustering has no ground-truth label against which to validate, so the choice of objective, dissimilarity, and number of clusters is criterion-relative and the resulting partition is a hypothesis rather than a verified answer.

Bibliography Master

@article{lloyd1982quantization,
  author  = {Lloyd, Stuart P.},
  title   = {Least squares quantization in {PCM}},
  journal = {IEEE Transactions on Information Theory},
  volume  = {28},
  number  = {2},
  pages   = {129--137},
  year    = {1982},
  note    = {Bell Laboratories technical note, 1957}
}

@inproceedings{arthur2007kmeanspp,
  author    = {Arthur, David and Vassilvitskii, Sergei},
  title     = {k-means++: The Advantages of Careful Seeding},
  booktitle = {Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA)},
  pages     = {1027--1035},
  year      = {2007}
}

@article{aloise2009nphard,
  author  = {Aloise, Daniel and Deshpande, Amit and Hansen, Pierre and Popat, Preyas},
  title   = {NP-hardness of {Euclidean} sum-of-squares clustering},
  journal = {Machine Learning},
  volume  = {75},
  number  = {2},
  pages   = {245--248},
  year    = {2009}
}

@article{macqueen1967kmeans,
  author  = {MacQueen, James},
  title   = {Some methods for classification and analysis of multivariate observations},
  journal = {Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability},
  volume  = {1},
  pages   = {281--297},
  year    = {1967}
}

@article{lindebuzogray1980,
  author  = {Linde, Yoseph and Buzo, Andr{\'e}s and Gray, Robert M.},
  title   = {An algorithm for vector quantizer design},
  journal = {IEEE Transactions on Communications},
  volume  = {28},
  number  = {1},
  pages   = {84--95},
  year    = {1980}
}

@article{fraleyraftery2002mclust,
  author  = {Fraley, Chris and Raftery, Adrian E.},
  title   = {Model-based clustering, discriminant analysis, and density estimation},
  journal = {Journal of the American Statistical Association},
  volume  = {97},
  number  = {458},
  pages   = {611--631},
  year    = {2002}
}

@article{vonluxburg2007spectral,
  author  = {von Luxburg, Ulrike},
  title   = {A tutorial on spectral clustering},
  journal = {Statistics and Computing},
  volume  = {17},
  number  = {4},
  pages   = {395--416},
  year    = {2007}
}

@article{shimalik2000ncut,
  author  = {Shi, Jianbo and Malik, Jitendra},
  title   = {Normalized cuts and image segmentation},
  journal = {IEEE Transactions on Pattern Analysis and Machine Intelligence},
  volume  = {22},
  number  = {8},
  pages   = {888--905},
  year    = {2000}
}

@book{hastie2009esl,
  author    = {Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome},
  title     = {The Elements of Statistical Learning: Data Mining, Inference, and Prediction},
  edition   = {2},
  publisher = {Springer},
  year      = {2009}
}

Prerequisites

45.08.07

Tier anchors

beginner: Hastie, Tibshirani & Friedman 2009 The Elements of Statistical Learning 2e (Springer) §14.3.6 (k-means as the alternating assign-to-nearest-centroid / recompute-centroid loop, with the human-tumour and simulated examples); MacKay 2003 Information Theory, Inference, and Learning Algorithms (Cambridge) §20 (the K-means algorithm presented with the soft-assignment / responsibility picture and the worked two-cluster runs)
intermediate: Hastie, Tibshirani & Friedman 2009 The Elements of Statistical Learning 2e (Springer) §14.3.6-14.3.9 (k-means, the within-cluster point-scatter objective W(C), the two-stage descent on (assignment, centroids), k-medoids, choosing K via the gap statistic, and the soft k-means / mixture connection); Lloyd 1982 IEEE Trans. Inform. Theory 28:129-137 (least-squares quantization in PCM: the centroid/nearest-neighbour fixed-point conditions and the descent algorithm)
master: Hastie, Tibshirani & Friedman 2009 The Elements of Statistical Learning 2e (Springer) §14.3.6-14.3.12 (model-based clustering with Gaussian mixtures and BIC, hierarchical/agglomerative clustering and the dendrogram, spectral clustering via the graph Laplacian); Arthur & Vassilvitskii 2007 Proc. SODA pp. 1027-1035 (k-means++: the $D^{2}$ seeding and its $O (lo g k)$ -competitive expected approximation guarantee); Aloise, Deshpande, Hansen & Popat 2009 Machine Learning 75:245-248 (NP-hardness of Euclidean sum-of-squares clustering for general dimension)

References

Hastie, T., Tibshirani, R. & Friedman, J. — The Elements of Statistical Learning (2nd ed.) · Springer Series in Statistics (2009). Chapter 14 'Unsupervised Learning', §14.3 'Cluster Analysis'. §14.3.6 'K-means' defines the within-cluster point-scatter / pooled within-cluster sum of squares W(C) = (1/2) sum_k sum_{C(i)=k} sum_{C(i')=k} ||x_i - x_i'||^2 = sum_k N_k sum_{C(i)=k} ||x_i - bar x_k||^2, shows that for squared Euclidean distance the optimal cluster representative is the cluster mean, and gives the two-step descent: with centroids fixed, assign each point to the nearest centroid; with the assignment fixed, recompute each centroid as the cluster mean. The text states that each step does not increase W(C) and the algorithm converges to a local optimum, and that the result depends on the random starts so multiple initialisations are used. §14.3.7 covers Gaussian mixtures as soft K-means (the responsibilities, the temperature parameter, and the hard limit). §14.3.8 'Example: Human Tumor Microarray Data' and §14.3.9 'Vector Quantization' give applications. §14.3.9-14.3.11 cover K-medoids (PAM, using a medoid and a general dissimilarity), the practical issue of choosing K via the gap statistic Gap(k) = E*[log W_k] - log W_k against a uniform reference distribution, and the silhouette. §14.3.12 'Hierarchical Clustering' develops agglomerative clustering with single, complete, and group-average (and Ward) linkage and the dendrogram as a representation of the nested merge sequence with monotone (cophenetic) heights, and divisive clustering. §14.5.3 'Spectral Clustering' constructs the graph from a similarity matrix, the graph Laplacian L = D - W (and the normalised versions), the eigenvector embedding using the smallest nonzero-eigenvalue eigenvectors of L, and the link to graph-cut (RatioCut / Ncut) relaxations; §8.5.1 and §6.8 relate the mixture-model and BIC machinery (BIC = -2 loglik + (log N) d) used to select the number of mixture components. The chapter also discusses the absence of a ground truth and the difficulty of validating a clustering.
Lloyd, S. P. — Least Squares Quantization in PCM · IEEE Transactions on Information Theory 28(2):129-137 (1982); the work circulated as an unpublished Bell Laboratories technical note from 1957. Treats the design of a scalar quantizer minimising mean squared error and derives the two necessary conditions for an optimal quantizer: (i) the nearest-neighbour / minimum-distortion partition condition — each point is encoded to the representative (quantization level) nearest it; and (ii) the centroid condition — each representative is the conditional mean (centroid) of the region assigned to it. The 'Lloyd algorithm' (Lloyd's Method I) alternates the two conditions and is the continuous / population ancestor of the finite-sample k-means iteration; the paper establishes the fixed-point characterisation that the alternating algorithm targets and the monotone decrease of the distortion. The vector / k-level generalisation (the Linde-Buzo-Gray generalised Lloyd algorithm, 1980) is the vector-quantization form used in signal coding and is identical in form to k-means.
Arthur, D. & Vassilvitskii, S. — k-means++: The Advantages of Careful Seeding · Proceedings of the eighteenth annual ACM-SIAM Symposium on Discrete Algorithms (SODA '07), pp. 1027-1035 (2007). Introduces the k-means++ seeding (D^2 weighting): choose the first centre uniformly at random from the data, then choose each subsequent centre x' with probability proportional to D(x')^2, the squared distance from x' to the nearest centre already chosen, and run Lloyd's algorithm from these seeds. The main theorem proves that the expected k-means cost of the resulting clustering is O(log k) times the optimal cost, E[phi] <= 8(ln k + 2) phi_OPT, so the seeding is O(log k)-competitive in expectation before any Lloyd iteration — a guarantee plain uniform/Forgy seeding does not enjoy. The paper also reports large empirical speed and quality improvements. (Background: Aloise, Deshpande, Hansen & Popat 2009, Machine Learning 75:245-248, prove that minimising the Euclidean sum-of-squares clustering objective is NP-hard even for k=2 clusters in general dimension; Mahajan, Nimbhorkar & Varadarajan 2009 prove NP-hardness for general k in the plane; Dasgupta 2008 for general dimension with k part of the input — so the global k-means optimum is intractable and the descent/seeding heuristics are the practical response.)

Estimated time

beginner: 20m
intermediate: 50m
master: 90m