Eigenfaces

Eigenfaces is the face recognition method in which every aligned face image is treated as one point in a space with one coordinate per pixel, the principal components of a training set of faces are computed, and a face is represented by its few coordinates along those components. Matthew Turk and Alex Pentland (1991) called the components eigenfaces, because they are the eigenvectors (principal components) of the set of faces; according to their abstract, they need not correspond to features such as eyes, ears and noses. A new face is recognized by projecting it onto the eigenfaces and comparing its weights with those of known people, and the same projection gives a test of whether an image region is a face at all: its distance from the subspace the eigenfaces span, which later papers by the same group call the distance from face space. The method followed work by Lawrence Sirovich and Michael Kirby (1987, 1990), who used the same decomposition, which they called eigenpictures, to represent and reconstruct face images compactly. Peter Belhumeur, João Hespanha and David Kriegman (1997) replaced principal components by Fisher’s linear discriminant, a variant they called Fisherfaces, and showed how much the eigenface coordinates are moved by lighting. The method is also called the Karhunen–Loève or eigenspace approach to face recognition.

The recognition task itself, its verification and identification protocols and today’s deep systems are covered in the face recognition article, and the decomposition in general in principal component analysis.

A note on sources. Turk and Pentland’s journal paper and its conference version (CVPR 1991) could not be read for this article: the publisher’s site refused automated access, no authors’ copy was found, and the Media Laboratory’s technical report of the same title is listed in its index as not online. Their method is described here from their abstract and from the papers that describe or reimplement it: Belhumeur et al. (1997), two papers by Pentland’s group (Pentland, Moghaddam and Starner, 1994; Moghaddam and Pentland, 1995), the U.S. National Institute of Standards and Technology (NIST) report on the FERET evaluation (Phillips et al., 1998) and Szeliski’s textbook (2022). Sirovich and Kirby’s papers could not be read either, and are attributed through the same sources.

Problem

Given images of known people, the gallery, and a new image of a face, decide who it is. Turk and Pentland treated face recognition as a two-dimensional problem rather than one requiring recovery of three-dimensional shape, relying, in the words of their abstract, on faces normally being upright and so describable by a small set of two-dimensional characteristic views. An image of $h \times w$ pixels becomes a vector of $P = hw$ numbers, so a $100 \times 100$ face is a point in a 10,000-dimensional space. Comparing such vectors directly, nearest-neighbour matching on the raw pixels, is the baseline that Belhumeur et al. call correlation. They list its disadvantages: it needs a gallery that densely samples the lighting conditions, every probe must be compared with every stored image, and every stored image must be kept. Eigenfaces address the cost and storage by comparing a few dozen coefficients instead of thousands of pixels, and add a way to say whether an image is a face.

Inputs and Outputs

Training. A set of $N$ face images, cropped and aligned so that the eyes and mouth fall at about the same pixels, all of the same size; in most uses the training set is the gallery itself or a representative sample of it. The output is the mean face $\boldsymbol{\mu}$ and $M$ eigenfaces $\mathbf{u}_1, \dots, \mathbf{u}_M$, each an image of $P$ pixels, together with the projection of every gallery image, a vector of $M$ weights.

Recognition. An aligned probe image goes in; out comes the identity of the nearest gallery entry in eigenface coordinates and the distance to it, which can be thresholded to reject unknown faces.

Detection. An image window goes in; out comes its distance from face space, low for faces. Scanning it over positions and scales, as in sliding window detection, gives a map whose minima are face candidates, an early form of face detection.

Intuition

Face images are far from uniformly spread in pixel space. All of them have two eyes above a nose above a mouth, the same overall shading, and pixels that vary together. Principal component analysis finds the few directions in which a set of face images varies most; most of the variation of the set lies in a subspace of much lower dimension than the number of pixels. Projecting onto that subspace keeps what distinguishes one training face from another, as far as variance goes, and discards directions in which the faces hardly vary. Two images of the same person should then land near each other, and an image that is not a face should land far from the subspace altogether, because no combination of eigenfaces reconstructs it well.

The weakness is in “as far as variance goes”. The components are chosen to maximize the scatter of all images, whatever its cause. If the training set contains the same people under different lighting, the largest variation can be the lighting, and the leading eigenfaces then encode where the light comes from rather than who is in the picture. Belhumeur et al. built Fisherfaces on exactly this observation: directions should be chosen by how well they separate people, not by how much they vary.

Algorithm

Training:

1. Crop, align and resize every training image to $h \times w$ pixels, optionally normalize its brightness and contrast, and stack each as a column vector $\mathbf{x}_i$ of length $P$.
2. Compute the mean face $\boldsymbol{\mu}$ and the centred vectors $\mathbf{a}_i = \mathbf{x}_i - \boldsymbol{\mu}$; collect them as the columns of $A$ ($P \times N$).
3. Form the $N \times N$ matrix $A^\top A$ of inner products between centred images and compute its eigenvectors $\mathbf{v}_j$ and eigenvalues $\lambda_j$.
4. Map each of the $M$ leading eigenvectors back to image space, $\mathbf{u}_j = A \mathbf{v}_j / \sqrt{\lambda_j}$: these are the eigenfaces, orthonormal images of $P$ pixels.
5. Project every gallery image: $\boldsymbol{\omega} = U^\top (\mathbf{x} - \boldsymbol{\mu})$, with $U = [\mathbf{u}_1 \cdots \mathbf{u}_M]$.

Recognition:

6. Align and preprocess the probe exactly as the training images, project it with step 5 and find the gallery vector at the smallest distance. Accept the identity if that distance is below a threshold, otherwise report an unknown face.

Detection:

7. For each window of the image, preprocess it, project it and compute the residual between the window and its reconstruction from the $M$ eigenfaces, the distance from face space. Windows whose residual is below a threshold, or local minima of the residual map, are face candidates.

