k-Means: A Bad Fixed Point

Math 123: Mathematical Aspects of Data Analysis — Fall 2026

Published

October 15, 2026

This example takes apart the bad fixed point of Lloyd’s algorithm: four points at the corners of an $L\times1$ rectangle and $k=2$. The optimal clustering is left versus right, with cost $1$, but Lloyd’s algorithm can stop at bottom versus top, with cost $L^2$. We look at both fixed points, at which starting centers lead to which one, and at how often each seeding rule falls into the trap.

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

L = 3.0                                 # width of the rectangle; any L > 1 works
X = np.array([[0, 0], [0, 1], [L, 0], [L, 1]], dtype=float)   # x_1, x_2, x_3, x_4
n, k = X.shape[0], 2
NAMES = ["x1", "x2", "x3", "x4"]

Figure style

# A muted ink colour and a translucent grid read correctly against both the
# light and the dark version of this site, so the figures need no theme switch.
INK = "#7c7c78"
GRID = "rgba(128, 128, 128, 0.28)"
# Dark bold text with a white halo stays legible over the points in both themes.
LABEL_FONT = dict(color="#1a1a19", size=13, weight="bold",
                  shadow="0 0 3px white, 0 0 3px white, 0 0 2px white")
CLUSTER_COLOR = ["#2a78d6", "#eb6834"]
# One colour per place Lloyd's algorithm can end up.
OUTCOME_COLOR = {"optimum": "#1baf7a", "bad": "#e34948", "one cluster": "#6250d6"}
OUTCOME_NAME = {"optimum": "left / right (cost 1)",
                "bad": f"bottom / top (cost L² = {L**2:g})",
                "one cluster": f"one empty cluster (cost L² + 1 = {L**2 + 1:g})"}


def style(fig, height=420, title=None, showlegend=None):
    """Apply the shared layout to a plot with axes."""
    fig.update_layout(
        template="simple_white",
        paper_bgcolor="rgba(0,0,0,0)",
        plot_bgcolor="rgba(0,0,0,0)",
        font=dict(color=INK, size=13),
        title=dict(text=title, x=0, xanchor="left", font=dict(size=15)),
        margin=dict(l=70, r=30, t=60 if title else 30, b=55),
        height=height,
        showlegend=showlegend,
        hovermode="closest",
    )
    fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    return fig


def equal_aspect(fig, n_panels):
    """Make one unit the same length on both axes of every subplot."""
    for i in range(1, n_panels + 1):
        suffix = "" if i == 1 else str(i)
        fig.layout[f"yaxis{suffix}"].scaleanchor = f"x{suffix}"
        fig.layout[f"xaxis{suffix}"].constrain = "domain"

Two fixed points

This is the lloyd function from the “k-means in practice” box in Lecture 12, changed only to also report how many update steps it took.

