B-plot writes when the script runsDictionary of Applied Machine Learning · clustering
Numerical companion to the entry clustering: it recomputes what the entry states and prints one line per check
Numerical companion to the glossary entry 'clustering'. A photograph taken in the Wachau (clustering_wachau_photo.jpg, the author's own, 640 x 480 px) is cut into 96 x 128 patches of 5 x 5 pixels. The 12288 patches are the data points, and the feature vector of a patch holds its mean red, green and blue values. The clusters are surface types (sky, cloud, shaded forest, sunlit slopes, rock), and painting each patch by its cluster segments the photograph.
Run it without installing anything:uv run https://dictionaryofml.org/terms/clustering.py
uv downloads this script and the pinned NumPy and Matplotlib it needs, then runs it; the script fetches any input file it uses. To keep the output files, download clustering.py into a folder and run uv run clustering.py there. With NumPy and Matplotlib already installed, python3 clustering.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download clustering.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""Five clustering methods segment a photograph of the Wachau by patch
color: k-means (also on a coreset), a GMM, DBSCAN, spectral clustering on
the nearest-neighbor graph of the colors, and convex clustering.
Purpose
-------
Numerical companion to the glossary entry 'clustering'. A photograph
taken in the Wachau (clustering_wachau_photo.jpg, the author's own,
640 x 480 px) is cut into 96 x 128 patches of 5 x 5 pixels. The 12288
patches are the data points, and the feature vector of a patch holds its
mean red, green and blue values. The clusters are surface types (sky,
cloud, shaded forest, sunlit slopes, rock), and painting each patch by
its cluster segments the photograph.
The demo checks the entry's claims, paragraph by paragraph: k-means
with k = 5 assigns every patch to the nearest of five cluster centroids
(five representative colors) and painting compresses the photograph to
five colors; k-means on a coreset, a small weighted subset of the
patches, gives nearly the same clustering error as k-means on all
patches; the degrees of belonging of a GMM sum to one and about one in
ten patches has none above 0.8; DBSCAN with eps = 0.03 and m* = 32 joins
sky and cloud through a chain of intermediate colors, so its two largest
clusters hold most of the patches, and the rare mixed colors along the
horizon are noise; on the 10-nearest-neighbor graph of the colors,
spectral clustering yields clusters that are connected in the graph; the
number of clusters of convex clustering falls along the path of its
penalty parameter, and one value on the path leaves about five clusters
of 80 or more patches.
Deterministic: k-means and the GMM are initialized by splitting the
patches into equal parts ordered by brightness, the k-means step of
spectral clustering by farthest-first seeding from the first patch, and
the coreset is every tenth patch in brightness order (no randomness).
Kept light on purpose: the distances, the neighbor graph, the new
feature vectors of spectral clustering and the convex-clustering solves
all work on blocks of rows or on the edge list, so no matrix with one row
and one column per patch is ever formed.
Self-contained: numpy + matplotlib only.
Blocks
------
[B-patches] Cut the photograph into 96 x 128 patches with mean-RGB
feature vectors; check the counts.
[B-families] The three design choices need no computation; this block
only records the number of clusters handed to the methods.
[B-kmeans] k-means with k = 5: check the nearest-centroid assignment,
the five-color compression, and that k-means on a coreset
of 1229 weighted patches nearly matches the full run.
[B-gmm] A GMM with k = 5, hard-assigned by the largest degree of
belonging: check that the degrees of belonging sum to one
and that 5-20% of the patches have none above 0.8.
[B-dbscan] DBSCAN in color space with eps = 0.03 and m* = 32: check
that every cluster is a chain of core points at most eps
apart, that noise patches are not within eps of a core
point, and that the two largest clusters hold most patches.
[B-graph] The 10-nearest-neighbor graph of the colors: spectral
clustering with k = 5 (eigenvectors of the Laplacian
matrix, then k-means); check that the clusters are
connected in the graph.
[B-convex] Convex clustering with the edge weights of that graph along
a path of penalty parameters: check that the number of
clusters falls along the path and that one value leaves
about five clusters of 80 or more patches.
[B-plot] One painted copy of the photograph per method (each patch in
the mean color of its cluster, white lines between clusters,
DBSCAN noise in black), plus the preview.
Outputs
-------
clustering_wachau_photo.jpg : input -- the photograph (640 x 480 px)
clustering_photo_kmeans.png, _gmm.png, _dbscan.png, _spectral.png,
_convex.png : the patches painted by cluster, for the entry's figure
clustering.png : preview (checking only) -- the photograph and its five
segmentations
"""
from pathlib import Path
import urllib.request
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.image import imread
OUT_DIR = Path(__file__).parent
ASSET_URL = "https://dictionaryofml.org/demoassets/"
def demo_asset(name):
"""Path to the input file `name`, fetched once if it is not already here.
Inside the repository the committed copy is used (pythondemos/ first,
then assets/), so a run there reproduces the committed outputs exactly.
Run from anywhere else -- from the copy published on dictionaryofml.org,
say -- the file is downloaded and cached beside the script, so the demo
is runnable without the repository.
"""
for cand in (OUT_DIR / name, OUT_DIR.parent / "assets" / name):
if cand.exists():
return cand
cached = OUT_DIR / name
print(f" fetching {name} from {ASSET_URL}{name}")
with urllib.request.urlopen(ASSET_URL + name, timeout=180) as resp:
cached.write_bytes(resp.read())
return cached
report = [] # collects (check name, pass/fail) pairs
def check(name, ok): # records and prints one verification
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
photograph cut into 96 x 128 = 12288 patches of 5 px; feature vector = mean (R, G, B) in [0, 1]
Cut the photograph into 96 x 128 patches with mean-RGB feature vectors; check the counts.
PHOTO_PATCH = 5 # side length of a patch in pixels
PER_PATCH = (10 // PHOTO_PATCH) ** 2 # 4: patches per former 10 x 10 patch
BIG = 20 * PER_PATCH # a cluster of this many patches counts as large
snap = np.asarray(imread(demo_asset("clustering_wachau_photo.jpg")),
dtype=float) / 255.0
p_rows, p_cols = snap.shape[0] // PHOTO_PATCH, snap.shape[1] // PHOTO_PATCH
colors = snap[:p_rows * PHOTO_PATCH, :p_cols * PHOTO_PATCH].reshape(
p_rows, PHOTO_PATCH, p_cols, PHOTO_PATCH, 3).mean(axis=(1, 3))
C = colors.reshape(-1, 3) # feature vectors: mean RGB
print(f" photograph cut into {p_rows} x {p_cols} = {len(C)} patches of "
f"{PHOTO_PATCH} px; feature vector = mean (R, G, B) in [0, 1]")
check("[B-patches] 96 x 128 = 12288 patches, each with a three-number feature vector",
C.shape == (12288, 3) and (p_rows, p_cols) == (96, 128))
[ok] [B-patches] 96 x 128 = 12288 patches, each with a three-number feature vector k = 5 clusters for k-means, the GMM and spectral clustering; DBSCAN and convex clustering determine their number of clusters themselves k-means with k = 5: cluster sizes [2453, 3201, 1996, 2149, 2489], centroid brightness [0.15, 0.32, 0.47, 0.62, 0.79], clustering error 0.0165
The three design choices need no computation; this block only records the number of clusters handed to the methods.
K = 5
print(f" k = {K} clusters for k-means, the GMM and spectral clustering; DBSCAN "
f"and convex clustering determine their number of clusters themselves")
k-means with k = 5: check the nearest-centroid assignment, the five-color compression, and that k-means on a coreset of 1229 weighted patches nearly matches the full run.
def assign(A, cents):
return ((A[:, None, :] - cents) ** 2).sum(axis=2).argmin(axis=1)
def lloyd(A, cents, weights=None):
"""Lloyd's algorithm: alternate nearest-centroid assignment and
(weighted) averaging until the assignment no longer changes."""
w = np.ones(len(A)) if weights is None else weights
lab = assign(A, cents)
for _ in range(300):
new = np.stack([(w[lab == c, None] * A[lab == c]).sum(axis=0) / w[lab == c].sum()
if (lab == c).any() else cents[c] for c in range(len(cents))])
new_lab = assign(A, new)
if np.array_equal(new_lab, lab) and np.allclose(new, cents):
return new_lab, new
lab, cents = new_lab, new
return lab, cents
def brightness_start(A, k):
"""Deterministic start: k equal parts of the points ordered by brightness."""
parts = np.array_split(np.argsort(A.sum(axis=1)), k)
return np.stack([A[g].mean(axis=0) for g in parts])
def clustering_error(A, cents):
return float(((A - cents[assign(A, cents)]) ** 2).sum(axis=1).mean())
lab_km, cents_km = lloyd(C, brightness_start(C, K))
bright = cents_km.mean(axis=1)
err_km = clustering_error(C, cents_km)
print(f" k-means with k = {K}: cluster sizes {np.bincount(lab_km).tolist()}, "
f"centroid brightness {np.round(np.sort(bright), 2).tolist()}, "
f"clustering error {err_km:.4f}")
check("[B-kmeans] every patch is assigned to the cluster centroid nearest in "
"Euclidean distance", np.array_equal(lab_km, assign(C, cents_km)))
check("[B-kmeans] the five cluster centroids span bright (sky) to dark (forest)",
bright.max() - bright.min() > 0.4)
check("[B-kmeans] painting each patch in the color of its centroid leaves five colors",
len(np.unique(cents_km[lab_km], axis=0)) == K)
# a coreset: every tenth patch in brightness order, each weighted by 10
order = np.argsort(C.sum(axis=1))
core_idx = order[::10]
S, w_S = C[core_idx], np.full(len(core_idx), 10.0)
lab_S, cents_S = lloyd(S, brightness_start(S, K), w_S)
err_S = clustering_error(C, cents_S)
worst = 0.0
for trial in range(20): # random centroid choices
cents = C[(trial * 613 + np.arange(K) * 997) % len(C)]
full = clustering_error(C, cents)
sub = float((w_S * ((S - cents[assign(S, cents)]) ** 2).sum(axis=1)).sum()
/ w_S.sum())
worst = max(worst, abs(sub - full) / full)
print(f" coreset of {len(S)} weighted patches: weighted clustering error within "
f"{100 * worst:.1f}% of the full one over 20 centroid choices; k-means on "
f"the coreset gives clustering error {err_S:.4f} on all patches")
check("[B-kmeans] the coreset's weighted clustering error is within 15% of the "
"full one for every tried choice of the centroids", worst < 0.15)
check("[B-kmeans] k-means on the coreset is within 5% of k-means on all patches",
err_S < 1.05 * err_km)
[ok] [B-kmeans] every patch is assigned to the cluster centroid nearest in Euclidean distance [ok] [B-kmeans] the five cluster centroids span bright (sky) to dark (forest) [ok] [B-kmeans] painting each patch in the color of its centroid leaves five colors coreset of 1229 weighted patches: weighted clustering error within 2.9% of the full one over 20 centroid choices; k-means on the coreset gives clustering error 0.0166 on all patches [ok] [B-kmeans] the coreset's weighted clustering error is within 15% of the full one for every tried choice of the centroids [ok] [B-kmeans] k-means on the coreset is within 5% of k-means on all patches GMM with k = 5: cluster sizes [2408, 3968, 1460, 1942, 2510]; 12% of the patches have no degree of belonging above 0.8
A GMM with k = 5, hard-assigned by the largest degree of belonging: check that the degrees of belonging sum to one and that 5-20% of the patches have none above 0.8.
def normal_pdf(A, mean, cov):
d = A.shape[1]
diff = A - mean
quad = np.einsum("ij,jk,ik->i", diff, np.linalg.inv(cov), diff)
return np.exp(-0.5 * quad) / np.sqrt(((2.0 * np.pi) ** d)
* np.linalg.det(cov))
def em_gmm(A, k, iters=300):
"""EM for a GMM with k components; returns hard labels and posteriors."""
d = A.shape[1]
parts = np.array_split(np.argsort(A.sum(axis=1)), k)
mu = np.stack([A[g].mean(axis=0) for g in parts])
sig = np.stack([np.cov(A[g].T) + 1e-5 * np.eye(d) for g in parts])
pk = np.full(k, 1.0 / k)
for _ in range(iters):
joint = np.stack([pk[c] * normal_pdf(A, mu[c], sig[c])
for c in range(k)])
post = joint / np.maximum(joint.sum(axis=0), 1e-300)
wgt = post.sum(axis=1)
pk = wgt / len(A)
mu = (post @ A) / wgt[:, None]
sig = np.stack([(post[c][:, None] * (A - mu[c])).T @ (A - mu[c])
/ wgt[c] + 1e-5 * np.eye(d) for c in range(k)])
return post.argmax(axis=0), post
lab_gmm, post = em_gmm(C, K)
unsure = float((post.max(axis=0) < 0.8).mean())
print(f" GMM with k = {K}: cluster sizes {np.bincount(lab_gmm).tolist()}; "
f"{unsure:.0%} of the patches have no degree of belonging above 0.8")
check("[B-gmm] the degrees of belonging of each patch sum to one",
np.allclose(post.sum(axis=0), 1.0))
check("[B-gmm] between 5% and 20% of the patches have no degree of belonging above 0.8",
0.05 <= unsure <= 0.20)
[ok] [B-gmm] the degrees of belonging of each patch sum to one [ok] [B-gmm] between 5% and 20% of the patches have no degree of belonging above 0.8 DBSCAN with eps = 0.03, m* = 32: 6 clusters of sizes [6445, 4360, 160, 116, 41, 26], 9% noise
DBSCAN in color space with eps = 0.03 and m* = 32: check that every cluster is a chain of core points at most eps apart, that noise patches are not within eps of a core point, and that the two largest clusters hold most patches.
def dbscan(within, core):
"""Clusters as chains of core points at most eps apart; -1 marks noise."""
lab = -np.ones(len(core), dtype=int)
current = 0
for seed in range(len(core)):
if lab[seed] != -1 or not core[seed]:
continue
lab[seed] = current
stack = [seed]
while stack:
u = stack.pop()
if not core[u]:
continue
for v in np.flatnonzero(within[u]):
if lab[v] == -1:
lab[v] = current
stack.append(v)
current += 1
return lab
def components(mask, adjacency):
"""Number of connected components of the subgraph induced by mask."""
seen = np.zeros(len(mask), dtype=bool)
count = 0
for s in np.flatnonzero(mask):
if seen[s]:
continue
count += 1
stack = [s]
seen[s] = True
while stack:
u = stack.pop()
for v in np.flatnonzero(adjacency[u] & mask & ~seen):
seen[v] = True
stack.append(v)
return count
EPS, MIN_PTS = 0.03, 8 * PER_PATCH # m* = 32: 8 per 10 x 10 patch, scaled
NEIGHBORS = 10 # for the graph built below
BLOCK = 1024 # patches per block of distances
def neighborhoods(A, eps, neighbors, block=BLOCK):
"""The eps-neighborhood of every patch and its nearest neighbors, one
block of rows at a time: the distances of a block are formed, used and
dropped, so the full distance matrix is never held."""
m = len(A)
sq = (A ** 2).sum(axis=1) # |a - b|^2 = |a|^2 + |b|^2 - 2 a.b
near = np.zeros((m, m), dtype=bool) # which patches lie within eps
idx = np.empty((m, neighbors), dtype=np.int64)
knn_d = np.empty((m, neighbors))
for start in range(0, m, block):
stop = min(start + block, m)
d = np.sqrt(np.maximum(sq[start:stop, None] + sq[None, :]
- 2.0 * (A[start:stop] @ A.T), 0.0))
near[start:stop] = d <= eps
d[np.arange(stop - start), np.arange(start, stop)] = np.inf # not itself
# a full row sort, not a partial selection: patches of uniform sky
# have exactly equal distances, and the two rules then keep different
# neighbors, which moves a few patches between spectral clusters
take = np.argsort(d, axis=1)[:, :neighbors]
idx[start:stop] = take
knn_d[start:stop] = np.take_along_axis(d, take, axis=1)
return near, idx, knn_d
within, knn_idx, knn_d = neighborhoods(C, EPS, NEIGHBORS)
core = within.sum(axis=1) >= MIN_PTS # counting the point itself
lab_db = dbscan(within, core)
n_db = lab_db.max() + 1
sizes_db = np.bincount(lab_db[lab_db >= 0])
noise = lab_db < 0
print(f" DBSCAN with eps = {EPS}, m* = {MIN_PTS}: {n_db} clusters of sizes "
f"{sorted(sizes_db.tolist(), reverse=True)}, {noise.mean():.0%} noise")
check("[B-dbscan] every cluster is connected through a chain of core points at "
"most eps apart",
all(components((lab_db == c) & core, within) == 1 for c in range(n_db)))
check("[B-dbscan] no noise patch is a core point or within eps of one",
not core[noise].any() and not within[noise][:, core].any())
check("[B-dbscan] sky with cloud and the land: the two largest clusters hold over "
"80% of the patches, 1-30% are noise",
np.sort(sizes_db)[-2:].sum() > 0.8 * len(C) and 0.01 <= noise.mean() <= 0.30)
[ok] [B-dbscan] every cluster is connected through a chain of core points at most eps apart [ok] [B-dbscan] no noise patch is a core point or within eps of one [ok] [B-dbscan] sky with cloud and the land: the two largest clusters hold over 80% of the patches, 1-30% are noise spectral clustering on the 10-nearest-neighbor color graph (75830 edges, 1 connected component(s)): cluster sizes [1831, 2671, 4826, 2, 2958], 5 of 5 clusters connected in the graph
The 10-nearest-neighbor graph of the colors: spectral clustering with k = 5 (eigenvectors of the Laplacian matrix, then k-means); check that the clusters are connected in the graph.
def knn_graph(A, idx, knn_dist):
"""Symmetrized nearest-neighbor graph with Gaussian edge weights, held as
an edge list (plus the adjacency indicator the connectivity checks read)."""
m = len(A)
rows = np.repeat(np.arange(m), idx.shape[1])
pairs = np.unique(np.stack([np.minimum(rows, idx.ravel()),
np.maximum(rows, idx.ravel())], axis=1), axis=0)
ei, ej = pairs[:, 0], pairs[:, 1]
bw = float(np.median(knn_dist))
w = np.exp(-((A[ei] - A[ej]) ** 2).sum(axis=1) / (2.0 * bw ** 2))
adjacency = np.zeros((m, m), dtype=bool)
adjacency[ei, ej] = True
adjacency |= adjacency.T
return adjacency, ei, ej, w
def farthest_first(A, k):
"""Deterministic start: the first point, then always the farthest one."""
chosen = [0]
for _ in range(k - 1):
d2 = ((A[:, None, :] - A[chosen]) ** 2).sum(axis=2).min(axis=1)
chosen.append(int(d2.argmax()))
return A[chosen]
adj, edge_i, edge_j, edge_w = knn_graph(C, knn_idx, knn_d)
def spectral_embedding(ei, ej, w, m, k, steps=600):
"""The k eigenvectors of the normalized Laplacian with the smallest
eigenvalues. They are the k largest of the normalized adjacency operator,
found by the Lanczos method with full reorthogonalization: every step is
one pass over the edge list, so the normalized Laplacian is never written
out as a matrix. Six hundred steps bring the residual of every returned
eigenvector below 1e-4; at three hundred it is 1e-3, which is not enough
to separate these five eigenvalues, all of which lie within 1e-4 of one."""
deg = (np.bincount(ei, weights=w, minlength=m)
+ np.bincount(ej, weights=w, minlength=m))
scale = 1.0 / np.sqrt(deg)
def normalized_adjacency(x):
y = scale * x
return scale * (np.bincount(ei, weights=w * y[ej], minlength=m)
+ np.bincount(ej, weights=w * y[ei], minlength=m))
gen = np.random.default_rng(0)
Q = np.zeros((m, steps))
Q[:, 0] = gen.standard_normal(m)
Q[:, 0] /= np.linalg.norm(Q[:, 0])
diag, offdiag, used = np.zeros(steps), np.zeros(steps), steps
for i in range(steps):
v = normalized_adjacency(Q[:, i])
diag[i] = Q[:, i] @ v
v -= Q[:, :i + 1] @ (Q[:, :i + 1].T @ v) # full reorthogonalization
v -= Q[:, :i + 1] @ (Q[:, :i + 1].T @ v) # twice is enough
nrm = np.linalg.norm(v)
if i + 1 == steps or nrm < 1e-10:
used = i + 1
break
offdiag[i] = nrm
Q[:, i + 1] = v / nrm
tri = (np.diag(diag[:used]) + np.diag(offdiag[:used - 1], 1)
+ np.diag(offdiag[:used - 1], -1))
vals, rot = np.linalg.eigh(tri)
return Q[:, :used] @ rot[:, np.argsort(vals)[::-1][:k]]
emb = spectral_embedding(edge_i, edge_j, edge_w, len(C), K) # k new feature entries per node
emb /= np.linalg.norm(emb, axis=1, keepdims=True)
lab_sp, _ = lloyd(emb, farthest_first(emb, K))
n_comp = components(np.ones(len(C), bool), adj)
connected_sp = sum(components(lab_sp == c, adj) == 1 for c in range(K))
print(f" spectral clustering on the {NEIGHBORS}-nearest-neighbor color graph "
f"({int(adj.sum() // 2)} edges, {n_comp} connected component(s)): cluster "
f"sizes {np.bincount(lab_sp).tolist()}, {connected_sp} of {K} clusters "
f"connected in the graph")
check("[B-graph] the new feature vectors have k entries, one row per node",
emb.shape == (len(C), K))
check("[B-graph] at least four of the five spectral clusters are connected in the graph",
connected_sp >= 4)
[ok] [B-graph] the new feature vectors have k entries, one row per node [ok] [B-graph] at least four of the five spectral clusters are connected in the graph convex clustering path (alpha: clusters / clusters of 80+ patches): 0.05: 4551/0, 0.0901: 3062/3, 0.162: 2045/16, 0.292: 1361/34, 0.527: 863/43, 0.949: 557/41, 1.71: 317/35, 3.08: 191/16, 5.55: 124/10, 10: 94/7, 18: 65/4, 32.5: 49/3, 58.5: 38/3, 105: 32/3, 190: 24/2, 342: 19/2, 616: 15/2, 1.11e+03: 12/2, 2e+03: 11/2 chosen alpha = 18: 65 clusters of sizes [6644, 4471, 602, 397, 20, 17, 15, 13, 10, 7, 7, 6, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
Convex clustering with the edge weights of that graph along a path of penalty parameters: check that the number of clusters falls along the path and that one value leaves about five clusters of 80 or more patches.
def label_components(m, ei, ej):
"""Connected-component label of every node, from a list of edges, by
union-find: no node-by-node adjacency matrix is formed."""
parent = list(range(m))
def find(x):
while parent[x] != x:
parent[x] = parent[parent[x]]
x = parent[x]
return x
for a, b in zip(ei.tolist(), ej.tolist()):
ra, rb = find(a), find(b)
if ra != rb:
parent[ra] = rb
roots = np.array([find(x) for x in range(m)])
_, first, inverse = np.unique(roots, return_index=True, return_inverse=True)
order = np.empty(len(first), dtype=int) # label by first appearance
order[np.argsort(first)] = np.arange(len(first))
return order[inverse]
def convex_clustering(A, i_idx, j_idx, a_edge, alphas, rho=1.0, iters=200,
cg_iters=25, cg_tol=1e-6):
"""Minimize sum_r ||x_r - w_r||^2 + alpha sum_edges W_rr' ||w_r - w_r'||
by alternating updates with one splitting variable per edge (Chi and
Lange, 2015), warm-started along the path of alphas. Returns, per
alpha, the labels of the clusters of patches that share a centroid.
The centroid update solves (2 I + rho E^T E) w = b by conjugate
gradients on the edge list, so the matrix is never formed or inverted."""
m, d = A.shape
bins_i = (i_idx[:, None] * d + np.arange(d)).ravel() # (node, channel) bins
bins_j = (j_idx[:, None] * d + np.arange(d)).ravel()
def edge_gather(X):
"""E^T (X_i - X_j) for a value per edge, as a value per node. One
bincount over a flattened (edge, channel) index does all channels at
once, which halves the cost of the inner loop."""
flat = X.ravel()
return (np.bincount(bins_i, weights=flat, minlength=m * d).reshape(m, d)
- np.bincount(bins_j, weights=flat, minlength=m * d).reshape(m, d))
degree = (np.bincount(i_idx, minlength=m)
+ np.bincount(j_idx, minlength=m)).astype(float)
precond = 1.0 / (2.0 + rho * degree)[:, None] # the system's diagonal
def system(X):
"""(2 I + rho E^T E) X."""
return 2.0 * X + rho * edge_gather(X[i_idx] - X[j_idx])
def solve(B, X):
"""Conjugate gradients with the diagonal preconditioner, warm-started
at X (one run per color channel). The warm start is the previous
centroid update, which is close, so the relative residual falls below
the tolerance in a handful of steps."""
R = B - system(X)
Z = precond * R
Pdir = Z.copy()
rz = (R * Z).sum(axis=0)
goal = cg_tol ** 2 * np.maximum((B * B).sum(axis=0), 1e-300)
for _ in range(cg_iters):
if np.all((R * R).sum(axis=0) <= goal):
break
Q = system(Pdir)
step = rz / np.maximum((Pdir * Q).sum(axis=0), 1e-300)
X = X + step * Pdir
R = R - step * Q
Z = precond * R
rz_next = (R * Z).sum(axis=0)
Pdir = Z + (rz_next / np.maximum(rz, 1e-300)) * Pdir
rz = rz_next
return X
V = A[i_idx] - A[j_idx]
U = np.zeros_like(V)
Wc = A.copy()
out = []
for alpha in alphas:
for _ in range(iters):
Wc = solve(2.0 * A + rho * edge_gather(V - U), Wc) # centroid update
Z = Wc[i_idx] - Wc[j_idx] + U
nz = np.maximum(np.linalg.norm(Z, axis=1), 1e-12)
V = np.maximum(0.0, 1.0 - alpha * a_edge / (rho * nz))[:, None] * Z
U = U + Wc[i_idx] - Wc[j_idx] - V
keep = np.linalg.norm(V, axis=1) == 0.0 # fused edges
out.append((float(alpha), label_components(m, i_idx[keep], j_idx[keep])))
return out
ALPHAS = np.geomspace(0.05, 2000.0, 19)
path = convex_clustering(C, edge_i, edge_j, edge_w, ALPHAS)
totals = [int(lab.max() + 1) for _, lab in path]
big = [int((np.bincount(lab) >= BIG).sum()) for _, lab in path]
# the point on the path where about k clusters of BIG or more patches remain;
# the rest are small clusters of a few patches each (ties: the larger alpha)
alpha_cc, lab_cc = min(path, key=lambda t: (
abs(int((np.bincount(t[1]) >= BIG).sum()) - K), -t[0]))
sizes_cc = sorted(np.bincount(lab_cc).tolist(), reverse=True)
print(f" convex clustering path (alpha: clusters / clusters of {BIG}+ patches): "
+ ", ".join(f"{a:.3g}: {t}/{b}" for (a, _), t, b in zip(path, totals, big)))
print(f" chosen alpha = {alpha_cc:.3g}: {len(sizes_cc)} clusters of sizes {sizes_cc}")
check("[B-convex] the number of clusters falls along the path from thousands "
"to at most twenty",
totals[0] > 2000 and totals[-1] <= 20 and all(
a >= b for a, b in zip(totals, totals[1:])))
check(f"[B-convex] a penalty parameter on the path leaves about k clusters of {BIG} "
"or more patches",
abs(sum(sz >= BIG for sz in sizes_cc) - K) <= 1)
[ok] [B-convex] the number of clusters falls along the path from thousands to at most twenty [ok] [B-convex] a penalty parameter on the path leaves about k clusters of 80 or more patches
One painted copy of the photograph per method (each patch in the mean color of its cluster, white lines between clusters, DBSCAN noise in black), plus the preview.
def paint_photo(lab, name, zoom=4):
"""Each patch in the mean color of its cluster, white lines between
clusters, DBSCAN noise in black."""
lab2 = lab.reshape(p_rows, p_cols)
img = np.zeros((p_rows, p_cols, 3))
for c in np.unique(lab):
img[lab2 == c] = C[lab == c].mean(axis=0) if c >= 0 else 0.0
big = img.repeat(zoom, axis=0).repeat(zoom, axis=1)
above = np.zeros((p_rows, p_cols), dtype=bool)
above[1:] = lab2[1:] != lab2[:-1]
left = np.zeros((p_rows, p_cols), dtype=bool)
left[:, 1:] = lab2[:, 1:] != lab2[:, :-1]
for r, c in zip(*np.nonzero(above)):
big[r * zoom, c * zoom:(c + 1) * zoom] = 1.0
for r, c in zip(*np.nonzero(left)):
big[r * zoom:(r + 1) * zoom, c * zoom] = 1.0
plt.imsave(OUT_DIR / name, np.clip(big, 0.0, 1.0))
return big
painted = {m: paint_photo(lab, f"clustering_photo_{m}.png")
for m, lab in (("kmeans", lab_km), ("gmm", lab_gmm),
("dbscan", lab_db), ("spectral", lab_sp),
("convex", lab_cc))}
fig, axes = plt.subplots(1, 6, figsize=(23, 3.6))
axes[0].imshow(snap)
axes[0].set_title("photograph (Wachau)")
for ax, m, title in zip(axes[1:], ("kmeans", "gmm", "dbscan", "spectral", "convex"),
("k-means, k = 5", "GMM, k = 5",
f"DBSCAN, eps = {EPS}, m* = {MIN_PTS} (noise black)",
"spectral clustering, k = 5",
f"convex clustering, alpha = {alpha_cc:.3g}")):
ax.imshow(np.clip(painted[m], 0, 1))
ax.set_title(title)
for ax in axes:
ax.set_xlabel("pixel column")
ax.set_ylabel("pixel row")
fig.tight_layout()
fig.savefig(OUT_DIR / "clustering.png", dpi=150)
plt.close(fig)
check("[B-plot] one painted copy of the photograph per method",
all((OUT_DIR / f"clustering_photo_{m}.png").exists() for m in painted))
passed = sum(1 for _, ok in report if ok)
print(f"\n{passed}/{len(report)} checks pass")
[ok] [B-plot] one painted copy of the photograph per method 16/16 checks pass
B-plot writes when the script runs