Mathematical Formulation

With the centred images as columns of $A$, the sample covariance of the training set is

$$C = \frac{1}{N} A A^\top = \frac{1}{N} \sum_{i=1}^{N} (\mathbf{x}_i - \boldsymbol{\mu})(\mathbf{x}_i - \boldsymbol{\mu})^\top ,$$

a $P \times P$ matrix, 10,000 by 10,000 for $100 \times 100$ images. Its eigenvectors with the largest eigenvalues are the eigenfaces (Szeliski, 2022, §5.2.3), and Belhumeur et al. write the same criterion as the projection $W$ with orthonormal columns that maximizes the determinant of the projected total scatter, $W_{\text{pca}} = \arg\max_W |W^\top S_T W|$, where $S_T = N C$.

The small-matrix trick. With $N$ much smaller than $P$, $C$ has rank at most $N - 1$ and it never needs to be formed. If $A^\top A \mathbf{v} = \lambda \mathbf{v}$, multiplying by $A$ gives $A A^\top (A\mathbf{v}) = \lambda (A\mathbf{v})$, so $A\mathbf{v}$ is an eigenvector of $A A^\top$ with the same eigenvalue, and since $\lVert A\mathbf{v} \rVert^2 = \mathbf{v}^\top A^\top A \mathbf{v} = \lambda$, the unit-norm eigenface is

$$\mathbf{u} = \frac{A \mathbf{v}}{\sqrt{\lambda}} .$$

Szeliski (2022, Appendix A.1.2) states this in terms of the singular value decomposition $A = U \Sigma V^\top$: the eigenvectors of the $N \times N$ matrix $A^\top A$ are the right singular vectors $V$, the eigenfaces are the left singular vectors $U = A V \Sigma^{-1}$, and computing them this way amounts to forming the eigenfaces as linear combinations of the difference images, which he attributes to Turk and Pentland (1991). The same identity is derived for PCA in general in principal component analysis, and the singular value decomposition of $A$ contains it, working with whichever dimension is smaller.

Projection and reconstruction. The weights of an image are $\omega_j = \mathbf{u}_j^\top (\mathbf{x} - \boldsymbol{\mu})$, and its reconstruction from $M$ eigenfaces is $\tilde{\mathbf{x}} = \boldsymbol{\mu} + \sum_{j=1}^{M} \omega_j \mathbf{u}_j$. Truncating after $M$ terms gives the best approximation in the least-squares sense among all $M$-term expansions (Szeliski, 2022, §5.2.3).

Recognition. The probe is assigned to the gallery entry $k$ that minimizes $\lVert \boldsymbol{\omega} - \boldsymbol{\omega}_k \rVert$; the k-nearest neighbours article covers the classifier and the similarity and distance measures article the alternatives to the Euclidean norm. FERET’s baseline implementation, for instance, used the L1 distance (Phillips et al., 1998). When the eigenfaces are computed from the gallery itself and all $N - 1$ components are kept, distances between projected gallery images equal their distances in pixel space, and the projection of a probe differs from the probe only by its component orthogonal to the training images, which is the same for every gallery entry; nearest-neighbour matching on all components therefore gives the same answers as matching on raw pixels. Belhumeur et al. state this equivalence with correlation and confirm it experimentally.

Distance from face space. Because the eigenfaces are orthonormal, the squared distance between an image and its reconstruction is

$$\varepsilon^2(\mathbf{x}) = \lVert \mathbf{x} - \tilde{\mathbf{x}} \rVert^2 = \lVert \mathbf{x} - \boldsymbol{\mu} \rVert^2 - \sum_{j=1}^{M} \omega_j^2 ,$$

computable from the $M$ weights and the energy of the centred image. Moghaddam and Pentland (1995, Eq. 2) write it in this form and say that it was called the “distance-from-face-space” in the context of the work with eigenfaces; Pentland, Moghaddam and Starner (1994) also credit the term to Turk and Pentland (1991). For a general principal subspace, Moghaddam and Pentland call the residual the distance-from-feature-space (DFFS) and the complementary part, the position of the projection within the subspace, the distance-in-feature-space (DIFS). Szeliski (2022), citing Moghaddam and Pentland’s later journal paper, calls the two the distance from face space (DFFS) and the distance in face space (DIFS), and describes the first as measuring the “faceness” of an image.

Fisherfaces. With $c$ people, class means $\boldsymbol{\mu}_i$ and $N_i$ images of person $i$, Belhumeur et al. define the between-class and within-class scatter

$$S_B = \sum_{i=1}^{c} N_i (\boldsymbol{\mu}_i - \boldsymbol{\mu})(\boldsymbol{\mu}_i - \boldsymbol{\mu})^\top , \qquad S_W = \sum_{i=1}^{c} \sum_{\mathbf{x}_k \in X_i} (\mathbf{x}_k - \boldsymbol{\mu}_i)(\mathbf{x}_k - \boldsymbol{\mu}_i)^\top ,$$

and Fisher’s linear discriminant chooses $W = \arg\max_W |W^\top S_B W| / |W^\top S_W W|$, whose columns are the generalized eigenvectors $S_B \mathbf{w}_i = \lambda_i S_W \mathbf{w}_i$ with the largest eigenvalues. At most $c - 1$ eigenvalues are nonzero, so the projection has at most $c - 1$ dimensions. For images, $S_W$ has rank at most $N - c$ and is singular whenever $N < P$, the usual case. Their solution is to apply PCA first, reducing the dimension to $N - c$ so that the within-class scatter becomes nonsingular, and then the discriminant, reducing it to $c - 1$:

$$W_{\text{opt}}^\top = W_{\text{fld}}^\top W_{\text{pca}}^\top , \qquad W_{\text{fld}} = \arg\max_W \frac{|W^\top W_{\text{pca}}^\top S_B W_{\text{pca}} W|}{|W^\top W_{\text{pca}}^\top S_W W_{\text{pca}} W|} .$$

