k-Means

The previous lesson froze the procedure deliberately: the starting points were chosen with a single rule, two refinement rounds were run, the number of clusters was held at 4. The only thing that changed was the similarity definition. This lesson fixes the definition — standard score with squared distance — and frees two of the three frozen constants: the iteration continues until it stops, and the starting points are chosen at random.

The resulting procedure is short: every point is assigned to the nearest centroid, every centroid moves to the mean of its own cluster, and this repeats until nothing changes. The lesson measures that this produces a decision in two separate places. The first is the start: how many distinct partitions come out when the same number of clusters is run from different starting points. The second is the number of clusters itself; k is a hyperparameter given from outside, and choosing it is a search.

• CL9. The dataset, columns, split and measure are as in the previous lesson: the fictional subscriber table, seed 20260218, the 756 subscribers of the training set, six numeric columns, standard score, squared distance. The label is again not generated.
• CL10. The measure is again the within-cluster sum, per point, smaller is better.
• CL11. The baseline is random assignment and is computed separately for every k; because the baseline changes with the number of clusters, a single number cannot be written for the whole table.
• CL12. The starting points are chosen from distinct points of the dataset; if a cluster ends up empty, its centroid is left where it was.
• CL13. The iteration stops when the assignment does not change for two rounds in a row; the upper limit is 60 rounds and no run reaches it. A floating-point number is never compared with == — the stopping criterion is equality of the assignment lists, that is, an integer comparison.
• CL14. The difference between two partitions is measured with the differing point pair ratio built in the previous lesson.
• CL15. In the sweep over the number of clusters, 12 starts are tried for every k and the best is taken; the number of candidates tried is written. The reported number is also a function of the search budget.

Where the Iteration Stops

The procedure is the alternating repetition of two steps, and each step lowers the measure from a different direction. The assignment step holds the centroids fixed and gives every point to its nearest one. The update step holds the assignment fixed and moves every centroid to the mean of its own cluster; because the point that minimizes the sum of squared distances is the mean, this step too lowers the measure or leaves it the same. Since the measure never rises and is bounded below, the iteration must stop. About where it stops, this argument says nothing.

# k_means.py — MODEL. The FICTIONAL set and the standard score definition of lesson 01
# are rebuilt with the same seed; the label is still not generated.
import math
import statistics

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [(0.28, 21), (0.22, 17), (0.18, 26), (0.14, 14), (0.18, 23)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value():
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u, w):
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


DATA = []
for i in range(RAW):
    r = rng(SEED + i)
    base = REGION[pick(r(), [x[0] for x in REGION])][1]
    household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    rating = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    if r() < 0.046:
        continue
    r, v = rng(SEED + 17001 + i), []
    for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.append(0.0 if r() < 0.038 else math.floor(
            base * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - d * 0.05) * 100 + 0.5) / 100)
        r()
    if v:
        x = {"mean_consumption": round(sum(v) / len(v), 2), "volatility": round(max(v) - min(v), 2),
             "household": household, "satisfaction": rating, "period": len(v)}
        x["per_person"] = round(x["mean_consumption"] / household, 3)
        DATA.append(x)

FIELDS = ["mean_consumption", "volatility", "per_person", "household", "satisfaction", "period"]


def shuffle(data, seed):
    r, s = rng(seed), list(range(len(data)))
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return [data[i] for i in s]


TRAIN = shuffle(DATA, SEED + 90000)[:756]
RAW_VALUES = [[x[a] for a in FIELDS] for x in TRAIN]
MEANS = [statistics.fmean(c) for c in zip(*RAW_VALUES)]
STD_DEVS = [statistics.pstdev(c) for c in zip(*RAW_VALUES)]
V = [[(v[j] - MEANS[j]) / STD_DEVS[j] for j in range(len(FIELDS))] for v in RAW_VALUES]
N = len(V)

squared = lambda a, b: sum((p - q) * (p - q) for p, q in zip(a, b))
centroid = lambda g: [statistics.fmean(c) for c in zip(*g)]


def within_sum(assign, k):                  # within-cluster sum, per point
    t = 0.0
    for c in range(k):
        g = [V[i] for i in range(N) if assign[i] == c]
        if g:
            t += sum(squared(v, centroid(g)) for v in g)
    return t / N


