Unsupervised Learning

Unsupervised learning fits a model to inputs that come with no labels: a collection of images, patches, pixel colours or descriptors, and nothing else. Instead of learning a mapping from an image to a class, it learns something about the images themselves, such as how they are distributed, which of them belong together, or which few numbers describe them. Methods that work without labels supply visual vocabularies, colour models, background models, compressed descriptors and pretrained features.

Definition

In unsupervised learning the training set is $\mathcal{D} = \{\mathbf{x}_n : n = 1, \dots, N\}$, inputs without the targets $y_n$ of supervised learning. Murphy (2022, §1.3) puts it in probabilistic terms: supervised learning fits a conditional model $p(y \mid \mathbf{x})$, while unsupervised learning fits an unconditional model $p(\mathbf{x})$, one that can also generate new data. He notes (§1.3, footnote 6; §4.2.1, footnote 2) that statisticians write the modelled variable as $y$, so for them unsupervised learning is the special case of fitting $p(y)$ with no inputs; machine learning switches to $\mathbf{x}$.

Torralba, Isola and Freeman (2024, §9.3) give a definition that does not need a density: the learner is given only input data and must come up with a model or representation of it that has useful properties, as measured by some objective function. Compressing the data into fewer dimensions while keeping its information is their example of such an objective.

Intuition

Take a large set of $8 \times 8$ patches cut from photographs. No one has said what any patch shows, yet they are far from random: neighbouring pixels have similar values, edges and smooth ramps recur, and patches that look like uniformly random static are, in our reading, vanishingly rare. An unsupervised method captures this regularity. A density model learns which patches are likely, a clustering groups similar patches into a vocabulary, and a dimensionality reduction finds a few coordinates, such as edge orientation and contrast, along which the patches vary.

Murphy (§1.3) gives three reasons to learn this way: it avoids collecting large labelled datasets, which can be slow and expensive (he cites doctors labelling medical images); it avoids forcing the world into categories that people themselves disagree on, such as where an action in a video starts; and it makes the model explain the high-dimensional input rather than a low-dimensional label.

Formal Definition

A density model $q_\theta(\mathbf{x})$ can be fitted by maximum likelihood, which Murphy (§4.2.2) justifies as minimizing the Kullback–Leibler divergence $D_{\mathrm{KL}}(p_{\mathcal{D}} \,\|\, q_\theta)$ from the empirical distribution $p_{\mathcal{D}}$ of the training images to the model; the generative model article gives the identity and cross-entropy loss owns its terms. Torralba and colleagues (§32.4.1, eq. 32.1) derive it from the true distribution $p_{\text{data}}$ and note that the reverse divergence, $\mathrm{KL}(q_\theta \,\|\, p_{\text{data}})$, usually cannot be computed, because $p_{\text{data}}$ itself is unknown.

The training log-likelihood cannot judge the model, because it can keep rising as the model grows: a model that memorizes its training images assigns them high probability and has learned nothing general, the generative counterpart of overfitting a classifier (Torralba et al., §32.3.1). A common yardstick is instead the negative log-likelihood of images not used in fitting (Murphy, §1.3.4, eq. 1.33):

$$\mathcal{L}(\theta; \mathcal{D}_{\text{test}}) = -\frac{1}{|\mathcal{D}_{\text{test}}|} \sum_{\mathbf{x} \in \mathcal{D}_{\text{test}}} \log q_\theta(\mathbf{x}).$$

For a continuous image, $q_\theta$ is a density over intensity values, so the log-likelihood depends on the units of the pixels and can be positive; two models are compared only on the same data in the same units.

The Tasks

The scikit-learn User Guide (§2) lists mixture models, manifold learning, clustering, biclustering, matrix factorization, covariance estimation, novelty and outlier detection, density estimation and unsupervised neural networks under unsupervised learning. For vision the tasks group as follows, and each has its own article.

Task	What it returns	Articles
Clustering	a partition, soft memberships or a hierarchy	clustering, k-means, mean shift, spectral clustering
Density estimation	a function $q(\mathbf{x})$ that can be evaluated	Gaussian mixture model, expectation–maximization, normalizing flow
Dimensionality reduction	a few coordinates per item	dimensionality reduction, principal component analysis, kernel methods, t-SNE, autoencoder
Generative modelling	a sampler, with or without a density	generative model, variational autoencoder, generative adversarial network, diffusion model
Representation learning	an encoder whose features serve later tasks	representation learning, self-supervised learning, contrastive learning