In computing $W_{\text{pca}}$ only the smallest $c - 1$ principal components are discarded (Belhumeur et al., 1997, §2.4).

Parameters

• Number of eigenfaces $M$. Recognition with too few components loses identity. On the Yale database, Belhumeur et al. found that the error of eigenfaces fell as components were added; with the first three components removed, their text says it levelled off at about 45 components, while in their Figure 8 the lowest point of that curve is at about 60 components and the error rises slightly beyond it. They note that Sirovich and Kirby had found a similar point of diminishing returns when representing face images. Pentland et al. (1994) used 20 eigenvectors, computed from a sample of 128 faces, to search a database of 7,562 images; FERET’s baseline used 200, computed from 500 faces (Phillips et al., 1998).
• Discarding leading components. Belhumeur et al. report the suggestion, which they do not attribute, that the first three principal components are mainly due to lighting and that dropping them improves recognition; they tested components 4 to 13 as an alternative to 1 to 10.
• Preprocessing. Belhumeur et al. normalized each image to zero mean and unit variance for the eigenface and correlation methods because it improved their performance, which also makes the result independent of light source intensity and a camera’s automatic gain. Moghaddam and Pentland (1995) obtained illumination invariance in their detector by normalizing each subimage to a zero-mean, unit-norm vector. FERET’s baseline translated, rotated and scaled each image to put the eye centres on fixed pixels, masked out background and hair, and histogram-equalized the remaining pixels (Phillips et al., 1998).
• Crop. What the crop includes changes the result: on the Yale database every method had lower error on full-face images, which included the face’s outline and some background, than on images cropped to the brow, eyes, nose, mouth and chin, and Belhumeur et al. warn that the full-face advantage might disappear if backgrounds or hair styles varied.
• Distance and thresholds. The distance used for matching (Euclidean, L1, or a Mahalanobis distance that weights each coefficient by its eigenvalue), the threshold for rejecting unknown faces, and the threshold on the distance from face space for detection all have to be set on held-out data; the sources read for this article do not give Turk and Pentland’s values.

Complexity

The following counts are our own, for $N$ training images of $P$ pixels, $M$ eigenfaces and a gallery of $K$ images. Forming $A^\top A$ costs $O(N^2 P)$ and its eigendecomposition $O(N^3)$, against $O(P^3)$ for the eigendecomposition of the full covariance; mapping $M$ eigenvectors back costs $O(MNP)$. With $N = 180$ and $P = 2{,}304$ in the example below, the eigenproblem is $180 \times 180$ instead of $2{,}304 \times 2{,}304$. Projecting a probe costs $O(MP)$, and matching it against the gallery $O(MK)$, against $O(PK)$ for raw pixels; the gallery takes $MK$ numbers instead of $PK$. For detection, each window costs one projection, so the residual map of an image needs $M$ correlations of the image with the eigenfaces plus a local energy computation, which Pentland et al. (1994) note can be computed efficiently with the FFT.

Implementation

The example renders synthetic face-like images and reproduces, on a small scale, the effects the sources describe. Each image is $48 \times 48$ pixels and is rendered as a Lambertian surface: an ellipsoidal head with a nose bump gives the surface normals, an albedo map gives darker eyes, brows and mouth, and the brightness is the albedo times the clipped cosine between the normal and the light direction, plus a small ambient term and noise. Twenty synthetic people differ in head shape, eye spacing, height and size, nose size, mouth position and width, and skin albedo; each rendering adds a jitter of up to half a pixel to the features. Light directions lie on a 15° grid of azimuth and elevation from −60° to 60°, grouped like Belhumeur et al.’s Harvard subsets by the larger of the two angles. The gallery holds nine renderings per person with near-frontal light (both angles within 15°), 180 images in all.

The script checks the small-matrix trick against the SVD (A), compares distances due to lighting and to identity (B), trains on the near-frontal gallery and recognizes fresh renderings under every light with five methods (C), measures what the leading components encode when the training set spans all lights (D), shifts the probes to test alignment (E), uses the distance from face space to separate faces from patches of scikit-image’s bundled camera photograph and to search that photograph for a pasted synthetic face (F), and repeats part C over five seeds (G). Belhumeur et al. normalized images for eigenfaces and correlation, and Turk and Pentland’s preprocessing could not be read, so parts C and F print every method with and without normalization. It runs in about 8 seconds with NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2.

import numpy as np
from scipy.linalg import eigh
from numpy.lib.stride_tricks import sliding_window_view
from skimage import data

H = W = 48
yy, xx = np.mgrid[0:H, 0:W].astype(float)
ANGLES = np.arange(-60, 61, 15)
LIGHTS = [(az, el) for el in ANGLES for az in ANGLES]
SUBSET = np.array([max(abs(az), abs(el)) for az, el in LIGHTS])  # 0,15,30,45,60


def identity(rng):
    """Shape and albedo parameters of one synthetic person."""
    return dict(a=rng.uniform(13, 17), b=rng.uniform(17, 21),
                eye_dx=rng.uniform(6, 9), eye_y=rng.uniform(-8, -4),
                eye_r=rng.uniform(1.6, 2.8), nose_h=rng.uniform(3, 7),
                nose_len=rng.uniform(3, 6), mouth_y=rng.uniform(7, 11),
                mouth_w=rng.uniform(4, 8), albedo=rng.uniform(0.6, 1.0),
                brow=rng.uniform(0.3, 0.7))


