K-Means Clustering

Imagine a retailer with ten million customer transaction histories. It has no labels telling it who is a "bargain hunter" versus a "premium loyalist." It still needs to group similar customers to target marketing campaigns. This is the central problem of The Unsupervised Learning Problem: discover structure in data when no ground truth exists. K-means clustering is the most widely taught and most widely deployed answer to this problem. It is simple to state, cheap to compute, and surprisingly effective. It works well on data that forms roughly round, well-separated blobs. K-means partitions a dataset of points in Euclidean space into $K$ groups. Points within a group should be close to each other. They should also be far from points in other groups. It does this by alternating between two simple steps. First, assign each point to its nearest cluster center. Then, move each center to the average of the points assigned to it. This alternation, known as Lloyd's algorithm, is the workhorse implementation. It is almost always what people mean when they say "k-means." The idea matters because it appears everywhere real systems need fast, interpretable grouping. Examples include compressing images by quantizing colors, and initializing more sophisticated probabilistic models. It also helps in building codebooks for vector quantization in speech and video compression. It is used for segmenting customers as well. It even serves as a subroutine inside more advanced pipelines, like feature learning and anomaly detection. Its appeal combines computational cheapness with conceptual transparency. It runs near-linear in the number of points per iteration. Every cluster is summarized by a single interpretable vector, its centroid. This article develops the within-cluster sum-of-squares objective that k-means minimizes. It derives Lloyd's algorithm from that objective. It proves why each step of the algorithm can only decrease the objective. It is also honest about the algorithm's central weakness. K-means converges only to a local optimum. The optimum it finds depends heavily on where you start. This last point sets up the natural next topics. These are K-Means++ Initialization and Choosing the Number of Clusters. It also leads to the probabilistic generalization known as Gaussian Mixture Models.

Key takeaways

K-means partitions $n$ points into $K$ clusters by minimizing the within-cluster sum of squared distances to cluster centroids.
Lloyd's algorithm alternates an assignment step (nearest centroid, by Distance Metrics and Similarity Measures) and an update step (recompute centroid as the mean).
Each Lloyd iteration provably never increases the objective, so the algorithm converges to a local minimum in finitely many steps.
Per-iteration cost is $O(nKd)$ time and $O(nd + Kd)$ space, making it scalable to large $n$ .
The final clustering depends strongly on initialization; poor starts give poor local optima — motivating K-Means++ Initialization.
K-means implicitly assumes roughly spherical, similarly-sized clusters, which is its main modeling limitation.

The core insight

K-means alternately assigns points to their nearest centroid and recomputes centroids as cluster means, and each of these two moves can only decrease the same sum-of-squared-distances objective, guaranteeing convergence to a local minimum.

The idea

Given a dataset $X = \{x_1, \dots, x_n\}$ with $x_i \in \mathbb{R}^d$ , and a chosen number of clusters $K$ , k-means seeks a partition of the points into sets $C_1, \dots, C_K$ and a set of centroids $\mu_1, \dots, \mu_K \in \mathbb{R}^d$ that minimize the within-cluster sum of squares (WCSS), also called inertia or the distortion:

J(C, \mu) = \sum_{k=1}^{K} \sum_{x_i \in C_k} \lVert x_i - \mu_k \rVert_2^2 .

This is a joint optimization over two coupled unknowns: a discrete assignment of every point to one of $K$ labels, and $K$ continuous centroid vectors. Solving it exactly is NP-hard even for $K=2$ in general dimension, because the number of possible partitions of $n$ points into $K$ groups grows combinatorially (a Stirling number of the second kind). K-means clustering, in practice, refers to a heuristic — Lloyd's algorithm — that finds a local minimum of $J$ efficiently rather than searching over all partitions. The objective $J$ has an important algebraic reading: for a fixed set of points assigned to a cluster, the point that minimizes the sum of squared Euclidean distances to all of them is exactly their mean. This is a direct consequence of the Distance Metrics and Similarity Measures used — squared Euclidean distance decomposes cleanly around the mean — and it is the reason the "update" half of the algorithm is simply an averaging operation rather than a general-purpose optimization.

How it works

Lloyd's algorithm is an iterative two-step procedure. It requires an initial guess for the $K$ centroids (commonly $K$ points chosen at random from the dataset, though smarter schemes exist), and it alternates: 1. Assignment step. Assign every point to the cluster whose centroid is nearest, using squared Euclidean distance. 2. Update step. Recompute each centroid as the mean of the points currently assigned to it. These two steps repeat until assignments stop changing (or change negligibly), at which point the algorithm has converged to a local minimum. ``` function LLOYDS-KMEANS(X, K, max_iters, tol):