def starting_points(r, k):                  # pick k distinct points as initial centroids
    s = []
    while len(s) < k:
        j = int(r() * N)
        if j not in s:
            s.append(j)
    return [V[j] for j in s]


def k_means(M, limit=60):                   # MODEL: assign, update centroids, repeat
    k, previous = len(M), None
    for it in range(1, limit + 1):
        a = [min(range(k), key=lambda c: squared(v, M[c])) for v in V]
        if a == previous:
            return a, it
        previous = a
        M = [centroid([V[i] for i in range(N) if a[i] == c]) or M[c] for c in range(k)]
    return a, limit


r = rng(SEED + 5100)
M0 = starting_points(r, 4)
print("convergence trace from a single start (k = 4)")
for it in range(1, 9):
    a = [min(range(4), key=lambda c: squared(v, M0[c])) for v in V]
    print(f"  round {it}  within-cluster sum {within_sum(a, 4):.4f}"
          f"  cluster sizes {sorted(a.count(c) for c in range(4))}")
    M0 = [centroid([V[i] for i in range(N) if a[i] == c]) or M0[c] for c in range(4)]

convergence trace from a single start (k = 4)
  round 1  within-cluster sum 4.2484  cluster sizes [109, 159, 233, 255]
  round 2  within-cluster sum 3.9366  cluster sizes [148, 166, 213, 229]
  round 3  within-cluster sum 3.7729  cluster sizes [153, 180, 198, 225]
  round 4  within-cluster sum 3.6519  cluster sizes [145, 175, 201, 235]
  round 5  within-cluster sum 3.6125  cluster sizes [140, 170, 215, 231]
  round 6  within-cluster sum 3.6058  cluster sizes [137, 172, 221, 226]
  round 7  within-cluster sum 3.6050  cluster sizes [138, 171, 222, 225]
  round 8  within-cluster sum 3.6049  cluster sizes [138, 170, 223, 225]

The number that was 4.2484 in the first round falls to 3.6049 by the eighth, and the fall slows quickly: from round one to round three 0.4755 is gained, from round five to round eight only 0.0076. The cluster sizes settle too, the 109–255 range of the first round narrowing to 138–225. The two-round procedure of the previous lesson had been cut off at 3.8910; here even the third round goes below it.

What this trace does not say matters more. The number falling and stopping does not show that where it stopped is the lowest number attainable. Every round tries partitions neighboring the previous one; no round tries a placement far away. The procedure therefore stops at a solution near where it started.

The Same k, a Different Start

How decisive the start is gets measured by running thirty different sets of starting points at the same number of clusters. The identity of a partition, independent of cluster numbers, is which points stand together.

def identity(a, k):                         # the label-independent identity of a partition
    return frozenset(frozenset(i for i in range(N) if a[i] == c) for c in range(k))


def disagreement(a, b):                     # on how many point pairs the two partitions differ
    differing = sum((a[i] == a[j]) != (b[i] == b[j])
                    for i in range(N) for j in range(i + 1, N))
    return 2 * differing / (N * (N - 1))


r, RUNS = rng(SEED + 7700), 30
RESULTS = []
for s in range(RUNS):
    a, it = k_means(starting_points(r, 4))
    RESULTS.append((within_sum(a, 4), it, a))
RESULTS.sort(key=lambda x: x[0])
DISTINCT = {identity(a, 4) for _, _, a in RESULTS}
print(f"{RUNS} random starts, k = 4 -> distinct partitions {len(DISTINCT)}")
print("  sorted within-cluster sums")
for b in range(0, RUNS, 10):
    print("   " + " ".join(f"{u:.4f}" for u, _, _ in RESULTS[b:b + 10]))
print(f"  best {RESULTS[0][0]:.4f}, median {RESULTS[RUNS // 2][0]:.4f}, "
      f"worst {RESULTS[-1][0]:.4f}, range {RESULTS[-1][0] - RESULTS[0][0]:.4f}")
print(f"  convergence rounds at least {min(t for _, t, _ in RESULTS)}, "
      f"at most {max(t for _, t, _ in RESULTS)}")
print(f"  lesson 01's single-pass procedure gave 3.8910; improvement "
      f"{3.8910 - RESULTS[0][0]:.4f}")
print(f"  best - second best: number difference {RESULTS[1][0] - RESULTS[0][0]:.6f}, "
      f"differing point pairs {disagreement(RESULTS[0][2], RESULTS[1][2]):.4f}")
