11 Clustering — First Principles¶

Status: ✅ Complete | Phase: 2 | Prerequisites: Linear Algebra, Probability & Statistics

Goal. Build K-Means, DBSCAN, and GMM/EM from scratch. Visualise convergence, compare with sklearn, and demonstrate failure cases on non-spherical data.

---

Prerequisites¶

• See theory.md for all derivations (K-Means objective, coordinate descent, DBSCAN definitions, EM updates)
• Required foundations:
  • Linear Algebra: Norms & Distances
  • Probability & Statistics

In [1]:

import random

import matplotlib.pyplot as plt
import numpy as np

%load_ext autoreload
%autoreload 2

SEED = 42
random.seed(SEED)
rng = np.random.default_rng(SEED)

plt.rcParams["figure.dpi"] = 90

import random import matplotlib.pyplot as plt import numpy as np %load_ext autoreload %autoreload 2 SEED = 42 random.seed(SEED) rng = np.random.default_rng(SEED) plt.rcParams["figure.dpi"] = 90

1. Problem Setup — WHY¶

We have data with no labels but suspect there are natural groups. Clustering is the task of discovering these groups automatically. Let's create a dataset where the structure is visually obvious.

In [2]:

# Generate 3 well-separated blobs in 2D
centers_true = np.array([[2.0, 2.0], [-2.0, -2.0], [2.0, -2.0]])
n_per_cluster = 100
blobs = []
for center in centers_true:
    blobs.append(center + rng.normal(0, 0.6, size=(n_per_cluster, 2)))
X_blobs = np.vstack(blobs)

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

# Left: unlabelled view (what the algorithm sees)
axes[0].scatter(X_blobs[:, 0], X_blobs[:, 1], alpha=0.5, edgecolor="k", s=20)
axes[0].set_title("What the algorithm sees (no labels)")
axes[0].set_xlabel("$x_1$")
axes[0].set_ylabel("$x_2$")

# Right: true cluster identity
colors_true = np.repeat([0, 1, 2], n_per_cluster)
axes[1].scatter(X_blobs[:, 0], X_blobs[:, 1], c=colors_true, cmap="Set1",
                alpha=0.6, edgecolor="k", s=20)
for i, c in enumerate(centers_true):
    axes[1].plot(*c, "kX", markersize=12)
axes[1].set_title("Ground truth (hidden from algorithm)")
axes[1].set_xlabel("$x_1$")
axes[1].set_ylabel("$x_2$")

plt.tight_layout()
plt.show()

# Generate 3 well-separated blobs in 2D centers_true = np.array([[2.0, 2.0], [-2.0, -2.0], [2.0, -2.0]]) n_per_cluster = 100 blobs = [] for center in centers_true: blobs.append(center + rng.normal(0, 0.6, size=(n_per_cluster, 2))) X_blobs = np.vstack(blobs) fig, axes = plt.subplots(1, 2, figsize=(11, 4.5)) # Left: unlabelled view (what the algorithm sees) axes[0].scatter(X_blobs[:, 0], X_blobs[:, 1], alpha=0.5, edgecolor="k", s=20) axes[0].set_title("What the algorithm sees (no labels)") axes[0].set_xlabel("$x_1$") axes[0].set_ylabel("$x_2$") # Right: true cluster identity colors_true = np.repeat([0, 1, 2], n_per_cluster) axes[1].scatter(X_blobs[:, 0], X_blobs[:, 1], c=colors_true, cmap="Set1", alpha=0.6, edgecolor="k", s=20) for i, c in enumerate(centers_true): axes[1].plot(*c, "kX", markersize=12) axes[1].set_title("Ground truth (hidden from algorithm)") axes[1].set_xlabel("$x_1$") axes[1].set_ylabel("$x_2$") plt.tight_layout() plt.show()

2. Mathematical Core — WHAT¶

K-Means objective¶

Minimise the within-cluster sum of squares (inertia):

$$J(c, \mu) = \sum_{k=1}^K \sum_{i \in C_k} \|x_i - \mu_k\|^2$$

DBSCAN¶

Clusters = maximal sets of density-connected core points, plus reachable border points. Parameters: $\varepsilon$ (radius) and minPts (minimum neighbours).

GMM¶

$$p(x) = \sum_{k=1}^K \pi_k \, \mathcal{N}(x; \mu_k, \Sigma_k)$$

EM alternates between computing responsibilities (E-step) and updating $\mu_k, \Sigma_k, \pi_k$ (M-step).

See theory.md for full derivations.

3. Solution Method — HOW¶

Lloyd's algorithm (K-Means)¶

Coordinate descent on $J$:

1. Assign each $x_i$ to its nearest centroid.
2. Update each $\mu_k$ to the mean of its assigned points.
3. Repeat until convergence.

Let's watch this process unfold step by step.

In [3]:

def kmeans_step_by_step(X, K, max_iter=20, seed=42):
    """Run K-Means and record centroids + labels at every iteration."""
    rng_local = np.random.default_rng(seed)
    idx = rng_local.choice(X.shape[0], K, replace=False)
    centers = X[idx].copy()
    history = [{"centers": centers.copy(), "labels": None}]

    for _ in range(max_iter):
        # Assign
        dists = np.sum((X[:, None, :] - centers[None, :, :]) ** 2, axis=2)
        labels = np.argmin(dists, axis=1)
        # Update
        new_centers = np.array([X[labels == k].mean(axis=0) for k in range(K)])
        history.append({"centers": new_centers.copy(), "labels": labels.copy()})
        if np.allclose(new_centers, centers, atol=1e-6):
            break
        centers = new_centers
    return history

history = kmeans_step_by_step(X_blobs, K=3)
print(f"Converged in {len(history) - 1} iterations")

