III-26 · k-Means

Give every assignment a distance and every center a mean

Drag a centroid, inspect its distance column, then separate the assignment and mean update that define Lloyd's algorithm.

Distance → assignment → mean

Hand-authored, unitless toy observations. Labels identify groups; they have no numerical order.

Filled circles: observations. Hollow diamonds: centroids. Dashed links use held labels. Drag snaps to 0.1; arrow keys preserve full input precision. Only handles capture touch. The selector brings an overlapping handle forward; Enter or Space does the same for a focused handle.

Edit handlesSelected handle

Handle x-5 to 5 · Enter to applyHandle y-5 to 5 · Enter to applyScale x0.1 to 10 · Enter to applyScale y0.1 to 10 · Enter to apply

k1 to 4 · Enter to applySeed0 to 4294967295 · Enter to apply

Frame 1 / 1 · edit begins at the displayed state.

X · samples in rowsEach mathematical feature vector is a 2 × 1 column.: 4 × 2
Centroids / squared distances / assignments: 2 × 2 / 4 × 2 / 4 × 1
Held-assignment WCSSNearest-label preview; press Assign to commit.: ≈ 8
Next nearest-assignment WCSSAllows labels to change; same current centroids.: ≈ 8
Mean update · current → proposed WCSSHeld assignments; empty centroids retained. Inspect before applying Update means.: ≈ 8 → ≈ 4
Phase and stopNot stopped; replay is capped at 40 half-steps.: ready · 0 mean updates
Empty-cluster policyEmpty mean is undefined. Retained coordinates are a policy, not a computed empty mean.: Retain previous centroid
Input provenance / initializationSeed affects partial Fisher–Yates observation selection only. Current centroids may be edited.: Hand-authored observations / seed 7

d_{ij}^2=\sum_{r=1}^2((x_{ir}-c_{jr})/s_r)^2,\quad a_i=\operatorname*{arg\,min}_j d_{ij}^2

\mu_j=\frac{1}{|C_j|}\sum_{i\in C_j}x_i,\quad J=\sum_i d_{i,a_i}^2

Actual scaled squared distances. Lowest centroid index wins exact ties. Preview labels; ≈ marks display rounding.
Observation	x₁	d² to c0	d² to c1	Held label	Nearest / tie
s1	-3	≈ 0	≈ 16	0	0 / unique
s2	-1	≈ 4	≈ 4	0	0 / tie
s3	1	≈ 16	≈ 0	1	1 / unique
s4	3	≈ 36	≈ 4	1	1 / unique

Held membership → sums → proposed mean. Empty clusters retain previous coordinates; their mean is undefined.
Cluster	IDs / n	Sum	Mean / policy	Current centroid
c0	s1, s2 / 2	[-4, 0]	≈ [-2, 0]	≈ [-3, 0]
c1	s3, s4 / 2	[4, 0]	≈ [2, 0]	≈ [1, 0]

Drag a hollow centroid diamond or use its arrow keys and numeric coordinates. Read distance → assignment → sum → mean. The proposed objective is visible before the update; every display uses the active numerical state.

Let the center make a claim

Reset, then select centroid 0. Move its x coordinate from −3 to −2. Only one distance column changes, but a nearest-center assignment can change anywhere its competitors were close. The plotted dashed links and labels show the current assignment; the table records every squared distance. The original point at −1 is tied between centers −3 and 1. Our exact tie chooses center 0.

Press Assign nearest to commit the displayed preview. Then read the cluster sums and the current → proposed WCSS before pressing Update means. Assigning points and moving centers are different operations. Repeat them, or replay their immutable half-step states. There is no animation that secretly fits a different model.

Define the unsupervised problem

We observe feature vectors but have no target labels. A cluster label is an algorithm’s grouping decision. It is not an observed class, a predicted regression target, or an ordered numerical quantity.

For $n$ observations and chosen $k$ , k-means uses centers $c_j$ and assignments $a_i$ . This sheet measures squared scaled Euclidean distance:

d_{ij}^2=\sum_{r=1}^2\left(\frac{x_{ir}-c_{jr}}{s_r}\right)^2, \qquad s_r>0,\qquad J(a,c)=\sum_{i=1}^n d_{i,a_i}^2.

