← All notes · ML revision (DA5400)

L17: Initialization of Centroids in K-means (Random Restarts & K-means++)

Source: https://www.youtube.com/watch?v=Cbrdtxq6bdk (≈19 min) Previous: L13 · L14 · L15 · L16

The lecture explains K-means++ and states its guarantee without proof ("constant times the optimum"). Section 5 adds: the exact theorem, a full d-dimensional proof of the key lemma (random seeding costs 2× the optimum in expectation), why the distance is squared, why the deterministic version breaks, efficient implementation, scalable variants, and a tested comparison.

Timestamped outline

1. Why initialization matters

• The convergence proof (L15) works for any starting point. But where the algorithm ends up, which local minimum, depends heavily on where it starts.
• Goal: start in a way that helps avoid bad local minima, like the rectangle example in L15 §7.5.

2. Simple strategies

(a) Random partition. Throw each point into one of the K boxes uniformly at random.

• It works, but it ignores the structure of the data entirely.
• Side effect: every initial mean comes out close to the global mean, so the starting means are bunched together.

(b) Random data points as means (the "Forgy" method). Pick K data points uniformly at random and treat them as the initial means.

• This initializes the means rather than the assignments. The first reassignment step then sends every point to its nearest chosen point.
• That isn't convergence yet: the means get recomputed and Lloyd's continues as usual.

(c) Random restarts. Repeat (b) several times with fresh random draws, run Lloyd's each time, and keep the partition with the lowest objective. This is very common in practice.

3. K-means++

Intuition: at the end, the K means act as K "representatives" of the data. If two of them were close together, the two clusters would be redundant. So start with means that are spread apart, but choose them with randomness.

Algorithm:

1. Choose μ₁⁰ uniformly at random from {x₁, …, xₙ}.

2. For l = 2, …, K:

   • Give every point x a score:

     $$ \text{score}(x)=\min_{j=1,\dots,l-1}\lVert x-\mu_j^0\rVert^2 $$

     This is the squared distance to the nearest mean chosen so far. If even the nearest mean is far away, then every chosen mean is far away.

   • Sample μₗ⁰ with probability ∝ score:

     $$ P(\mu_l^0=x_i)=\frac{\text{score}(x_i)}{\sum_{j}\text{score}(x_j)} . $$

3. Run ordinary Lloyd's algorithm from these K means.

Sampling example: three candidates with scores 10, 20 and 30. The total is 60, so the probabilities are 1/6, 1/3 and 1/2. The farthest point is the most likely pick, but not a certain one.

Note: points that have already been chosen have score 0, so they can never be picked again.

4. Why randomise? And the guarantee

• Deterministic alternative: always pick the point with the largest score. For any fixed deterministic rule, you can construct an adversarial dataset that defeats it.
• Randomising keeps the "go far away" preference, because the max is the most likely pick, while still allowing other choices. That makes an average-case (expected) guarantee possible.
• The guarantee: the randomness comes only from the seeding; Lloyd's itself is deterministic. Over that randomness,

$$ \mathbb{E}\Big[\sum_i\lVert x_i-\mu_{z_i}\rVert^2\Big]\;\le\;c\cdot\min_{z_1..z_n}\sum_i\lVert x_i-\mu_{z_i}\rVert^2 . $$

The lecture's illustration: if the best possible objective is 10, then on average K-means++ might give something like ≤ 5 × 10. It won't be optimal, but it also won't be terrible. (The exact factor is in §5.1.)

Trade-offs:

• K-means++ needs K sequential passes over the data, each computing scores. That's expensive for something like a billion points.
• Uniform random seeding with restarts is cheaper and still widely used, especially for high-dimensional data.
• K-means++ is a solid default. The right choice depends on your compute budget.

5. Going further: d-dimensional analysis

This section goes beyond the video.

5.1 The exact theorem

Arthur & Vassilvitskii (2007). Let OPT be the optimal K-means cost. Then for the K-means++ seeds alone, before any Lloyd iterations:

$$ \mathbb{E}[F_{\text{seed}}]\;\le\;8(\ln K+2)\cdot\text{OPT}. $$

• Lloyd's iterations can only lower the cost (L15), so the final result satisfies the same bound.
• The factor is O(log K). It depends on K, not on n or d. For K = 10 that's ≈ 34, a loose worst-case bound. In practice it does much better.
• The paper also shows the bound is tight up to constants: some datasets really do reach Ω(log K).