def render(p, az, el, rng, shift=0.0):
    """Lambertian rendering: ellipsoid head + nose bump, albedo features."""
    j = lambda: rng.uniform(-0.5, 0.5)            # small 'expression' jitter
    X, Y = xx - 23.5 - shift, yy - 23.5
    r2 = (X / p['a'])**2 + (Y / p['b'])**2
    inside = r2 < 1
    z = 12 * np.sqrt(np.clip(1 - r2, 0, None))
    z += p['nose_h'] * np.exp(-X**2 / 4.5 - (Y - 1)**2 / (2 * p['nose_len']**2))
    gy, gx = np.gradient(z)
    n = np.stack([-gx, -gy, np.ones_like(z)])
    n /= np.linalg.norm(n, axis=0)
    a, e = np.radians(az), np.radians(el)
    s = np.array([np.sin(a) * np.cos(e), -np.sin(e), np.cos(a) * np.cos(e)])
    shade = np.clip(np.tensordot(s, n, 1), 0, None)
    alb = np.full((H, W), p['albedo'])
    for side in (-1, 1):
        ex, ey = side * p['eye_dx'] + j(), p['eye_y'] + j()
        alb[(X - ex)**2 + (Y - ey)**2 < p['eye_r']**2] = 0.15
        brow = (abs(X - ex) < p['eye_r'] + 1) & (abs(Y - ey + 3.5) < 0.8)
        alb[brow] *= p['brow']
    mx, my = j(), p['mouth_y'] + j()
    alb[((X - mx) / p['mouth_w'])**2 + ((Y - my) / 1.2)**2 < 1] = 0.3
    img = inside * alb * (0.1 + 0.9 * shade)
    return (img + rng.normal(0, 0.01, img.shape)).ravel()


def pca(X, k):
    """Eigenfaces by the small-matrix trick: eigenvectors of the N x N Gram matrix."""
    mu = X.mean(0)
    A = (X - mu).T                                # P x N, columns are centred images
    lam, V = np.linalg.eigh(A.T @ A)              # N x N instead of P x P
    order = np.argsort(lam)[::-1][:k]
    lam, V = lam[order], V[:, order]
    U = A @ V / np.sqrt(lam)                      # unit-norm eigenfaces, P x k
    return mu, U, lam


def normalize(X):
    X = X - X.mean(1, keepdims=True)
    return X / X.std(1, keepdims=True)


def features(G, labels, c):
    """Projection functions for each method, fitted on the gallery G."""
    N = len(G)
    mu, U, _ = pca(G, N - 1)
    Up = U[:, :N - c]                             # Fisherfaces step 1: PCA to N - c
    Y = (G - mu) @ Up
    m = Y.mean(0)
    Sb = np.zeros((N - c, N - c)); Sw = np.zeros_like(Sb)
    for k in range(c):
        Yk = Y[labels == k]; mk = Yk.mean(0)
        Sb += len(Yk) * np.outer(mk - m, mk - m)
        Sw += (Yk - mk).T @ (Yk - mk)
    w, V = eigh(Sb, Sw)                           # generalized eigenproblem
    Wf = Up @ V[:, np.argsort(w)[::-1][:c - 1]]   # Fisherfaces step 2: LDA to c - 1
    return {'pixels': lambda X: X,
            'eigen M=10': lambda X: (X - mu) @ U[:, :10],
            'eigen 4-13': lambda X: (X - mu) @ U[:, 3:13],
            f'eigen M={N-1}': lambda X: (X - mu) @ U,
            f'fisher M={c-1}': lambda X: (X - mu) @ Wf}


def nn_error(f, G, gl, P, pl):
    g, q = f(G), f(P)
    d = (q**2).sum(1)[:, None] - 2 * q @ g.T + (g**2).sum(1)[None, :]
    return (gl[d.argmin(1)] != pl).mean()


def experiment(seed, c=20, verbose=True):
    rng = np.random.default_rng(seed)
    people = [identity(rng) for _ in range(c)]
    gal_idx = np.flatnonzero(SUBSET <= 15)        # 9 near-frontal lights
    G = np.array([render(p, *LIGHTS[i], rng) for p in people for i in gal_idx])
    gl = np.repeat(np.arange(c), len(gal_idx))
    P = np.array([render(p, *L, rng) for p in people for L in LIGHTS])
    pl = np.repeat(np.arange(c), len(LIGHTS))
    psub = np.tile(SUBSET, c)
    out = {}
    for norm in (False, True):
        Gx, Px = (normalize(G), normalize(P)) if norm else (G, P)
        for name, f in features(Gx, gl, c).items():
            out[(norm, name)] = [nn_error(f, Gx, gl, Px[psub == s], pl[psub == s])
                                 for s in (0, 15, 30, 45, 60)]
    return rng, people, G, gl, P, pl, psub, out


rng, people, G, gl, P, pl, psub, out = experiment(seed=0)
c, N = len(people), len(G)

# A. small-matrix trick against the SVD of the data
mu, U, lam = pca(G, 20)
_, sv, Vt = np.linalg.svd(G - mu, full_matrices=False)
print(f"A. images {G.shape[1]} px, N = {N}: Gram matrix {N}x{N} "
      f"instead of covariance {G.shape[1]}x{G.shape[1]}")
print(f"   max |eigenvalue - sigma^2| / sigma_1^2 = {np.max(abs(lam - sv[:20]**2)) / sv[0]**2:.1e}, "
      f"max |u_gram| - |u_svd| = {np.max(abs(abs(U) - abs(Vt[:20].T))):.1e}")

# B. lighting versus identity in pixel space (Moses et al.)
front = np.array([LIGHTS.index((0, 0)), LIGHTS.index((45, 0)), LIGHTS.index((0, 45))])
Pf = P.reshape(c, len(LIGHTS), -1)
same = np.mean([np.linalg.norm(Pf[k, front[0]] - Pf[k, front[1:]], axis=1).mean() for k in range(c)])
diff = np.mean([np.linalg.norm(Pf[k, front[0]] - Pf[l, front[0]]) for k in range(c) for l in range(c) if k != l])
print(f"B. mean distance, same person frontal vs 45-degree light: {same:.2f}; "
      f"different people, same frontal light: {diff:.2f}")

