Chapter 3 — Unsupervised Learning¶

import numpy as np
from typing import NamedTuple


class KMeansResult(NamedTuple):
    labels: np.ndarray
    centroids: np.ndarray
    inertia: float
    n_iter: int


def kmeans_plusplus_init(X: np.ndarray, K: int, rng: np.random.Generator) -> np.ndarray:
    """K-Means++ centroid initialisation."""
    n = X.shape[0]
    first_idx = rng.integers(n)
    centroids = [X[first_idx]]

    for _ in range(K - 1):
        dists = np.array([min(np.sum((x - c) ** 2) for c in centroids) for x in X])
        probs = dists / dists.sum()
        idx = rng.choice(n, p=probs)
        centroids.append(X[idx])

    return np.array(centroids)


def kmeans(
    X: np.ndarray,
    K: int,
    max_iter: int = 300,
    tol: float = 1e-4,
    random_state: int = 42,
) -> KMeansResult:
    """K-Means with K-Means++ initialisation."""
    rng = np.random.default_rng(random_state)
    centroids = kmeans_plusplus_init(X, K, rng)

    for iteration in range(max_iter):
        # Assignment step
        dists = np.linalg.norm(X[:, None, :] - centroids[None, :, :], axis=2)
        labels = dists.argmin(axis=1)

        # Update step
        new_centroids = np.array([
            X[labels == k].mean(axis=0) if (labels == k).any() else centroids[k]
            for k in range(K)
        ])

        shift = np.linalg.norm(new_centroids - centroids)
        centroids = new_centroids
        if shift < tol:
            break

    inertia = sum(np.sum((X[labels == k] - centroids[k]) ** 2) for k in range(K))
    return KMeansResult(labels=labels, centroids=centroids, inertia=float(inertia), n_iter=iteration + 1)

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
import matplotlib.pyplot as plt

K_range = range(2, 15)
inertias: list[float] = []
silhouettes: list[float] = []

for k in K_range:
    km = KMeans(n_clusters=k, init="k-means++", n_init=10, random_state=42)
    labels = km.fit_predict(X)
    inertias.append(km.inertia_)
    silhouettes.append(silhouette_score(X, labels))

best_k_silhouette = K_range[silhouettes.index(max(silhouettes))]
print(f"Best K by silhouette: {best_k_silhouette} (score={max(silhouettes):.4f})")

from sklearn.cluster import DBSCAN
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

dbscan = DBSCAN(eps=0.5, min_samples=5)
labels = dbscan.fit_predict(X_scaled)

n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
n_noise = (labels == -1).sum()
print(f"Clusters: {n_clusters}, Noise points: {n_noise} ({100*n_noise/len(labels):.1f}%)")

from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram, linkage
import matplotlib.pyplot as plt

# Compute full linkage matrix for dendrogram
Z = linkage(X_scaled, method="ward")

fig, ax = plt.subplots(figsize=(12, 5))
dendrogram(Z, ax=ax, truncate_mode="level", p=5, leaf_font_size=10)
ax.set_title("Dendrogram (Ward linkage, top 5 levels)")
ax.set_xlabel("Sample index or cluster size")
ax.set_ylabel("Ward distance")
plt.tight_layout()

# Cut the tree to get cluster labels
agg = AgglomerativeClustering(n_clusters=4, linkage="ward")
labels = agg.fit_predict(X_scaled)

import numpy as np


class PCAFromScratch:
    """PCA implemented via eigendecomposition of the covariance matrix."""

    def __init__(self, n_components: int) -> None:
        self.n_components = n_components
        self.components_: np.ndarray | None = None       # shape (n_components, d)
        self.explained_variance_: np.ndarray | None = None
        self.explained_variance_ratio_: np.ndarray | None = None
        self.mean_: np.ndarray | None = None

    def fit(self, X: np.ndarray) -> "PCAFromScratch":
        self.mean_ = X.mean(axis=0)
        X_centered = X - self.mean_

        cov = X_centered.T @ X_centered / (X.shape[0] - 1)  # unbiased estimate
        eigenvalues, eigenvectors = np.linalg.eigh(cov)

        # eigh returns ascending order; reverse for descending
        idx = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx]
        eigenvectors = eigenvectors[:, idx]

        self.components_ = eigenvectors[:, : self.n_components].T
        self.explained_variance_ = eigenvalues[: self.n_components]
        self.explained_variance_ratio_ = self.explained_variance_ / eigenvalues.sum()
        return self

    def transform(self, X: np.ndarray) -> np.ndarray:
        return (X - self.mean_) @ self.components_.T

    def fit_transform(self, X: np.ndarray) -> np.ndarray:
        return self.fit(X).transform(X)