def lloyd(X, centers, max_iter=100):
    centers = np.array(centers, dtype=float)
    labels = None
    steps = 0
    for _ in range(max_iter):
        dists = ((X[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)  # (n, k)
        new_labels = dists.argmin(axis=1)          # assignment step
        if labels is not None and np.array_equal(new_labels, labels):
            break
        labels = new_labels
        steps += 1
        for l in range(len(centers)):              # update step
            if np.any(labels == l):
                centers[l] = X[labels == l].mean(axis=0)
    cost = ((X - centers[labels]) ** 2).sum()
    return labels, centers, cost, steps

Start once from the two bad centers $(\tfrac L2,0)$ and $(\tfrac L2,1)$ of the lecture, and once from two points near the middles of the short sides.

starts = {
    "bad": np.array([[L / 2, 0], [L / 2, 1]]),
    "optimum": np.array([[0.4, 0.5], [L - 0.4, 0.5]]),
}
runs = {name: lloyd(X, c) for name, c in starts.items()}

for name, (labels, centers, cost, steps) in runs.items():
    print(f"{name:8s} labels {labels}, centers {centers.tolist()}, "
          f"cost {cost:g} after {steps} update step(s)")
print(f"\nratio of the two costs: {runs['bad'][2] / runs['optimum'][2]:g} = L^2")

bad      labels [0 1 0 1], centers [[1.5, 0.0], [1.5, 1.0]], cost 9 after 1 update step(s)
optimum  labels [0 0 1 1], centers [[0.0, 0.5], [3.0, 0.5]], cost 1 after 1 update step(s)

ratio of the two costs: 9 = L^2

From the bad start the first update step reproduces the starting centers exactly, so the assignment never changes and the algorithm stops with cost $L^2$. Both final configurations satisfy part 3 of the termination theorem: every center is the mean of its cluster and every point is at a nearest center. Nothing in the stopping rule distinguishes the good one from the bad one.

Plot both fixed points

xlim, ylim = (-0.75, L + 0.75), (-0.75, 1.75)
gx, gy = np.meshgrid(np.linspace(*xlim, 300), np.linspace(*ylim, 150))
grid = np.c_[gx.ravel(), gy.ravel()]
two_colors = [[0, CLUSTER_COLOR[0]], [0.5, CLUSTER_COLOR[0]],
              [0.5, CLUSTER_COLOR[1]], [1, CLUSTER_COLOR[1]]]

fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.08, subplot_titles=[
    f"Optimum: left / right, cost {runs['optimum'][2]:g}",
    f"Bad fixed point: bottom / top, cost {runs['bad'][2]:g}"])
for col, name in enumerate(["optimum", "bad"], start=1):
    labels, centers, cost, _ = runs[name]
    cells = ((grid[:, None, :] - centers[None]) ** 2).sum(axis=2).argmin(axis=1)
    fig.add_trace(go.Heatmap(x=gx[0], y=gy[:, 0], z=cells.reshape(gx.shape),
                             colorscale=two_colors, zmin=0, zmax=1, opacity=0.18,
                             showscale=False, hoverinfo="skip"), row=1, col=col)
    fig.add_trace(go.Scatter(x=X[:, 0], y=X[:, 1], mode="markers+text", showlegend=False,
                             text=NAMES, textposition="middle left", textfont=LABEL_FONT,
                             marker=dict(size=14, color=[CLUSTER_COLOR[l] for l in labels],
                                         line=dict(color="white", width=2)),
                             customdata=labels,
                             hovertemplate="%{text}: cluster %{customdata}<extra></extra>"),
                  row=1, col=col)
    fig.add_trace(go.Scatter(x=centers[:, 0], y=centers[:, 1], mode="markers", showlegend=False,
                             marker=dict(symbol="x", size=14, color=CLUSTER_COLOR,
                                         line=dict(color="white", width=1.5)),
                             hovertemplate="center (%{x:.2f}, %{y:.2f})<extra></extra>"),
                  row=1, col=col)
fig.update_xaxes(range=xlim, title="first coordinate")
fig.update_yaxes(range=ylim)
fig.update_yaxes(title="second coordinate", col=1)
equal_aspect(fig, 2)
style(fig, height=380)
fig.show()

Figure 1: The two fixed points of Lloyd’s algorithm for $L=3$. Shading shows the nearest-center regions, split by the perpendicular bisector of the two centers; crosses are the centers. Left: the optimum, cost $1$. Right: the bad fixed point, cost $L^2=9$.

Stuck in a basin

A fixed point that the slightest nudge escapes would be harmless in practice. This one is not like that. Move both bad centers by up to $0.1$ in each coordinate, at random, and run Lloyd’s algorithm again.

rng = np.random.default_rng(0)
nudged_costs = [lloyd(X, starts["bad"] + rng.uniform(-0.1, 0.1, size=(2, 2)))[2]
                for _ in range(1000)]
print("final costs after 1000 nudges:", np.unique(np.round(nudged_costs, 9)))

final costs after 1000 nudges: [9.]

Every run returns to cost $L^2$. As long as both bottom corners stay closer to the lower center and both top corners to the upper one, the first assignment is still bottom versus top, and the first update step snaps the centers back to $(\tfrac L2,0)$ and $(\tfrac L2,1)$. The bad fixed point is a genuine local minimum with a whole region of starting centers around it, its basin of attraction.

Two initial centers make four numbers, so we cannot draw the basin directly, but we can draw slices of it: fix the first center $\bfmu_1$, and color each possible position of the second center $\bfmu_2$ by where Lloyd’s algorithm ends up.

def outcome(cost):
    """Name the fixed point that a final cost belongs to."""
    if np.isclose(cost, 1):
        return "optimum"
    if np.isclose(cost, L**2):
        return "bad"
    if np.isclose(cost, L**2 + 1):
        return "one cluster"
    return "other"


bx, by = np.meshgrid(np.linspace(-1.5, L + 1.5, 161), np.linspace(-2, 3, 121))
fixed_mu1 = {"μ1 = (L/2, 0), a bad center": (L / 2, 0),
             "μ1 = x1 = (0, 0), a corner": (0, 0)}
basins = {}
for title, mu1 in fixed_mu1.items():
    basins[title] = np.array([outcome(lloyd(X, [mu1, (a, b)])[2])
                              for a, b in zip(bx.ravel(), by.ravel())]).reshape(bx.shape)
    names, counts = np.unique(basins[title], return_counts=True)
    print(title, {str(nm): f"{c / counts.sum():.1%}" for nm, c in zip(names, counts)})

μ1 = (L/2, 0), a bad center {'bad': '2.2%', 'one cluster': '23.7%', 'optimum': '74.0%'}
μ1 = x1 = (0, 0), a corner {'bad': '2.1%', 'one cluster': '8.9%', 'optimum': '88.9%'}

Plot the basin slices

order = ["optimum", "bad", "one cluster"]
colorscale = [[edge, OUTCOME_COLOR[nm]] for i, nm in enumerate(order)
              for edge in (i / len(order), (i + 1) / len(order))]

fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.08, subplot_titles=list(fixed_mu1))
for col, (title, mu1) in enumerate(fixed_mu1.items(), start=1):
    z = np.vectorize(order.index)(basins[title])
    fig.add_trace(go.Heatmap(x=bx[0], y=by[:, 0], z=z, colorscale=colorscale,
                             zmin=-0.5, zmax=len(order) - 0.5, opacity=0.55, showscale=False,
                             text=np.vectorize(OUTCOME_NAME.get)(basins[title]),
                             hovertemplate="μ2 = (%{x:.2f}, %{y:.2f})<br>%{text}<extra></extra>"),
                  row=1, col=col)
    fig.add_trace(go.Scatter(x=X[:, 0], y=X[:, 1], mode="markers+text", showlegend=False,
                             text=NAMES, textposition="middle left", textfont=LABEL_FONT,
                             marker=dict(size=10, color=INK, line=dict(color="white", width=2)),
                             hoverinfo="skip"), row=1, col=col)
    fig.add_trace(go.Scatter(x=[mu1[0]], y=[mu1[1]], mode="markers", showlegend=False,
                             marker=dict(symbol="x", size=16, color="#1a1a19",
                                         line=dict(color="white", width=1.5)),
                             hovertemplate="μ1 (fixed)<extra></extra>"), row=1, col=col)
for nm in order:                                   # legend entries for the heatmap colors
    fig.add_trace(go.Scatter(x=[None], y=[None], mode="markers", name=OUTCOME_NAME[nm],
                             marker=dict(symbol="square", size=12, color=OUTCOME_COLOR[nm],
                                         opacity=0.7)))
fig.update_xaxes(range=[bx.min(), bx.max()], title="first coordinate of μ2")
fig.update_yaxes(range=[by.min(), by.max()])
fig.update_yaxes(title="second coordinate of μ2", col=1)
equal_aspect(fig, 2)
style(fig, height=470, showlegend=True)
fig.update_layout(legend=dict(orientation="h", x=0, y=-0.18, xanchor="left", yanchor="top"),
                  margin=dict(b=110))