# X: array of n points in R^d n, d = shape(X)

# 1. Initialize K centroids (e.g., K distinct points sampled from X) centroids = INITIALIZE(X, K) # shape (K, d) assignments = array of size n, filled with -1 for iter in 1 .. max_iters:

# --- Assignment step --- new_assignments = array of size n for i in 1 .. n: best_k = argmin over k in 1..K of squared_distance(X[i], centroids[k]) new_assignments[i] = best_k

# --- Update step --- new_centroids = array of shape (K, d), filled with 0 counts = array of size K, filled with 0 for i in 1 .. n: k = new_assignments[i] new_centroids[k] += X[i] counts[k] += 1 for k in 1 .. K: if counts[k] == 0:

# Empty cluster: re-seed, e.g. to the farthest point from

# its current centroid, to avoid a degenerate cluster new_centroids[k] = REINITIALIZE_EMPTY_CLUSTER(X, centroids) else: new_centroids[k] = new_centroids[k] / counts[k]

# --- Convergence check --- shift = sum over k of squared_distance(new_centroids[k], centroids[k]) centroids = new_centroids if new_assignments == assignments or shift < tol: assignments = new_assignments break assignments = new_assignments return assignments, centroids ``` Two implementation details matter for correctness. First, the empty cluster problem: if no point is assigned to some centroid $k$ , the mean is undefined, so an empty cluster must be explicitly re-seeded (common choices: the point farthest from its current centroid, or a random point). Second, the stopping rule: implementations may stop on "assignments unchanged" (exact convergence, guaranteed to occur since there are finitely many partitions), or on a small centroid shift / objective decrease (an approximate, faster-terminating criterion used in practice with a `max_iters` cap as a safety net). The figure below shows the state of a small 2D dataset through one full assignment–update cycle, using $K=3$ .

Figure

One Lloyd iteration: points colored by nearest current centroid, then centroids re-centered on their assigned points

Complexity & correctness

Time complexity. Each iteration performs the assignment step, which computes the distance from every one of $n$ points to every one of $K$ centroids in $d$ dimensions: $O(nKd)$ . The update step sums points into $K$ bins, also $O(nd)$ , dominated by the assignment step. If the algorithm runs for $T$ iterations, total time is

O(nKdT).

$T$ is not bounded by a nice closed-form function of $n$ in the worst case — in fact, pathological inputs exist on which Lloyd's algorithm requires a super-polynomial number of iterations — but in practice $T$ is small (often under a few dozen) and the algorithm is treated as near-linear in $n$ for fixed $K$ and $d$ .

Space complexity. Storing the dataset needs $O(nd)$ , the centroids need $O(Kd)$ , and the assignment vector needs $O(n)$ , giving $O(nd + Kd)$ total — dominated by the data itself, so k-means adds negligible memory overhead beyond holding the dataset.

Correctness: why Lloyd's algorithm converges. The key fact is that both steps of the algorithm are individually optimal descent moves on the same objective $J$ , when the other variable is held fixed. Assignment step is optimal given fixed centroids. For a fixed set of centroids $\mu_1, \dots, \mu_K$ , the objective decomposes point-by-point:

J = \sum_{i=1}^n \lVert x_i - \mu_{c(i)} \rVert^2,

where $c(i)$ is the cluster assigned to $x_i$ . Each term depends only on that point's own assignment, so $J$ is minimized by choosing, independently for each $i$ , the centroid $\mu_k$ that minimizes $\lVert x_i - \mu_k \rVert^2$ — exactly the nearest-centroid rule. No other assignment can produce a smaller $J$ for these centroids. Update step is optimal given fixed assignments. For a fixed partition, $J$ splits into $K$ independent sums, one per cluster:

J = \sum_{k=1}^K \sum_{x_i \in C_k} \lVert x_i - \mu_k \rVert^2 .

For a fixed set of points $C_k$ , the function $f(\mu) = \sum_{x_i \in C_k} \lVert x_i - \mu \rVert^2$ is a strictly convex quadratic in $\mu$ . Differentiating and setting the gradient to zero,

\nabla_\mu f(\mu) = \sum_{x_i \in C_k} -2(x_i - \mu) = 0 \quad \Longrightarrow \quad \mu = \frac{1}{|C_k|}\sum_{x_i \in C_k} x_i,

