Unsupervised Learning

Classification: you have $(x_i, y_i)$ pairs. The model learns to predict $y$ from $x$. Clustering: you have only $x_i$. The model learns to assign each $x_i$ to a group, where the groups are defined by similarity rather than by external labels.

The Iris dataset shows both sides cleanly. With labels: three species, separable in petal-length × petal-width. Without labels: still three clusters, recoverable by k-means with k=3. The cluster labels k-means returns will be permutations of the true species labels, but the partition is the same up to renaming.

from sklearn.datasets import load_iris
from sklearn.cluster import KMeans
import numpy as np

X = load_iris().data
y_true = load_iris().target

km = KMeans(n_clusters=3, n_init=10, random_state=42)
y_pred = km.fit_predict(X)
# y_pred contains values {0, 1, 2} that are some permutation of y_true.

Three reasons unsupervised learning matters even when labels are available.

Exploratory analysis. You do not know what classes exist until you have looked. Clustering finds the structure before you ask "what should I label these as?"

Semi-supervised learning. Labels are expensive; raw data is cheap. Cluster the data, hand-label one example per cluster, propagate to the rest. §10 makes this concrete.

Anomaly detection. A trained clustering / density model assigns low probability to points that do not fit any cluster. Those are your anomalies.

K-means (Lloyd 1957, MacQueen 1967) clusters $m$ points into $k$ groups by alternating two steps until convergence.

1. Assignment step. Each point is assigned to its nearest centroid (Euclidean distance).
2. Update step. Each centroid is moved to the mean of the points assigned to it.

The objective is the inertia (within-cluster sum of squares):

where $z_i \in \{1, \ldots, k\}$ is the cluster assignment and $\boldsymbol{\mu}_j$ is the centroid of cluster $j$. Lloyd's algorithm is coordinate descent on this objective: fix centroids, optimize assignments; fix assignments, optimize centroids; repeat. Each step reduces (or holds) $J$, so the algorithm converges. Not necessarily to the global minimum, which is NP-hard.

km = KMeans(n_clusters=k, init="k-means++", n_init=10, random_state=42)
labels = km.fit_predict(X)
centroids = km.cluster_centers_
inertia = km.inertia_

def kmeans_pp_init(X, k, rng):
    m, d = X.shape
    centroids = np.empty((k, d))
    centroids[0] = X[rng.integers(0, m)]
    for j in range(1, k):
        diffs = X[:, None, :] - centroids[None, :j, :]
        d2 = (diffs ** 2).sum(axis=-1).min(axis=1)
        probs = d2 / d2.sum() if d2.sum() > 0 else np.ones(m) / m
        centroids[j] = X[rng.choice(m, p=probs)]
    return centroids

def kmeans(X, k, n_iters=100, seed=42, init="kmeans++", n_init=10):
    rng = np.random.default_rng(seed); m, d = X.shape; best = None
    for _ in range(n_init):
        centroids = kmeans_pp_init(X, k, rng)
        labels = np.full(m, -1, dtype=int)
        for _ in range(n_iters):
            dists = ((X[:, None, :] - centroids[None, :, :]) ** 2).sum(axis=-1)
            new_labels = dists.argmin(axis=1)
            if np.array_equal(new_labels, labels): break
            labels = new_labels
            for j in range(k):
                mask = labels == j
                if mask.any(): centroids[j] = X[mask].mean(axis=0)
        inertia = float(((X - centroids[labels]) ** 2).sum(axis=-1).sum())
        if best is None or inertia < best[2]:
            best = (centroids.copy(), labels.copy(), inertia)
    return best

Three failure modes worth knowing.

Failure 1: bad initialization. Random initialization can land all centroids near each other; the algorithm then converges to a partition where some clusters are empty or near-empty. n_init=10 (run from 10 random starts and keep the best) addresses this partially. K-means++ (§3) addresses it more cleanly.

Failure 2: spherical-cluster assumption. K-means assumes clusters are roughly spherical and equally-sized. For elongated or differently-sized clusters, k-means systematically misassigns the boundary points. The "spherical" assumption comes from the Euclidean distance + equal-weight inertia objective.