The rows overlap. A Gaussian mixture is a density and a soft clustering; probabilistic PCA is a dimensionality reduction and a density, as the PCA article describes. Torralba and colleagues (§30.3.3, Table 30.1) file clustering under representation learning as compression: it quantizes continuous data into discrete categories, as autoencoders compress by removing redundant information. Murphy (§1.3.3) presents self-supervised learning, which builds proxy supervised tasks from unlabelled data, as an approach to unsupervised learning; semi-supervised learning adds a few labels to many unlabelled images.

Evaluating Without Labels

Evaluation is hard, because there is no ground truth to compare with (Murphy, §1.3.4). Murphy describes three routes.

• Held-out likelihood. A good density model is not surprised by new data. Because the density must integrate to one, putting mass where the data are takes it away from where they are not. Murphy adds two caveats: density estimation is difficult in high dimensions, and a model that assigns high probability to the data need not have learned useful high-level patterns. The generative model article shows a related gap: likelihood and sample quality can disagree.
• Downstream performance. The learned representation feeds a supervised method, and a good representation reaches the same accuracy with fewer labelled examples, a gain in sample efficiency. Representation learning compares encoders this way.
• Understanding. In science the aim is often an interpretable model that explains the observed patterns; checking that it found the real structure often needs experiments.

Clusterings have their own measures: external indices such as the adjusted Rand index against a reference partition, and internal indices such as the silhouette, which the clustering article defines. Choosing a model’s size, the number of components or of dimensions, is part of the same problem; held-out likelihood is one answer, and BIC and AIC another, as the Gaussian mixture model article explains.

Kernel Density Estimation

A density estimate that makes no parametric assumption is the kernel density estimate, or Parzen window: a kernel of width $h$ on every training sample, summed and normalized. The scikit-learn User Guide (§2.8) introduces it as a histogram whose blocks are centred on the data instead of on a grid, notes that the bandwidth trades bias (large $h$, smooth) against variance (small $h$, rough), and warns that the curse of dimensionality degrades it in high dimensions. Murphy (§16.3) calls it a generative model that needs no fitting beyond choosing $h$, at the cost of storing all the data and of slow evaluation. Mean shift owns the estimator’s formula and climbs its gradient to the modes. The GAN paper article’s worked example owns the measurement of its bias: a Parzen window fitted to samples from a known density of $8 \times 8$ images, scored on held-out images, falls far below the true log-likelihood. In scikit-learn 1.7.2, KernelDensity.score returns the total log-likelihood of its input, while PCA.score and GaussianMixture.score return the mean per sample (installed source), so a comparison divides the first by the number of samples.

Examples

One yardstick for three model families

The script builds $8 \times 8$ patches whose density is known exactly: four basis images (vertical stripes, a diagonal grating and two intensity ramps, orthonormalized) with independent Gaussian coefficients of standard deviations 0.8, 0.6, 0.4 and 0.2, plus white noise $\sigma = 0.05$, around a mean of 0.5. The patches are therefore Gaussian with covariance $B\,\mathrm{diag}(s^2)\,B^\top + \sigma^2 I$, which has four large eigenvalues and 60 equal to $\sigma^2$. Three families are fitted to 2,000 patches and scored on 5,000 others: probabilistic PCA with $q = 1, \dots, 10$ components (scikit-learn’s PCA.score_samples), full-covariance Gaussian mixtures with $k = 1, \dots, 4$ components, and a Gaussian with independent pixels. Each model is compared with the exact density on the same test patches, with the standard error of that paired difference. The parameter counts are ours. scikit-learn’s PCA divides the sample covariance by $n - 1$, and GaussianMixture by $n$ before adding $10^{-6}$ to the diagonal (installed source); the hand cross-checks follow each. Log-likelihoods are in nats per patch, for intensities in the units of the construction. Seed 0 prints the full table and a summary covers data seeds 0–4 (NumPy 2.3.5, SciPy 1.16.3, scikit-learn 1.7.2; about 3 s; two runs give identical output).

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v] = "1"                      # one BLAS thread, set before NumPy is imported
import numpy as np
from scipy.stats import multivariate_normal
from sklearn.decomposition import PCA
from sklearn.mixture import GaussianMixture

# 8x8 patches (d = 64) built from q0 = 4 known basis images plus white noise, so the
# patch density is exactly Gaussian: N(0.5, B diag(s^2) B^T + sigma^2 I).
H = W = 8; D = H * W; Q0 = 4; SIGMA = 0.05
S_STD = np.array([0.8, 0.6, 0.4, 0.2])      # coefficient standard deviations
N_TRAIN, N_TEST = 2000, 5000