# C. extrapolation from near-frontal lighting
print("C. nearest-neighbour error rate (%) by light subset; gallery = subset <=15 deg")
print(f"   {'method':14s} {'norm':>5s}" + "".join(f"{s:>7d}" for s in (0, 15, 30, 45, 60)))
for (norm, name), e in out.items():
    print(f"   {name:14s} {'yes' if norm else 'no':>5s}" + "".join(f"{100*v:7.1f}" for v in e))

# D. what the leading components encode when training spans all lights
Xall = P                                          # one render of each person at each light
mu_a, U_a, lam_a = pca(Xall, 8)
coef = (Xall - mu_a) @ U_a
total_var = ((Xall - mu_a)**2).sum()
s_xy = np.array([[np.sin(np.radians(az)) * np.cos(np.radians(el)), np.sin(np.radians(el))]
                 for az, el in LIGHTS])
sx, sy = np.tile(s_xy[:, 0], c), np.tile(s_xy[:, 1], c)
print("D. PCA on all 81 lights x 20 people: leading components")
print("   comp  var%  identity-share  corr(light x)  corr(light y)")
for k in range(8):
    ck = coef[:, k]
    share = np.var([ck[pl == i].mean() for i in range(c)]) / ck.var()
    print(f"   {k+1:4d} {100*lam_a[k]/total_var:5.1f} {share:15.2f}"
          f" {np.corrcoef(ck, sx)[0,1]:14.2f} {np.corrcoef(ck, sy)[0,1]:14.2f}")

# E. misalignment: frontal-light probes shifted horizontally
print("E. error rate (%) for frontal-light probes shifted by 0-4 px (normalized images)")
fs = features(normalize(G), gl, c)
for name in ('pixels', 'eigen M=10', f'fisher M={c-1}'):
    errs = []
    for sh in range(5):
        Ps = normalize(np.array([render(p, 0, 0, rng, shift=sh) for p in people for _ in range(5)]))
        errs.append(nn_error(fs[name], normalize(G), gl, Ps, np.repeat(np.arange(c), 5)))
    print(f"   {name:14s}" + "".join(f"{100*v:7.1f}" for v in errs))

# F. distance from face space (DFFS) as a face detector
mask = (((xx - 23.5) / 12)**2 + ((yy - 23.5) / 16)**2 < 1).ravel()
def unit(X):
    X = X - X.mean(1, keepdims=True)
    return X / np.linalg.norm(X, axis=1, keepdims=True)
cam = data.camera().astype(float) / 255
win = sliding_window_view(cam, (H, W))[::4, ::4]
new = [identity(rng) for _ in range(10)]
faces = np.array([render(p, *LIGHTS[i], rng) for p in new for i in np.flatnonzero(SUBSET <= 15)])
target = render(new[0], 15, 15, rng).reshape(H, W)
scene = cam.copy(); r0, c0 = 200, 300
patch = scene[r0:r0+H, c0:c0+W]
scene[r0:r0+H, c0:c0+W] = np.where(target > 0.02, target, patch)
swin = sliding_window_view(scene, (H, W))[::4, ::4]
print("F. face/non-face scores with 10 eigenfaces (pixels inside a fixed ellipse)")
for norm, ml in ((False, False), (True, False), (True, True)):
    prep = unit if norm else (lambda X: X)
    Gm = prep(G[:, mask])
    mu_f, U_f, lam_f = pca(Gm, 10)
    var = lam_f / len(Gm)                         # covariance eigenvalues
    rho = (((Gm - mu_f)**2).sum() / len(Gm) - var.sum()) / (mask.sum() - 10)
    def score(X):
        Z = prep(X[:, mask]) - mu_f
        y = Z @ U_f
        dffs2 = ((Z - y @ U_f.T)**2).sum(1)       # squared distance from face space
        # Moghaddam & Pentland (1995): add the distance-in-feature-space (DIFS),
        # Mahalanobis-weighted, and weight the DFFS by the mean residual variance
        return (y**2 / var).sum(1) + dffs2 / rho if ml else np.sqrt(dffs2)
    df, dc = score(faces), score(win.reshape(-1, H * W))
    thr = np.quantile(df, 0.95)
    smap = score(swin.reshape(-1, H * W)).reshape(swin.shape[:2])
    r, cc = np.unravel_index(smap.argmin(), smap.shape)
    name = 'DIFS+DFFS' if ml else 'DFFS'
    print(f"   {name:9s} normalized={str(norm):5s} camera patches scoring below 95% of faces: "
          f"{100*(dc < thr).mean():5.1f}%; scan minimum at ({4*r},{4*cc}), face at ({r0},{c0})")

# G. spread over seeds of the 45-degree error (normalized images)
print("G. error at 45 deg over seeds 0-4 (normalized), min-max %")
res = {}
for seed in range(5):
    o = experiment(seed)[-1]
    for (norm, name), e in o.items():
        if norm:
            res.setdefault(name, []).append(100 * e[3])
for name, v in res.items():
    print(f"   {name:14s} {min(v):5.1f} - {max(v):5.1f}")

Output:

A. images 2304 px, N = 180: Gram matrix 180x180 instead of covariance 2304x2304
   max |eigenvalue - sigma^2| / sigma_1^2 = 9.8e-16, max |u_gram| - |u_svd| = 3.1e-15
