Lloyd's Algorithm and the Local Minima of k-Means

The k-means objective takes a set of n points and a fixed integer k and asks for the partition that minimizes the sum of squared distances from each point to the mean of its own group:

$$J = \sum_{i=1}^{n} \lVert x_i - c_{a(i)} \rVert^2$$

where a(i) is the cluster index assigned to point i and c_j is the mean of the points assigned to cluster j. The same quantity can be written as a sum over clusters of the pairwise squared distances divided by twice the cluster size, which shows the objective is a measure of within-group scatter rather than of distance to an external label. Minimizing J over all assignments is NP-hard for general k and dimension, so no algorithm computes the optimum in reasonable time. Lloyd’s algorithm instead descends to a fixed point of a two-step iteration, and each step is an exact minimization of J over one of the two unknowns.

Assignment and update are two exact minimizations of one objective

Holding the centroids fixed, the best assignment for point i is the index of the nearest centroid, because each point contributes only its own squared distance to the sum. Assigning every point to its nearest centroid partitions the space into the Voronoi cells of the centroid set, so the regions are convex polytopes, unbounded by default, and adjacent cells meet on hyperplanes. Ties are possible and any consistent tie-breaking rule is equally optimal for J; different choices give different subsequent iterations.

Holding the assignment fixed, the contribution of cluster j is sum over its members of ||x_i - c_j||^2, which is a strictly convex function of c_j with gradient 2 * sum (c_j - x_i). Setting the gradient to zero gives c_j = (1/n_j) * sum x_i, the arithmetic mean, and the Hessian is 2 n_j I, positive definite because the cluster is non-empty at that moment. So the update step is not a heuristic; it is the closed-form minimizer.

The consequence is a monotonicity guarantee. The assignment step cannot increase J because each point’s cost is individually minimized; the update step cannot increase it because each centroid is replaced by the argument that minimizes its own term. The sequence of objective values is non-increasing and bounded below by zero. Because there are finitely many partitions of n points into k groups, each strictly decreasing step lands on a partition that has not been visited before in the run, so the loop must terminate. If both steps produce no change, the state is a fixed point: no single point wants a different cluster, and no centroid wants to move.

Termination is not optimality

A fixed point is only a local minimum. Two runs on the same data with different starting centroids typically give different partitions with different objective values, and either may be far from the global optimum. The iteration count is bounded by the number of partitions, which is astronomically large, and the bound is not merely slack: there are constructions in the plane where the algorithm passes through a number of states exponential in the number of points before halting, so no polynomial iteration bound exists. Practically, the number of iterations is usually small, but that is an observation about typical data rather than a guarantee.

The remedy used in every implementation is restarts. The loop is run several times from different initial centroids and the run with the smallest final objective is kept. Since the objective of the best of several runs is a random variable, the number of restarts trades compute for the chance of reaching a better basin. A second, cheaper remedy is a bound: cluster assignment decisions that cannot change the objective can be skipped, and a run can be abandoned once its partial cost already exceeds the best complete cost found so far.

Seeding scheme	How centroids are chosen	Guarantee
Uniform sample	k distinct data points at random	None
k-means++	First point uniform, each later point with probability proportional to its squared distance to the nearest chosen centroid	Expected O(log k) approximation
Greedy	Point that most reduces the current objective	Constant factor, O(k) passes

Seeding changes which basin the loop reaches

Uniformly sampling data points as initial centroids is the simplest option and has no approximation guarantee at all; adversarial configurations make it arbitrarily bad. The k-means++ scheme samples the first centroid uniformly from the data, then samples each subsequent centroid from the remaining points with probability proportional to D(x)^2, where D(x) is the distance from x to the closest centroid chosen so far. Spreading the seed set this way competes with the cost of the optimal solution: the resulting clustering has objective value within a factor of 8(ln k + 2) of optimal with the standard analysis, and O(log k) in expectation.

The proof idea is worth knowing because it explains why the squared distance appears. Any cluster of the optimal solution that contains no chosen centroid must have some point, and the probability of selecting a point from that cluster on a given step is at least proportional to the fraction of the total potential that cluster contributes. Clusters that are expensive relative to the current potential are therefore likely to be seeded, and after O(log k) careful steps the expected excess potential has been amortized down.