yy, xx = np.mgrid[0:H, 0:W].astype(float)
raw = np.stack([np.cos(2 * np.pi * 2 * xx / W),              # vertical stripes, 2 cycles per patch
                np.cos(2 * np.pi * (xx + yy) / W),            # diagonal grating
                xx - (W - 1) / 2,                              # horizontal ramp
                yy - (H - 1) / 2]).reshape(Q0, D).T            # vertical ramp
B, _ = np.linalg.qr(raw)                                       # orthonormalized, in this order
MU = np.full(D, 0.5)
COV = B @ np.diag(S_STD**2) @ B.T + SIGMA**2 * np.eye(D)

def make_data(seed):
    rng = np.random.default_rng(seed)
    def draw(n):
        a = rng.standard_normal((n, Q0)) * S_STD
        return MU + a @ B.T + SIGMA * rng.standard_normal((n, D))
    return draw(N_TRAIN), draw(N_TEST)

def gauss_ll(X, mu, C):
    """Per-patch log-density (nats) of N(mu, C)."""
    R = X - mu
    _, logdet = np.linalg.slogdet(C)
    return -0.5 * (np.einsum("ij,ij->i", R @ np.linalg.inv(C), R) + logdet + D * np.log(2 * np.pi))

def ppca_cov(Xtr, q):
    """Tipping-Bishop covariance from our own eigendecomposition (n - 1 denominator, as sklearn)."""
    lam, U = np.linalg.eigh(np.cov(Xtr, rowvar=False))
    lam, U = lam[::-1], U[:, ::-1]
    s2 = lam[q:].mean()
    return U[:, :q] @ np.diag(lam[:q] - s2) @ U[:, :q].T + s2 * np.eye(D)

