Unsupervised Learning Lab
Unsupervised learning finds structure in unlabeled data, and a self-stretching map can even solve the traveling salesman problem.
Based on Chapter 11, “Supervised and unsupervised learning,” Sections 11.2, 11.4 and 11.8, in Alaa Khamis, Optimization Algorithms: AI techniques for design, planning, and control problems, Manning Publications, 2024 (ISBN 978-1633438835).
A self-organizing ring of neurons stretching over 194 cities in Qatar.
Unsupervised learning
Unsupervised learning works with unlabeled data: clustering groups similar objects, and dimensionality reduction maps data to fewer dimensions while keeping its structure.
Section 11.2 of the book surveys the main families of clustering and dimensionality-reduction algorithms.
Similarity measures
Pick a measure to see its contours of equal distance; click to place a point.
Chapter 11 of the book gives the formula for each measure and works a haversine example.
Distance explorer
Click to place a pointOn the globe, the haversine formula gives great-circle distances. Pick two cities to compare it with a flat-map straight line.
Haversine
K-means clustering
K-means alternates two steps: assign every point to its nearest center, then move each center to the mean of its points. It converges, but to a local minimum that depends on the seeds. Click Next step, or try Bad seeds.
Chapter 11 of the book covers k-means++ seeding, the cost of k-means, and how it compares with other algorithms.
K-means step by step
Click to add pointsValidating clusters
Internal measures (silhouette, Davies-Bouldin) judge a clustering from the data alone; external ones (ARI, NMI) compare it with known classes. The chart scores each k for the data above.
Chapter 11 of the book defines each validation index and explains the elbow method for choosing k.
Choosing k for the data above
Beyond k-means: DBSCAN and agglomerative clustering
DBSCAN grows clusters of any shape from dense core points; agglomerative clustering keeps merging the closest clusters. Change ε and the linkage and compare ARI.
Three algorithms, one dataset
Self-organizing maps
In a self-organizing map, the neuron closest to each input wins, and it and its lattice neighbors move toward the input. Change σ to compare neighborhood shapes, then train the 2D map below.
Section 11.4 of the book gives Algorithm 11.1 and the update rule, and explains why a SOM preserves topology.
Neighborhood functions (equation 11.6)
A 2D map unfolding over data
Neurons drawn at their weightsSolving the TSP with a self-organizing map
Make the lattice a ring: it stretches like an elastic band over the cities, and reading them in ring order gives a tour. Press Run, then polish with 2-opt.
Section 11.8 of the book walks through Listing 11.5, its settings and results on qa194.
194 cities in Qatar
TSPLIB qa194, from the book’s repoPython in Google Colab
Open the book’s notebooks or the companion notebook with the from-scratch code below in Colab.
Notebooks
From-scratch implementations
# Similarity measures for clustering import math import numpy as np x, y = np.array([1.0, 2.0]), np.array([4.0, 6.0]) print("Euclidean:", np.sqrt(((x - y) ** 2).sum())) # 5.0 print("Manhattan:", np.abs(x - y).sum()) # 7.0 print("Chebyshev:", np.abs(x - y).max()) # 4.0 print("cosine similarity:", round(x @ y / (np.linalg.norm(x) * np.linalg.norm(y)), 4)) # Mahalanobis distance accounts for the spread and correlation of the data rng = np.random.default_rng(0) data = rng.multivariate_normal([0, 0], [[4, 3], [3, 4]], 500) S_inv = np.linalg.inv(np.cov(data.T)) maha = lambda a, b: math.sqrt((a - b) @ S_inv @ (a - b)) p, q = np.array([2.0, 2.0]), np.array([2.0, -2.0]) print(f"same Euclidean distance from the origin ({np.linalg.norm(p):.2f}), Mahalanobis: along the trend {maha(p, np.zeros(2)):.2f}, across it {maha(q, np.zeros(2)):.2f}") def haversine(lat1, lon1, lat2, lon2, R=6371.0): """Great-circle distance in km.""" la1, lo1, la2, lo2 = map(math.radians, (lat1, lon1, lat2, lon2)) a = math.sin((la2 - la1) / 2) ** 2 + math.cos(la1) * math.cos(la2) * math.sin((lo2 - lo1) / 2) ** 2 return 2 * R * math.atan2(math.sqrt(a), math.sqrt(1 - a)) print(f"Los Angeles to Madrid: {haversine(34.0522, -118.2437, 40.4168, -3.7038):,.1f} km") # about 9,363 km
# K-means from scratch, with internal (silhouette, Davies-Bouldin) and external (ARI, NMI) validation import math import numpy as np rng = np.random.default_rng(1) centers = np.array([[-2, -1], [2, -1.5], [0, 2]]) X = np.vstack([c + 0.6 * rng.standard_normal((60, 2)) for c in centers]) truth = np.repeat([0, 1, 2], 60) def kmeans(X, k, rng, plus_plus=True, iters=100): if plus_plus: # k-means++: spread the initial centers out C = [X[rng.integers(len(X))]] for _ in range(k - 1): d2 = np.min(((X[:, None] - np.array(C)) ** 2).sum(-1), axis=1) C.append(X[rng.choice(len(X), p=d2 / d2.sum())]) C = np.array(C) else: C = X[rng.choice(len(X), k, replace=False)] for it in range(iters): labels = ((X[:, None] - C) ** 2).sum(-1).argmin(1) # assign each point to its nearest center newC = np.array([X[labels == j].mean(0) if (labels == j).any() else C[j] for j in range(k)]) if np.allclose(newC, C): break C = newC # move each center to the mean of its points inertia = ((X - C[labels]) ** 2).sum() return labels, C, inertia, it + 1 def silhouette(X, labels): D = np.sqrt(((X[:, None] - X) ** 2).sum(-1)); s = [] for i in range(len(X)): same = labels == labels[i] a = D[i, same].sum() / max(same.sum() - 1, 1) b = min(D[i, labels == j].mean() for j in set(labels) if j != labels[i]) s.append((b - a) / max(a, b)) return np.mean(s) def davies_bouldin(X, labels): ks = sorted(set(labels)); C = np.array([X[labels == k].mean(0) for k in ks]) S = np.array([np.sqrt(((X[labels == k] - C[i]) ** 2).sum(1)).mean() for i, k in enumerate(ks)]) return np.mean([max((S[i] + S[j]) / np.linalg.norm(C[i] - C[j]) for j in range(len(ks)) if j != i) for i in range(len(ks))]) def contingency(a, b): A, B = sorted(set(a)), sorted(set(b)) return np.array([[np.sum((a == i) & (b == j)) for j in B] for i in A]) def ari(a, b): n = contingency(a, b); comb = lambda v: v * (v - 1) / 2 idx, ra, rb, tot = comb(n).sum(), comb(n.sum(1)).sum(), comb(n.sum(0)).sum(), comb(n.sum()) exp = ra * rb / tot return (idx - exp) / (0.5 * (ra + rb) - exp) def nmi(a, b): n = contingency(a, b) / len(a); pa, pb = n.sum(1), n.sum(0) mi = sum(n[i, j] * math.log(n[i, j] / (pa[i] * pb[j])) for i in range(n.shape[0]) for j in range(n.shape[1]) if n[i, j] > 0) h = lambda p: -sum(v * math.log(v) for v in p if v > 0) return mi / ((h(pa) + h(pb)) / 2) best_of = lambda k, runs=5: min((kmeans(X, k, rng) for _ in range(runs)), key=lambda r: r[2]) # keep the best of 5 starts print(" k inertia silhouette Davies-Bouldin") for k in range(2, 7): labels, C, inertia, its = best_of(k) print(f" {k} {inertia:8.1f} {silhouette(X, labels):.3f} {davies_bouldin(X, labels):.3f}") labels, C, inertia, its = kmeans(X, 3, rng) print(f"k = 3 converged in {its} iterations; ARI = {ari(truth, labels):.3f}, NMI = {nmi(truth, labels):.3f}")
# K-means vs. DBSCAN vs. agglomerative clustering on two interleaving moons import numpy as np rng = np.random.default_rng(0) t = rng.uniform(0, np.pi, 150) X = np.vstack([np.c_[np.cos(t), np.sin(t)], np.c_[1 - np.cos(t), 0.5 - np.sin(t)]]) + 0.06 * rng.standard_normal((300, 2)) truth = np.repeat([0, 1], 150) D = np.sqrt(((X[:, None] - X) ** 2).sum(-1)) def dbscan(D, eps, min_pts): """Label -1 is noise. Core points have at least min_pts neighbors within eps; clusters grow from core points.""" n, labels, c = len(D), np.full(len(D), -2), -1 for i in range(n): if labels[i] != -2: continue nb = np.flatnonzero(D[i] <= eps) if len(nb) < min_pts: labels[i] = -1; continue c += 1; labels[i] = c; queue = list(nb) while queue: j = queue.pop() if labels[j] == -1: labels[j] = c # border point if labels[j] != -2: continue labels[j] = c nb_j = np.flatnonzero(D[j] <= eps) if len(nb_j) >= min_pts: queue.extend(nb_j) # j is a core point: keep expanding return labels def single_linkage(D, k): """Agglomerative clustering: repeatedly merge the two closest clusters until k remain.""" labels = np.arange(len(D)) order = np.dstack(np.unravel_index(np.argsort(D, axis=None), D.shape))[0] groups = len(D) for i, j in order: if groups == k: break if i < j and labels[i] != labels[j]: labels[labels == labels[j]] = labels[i]; groups -= 1 return np.unique(labels, return_inverse=True)[1] def kmeans2(X, rng): C = X[rng.choice(len(X), 2, replace=False)] for _ in range(100): lab = ((X[:, None] - C) ** 2).sum(-1).argmin(1) C = np.array([X[lab == j].mean(0) for j in range(2)]) return lab accuracy = lambda lab: max(np.mean(lab == truth), np.mean(lab == 1 - truth)) db = dbscan(D, 0.2, 5) print(f"K-means (k = 2): agreement with the true moons {accuracy(kmeans2(X, rng)):.0%}") print(f"single linkage (k = 2): {accuracy(single_linkage(D, 2)):.0%}") print(f"DBSCAN (eps = 0.2, minPts = 5): {db.max() + 1} clusters, {np.sum(db == -1)} noise points, agreement {accuracy(np.where(db < 0, 0, db)):.0%}")
# A self-organizing map from scratch (Section 11.4, Algorithm 11.1) import numpy as np rng = np.random.default_rng(0) # data: random colors (3D inputs) mapped onto a 12 x 12 lattice of neurons X = rng.random((500, 3)) rows = cols = 12 W = rng.random((rows * cols, 3)) # randomly initialize the weights grid = np.array([(r, c) for r in range(rows) for c in range(cols)], float) def neighborhood(u, sigma): """Gaussian neighborhood of equation 11.6, on the lattice distance to the winner u.""" d2 = ((grid - grid[u]) ** 2).sum(1) return np.exp(-d2 / (2 * sigma ** 2)) steps, alpha0, sigma0 = 3000, 0.5, 6.0 for s in range(steps): x = X[rng.integers(len(X))] # pick a random input u = ((W - x) ** 2).sum(1).argmin() # best matching unit (winning neuron) frac = s / steps alpha, sigma = alpha0 * (1 - frac), max(0.5, sigma0 * (1 - frac)) # both shrink over time W += alpha * neighborhood(u, sigma)[:, None] * (x - W) # pull the winner and its neighbors toward x # Topology preservation: neighboring neurons should hold similar colors diff = np.mean([np.linalg.norm(W[i] - W[i + 1]) for i in range(len(W) - 1) if (i + 1) % cols]) rand = np.mean([np.linalg.norm(W[i] - W[j]) for i, j in rng.integers(len(W), size=(500, 2))]) print(f"mean color difference between lattice neighbors {diff:.3f} vs. random pairs {rand:.3f}") q_err = np.mean([np.sqrt(((W - x) ** 2).sum(1)).min() for x in X]) print(f"quantization error (mean distance from an input to its winner): {q_err:.3f}")
# Solving the TSP with a self-organizing map: 194 cities in Qatar (Section 11.8, Listing 11.5) import urllib.request import numpy as np URL = ("https://raw.githubusercontent.com/Optimization-Algorithms-Book/Code-Listings/" "256207c4a8badc0977286c48a6e1cfd33237a51d/Appendix%20B/data/TSP/qa194.tsp") lines = urllib.request.urlopen(URL).read().decode().splitlines() start = lines.index("NODE_COORD_SECTION") + 1 cities = np.array([[float(v) for v in l.split()[1:3]] for l in lines[start:start + 194]]) span = cities.max(0) - cities.min(0) norm = (cities - cities.min(0)) / span.max() # scale into [0, 1], keeping the aspect ratio rng = np.random.default_rng(1) n = len(norm) * 8 # neurons: 8 per city, on a ring network = rng.random((n, 2)) lr, radius = 0.8, n / 10 def route_length(order): pts = cities[order] return np.sqrt(((pts - np.roll(pts, 1, axis=0)) ** 2).sum(1)).sum() for i in range(12000): city = norm[rng.integers(len(norm))] winner = ((network - city) ** 2).sum(1).argmin() d = np.abs(winner - np.arange(n)); d = np.minimum(d, n - d) # distance along the ring g = np.exp(-d ** 2 / (2 * max(radius, 1) ** 2)) network += lr * g[:, None] * (city - network) lr *= 0.99997; radius *= 0.9997 # decay, as in Listing 11.5 if radius < 1 or lr < 0.001: break # read the tour: sort the cities by the ring position of their winning neuron winners = [((network - c) ** 2).sum(1).argmin() for c in norm] order = np.argsort(winners) L = route_length(order) print(f"SOM route length: {L:,.0f} optimal: 9,352 gap: {(L / 9352 - 1):.1%}")
What to remember
Distance matters
The measure decides what “close” means.
K-means is local search
It lowers inertia every step but depends on its seeds.
Validate
Internal indexes need no labels; external ones do.
SOMs
A ring-shaped SOM gives good TSP tours.
References: T. Kohonen, “Self-organized formation of topologically correct feature maps,” 1982. M. Ester, H.-P. Kriegel, J. Sander and X. Xu, “A density-based algorithm for discovering clusters in large spatial databases with noise,” 1996. D. Arthur and S. Vassilvitskii, “k-means++: the advantages of careful seeding,” 2007. TSPLIB national instances (qa194). A. Khamis, Optimization Algorithms: AI techniques for design, planning, and control problems, Manning, 2024, Chapter 11.
Check your understanding
Two quick questions on Chapter 11. Pick an answer to check it.
Q1. Two runs of k-means on the same data with the same k give different clusterings. Why?
Each k-means step lowers the within-cluster sum of squares, so it always converges, but only to a local minimum that depends on where the centers started; k-means++ seeding reduces this effect.Q2. When a ring-shaped self-organizing map is used to solve the TSP, how is the tour obtained?
Training stretches the ring of neurons like an elastic band toward the cities; reading the cities in the ring order of their best matching units gives a tour.For more exercises with full solutions, see Appendix C, “Exercises and solutions,” in Optimization Algorithms. The coding solutions are also on GitHub (Chapter 11).