def kmeans_step_by_step(X, K, max_iter=20, seed=42): """Run K-Means and record centroids + labels at every iteration.""" rng_local = np.random.default_rng(seed) idx = rng_local.choice(X.shape[0], K, replace=False) centers = X[idx].copy() history = [{"centers": centers.copy(), "labels": None}] for _ in range(max_iter): # Assign dists = np.sum((X[:, None, :] - centers[None, :, :]) ** 2, axis=2) labels = np.argmin(dists, axis=1) # Update new_centers = np.array([X[labels == k].mean(axis=0) for k in range(K)]) history.append({"centers": new_centers.copy(), "labels": labels.copy()}) if np.allclose(new_centers, centers, atol=1e-6): break centers = new_centers return history history = kmeans_step_by_step(X_blobs, K=3) print(f"Converged in {len(history) - 1} iterations")

Converged in 3 iterations

In [4]:

# Visualise first 4 iterations
n_show = min(4, len(history) - 1)
fig, axes = plt.subplots(1, n_show, figsize=(4 * n_show, 4), sharey=True)
if n_show == 1:
    axes = [axes]

for step_idx in range(n_show):
    ax = axes[step_idx]
    entry = history[step_idx + 1]
    ax.scatter(X_blobs[:, 0], X_blobs[:, 1], c=entry["labels"],
              cmap="Set1", alpha=0.4, s=15, edgecolor="none")
    # Show previous centroids (faded) and current (solid)
    prev_c = history[step_idx]["centers"]
    curr_c = entry["centers"]
    ax.scatter(prev_c[:, 0], prev_c[:, 1], marker="o", s=100,
              facecolors="none", edgecolors="gray", linewidths=2, label="prev")
    ax.scatter(curr_c[:, 0], curr_c[:, 1], marker="X", s=120,
              c="black", zorder=5, label="current")
    for k in range(3):
        ax.annotate("", xy=curr_c[k], xytext=prev_c[k],
                    arrowprops=dict(arrowstyle="->", color="red", lw=1.5))
    ax.set_title(f"Iteration {step_idx + 1}")
    ax.set_xlabel("$x_1$")
axes[0].set_ylabel("$x_2$")
axes[0].legend(fontsize=8)
plt.suptitle("K-Means convergence: centroids moving to cluster means", y=1.02)
plt.tight_layout()
plt.show()

# Visualise first 4 iterations n_show = min(4, len(history) - 1) fig, axes = plt.subplots(1, n_show, figsize=(4 * n_show, 4), sharey=True) if n_show == 1: axes = [axes] for step_idx in range(n_show): ax = axes[step_idx] entry = history[step_idx + 1] ax.scatter(X_blobs[:, 0], X_blobs[:, 1], c=entry["labels"], cmap="Set1", alpha=0.4, s=15, edgecolor="none") # Show previous centroids (faded) and current (solid) prev_c = history[step_idx]["centers"] curr_c = entry["centers"] ax.scatter(prev_c[:, 0], prev_c[:, 1], marker="o", s=100, facecolors="none", edgecolors="gray", linewidths=2, label="prev") ax.scatter(curr_c[:, 0], curr_c[:, 1], marker="X", s=120, c="black", zorder=5, label="current") for k in range(3): ax.annotate("", xy=curr_c[k], xytext=prev_c[k], arrowprops=dict(arrowstyle="->", color="red", lw=1.5)) ax.set_title(f"Iteration {step_idx + 1}") ax.set_xlabel("$x_1$") axes[0].set_ylabel("$x_2$") axes[0].legend(fontsize=8) plt.suptitle("K-Means convergence: centroids moving to cluster means", y=1.02) plt.tight_layout() plt.show()

4. Implementation — BUILD¶

4.1 K-Means from scratch¶

In [5]:

class KMeansScratch:
    """K-Means with random or K-Means++ initialisation."""

    def __init__(self, n_clusters=3, max_iter=100, tol=1e-4,
                 init="kmeans++", random_state=42):
        self.n_clusters = n_clusters
        self.max_iter = max_iter
        self.tol = tol
        self.init = init
        self.random_state = random_state
        self.cluster_centers_ = None
        self.labels_ = None
        self.inertia_ = None
        self.n_iter_ = 0

    def _init_centers(self, X):
        rng_local = np.random.default_rng(self.random_state)
        n = X.shape[0]
        if self.init == "random":
            idx = rng_local.choice(n, self.n_clusters, replace=False)
            return X[idx].copy()
        # K-Means++
        centers = [X[rng_local.integers(n)]]
        for _ in range(1, self.n_clusters):
            dists = np.min([np.sum((X - c) ** 2, axis=1) for c in centers], axis=0)
            probs = dists / dists.sum()
            idx = rng_local.choice(n, p=probs)
            centers.append(X[idx])
        return np.array(centers)

    def fit(self, X):
        X = np.asarray(X, dtype=float)
        centers = self._init_centers(X)

        for iteration in range(1, self.max_iter + 1):
            # Assign
            sq_dists = np.sum((X[:, None, :] - centers[None, :, :]) ** 2, axis=2)
            labels = np.argmin(sq_dists, axis=1)

            # Update
            new_centers = np.empty_like(centers)
            for k in range(self.n_clusters):
                members = X[labels == k]
                if len(members) > 0:
                    new_centers[k] = members.mean(axis=0)
                else:
                    new_centers[k] = centers[k]  # keep old center

            shift = np.max(np.linalg.norm(new_centers - centers, axis=1))
            centers = new_centers
            self.n_iter_ = iteration
            if shift <= self.tol:
                break

        # Final assignment
        sq_dists = np.sum((X[:, None, :] - centers[None, :, :]) ** 2, axis=2)
        self.labels_ = np.argmin(sq_dists, axis=1)
        self.cluster_centers_ = centers
        self.inertia_ = float(np.sum(sq_dists[np.arange(X.shape[0]), self.labels_]))
        return self

    def predict(self, X):
        X = np.asarray(X, dtype=float)
        sq_dists = np.sum((X[:, None, :] - self.cluster_centers_[None, :, :]) ** 2, axis=2)
        return np.argmin(sq_dists, axis=1)