5.2 Key lemma (proved in ℝᵈ): one random point costs 2× the optimum

Take a single optimal cluster C with mean μ, and suppose its centre is a uniformly random point a ∈ C. What's the expected cost?

Use the L15 identity, Σ_{x∈C}‖x − v‖² = Σ_{x∈C}‖x − μ‖² + |C|·‖μ − v‖², with v = a:

$$ \mathbb{E}_a\Big[\sum_{x\in C}\lVert x-a\rVert^2\Big] =\underbrace{\sum_{x\in C}\lVert x-\mu\rVert^2}_{\text{OPT}(C)} +|C|\cdot\underbrace{\frac{1}{|C|}\sum_{a\in C}\lVert a-\mu\rVert^2}_{=\,\text{OPT}(C)/|C|} =2\,\text{OPT}(C). $$

So choosing a centre uniformly from the right cluster costs, on average, exactly twice the optimum, in any dimension. That's why step 1 of K-means++ (uniform) is fine.

The hard part of the proof (stated, not proved here): when a later centre is drawn with D² weighting and it lands in a not-yet-covered optimal cluster, the expected cost of that cluster is at most 8·OPT(C). The D² weighting makes "landing in an uncovered cluster" likely, because uncovered clusters contain most of the remaining score. A harmonic-series argument over the K rounds gives the (ln K + 2) factor.

5.3 Why squared distance? (D^α sampling)

You could sample with probability ∝ D(x)^α for any α ≥ 0.

• α = 0 gives uniform sampling, the Forgy method. This has no constant-factor guarantee: K well-separated clusters, one of them tiny, and uniform sampling will usually miss it.
• α = 2 matches the objective, which is also a squared distance. A point's score then is its contribution to the cost if no further centre were added. So you sample in proportion to "how much cost this point is currently causing". That alignment is what makes the proof work.
• α → ∞ concentrates all the probability on the farthest point: deterministic farthest-first (§5.4).

General result: D^α seeding is an O(2^{2α} log K) approximation for the cost Σ‖x − μ‖^α (K-medians uses α = 1).

5.4 Why the deterministic version (farthest-first) fails

Farthest-first traversal (Gonzalez, 1985) always picks argmax score. It is a provable 2-approximation for a different problem, K-center (minimise the maximum radius). But for K-means it is fragile to outliers:

• One outlier at distance 1000 has score 10⁶, while every "normal" point has score around 1.
• Farthest-first picks the outlier with certainty and wastes a centre on a single point.
• K-means++ gives the outlier probability ≈ 10⁶ / (10⁶ + n·1). That's still high, but if there are many genuine far-apart clusters, their combined score competes. And over restarts it won't make the same mistake every time.

The tested demo in §5.7 shows exactly this.

5.5 Efficient implementation: O(nKd) total, the same as one Lloyd iteration

The naive approach recomputes the distance from every point to every chosen centre in every round, giving O(nK²d) work. Instead, keep a running minimum:

$$ D^2_{(l)}(x)=\min\big(D^2_{(l-1)}(x),\;\lVert x-\mu_l\rVert^2\big)\quad\Rightarrow\quad O(nd)\text{ per new centre}. $$

Sampling ∝ D²: compute the cumulative sum of the scores, draw u ~ Uniform(0, ΣD²), and binary-search for it. That's O(n) for the cumulative sum and O(log n) per draw. So seeding costs about one extra Lloyd iteration, which is cheap compared with the iterations that follow. The lecture's "slow" concern applies mainly at billion-point scale, where K sequential passes are the bottleneck.

Greedy K-means++ (the scikit-learn default): in each round, draw 2 + ⌊ln K⌋ candidates by D² sampling and keep the one that lowers the current cost the most. This is empirically better and costs a little more.

5.6 Scaling up: K-means|| (Bahmani et al., 2012)

This fixes the "K sequential passes" problem:

1. Start with one uniform centre.
2. For about 5 rounds (O(log ψ) in theory): sample each point independently with probability ∝ ℓ·D²(x)/ΣD², where ℓ ≈ 2K. Each round adds about ℓ centres and is a single parallel pass, which suits map-reduce.
3. You now have roughly ℓ × rounds candidate centres. Weight each by the number of points closest to it, then run ordinary K-means++ on this small weighted set to get K seeds.

The guarantee is still O(log K). This is the default seeding in Spark MLlib.

5.7 Tested comparison (n = 1000, K = 10, d = 5, 50 seeds each)