Failure 3: you have to choose k. K-means with k=5 will partition any dataset into 5 clusters, including datasets that have 2 true clusters or 10 true clusters. Diagnostic methods (§4) help but do not eliminate this problem.

Random initialization can pick all $k$ centroids from the same cluster, which means the algorithm has to migrate them apart over many iterations and often gets stuck. K-means++ (Arthur and Vassilvitskii 2007) picks them strategically.

The recipe:

1. Choose the first centroid uniformly at random from the data.
2. For each subsequent centroid: pick a data point with probability proportional to its squared distance from the nearest already-chosen centroid.

This biases the initialization toward spreading the centroids out, which dramatically improves convergence quality. The expected approximation ratio (against the optimal k-means objective) is $O(\log k)$, which is a real theoretical guarantee.

def kmeans_plusplus(X: np.ndarray, k: int, rng) -> np.ndarray:
    """k-means++ init. Returns (k, d) initial centroids."""
    m, d = X.shape
    centroids = np.empty((k, d))
    # First centroid uniformly at random
    centroids[0] = X[rng.integers(0, m)]
    for j in range(1, k):
        # Distance from each point to its nearest current centroid
        dists = ((X[:, None, :] - centroids[None, :j, :]) ** 2).sum(axis=-1).min(axis=1)
        probs = dists / dists.sum()
        idx = rng.choice(m, p=probs)
        centroids[j] = X[idx]
    return centroids

sklearn's KMeans uses k-means++ by default (init="k-means++"). The only reason to ever set init="random" is to demonstrate how bad the alternative is.

Inertia decreases monotonically as $k$ increases. With $k = m$ (one cluster per point) inertia is zero. So "minimize inertia" cannot pick $k$.

The elbow method. Plot inertia as a function of $k$ from 1 to 10. Look for the "elbow" where the curve transitions from steep to flat. Pick the $k$ at the elbow.

inertias = []
for k in range(1, 11):
    km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(X)
    inertias.append(km.inertia_)
# Plot inertias[k] for k=1..10; visually identify the elbow.

This is subjective. The elbow is sometimes sharp, often not.

Silhouette score. For each point, define $a(i)$ as its average distance to other points in its own cluster, and $b(i)$ as its average distance to points in the nearest other cluster. The silhouette for point $i$ is:

High silhouette = well-clustered. Low = ambiguous. Negative = probably misassigned. The dataset-level silhouette score is the mean over points.

from sklearn.metrics import silhouette_score

scores = []
for k in range(2, 11):
    km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(X)
    scores.append(silhouette_score(X, km.labels_))
best_k = np.argmax(scores) + 2

Pick the $k$ that maximizes silhouette. Silhouette is more objective than the elbow but slower to compute ($O(m^2)$ for the distance matrix).

Silhouette diagrams. sklearn's silhouette score is a single scalar. The diagram (per-cluster silhouette distributions, stacked) is much more informative. Géron's Figure 9-9 shows what they look like. An evenly-distributed silhouette diagram with all clusters having similar widths and high mean silhouettes is "clean clustering". Long horizontal protrusions are warning signs.

A classical k-means application: color quantization. Treat each pixel as a 3D point in (R, G, B) space, cluster into $k$ colors, replace each pixel with its cluster centroid. Result: an image with at most $k$ unique colors.

from sklearn.cluster import KMeans
from PIL import Image
import numpy as np

image = np.asarray(Image.open("ladybug.jpg")) / 255.0   # shape (H, W, 3)
H, W, _ = image.shape
pixels = image.reshape(-1, 3)

km = KMeans(n_clusters=8, n_init=10, random_state=42).fit(pixels)
segmented = km.cluster_centers_[km.labels_].reshape(H, W, 3)
Image.fromarray((segmented * 255).astype(np.uint8)).save("ladybug_k8.png")

The compressed image (8 unique colors) is recognizable but cartoonish. At $k = 64$ it is nearly indistinguishable from the original. The compression ratio is real: storing the centroid table plus the per-pixel cluster ID is much smaller than 24-bit color per pixel.

