Unsupervised Learning: Structure Without Labels

Unsupervised learning asks a fundamentally harder question than supervised learning: given data X with no labels, what structure exists? The difficulty is that “structure” has no ground truth — you can’t compute a test accuracy against reality. This is why unsupervised methods require more mathematical discipline, not less. You need to understand exactly what objective each algorithm is optimizing and what assumptions it bakes in, because nothing will tell you if you’re wrong.

This file covers five methods in depth: k-means, hierarchical clustering, PCA (revisited from the ML application angle), t-SNE/UMAP (with honest caveats about what they lie about), and autoencoders (the bridge to deep learning). Each connects directly to the Phase 0 mathematical machinery.

---

1. K-Means Clustering

What K-Means Is Actually Optimizing

K-means minimizes the within-cluster sum of squares (WCSS):

J = Σₖ Σᵢ∈Cₖ ‖xᵢ - μₖ‖²

Where Cₖ is cluster k and μₖ is its centroid. This is a non-convex optimization problem — finding the global minimum is NP-hard. K-means finds a local minimum via coordinate descent:

E-step: Assign each point to its nearest centroid
M-step: Update each centroid to the mean of its assigned points
Repeat until convergence (assignments stop changing)

This is the Expectation-Maximization (EM) structure. E-step = soft assignment to “expectation”; M-step = “maximization” of likelihood. K-means is a hard-assignment EM algorithm. Gaussian Mixture Models (GMMs) are the soft-assignment version — each point has a probability of belonging to each cluster. Understanding this connection is essential for Phase 3 (generative models).

Convergence Proof Intuition

Each step decreases J:

• E-step: reassigning each point to its nearest centroid can only decrease J (it’s the definition of nearest)

• M-step: updating centroid to the mean minimizes the sum of squared distances to that centroid (the mean is the unique minimizer of squared Euclidean distance)

Therefore J is non-increasing at each step. Since there are finitely many possible assignments (K^n partitions), the algorithm must terminate. But it may terminate at a local minimum — hence multiple random restarts (n_init=10 in sklearn).

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs

class KMeansScratch:
    """K-means via E-M coordinate descent."""
    
    def __init__(self, k, n_init=10, max_iters=300, tol=1e-4):
        self.k = k
        self.n_init = n_init
        self.max_iters = max_iters
        self.tol = tol
    
    def _run_once(self, X):
        """Single run with random initialization."""
        n, d = X.shape
        # Random initialization: pick k data points as initial centroids
        idx = np.random.choice(n, self.k, replace=False)
        centroids = X[idx].copy()
        
        for iteration in range(self.max_iters):
            # E-step: assign to nearest centroid
            dists = np.linalg.norm(X[:, None, :] - centroids[None, :, :], axis=2)
            labels = np.argmin(dists, axis=1)
            
            # M-step: recompute centroids
            new_centroids = np.array([
                X[labels == k].mean(axis=0) if (labels == k).sum() > 0 else centroids[k]
                for k in range(self.k)
            ])
            
            # Convergence check
            if np.linalg.norm(new_centroids - centroids) < self.tol:
                break
            
            centroids = new_centroids
        
        # Compute WCSS for this run
        wcss = sum(
            np.sum((X[labels == k] - centroids[k])**2)
            for k in range(self.k)
        )
        return centroids, labels, wcss
    
    def fit(self, X):
        """Multiple restarts, keep best (lowest WCSS)."""
        best_wcss = np.inf
        for _ in range(self.n_init):
            centroids, labels, wcss = self._run_once(X)
            if wcss < best_wcss:
                best_wcss = wcss
                self.centroids_ = centroids
                self.labels_ = labels
                self.inertia_ = wcss
        return self

# The Elbow Method: how to choose K
X, y_true = make_blobs(n_samples=500, centers=4, cluster_std=0.8, random_state=42)

wcss_values = []
k_range = range(1, 11)
for k in k_range:
    km = KMeansScratch(k=k, n_init=5)
    km.fit(X)
    wcss_values.append(km.inertia_)

# The "elbow" is where marginal WCSS decrease starts to slow — heuristic for K
# Note: elbow method is a heuristic. For rigorous K selection, use silhouette score.
print("WCSS by K:", [f"K={k}: {w:.1f}" for k, w in zip(k_range, wcss_values)])

# Verify against sklearn
from sklearn.cluster import KMeans
km_sk = KMeans(n_clusters=4, random_state=42)
km_sk.fit(X)
print(f"\nScratch inertia: {KMeansScratch(k=4).fit(X).inertia_:.1f}")
print(f"sklearn inertia: {km_sk.inertia_:.1f}")