so the unique minimizer is the sample mean. This is why the update step recomputes centroids as arithmetic means — it is literally solving the exact minimization for that half of the problem. Putting it together. Let $J_t$ be the objective value after the $t$ -th assignment step and update step. Because the assignment step strictly minimizes $J$ over all assignments given the previous centroids, and the update step strictly minimizes $J$ over all centroids given the new assignments,

J_0 \ge J_1 \ge J_2 \ge \cdots

is a monotonically non-increasing sequence. Since there are only finitely many ways to partition $n$ points into $K$ non-empty groups (a finite number bounded by $K^n$ ), and each distinct partition has a fixed, well-defined objective value, the sequence of partitions visited cannot decrease $J$ forever without repeating — and since $J$ strictly decreases whenever the partition changes (excluding pathological ties), the algorithm must reach a partition that repeats, i.e., converge, in a finite number of steps. In practice this bound is astronomically loose; convergence is typically fast because the descent is steep early on.

What convergence does not guarantee. Lloyd's algorithm converges to a local minimum of $J$ , not the global minimum. Because $J$ is a non-convex function of the joint (assignment, centroid) variables, different initializations can lead to different local optima with different final objective values. This single fact is the central limitation of k-means and is the reason smarter seeding strategies — see K-Means++ Initialization — exist: they don't change what each iteration does, only where the search starts, but that alone substantially improves the typical quality of the local optimum reached.

Modeling assumptions and honest tradeoffs. K-means implicitly assumes clusters are:

roughly spherical (isotropic) in the chosen feature space, because squared Euclidean distance treats all directions equally,
of comparable size and density, because the mean is dominated by whichever cluster has more mass,
separable by a Voronoi-like partition, since the assignment step always produces convex polyhedral regions around centroids.

When clusters are elongated, of very different sizes, or non-convex (e.g., two concentric rings), k-means will confidently produce a wrong-looking but locally optimal answer. This limitation motivates the probabilistic generalization in Gaussian Mixture Models, which allows elliptical, differently-shaped clusters via a full covariance matrix per component, at the cost of a more expensive EM-based fit.

Worked example

Consider six 1D points (for arithmetic that is easy to check by hand, using scalar distances):

X = \{1, 2, 4, 5, 8, 9\}, \qquad K = 2.

Initialize centroids at $\mu_1 = 1$ , $\mu_2 = 9$ (the two extreme points — a deliberately "lucky" start for intuition first, then we will also examine a bad start).

Iteration 1 — assignment. Distance to $\mu_1=1$ vs. $\mu_2=9$ for each point:

Point	$\|x-\mu_1\|$	$\|x-\mu_2\|$	Assigned to
1	0	8	cluster 1
2	1	7	cluster 1
4	3	5	cluster 1
5	4	4	cluster 1 (tie broken to first, cluster 1)
8	7	1	cluster 2
9	8	0	cluster 2

Iteration 1 — update. Cluster 1 = $\{1,2,4,5\}$ , mean $= 12/4 = 3$ . Cluster 2 = $\{8,9\}$ , mean $= 17/2 = 8.5$ . New centroids: $\mu_1 = 3$ , $\mu_2 = 8.5$ .

Iteration 2 — assignment.

Point	$\|x-3\|$	$\|x-8.5\|$	Assigned to
1	2	7.5	cluster 1
2	1	6.5	cluster 1
4	1	4.5	cluster 1
5	2	3.5	cluster 1
8	5	0.5	cluster 2
9	6	0.5	cluster 2

Assignments are unchanged from iteration 1 ( $\{1,2,4,5\}$ vs $\{8,9\}$ ), so the algorithm has converged. Final objective:

J = (1-3)^2+(2-3)^2+(4-3)^2+(5-3)^2 + (8-8.5)^2+(9-8.5)^2 = 4+1+1+4+0.25+0.25 = 10.5.

