B-fixedpoint writes when the script runsDictionary of Applied Machine Learning · $k$-means
Numerical companion to the entry $k$-means: it recomputes what the entry states and prints one line per check
Numerical companion to the glossary entry 'kmeans' ($k$-means). The GeoSphere Austria weather station Krems (station id 3805) records the minimum and maximum air temperature of each day; this script downloads the records for 2024 from the GeoSphere data hub (dataset klima-v2-1d) and writes them to kmeans_weather.csv. Each day is a data point with feature vector (tmin, tmax), and $k$-means with $k = 2$ partitions the 366 days into a cold-season and a warm-season cluster.
Run it with python3 kmeans.py, from any directory — it writes its output files into the current directory. Requires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Download kmeans.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""Lloyd's algorithm on a year of days at Krems: the clustering error
never increases, and the iteration stops at a fixed point.
Purpose
-------
Numerical companion to the glossary entry 'kmeans' ($k$-means). The
GeoSphere Austria weather station Krems (station id 3805) records the
minimum and maximum air temperature of each day; this script downloads
the records for 2024 from the GeoSphere data hub (dataset klima-v2-1d)
and writes them to kmeans_weather.csv. Each day is a data point with
feature vector (tmin, tmax), and $k$-means with $k = 2$ partitions the
366 days into a cold-season and a warm-season cluster.
The demo checks the entry's central claims about Lloyd's algorithm:
each iteration -- assign every data point to its nearest cluster
centroid, then recompute each centroid as the mean of its assigned
data points -- never increases the clustering error; after finitely
many iterations nothing changes anymore, i.e., the iteration reaches a
fixed point, where each centroid coincides with the mean of the data
points assigned to it.
Deterministic: the initial centroids are the coldest and the warmest
day of the year (no randomness). Self-contained: numpy + matplotlib
only (stdlib urllib for the download).
Blocks
------
[B-fetch] Download the 366 daily temperature pairs at Krems for 2024
and write them to kmeans_weather.csv; check the count and
one pinned value against the archive.
[B-lloyd] Run Lloyd's algorithm with k = 2 from the coldest/warmest
day; check that the clustering error never increases and
that assignments stop changing after finitely many
iterations.
[B-fixedpoint] Verify the fixed-point property of the result: each
final cluster centroid equals the mean of the data points
assigned to it, so one more iteration changes nothing.
[B-image] Image compression and image segmentation on a subsampled
photo of the Oetscher massif (assets/oetscher.jpg): k-means
on the pixel colors with k = 4 replaces each pixel's color by
the nearest palette color (compression factor ~12), and with
k = 2 partitions the pixels into a sky-and-mountain region
and a vegetation region.
Outputs
-------
kmeans_weather.csv : date, tmin, tmax for the 366 days of 2024
kmeans_cluster1.csv : tmin, tmax of the days in the cold cluster
kmeans_cluster2.csv : tmin, tmax of the days in the warm cluster
kmeans_centroids.csv : tmin, tmax of the two final cluster centroids
kmeans_error.csv : iter, error -- clustering error per iteration
kmeans_oetscher_original.png : the subsampled photo
kmeans_oetscher_compressed.png : the photo quantized to 4 palette colors
kmeans_oetscher_mask.png : sky-and-mountain/vegetation mask (k = 2)
kmeans.png : preview (checking only) -- the two clusters
with their centroids, and the monotone
clustering error
"""
import json
import urllib.request
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
OUT_DIR = Path(__file__).parent
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}")
Download the 366 daily temperature pairs at Krems for 2024 and write them to kmeans_weather.csv; check the count and one pinned value against the archive.
URL = ("https://dataset.api.hub.geosphere.at/v1/station/historical/"
"klima-v2-1d?parameters=tlmin,tlmax&station_ids=3805"
"&start=2024-01-01&end=2024-12-31")
with urllib.request.urlopen(URL, timeout=120) as resp:
payload = json.load(resp)
params = payload["features"][0]["properties"]["parameters"]
stamps = [t[:10] for t in payload["timestamps"]]
X = np.stack([np.array(params["tlmin"]["data"], dtype=float),
np.array(params["tlmax"]["data"], dtype=float)], 1)
with open(OUT_DIR / "kmeans_weather.csv", "w") as f:
f.write("date,tmin,tmax\n")
for day, (lo, hi) in zip(stamps, X):
f.write(f"{day},{lo},{hi}\n")
check("[B-fetch] 366 daily temperature pairs downloaded for 2024",
len(X) == 366)
check("[B-fetch] the record matches the archive (Feb 1: -3.8 to 10.4)",
stamps[31] == "2024-02-01" and np.allclose(X[31], [-3.8, 10.4]))
[ok] [B-fetch] 366 daily temperature pairs downloaded for 2024 [ok] [B-fetch] the record matches the archive (Feb 1: -3.8 to 10.4) fixed point reached after 5 iterations, clustering error 15380
Run Lloyd's algorithm with k = 2 from the coldest/warmest day; check that the clustering error never increases and that assignments stop changing after finitely many iterations.
def assign(X, centroids):
"""Index of the nearest cluster centroid for every data point."""
d = ((X[:, None, :] - centroids[None, :, :]) ** 2).sum(axis=2)
return d.argmin(axis=1)
def clustering_error(X, centroids, labels):
return float(((X - centroids[labels]) ** 2).sum())
mean_temp = X.mean(axis=1)
centroids = np.stack([X[mean_temp.argmin()], X[mean_temp.argmax()]])
labels = assign(X, centroids)
errors = [clustering_error(X, centroids, labels)]
iterations = 0
while True:
iterations += 1
centroids = np.stack([X[labels == c].mean(axis=0) for c in (0, 1)])
new_labels = assign(X, centroids)
errors.append(clustering_error(X, centroids, new_labels))
if np.array_equal(new_labels, labels):
break
labels = new_labels
errors = np.array(errors)
print(f" fixed point reached after {iterations} iterations, "
f"clustering error {errors[-1]:.0f}")
check("[B-lloyd] the clustering error never increases",
bool(np.all(np.diff(errors) <= 1e-9)))
check("[B-lloyd] assignments stop changing after finitely many iterations",
iterations < 50)
check("[B-lloyd] both clusters are nonempty",
0 < int(labels.sum()) < len(X))
[ok] [B-lloyd] the clustering error never increases [ok] [B-lloyd] assignments stop changing after finitely many iterations [ok] [B-lloyd] both clusters are nonempty
Verify the fixed-point property of the result: each final cluster centroid equals the mean of the data points assigned to it, so one more iteration changes nothing.
recomputed = np.stack([X[labels == c].mean(axis=0) for c in (0, 1)])
check("[B-fixedpoint] each centroid equals the mean of its assigned "
"data points", np.allclose(recomputed, centroids))
check("[B-fixedpoint] one more iteration changes no assignment",
np.array_equal(assign(X, recomputed), labels))
header = "tmin,tmax"
np.savetxt(OUT_DIR / "kmeans_cluster1.csv", X[labels == 0],
delimiter=",", header=header, comments="", fmt="%.1f")
np.savetxt(OUT_DIR / "kmeans_cluster2.csv", X[labels == 1],
delimiter=",", header=header, comments="", fmt="%.1f")
np.savetxt(OUT_DIR / "kmeans_centroids.csv", centroids,
delimiter=",", header=header, comments="", fmt="%.2f")
np.savetxt(OUT_DIR / "kmeans_error.csv",
np.stack([np.arange(len(errors)), errors], 1),
delimiter=",", header="iter,error", comments="", fmt="%.1f")
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))
ax = axes[0]
ax.plot(X[labels == 0, 0], X[labels == 0, 1], "o", color="tab:blue",
ms=3, label="cold-season cluster")
ax.plot(X[labels == 1, 0], X[labels == 1, 1], "^", mfc="none",
mec="tab:red", ms=4, label="warm-season cluster")
ax.plot(centroids[:, 0], centroids[:, 1], "kx", ms=10, mew=2,
label="cluster centroids")
ax.set_xlabel("minimum temperature of the day (°C)")
ax.set_ylabel("maximum temperature of the day (°C)")
ax.set_title("k-means with k = 2 on the 366 days of 2024 (Krems)")
ax.legend(frameon=False, fontsize=8)
ax = axes[1]
ax.plot(np.arange(len(errors)), errors, "k.-")
ax.set_xlabel("Lloyd iteration")
ax.set_ylabel("clustering error")
ax.set_title("monotone descent to a fixed point")
fig.tight_layout()
fig.savefig(OUT_DIR / "kmeans.png", dpi=150)
[ok] [B-fixedpoint] each centroid equals the mean of its assigned data points [ok] [B-fixedpoint] one more iteration changes no assignment
B-fixedpoint writes when the script runsImage compression and image segmentation on a subsampled photo of the Oetscher massif (assets/oetscher.jpg): k-means on the pixel colors with k = 4 replaces each pixel's color by the nearest palette color (compression factor ~12), and with k = 2 partitions the pixels into a sky-and-mountain region and a vegetation region.
# A photo of the Oetscher massif, subsampled by keeping every 36th pixel
# in each direction: the blue sky and gray summit, the dark forest and
# the bright meadow give a small image whose pixel colors k-means can
# quantize (compression) and partition (segmentation).
from matplotlib.image import imread
SUBSAMPLE = 36
photo = imread(OUT_DIR.parent / "assets" / "oetscher.jpg") / 255.0
sub = photo[::SUBSAMPLE, ::SUBSAMPLE]
colors = sub.reshape(-1, 3) # one RGB vector per pixel
luminance = colors.mean(axis=1)
check("[B-image] subsampling keeps every 36th pixel",
sub.shape == (97, 129, 3))
print(f" photo subsampled to {sub.shape[0]} x {sub.shape[1]} pixels")
def lloyd_colors(k):
"""Lloyd's algorithm on the pixel colors; deterministic init at the
colors of the pixels whose luminances sit at k evenly spaced
quantiles."""
qs = np.quantile(luminance, np.linspace(0.0, 1.0, k))
cents = np.array([colors[np.abs(luminance - q).argmin()] for q in qs])
labs = assign(colors, cents)
while True:
cents = np.stack([colors[labs == c].mean(axis=0) for c in range(k)])
new = assign(colors, cents)
if np.array_equal(new, labs):
return cents, labs
labs = new
K_PALETTE = 4
palette, plabels = lloyd_colors(K_PALETTE)
compressed = palette[plabels].reshape(sub.shape)
npix = colors.shape[0]
bits_orig = 24 * npix # 8-bit RGB per pixel
bits_comp = 2 * npix + K_PALETTE * 24 # 2-bit index + palette
factor = bits_orig / bits_comp
check("[B-image] 4-color palette compresses by a factor of about 12",
11.5 < factor < 12.0)
mask_cents, mlabels = lloyd_colors(2)
# the bluer of the two clusters collects the sky and the gray summit,
# the other the vegetation
blueness = [float((colors[mlabels == c][:, 2]
- colors[mlabels == c][:, 0]).mean()) for c in (0, 1)]
sky = int(np.argmax(blueness))
mask = (mlabels == sky).reshape(sub.shape[:2])
check("[B-image] the sky-and-mountain cluster is markedly bluer",
blueness[sky] > 0.3 > 0.1 > blueness[1 - sky])
check("[B-image] both regions are present",
0.1 < mask.mean() < 0.9)
UPSCALE = 6 # keep the pixels crisp
def save_pixels(img01, path):
from PIL import Image
rgb8 = (np.clip(img01, 0.0, 1.0) * 255.0 + 0.5).astype(np.uint8)
im = Image.fromarray(rgb8)
im = im.resize((im.width * UPSCALE, im.height * UPSCALE), Image.NEAREST)
im.save(path)
save_pixels(sub, OUT_DIR / "kmeans_oetscher_original.png")
save_pixels(compressed, OUT_DIR / "kmeans_oetscher_compressed.png")
save_pixels(np.repeat(mask[:, :, None], 3, axis=2).astype(float),
OUT_DIR / "kmeans_oetscher_mask.png")
print(f" wrote kmeans_oetscher_original/compressed/mask.png "
f"(compression factor {factor:.1f})")
failed = [name for name, ok in report if not ok]
print(f"{len(report) - len(failed)}/{len(report)} checks passed"
+ (f", FAILED: {failed}" if failed else ""))
[ok] [B-image] subsampling keeps every 36th pixel photo subsampled to 97 x 129 pixels [ok] [B-image] 4-color palette compresses by a factor of about 12 [ok] [B-image] the sky-and-mountain cluster is markedly bluer [ok] [B-image] both regions are present wrote kmeans_oetscher_original/compressed/mask.png (compression factor 12.0) 11/11 checks passed