🌀 k-means: Lloyd's Algorithm and k-means++ Initialization

By Better Wise

Last lesson, PCA gave you a lower-dimensional space via SVD. Clustering answers a different question: given points, which ones belong together? The discipline is the same one you applied to the linear models: no NumPy until the objective and both update rules are on paper. k-means is a two-step coordinate descent on one scalar, and once you can write that scalar you can debug every failure mode you will hit in the capstone.

What you will be able to do by the end

  • Write the k-means objective and derive both Lloyd steps as exact minimizations, not heuristics.
  • Implement k-means++ seeding in roughly fifteen lines of NumPy.
  • Justify a restart budget and a choice of k using inertia plus silhouette, not inertia alone.
  • Wrap the result in the same estimator shape you already use for the linear models, so it drops into the baseline comparison.

The inertia objective, derived on paper

$$ J = \sum_{n=1}^{N} \sum_{k=1}^{K} r_{nk} \lVert x_n - \mu_k \rVert^2 , \qquad r_{nk} \in {0,1}, \quad \sum_k r_{nk} = 1 $$

Lloyd's algorithm is block coordinate descent on $J$:

  1. Assignment. Fix the centroids $\mu$. The objective separates over points, so each point independently chooses its nearest centroid. Exact, closed form.
  2. Update. Fix the assignments $r$. Setting $\partial J / \partial \mu_k = -2 \sum_{n \in C_k} (x_n - \mu_k) = 0$ gives $\mu_k = \frac{1}{|C_k|}\sum_{n \in C_k} x_n$, the mean of the assigned points. Also exact.

Neither step can increase $J$, and $J \ge 0$, so the algorithm decreases the inertia _____ and stops at a fixed point. There are finitely many partitions, so it terminates. What it does not give you is a global minimum: the fixed point it lands on is a _____. That single sentence is why the second half of this lesson exists.

Important: k-means assumes clusters are roughly isotropic (spherical), similar in size, and comparable in scale. Standardize features before you run it, or the feature with the largest variance silently becomes the only feature.

Worked example: one Lloyd iteration and one k-means++ draw

Take six points: A(0,0), B(1,0), C(0,1), D(9,9), E(10,9), F(9,10), with $k=2$ and a deliberately bad seed from the first two rows: $\mu_1=(0,0)$, $\mu_2=(1,0)$.

Assignment. A goes to 1 (distance 0). B goes to 2 (0 versus 2). C goes to 1 (1 versus 2). D, E, F go to 2.

Update. $\mu_1 = (\tfrac13, \tfrac13)$, $\mu_2 = (\tfrac{28}{3}, \tfrac{28}{3}) \approx (9.33, 9.33)$.

Inertia. Cluster 1: $\tfrac29 + \tfrac59 + \tfrac59 = \tfrac43$. Cluster 2: same by symmetry. Total $J = \tfrac83 \approx 2.67$. Reassignment changes nothing, so this is a fixed point in one iteration.

Seeding. Now the 1-D version, $x = [0,1,2,10,11,12]$, $k=2$. The first center is drawn uniformly, and suppose it lands on 0. Squared distances to the nearest chosen center are $[0,1,4,100,121,144]$, summing to 370. The second center is drawn with probability proportional to those values: $p(1)=1/370 \approx 0.003$, $p(10)=100/370 \approx 0.27$, $p(12)=144/370 \approx 0.39$. The seed almost surely lands in the far group, which is exactly the failure the naive "take k random rows" seeding produces all the time.

For a fixed assignment the inertia-minimizing center is the _____; k-means++ draws each center after the first with probability proportional to _____.

Guided practice: build and judge the clusterer

Exercise 1 (paper, five minutes). Differentiate $J$ with respect to $\mu_k$ and recover the mean. Then write one sentence on why the assignment step cannot increase $J$.

Answer

$\partial J/\partial \mu_k = -2\sum_{n \in C_k}(x_n - \mu_k) = 0 \Rightarrow \mu_k = \frac{1}{|C_k|}\sum_{n \in C_k} x_n$. The Hessian is $2|C_k| I$, positive definite, so this is the minimizer. For the assignment step: given fixed centroids, moving any point to a different cluster replaces its term with a distance at least as large, and every point is optimized independently.

Exercise 2. Implement k-means++ seeding for the project dataset, returning (k, d) initial centers.

Answer

import numpy as np

def kmeans_pp(X, k, rng):
    n = X.shape[0]
    centers = np.empty((k, X.shape[1]))
    centers[0] = X[rng.integers(n)]
    d2 = ((X - centers[0]) ** 2).sum(axis=1)
    for j in range(1, k):
        probs = d2 / d2.sum()
        idx = rng.choice(n, p=probs)
        centers[j] = X[idx]
        d2 = np.minimum(d2, ((X - centers[j]) ** 2).sum(axis=1))
    return centers