Now the cautionary case — a bad initialization. Start instead with $\mu_1 = 1$ , $\mu_2 = 2$ (both centroids inside the left cluster). Assignment: points $1,2$ split trivially near their own centroid; $4$ is closer to $\mu_2=2$ (distance 2) than $\mu_1=1$ (distance 3); every remaining point $4,5,8,9$ is closer to $\mu_2$ . So cluster 1 $=\{1\}$ , cluster 2 $=\{2,4,5,8,9\}$ . Update: $\mu_1 = 1$ (unchanged, singleton), $\mu_2 = (2+4+5+8+9)/5 = 28/5 = 5.6$ . Next assignment: point 2 is distance 1 from $\mu_1{=}1$ and 3.6 from $\mu_2$ , stays cluster 1; point 4 is distance 3 from $\mu_1$ , 1.6 from $\mu_2$ , stays cluster 2; nothing changes further. Converged partition: $\{1,2\}$ vs $\{4,5,8,9\}$ , means $1.5$ and $6.5$ . Objective here:

J = (1-1.5)^2+(2-1.5)^2 + (4-6.5)^2+(5-6.5)^2+(8-6.5)^2+(9-6.5)^2 = 0.25+0.25+6.25+2.25+2.25+6.25 = 17.5.

This is a worse (higher) local optimum than the $J=10.5$ found from the first initialization, even though both runs fully converged, and both used the identical algorithm on the identical data. This trace is the concrete demonstration of the theoretical point above: Lloyd's algorithm guarantees monotonic decrease and convergence, but not which local minimum you land in.

Common pitfalls & traps

Pitfall	Why it happens	Consequence
Treating k-means output as globally optimal	The algorithm only guarantees local minimization of a non-convex objective	Reported clusters can vary run to run; results can be misleadingly "confident"
Forgetting to standardize features	Squared Euclidean distance is scale-sensitive; a feature measured in thousands dominates one measured in fractions	Clustering is driven almost entirely by the largest-scale feature
Ignoring the empty-cluster case	A centroid can end up with zero assigned points after an update	Division by zero / undefined mean if not explicitly handled
Choosing $K$ by eyeballing inertia alone	$J$ is monotonically non-increasing in $K$ by construction ( $K=n$ gives $J=0$ )	Naively minimizing $J$ over $K$ always favors more clusters; needs a principled method, see Choosing the Number of Clusters
Applying k-means to non-convex or unequal-density clusters	The assignment step always produces convex Voronoi regions around centroids	Rings, crescents, or clusters of very different density are split incorrectly
Running only once with random initialization	A single unlucky start can land in a much worse local optimum	Practitioners typically run many restarts and keep the lowest- $J$ result

Pioneers & history

The algorithm's name and its most common implementation trace to Stuart Lloyd, who developed it in 1957 at Bell Labs for pulse-code modulation (signal quantization), though the work was not published until 1982. Independently, E.W. Forgy described essentially the same iterative procedure in 1965 in the statistics literature, and James MacQueen coined the term "k-means" in a 1967 paper that also proposed an online variant updating centroids after each point. The convergence of the algorithm — that alternating exact minimization over two blocks of variables monotonically decreases a joint objective — is a special case of a much older and more general idea now called block-coordinate descent or, in its EM incarnation, the ideas behind Gaussian Mixture Models fitting.

Origins · 1 minStuart P. Lloyd and how this was discovered

Stuart P. Lloyd

1923–present · Bell Labs, quantization theory

Lloyd worked on the problem of quantizing an analog signal into discrete levels for transmission — effectively clustering a continuum of possible amplitudes into $K$ representative values that minimize expected squared quantization error. The iterative centroid-update rule he devised in an unpublished 1957 memorandum is mathematically identical to modern k-means; his original motivation was engineering signal compression, not statistical pattern discovery, and the technique migrated into general-purpose clustering only later.

History

1957

Stuart Lloyd devises the iterative centroid algorithm for scalar quantization at Bell Labs (published 1982).

1965

E.W. Forgy independently describes an equivalent iterative clustering procedure in statistics.

1967

James MacQueen names the method "k-means" and proposes an online, one-point-at-a-time variant.

1982

Lloyd's original 1957 memorandum is finally published, cementing the algorithm's name association.

2007

Arthur and Vassilvitskii introduce k-means++, a seeding strategy with provable approximation guarantees.

Where this matters