fig.show()

Figure 2: Two slices of the basins of attraction for $L=3$. The first center $\bfmu_1$ is fixed (black cross) and each pixel is a choice of the second center $\bfmu_2$, colored by the fixed point that Lloyd’s algorithm reaches from $(\bfmu_1,\bfmu_2)$. The four data points are the gray dots.

Three things show up.

The red region has positive area, so the bad fixed point is reached from a whole region of starts, not just one exact position. It is only a thin sliver in these slices, though, so centers placed uniformly at random in the picture would rarely land in it.
There is a third fixed point, in purple: if $\bfmu_2$ starts far enough away that no data point is closest to it, its cluster is empty, $\bfmu_2$ stays where it is, $\bfmu_1$ moves to the mean of all four points, and the algorithm stops with cost $L^2+1$. This is why practical initializations put the centers at data points.
In the right panel, $\bfmu_1$ sits on the corner $\bfx_1$, as in the first step of k-means++. Of the other three corners, only $\bfx_2$, directly above $\bfx_1$, lies in the red region, while $\bfx_3$ and $\bfx_4$ lie in the green one. The sliver is thin, but a data point sits right in it. This answers the “(check this!)” in the lecture: from $\bfmu_1=\bfx_1$ and $\bfmu_2\in\{\bfx_3,\bfx_4\}$, Lloyd’s algorithm finds the left/right optimum.

How often does seeding fall in?

Now restrict the starting centers to two distinct data points, as both standard seedings do. There are only $\binom42=6$ such pairs, so we can try them all.

pairs = [(0, 1), (2, 3), (0, 2), (1, 3), (0, 3), (1, 2)]
pair_runs = {p: lloyd(X, X[list(p)]) for p in pairs}
for (i, j), (labels, centers, cost, steps) in pair_runs.items():
    print(f"start at x{i + 1}, x{j + 1}: final cost {cost:g} ({outcome(cost)})")

start at x1, x2: final cost 9 (bad)
start at x3, x4: final cost 9 (bad)
start at x1, x3: final cost 1 (optimum)
start at x2, x4: final cost 1 (optimum)
start at x1, x4: final cost 1 (optimum)
start at x2, x3: final cost 1 (optimum)

Plot all six starts

fig = make_subplots(rows=2, cols=3, horizontal_spacing=0.05, vertical_spacing=0.16,
                    subplot_titles=[f"start x{i + 1}, x{j + 1}: cost {pair_runs[(i, j)][2]:g}"
                                    for i, j in pairs])
for idx, p in enumerate(pairs):
    row, col = idx // 3 + 1, idx % 3 + 1
    labels, centers, cost, _ = pair_runs[p]
    start = X[list(p)]
    for l in range(k):                             # path of each center
        fig.add_trace(go.Scatter(x=[start[l, 0], centers[l, 0]], y=[start[l, 1], centers[l, 1]],
                                 mode="lines", line=dict(color=CLUSTER_COLOR[l], width=2, dash="dot"),
                                 showlegend=False, hoverinfo="skip"), row=row, col=col)
    fig.add_trace(go.Scatter(x=X[:, 0], y=X[:, 1], mode="markers", showlegend=False,
                             marker=dict(size=11, color=[CLUSTER_COLOR[l] for l in labels],
                                         line=dict(color="white", width=2)),
                             text=NAMES, customdata=labels,
                             hovertemplate="%{text}: cluster %{customdata}<extra></extra>"),
                  row=row, col=col)
    fig.add_trace(go.Scatter(x=start[:, 0], y=start[:, 1], mode="markers", showlegend=False,
                             marker=dict(symbol="circle-open", size=22, color=CLUSTER_COLOR,
                                         line=dict(width=2.5)),
                             hovertemplate="initial center (%{x:g}, %{y:g})<extra></extra>"),
                  row=row, col=col)
    fig.add_trace(go.Scatter(x=centers[:, 0], y=centers[:, 1], mode="markers", showlegend=False,
                             marker=dict(symbol="x", size=13, color=CLUSTER_COLOR,
                                         line=dict(color="white", width=1.5)),
                             hovertemplate="final center (%{x:.2f}, %{y:.2f})<extra></extra>"),
                  row=row, col=col)
fig.update_xaxes(range=[-0.5, L + 0.5], showticklabels=False)
fig.update_yaxes(range=[-0.5, 1.5], showticklabels=False)
equal_aspect(fig, 6)
style(fig, height=420)
fig.update_annotations(font=dict(size=13))
fig.show()

Figure 3: Lloyd’s algorithm from each of the six pairs of data points as initial centers, for $L=3$. Circles are the initial centers, crosses the final ones, and the dashed lines show how far each center moves. Starting from two corners on the same short side (top row, left and middle) leads to the bad fixed point.

Only the two vertical pairs, $\{\bfx_1,\bfx_2\}$ and $\{\bfx_3,\bfx_4\}$, are bad. With both centers on the same short side, their bisector is horizontal, so it splits bottom from top. This gives the failure probabilities of the two seedings exactly.

Random data points: the pair is uniform among the $6$, so Lloyd’s algorithm fails with probability $\tfrac26=\tfrac13$, whatever $L$ is. As $L$ grows the failure gets worse (its cost is $L^2$) but no rarer, and the expected cost $\tfrac23\cdot1+\tfrac13L^2$ grows without bound.
k-means++: whichever corner is chosen first, the bad second choice is its vertical neighbor, at squared distance $1$, against $L^2$ and $L^2+1$ for the other two corners. So it fails with probability $1/(2L^2+2)$, and the expected cost is $1+(L^2-1)/(2L^2+2)<\tfrac32$ for every $L$.

To check these numbers, here is k-means++ seeding written out, followed by a simulation of both seedings for several widths $L$.