$J$ is the within-cluster sum of squares, WCSS. It is a sum, without $1/n$ , one half, or a square root. Each observation has weight one, including coincident observations with different IDs. Reference scales remain fixed during a Lloyd run. With unit scales, this is ordinary Euclidean k-means. There is no intercept or target column.

Mathematical observations are $2\times1$ columns; $X$ stores samples in its $n\times2$ rows. The center matrix is $k\times2$ , the squared-distance matrix $n\times k$ , and assignments an $n\times1$ integer column. Zero-based cluster indices run from 0 to $k-1$ .

Two exact subproblems give the mechanics

With centers fixed, choose the least costly center independently for every row:

a_i\leftarrow\operatorname*{arg\,min}_{0\le j<k}d_{ij}^2.

The browser and NumPy choose the lowest center index on an exact computed tie. Rounding the table to six significant digits does not create a tie in the model.

With assignments fixed, let $C_j$ contain the observations assigned to center $j$ . For a nonempty cluster, differentiate each coordinate of its quadratic cost:

\frac{\partial J_j}{\partial c_{jr}} =\frac{2}{s_r^2}\left(|C_j|c_{jr}-\sum_{i\in C_j}x_{ir}\right)=0, \qquad c_{jr}\leftarrow\frac{1}{|C_j|}\sum_{i\in C_j}x_{ir}.

Thus the coordinatewise arithmetic mean minimizes that cluster’s cost. Fixed feature scaling changes memberships but does not change the raw-coordinate mean of a fixed membership. The positive quadratic curvature gives the unique mean for a nonempty cluster.

The equivalent identity is useful for an independent proof:

\sum_{i\in C_j}\left\|S^{-1}(x_i-c_j)\right\|^2 =\sum_{i\in C_j}\left\|S^{-1}(x_i-\mu_j)\right\|^2 +|C_j|\left\|S^{-1}(c_j-\mu_j)\right\|^2.

The mean removes the last nonnegative term. Assignment and update cannot increase $J$ in exact arithmetic under these fixed data, scales and declared empty-cluster conventions. Arbitrary dragging changes centers and can increase it.

Calculate the first full iteration literally

The four hand-authored observations have x coordinates −3, −1, 1 and 3, all with y=0. Initial centers are (−3,0) and (1,0). Neither the observations nor these default centers were sampled by a seed.

ID	x₁	d² to c0	d² to c1	Assigned label
s1	−3	0	16	0
s2	−1	4	4	0 by tie
s3	1	16	0	1
s4	3	36	4	1

The objective is $0+4+0+4=8$ . Cluster 0 contains s1/s2, sum (−4,0), count 2, mean (−2,0). Cluster 1 contains s3/s4, sum (4,0), count 2, mean (2,0). At those proposed centers each assigned squared distance is 1, so proposed WCSS is 4. Pressing Update means commits exactly those values.

The next assignment is unchanged. A second mean update also leaves the centers unchanged, so the sheet stops at binary64 equality. Its trace objectives are 8 → 8 → 4 → 4 → 4, including the initial preview and each assignment/update phase. Reaching this fixed point does not generally prove global optimality.

Explore without changing history

The centroid selector, numeric fields, pointer dragging and arrow keys all call the same coordinate edit. Pointer input snaps to a stated 0.1 lattice; keyboard movement adds 0.1, or 0.5 with Shift, to the unrounded displayed coordinate. Numeric Enter preserves fine precision. Values outside −5…5 are visibly rejected; pointer movement stops at that boundary.

Switch Edit handles to observations to move the data themselves. Editing data, centers or scales clears committed assignments and returns to a nearest-label preview. This prevents a historical assignment from masquerading as a new optimization run. Seed changes alone leave the current centers unchanged; Initialize from seed selects new initial observation rows.

Replay records deeply detached, frozen complete states. First, Previous and Next select them; they do not reverse mutations. Edit a center while viewing a historical frame and the new experiment begins with that frame’s centers, observations and scales. Reset restores the literal original data and selected handle.

Empty means, precision and evaluation

An empty cluster has no mean. The policy is retain its previous center. The membership table says “Undefined; retain” and shows the coordinates; it never reports the retained value as an observed empty-cluster mean. Such a center contributes zero to held WCSS. Part V deliberately does not reseed it behind the learner’s back.