Keep d2 as a running minimum. Recomputing all distances to all chosen centers each round is correct but quadratic in $k$, and you will notice on the capstone dataset.

Exercise 3. Your Lloyd loop converges in under 30 iterations on this dataset with k=5, but the inertia varies by 40 percent across seeds. Decide what to report.

Guidance

Do not report the best run as if it were typical. Run 10 restarts with k-means++ seeding, keep the lowest-inertia model, and report the median and spread of final inertia alongside it. Then check the silhouette of the chosen model. If the best inertia comes with a low silhouette, the extra tightness is probably a pair of clusters carved out of one true group.

Exercise 4. Compute the silhouette for the chosen model.

Answer

For point $i$ in cluster $C$: $a(i)$ is the mean distance to the other members of $C$; $b(i)$ is the minimum over other clusters of the mean distance to that cluster. Then $s(i) = (b-a)/\max(a,b)$, and the score is the mean over points. Silhouette ranges from _____ to 1; near 1 means well inside its cluster, near 0 means on a boundary, negative means likely misassigned. Cost is $O(n^2)$ distances, so subsample to a few thousand points on large data.

Where this breaks: mistakes and fixes

Symptom Cause Fix
Two centroids on one true group, one group split poor seeding k-means++ plus 10 restarts, tracked by inertia distribution
RuntimeWarning: Mean of empty slice a cluster lost all points re-seed the empty centroid to the point farthest from its centroid
Iteration count explodes comparing full centroid arrays with == stop on centroid shift norm or assignment change, tolerance 1e-6 to 1e-4
Best k is always k_max selecting by inertia inertia decreases with k and hits 0 at k = n; use silhouette or an elbow plus domain knowledge
Clusters dominated by one feature unscaled inputs standardize; the distance is Euclidean and unit-blind

Seeding and restart rules of thumb

  • Seed with k-means++, then spend your compute on restarts, not iterations. Ten restarts beats one run with a hundred extra iterations.
  • Track quality with both numbers: inertia for "did this run land in a better basin", silhouette for "is this partition believable".
  • Set n_init=10 in your library baseline too. A baseline crippled by n_init=1 is not a fair comparison, and the capstone protocol will call that out.
  • When silhouette is flat across k and inertia shows no elbow, your data probably is not spherical. That is the signal to reach for the next lesson, not to keep tuning k.
  • Fix the RNG seed per experiment and reuse it across your implementation and the library baseline so restart budgets are comparable.

Closing notes and what EM adds next

You now have the full loop: derive $J$, alternate assignment and mean update, seed with k-means++ and justify restarts, judge with inertia plus silhouette. The one assumption to carry forward is the hard assignment: every point belongs to exactly one centroid. Next lesson replaces that with soft responsibilities and a full covariance per component, which is the Gaussian mixture model fit by expectation maximization. Same alternating structure, one extra derivation.

Front: What does k-means minimize? Back: Inertia: the sum of squared distances from each point to its assigned centroid.

Front: What is the assignment step of Lloyd's algorithm? Back: Assign every point to its nearest centroid, which is the exact minimizer for fixed centroids.

Front: What is the update step of Lloyd's algorithm? Back: Set each centroid to the mean of its assigned points, the exact minimizer for a fixed assignment.

Front: Why does Lloyd's algorithm always terminate? Back: Inertia never increases and there are finitely many partitions, so it reaches a fixed point.

Front: What kind of optimum does k-means return? Back: A local minimum whose location depends on the initialization.

Front: How does k-means++ choose the first center? Back: Uniformly at random from the data points.

Front: How does k-means++ choose later centers? Back: With probability proportional to the squared distance to the nearest already chosen center.

Front: Why does k-means++ reduce bad runs? Back: It spreads seeds apart, so two seeds rarely land in the same true cluster.

Front: What range does the silhouette score take? Back: From -1 to 1; near 1 is well inside a cluster, near 0 is on a boundary, negative suggests misassignment.

Front: What is the cost of an exact silhouette score? Back: O(n^2) pairwise distances, so subsample for large datasets.

Front: Why can't you select k by minimum inertia? Back: Inertia decreases monotonically with k and reaches 0 when k equals the number of points.

Front: How do you handle an empty cluster? Back: Re-seed that centroid, typically at the point farthest from its own centroid, then continue.

Front: Why standardize features before k-means? Back: The objective is Euclidean, so a high-variance feature dominates the distance.

Front: What should the convergence check compare? Back: Centroid movement or assignment changes, not exact array equality.