k-Means Clustering

k-Means clustering is a foundational unsupervised learning, data partitioning, and vector quantization algorithm in computation and algorithms, multivariate statistics, and signal processing. Given a dataset of n observations X = {x₁, x₂, …, x_n} in d-dimensional Euclidean space ℝ^d and a user-specified integer k ≤ n, k-means partitions the data into k disjoint clusters S = {S₁, S₂, …, S_k} with cluster centers (centroids) μ₁, …, μ_k ∈ ℝ^d so as to minimize the within-cluster sum of squares (WCSS), also known as the distortion or inertia objective.

The mathematical problem and its standard alternating minimization heuristic were discovered independently across four disciplines: by Hugo Steinhaus in 1956 (spatial mechanics), Stuart P. Lloyd at Bell Labs in 1957 (pulse-code modulation quantization, published in 1982), Edward W. Forgy in 1965, and James MacQueen in 1967, who coined the term "k-means." Geometrically, k-means tessellates ℝ^d into a centroidal Voronoi diagram whose cells are convex polyhedra, contrasting with density-based clustering methods such as DBSCAN, which discover arbitrarily shaped clusters and reject background noise without fixing k in advance.

Given X ⊂ ℝ^d, the k-means optimization problem seeks a partition S = {S₁, …, S_k} and centroid vectors μ₁, …, μ_k ∈ ℝ^d minimizing the distortion function:

J(S, μ) = ∑_j=1^k ∑_xi∈Sj ‖x_i − μ_j‖₂².

By the variance decomposition identity, minimizing the within-cluster sum of squares J(S, μ) is algebraically equivalent to maximizing the between-cluster weighted sum of squared distances ∑_j=1^k |S_j| ‖μ_j − x̄‖₂², and to a relaxed spectral trace maximization linked to principal component analysis (Ding & He, 2004). Despite its simple quadratic form, finding the global minimum of J(S, μ) is NP-hard even for k = 2 clusters in general dimension d (Aloise et al., 2009) and for arbitrary k in the Euclidean plane d = 2 (Mahajan, Nimbhorkar, & Varadarajan, 2012).

Lloyd's Two-Step Coordinate Descent Iteration

Lloyd's algorithm finds a local minimum of J(S, μ) by alternating between two coordinate-descent steps from an initial set of centroids {μ₁⁽⁰⁾, …, μ_k⁽⁰⁾}:

1. Assignment (Voronoi) step: Holding the centroids μ₁^(t), …, μ_k^(t) fixed, assign each data point x_i to the cluster whose centroid is nearest in Euclidean distance (breaking ties deterministically by lowest index):

S_j^(t) = { x_i ∈ X : ‖x_i − μ_j^(t)‖₂² ≤ ‖x_i − μ_l^(t)‖₂² for all 1 ≤ l ≤ k }.

1. Update (Centroid) step: Holding the cluster assignments S₁^(t), …, S_k^(t) fixed, differentiate ∑_xi∈Sj ‖x_i − μ_j‖₂² with respect to μ_j to obtain the unique minimizer—the arithmetic mean (center of mass) of the points in S_j^(t):

μ_j^(t+1) = (1 / |S_j^(t)|) ∑_xi∈Sj^(t) x_i.

Because both steps monotonically decrease the non-negative objective J(S, μ) and there are at most kⁿ finite partitions of n points, Lloyd's iteration always terminates in a finite number of steps at a stationary centroidal Voronoi configuration. Each iteration requires O(nkd) floating-point operations.

When initial centroids are chosen uniformly at random from X (Forgy initialization), two centroids frequently land inside the same true cluster while another cluster receives none, trapping Lloyd's algorithm in arbitrarily poor local minima. In 2007, David Arthur and Sergei Vassilvitskii introduced k-means++ initialization, which spreads out initial centers via randomized D²-weighting:

1. Choose the first center μ₁ uniformly at random from X.
2. For each subsequent center m = 2, …, k, compute for every point x_i ∈ X its squared distance D(x_i)² = min_1≤j<m ‖x_i − μ_j‖₂² to the closest center already chosen, and sample the next center μ_m = x_i with probability proportional to D(x_i)² / ∑_r=1ⁿ D(x_r)².

Arthur and Vassilvitskii proved that k-means++ initialization alone guarantees an expected distortion 𝔼[J] ≤ 8(ln k + 2) J_opt before running a single Lloyd step, while accelerating subsequent Lloyd convergence by two- to tenfold.

Elkan's Triangle-Inequality Acceleration and Mini-Batch k-Means

In later iterations of Lloyd's algorithm, centroids move very little, so most point-to-centroid distance calculations ‖x_i − μ_j‖₂ are redundant. Elkan's algorithm (2003) applies the metric triangle inequality ‖x_i − μ_j‖₂ ≥ ‖x_i − c(x_i)‖₂ whenever ‖c(x_i) − μ_j‖₂ ≥ 2 ‖x_i − c(x_i)‖₂ to skip up to 90% of distance evaluations while returning the same partition as standard Lloyd's algorithm. For web-scale datasets where n > 10⁸, Mini-Batch k-Means (Sculley, 2010) samples random mini-batches of size b (e.g., b = 1,024) via a Monte Carlo method and applies stochastic gradient descent updates to the centroids.

Because the raw distortion J(S, μ) decreases monotonically as k increases (reaching J = 0 when k = n), selecting an appropriate cluster count k requires external criteria such as the Silhouette coefficient (Rousseeuw, 1987), which compares each point's mean intra-cluster distance a(i) against its nearest neighbor-cluster distance b(i) via s(i) = [b(i) − a(i)] / max(a(i), b(i)) ∈ [−1, +1]; the Calinski–Harabasz variance-ratio criterion; or Tibshirani's Gap statistic, which compares log J_k against a Monte Carlo null reference distribution.

Probabilistically, k-means is the small-variance zero-temperature limit (σ² → 0⁺) of the Expectation–Maximization (EM) algorithm for an isotropic Gaussian Mixture Model (GMM) ∑_j=1^k π_j N(x | μ_j, σ²I) with equal spherical covariances and equal mixing weights, where soft posterior responsibilities collapse to hard 0–1 cluster assignments.

When does k-means clustering fail?

Because k-means uses isotropic Euclidean distance from cluster means, it fails in three common scenarios: (1) non-convex or intertwined manifolds (such as concentric rings or crescent moons, where DBSCAN or spectral clustering succeeds); (2) elongated anisotropic elliptical clusters (where full-covariance Gaussian Mixture Models succeed); and (3) datasets with strong outliers, where squared Euclidean loss pulls centroids away from the cluster core (where k-medians or k-medoids succeeds).

Why should features be standardized before running k-means?

Because ‖x_i − μ_j‖₂² = ∑_m=1^d (x_i,m − μ_j,m)² sums squared coordinate differences across all dimensions, a single variable measured in large numerical units (such as income in dollars versus age in years) dominates the Euclidean distance unless features are standardized to zero mean and unit variance (z-scores).

1. Lloyd, S. P. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2), 129–137. https://doi.org/10.1109/TIT.1982.1056489
2. MacQueen, J. (1967). Some methods for classification and analysis of multivariate observations. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, 1, 281–297. https://doi.org/10.1525/9780520313880-018
3. Arthur, D., & Vassilvitskii, S. (2007). k-means++: The advantages of careful seeding. Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, 1027–1035. https://doi.org/10.1145/1283383.1283494
4. Elkan, C. (2003). Using the triangle inequality to accelerate k-means. Proceedings of the Twentieth International Conference on Machine Learning (ICML), 147–153. https://doi.org/10.5555/3041838.3041857