WCSS has squared normalized-feature units. Smaller WCSS is not classification accuracy, truth, fairness, or usefulness. With distinct points, allowing $k=n$ can make WCSS zero by assigning each point its own center; that alone does not justify so many clusters. Choosing k requires the intended task, stability and an appropriate external or structural evaluation.

The physical plot uses equal raw x/y scales. Scaling changes the metric without redrawing a misleading transformed coordinate system. Numeric outputs are authoritative where observations coincide. Visible labels are sized for the measured viewport; wide tables and equations scroll in named bounded regions.

Read the actual canonical code

The source below implements inspect_lloyd, lloyd_step, lloyd_trace and deterministic initialization directly with NumPy arrays and explicit boundary checks. From the book root, the default invocation prints the literal first mean update:

.venv/bin/python -B code/python/p05/clustering.py

code/python/p05/clustering.py

"""Canonical from-scratch NumPy clustering. X is n×2; centers are k×2.
No target column, fitting library, hidden reseeding, or optimum guarantee.
"""
import json
import math
import sys
from copy import deepcopy
import numpy as np


def number(v, name, lo=None, hi=None, integer=False):
    try:
        valid = isinstance(v, (int, float)) and not isinstance(v, bool) and math.isfinite(float(v))
    except (OverflowError, TypeError):
        valid = False
    if not valid:
        raise ValueError(f"{name} must be a finite number.")
    v = float(v)
    if (lo is not None and v < lo) or (hi is not None and v > hi) or (integer and v != math.floor(v)):
        raise ValueError(f"{name} is outside its supported range.")
    return int(v) if integer else v


def vector(v, name):
    if not isinstance(v, (list, tuple)) or len(v) != 2:
        raise ValueError(f"{name} needs a dense two-feature row.")
    return [number(x, name) for x in v]


def observations(points):
    if not isinstance(points, list) or not 1 <= len(points) <= 32:
        raise ValueError("Use 1–32 observations; empty data have no objective.")
    ids, rows = set(), []
    for p in points:
        if not isinstance(p, dict) or not isinstance(p.get("id"), str) or not p["id"].strip() or p["id"] in ids:
            raise ValueError("Observation IDs must be unique nonempty strings.")
        ids.add(p["id"])
        rows.append(vector(p.get("x"), "Observation"))
    return np.asarray(rows, dtype=np.float64)


def scales_vector(scales):
    values = vector(scales, "Scales")
    if any(x <= 0 for x in values):
        raise ValueError("Reference scales must be positive.")
    return values


def squared_distance(a, b, scales=(1, 1)):
    a, b, scales = vector(a, "First vector"), vector(b, "Second vector"), scales_vector(scales)
    total = 0.0
    for j in range(2):
        difference = number(a[j]-b[j], "Coordinate difference")
        delta = number(difference/scales[j], "Scaled difference")
        square = number(delta*delta, "Squared difference")
        if difference != 0 and (delta == 0 or square == 0):
            raise ValueError("Distance underflow: rescale coordinates or reference scales so distinct points remain distinguishable.")
        total = number(total+square, "Squared distance")
    return total


def distance_matrix(X, centers, scales):
    # Explicit scalar boundary checks keep NumPy warning settings from changing semantics.
    return np.asarray([[squared_distance(x.tolist(), c.tolist(), scales) for c in centers] for x in X], dtype=np.float64)