import numpy as np

def init_random(X, K, rng):
    return X[rng.choice(len(X), K, replace=False)]

def init_kmeanspp(X, K, rng):
    n = len(X)
    M = [X[rng.integers(n)]]                       # first mean: uniform
    D2 = ((X - M[0])**2).sum(1)                    # score = min squared distance
    for _ in range(1, K):
        p = D2 / D2.sum()                          # normalise scores
        M.append(X[rng.choice(n, p=p)])            # sample ∝ score
        D2 = np.minimum(D2, ((X - M[-1])**2).sum(1))   # O(nd) incremental update
    return np.array(M)

def init_farthest(X, K, rng):
    n = len(X)
    M = [X[rng.integers(n)]]
    D2 = ((X - M[0])**2).sum(1)
    for _ in range(1, K):
        M.append(X[D2.argmax()])                   # deterministic: take the max
        D2 = np.minimum(D2, ((X - M[-1])**2).sum(1))
    return np.array(M)

def lloyd(X, M, iters=300):
    for _ in range(iters):
        z = ((X[:, None] - M[None])**2).sum(-1).argmin(1)
        M_new = np.stack([X[z == k].mean(0) if np.any(z == k) else M[k] for k in range(len(M))])
        if np.allclose(M_new, M): break
        M = M_new
    return ((X - M[z])**2).sum()

rng = np.random.default_rng(0)
K, d = 10, 5
centres = rng.uniform(-20, 20, (K, d))
X = np.vstack([rng.normal(c, 1.0, (100, d)) for c in centres])
X_out = np.vstack([X, rng.uniform(-200, 200, (5, d))])      # same data + 5 far outliers

for name, data in [("clean", X), ("with outliers", X_out)]:
    print(f"--- {name} (n={len(data)}, K={K}, d={d}) ---")
    for init in (init_random, init_kmeanspp, init_farthest):
        F = np.array([lloyd(data, init(data, K, np.random.default_rng(s))) for s in range(50)])
        print(f"{init.__name__:14s} mean F = {F.mean():9.0f}   best = {F.min():9.0f}   worst = {F.max():9.0f}")

Output when I ran it (final objective after Lloyd's, over 50 runs):

How to read it:

• On clean, well-separated blobs, farthest-first is perfect every time. Each blob's farthest point is in a new blob.
• Add just 5 outliers and farthest-first gets worse than K-means++ on average. It spends 5 of its 10 centres on the outliers every time. Even its best run is worse than K-means++'s mean.
• Random seeding has the widest spread. Its best run matches the others, which is why many restarts can rescue it.
• K-means++ is the most robust of the three: never the worst, and the best on average with outliers present.

Key takeaways

1. Lloyd's converges from any start, but the start decides which local minimum it reaches.
2. Forgy + restarts: pick K random data points, run, repeat, and keep the lowest F. It's cheap and common.
3. K-means++: choose the first centre uniformly, then each later centre ∝ D²(x), the squared distance to the nearest centre so far. It spreads the centres while staying randomised.
4. Guarantee: E[cost] ≤ 8(ln K + 2)·OPT. The key lemma is that a uniform random point from a cluster costs exactly 2·OPT in expectation.
5. Deterministic farthest-first is great for K-center but breaks with outliers. Randomisation gives robustness and makes an expected guarantee possible.
6. With a running minimum, seeding costs O(nKd), about one Lloyd iteration. K-means|| makes it parallel.

Self-check questions

1. Scores are 4, 4, 2 and 0 for four points. What is the sampling probability of each? Why is the 0-score point special? (0.4, 0.4, 0.2, 0. Score 0 means it's already a centre.)
2. Prove that E_a[Σ_{x∈C}‖x − a‖²] = 2·OPT(C) for a uniform over C, using the L15 identity.
3. Why does sampling ∝ D (not D²) fit K-medians better than K-means?
4. Construct a 1-D dataset with K = 2 where farthest-first gives an arbitrarily bad K-means cost. (Hint: one far outlier plus two tight, well-separated clusters.)
5. How does the running-minimum trick bring seeding from O(nK²d) down to O(nKd)?
6. For farthest-first specifically, explain how an adversary who knows the rule can build a dataset that defeats it. Then explain why the same trick doesn't work against K-means++. (Note: deterministic constant-factor algorithms for K-means do exist, e.g. local search, but they are much slower. The lecture's point is about simple greedy seeding rules.)