This is not how production image codecs work. JPEG and PNG use DCT and run-length encoding respectively, both of which exploit spatial correlations between pixels in ways k-means does not. K-means color quantization is a teaching example; it is not the modern compression algorithm.

K-means fails on the canonical "moons" and "concentric circles" datasets because its objective assumes spherical clusters. DBSCAN (Density-Based Spatial Clustering of Applications with Noise; Ester et al. 1996) makes a different assumption: clusters are dense regions separated by sparse regions.

The algorithm has two hyperparameters: eps (a distance threshold) and min_samples (a count). Each point falls into one of three categories.

Core point: has at least min_samples points (including itself) within distance eps.

Border point: not a core point, but within eps of one.

Noise point: neither core nor border. Gets label -1.

A cluster is a connected component of core points (plus their associated border points). Two core points are in the same cluster if one is within eps of the other; transitivity does the rest.

from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons

X, _ = make_moons(n_samples=1000, noise=0.05, random_state=42)
dbscan = DBSCAN(eps=0.05, min_samples=5)
labels = dbscan.fit_predict(X)
# labels is in {-1, 0, 1}: noise + two clusters

DBSCAN's three superpowers:

1. Finds non-spherical clusters (moons, rings, snakes).
2. Auto-determines the number of clusters (no $k$ to choose).
3. Identifies outliers natively (label -1).

DBSCAN's weaknesses:

1. Hyperparameter sensitivity. eps must be tuned. Too small: every point is noise. Too large: everything is one cluster.
2. Variable density. If your dataset has clusters of very different densities, a single eps cannot capture all of them. HDBSCAN (the hierarchical version) fixes this.
3. High-dimensional spaces. DBSCAN's distance assumptions break in high dimensions (the curse strikes again). Use it after dim reduction.

Note on eps: the standard heuristic is to plot the $k$-distance graph (distance from each point to its $k$-th nearest neighbor, sorted), look for an elbow, and use that distance as eps. sklearn.neighbors.NearestNeighbors does the heavy lifting.

K-means and DBSCAN are hard clustering: each point gets one label. Gaussian Mixture Models (GMMs) are soft clustering: each point gets a probability distribution over clusters.

The generative story: data is drawn from a mixture of $k$ Gaussian distributions, each with its own mean $\boldsymbol{\mu}_j$, covariance $\boldsymbol{\Sigma}_j$, and prior weight $\pi_j$ (with $\sum_j \pi_j = 1$). To generate a point: sample a cluster $j$ from the prior $\boldsymbol{\pi}$, then sample $\mathbf{x}$ from $\mathcal{N}(\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)$.

Here $\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) = (2\pi)^{-d/2} |\boldsymbol{\Sigma}|^{-1/2} \exp\!\big(-\tfrac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\big)$ is the multivariate Gaussian density in $d$ dimensions. The covariance matrix $\boldsymbol{\Sigma}$ is the $d \times d$ symmetric matrix whose diagonal holds each feature's variance and whose off-diagonals hold the pairwise covariances; it is what sets the size, shape, and tilt of each cluster's ellipse (a diagonal $\boldsymbol{\Sigma}$ gives axis-aligned ellipses, a full $\boldsymbol{\Sigma}$ lets them rotate). The quadratic form in the exponent, $(\mathbf{x}-\boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})$, is the squared Mahalanobis distance from $\mathbf{x}$ to the center $\boldsymbol{\mu}$ — Euclidean distance rescaled by the covariance, so it counts "how many standard deviations out" along each direction.

The likelihood:

Fitting a GMM means estimating $\{\pi_j, \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j\}_{j=1}^{k}$ that maximize the data log-likelihood.

The optimization is solved by Expectation-Maximization (EM), alternating two steps.

E step (Expectation). For each point $\mathbf{x}_i$ and cluster $j$, compute the responsibility $\gamma_{ij}$: the posterior probability that point $i$ came from cluster $j$.