def n_params(name, q):
    if name == "PPCA": return D + D * q - q * (q - 1) // 2 + 1
    if name == "GMM full": return (q - 1) + q * (D + D * (D + 1) // 2)
    return 2 * D                               # independent pixels: a mean and a variance each

print(f"C2 basis orthonormal: max |B^T B - I| = {np.abs(B.T @ B - np.eye(Q0)).max():.1e}")
summary = []
for seed in range(5):
    Xtr, Xte = make_data(seed)
    exact = gauss_ll(Xte, MU, COV)
    c3 = np.abs(exact - multivariate_normal(MU, COV).logpdf(Xte)).max()
    ev = np.sort(np.linalg.eigvalsh(np.cov(Xtr, rowvar=False)))[::-1]
    c1 = int((ev > 2 * SIGMA**2).sum())

    rows, c4, c5, conv = [], 0.0, 0.0, []
    for q in range(1, 11):
        m = PCA(n_components=q).fit(Xtr)
        ll = m.score_samples(Xte)
        c4 = max(c4, abs(m.score(Xte) - gauss_ll(Xte, Xtr.mean(0), ppca_cov(Xtr, q)).mean()))
        rows.append(("PPCA", q, m.score(Xtr), ll))
    for k in range(1, 5):
        m = GaussianMixture(n_components=k, covariance_type="full", random_state=0).fit(Xtr)
        conv.append(bool(m.converged_))
        rows.append(("GMM full", k, m.score(Xtr), m.score_samples(Xte)))
    m = GaussianMixture(n_components=1, covariance_type="diag", random_state=0).fit(Xtr)
    conv.append(bool(m.converged_))
    rows.append(("indep. pixels", 1, m.score(Xtr), m.score_samples(Xte)))

    mu_hat = Xtr.mean(0)
    full_hand = gauss_ll(Xte, mu_hat, np.cov(Xtr, rowvar=False, bias=True) + 1e-6 * np.eye(D))
    diag_hand = gauss_ll(Xte, mu_hat, np.diag(Xtr.var(0) + 1e-6))
    c5 = max(abs(rows[10][3].mean() - full_hand.mean()), abs(rows[14][3].mean() - diag_hand.mean()))

    def paired(a, b):                          # mean of a - b per patch, and its standard error
        d = a - b
        return d.mean(), d.std(ddof=1) / np.sqrt(len(d))

    if seed == 0:
        print(f"\nseed 0: exact density, held-out mean log-likelihood {exact.mean():.4f} nats per patch "
              f"(SE {exact.std(ddof=1) / np.sqrt(N_TEST):.4f})")
        print(f"C1 training-covariance eigenvalues 1-6: {np.round(ev[:6], 5)}; above 2 sigma^2: {c1}")
        print("model           q/k  params  train LL  held-out LL  minus exact  SE(diff)")
        for name, q, tr, ll in rows:
            dm, ds = paired(ll, exact)
            print(f"{name:14s} {q:3d} {n_params(name, q):7d} {tr:9.4f} {ll.mean():12.4f} {dm:12.4f} {ds:9.4f}")
        print(f"C3 max |hand - scipy| exact log-density: {c3:.1e}; C4 max |hand PPCA - PCA.score|: {c4:.1e}; "
              f"C5 max |hand Gaussian - GaussianMixture.score|: {c5:.1e}")

    ppca = [r for r in rows if r[0] == "PPCA"]; gmm = [r for r in rows if r[0] == "GMM full"]
    best_q = ppca[int(np.argmax([r[3].mean() for r in ppca]))][1]
    best_k = gmm[int(np.argmax([r[3].mean() for r in gmm]))][1]
    z_max = max(paired(r[3], exact)[0] / paired(r[3], exact)[1] for r in rows)
    d3, s3 = paired(ppca[Q0 - 1][3], rows[14][3])          # PPCA q0 minus independent pixels
    d5, s5 = paired(ppca[Q0 - 1][3], gmm[0][3])            # PPCA q0 minus full Gaussian
    summary.append((seed, c1, best_q, best_k, z_max, d3, s3, d5, s5, all(conv), c3, c4, c5))

print("\nseed  C1  best q  best k  max (model - exact)/SE  PPCA4 - indep (SE)  PPCA4 - full (SE)  "
      "all GMM converged  C3, C4, C5")
for s, c1, bq, bk, z, d3, s3, d5, s5, cv, c3, c4, c5 in summary:
    print(f"{s:4d} {c1:3d} {bq:7d} {bk:7d} {z:22.4f}   {d3:8.4f} ({s3:.4f})  {d5:7.4f} ({s5:.4f})  "
          f"{str(cv):>17s}  {c3:.1e}, {c4:.1e}, {c5:.1e}")

Output:

C2 basis orthonormal: max |B^T B - I| = 6.7e-16

seed 0: exact density, held-out mean log-likelihood 92.0757 nats per patch (SE 0.0801)
C1 training-covariance eigenvalues 1-6: [0.6493  0.34283 0.16673 0.04347 0.00341 0.00339]; above 2 sigma^2: 4
model           q/k  params  train LL  held-out LL  minus exact  SE(diff)
PPCA             1     129   51.0185      50.0266     -42.0490    0.3590
PPCA             2     192   69.5326      69.8688     -22.2069    0.2773
PPCA             3     254   86.3009      86.2200      -5.8557    0.1329
PPCA             4     315   92.1216      91.9981      -0.0775    0.0053
PPCA             5     375   92.1489      91.9799      -0.0957    0.0060
PPCA             6     434   92.1759      91.9569      -0.1188    0.0066
PPCA             7     492   92.1996      91.9337      -0.1419    0.0072
PPCA             8     549   92.2198      91.9145      -0.1611    0.0076
PPCA             9     605   92.2402      91.8997      -0.1759    0.0081
PPCA            10     660   92.2592      91.8816      -0.1941    0.0084
GMM full         1    2144   92.5987      91.5114      -0.5643    0.0156
GMM full         2    4289   93.4251      90.7381      -1.3376    0.0240
GMM full         3    6434   94.0908      90.1264      -1.9492    0.0294
GMM full         4    8579   94.7283      89.2795      -2.7961    0.0367
indep. pixels    1     128   38.3108      37.7002     -54.3754    0.2992
C3 max |hand - scipy| exact log-density: 3.1e-13; C4 max |hand PPCA - PCA.score|: 1.3e-13; C5 max |hand Gaussian - GaussianMixture.score|: 7.1e-14

seed  C1  best q  best k  max (model - exact)/SE  PPCA4 - indep (SE)  PPCA4 - full (SE)  all GMM converged  C3, C4, C5
   0   4       4       1               -14.6489    54.2979 (0.2984)   0.4868 (0.0146)               True  3.1e-13, 1.3e-13, 7.1e-14
   1   4       4       1               -14.8985    54.4303 (0.3032)   0.4828 (0.0142)               True  3.3e-13, 1.1e-13, 7.1e-14
   2   4       4       1               -11.2098    53.9672 (0.2911)   0.4544 (0.0142)               True  4.0e-13, 1.1e-13, 3.6e-14
   3   4       4       1               -13.3940    54.0090 (0.2819)   0.4687 (0.0143)               True  3.7e-13, 1.1e-13, 2.1e-14
   4   4       4       1               -12.7442    53.7771 (0.3081)   0.4871 (0.0142)               True  2.8e-13, 1.7e-13, 1.3e-13

Every check and hypothesis was written down, with its tolerance, before any script ran, so each is a test rather than an illustration. The checks held on every seed: the basis is orthonormal to $6.7 \times 10^{-16}$; exactly four training-covariance eigenvalues exceed $2\sigma^2 = 0.005$ (the fifth, 0.00341 on seed 0, is above $\sigma^2 = 0.0025$, so a bar at $\sigma^2$ would not separate signal from noise); our Gaussian log-density, probabilistic PCA and single Gaussians agree with SciPy, PCA.score and GaussianMixture.score to $4 \times 10^{-13}$ or better, against tolerances of $10^{-9}$ and $10^{-6}$; and every mixture fit converged.

All five hypotheses held on all five seeds. Probabilistic PCA scored best at $q = 4$, the true number of components (H1). No model scored above the exact density; even the closest model on each seed was 11.2 to 14.9 paired standard errors below it (H2). It follows from the definition of the KL divergence that this must hold in expectation: the expected held-out difference is $-D_{\mathrm{KL}}(p \,\|\, q) \le 0$. The independent-pixel Gaussian fell 53.8 to 54.4 nats below probabilistic PCA with four components, with standard errors near 0.3 (H3): treating pixels as independent throws away the correlations the basis images create. The mixtures scored best with one component (H4), the true family. Probabilistic PCA with four components beat the single full-covariance Gaussian by 0.45 to 0.49 nats, standard errors 0.014 to 0.015 (H5), although both families contain the true density; in our reading the full covariance’s 2,144 free parameters, against 315, cost it more estimation error with 2,000 patches.

The training column shows why held-out data are needed. On seed 0 the training log-likelihood rises at every step in both families, while the held-out value peaks at $q = 4$ and falls with every added mixture component. The training likelihood would choose the largest model in both.

Seven changes from the batch plan were made before the pilot: a kernel density estimate over a bandwidth grid was left out, with its hypothesis and its best-bandwidth measure, because the GAN paper article’s worked example already measures it against an exact density on $8 \times 8$ images; the independent-pixel Gaussian replaced it as the third family; hypotheses H4 and H5 and the training-likelihood and parameter-count columns were added; data seeds 0–4 were added (the plan named none); the eigenvalue check “above $\sigma^2$” was made precise as “above $2\sigma^2$”; the library cross-checks were added; and a standard error is reported both for the exact density’s mean and for each paired difference. The two library cross-checks were first written for seed 0 only, and the exact-density check for no stated seed; all three were run on every seed, with the same tolerances, before the pilot. The approved pilot and the first run are the same script with the same output.

Common Misconceptions

• “Unsupervised means assumption-free.” The choice of features, distance, scale and model family does the work that labels do in supervised learning; the clustering article gives Kleinberg’s theorem that no clustering function has all three of the properties he asks for.
• “A higher likelihood on the training images means a better model.” Training likelihood rises with every component, as in the example; held-out likelihood or a penalized criterion such as BIC must choose the model’s size.
• “A log-likelihood is a probability, so it cannot be positive.” For continuous pixels it is a log density: the example’s values are near +92 nats per patch, and they change if the intensities are rescaled.
• “Self-supervised learning is something else.” It learns from unlabelled images and Murphy (§1.3.3) treats it as an approach to unsupervised learning; what distinguishes it is the proxy supervised task, as the self-supervised learning article describes.
• “A good density model gives good features.” A model can assign high probability to the data without learning useful high-level patterns (Murphy, §1.3.4), which is why downstream evaluation exists.

Where It Is Used

• Visual vocabularies and descriptor encoding: clustering local descriptors into the words of a bag of visual words.
• Colour models and segmentation: colour Gaussian mixtures in segmentation energies, k-means colour quantization, mean shift segmentation and superpixels.
• Background modelling: per-pixel mixtures in background subtraction.
• Anomaly detection: a model of normal images that flags unlikely ones, as in visual anomaly detection.
• Compression and visualization: principal component analysis for descriptors and patches, t-SNE for inspecting feature spaces.
• Pretraining and generation: self-supervised encoders that serve as starting points for later tasks, and generative models that synthesize images.