B. mean distance, same person frontal vs 45-degree light: 7.87; different people, same frontal light: 7.42
C. nearest-neighbour error rate (%) by light subset; gallery = subset <=15 deg
   method          norm      0     15     30     45     60
   pixels            no    0.0    0.6    2.5   37.3   71.7
   eigen M=10        no    0.0    0.6   10.0   64.0   82.2
   eigen 4-13        no    0.0    0.6    2.5    7.3   26.2
   eigen M=179       no    0.0    0.6    2.5   37.3   71.7
   fisher M=19       no    0.0    2.5    9.4   30.8   51.7
   pixels           yes    5.0    1.9    0.0    0.8    1.6
   eigen M=10       yes    5.0    2.5    5.3   12.3   19.5
   eigen 4-13       yes    0.0    1.9    2.8   20.8   34.4
   eigen M=179      yes    5.0    1.9    0.0    0.8    1.6
   fisher M=19      yes    0.0    0.0    1.9   18.1   41.2
D. PCA on all 81 lights x 20 people: leading components
   comp  var%  identity-share  corr(light x)  corr(light y)
      1  28.2            0.00           0.98           0.00
      2  25.6            0.00           0.00          -0.98
      3  17.7            0.41          -0.00           0.01
      4   5.5            0.94           0.00           0.01
      5   3.6            0.93           0.00          -0.02
      6   1.8            0.00           0.01           0.00
      7   1.4            0.61          -0.00           0.00
      8   1.3            0.30          -0.00          -0.02
E. error rate (%) for frontal-light probes shifted by 0-4 px (normalized images)
   pixels            0.0    4.0   32.0   50.0   76.0
   eigen M=10        1.0    5.0   19.0   35.0   62.0
   fisher M=19       0.0   74.0   78.0   87.0   93.0
F. face/non-face scores with 10 eigenfaces (pixels inside a fixed ellipse)
   DFFS      normalized=False camera patches scoring below 95% of faces:  75.9%; scan minimum at (192,188), face at (200,300)
   DFFS      normalized=True  camera patches scoring below 95% of faces:   9.9%; scan minimum at (20,108), face at (200,300)
   DIFS+DFFS normalized=True  camera patches scoring below 95% of faces:   9.7%; scan minimum at (20,108), face at (200,300)
G. error at 45 deg over seeds 0-4 (normalized), min-max %
   pixels           0.0 -   0.8
   eigen M=10       5.4 -  12.3
   eigen 4-13       3.1 -  20.8
   eigen M=179      0.0 -   0.8
   fisher M=19     18.1 -  78.3

Two changes were made after the first run. Part D first recomputed a $1{,}620 \times 1{,}620$ matrix product in every loop iteration to obtain the total variance; it now computes the total once, which cut the running time from about half a minute and left every printed value of parts A to E unchanged. In part F, the first run printed only the two distance-from-face-space lines, with the median distances of faces and patches. Its scan minimum fell on the sky of the camera image instead of on the pasted face, against the plan’s expectation that it would land “within a few pixels of the pasted face, at least with normalized windows”. The likelihood score of Moghaddam and Pentland (1995) was then added as a third line, and the print format was changed to label each line and drop the medians; no setting was changed after the added line was run.

A. The eigenvalues and eigenvectors from the $180 \times 180$ Gram matrix agree with the SVD of the data to rounding error, so the trick is exact, not an approximation.

B. Changing the light from frontal to 45° moves an image further, on average, than replacing the person by another under the same light. This is the observation of Moses, Adini and Ullman that Belhumeur et al. quote: variations due to illumination and viewing direction are almost always larger than variations due to identity. Any method that compares raw distances will be confused by lighting it has not seen.

C. With raw images, every method is nearly perfect on light like the gallery’s and degrades as the light moves away. Ten eigenfaces do worse than raw pixels at 30° and beyond, because, by our reading, the leading components of even a near-frontal gallery describe the small lighting changes in it and give probe lighting a large weight in the matching. Dropping the first three components has the largest effect in the table, in the same direction as in Belhumeur et al.’s Harvard experiments (see Reported Results). Eigenfaces with all 179 components reproduce the pixel results exactly, the equivalence derived above. Fisherfaces beat raw pixels at 45° and 60° but not the eigenfaces without the first three.

Normalization changes the picture more than any choice of subspace. Normalized pixels are nearly perfect at every angle, because in this renderer much of what a light change does to a convex, mostly smooth head is to scale and tilt the brightness, which subtracting the mean and dividing by the standard deviation largely removes (our reading of the renderer; real faces cast shadows from the nose and brows, which this one mostly lacks). With normalization, dropping the first three eigenfaces now hurts at 45° and 60°: the trick depends on what the leading components happen to encode, which Belhumeur et al. caution against when they say it is unlikely that the first components correspond solely to lighting, so useful information may be lost. The single error on frontal probes in some normalized rows is one probe of twenty.

Fisherfaces did not do best, against the plan’s hypothesis and against Belhumeur et al.’s Harvard results, where they had the lowest or joint-lowest error on every subset of the extrapolation test. The seed loop in part G shows how unstable they are here: their error at 45° ranges from 18% to 78% over five seeds, while normalized pixels never exceed 0.8%. Our explanation, which the example does not test separately: the gallery offers nine images per person under nearly the same light, and the discriminant weights most heavily the directions in which each person’s gallery images hardly vary. Any change absent from the gallery that has a component along those directions, here the attached shadows and larger shading changes of oblique light, is magnified rather than suppressed. Belhumeur et al.’s data were real faces and their method was trained on a lighting subset with the same 15° range; the toy differs in its rendering and its small, nearly noise-free variation within each class.

D. When the training set contains every person under every light, the first two components have no identity content at all and correlate almost perfectly with the horizontal and vertical components of the light direction; the third mixes the two, and identity appears from the fourth on, interrupted at the sixth by another component with no identity content. This is the effect Belhumeur et al. describe: principal components that retain the variation due to lighting, so that the projected points are not well clustered and the classes may be smeared together.

E. With normalized images, a shift of two pixels raises the error of raw pixels to about a third, and three pixels to half; ten eigenfaces degrade somewhat more slowly, but not enough to make alignment optional. Fisherfaces fail after a single pixel, the same magnification of nuisances that are absent from the gallery as in part C.