def kmeans_pp(X, k, rng):
    """Choose k initial centers from the rows of X by k-means++ seeding."""
    centers = [X[rng.integers(len(X))]]
    for _ in range(1, k):
        D2 = ((X[:, None, :] - np.array(centers)[None]) ** 2).sum(axis=2).min(axis=1)
        centers.append(X[rng.choice(len(X), p=D2 / D2.sum())])
    return np.array(centers)


def random_points(X, k, rng):
    """Choose k distinct rows of X uniformly at random."""
    return X[rng.choice(len(X), size=k, replace=False)]


def failure_rate(seeding, L, trials, rng):
    """Fraction of seeded Lloyd runs on the L x 1 rectangle that miss the optimum."""
    XL = np.array([[0, 0], [0, 1], [L, 0], [L, 1]], dtype=float)
    return np.mean([lloyd(XL, seeding(XL, 2, rng))[2] > 1 + 1e-9 for _ in range(trials)])


rng = np.random.default_rng(1)
widths = np.array([1.25, 1.5, 2, 3, 4, 6, 8])
trials = 20_000
sim = {"random data points": [failure_rate(random_points, w, trials, rng) for w in widths],
       "k-means++": [failure_rate(kmeans_pp, w, trials, rng) for w in widths]}

print("    L   random (sim / exact)   k-means++ (sim / exact)")
for w, r, p in zip(widths, sim["random data points"], sim["k-means++"]):
    print(f"{w:5.2f}   {r:.4f} / {1/3:.4f}        {p:.4f} / {1 / (2 * w**2 + 2):.4f}")

    L   random (sim / exact)   k-means++ (sim / exact)
 1.25   0.3361 / 0.3333        0.1989 / 0.1951
 1.50   0.3329 / 0.3333        0.1515 / 0.1538
 2.00   0.3362 / 0.3333        0.0980 / 0.1000
 3.00   0.3311 / 0.3333        0.0512 / 0.0500
 4.00   0.3313 / 0.3333        0.0292 / 0.0294
 6.00   0.3371 / 0.3333        0.0136 / 0.0135
 8.00   0.3318 / 0.3333        0.0081 / 0.0077

Plot the failure probabilities

Ls = np.linspace(1, widths.max(), 200)
exact = {"random data points": np.full_like(Ls, 1 / 3), "k-means++": 1 / (2 * Ls**2 + 2)}

fig = go.Figure()
for color, name in zip(CLUSTER_COLOR, exact):
    fig.add_trace(go.Scatter(x=Ls, y=exact[name], mode="lines", name=f"{name} (exact)",
                             line=dict(color=color, width=2), legendgroup=name,
                             hovertemplate="L = %{x:.2f}<br>exact %{y:.4f}<extra>" + name + "</extra>"))
    fig.add_trace(go.Scatter(x=widths, y=sim[name], mode="markers", name=f"{name} (simulated)",
                             marker=dict(color=color, size=9, line=dict(color="white", width=2)),
                             legendgroup=name,
                             hovertemplate="L = %{x:.2f}<br>simulated %{y:.4f}<extra>" + name + "</extra>"))
    fig.add_annotation(x=Ls[-1], y=np.log10(exact[name][-1]), text=name, showarrow=False,
                       xanchor="right", yanchor="bottom", yshift=4, font=dict(color=INK, size=13))
fig.update_xaxes(title="width L of the rectangle", range=[1, widths.max() + 0.2])
ticks = [0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5]
fig.update_yaxes(title="probability of missing the optimum", type="log",
                 range=[np.log10(0.004), np.log10(0.6)],
                 tickvals=ticks, ticktext=[f"{t:g}" for t in ticks])
style(fig, height=430, showlegend=True)
fig.update_layout(legend=dict(x=0.02, y=0.02, xanchor="left", yanchor="bottom"))
fig.show()

Figure 4: Probability that Lloyd’s algorithm misses the optimum on the $L\times1$ rectangle, for the two seedings. Lines are the exact values $\tfrac13$ and $1/(2L^2+2)$; dots are the simulated rates from $20{,}000$ runs each. The vertical axis is logarithmic.

The simulations agree with the exact values to within sampling error. k-means++ helps exactly where the problem is worst: the wider the rectangle, the bigger the disaster at the bad fixed point, and the less likely k-means++ is to land there.

Restarts in scikit-learn

sklearn.cluster.KMeans uses both remedies at once. With init="random" it seeds from random data points, and with n_init it keeps the best of that many runs, so it fails only if every run fails, with probability $(\tfrac13)^{\texttt{n\_init}}$ here.

from sklearn.cluster import KMeans

seeds = range(300)
for init, n_init in [("random", 1), ("random", 3), ("random", 10), ("k-means++", 1)]:
    costs = [KMeans(n_clusters=2, init=init, n_init=n_init, random_state=s).fit(X).inertia_
             for s in seeds]
    rate = np.mean(np.array(costs) > 1 + 1e-9)
    print(f'init="{init}", n_init={n_init:2d}: missed the optimum in {rate:.1%} of {len(seeds)} seeds')

init="random", n_init= 1: missed the optimum in 35.3% of 300 seeds
init="random", n_init= 3: missed the optimum in 2.3% of 300 seeds
init="random", n_init=10: missed the optimum in 0.0% of 300 seeds
init="k-means++", n_init= 1: missed the optimum in 0.3% of 300 seeds

The rates for init="random" fall roughly like $\tfrac13$, $\tfrac1{27}$, and $3^{-10}$. scikit-learn’s "k-means++" fails even less often than the $1/(2L^2+2)=5\%$ computed above, because it is a greedy variant: for each new center it draws a few candidates by the $D^2$ rule (two, when $k=2$) and keeps the one that lowers the cost the most. It picks the bad vertical neighbor only if every candidate is that neighbor.

Even so, no seeding rule makes Lloyd’s algorithm exact. The bad fixed point is still there, and its basin (Figure 2) still contains two of the six pairs of data points; good seeding just makes it unlikely that we start inside it.