print("KMeansScratch defined ✓")

class KMeansScratch: """K-Means with random or K-Means++ initialisation.""" def __init__(self, n_clusters=3, max_iter=100, tol=1e-4, init="kmeans++", random_state=42): self.n_clusters = n_clusters self.max_iter = max_iter self.tol = tol self.init = init self.random_state = random_state self.cluster_centers_ = None self.labels_ = None self.inertia_ = None self.n_iter_ = 0 def _init_centers(self, X): rng_local = np.random.default_rng(self.random_state) n = X.shape[0] if self.init == "random": idx = rng_local.choice(n, self.n_clusters, replace=False) return X[idx].copy() # K-Means++ centers = [X[rng_local.integers(n)]] for _ in range(1, self.n_clusters): dists = np.min([np.sum((X - c) ** 2, axis=1) for c in centers], axis=0) probs = dists / dists.sum() idx = rng_local.choice(n, p=probs) centers.append(X[idx]) return np.array(centers) def fit(self, X): X = np.asarray(X, dtype=float) centers = self._init_centers(X) for iteration in range(1, self.max_iter + 1): # Assign sq_dists = np.sum((X[:, None, :] - centers[None, :, :]) ** 2, axis=2) labels = np.argmin(sq_dists, axis=1) # Update new_centers = np.empty_like(centers) for k in range(self.n_clusters): members = X[labels == k] if len(members) > 0: new_centers[k] = members.mean(axis=0) else: new_centers[k] = centers[k] # keep old center shift = np.max(np.linalg.norm(new_centers - centers, axis=1)) centers = new_centers self.n_iter_ = iteration if shift <= self.tol: break # Final assignment sq_dists = np.sum((X[:, None, :] - centers[None, :, :]) ** 2, axis=2) self.labels_ = np.argmin(sq_dists, axis=1) self.cluster_centers_ = centers self.inertia_ = float(np.sum(sq_dists[np.arange(X.shape[0]), self.labels_])) return self def predict(self, X): X = np.asarray(X, dtype=float) sq_dists = np.sum((X[:, None, :] - self.cluster_centers_[None, :, :]) ** 2, axis=2) return np.argmin(sq_dists, axis=1) print("KMeansScratch defined ✓")

KMeansScratch defined ✓

In [6]:

# Quick test on blobs
km = KMeansScratch(n_clusters=3, init="kmeans++", random_state=42)
km.fit(X_blobs)
print(f"Iterations: {km.n_iter_}, Inertia: {km.inertia_:.2f}")

fig, ax = plt.subplots(figsize=(5, 4.5))
ax.scatter(X_blobs[:, 0], X_blobs[:, 1], c=km.labels_, cmap="Set1",
           alpha=0.5, s=20, edgecolor="none")
ax.scatter(km.cluster_centers_[:, 0], km.cluster_centers_[:, 1],
           c="black", marker="X", s=150, zorder=5, label="centroids")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.set_title("K-Means++ on blobs")
ax.legend()
plt.tight_layout()
plt.show()

# Quick test on blobs km = KMeansScratch(n_clusters=3, init="kmeans++", random_state=42) km.fit(X_blobs) print(f"Iterations: {km.n_iter_}, Inertia: {km.inertia_:.2f}") fig, ax = plt.subplots(figsize=(5, 4.5)) ax.scatter(X_blobs[:, 0], X_blobs[:, 1], c=km.labels_, cmap="Set1", alpha=0.5, s=20, edgecolor="none") ax.scatter(km.cluster_centers_[:, 0], km.cluster_centers_[:, 1], c="black", marker="X", s=150, zorder=5, label="centroids") ax.set_xlabel("$x_1$") ax.set_ylabel("$x_2$") ax.set_title("K-Means++ on blobs") ax.legend() plt.tight_layout() plt.show()

Iterations: 2, Inertia: 203.95

4.2 Elbow method for choosing K¶

In [7]:

K_range = range(1, 9)
inertias = []
for K in K_range:
    model = KMeansScratch(n_clusters=K, random_state=42)
    model.fit(X_blobs)
    inertias.append(model.inertia_)

fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(list(K_range), inertias, "bo-", linewidth=2)
ax.axvline(3, color="red", linestyle="--", label="K=3 (elbow)")
ax.set_xlabel("Number of clusters K")
ax.set_ylabel("Inertia (WCSS)")
ax.set_title("Elbow method — inertia vs K")
ax.legend()
plt.tight_layout()
plt.show()

K_range = range(1, 9) inertias = [] for K in K_range: model = KMeansScratch(n_clusters=K, random_state=42) model.fit(X_blobs) inertias.append(model.inertia_) fig, ax = plt.subplots(figsize=(6, 4)) ax.plot(list(K_range), inertias, "bo-", linewidth=2) ax.axvline(3, color="red", linestyle="--", label="K=3 (elbow)") ax.set_xlabel("Number of clusters K") ax.set_ylabel("Inertia (WCSS)") ax.set_title("Elbow method — inertia vs K") ax.legend() plt.tight_layout() plt.show()

4.3 DBSCAN from scratch¶

In [8]:

class DBSCANScratch:
    """DBSCAN clustering with brute-force neighbour search."""

    def __init__(self, eps=0.5, min_samples=5):
        self.eps = eps
        self.min_samples = min_samples
        self.labels_ = None

    def fit(self, X):
        X = np.asarray(X, dtype=float)
        n = X.shape[0]
        labels = np.full(n, -1, dtype=int)  # -1 = unvisited/noise
        cluster_id = 0

        # Precompute pairwise distances
        dist_matrix = np.sqrt(np.sum((X[:, None, :] - X[None, :, :]) ** 2, axis=2))

        def region_query(i):
            return np.where(dist_matrix[i] <= self.eps)[0]

        visited = np.zeros(n, dtype=bool)

        for i in range(n):
            if visited[i]:
                continue
            visited[i] = True
            neighbours = region_query(i)

            if len(neighbours) < self.min_samples:
                labels[i] = -1  # noise (may be relabelled as border)
                continue

            # Start new cluster
            labels[i] = cluster_id
            seed_set = list(neighbours)
            j = 0
            while j < len(seed_set):
                q = seed_set[j]
                if not visited[q]:
                    visited[q] = True
                    q_neighbours = region_query(q)
                    if len(q_neighbours) >= self.min_samples:
                        seed_set.extend(q_neighbours.tolist())
                if labels[q] == -1:  # was noise or unvisited
                    labels[q] = cluster_id
                j += 1

            cluster_id += 1

        self.labels_ = labels
        return self

print("DBSCANScratch defined ✓")

class DBSCANScratch: """DBSCAN clustering with brute-force neighbour search.""" def __init__(self, eps=0.5, min_samples=5): self.eps = eps self.min_samples = min_samples self.labels_ = None def fit(self, X): X = np.asarray(X, dtype=float) n = X.shape[0] labels = np.full(n, -1, dtype=int) # -1 = unvisited/noise cluster_id = 0 # Precompute pairwise distances dist_matrix = np.sqrt(np.sum((X[:, None, :] - X[None, :, :]) ** 2, axis=2)) def region_query(i): return np.where(dist_matrix[i] <= self.eps)[0] visited = np.zeros(n, dtype=bool) for i in range(n): if visited[i]: continue visited[i] = True neighbours = region_query(i) if len(neighbours) < self.min_samples: labels[i] = -1 # noise (may be relabelled as border) continue # Start new cluster labels[i] = cluster_id seed_set = list(neighbours) j = 0 while j < len(seed_set): q = seed_set[j] if not visited[q]: visited[q] = True q_neighbours = region_query(q) if len(q_neighbours) >= self.min_samples: seed_set.extend(q_neighbours.tolist()) if labels[q] == -1: # was noise or unvisited labels[q] = cluster_id j += 1 cluster_id += 1 self.labels_ = labels return self print("DBSCANScratch defined ✓")

DBSCANScratch defined ✓

In [9]:

db = DBSCANScratch(eps=1.0, min_samples=5)
db.fit(X_blobs)
n_clusters_found = len(set(db.labels_) - {-1})
n_noise = np.sum(db.labels_ == -1)
print(f"Clusters found: {n_clusters_found}, Noise points: {n_noise}")

fig, ax = plt.subplots(figsize=(5, 4.5))
mask = db.labels_ >= 0
ax.scatter(X_blobs[mask, 0], X_blobs[mask, 1], c=db.labels_[mask],
           cmap="Set1", alpha=0.5, s=20, edgecolor="none", label="clustered")
if n_noise > 0:
    ax.scatter(X_blobs[~mask, 0], X_blobs[~mask, 1], c="gray",
              marker="x", s=30, label="noise")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.set_title(f"DBSCAN (eps={db.eps}, minPts={db.min_samples})")
ax.legend()
plt.tight_layout()
plt.show()

db = DBSCANScratch(eps=1.0, min_samples=5) db.fit(X_blobs) n_clusters_found = len(set(db.labels_) - {-1}) n_noise = np.sum(db.labels_ == -1) print(f"Clusters found: {n_clusters_found}, Noise points: {n_noise}") fig, ax = plt.subplots(figsize=(5, 4.5)) mask = db.labels_ >= 0 ax.scatter(X_blobs[mask, 0], X_blobs[mask, 1], c=db.labels_[mask], cmap="Set1", alpha=0.5, s=20, edgecolor="none", label="clustered") if n_noise > 0: ax.scatter(X_blobs[~mask, 0], X_blobs[~mask, 1], c="gray", marker="x", s=30, label="noise") ax.set_xlabel("$x_1$") ax.set_ylabel("$x_2$") ax.set_title(f"DBSCAN (eps={db.eps}, minPts={db.min_samples})") ax.legend() plt.tight_layout() plt.show()

Clusters found: 2, Noise points: 1

4.4 Gaussian Mixture Model with EM from scratch¶

In [10]:

class GMMScratch:
    """Gaussian Mixture Model fitted by EM."""

    def __init__(self, n_components=3, max_iter=100, tol=1e-6,
                 reg_covar=1e-6, random_state=42):
        self.n_components = n_components
        self.max_iter = max_iter
        self.tol = tol
        self.reg_covar = reg_covar
        self.random_state = random_state

    @staticmethod
    def _multivariate_gaussian(X, mu, sigma):
        """Evaluate N(x; mu, sigma) for each row of X."""
        d = X.shape[1]
        diff = X - mu  # (n, d)
        sigma_inv = np.linalg.inv(sigma)
        det_sigma = np.linalg.det(sigma)
        # Mahalanobis: each row diff[i] @ sigma_inv @ diff[i]
        maha = np.sum(diff @ sigma_inv * diff, axis=1)
        norm_const = np.sqrt((2 * np.pi) ** d * det_sigma)
        return np.exp(-0.5 * maha) / norm_const

    def fit(self, X):
        X = np.asarray(X, dtype=float)
        n, d = X.shape
        K = self.n_components
        rng_local = np.random.default_rng(self.random_state)

        # Initialise with K-Means
        km_init = KMeansScratch(n_clusters=K, random_state=self.random_state)
        km_init.fit(X)
        self.means_ = km_init.cluster_centers_.copy()
        self.covariances_ = np.array([np.eye(d) for _ in range(K)])
        self.weights_ = np.ones(K) / K

        self.log_likelihoods_ = []

        for iteration in range(self.max_iter):
            # === E-step: compute responsibilities ===
            resp = np.zeros((n, K))
            for k in range(K):
                resp[:, k] = self.weights_[k] * self._multivariate_gaussian(
                    X, self.means_[k], self.covariances_[k]
                )
            resp_sum = resp.sum(axis=1, keepdims=True)
            resp_sum = np.maximum(resp_sum, 1e-300)  # avoid division by zero
            resp /= resp_sum

            # Log-likelihood
            ll = float(np.sum(np.log(resp_sum.ravel())))
            self.log_likelihoods_.append(ll)
            if len(self.log_likelihoods_) > 1:
                if abs(self.log_likelihoods_[-1] - self.log_likelihoods_[-2]) < self.tol:
                    break

            # === M-step: update parameters ===
            Nk = resp.sum(axis=0)  # (K,)
            for k in range(K):
                # Mean
                self.means_[k] = (resp[:, k] @ X) / Nk[k]
                # Covariance
                diff = X - self.means_[k]  # (n, d)
                self.covariances_[k] = (
                    (diff * resp[:, k:k+1]).T @ diff / Nk[k]
                    + self.reg_covar * np.eye(d)
                )
                # Weight
                self.weights_[k] = Nk[k] / n

        self.responsibilities_ = resp
        self.labels_ = np.argmax(resp, axis=1)
        self.n_iter_ = iteration + 1
        return self

    def predict(self, X):
        X = np.asarray(X, dtype=float)
        K = self.n_components
        resp = np.zeros((X.shape[0], K))
        for k in range(K):
            resp[:, k] = self.weights_[k] * self._multivariate_gaussian(
                X, self.means_[k], self.covariances_[k]
            )
        return np.argmax(resp, axis=1)

print("GMMScratch defined ✓")

class GMMScratch: """Gaussian Mixture Model fitted by EM.""" def __init__(self, n_components=3, max_iter=100, tol=1e-6, reg_covar=1e-6, random_state=42): self.n_components = n_components self.max_iter = max_iter self.tol = tol self.reg_covar = reg_covar self.random_state = random_state @staticmethod def _multivariate_gaussian(X, mu, sigma): """Evaluate N(x; mu, sigma) for each row of X.""" d = X.shape[1] diff = X - mu # (n, d) sigma_inv = np.linalg.inv(sigma) det_sigma = np.linalg.det(sigma) # Mahalanobis: each row diff[i] @ sigma_inv @ diff[i] maha = np.sum(diff @ sigma_inv * diff, axis=1) norm_const = np.sqrt((2 * np.pi) ** d * det_sigma) return np.exp(-0.5 * maha) / norm_const def fit(self, X): X = np.asarray(X, dtype=float) n, d = X.shape K = self.n_components rng_local = np.random.default_rng(self.random_state) # Initialise with K-Means km_init = KMeansScratch(n_clusters=K, random_state=self.random_state) km_init.fit(X) self.means_ = km_init.cluster_centers_.copy() self.covariances_ = np.array([np.eye(d) for _ in range(K)]) self.weights_ = np.ones(K) / K self.log_likelihoods_ = [] for iteration in range(self.max_iter): # === E-step: compute responsibilities === resp = np.zeros((n, K)) for k in range(K): resp[:, k] = self.weights_[k] * self._multivariate_gaussian( X, self.means_[k], self.covariances_[k] ) resp_sum = resp.sum(axis=1, keepdims=True) resp_sum = np.maximum(resp_sum, 1e-300) # avoid division by zero resp /= resp_sum # Log-likelihood ll = float(np.sum(np.log(resp_sum.ravel()))) self.log_likelihoods_.append(ll) if len(self.log_likelihoods_) > 1: if abs(self.log_likelihoods_[-1] - self.log_likelihoods_[-2]) < self.tol: break # === M-step: update parameters === Nk = resp.sum(axis=0) # (K,) for k in range(K): # Mean self.means_[k] = (resp[:, k] @ X) / Nk[k] # Covariance diff = X - self.means_[k] # (n, d) self.covariances_[k] = ( (diff * resp[:, k:k+1]).T @ diff / Nk[k] + self.reg_covar * np.eye(d) ) # Weight self.weights_[k] = Nk[k] / n self.responsibilities_ = resp self.labels_ = np.argmax(resp, axis=1) self.n_iter_ = iteration + 1 return self def predict(self, X): X = np.asarray(X, dtype=float) K = self.n_components resp = np.zeros((X.shape[0], K)) for k in range(K): resp[:, k] = self.weights_[k] * self._multivariate_gaussian( X, self.means_[k], self.covariances_[k] ) return np.argmax(resp, axis=1) print("GMMScratch defined ✓")

GMMScratch defined ✓

In [11]:

gmm = GMMScratch(n_components=3, random_state=42)
gmm.fit(X_blobs)
print(f"EM iterations: {gmm.n_iter_}")
print(f"Final log-likelihood: {gmm.log_likelihoods_[-1]:.2f}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

# Left: cluster assignments
axes[0].scatter(X_blobs[:, 0], X_blobs[:, 1], c=gmm.labels_, cmap="Set1",
                alpha=0.5, s=20, edgecolor="none")
for k in range(3):
    axes[0].plot(*gmm.means_[k], "kX", markersize=12)
axes[0].set_title("GMM cluster assignments")
axes[0].set_xlabel("$x_1$")
axes[0].set_ylabel("$x_2$")

# Right: log-likelihood convergence
axes[1].plot(gmm.log_likelihoods_, "b-o", markersize=3)
axes[1].set_xlabel("EM iteration")
axes[1].set_ylabel("Log-likelihood")
axes[1].set_title("EM convergence (monotonically increasing)")

plt.tight_layout()
plt.show()

gmm = GMMScratch(n_components=3, random_state=42) gmm.fit(X_blobs) print(f"EM iterations: {gmm.n_iter_}") print(f"Final log-likelihood: {gmm.log_likelihoods_[-1]:.2f}") fig, axes = plt.subplots(1, 2, figsize=(11, 4.5)) # Left: cluster assignments axes[0].scatter(X_blobs[:, 0], X_blobs[:, 1], c=gmm.labels_, cmap="Set1", alpha=0.5, s=20, edgecolor="none") for k in range(3): axes[0].plot(*gmm.means_[k], "kX", markersize=12) axes[0].set_title("GMM cluster assignments") axes[0].set_xlabel("$x_1$") axes[0].set_ylabel("$x_2$") # Right: log-likelihood convergence axes[1].plot(gmm.log_likelihoods_, "b-o", markersize=3) axes[1].set_xlabel("EM iteration") axes[1].set_ylabel("Log-likelihood") axes[1].set_title("EM convergence (monotonically increasing)") plt.tight_layout() plt.show()

EM iterations: 6
Final log-likelihood: -847.79

5. Library Comparison¶

Compare our from-scratch implementations with:

• ml_first_principles.distance_models.KMeans (src library)
• sklearn.cluster.KMeans and sklearn.cluster.DBSCAN

In [12]:

from ml_first_principles.distance_models import KMeans as LibKMeans
from sklearn.cluster import KMeans as SklearnKMeans

# Our from-scratch
ours = KMeansScratch(n_clusters=3, init="random", random_state=42)
ours.fit(X_blobs)

# Library
lib = LibKMeans(n_clusters=3, random_state=42)
lib.fit(X_blobs)

# Sklearn
skl = SklearnKMeans(n_clusters=3, random_state=42, n_init=1)
skl.fit(X_blobs)

# Since cluster label indices may differ, compare by matching centroids
def match_labels(labels_a, labels_b, K):
    """Find label permutation that maximises agreement."""
    from itertools import permutations
    best_acc = 0
    best_perm = None
    for perm in permutations(range(K)):
        remapped = np.array([perm[l] for l in labels_b])
        acc = np.mean(labels_a == remapped)
        if acc > best_acc:
            best_acc = acc
            best_perm = perm
    return best_acc, best_perm

acc_lib, _ = match_labels(ours.labels_, lib.labels_, 3)
acc_skl, _ = match_labels(ours.labels_, skl.labels_, 3)

print(f"Agreement (ours vs lib):     {acc_lib:.4f}")
print(f"Agreement (ours vs sklearn): {acc_skl:.4f}")
print(f"\nInertia — ours: {ours.inertia_:.2f}, lib: {lib.inertia_:.2f}, "
      f"sklearn: {skl.inertia_:.2f}")

from ml_first_principles.distance_models import KMeans as LibKMeans from sklearn.cluster import KMeans as SklearnKMeans # Our from-scratch ours = KMeansScratch(n_clusters=3, init="random", random_state=42) ours.fit(X_blobs) # Library lib = LibKMeans(n_clusters=3, random_state=42) lib.fit(X_blobs) # Sklearn skl = SklearnKMeans(n_clusters=3, random_state=42, n_init=1) skl.fit(X_blobs) # Since cluster label indices may differ, compare by matching centroids def match_labels(labels_a, labels_b, K): """Find label permutation that maximises agreement.""" from itertools import permutations best_acc = 0 best_perm = None for perm in permutations(range(K)): remapped = np.array([perm[l] for l in labels_b]) acc = np.mean(labels_a == remapped) if acc > best_acc: best_acc = acc best_perm = perm return best_acc, best_perm acc_lib, _ = match_labels(ours.labels_, lib.labels_, 3) acc_skl, _ = match_labels(ours.labels_, skl.labels_, 3) print(f"Agreement (ours vs lib): {acc_lib:.4f}") print(f"Agreement (ours vs sklearn): {acc_skl:.4f}") print(f"\nInertia — ours: {ours.inertia_:.2f}, lib: {lib.inertia_:.2f}, " f"sklearn: {skl.inertia_:.2f}")

Agreement (ours vs lib):     1.0000
Agreement (ours vs sklearn): 1.0000

Inertia — ours: 203.95, lib: 203.95, sklearn: 203.95

In [13]:

from sklearn.cluster import DBSCAN as SklearnDBSCAN

db_ours = DBSCANScratch(eps=1.0, min_samples=5)
db_ours.fit(X_blobs)

db_skl = SklearnDBSCAN(eps=1.0, min_samples=5)
db_skl.fit(X_blobs)

n_ours = len(set(db_ours.labels_) - {-1})
n_skl = len(set(db_skl.labels_) - {-1})
noise_ours = np.sum(db_ours.labels_ == -1)
noise_skl = np.sum(db_skl.labels_ == -1)

# Compare cluster counts and noise points
print(f"Clusters — ours: {n_ours}, sklearn: {n_skl}")
print(f"Noise pts — ours: {noise_ours}, sklearn: {noise_skl}")

# Agreement (handle label permutation)
mask = (db_ours.labels_ >= 0) & (db_skl.labels_ >= 0)
if mask.sum() > 0:
    acc_db, _ = match_labels(db_ours.labels_[mask], db_skl.labels_[mask],
                             max(n_ours, n_skl))
    print(f"Label agreement (non-noise): {acc_db:.4f}")

from sklearn.cluster import DBSCAN as SklearnDBSCAN db_ours = DBSCANScratch(eps=1.0, min_samples=5) db_ours.fit(X_blobs) db_skl = SklearnDBSCAN(eps=1.0, min_samples=5) db_skl.fit(X_blobs) n_ours = len(set(db_ours.labels_) - {-1}) n_skl = len(set(db_skl.labels_) - {-1}) noise_ours = np.sum(db_ours.labels_ == -1) noise_skl = np.sum(db_skl.labels_ == -1) # Compare cluster counts and noise points print(f"Clusters — ours: {n_ours}, sklearn: {n_skl}") print(f"Noise pts — ours: {noise_ours}, sklearn: {noise_skl}") # Agreement (handle label permutation) mask = (db_ours.labels_ >= 0) & (db_skl.labels_ >= 0) if mask.sum() > 0: acc_db, _ = match_labels(db_ours.labels_[mask], db_skl.labels_[mask], max(n_ours, n_skl)) print(f"Label agreement (non-noise): {acc_db:.4f}")

Clusters — ours: 2, sklearn: 2
Noise pts — ours: 1, sklearn: 1
Label agreement (non-noise): 1.0000

6. Experiments and Failures — VERIFY¶

6.1 Non-spherical data: moons and circles¶

K-Means assumes spherical clusters (Voronoi cells). Let's see what happens on non-convex shapes.

In [14]:

from sklearn.datasets import make_moons, make_circles

X_moons, y_moons = make_moons(n_samples=300, noise=0.08, random_state=42)
X_circles, y_circles = make_circles(n_samples=300, noise=0.05,
                                     factor=0.5, random_state=42)

from sklearn.datasets import make_moons, make_circles X_moons, y_moons = make_moons(n_samples=300, noise=0.08, random_state=42) X_circles, y_circles = make_circles(n_samples=300, noise=0.05, factor=0.5, random_state=42)

In [15]:

datasets = [("Moons", X_moons, y_moons), ("Circles", X_circles, y_circles)]
methods = [
    ("K-Means", lambda X: KMeansScratch(n_clusters=2, random_state=42).fit(X).labels_),
    ("DBSCAN", lambda X: DBSCANScratch(eps=0.3, min_samples=5).fit(X).labels_),
    ("GMM", lambda X: GMMScratch(n_components=2, random_state=42).fit(X).labels_),
]

fig, axes = plt.subplots(2, 4, figsize=(16, 7))

for row, (data_name, X_data, y_data) in enumerate(datasets):
    # Ground truth
    axes[row, 0].scatter(X_data[:, 0], X_data[:, 1], c=y_data,
                         cmap="Set1", s=15, alpha=0.7)
    axes[row, 0].set_title(f"{data_name} — Ground truth")
    axes[row, 0].set_ylabel("$x_2$")

    for col, (method_name, method_fn) in enumerate(methods, start=1):
        labels = method_fn(X_data)
        axes[row, col].scatter(X_data[:, 0], X_data[:, 1], c=labels,
                               cmap="Set1", s=15, alpha=0.7)
        axes[row, col].set_title(f"{data_name} — {method_name}")

for ax in axes.flat:
    ax.set_xlabel("$x_1$")

plt.suptitle("Clustering comparison on non-spherical data", fontsize=13, y=1.01)
plt.tight_layout()
plt.show()

print("Observation: K-Means and GMM fail on moons/circles because they")
print("assume convex (spherical/ellipsoidal) cluster shapes.")
print("DBSCAN succeeds because it follows density connectivity.")

datasets = [("Moons", X_moons, y_moons), ("Circles", X_circles, y_circles)] methods = [ ("K-Means", lambda X: KMeansScratch(n_clusters=2, random_state=42).fit(X).labels_), ("DBSCAN", lambda X: DBSCANScratch(eps=0.3, min_samples=5).fit(X).labels_), ("GMM", lambda X: GMMScratch(n_components=2, random_state=42).fit(X).labels_), ] fig, axes = plt.subplots(2, 4, figsize=(16, 7)) for row, (data_name, X_data, y_data) in enumerate(datasets): # Ground truth axes[row, 0].scatter(X_data[:, 0], X_data[:, 1], c=y_data, cmap="Set1", s=15, alpha=0.7) axes[row, 0].set_title(f"{data_name} — Ground truth") axes[row, 0].set_ylabel("$x_2$") for col, (method_name, method_fn) in enumerate(methods, start=1): labels = method_fn(X_data) axes[row, col].scatter(X_data[:, 0], X_data[:, 1], c=labels, cmap="Set1", s=15, alpha=0.7) axes[row, col].set_title(f"{data_name} — {method_name}") for ax in axes.flat: ax.set_xlabel("$x_1$") plt.suptitle("Clustering comparison on non-spherical data", fontsize=13, y=1.01) plt.tight_layout() plt.show() print("Observation: K-Means and GMM fail on moons/circles because they") print("assume convex (spherical/ellipsoidal) cluster shapes.") print("DBSCAN succeeds because it follows density connectivity.")

Observation: K-Means and GMM fail on moons/circles because they
assume convex (spherical/ellipsoidal) cluster shapes.
DBSCAN succeeds because it follows density connectivity.

6.2 Failure case: K-Means on elongated clusters¶

In [16]:

# Create two elongated, angled clusters
n_elong = 150
t = rng.uniform(0, 1, n_elong)
cluster_a = np.column_stack([t * 4, t * 4 + rng.normal(0, 0.2, n_elong)])
cluster_b = np.column_stack([t * 4 + 1, -t * 4 + 4 + rng.normal(0, 0.2, n_elong)])
X_elong = np.vstack([cluster_a, cluster_b])
y_elong = np.array([0] * n_elong + [1] * n_elong)

km_elong = KMeansScratch(n_clusters=2, random_state=42).fit(X_elong)
db_elong = DBSCANScratch(eps=0.5, min_samples=5).fit(X_elong)
gmm_elong = GMMScratch(n_components=2, random_state=42).fit(X_elong)

fig, axes = plt.subplots(1, 4, figsize=(16, 3.5))

for ax, (title, labels) in zip(axes, [
    ("Ground truth", y_elong),
    ("K-Means", km_elong.labels_),
    ("DBSCAN", db_elong.labels_),
    ("GMM", gmm_elong.labels_),
]):
    ax.scatter(X_elong[:, 0], X_elong[:, 1], c=labels, cmap="Set1",
              s=15, alpha=0.7)
    ax.set_title(title)
    ax.set_xlabel("$x_1$")
    ax.set_ylabel("$x_2$")

plt.suptitle("Elongated clusters — K-Means splits by centroid, misses the shape",
             fontsize=11, y=1.02)
plt.tight_layout()
plt.show()

# Create two elongated, angled clusters n_elong = 150 t = rng.uniform(0, 1, n_elong) cluster_a = np.column_stack([t * 4, t * 4 + rng.normal(0, 0.2, n_elong)]) cluster_b = np.column_stack([t * 4 + 1, -t * 4 + 4 + rng.normal(0, 0.2, n_elong)]) X_elong = np.vstack([cluster_a, cluster_b]) y_elong = np.array([0] * n_elong + [1] * n_elong) km_elong = KMeansScratch(n_clusters=2, random_state=42).fit(X_elong) db_elong = DBSCANScratch(eps=0.5, min_samples=5).fit(X_elong) gmm_elong = GMMScratch(n_components=2, random_state=42).fit(X_elong) fig, axes = plt.subplots(1, 4, figsize=(16, 3.5)) for ax, (title, labels) in zip(axes, [ ("Ground truth", y_elong), ("K-Means", km_elong.labels_), ("DBSCAN", db_elong.labels_), ("GMM", gmm_elong.labels_), ]): ax.scatter(X_elong[:, 0], X_elong[:, 1], c=labels, cmap="Set1", s=15, alpha=0.7) ax.set_title(title) ax.set_xlabel("$x_1$") ax.set_ylabel("$x_2$") plt.suptitle("Elongated clusters — K-Means splits by centroid, misses the shape", fontsize=11, y=1.02) plt.tight_layout() plt.show()

6.3 K-Means++ vs random initialisation¶

In [17]:

inertias_random = []
inertias_pp = []
for seed in range(50):
    km_r = KMeansScratch(n_clusters=3, init="random", random_state=seed)
    km_r.fit(X_blobs)
    inertias_random.append(km_r.inertia_)

    km_p = KMeansScratch(n_clusters=3, init="kmeans++", random_state=seed)
    km_p.fit(X_blobs)
    inertias_pp.append(km_p.inertia_)

fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(inertias_random, bins=15, alpha=0.6, label="Random init", color="steelblue")
ax.hist(inertias_pp, bins=15, alpha=0.6, label="K-Means++ init", color="coral")
ax.axvline(np.median(inertias_random), color="steelblue", linestyle="--",
           label=f"Random median: {np.median(inertias_random):.0f}")
ax.axvline(np.median(inertias_pp), color="coral", linestyle="--",
           label=f"K-Means++ median: {np.median(inertias_pp):.0f}")
ax.set_xlabel("Final inertia")
ax.set_ylabel("Count")
ax.set_title("K-Means++: more consistent, lower inertia (50 runs)")
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

inertias_random = [] inertias_pp = [] for seed in range(50): km_r = KMeansScratch(n_clusters=3, init="random", random_state=seed) km_r.fit(X_blobs) inertias_random.append(km_r.inertia_) km_p = KMeansScratch(n_clusters=3, init="kmeans++", random_state=seed) km_p.fit(X_blobs) inertias_pp.append(km_p.inertia_) fig, ax = plt.subplots(figsize=(7, 4)) ax.hist(inertias_random, bins=15, alpha=0.6, label="Random init", color="steelblue") ax.hist(inertias_pp, bins=15, alpha=0.6, label="K-Means++ init", color="coral") ax.axvline(np.median(inertias_random), color="steelblue", linestyle="--", label=f"Random median: {np.median(inertias_random):.0f}") ax.axvline(np.median(inertias_pp), color="coral", linestyle="--", label=f"K-Means++ median: {np.median(inertias_pp):.0f}") ax.set_xlabel("Final inertia") ax.set_ylabel("Count") ax.set_title("K-Means++: more consistent, lower inertia (50 runs)") ax.legend(fontsize=8) plt.tight_layout() plt.show()

7. Connections¶

• Prerequisites: Linear Algebra, Probability & Statistics
• Theory: theory.md — all derivations live there
• Related Topics: 10 PCA (reduce dimensions before clustering), 08 KNN (distance-based supervised)
• Synthesis: Supervised vs. Unsupervised
• Graph Map: See INDEX.md

Takeaway¶

• K-Means minimises within-cluster sum of squares via coordinate descent (assign → update). Fast and simple, but assumes spherical clusters and requires choosing $K$.
• K-Means++ initialisation consistently finds better local minima than random init.
• DBSCAN discovers clusters of arbitrary shape by following density connectivity. No need to specify $K$, but sensitive to the $\varepsilon$ parameter and cannot handle varying density.
• GMM/EM provides soft probabilistic assignments and models ellipsoidal clusters. EM monotonically increases log-likelihood but converges to local maxima.
• No single method dominates — the right choice depends on cluster shape, noise, and whether you need soft assignments.

Exercises¶

See exercises.ipynb for practice problems.