F. On raw windows the distance from face space hardly separates faces from the photograph: three quarters of the camera patches score as well as the 95th percentile of faces. Normalizing each window to zero mean and unit norm, as Moghaddam and Pentland did, brings that to about one in ten. No line finds the pasted face: the raw distance puts its minimum elsewhere in the photograph, and the normalized distance and the likelihood score that adds the distance-in-feature-space both put theirs on the same window in the sky. By our reading, a normalized smooth patch is mostly a gentle brightness gradient, which the eigenfaces of a lighting-varied training set can reconstruct well, while ten eigenfaces of twenty people reconstruct a new person’s face only roughly. The detectors in the sources were built from real face images, and Moghaddam and Pentland’s searched over scale as well as position.

Properties and Behavior

Lighting. Under a Lambertian model without shadowing, the image intensity at a point is linear in the light source vector, so images of one face under any light from a fixed viewpoint lie in a three-dimensional linear subspace of image space; real faces deviate from it because of self-shadowing, specularities and expression (Belhumeur et al., 1997, §2.3). The variation this produces is not separated from identity by PCA, which maximizes the total scatter. Belhumeur et al.’s “Linear Subspace” method used the observation directly, representing each person by three or more images under different lights and classifying a probe by its distance to each person’s subspace; Fisherfaces instead learn a projection in which the within-class scatter is small, without modelling the subspace.

Pose. A single eigenspace for all views has to describe face images under rotation, which Pentland et al. (1994) note form non-convex, if connected, regions of image space. Their view-based alternative builds a separate eigenspace for each view, chooses the view whose eigenspace best describes the input by the distance from face space, and recognizes in that eigenspace. Belhumeur et al. made no attempt to deal with pose and assumed the face had been located and aligned.

Alignment and scale. Because each pixel is one coordinate, a face shifted by a few pixels, scaled or rotated in the image plane is a different vector. The example above measures how quickly nearest-neighbour matching fails, and Pentland et al. (1994) describe registration, locating the face and placing it in a standard-size frame, as an important preprocessing step in an eigenvector recognition system.

Occlusion and local change. Pentland et al. (1994) show test views with a hand near the face, a painted face and a beard, for which the nearest eigenface matches were all wrong while matches on the eyes and nose alone were correct. Their explanation is that a pure eigenface system can be fooled by gross variations in the input image.

Learning new people. According to Turk and Pentland’s abstract, the method can learn and later recognize new faces in an unsupervised manner. As our reasoning, adding a person needs only the projection of their images onto existing eigenfaces, as long as the eigenfaces still describe the new faces well; FERET’s design rules note that PCA-based algorithms are gallery-insensitive in this sense, while a straightforward Fisher discriminant is tuned to the gallery it was trained on (Phillips et al., 1998).

Reported Results

Results from the 1990s are on small private databases with different protocols and are not comparable with one another or with the deep-learning results in the face recognition article. Turk and Pentland’s own experiments are not reported here because their paper could not be read.

Belhumeur et al. (1997), Harvard database. 330 images of five people (66 each), collected by Hallinan at the Harvard Robotics Laboratory, with the light source direction sampled in 15° steps of longitude and latitude and grouped into subsets by the larger angle from the camera axis: Subset 1 within 15° (30 images), Subset 2 at 30° (45), Subset 3 at 45° (65), Subset 4 at 60° (85) and Subset 5 at 75° (105). The images were cropped within the face, excluding the head’s outline, and classified by a nearest-neighbour rule. Trained on Subset 1 and tested on Subsets 1 to 3 (their Figure 5), the error rates on Subset 3 were 41.5% for ten eigenfaces, 27.7% for ten eigenfaces without the first three, 33.9% for correlation, 9.2% for the Linear Subspace method and 4.6% for Fisherfaces with four dimensions; all methods had 0% error on Subset 1, tested by leaving each image out of the training set. Trained on Subsets 1 and 5 and tested on Subsets 2 to 4 (their Figure 6), Fisherfaces had 0.0%, 0.0% and 1.2% error and ten eigenfaces 11.11%, 33.9% and about 20%.

Belhumeur et al. (1997), Yale database. 160 frontal images of 16 people, ten each, taken in one session: females and males, some with facial hair and some with glasses, under ambient light with and without glasses, three positions of a lamp and five expressions (happy, sad, winking, sleepy, surprised). Each image was classified by leaving it out, recomputing the projection and taking the nearest neighbour. With closely cropped images the error rates were 24.4% for 30 eigenfaces, 15.3% for 30 eigenfaces without the first three, 23.9% for correlation, 21.6% for the Linear Subspace method and 7.3% for Fisherfaces (15 dimensions); with full-face images, which included the outline and part of the background, they were 19.4%, 10.8%, 20.0%, 15.6% and 0.6% (their Figure 9). In a further test on 36 images from a superset of the Yale database, half with glasses, the task was to decide whether a person wore glasses: PCA with ten components had 52.6% error, near chance, and a one-dimensional Fisherface 5.3% (their Table 1).

Pentland, Moghaddam and Starner (1994). On a database of 7,562 images of about 3,000 people, each image annotated by hand for sex, race, approximate age and expression, 200 faces chosen at random and matched by nearest neighbour with 20 eigenvectors were recognized with 95% accuracy. The authors report that the accuracy was maintained across race and sex categories, with a possible trend, not statistically significant, toward lower performance on the faces they categorize as “Oriental”. On 189 images of 21 people in nine views from −90° to +90°, view-based eigenspaces recognized 90% of intermediate views against 88% for a single parametric eigenspace, and 83% against 78% for views 23° outside the training range. For eye detection by the distance from feature space with ten eigenvectors on the large database, the best operating point detected 94% of left eyes at a false alarm rate of 6%, and the authors report an order of magnitude fewer false alarms than template matching by the sum of squared differences at the same detection rate.