This is just Bayes' rule — $\text{posterior} = \frac{\text{prior} \times \text{likelihood}}{\text{evidence}}$, i.e. $p(j \mid \mathbf{x}_i) = \pi_j \, \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j) \,/\, p(\mathbf{x}_i)$ — applied to the latent cluster assignment, with the denominator just normalizing the responsibilities to sum to one across clusters.

M step (Maximization). Update parameters using the responsibilities as soft weights.

Alternate E and M until convergence. The log-likelihood is monotonically increasing under EM (Dempster, Laird, Rubin 1977).

gmm = GaussianMixture(n_components=k, covariance_type="full", reg_covar=1e-6, random_state=42).fit(X)
labels = gmm.predict(X)          # argmax responsibility
resp = gmm.predict_proba(X)      # (m, k) responsibilities
logp = gmm.score_samples(X)      # log p(x_i)

def gaussian_pdf(X, mu, Sigma):
    d = X.shape[1]; diff = X - mu; inv = np.linalg.inv(Sigma)
    exponent = -0.5 * np.einsum("ij,jk,ik->i", diff, inv, diff)
    det = np.linalg.det(Sigma)
    norm_const = 1.0 / np.sqrt((2 * np.pi) ** d * max(det, 1e-300))
    return norm_const * np.exp(exponent)

def fit_gmm_em(X, k, n_iters=100, seed=42, reg_covar=1e-6):
    rng = np.random.default_rng(seed); m, d = X.shape
    pi = np.ones(k) / k
    mu = X[rng.choice(m, k, replace=False)].astype(np.float64)
    Sigma = np.tile(np.eye(d), (k, 1, 1))
    for it in range(n_iters):
        # E step
        weighted = np.zeros((m, k))
        for j in range(k):
            weighted[:, j] = pi[j] * gaussian_pdf(X, mu[j], Sigma[j])
        denom = weighted.sum(axis=1, keepdims=True) + 1e-300
        gamma = weighted / denom
        # M step
        Nk = gamma.sum(axis=0) + 1e-12
        pi = Nk / m
        for j in range(k):
            mu[j] = (gamma[:, j:j+1] * X).sum(axis=0) / Nk[j]
            diff = X - mu[j]
            Sigma[j] = (gamma[:, j:j+1] * diff).T @ diff / Nk[j]
            Sigma[j] = Sigma[j] + reg_covar * np.eye(d)
    return pi, mu, Sigma, gamma

GaussianMixture(n_components=k) requires you to set $k$. BayesianGaussianMixture uses a Dirichlet process prior over component weights to effectively learn $k$ from the data. Treat that prior as a black box for now: it is a prior over the mixture weights that favors using as few components as possible, so unused ones get pushed toward zero. You set an upper bound; the algorithm shrinks unused components' weights to near-zero.

from sklearn.mixture import BayesianGaussianMixture

bgmm = BayesianGaussianMixture(n_components=10, n_init=10, random_state=42).fit(X)
print(bgmm.weights_.round(3))
# e.g., [0.34, 0.33, 0.33, 0.00, 0.00, ..., 0.00]
#       Three clusters auto-detected; the rest dropped.

The "effective" number of clusters is the count of components with non-negligible weight (typically > 0.01).

bics = [GaussianMixture(n_components=k, n_init=10, random_state=42).fit(X).bic(X)
        for k in range(1, 11)]
best_k = np.argmin(bics) + 1

The BIC (Bayesian Information Criterion) is $\log(m)\,p - 2\log\hat{L}$, where $m$ is the number of points, $p$ the number of fitted parameters, and $\hat{L}$ the maximized likelihood. The first term penalizes complexity (more components means more parameters) and the second rewards fit, so minimizing BIC auto-prefers smaller $k$ when extra components do not improve the likelihood much.

Once you have a GMM, anomaly detection is a one-liner: compute the log-likelihood of each point under the fitted model, and threshold.

High score = low density = unusual point.

import numpy as np

gmm = GaussianMixture(n_components=4, n_init=10, random_state=42).fit(X_train)
log_densities = gmm.score_samples(X_test)
threshold = np.percentile(log_densities, 4)   # bottom 4% as anomalies
anomalies = X_test[log_densities < threshold]