--- title: "k-Means: A Bad Fixed Point" subtitle: "Math 123: Mathematical Aspects of Data Analysis — Fall 2026" date: 2026-10-15 number-sections: false jupyter: python3 format: html: toc: true toc-depth: 2 code-line-numbers: true code-tools: true out-width: "100%" column: page execute: echo: true draft: false --- This example takes apart the [bad fixed point](../../lectures/lecture-12.qmd#exm-kmeans-bad-local-min) of [Lloyd's algorithm](../../lectures/lecture-12.qmd#alg-lloyd): four points at the corners of an $L\times1$ rectangle and $k=2$. The optimal clustering is left versus right, with cost $1$, but Lloyd's algorithm can stop at bottom versus top, with cost $L^2$. We look at both fixed points, at which starting centers lead to which one, and at how often each seeding rule falls into the trap. ```{python} import numpy as np import plotly.graph_objects as go from plotly.subplots import make_subplots L = 3.0 # width of the rectangle; any L > 1 works X = np.array([[0, 0], [0, 1], [L, 0], [L, 1]], dtype=float) # x_1, x_2, x_3, x_4 n, k = X.shape[0], 2 NAMES = ["x1", "x2", "x3", "x4"] ``` ```{python} #| code-fold: true #| code-summary: "Figure style" # A muted ink colour and a translucent grid read correctly against both the # light and the dark version of this site, so the figures need no theme switch. INK = "#7c7c78" GRID = "rgba(128, 128, 128, 0.28)" # Dark bold text with a white halo stays legible over the points in both themes. LABEL_FONT = dict(color="#1a1a19", size=13, weight="bold", shadow="0 0 3px white, 0 0 3px white, 0 0 2px white") CLUSTER_COLOR = ["#2a78d6", "#eb6834"] # One colour per place Lloyd's algorithm can end up. OUTCOME_COLOR = {"optimum": "#1baf7a", "bad": "#e34948", "one cluster": "#6250d6"} OUTCOME_NAME = {"optimum": "left / right (cost 1)", "bad": f"bottom / top (cost L² = {L**2:g})", "one cluster": f"one empty cluster (cost L² + 1 = {L**2 + 1:g})"} def style(fig, height=420, title=None, showlegend=None): """Apply the shared layout to a plot with axes.""" fig.update_layout( template="simple_white", paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)", font=dict(color=INK, size=13), title=dict(text=title, x=0, xanchor="left", font=dict(size=15)), margin=dict(l=70, r=30, t=60 if title else 30, b=55), height=height, showlegend=showlegend, hovermode="closest", ) fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) return fig def equal_aspect(fig, n_panels): """Make one unit the same length on both axes of every subplot.""" for i in range(1, n_panels + 1): suffix = "" if i == 1 else str(i) fig.layout[f"yaxis{suffix}"].scaleanchor = f"x{suffix}" fig.layout[f"xaxis{suffix}"].constrain = "domain" ``` ## Two fixed points {#sec-fixed-points} This is the `lloyd` function from the "k-means in practice" box in [Lecture 12](../../lectures/lecture-12.qmd#alg-lloyd), changed only to also report how many update steps it took. ```{python} def lloyd(X, centers, max_iter=100): centers = np.array(centers, dtype=float) labels = None steps = 0 for _ in range(max_iter): dists = ((X[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2) # (n, k) new_labels = dists.argmin(axis=1) # assignment step if labels is not None and np.array_equal(new_labels, labels): break labels = new_labels steps += 1 for l in range(len(centers)): # update step if np.any(labels == l): centers[l] = X[labels == l].mean(axis=0) cost = ((X - centers[labels]) ** 2).sum() return labels, centers, cost, steps ``` Start once from the two bad centers $(\tfrac L2,0)$ and $(\tfrac L2,1)$ of the lecture, and once from two points near the middles of the short sides. ```{python} starts = { "bad": np.array([[L / 2, 0], [L / 2, 1]]), "optimum": np.array([[0.4, 0.5], [L - 0.4, 0.5]]), } runs = {name: lloyd(X, c) for name, c in starts.items()} for name, (labels, centers, cost, steps) in runs.items(): print(f"{name:8s} labels {labels}, centers {centers.tolist()}, " f"cost {cost:g} after {steps} update step(s)") print(f"\nratio of the two costs: {runs['bad'][2] / runs['optimum'][2]:g} = L^2") ``` From the bad start the first update step reproduces the starting centers exactly, so the assignment never changes and the algorithm stops with cost $L^2$. Both final configurations satisfy part 3 of [the termination theorem](../../lectures/lecture-12.qmd#thm-lloyd): every center is the mean of its cluster and every point is at a nearest center. Nothing in the stopping rule distinguishes the good one from the bad one. ```{python} #| label: fig-fixed-points #| fig-cap: "The two fixed points of Lloyd's algorithm for $L=3$. Shading shows the nearest-center regions, split by the perpendicular bisector of the two centers; crosses are the centers. Left: the optimum, cost $1$. Right: the bad fixed point, cost $L^2=9$." #| fig-alt: "Two plots of the same four points at the corners of a 3 by 1 rectangle. On the left, a vertical line splits the plane into a blue left half and an orange right half; the two left corners are blue, the two right corners are orange, and the centers sit at the midpoints of the short sides. On the right, a horizontal line splits the plane into a blue bottom half and an orange top half; the two bottom corners are blue, the two top corners are orange, and the centers sit at the midpoints of the long sides." #| code-fold: true #| code-summary: "Plot both fixed points" xlim, ylim = (-0.75, L + 0.75), (-0.75, 1.75) gx, gy = np.meshgrid(np.linspace(*xlim, 300), np.linspace(*ylim, 150)) grid = np.c_[gx.ravel(), gy.ravel()] two_colors = [[0, CLUSTER_COLOR[0]], [0.5, CLUSTER_COLOR[0]], [0.5, CLUSTER_COLOR[1]], [1, CLUSTER_COLOR[1]]] fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.08, subplot_titles=[ f"Optimum: left / right, cost {runs['optimum'][2]:g}", f"Bad fixed point: bottom / top, cost {runs['bad'][2]:g}"]) for col, name in enumerate(["optimum", "bad"], start=1): labels, centers, cost, _ = runs[name] cells = ((grid[:, None, :] - centers[None]) ** 2).sum(axis=2).argmin(axis=1) fig.add_trace(go.Heatmap(x=gx[0], y=gy[:, 0], z=cells.reshape(gx.shape), colorscale=two_colors, zmin=0, zmax=1, opacity=0.18, showscale=False, hoverinfo="skip"), row=1, col=col) fig.add_trace(go.Scatter(x=X[:, 0], y=X[:, 1], mode="markers+text", showlegend=False, text=NAMES, textposition="middle left", textfont=LABEL_FONT, marker=dict(size=14, color=[CLUSTER_COLOR[l] for l in labels], line=dict(color="white", width=2)), customdata=labels, hovertemplate="%{text}: cluster %{customdata}<extra></extra>"), row=1, col=col) fig.add_trace(go.Scatter(x=centers[:, 0], y=centers[:, 1], mode="markers", showlegend=False, marker=dict(symbol="x", size=14, color=CLUSTER_COLOR, line=dict(color="white", width=1.5)), hovertemplate="center (%{x:.2f}, %{y:.2f})<extra></extra>"), row=1, col=col) fig.update_xaxes(range=xlim, title="first coordinate") fig.update_yaxes(range=ylim) fig.update_yaxes(title="second coordinate", col=1) equal_aspect(fig, 2) style(fig, height=380) fig.show() ``` ## Stuck in a basin {#sec-basin} A fixed point that the slightest nudge escapes would be harmless in practice. This one is not like that. Move both bad centers by up to $0.1$ in each coordinate, at random, and run Lloyd's algorithm again. ```{python} rng = np.random.default_rng(0) nudged_costs = [lloyd(X, starts["bad"] + rng.uniform(-0.1, 0.1, size=(2, 2)))[2] for _ in range(1000)] print("final costs after 1000 nudges:", np.unique(np.round(nudged_costs, 9))) ``` Every run returns to cost $L^2$. As long as both bottom corners stay closer to the lower center and both top corners to the upper one, the first assignment is still bottom versus top, and the first update step snaps the centers back to $(\tfrac L2,0)$ and $(\tfrac L2,1)$. The bad fixed point is a genuine local minimum with a whole region of starting centers around it, its *basin of attraction*. Two initial centers make four numbers, so we cannot draw the basin directly, but we can draw slices of it: fix the first center $\bfmu_1$, and color each possible position of the second center $\bfmu_2$ by where Lloyd's algorithm ends up. ```{python} def outcome(cost): """Name the fixed point that a final cost belongs to.""" if np.isclose(cost, 1): return "optimum" if np.isclose(cost, L**2): return "bad" if np.isclose(cost, L**2 + 1): return "one cluster" return "other" bx, by = np.meshgrid(np.linspace(-1.5, L + 1.5, 161), np.linspace(-2, 3, 121)) fixed_mu1 = {"μ1 = (L/2, 0), a bad center": (L / 2, 0), "μ1 = x1 = (0, 0), a corner": (0, 0)} basins = {} for title, mu1 in fixed_mu1.items(): basins[title] = np.array([outcome(lloyd(X, [mu1, (a, b)])[2]) for a, b in zip(bx.ravel(), by.ravel())]).reshape(bx.shape) names, counts = np.unique(basins[title], return_counts=True) print(title, {str(nm): f"{c / counts.sum():.1%}" for nm, c in zip(names, counts)}) ``` ```{python} #| label: fig-basins #| fig-cap: "Two slices of the basins of attraction for $L=3$. The first center $\\bfmu_1$ is fixed (black cross) and each pixel is a choice of the second center $\\bfmu_2$, colored by the fixed point that Lloyd's algorithm reaches from $(\\bfmu_1,\\bfmu_2)$. The four data points are the gray dots." #| fig-alt: "Two maps of the plane around the 3 by 1 rectangle, each mostly green. In the left map the first center is at the middle of the bottom side; green fills two large overlapping disks around the left and right sides of the rectangle, purple fills the corners of the map outside them, and a thin red sliver runs vertically from the first center up through the middle of the top side. In the right map the first center is at the bottom-left corner x1; green covers everything to the right, purple fills the far left and the top-left corner, and a thin red wedge rises from x1 through the top-left corner x2, while the two right corners x3 and x4 sit in green." #| code-fold: true #| code-summary: "Plot the basin slices" order = ["optimum", "bad", "one cluster"] colorscale = [[edge, OUTCOME_COLOR[nm]] for i, nm in enumerate(order) for edge in (i / len(order), (i + 1) / len(order))] fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.08, subplot_titles=list(fixed_mu1)) for col, (title, mu1) in enumerate(fixed_mu1.items(), start=1): z = np.vectorize(order.index)(basins[title]) fig.add_trace(go.Heatmap(x=bx[0], y=by[:, 0], z=z, colorscale=colorscale, zmin=-0.5, zmax=len(order) - 0.5, opacity=0.55, showscale=False, text=np.vectorize(OUTCOME_NAME.get)(basins[title]), hovertemplate="μ2 = (%{x:.2f}, %{y:.2f})<br>%{text}<extra></extra>"), row=1, col=col) fig.add_trace(go.Scatter(x=X[:, 0], y=X[:, 1], mode="markers+text", showlegend=False, text=NAMES, textposition="middle left", textfont=LABEL_FONT, marker=dict(size=10, color=INK, line=dict(color="white", width=2)), hoverinfo="skip"), row=1, col=col) fig.add_trace(go.Scatter(x=[mu1[0]], y=[mu1[1]], mode="markers", showlegend=False, marker=dict(symbol="x", size=16, color="#1a1a19", line=dict(color="white", width=1.5)), hovertemplate="μ1 (fixed)<extra></extra>"), row=1, col=col) for nm in order: # legend entries for the heatmap colors fig.add_trace(go.Scatter(x=[None], y=[None], mode="markers", name=OUTCOME_NAME[nm], marker=dict(symbol="square", size=12, color=OUTCOME_COLOR[nm], opacity=0.7))) fig.update_xaxes(range=[bx.min(), bx.max()], title="first coordinate of μ2") fig.update_yaxes(range=[by.min(), by.max()]) fig.update_yaxes(title="second coordinate of μ2", col=1) equal_aspect(fig, 2) style(fig, height=470, showlegend=True) fig.update_layout(legend=dict(orientation="h", x=0, y=-0.18, xanchor="left", yanchor="top"), margin=dict(b=110)) fig.show() ``` Three things show up. - The red region has positive area, so the bad fixed point is reached from a whole region of starts, not just one exact position. It is only a thin sliver in these slices, though, so centers placed uniformly at random in the picture would rarely land in it. - There is a third fixed point, in purple: if $\bfmu_2$ starts far enough away that no data point is closest to it, its cluster is empty, $\bfmu_2$ [stays where it is](../../lectures/lecture-12.qmd#alg-lloyd), $\bfmu_1$ moves to the mean of all four points, and the algorithm stops with cost $L^2+1$. This is why practical initializations put the centers *at data points*. - In the right panel, $\bfmu_1$ sits on the corner $\bfx_1$, as in the first step of [k-means++](../../lectures/lecture-12.qmd#alg-kmeans-pp). Of the other three corners, only $\bfx_2$, directly above $\bfx_1$, lies in the red region, while $\bfx_3$ and $\bfx_4$ lie in the green one. The sliver is thin, but a data point sits right in it. This answers the "(check this!)" in the lecture: from $\bfmu_1=\bfx_1$ and $\bfmu_2\in\{\bfx_3,\bfx_4\}$, Lloyd's algorithm finds the left/right optimum. ## How often does seeding fall in? {#sec-seeding} Now restrict the starting centers to two distinct data points, as both standard seedings do. There are only $\binom42=6$ such pairs, so we can try them all. ```{python} pairs = [(0, 1), (2, 3), (0, 2), (1, 3), (0, 3), (1, 2)] pair_runs = {p: lloyd(X, X[list(p)]) for p in pairs} for (i, j), (labels, centers, cost, steps) in pair_runs.items(): print(f"start at x{i + 1}, x{j + 1}: final cost {cost:g} ({outcome(cost)})") ``` ```{python} #| label: fig-pairs #| fig-cap: "Lloyd's algorithm from each of the six pairs of data points as initial centers, for $L=3$. Circles are the initial centers, crosses the final ones, and the dashed lines show how far each center moves. Starting from two corners on the same short side (top row, left and middle) leads to the bad fixed point." #| fig-alt: "Six small plots of the 3 by 1 rectangle arranged in two rows of three. In the two plots where the initial centers are the two left corners or the two right corners, the final clusters are bottom and top, with the centers at the midpoints of the long sides. In the other four plots, where the initial centers are the two bottom corners, the two top corners, or a pair of opposite corners, the final clusters are left and right, with the centers at the midpoints of the short sides." #| code-fold: true #| code-summary: "Plot all six starts" fig = make_subplots(rows=2, cols=3, horizontal_spacing=0.05, vertical_spacing=0.16, subplot_titles=[f"start x{i + 1}, x{j + 1}: cost {pair_runs[(i, j)][2]:g}" for i, j in pairs]) for idx, p in enumerate(pairs): row, col = idx // 3 + 1, idx % 3 + 1 labels, centers, cost, _ = pair_runs[p] start = X[list(p)] for l in range(k): # path of each center fig.add_trace(go.Scatter(x=[start[l, 0], centers[l, 0]], y=[start[l, 1], centers[l, 1]], mode="lines", line=dict(color=CLUSTER_COLOR[l], width=2, dash="dot"), showlegend=False, hoverinfo="skip"), row=row, col=col) fig.add_trace(go.Scatter(x=X[:, 0], y=X[:, 1], mode="markers", showlegend=False, marker=dict(size=11, color=[CLUSTER_COLOR[l] for l in labels], line=dict(color="white", width=2)), text=NAMES, customdata=labels, hovertemplate="%{text}: cluster %{customdata}<extra></extra>"), row=row, col=col) fig.add_trace(go.Scatter(x=start[:, 0], y=start[:, 1], mode="markers", showlegend=False, marker=dict(symbol="circle-open", size=22, color=CLUSTER_COLOR, line=dict(width=2.5)), hovertemplate="initial center (%{x:g}, %{y:g})<extra></extra>"), row=row, col=col) fig.add_trace(go.Scatter(x=centers[:, 0], y=centers[:, 1], mode="markers", showlegend=False, marker=dict(symbol="x", size=13, color=CLUSTER_COLOR, line=dict(color="white", width=1.5)), hovertemplate="final center (%{x:.2f}, %{y:.2f})<extra></extra>"), row=row, col=col) fig.update_xaxes(range=[-0.5, L + 0.5], showticklabels=False) fig.update_yaxes(range=[-0.5, 1.5], showticklabels=False) equal_aspect(fig, 6) style(fig, height=420) fig.update_annotations(font=dict(size=13)) fig.show() ``` Only the two vertical pairs, $\{\bfx_1,\bfx_2\}$ and $\{\bfx_3,\bfx_4\}$, are bad. With both centers on the same short side, their bisector is horizontal, so it splits bottom from top. This gives the failure probabilities of the two seedings exactly. - *Random data points*: the pair is uniform among the $6$, so Lloyd's algorithm fails with probability $\tfrac26=\tfrac13$, *whatever* $L$ is. As $L$ grows the failure gets worse (its cost is $L^2$) but no rarer, and the expected cost $\tfrac23\cdot1+\tfrac13L^2$ grows without bound. - *k-means++*: whichever corner is chosen first, the bad second choice is its vertical neighbor, at squared distance $1$, against $L^2$ and $L^2+1$ for the other two corners. So it fails with probability $1/(2L^2+2)$, and the expected cost is $1+(L^2-1)/(2L^2+2)<\tfrac32$ for every $L$. To check these numbers, here is [k-means++ seeding](../../lectures/lecture-12.qmd#alg-kmeans-pp) written out, followed by a simulation of both seedings for several widths $L$. ```{python} def kmeans_pp(X, k, rng): """Choose k initial centers from the rows of X by k-means++ seeding.""" centers = [X[rng.integers(len(X))]] for _ in range(1, k): D2 = ((X[:, None, :] - np.array(centers)[None]) ** 2).sum(axis=2).min(axis=1) centers.append(X[rng.choice(len(X), p=D2 / D2.sum())]) return np.array(centers) def random_points(X, k, rng): """Choose k distinct rows of X uniformly at random.""" return X[rng.choice(len(X), size=k, replace=False)] def failure_rate(seeding, L, trials, rng): """Fraction of seeded Lloyd runs on the L x 1 rectangle that miss the optimum.""" XL = np.array([[0, 0], [0, 1], [L, 0], [L, 1]], dtype=float) return np.mean([lloyd(XL, seeding(XL, 2, rng))[2] > 1 + 1e-9 for _ in range(trials)]) rng = np.random.default_rng(1) widths = np.array([1.25, 1.5, 2, 3, 4, 6, 8]) trials = 20_000 sim = {"random data points": [failure_rate(random_points, w, trials, rng) for w in widths], "k-means++": [failure_rate(kmeans_pp, w, trials, rng) for w in widths]} print(" L random (sim / exact) k-means++ (sim / exact)") for w, r, p in zip(widths, sim["random data points"], sim["k-means++"]): print(f"{w:5.2f} {r:.4f} / {1/3:.4f} {p:.4f} / {1 / (2 * w**2 + 2):.4f}") ``` ```{python} #| label: fig-failure #| fig-cap: "Probability that Lloyd's algorithm misses the optimum on the $L\\times1$ rectangle, for the two seedings. Lines are the exact values $\\tfrac13$ and $1/(2L^2+2)$; dots are the simulated rates from $20{,}000$ runs each. The vertical axis is logarithmic." #| fig-alt: "A line chart of failure probability against the width L, from 1 to 8, on a logarithmic vertical axis. The random data points seeding is a flat line at one third. The k-means++ seeding starts near one quarter at L equal to 1 and falls steadily to below one percent at L equal to 8. Simulated dots sit on both lines." #| code-fold: true #| code-summary: "Plot the failure probabilities" Ls = np.linspace(1, widths.max(), 200) exact = {"random data points": np.full_like(Ls, 1 / 3), "k-means++": 1 / (2 * Ls**2 + 2)} fig = go.Figure() for color, name in zip(CLUSTER_COLOR, exact): fig.add_trace(go.Scatter(x=Ls, y=exact[name], mode="lines", name=f"{name} (exact)", line=dict(color=color, width=2), legendgroup=name, hovertemplate="L = %{x:.2f}<br>exact %{y:.4f}<extra>" + name + "</extra>")) fig.add_trace(go.Scatter(x=widths, y=sim[name], mode="markers", name=f"{name} (simulated)", marker=dict(color=color, size=9, line=dict(color="white", width=2)), legendgroup=name, hovertemplate="L = %{x:.2f}<br>simulated %{y:.4f}<extra>" + name + "</extra>")) fig.add_annotation(x=Ls[-1], y=np.log10(exact[name][-1]), text=name, showarrow=False, xanchor="right", yanchor="bottom", yshift=4, font=dict(color=INK, size=13)) fig.update_xaxes(title="width L of the rectangle", range=[1, widths.max() + 0.2]) ticks = [0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5] fig.update_yaxes(title="probability of missing the optimum", type="log", range=[np.log10(0.004), np.log10(0.6)], tickvals=ticks, ticktext=[f"{t:g}" for t in ticks]) style(fig, height=430, showlegend=True) fig.update_layout(legend=dict(x=0.02, y=0.02, xanchor="left", yanchor="bottom")) fig.show() ``` The simulations agree with the exact values to within sampling error. k-means++ helps exactly where the problem is worst: the wider the rectangle, the bigger the disaster at the bad fixed point, and the less likely k-means++ is to land there. ## Restarts in scikit-learn {#sec-sklearn} `sklearn.cluster.KMeans` uses both remedies at once. With `init="random"` it seeds from random data points, and with `n_init` it keeps the best of that many runs, so it fails only if *every* run fails, with probability $(\tfrac13)^{\texttt{n\_init}}$ here. ```{python} from sklearn.cluster import KMeans seeds = range(300) for init, n_init in [("random", 1), ("random", 3), ("random", 10), ("k-means++", 1)]: costs = [KMeans(n_clusters=2, init=init, n_init=n_init, random_state=s).fit(X).inertia_ for s in seeds] rate = np.mean(np.array(costs) > 1 + 1e-9) print(f'init="{init}", n_init={n_init:2d}: missed the optimum in {rate:.1%} of {len(seeds)} seeds') ``` The rates for `init="random"` fall roughly like $\tfrac13$, $\tfrac1{27}$, and $3^{-10}$. scikit-learn's `"k-means++"` fails even less often than the $1/(2L^2+2)=5\%$ computed above, because it is a *greedy* variant: for each new center it draws a few candidates by the $D^2$ rule (two, when $k=2$) and keeps the one that lowers the cost the most. It picks the bad vertical neighbor only if every candidate is that neighbor. Even so, no seeding rule makes Lloyd's algorithm exact. The bad fixed point is still there, and its basin (@fig-basins) still contains two of the six pairs of data points; good seeding just makes it unlikely that we start inside it.