K-Means Failure Modes

Problem	What Goes Wrong	Fix
Unequal cluster sizes	Large clusters dominate WCSS	Use GMMs or density-based clustering
Non-spherical clusters	K-means assumes Voronoi/spherical structure	Use DBSCAN or spectral clustering
Different densities	Low-density cluster absorbs high-density cluster	Normalize or use GMM
Local minima	Different runs give different results	n_init=10+ with k-means++ initialization
Wrong K	Forced K=5 when true K=3	Elbow + silhouette + domain knowledge

---

2. Hierarchical Clustering

Hierarchical clustering builds a dendrogram — a tree of nested clusters — without specifying K in advance. Cut the dendrogram at any height to get any number of clusters.

Agglomerative (bottom-up):

1. Start: each point is its own cluster

2. Merge the two closest clusters

3. Repeat until one cluster remains

4. Cut the resulting tree at the desired level

Distance between clusters (linkage criterion) is the key choice:

Linkage	Cluster distance definition	Tendency
Single	Min pairwise distance	Chaining (elongated clusters)
Complete	Max pairwise distance	Compact, roughly equal clusters
Average	Mean pairwise distance	Compromise
Ward	Minimizes WCSS increase on merge	Most similar to k-means; often best

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

X, _ = make_blobs(n_samples=50, centers=3, random_state=42)

# Scipy for dendrogram visualization
Z = linkage(X, method='ward')

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
dendrogram(Z, truncate_mode='level', p=5)
plt.title('Dendrogram (Ward linkage)')
plt.xlabel('Sample index')
plt.ylabel('Distance')

# Cut at K=3
plt.subplot(1, 2, 2)
agg = AgglomerativeClustering(n_clusters=3, linkage='ward')
labels = agg.fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis')
plt.title('Agglomerative Clustering (K=3)')
plt.tight_layout()
plt.savefig('hierarchical_clustering.png', dpi=150)

# Complexity: O(n² log n) time, O(n²) space — does not scale to large datasets
print(f"Unique clusters found: {len(set(labels))}")

When to use: When you need to explore cluster hierarchy, don’t know K in advance, or have small-to-medium data (<10K samples). Ward linkage + agglomerative is usually the default starting point.

---

3. PCA — Revisited from the ML Application Angle

PCA was covered mathematically in Phase 0. Here we focus on the applied ML usage patterns: when it helps, when it hurts, and what you should actually do in production pipelines.

The Core Mechanism (Brief Recap)

PCA finds directions of maximum variance. The k-th principal component is the k-th eigenvector of the covariance matrix Σ = (1/n)XᵀX, sorted by eigenvalue magnitude.

from sklearn.decomposition import PCA
from sklearn.datasets import load_digits
import matplotlib.pyplot as plt
import numpy as np

# Digits dataset: 64 features (8×8 pixel values)
X, y = load_digits(return_X_y=True)
print(f"Original shape: {X.shape}")  # (1797, 64)

# Fit PCA
pca = PCA()
pca.fit(X)

# Explained variance: how much information each component captures
cumvar = np.cumsum(pca.explained_variance_ratio_)
n_components_95 = np.argmax(cumvar >= 0.95) + 1
print(f"Components for 95% variance: {n_components_95}")  # Typically ~29 for digits

# Reconstruction quality at different compression levels
for n_comp in [2, 10, 30, 64]:
    pca_k = PCA(n_components=n_comp)
    X_reduced = pca_k.fit_transform(X)
    X_reconstructed = pca_k.inverse_transform(X_reduced)
    
    mse = np.mean((X - X_reconstructed)**2)
    compression = n_comp / 64
    print(f"n_components={n_comp:3d}: reconstruction MSE={mse:.2f}, compression={compression:.1%}")

# Visualization in 2D — the standard diagnostic
pca2 = PCA(n_components=2)
X_2d = pca2.fit_transform(X)

plt.figure(figsize=(10, 8))
scatter = plt.scatter(X_2d[:, 0], X_2d[:, 1], c=y, cmap='tab10', alpha=0.6)
plt.colorbar(scatter)
plt.title('Digits dataset: PCA 2D projection')
plt.xlabel(f'PC1 ({pca2.explained_variance_ratio_[0]:.1%} variance)')
plt.ylabel(f'PC2 ({pca2.explained_variance_ratio_[1]:.1%} variance)')
plt.savefig('pca_digits_2d.png', dpi=150)