def inspect_lloyd(state):
    X = observations(state["points"])
    scales = scales_vector(state["scales"])
    centers = state["centers"]
    if not isinstance(centers, list):
        raise ValueError("Centroids must be a dense array.")
    number(len(centers), "k", 1, len(X), True)
    C = np.asarray([vector(c, "Centroid") for c in centers], dtype=np.float64)
    phase, labels = state["phase"], state["labels"]
    if phase not in ("ready", "assign", "update") or type(state["stopped"]) is not bool:
        raise ValueError("Invalid phase or stop flag.")
    number(state["iteration"], "Iteration", 0, 1000, True)
    if phase == "ready":
        if labels is not None or state["stopped"] or state["iteration"] != 0:
            raise ValueError("Invalid ready state.")
    elif not isinstance(labels, list) or len(labels) != len(X):
        raise ValueError("Assignments must match rows.")
    if state["stopped"] and phase != "update":
        raise ValueError("Stop is reported only after mean update.")
    if labels is not None:
        labels = [number(v, "Assignment", 0, len(C)-1, True) for v in labels]
    D = distance_matrix(X, C, scales)
    nearest = np.argmin(D, axis=1).tolist()  # first index on computed equality
    held = nearest if labels is None else labels
    clusters = []
    for j, c in enumerate(C):
        indices = [i for i, label in enumerate(held) if label == j]
        sums = [0.0, 0.0]
        for i in indices:
            for d in range(2):
                sums[d] = number(sums[d]+float(X[i, d]), "Cluster sum")
        mean = [number(v/len(indices), "Mean") for v in sums] if indices else c.tolist()
        clusters.append({"index":j, "ids":[state["points"][i]["id"] for i in indices], "count":len(indices), "sums":sums, "mean":mean, "empty":not indices})
    def objective(assignment, matrix=D):
        total = 0.0
        for i, label in enumerate(assignment):
            total = number(total+float(matrix[i, label]), "WCSS")
        return total
    return {"distances":D.tolist(), "nearest":nearest, "labels":held, "clusters":clusters,
            "objective":objective(held), "nearestObjective":objective(nearest),
            "proposedObjective":objective(held, distance_matrix(X, np.asarray([c["mean"] for c in clusters]), scales)),
            "ties":[sum(v == row[nearest[i]] for v in row) > 1 for i, row in enumerate(D.tolist())]}


def lloyd_step(state):
    active = inspect_lloyd(state)
    result = deepcopy(state)
    if state["stopped"]:
        return result
    if state["phase"] != "assign":
        result.update(labels=active["nearest"], phase="assign")
    else:
        centers = [c["mean"] for c in active["clusters"]]
        result.update(centers=centers, phase="update", iteration=state["iteration"]+1, stopped=centers == state["centers"])
    inspect_lloyd(result)
    return result


def lloyd_trace(state, half_steps=40):
    number(half_steps, "Half steps", 1, 100, True)
    frames = [deepcopy(state)]
    inspect_lloyd(state)
    for _ in range(half_steps):
        if frames[-1]["stopped"]:
            break
        frames.append(lloyd_step(frames[-1]))
    return frames


def initialize(points, k, seed):
    observations(points)
    k = number(k, "k", 1, len(points), True)
    state = number(seed, "Seed", 0, 2**32-1, True)
    order = list(range(len(points)))
    for i in range(k):
        state = (1664525*state+1013904223) % 2**32
        j = i+math.floor((state/2**32)*(len(points)-i))
        order[i], order[j] = order[j], order[i]
    return [vector(points[i]["x"], "Selected centroid") for i in order[:k]]


def same_partition(a, b):
    if not isinstance(a, list) or not isinstance(b, list) or len(a) != len(b):
        raise ValueError("Partitions need equal row counts.")
    a = [number(x, "Label", 0, 31, True) for x in a]
    b = [number(x, "Label", 0, 31, True) for x in b]
    return all((a[i] == a[j]) == (b[i] == b[j]) for i in range(len(a)) for j in range(i))


def hierarchical(points, linkage, scales=(1, 1)):
    X = observations(points)
    scales_vector(scales)
    if linkage not in ("single", "complete", "average"):
        raise ValueError("Choose single, complete or average linkage.")
    D = np.sqrt(distance_matrix(X, X, scales)).tolist()
    groups = [{"id":i, "members":[i]} for i in range(len(X))]
    merges = []
    while len(groups) > 1:
        best = None
        for i in range(len(groups)):
            for j in range(i+1, len(groups)):
                values = [D[a][b] for a in groups[i]["members"] for b in groups[j]["members"]]
                if linkage == "single":
                    height = min(values)
                elif linkage == "complete":
                    height = max(values)
                else:
                    total = 0.0
                    for d in values:
                        total = number(total+d, "Linkage sum")
                    height = number(total/len(values), "Average linkage")
                if best is None or height < best[0]:
                    best = height, i, j
        height, i, j = best
        left, right = groups[i], groups[j]
        group = {"id":len(X)+len(merges), "members":sorted(left["members"]+right["members"])}
        merges.append({"left":deepcopy(left), "right":deepcopy(right), "group":deepcopy(group), "height":height})
        groups = sorted([g for g in groups if g["id"] not in (left["id"], right["id"])]+[group], key=lambda g:g["id"])
    return {"distances":D, "merges":merges}