from sklearn.decomposition import PCA
import numpy as np

pca_full = PCA(n_components=None)
pca_full.fit(X_scaled)

cumulative_variance = np.cumsum(pca_full.explained_variance_ratio_)
n_components_95 = np.searchsorted(cumulative_variance, 0.95) + 1
print(f"Components for 95% variance: {n_components_95}")

# Apply PCA
pca = PCA(n_components=n_components_95, random_state=42)
X_reduced = pca.fit_transform(X_scaled)
print(f"Original shape: {X_scaled.shape} → Reduced shape: {X_reduced.shape}")

from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, perplexity=30, n_iter=1000, learning_rate="auto",
            init="pca", random_state=42)
X_tsne = tsne.fit_transform(X_scaled)

import umap

reducer = umap.UMAP(
    n_components=2,
    n_neighbors=15,     # controls local vs global balance
    min_dist=0.1,       # controls cluster compactness
    metric="euclidean",
    random_state=42,
)
X_umap = reducer.fit_transform(X_scaled)

# UMAP supports transforming new points (t-SNE does not)
X_new_umap = reducer.transform(X_new_scaled)

import torch
import torch.nn as nn
from torch.utils.data import DataLoader, TensorDataset


class Autoencoder(nn.Module):
    """Simple fully-connected autoencoder."""

    def __init__(self, input_dim: int, latent_dim: int, hidden_dim: int = 128) -> None:
        super().__init__()
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, latent_dim),
        )
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, input_dim),
        )

    def forward(self, x: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor]:
        z = self.encoder(x)
        x_hat = self.decoder(z)
        return x_hat, z


def train_autoencoder(
    model: Autoencoder,
    X_train: torch.Tensor,
    n_epochs: int = 50,
    lr: float = 1e-3,
    batch_size: int = 256,
) -> list[float]:
    """Train autoencoder and return per-epoch reconstruction losses."""
    loader = DataLoader(TensorDataset(X_train), batch_size=batch_size, shuffle=True)
    optimiser = torch.optim.Adam(model.parameters(), lr=lr)
    criterion = nn.MSELoss()
    epoch_losses: list[float] = []

    model.train()
    for epoch in range(n_epochs):
        batch_losses: list[float] = []
        for (x_batch,) in loader:
            x_hat, _ = model(x_batch)
            loss = criterion(x_hat, x_batch)
            optimiser.zero_grad()
            loss.backward()
            optimiser.step()
            batch_losses.append(loss.item())
        epoch_loss = sum(batch_losses) / len(batch_losses)
        epoch_losses.append(epoch_loss)
        if (epoch + 1) % 10 == 0:
            print(f"Epoch {epoch+1}/{n_epochs}: loss = {epoch_loss:.6f}")

    return epoch_losses

from sklearn.ensemble import IsolationForest

iso_forest = IsolationForest(
    n_estimators=200,
    contamination=0.05,   # expected fraction of anomalies
    random_state=42,
    n_jobs=-1,
)
anomaly_labels = iso_forest.fit_predict(X_scaled)   # -1 = anomaly, 1 = normal
anomaly_scores = iso_forest.score_samples(X_scaled) # lower = more anomalous

print(f"Anomalies detected: {(anomaly_labels == -1).sum()}")

from sklearn.svm import OneClassSVM

oc_svm = OneClassSVM(kernel="rbf", nu=0.05, gamma="scale")
oc_svm.fit(X_train_normal)
labels = oc_svm.predict(X_test)  # 1 = normal, -1 = anomaly

import torch
import numpy as np


def detect_anomalies_autoencoder(
    model: Autoencoder,
    X: torch.Tensor,
    threshold_percentile: float = 95.0,
) -> tuple[np.ndarray, np.ndarray]:
    """Detect anomalies using autoencoder reconstruction error."""
    model.eval()
    with torch.no_grad():
        X_hat, _ = model(X)
        errors = ((X - X_hat) ** 2).mean(dim=1).cpu().numpy()

    threshold = np.percentile(errors, threshold_percentile)
    is_anomaly = errors > threshold
    return is_anomaly, errors