Moghaddam and Pentland (1995). Their maximum-likelihood detector, which combines the distance-from- and the distance-in-feature-space, located the left eye in about 7,000 test images of the Media Laboratory’s database with a false alarm rate nearly two orders of magnitude below that of sum-of-squared-differences template matching at the same detection rate, and a multiscale version detected 97% of 2,000 faces in the FERET database.

FERET, September 1996 (Phillips et al., 1998). NIST’s evaluation ran its own PCA baseline: images aligned by the eye centres, masked and histogram-equalized, 200 eigenvectors computed from 500 faces, nearest neighbour with the L1 distance. On the FB probes (of two frontal images taken in one session, the second with a different expression requested, one went at random into the gallery of 1,196 people and the other into the probe set), its rank-1 identification rate was about 0.80, read from the report’s Figure 3, against about 0.95 for the best partially automatic algorithm of September 1996, one of two from the MIT Media Laboratory. Ranking the algorithms on six disjoint galleries of about 200 people, the baseline PCA placed eighth to tenth of ten on FB probes. The report does not break down results by demographic group.

Limitations

• Lighting, pose and expression change the coordinates. Principal components are chosen for variance, not for identity, and the variance of real face images is dominated by lighting and viewing direction (Belhumeur et al., 1997). The method works when probes resemble the gallery; all methods compared by Belhumeur et al. performed well when a test image was similar to a training image.
• It needs precise alignment. Every pixel is a feature, so the face must be located, scaled and registered before projection (part E of the example).
• Holistic representation. A local change such as an occluding hand or a beard moves the whole coefficient vector (Pentland et al., 1994), which motivated modular eigenfeatures and, later, local descriptors.
• Linear and global. One linear subspace cannot describe the non-convex set of views of a face (Pentland et al., 1994), and the distance from a single subspace is a weak detector on its own, as part F of the example shows.
• Small evaluations. The early results were on tens to hundreds of people. FERET, the first large independent test, placed the PCA baseline near the bottom of the algorithms of 1996 (Phillips et al., 1998).

Variants

• Fisherfaces (Belhumeur, Hespanha and Kriegman, 1997): PCA to $N - c$ dimensions followed by Fisher’s linear discriminant to at most $c - 1$, choosing directions that separate people relative to the variation within each person. Belhumeur et al. note that the discriminant had been applied to faces before, by Cheng et al., on features from a polar quantization of shape, and that Cui et al.’s contemporaneous Most Discriminating Feature method applied it to hand gestures; Fisherfaces apply it to the image pixels. Because the projection depends on the gallery’s classes, it can also be trained for other partitions, such as glasses against no glasses.
• View-based and modular eigenspaces (Pentland, Moghaddam and Starner, 1994): a separate eigenspace per view for pose, and eigenfeatures (eigeneyes, eigennoses, eigenmouths) for facial features, detected by their distance from feature space. On 45 people with neutral images for training and smiling ones for testing, eigenfeatures alone matched the 95% of eigenfaces and the combination reached 98%.
• Probabilistic eigenspaces (Moghaddam and Pentland, 1995): the distance from face space is treated as one marginal of a Gaussian density in image space and combined with a Mahalanobis distance in face space, the residual weighted by the mean of the discarded eigenvalues, their Eq. 11; mixtures of Gaussians handle multimodal training sets. Szeliski (2022) lists a later Bayesian extension by Moghaddam, Jebara and Pentland (2000).
• The Linear Subspace method (Belhumeur et al., 1997): three or more images per person under different lights span that person’s illumination subspace, and the probe goes to the nearest subspace.
• Successors. Later methods moved away from holistic linear subspaces: local texture descriptors, active appearance models that model shape and texture together, and then deep embeddings such as DeepFace and those trained with the margin losses of ArcFace; the face recognition article follows that path.

History

Principal component analysis is much older than its use on faces; its history from Pearson and Hotelling is in the principal component analysis article. Sirovich and Kirby (1987) applied it to face images: according to Belhumeur et al. they called the eigenvectors eigenpictures, and found a point of diminishing returns in the number needed to represent face images. Szeliski (2022) attributes to Kirby and Sirovich (1990) the observation that an arbitrary face image can be compressed and reconstructed from a mean image plus a small number of scaled signed basis images, and the term eigenpictures; Belhumeur et al. cite the 1987 paper for the term and Szeliski the 1990 paper, and since neither paper could be checked, which of the two first used it is left open. Their aim was representation and compression.

Turk and Pentland (1991) turned the representation into a recognizer. Szeliski (2022) credits them with recognizing that the coefficients of the eigenface expansion could themselves be used to construct a fast image matching algorithm. Their abstract describes a near-real-time system that locates and tracks a subject’s head and recognizes the person by comparing the eigenface weights with those of known individuals, motivated by physiology and information theory as well as by the need for speed and accuracy, and easy to implement with a neural network architecture. The conference version, “Face Recognition Using Eigenfaces”, appeared at CVPR 1991.

The work that followed addressed its failures one at a time. Pentland’s group scaled it to thousands of people and extended it to multiple views and facial features (1994) and gave the subspace distances a probabilistic interpretation for detection (1995). Belhumeur, Hespanha and Kriegman (1997) took on lighting: noting that PCA keeps the variation due to illumination, they replaced it with a class-specific projection and reported lower error rates on the Harvard and Yale databases. NIST’s FERET evaluations made a PCA algorithm the baseline against which other methods were measured, and in September 1996 an MIT Media Laboratory algorithm based on more recent work by Moghaddam, Nastar and Pentland and by Moghaddam and Pentland was the best algorithm tested for identification (Phillips et al., 1998). The face recognition article continues the story from holistic subspaces to local descriptors and deep networks.