def dbscan(points, epsilon, min_points, scales=(1, 1)):
    X = observations(points)
    epsilon = number(epsilon, "Epsilon", 0)
    min_points = number(min_points, "Minimum points", 1, len(X), True)
    limit = number(epsilon*epsilon, "Squared epsilon")
    D = distance_matrix(X, X, scales).tolist()
    neighborhoods = [[i for i, d in enumerate(row) if d <= limit] for row in D]
    core = [len(ns) >= min_points for ns in neighborhoods]
    labels, visited, count = [-1]*len(X), [False]*len(X), 0
    for i in range(len(X)):
        if visited[i]:
            continue
        visited[i] = True
        if not core[i]:
            continue
        labels[i] = count
        queue, queued = list(neighborhoods[i]), set(neighborhoods[i])
        q = 0
        while q < len(queue):
            j = queue[q]
            if not visited[j]:
                visited[j] = True
                if core[j]:
                    for v in neighborhoods[j]:
                        if v not in queued:
                            queue.append(v); queued.add(v)
            if labels[j] == -1:
                labels[j] = count
            q += 1
        count += 1
    return {"distances":D, "neighborhoods":neighborhoods, "core":core, "labels":labels, "count":count,
            "kinds":["core" if core[i] else "noise" if label < 0 else "border" for i, label in enumerate(labels)]}


def dispatch(request):
    operation = request["operation"]
    if operation == "lloyd":
        state = request["state"]
        before = inspect_lloyd(state)
        after_state = lloyd_step(state)
        after = inspect_lloyd(after_state)
        return {"distances":before["distances"], "labels":before["labels"], "means":[c["mean"] for c in before["clusters"]],
                "empty":[c["empty"] for c in before["clusters"]], "objective":before["objective"],
                "proposedObjective":before["proposedObjective"], "afterCenters":after_state["centers"], "afterObjective":after["objective"]}
    if operation == "hierarchical":
        result = hierarchical(request["points"], request["linkage"], request.get("scales", [1, 1]))
        return {"merges":[{"left":m["left"]["members"], "right":m["right"]["members"], "height":m["height"]} for m in result["merges"]]}
    if operation == "dbscan":
        result = dbscan(request["points"], request["epsilon"], request["minPoints"], request.get("scales", [1, 1]))
        return {k:result[k] for k in ("neighborhoods", "core", "labels", "kinds", "count")}
    if operation == "initialize":
        return {"centers":initialize(request["points"], request["k"], request["seed"])}
    raise ValueError("Unknown clustering operation.")


if __name__ == "__main__":
    raw = "" if sys.stdin.isatty() else sys.stdin.read()
    default = {"operation":"lloyd", "state":{"points":[{"id":f"s{i+1}", "x":[x,0]} for i, x in enumerate([-3,-1,1,3])], "centers":[[-3,0],[1,0]], "scales":[1,1], "labels":[0,0,1,1], "phase":"assign", "iteration":0, "stopped":False}}
    try:
        print(json.dumps(dispatch(json.loads(raw) if raw.strip() else default), allow_nan=False, indent=2))
    except (ValueError, KeyError, TypeError, OverflowError) as error:
        print(f"Clustering input error: {error}", file=sys.stderr)
        raise SystemExit(2)

Independent test source evaluates expanded rational quadratics and the pairwise-distance identity for cluster SSE, and enumerates every tiny unlabeled partition. It imports and tests this actual canonical module; agreement between Python and TypeScript alone is insufficient evidence. Their purpose is to distinguish mathematical evidence from two implementations making the same mistake.

Five graded investigations

Connections and checkpoints

You can now calculate an entire distance matrix, distinguish preview from committed labels, derive a mean update and account for both objective values. III-27 tests the role of initialization and feature units; III-28 changes the definition of a cluster. Book II’s norms, means and partial derivatives supply the mathematical tools. A cluster number describes membership, so relabeling whole groups must not change an evaluation of the partition.