print(f"  best - worst:       number difference {RESULTS[-1][0] - RESULTS[0][0]:.4f}, "
      f"differing point pairs {disagreement(RESULTS[0][2], RESULTS[-1][2]):.4f}")

30 random starts, k = 4 -> distinct partitions 26
  sorted within-cluster sums
   3.5649 3.5649 3.5656 3.5656 3.5656 3.5658 3.5711 3.5711 3.5711 3.5711
   3.5711 3.5712 3.5712 3.5727 3.5903 3.5906 3.5936 3.5936 3.5941 3.5955
   3.5956 3.5957 3.5981 3.5995 3.6026 3.6077 3.6165 3.7467 3.7744 3.7763
  best 3.5649, median 3.5906, worst 3.7763, range 0.2115
  convergence rounds at least 11, at most 52
  lesson 01's single-pass procedure gave 3.8910; improvement 0.3261
  best - second best: number difference 0.000028, differing point pairs 0.0037
  best - worst:       number difference 0.2115, differing point pairs 0.2844

Thirty runs gave 26 distinct partitions. The same data, the same procedure, the same number of clusters; only the starting points changed. The numbers sit in a narrow band — between 3.5649 and 3.7763 — and twenty-seven of the thirty runs are below 3.62. Since the random assignment baseline is 5.9730 at this number of clusters, every run beats the baseline by a wide margin; even the worst is 2.1967 below it.

The bottom two lines are the real finding of the lesson. The difference between the numbers of the best two runs is 0.000028, the same number to four decimals; yet the two partitions disagree on 0.0037 of the point pairs. Between the best and the worst run the difference in numbers is 0.2115 and the disagreement is 0.2844 — more than a quarter of the pairs. The within-cluster sum reduces a partition to a single number, and what is preserved when that number stays the same is not who stands in the same cluster as whom.

The convergence rounds vary between 11 and 52, and no run reaches the 60-round limit. The iteration really does stop, but where it stops depends on where the starting points fall.

The Number of Clusters Is a Hyperparameter

k is not something the model learns from the data but a number given from outside — a hyperparameter in the terminology of the Supervised Learning course. Choosing it by looking at the measure does not work directly: as the number of clusters grows, the within-cluster sum falls in only one direction, and the smallest value comes from the partition in which every point is its own cluster. The reading is therefore made not from the number itself but from the successive differences; the place where the differences shrink sharply is called the elbow.

def random_baseline(k, seed):
    r = rng(seed)
    return within_sum([int(r() * k) for _ in range(N)], k)


STARTS, CANDIDATES, WITHIN = 12, 0, {}
print(f"{'k':>3}{'baseline':>9}{'method':>9}{'difference':>11}{'successive difference':>23}{'drop ratio':>13}")
for k in range(1, 11):
    r = rng(SEED + 8800 + k)
    WITHIN[k] = min(within_sum(k_means(starting_points(r, k))[0], k) for _ in range(STARTS))
    CANDIDATES += STARTS
    baseline = statistics.fmean(random_baseline(k, SEED + 9900 + 50 * k + t) for t in range(5))
    d = "" if k == 1 else f"{WITHIN[k - 1] - WITHIN[k]:23.4f}"
    o = "" if k == 1 else f"{(WITHIN[k - 1] - WITHIN[k]) / WITHIN[k - 1]:13.3f}"
    print(f"{k:>3}{baseline:>9.4f}{WITHIN[k]:>9.4f}{baseline - WITHIN[k]:>11.4f}{d}{o}")
print(f"candidates tried {CANDIDATES} runs ({STARTS} starts x 10 cluster counts)")

STEP = {k: WITHIN[k - 1] - WITHIN[k] for k in range(2, 11)}
biggest_drop = max(range(2, 10), key=lambda k: STEP[k] - STEP[k + 1])
threshold_k = min(k for k in range(2, 11) if STEP[k] / WITHIN[k - 1] < 0.10)
slope = (WITHIN[10] - WITHIN[1]) / 9
deviation = {k: (WITHIN[1] + slope * (k - 1)) - WITHIN[k] for k in range(2, 10)}
chord_k = max(deviation, key=deviation.get)
runner_up = sorted(deviation, key=deviation.get)[-2]
print("\nthree separate readings of the elbow")
print(f"  k where the successive difference falls the most  {biggest_drop}")
print(f"  k where the drop ratio first falls below 0.10     {threshold_k}")
print(f"  k farthest from the chord joining the endpoints   {chord_k}"
      f"   (deviation {deviation[chord_k]:.4f}; runner-up k = {runner_up}, deviation {deviation[runner_up]:.4f})")