K-means is a default building block across data-driven engineering. In image and video compression, it performs vector quantization: pixel colors or motion vectors are clustered into a small codebook, and each pixel is replaced by its nearest codeword, directly reducing storage. In computer vision pipelines, k-means historically clustered local feature descriptors (like SIFT) into "visual words" for bag-of-visual-words image representations, and is still used to build codebooks for approximate nearest-neighbor search structures. In natural language and recommendation systems, k-means clusters embedding vectors — of documents, users, or products — to power segment-based personalization and exploratory analysis at scale. In bioinformatics, it groups gene-expression profiles or single-cell measurements into candidate cell types as an exploratory first pass. It is also a standard initialization routine for more expensive iterative methods, most notably providing a starting partition for Gaussian Mixture Models fit via expectation-maximization, since a good k-means solution is a fast, cheap proxy for a good initial guess at cluster means and covariances. Nearly every major machine learning library (scikit-learn, Spark MLlib, TensorFlow, cuML for GPU acceleration) ships a production-grade k-means implementation because its simplicity, interpretability, and near-linear scaling make it the first tool reached for whenever raw grouping of vector data is needed.

Ask the tutor: How does k-means++ initialization change the guarantees on how close the found solution is to the global optimum?

Does k-means always converge, and how many iterations does it take?

Yes, Lloyd's algorithm always converges in a finite number of iterations, because there are only finitely many ways to partition $n$ points into $K$ groups, and the objective $J$ strictly decreases whenever the partition changes between iterations (ties aside). So the sequence of partitions cannot cycle or decrease forever; it must reach a fixed point. In the worst case the number of iterations can be very large (theoretical worst-case bounds are exponential-ish on adversarial inputs), but in practice, real datasets typically converge within a few dozen iterations, and most implementations cap iterations with a `max_iters` safeguard regardless.

Why does k-means use squared Euclidean distance instead of plain Euclidean distance?

Squared distance is what makes the update step solvable in closed form: minimizing the sum of squared distances to a fixed point set yields the arithmetic mean, a simple, differentiable, convex problem. Minimizing the sum of plain (unsquared) Euclidean distances instead yields the geometric median, which has no closed-form solution and requires iterative numerical methods (e.g., Weiszfeld's algorithm) even to compute for one cluster. Squared distance also connects k-means directly to variance: WCSS is essentially $K$ times the average within-cluster variance, tying the method to familiar statistical quantities.

How is K chosen if k-means doesn't tell us automatically?

Raw inertia $J$ is not a valid criterion by itself because it is monotonically non-increasing in $K$ : adding more clusters can only help fit the data more tightly, all the way down to $J=0$ when $K=n$ . Practical methods instead look for a "knee" in $J$ versus $K$ (the elbow method), use silhouette scores that reward well-separated, well-formed clusters, or use information-theoretic penalties (like a Bayesian Information Criterion adaptation) that explicitly penalize extra clusters. This selection problem is developed fully in Choosing the Number of Clusters.

Is k-means sensitive to feature scaling, and what should be done about it?

Yes, strongly. Squared Euclidean distance treats every feature dimension as equally important in absolute units, so a feature ranging over thousands (like annual income) will dominate a feature ranging over single digits (like number of purchases), effectively drowning out the smaller-scale feature in every distance computation. Standard practice, connecting back to Distance Metrics and Similarity Measures, is to standardize each feature (subtract mean, divide by standard deviation) before running k-means, unless there is a specific domain reason the raw scales should matter.

How does k-means relate to Gaussian Mixture Models?

K-means can be seen as a limiting, simplified special case of fitting a Gaussian Mixture Model with expectation-maximization, where each Gaussian component is forced to have identical, spherical (isotropic) covariance, and instead of soft probabilistic assignments, each point is hard-assigned to its single nearest component. Gaussian Mixture Models relax both restrictions: they allow elliptical clusters via full covariance matrices, and they allow a point to have fractional membership across clusters. This makes GMMs strictly more expressive, at the cost of more parameters to estimate and a more expensive fitting procedure.

Why can two runs of k-means on the same data give different clusters?

Because Lloyd's algorithm is a local search over a non-convex objective, and its trajectory is entirely determined by where the centroids start. Different random initializations can lead the assignment/update alternation down different paths, converging to different local minima with different final cluster boundaries and different values of $J$ . This is not a bug; it is an intrinsic property of local optimization on a non-convex landscape, and it is precisely the motivation for smarter initialization schemes such as K-Means++ Initialization, which chooses spread-out starting centroids to make good local optima much more likely.

What happens if a cluster becomes empty during the update step?

If no point is assigned to a centroid after the assignment step, the "mean of its points" is undefined, and the naive update formula divides by zero. Real implementations detect this and re-seed the empty centroid — common strategies include moving it to the point currently farthest from its own assigned centroid (attacking the largest source of remaining error), or restarting it at a random data point. This is a necessary edge-case handling step that a correct implementation must include explicitly, since it is not automatically resolved by the core algorithm.