The seeding pass itself costs O(n k d) arithmetic because each new centroid requires distances to all points, which is one extra iteration’s worth of work and is usually repaid by reaching a better fixed point. A cheaper variant samples a small number of candidate points per step and chooses the best by the potential reduction; another applies a few local search steps that swap centroids with data points while the objective improves.

Cost per iteration and the tricks that reduce it

One iteration computes n k squared distances in d dimensions, so the arithmetic is O(n k d). The identity ||x - c||^2 = ||x||^2 - 2 x . c + ||c||^2 lets the assignment step drop the ||x||^2 term, since it is constant for a fixed query point, and reduces the remaining work to a matrix product of shape n x d by d x k. That matters because the matrix product can run on cache-blocked and vectorized kernels, while the naive loop over pairs cannot.

Acceleration methods exploit the geometry rather than the arithmetic. If every point keeps a lower bound on its distance to every other centroid and centroids move by known amounts, the triangle inequality rules out most centroid comparisons without computing them; the fraction of comparisons skipped grows as clusters become separated. Mini-batch variants process a small random subset per iteration and update only the centroids that the batch touches, using per-cluster counts so that the running mean stays exact with respect to the points seen. That variant loses the monotone descent property of full-batch Lloyd because the assignment is computed on a sample, and its iterates are a stochastic process rather than a descent sequence.

Degenerate inputs deserve explicit handling. A centroid can lose all of its points, leaving its mean undefined; keeping the previous position is valid and keeps J unchanged, and re-seeding the centroid at the point farthest from its own centroid is the common repair. Duplicate points at identical coordinates produce ties in assignment, and a run that places two centroids at the same location is a fixed point for the assignment step but not a strict minimum. Points at identical coordinates in a cluster make the within-cluster term zero for that cluster.

The geometry the objective assumes

Squared Euclidean distance encodes a specific model of a cluster: an isotropic Gaussian, the same volume for every group, and convex cells. Elongated clusters, clusters of very different density, and concentric rings are all cut apart by the Voronoi structure, no matter how many restarts are used. The metric is a modelling decision rather than a detail, and alternatives change what the cells look like. A soft version of the same model, in which points carry membership probabilities and each group has its own covariance, is described in Gaussian Mixture Models and the EM Lower Bound, and its fixed points have the same local-minimum character. When group shapes are genuinely non-convex, the neighbourhood-based definition in DBSCAN, Density Reachability, and the Parameter Cliff avoids cells altogether.

Two comparisons are frequently attempted and are not valid. The first is comparing J across different values of k: J is non-increasing in k because extra centroids can only reduce the distance to the nearest one, so a larger k always looks better and an elbow in the curve is a heuristic, not a test. The second is comparing J between runs on different data or after different preprocessing, since scaling a feature changes distances directly. The metric question is treated in What a Distance Function Decides About Cluster Shape, and the model-selection question in Cluster Validity, Stability, and Choosing the Number of Clusters.

A compact implementation makes the two minimizations explicit:

import numpy as np

def lloyd(X, k, max_iter=100, seed=0):
    rng = np.random.default_rng(seed)
    C = X[rng.choice(len(X), k, replace=False)].copy()   # uniform seeding
    for _ in range(max_iter):
        # squared distances via ||x||^2 - 2 x.c + ||c||^2, argmin unaffected by ||x||^2
        d2 = (X ** 2).sum(1)[:, None] - 2.0 * X @ C.T + (C ** 2).sum(1)[None, :]
        a = d2.argmin(1)
        counts = np.bincount(a, minlength=k)
        # empty clusters keep their old position; the update is the exact mean otherwise
        new = np.array([X[a == j].mean(0) if counts[j] else C[j] for j in range(k)])
        if np.allclose(new, C):
            break
        C = new
    return C, a

The final partition depends on the seeding scheme, the tie-breaking order, and the number of restarts; the objective value of a single run says nothing about how close it is to the best partition that exists for that data.