print(f"  successive differences are not one-directional: k = 9 gives {STEP[9]:.4f}, k = 10 gives {STEP[10]:.4f}")

  k baseline   method difference  successive difference   drop ratio
  1   6.0000   6.0000     0.0000
  2   5.9928   4.8521     1.1407                 1.1479        0.191
  3   5.9894   3.9996     1.9897                 0.8524        0.176
  4   5.9730   3.5655     2.4075                 0.4342        0.109
  5   5.9721   3.2151     2.7570                 0.3504        0.098
  6   5.9635   2.8969     3.0666                 0.3182        0.099
  7   5.9486   2.6732     3.2754                 0.2237        0.077
  8   5.9384   2.5020     3.4365                 0.1712        0.064
  9   5.9443   2.3874     3.5568                 0.1146        0.046
 10   5.9181   2.2570     3.6611                 0.1304        0.055
candidates tried 120 runs (12 starts x 10 cluster counts)

three separate readings of the elbow
  k where the successive difference falls the most  3
  k where the drop ratio first falls below 0.10     5
  k farthest from the chord joining the endpoints   4   (deviation 1.1869; runner-up k = 3, deviation 1.1686)
  successive differences are not one-directional: k = 9 gives 0.1146, k = 10 gives 0.1304

The baseline column alone is a warning. The number for random assignment falls too as k grows: from 6.0000 to 5.9181. The fall is small but real, because every new cluster estimates its centroid from its own points and the estimation noise pulls the measure down. If even the number of a procedure that does not learn improves as k grows, the whole of the method’s fall cannot be attributed to structure. That is why the difference column is needed.

The elbow, on the other hand, does not give a single answer from the table. Three common readings were tried on the same ten rows and the three chose three separate numbers of clusters: the place where the successive difference falls the most is k = 3, the place where the drop ratio first falls below 0.10 is k = 5, and the point farthest from the chord joining the endpoints is k = 4. In the third, the winner’s deviation is 1.1869 and the runner-up’s 1.1686 — the gap is under 2%. The last line shuts the table altogether: the successive differences are not even one-directional, rising from 0.1146 at k = 9 to 0.1304 at k = 10.

The number of candidates tried is 120 runs, and this number stands beside every reported cell. The 3.5655 in the k = 4 row is the best of twelve starts; the search with thirty starts had found 3.5649 for the same k. The 0.0006 between them comes not from the data but from the search budget.

Summary

• The procedure is the alternating repetition of two steps and each step lowers the measure or leaves it the same; the iteration must therefore stop, but there is no guarantee that where it stops is the best.
• Thirty random starts gave 26 distinct partitions at k = 4; the numbers varied between 3.5649 and 3.7763, and convergence between 11 and 52 rounds.
• A difference in numbers and a difference in partitions are not the same thing: two runs that give the same number to four decimals disagree on 0.0037 of the point pairs, the best and the worst run on 0.2844.
• The random assignment baseline also falls with the number of clusters, from 6.0000 to 5.9181; once this baseline is subtracted, the method’s difference rises from 2.4075 to 3.6611.
• Three separate readings of the elbow chose three separate numbers of clusters — 3, 5 and 4 — and the successive differences are not even one-directional; the chosen number is the result of a search budget of 120 runs.

Next Step

The two constraints of this procedure are tied together: the number of clusters must be given before starting, and every run produces a single partition. Three separate searches were made for k = 3, k = 4 and k = 5 and the three results came out unrelated — the clusters of k = 3 are not unions of the clusters of k = 5. The next lesson takes a family that reverses the order: every point starts as its own cluster, the nearest two clusters are merged, and a nested family of clusters is produced up to a single cluster. The number of clusters is here not an input but the answer to the question of which height the tree is cut at. In return a new decision appears: “the nearest two clusters” cannot be derived from point distance, and the distance between clusters must be defined separately. The lesson counts on how many points three linkage criteria disagree on the same data and prints the dendrogram as text.