This is a generative anomaly detector: it learns the distribution of normal data and flags points unlikely under it. Alternatives include:

• Isolation Forest: anomalies are points easy to isolate via random trees.
• One-class SVM: learn a boundary around the bulk of the data; points outside are anomalies.
• Autoencoders: anomalies have high reconstruction error.

Each has its niche. GMM-based anomaly detection is interpretable (you can ask "which component is this point closest to?") and fast at inference. It struggles in high dimensions for the usual reasons.

You have 1000 unlabeled examples and 50 labeled ones. A fully-supervised classifier on 50 examples will overfit. A clustering of all 1050 examples, with the labels propagated through cluster membership, can dramatically beat that.

The recipe:

1. Cluster all examples (labeled + unlabeled) with k-means or GMM, choosing $k$ slightly larger than the number of classes.
2. For each cluster, look at the labels of the labeled examples that fell into it. Assign the cluster's label as the majority.
3. Propagate that label to every unlabeled point in the cluster.

from sklearn.cluster import KMeans
import numpy as np

km = KMeans(n_clusters=50, n_init=10, random_state=42).fit(X_train)
# For each cluster, find the labeled example closest to its centroid
representative_idx = np.argmin(((X_train[labeled_idx] - km.cluster_centers_[:, None, :]) ** 2)
                                .sum(axis=-1), axis=1)
# Use those labels as cluster labels, propagate to all unlabeled points
cluster_labels = y_train[labeled_idx][representative_idx]
y_propagated = cluster_labels[km.labels_]

Géron shows this on MNIST: training a logistic regression on 50 manually-labeled examples + the propagated labels of the other 950 gives ~90% accuracy, compared to ~75% for logistic regression on the 50 labeled alone.

GMMs are generative models: they model $p(\mathbf{x})$, the distribution of data. Classifiers like logistic regression are discriminative: they model $p(y \mid \mathbf{x})$, the conditional distribution given features.

Generative models are more flexible. You can sample new data, detect anomalies, do imputation of missing values, and combine with prior information via Bayes' rule. They are also harder to train and often less accurate when the only goal is classification.

The 2026 frontier of generative modeling is not GMMs. It is diffusion models, variational autoencoders, and GANs (another chapter). The connection back to this chapter: GANs are trained against a discriminator, but their objective can be reinterpreted as fitting a density, which is exactly what a GMM does. The continuum from "fit a mixture of Gaussians" to "train a deep generator network" is the lineage of modern generative AI. Lilian Weng's GAN post (18-lilian-weng/2017-08-20-gan) draws this line explicitly.

For tabular and small-scale data, GMMs are still a serious tool. For images, text, audio, you need the deep stack.

A decision tree of decisions, learned by years of practitioners falling into the same traps.

Do you know k?

• Yes, your application has a natural cluster count → k-means.
• No, you have to discover it → BayesianGaussianMixture, DBSCAN, or BIC-tuned GMM.

Are clusters spherical (roughly equal-density, blob-shaped)?

• Yes → k-means.
• No, but they are smooth and density-defined → DBSCAN.
• No, they are elongated or oddly-shaped → GMM with full covariance, or hierarchical clustering.

Do you need soft assignments / probabilities?

• Yes → GMM (provides predict_proba).
• No → k-means or DBSCAN.

Do you have outliers you want to identify?

• Yes → DBSCAN (gives noise label) or GMM-based anomaly detection.
• No → k-means.

Is your data high-dimensional?

• Yes → run PCA first, then cluster on the reduced output.
• No → cluster directly.

Is your dataset huge (m > 100k)?

• Yes → MiniBatchKMeans, or DBSCAN with appropriate indexes.
• No → standard sklearn fits.

Use 2D PCA or UMAP projections as exploratory views, not as ground truth. They can suggest candidate structure and outliers, but every projection discards or distorts information; validate any choice with metrics and diagnostics in the original feature space. Five minutes of visual inspection can narrow the search, but it cannot identify the right clustering algorithm by itself.