When PCA Helps vs. Hurts

Scenario	PCA Verdict
High-dimensional data with correlated features (e.g., pixels)	✅ Strong dimensionality reduction, speeds up downstream models
Tabular data with 10-50 features	⚠️ Often hurts — destroys interpretability, rarely helps tree models
Preprocessing for k-means	✅ Removes noise dimensions that distort Euclidean distance
Preprocessing before neural network	⚠️ Usually unnecessary — network learns its own projections
Data visualization (2D/3D)	✅ Always useful as a diagnostic tool

What most people get wrong: applying PCA as a default preprocessing step for tabular ML competitions. Tree-based models (XGBoost, Random Forest) are invariant to monotone feature transformations and rotation — PCA doesn’t help them. Apply PCA when you have evidence of redundant correlation structure, not by default.

---

4. t-SNE and UMAP — What They Do and What They Lie About

t-SNE and UMAP are visualization techniques, not dimensionality reduction techniques for machine learning pipelines. This distinction is critical and routinely misunderstood.

t-SNE: The Mechanism

t-SNE (t-distributed Stochastic Neighbor Embedding) preserves local structure at the cost of global structure:

1. Compute pairwise similarities in high-d space: P(i|j) proportional to exp(-‖xᵢ-xⱼ‖²/2σ²) (Gaussian)

2. Define a target distribution in 2D: Q proportional to (1 + ‖yᵢ-yⱼ‖²)⁻¹ (Student-t, heavier tails)

3. Minimize KL(P ‖ Q) via gradient descent

The Student-t distribution in low-d space is the key insight: it pushes dissimilar points far apart (t has heavier tails than Gaussian, so distant points need to be very far in 2D to match low P). This creates the characteristic “cluster” appearance.

from sklearn.manifold import TSNE
from sklearn.datasets import load_digits
import matplotlib.pyplot as plt
import time

X, y = load_digits(return_X_y=True)

# t-SNE is expensive: O(n² log n) with Barnes-Hut approximation
# Reduce to PCA first for speed (standard practice)
from sklearn.decomposition import PCA
X_pca50 = PCA(n_components=50).fit_transform(X)

start = time.time()
tsne = TSNE(n_components=2, perplexity=30, n_iter=1000, random_state=42)
X_tsne = tsne.fit_transform(X_pca50)
print(f"t-SNE time: {time.time() - start:.1f}s")

plt.figure(figsize=(10, 8))
scatter = plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=y, cmap='tab10', alpha=0.6)
plt.colorbar(scatter)
plt.title('Digits: t-SNE (perplexity=30)')
plt.savefig('tsne_digits.png', dpi=150)

What t-SNE Lies About (Critical Warnings)

⚠️ DO NOT trust t-SNE for these inferences:

1. Cluster sizes are meaningless: t-SNE distorts distances to pack points into 2D. A large cluster is NOT a more common class.

2. Distances between clusters are meaningless: The space between cluster A and cluster B in a t-SNE plot tells you nothing about their true similarity.

3. Perplexity changes the plot significantly: perplexity=5 and perplexity=50 can look completely different on the same data. There is no “correct” perplexity.

4. Random seed dependence: Two runs with different seeds produce different-looking (but equally valid) plots.

5. Cannot be applied to new points: t-SNE has no transform() method — it must be refit on the entire dataset. → Use UMAP for production visualization.

UMAP: The Better Production Tool

UMAP (Uniform Manifold Approximation and Projection) is faster, preserves more global structure, and has a transform() for new points.

# pip install umap-learn
import umap

reducer = umap.UMAP(n_components=2, n_neighbors=15, min_dist=0.1, random_state=42)
X_umap = reducer.fit_transform(X)

# Unlike t-SNE, you can embed NEW points:
# X_new_embedded = reducer.transform(X_new)

plt.figure(figsize=(10, 8))
scatter = plt.scatter(X_umap[:, 0], X_umap[:, 1], c=y, cmap='tab10', alpha=0.6)
plt.colorbar(scatter)
plt.title('Digits: UMAP (n_neighbors=15, min_dist=0.1)')
plt.savefig('umap_digits.png', dpi=150)

Practical rule: Use t-SNE for exploratory analysis in notebooks. Use UMAP when you need a persistent embedding or need to embed new data points.

---

5. Autoencoders — The Bridge to Deep Learning

An autoencoder is a neural network trained to reproduce its input. It has two parts:

• Encoder: f: X → Z (compresses input to latent code Z)

• Decoder: g: Z → X̂ (reconstructs from latent code)

Training objective: minimize ‖X - g(f(X))‖² (reconstruction error)

The bottleneck dimension of Z forces the network to learn a compressed representation. When Z is linear and the reconstruction loss is MSE, an autoencoder is equivalent to PCA — it learns the same subspace. When Z is nonlinear (hidden layers with activations), it learns a nonlinear manifold — strictly more expressive than PCA.

import torch
import torch.nn as nn
import torch.optim as optim
from sklearn.datasets import load_digits
from sklearn.preprocessing import MinMaxScaler
import numpy as np

# Data prep
X, y = load_digits(return_X_y=True)
X_scaled = MinMaxScaler().fit_transform(X)  # [0, 1]
X_tensor = torch.FloatTensor(X_scaled)

class Autoencoder(nn.Module):
    def __init__(self, input_dim=64, latent_dim=8):
        super().__init__()
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, 32),
            nn.ReLU(),
            nn.Linear(32, latent_dim),
        )
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, 32),
            nn.ReLU(),
            nn.Linear(32, input_dim),
            nn.Sigmoid()  # output in [0,1]
        )
    
    def forward(self, x):
        z = self.encoder(x)
        return self.decoder(z), z

model = Autoencoder(input_dim=64, latent_dim=8)
optimizer = optim.Adam(model.parameters(), lr=1e-3)
criterion = nn.MSELoss()

# Training loop
losses = []
for epoch in range(200):
    optimizer.zero_grad()
    x_hat, z = model(X_tensor)
    loss = criterion(x_hat, X_tensor)
    loss.backward()
    optimizer.step()
    if epoch % 20 == 0:
        losses.append(loss.item())
        print(f"Epoch {epoch}: loss={loss.item():.4f}")

# Extract latent representations
model.eval()
with torch.no_grad():
    _, Z = model(X_tensor)
    Z_np = Z.numpy()

# How much structure does the 8D latent space preserve?
# Measure: KNN classification accuracy in latent space vs. original space
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score

knn_original = KNeighborsClassifier(n_neighbors=5)
knn_latent = KNeighborsClassifier(n_neighbors=5)

acc_original = cross_val_score(knn_original, X_scaled, y, cv=5).mean()
acc_latent = cross_val_score(knn_latent, Z_np, y, cv=5).mean()

print(f"\nKNN accuracy in original space (64D): {acc_original:.4f}")
print(f"KNN accuracy in latent space (8D):   {acc_latent:.4f}")
# Good autoencoder should retain most structure

Variants That Bridge to Deep Learning

Variant	Change	What It Unlocks
Denoising AE	Train to reconstruct clean X from noisy input	Robust features; precursor to diffusion models
Sparse AE	L1 penalty on Z	Dictionary learning; interpretable features
Variational AE (VAE)	Z is a distribution N(μ, σ²); KL regularization	Generative modeling; latent space interpolation
Contractive AE	Penalize Jacobian of encoder	Smooth, robust manifold representation

The VAE is Phase 3 material. When you get there, this autoencoder implementation is the direct precursor — the only additions are the reparameterization trick and the KL term in the loss.

---

Practice Problems

1. K-means vs. GMM: Generate data from two overlapping Gaussian distributions. Apply k-means and sklearn’s GaussianMixture. Compare cluster assignments for the overlapping region. Explain why GMM is “softer” in terms of the probability distribution it assumes.

2. t-SNE perplexity sweep: On the digits dataset, run t-SNE with perplexity ∈ {5, 30, 100}. Plot all three. Write 3 sentences quantifying how the cluster structure changes. Confirm that none of these is “correct” — they are all valid views of different neighborhood scales.

3. Autoencoder compression benchmark: Train autoencoders with latent_dim ∈ {2, 4, 8, 16, 32} on digits. For each, record reconstruction MSE and 5-fold KNN accuracy in latent space. Plot both curves vs. latent_dim. Identify the latent_dim where you get ≥90% of original KNN accuracy with the smallest bottleneck.

4. PCA vs. Autoencoder: On the digits dataset, apply PCA to k components and train an autoencoder with latent_dim=k for k ∈ {2, 8, 32}. Compare reconstruction MSE. For k=2, plot both 2D projections colored by digit class. Quantify when the nonlinear autoencoder beats linear PCA.

---

Return to README.md · Next: 03_model_selection_and_evaluation.md