PCA of a collection of faces

Math 123: Mathematical Aspects of Data Analysis — Fall 2026

Published

September 24, 2026

This is a sequel to the SVD of a collection of faces, run on a harder dataset. The Labeled Faces in the Wild people dataset is $13233$ news photographs of $5749$ people, most of whom appear only once.

import base64
import io

import numpy as np
import plotly.graph_objects as go
from matplotlib import colormaps
from PIL import Image
from plotly.subplots import make_subplots
from sklearn.datasets import fetch_lfw_people
from sklearn.decomposition import PCA

lfw = fetch_lfw_people()
images = lfw.images                # (13233, 62, 47), brightness in [0, 1]
A = lfw.data.astype(float)         # (13233, 2914), the same images flattened
person = lfw.target                # which of the 5749 people is in each image
names = [name.split()[-1] for name in lfw.target_names]   # surnames, for labels
n, d = A.shape
H, W = images.shape[1:]

counts = np.bincount(person)
top = np.argsort(counts)[::-1]     # people, most photographed first

print(f"{n} images of {len(names)} people, each {H} x {W} pixels")
print(f"Data matrix A has shape {n} x {d}")
print(f"{(counts == 1).sum()} people appear exactly once")
print("Most photographed:", dict(zip([names[p] for p in top[:7]], counts[top[:7]].tolist())))

13233 images of 5749 people, each 62 x 47 pixels
Data matrix A has shape 13233 x 2914
4069 people appear exactly once
Most photographed: {'Bush': 530, 'Powell': 236, 'Blair': 144, 'Rumsfeld': 121, 'Schroeder': 109, 'Sharon': 77, 'Chavez': 71}

Figure style and image helpers

# The same muted palette as the Olivetti example, which reads correctly against
# both the light and the dark version of this site.
INK = "#7c7c78"
GRID = "rgba(128, 128, 128, 0.28)"
BLUE = "#2a78d6"
GAP = 3


def style(fig, height=420, title=None, showlegend=None):
    """Apply the shared layout to a plot with axes."""
    fig.update_layout(
        template="simple_white",
        paper_bgcolor="rgba(0,0,0,0)",
        plot_bgcolor="rgba(0,0,0,0)",
        font=dict(color=INK, size=13),
        title=dict(text=title, x=0, xanchor="left", font=dict(size=15)),
        margin=dict(l=70, r=30, t=60 if title else 30, b=55),
        height=height,
        showlegend=showlegend,
        hovermode="closest",
    )
    fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    return fig


def colour(tile, cmap="gray", vmin=0.0, vmax=1.0):
    """Colour one H x W image with a matplotlib colormap, as RGBA bytes."""
    scaled = np.clip((tile - vmin) / (vmax - vmin), 0, 1)
    return (colormaps[cmap](scaled) * 255).astype(np.uint8)


def montage(rows):
    """Tile a grid of coloured H x W images into one array; gaps and None stay transparent."""
    n_rows, n_cols = len(rows), max(len(row) for row in rows)
    out = np.zeros((n_rows * (H + GAP) - GAP, n_cols * (W + GAP) - GAP, 4), np.uint8)
    for i, row in enumerate(rows):
        for j, tile in enumerate(row):
            if tile is not None:
                top, left = i * (H + GAP), j * (W + GAP)
                out[top:top + H, left:left + W] = tile
    return out


def png_uri(rgba):
    """Encode an RGBA array as a PNG data URI, which is far smaller than a heatmap."""
    buffer = io.BytesIO()
    Image.fromarray(rgba, "RGBA").save(buffer, "PNG")
    return "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode()


def tick_centres(count, size):
    return [j * (size + GAP) + (size - 1) / 2 for j in range(count)]


def image_axes(fig, n_rows, n_cols, col_labels=None, row_labels=None):
    """Label the columns and rows of a montage and hide everything else."""
    fig.update_xaxes(
        tickvals=tick_centres(n_cols, W), ticktext=col_labels or [""] * n_cols,
        showgrid=False, zeroline=False, showline=False, ticks="",
    )
    fig.update_yaxes(
        tickvals=tick_centres(n_rows, H), ticktext=row_labels or [""] * n_rows,
        showgrid=False, zeroline=False, showline=False, ticks="",
    )


def image_figure(rows, col_labels=None, row_labels=None, title=None, height=None,
                 cmaps=None):
    """Show a grid of H x W images as one picture; cmaps is one (cmap, vmin, vmax) per row."""
    n_rows, n_cols = len(rows), max(len(row) for row in rows)
    cmaps = cmaps or [("gray", 0.0, 1.0)] * n_rows
    coloured = [[None if tile is None else colour(tile, *cmap) for tile in row]
                for row, cmap in zip(rows, cmaps)]
    fig = go.Figure(go.Image(source=png_uri(montage(coloured)), hoverinfo="skip"))
    image_axes(fig, n_rows, n_cols, col_labels, row_labels)
    style(fig, height=height or 170 * n_rows + 90, title=title)
    fig.update_layout(margin=dict(l=130 if row_labels else 30, r=30,
                                  t=60 if title else 20, b=60))
    return fig

The data

Each image is $62\times47$ with brightness between $0$ and $1$, so the data matrix is $\bfA\in\Rbb^{13233\times2914}$, one flattened photograph per row. Unlike the Olivetti faces, these were not taken in a studio: lighting, background, expression, and camera angle are all uncontrolled. The people are also very unevenly represented: $4069$ of them appear exactly once, while Bush alone appears $530$ times.

Plot one face each of the most photographed people

fig = image_figure(
    [[images[person == p][0] for p in top[:7]]],
    col_labels=[f"{names[p]}<br>({counts[p]})" for p in top[:7]],
    title="Every row of A is a picture in R<sup>2914</sup>",
    height=300,
)
fig.show()

Figure 1: One photograph of each of the seven most photographed people, with the number of photographs of each.

The mean face

Recall, PCA first subtracts the sample mean $\bar{\bfx}\in\Rbb^{2914}$ from every row, giving the centered data matrix $\bfX$ whose $i$-th row is $(\bfx_i-\bar{\bfx})^\top$. The mean is itself a face, and a centered row is no longer interpretable as a greyscale image. It is positive where the photograph is brighter than average and negative where it is darker.

x_bar = A.mean(axis=0)
X = A - x_bar

print(f"Column means of X are zero up to {np.abs(X.mean(axis=0)).max():.1e}")
print(f"Entries of A lie in [{A.min():.2f}, {A.max():.2f}], "
      f"entries of X in [{X.min():.2f}, {X.max():.2f}]")

Column means of X are zero up to 1.0e-15
Entries of A lie in [0.00, 1.00], entries of X in [-0.65, 0.72]

Plot the mean face and some centered rows

shown_mean = [np.flatnonzero(person == p)[2] for p in top[[5, 3, 4, 6]]]
centered = [X[i].reshape(H, W) for i in shown_mean]
limit = float(np.abs(centered).max())
fig = image_figure(
    [[x_bar.reshape(H, W)] + [images[i] for i in shown_mean], [None] + centered],
    col_labels=["mean face"] + [names[person[i]] for i in shown_mean],
    row_labels=["photograph", "minus the mean"],
    title="Centering subtracts the mean face from every row",
    height=430,
    cmaps=[("gray", 0.0, 1.0), ("RdBu_r", -limit, limit)],
)
fig.show()

Figure 2: The mean face, four photographs (top), and the same photographs with the mean subtracted (bottom). Red is brighter than average and blue is darker, on a scale symmetric about zero.

PCA of the data matrix

By definition, PCA with $k$ components is read off from the SVD $\bfX=\sum_j\sigma_j\bfu_j\bfv_j^\top$ of the centered matrix: the principal directions are $\bfv_1,\ldots,\bfv_k\in\Rbb^{2914}$, and the scores are $\bfZ=\bfX\bfV_k=\bfU_k\Sigma_k\in\Rbb^{13233\times k}$.

Computing PCA with scikit-learn

In the Olivetti example we asked NumPy for every singular vector and then threw most of them away. Here we ask sklearn.decomposition.PCA for only the first $k=300$. For a request that small relative to $\min\{n,d\}=2914$, scikit-learn switches to a randomized SVD, which never forms the full factorization. I hope to get to this topic later in the course!

The fitted object stores the pieces of definition under its own names:

`PCA` attribute	In our notation
`mean_`	$\bar{\bfx}$
`components_`	$\bfV_k^\top$, one principal direction per row
`singular_values_`	$\sigma_1,\ldots,\sigma_k$ of the centered $\bfX$
`transform(A)`	the scores $\bfZ=\bfX\bfV_k$
`explained_variance_`	$\sigma_j^2/(n-1)$, not $\lambda_j=\sigma_j^2/n$
`explained_variance_ratio_`	$\sigma_j^2/\sum_{i}\sigma_i^2$

k_max = 300
pca = PCA(n_components=k_max, random_state=0).fit(A)
V = pca.components_.T                      # (2914, 300), the principal directions
sigma = pca.singular_values_
lam = sigma**2 / n                         # eigenvalues of C = X^T X / n
Z = pca.transform(A)                       # (13233, 300), the scores

print(f"solver used: {pca._fit_svd_solver}")
print(f"mean_ = x_bar            up to {np.abs(pca.mean_ - x_bar).max():.1e}")
print(f"transform(A) = X V_k     up to {np.abs(Z - X @ V).max():.1e}")
print(f"explained_variance_ = sigma^2/(n-1)  up to "
      f"{np.abs(pca.explained_variance_ - sigma**2 / (n - 1)).max():.1e}")

# The scores are uncorrelated, with variances lambda_j (lecture 6).
score_cov = Z.T @ Z / n
print(f"Z^T Z / n = diag(lambda)  up to {np.abs(score_cov - np.diag(lam)).max():.1e}")

solver used: randomized
mean_ = x_bar            up to 0.0e+00
transform(A) = X V_k     up to 4.7e-14
explained_variance_ = sigma^2/(n-1)  up to 0.0e+00
Z^T Z / n = diag(lambda)  up to 1.5e-04

The randomized solver is approximate, so it is worth checking how well each $\bfv_j$ satisfies the eigenvalue equation $\bfC\bfv_j=\lambda_j\bfv_j$ from PCA diagonalizes the covariance matrix.

C = X.T @ X / n
residual = np.linalg.norm(C @ V - V * lam, axis=0) / lam
for j in (1, 10, 50, 100, 207, 300):
    print(f"j = {j:>3}:  ||C v_j - lambda_j v_j|| / lambda_j = {residual[j - 1]:.1e}")

j =   1:  ||C v_j - lambda_j v_j|| / lambda_j = 7.5e-15
j =  10:  ||C v_j - lambda_j v_j|| / lambda_j = 4.7e-11
j =  50:  ||C v_j - lambda_j v_j|| / lambda_j = 7.8e-07
j = 100:  ||C v_j - lambda_j v_j|| / lambda_j = 7.1e-05
j = 207:  ||C v_j - lambda_j v_j|| / lambda_j = 1.4e-02
j = 300:  ||C v_j - lambda_j v_j|| / lambda_j = 9.9e-02

The leading directions are exact to machine precision and the last few are only rough. That trade is usually fine, since each of the last directions carries well under a tenth of a percent of the variance, as we see next.

Explained variance

print("sigma_1 through sigma_6:", sigma[:6].round(1))
print(f"sigma_1 / sigma_2 = {sigma[0] / sigma[1]:.2f}")

cosine = V[:, 0] @ x_bar / np.linalg.norm(x_bar)
print(f"cos(angle between v_1 and the mean face) = {cosine:.3f}")

sigma_1 through sigma_6: [528.7 321.1 281.3 261.8 177.2 169.7]
sigma_1 / sigma_2 = 1.65
cos(angle between v_1 and the mean face) = 0.991

The total variance $\operatorname{TV}=\|\bfX\|_F^2/n$ needs no SVD at all, which is why scikit-learn can report the explained variance ratio \[ R_k=\frac{\sum_{j\leq k}\sigma_j^2}{\sum_{j}\sigma_j^2}=\frac{\sum_{j\leq k}\lambda_j}{\operatorname{TV}} \] after computing only $300$ singular values.

total_variance = np.linalg.norm(X) ** 2 / n
R = np.cumsum(pca.explained_variance_ratio_)            # R[k - 1] is R_k

print(f"TV = ||X||_F^2 / n = {total_variance:.3f},  trace(C) = {np.trace(C):.3f}")
print(f"R_k from sklearn matches sum(lambda) / TV up to "
      f"{np.abs(R - np.cumsum(lam) / total_variance).max():.1e}")

thresholds = (0.80, 0.90, 0.95)
for t in thresholds:
    k_t = int(np.searchsorted(R, t)) + 1
    print(f"smallest k with R_k >= {t:.0%}: " + (f"{k_t}" if k_t <= k_max else f"more than {k_max}"))
print(f"R_{k_max} = {R[-1]:.1%}")

TV = ||X||_F^2 / n = 78.979,  trace(C) = 78.979
R_k from sklearn matches sum(lambda) / TV up to 1.4e-15
smallest k with R_k >= 80%: 37
smallest k with R_k >= 90%: 101
smallest k with R_k >= 95%: 207
R_300 = 96.8%

Plot the explained variance

ks = np.arange(1, k_max + 1)
fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.1,
                    subplot_titles=["Share of variance along v<sub>j</sub>",
                                    "Explained variance ratio R<sub>k</sub>"])
fig.add_trace(
    go.Scatter(x=ks, y=pca.explained_variance_ratio_, mode="lines",
               line=dict(color=BLUE, width=2),
               hovertemplate="j = %{x}<br>%{y:.2%} of the variance<extra></extra>"),
    row=1, col=1,
)
fig.add_trace(
    go.Scatter(x=ks, y=R, mode="lines", line=dict(color=BLUE, width=2),
               hovertemplate="k = %{x}<br>R<sub>k</sub> = %{y:.1%}<extra></extra>"),
    row=1, col=2,
)
for t in thresholds:
    fig.add_hline(y=t, line=dict(color=INK, width=1, dash="dot"), row=1, col=2)
    k_t = int(np.searchsorted(R, t)) + 1
    text = f"{t:.0%}: k = {k_t}" if k_t <= k_max else f"{t:.0%}: not reached"
    fig.add_annotation(x=k_max, y=t, text=text, showarrow=False, xanchor="right",
                       yanchor="bottom", font=dict(color=INK, size=12), row=1, col=2)
    if k_t <= k_max:
        fig.add_trace(
            go.Scatter(x=[k_t], y=[R[k_t - 1]], mode="markers", hoverinfo="skip",
                       marker=dict(color=BLUE, size=9, line=dict(color="white", width=2))),
            row=1, col=2,
        )
fig.update_xaxes(title="Index j", row=1, col=1)
fig.update_xaxes(title="Number of components k", row=1, col=2)
fig.update_yaxes(type="log", tickformat=".1%", row=1, col=1)
fig.update_yaxes(tickformat=".0%", range=[0, 1], row=1, col=2)
style(fig, height=420, title="How much of the spread do the first k directions keep?",
      showlegend=False)
fig.show()

Figure 3: Left: the share of the total variance along each principal direction, on a logarithmic scale. Right: the explained variance ratio R_k of the first k directions, with the 80%, 90%, and 95% levels dotted.

The usual rule of thumb, $R_k\geq95\%$, takes $k=207$ directions, out of $2914$ pixels. Had we asked for only $150$ components we could not have met it, since $R_{150}\approx93\%$. Settling for $90\%$ takes $k=101$.

What do the principal directions look like?

Each $\bfv_j\in\Rbb^{2914}$ reshapes to a $62\times47$ image, often called an eigenface.

Plot the leading principal directions

tiles = [V[:, j].reshape(H, W) for j in range(10)]
limit = float(np.abs(tiles).max())
fig = image_figure(
    [tiles[:5], tiles[5:]],
    row_labels=["v<sub>1</sub> to v<sub>5</sub>", "v<sub>6</sub> to v<sub>10</sub>"],
    title="Principal directions v<sub>1</sub> through v<sub>10</sub>",
    height=430,
    cmaps=[("RdBu_r", -limit, limit)] * 2,
)
fig.show()

Figure 4: The first ten principal directions, each reshaped to 62 × 47. Red is positive and blue is negative, on a scale symmetric about zero.

Low-rank reconstructions

The reconstruction of photograph $i$ from $k$ components is $\widehat{\bfx}_i=\bar{\bfx}+\bfV_k\bfz_i$, where $\bfz_i$ is the $i$-th row of $\bfZ$. The mean is added back on, so $k=0$ already gives every photograph the mean face, rather than the blank image the SVD would give. Storing the reconstructions costs $d$ numbers for $\bar{\bfx}$, $kd$ for $\bfV_k$, and $kn$ for the scores, and the relative error of the centered reconstruction is $\|\bfX-\bfX_k\|_F/\|\bfX\|_F=\sqrt{1-R_k}$.

def reconstruct(rows, k):
    """x_bar + V_k z_i for the photographs with the given row indices."""
    return x_bar + Z[rows, :k] @ V[:, :k].T


print(f"{'k':>4}  {'R_k':>7}  {'rel. F error':>12}  {'storage / A':>11}")
for k in (0, 1, 5, 20, 37, 101, 207, 300):
    R_k = R[k - 1] if k else 0.0
    X_k = reconstruct(np.arange(n), k) - x_bar
    error = np.linalg.norm(X - X_k) / np.linalg.norm(X)
    storage = (d + k * (n + d)) / (n * d)
    print(f"{k:>4}  {R_k:>7.3f}  {error:>12.3f}  {storage:>11.2%}")

   k      R_k  rel. F error  storage / A
   0    0.000         1.000        0.01%
   1    0.267         0.856        0.05%
   5    0.538         0.680        0.22%
  20    0.729         0.521        0.85%
  37    0.802         0.445        1.56%
 101    0.901         0.315        4.24%
 207    0.950         0.223        8.68%
 300    0.968         0.179       12.57%

The error column matches $\sqrt{1-R_k}$: at $k=101$ we keep $90\%$ of the variance, but the Frobenius error is still over $30\%$, because variance is measured in squared norms. Move the slider to watch five faces come back into focus, with $R_k$ in the title.

Build the reconstruction slider

shown_people = top[[1, 3, 0, 4, 2]]
shown = np.array([np.flatnonzero(person == p)[1] for p in shown_people])
k_values = [0, 1, 2, 5, 10, 20, 37, 50, 101, 207, 300]


def slider_title(k):
    R_k = R[k - 1] if k else 0.0
    return f"Reconstructions from k principal components, R<sub>k</sub> = {R_k:.0%}"


fig = go.Figure()
for k in k_values:
    rows = [images[shown], reconstruct(shown, k).reshape(-1, H, W)]
    fig.add_trace(
        go.Image(source=png_uri(montage([[colour(t) for t in row] for row in rows])),
                 visible=k == k_values[0], hoverinfo="skip")
    )

fig.update_layout(
    sliders=[dict(
        active=0, currentvalue=dict(prefix="k = ", font=dict(color=INK)),
        pad=dict(t=40), len=0.9, x=0.05,
        steps=[
            dict(method="update", label=str(k),
                 args=[{"visible": [k == other for other in k_values]},
                       {"title.text": slider_title(k)}])
            for k in k_values
        ],
    )]
)
image_axes(fig, 2, len(shown), [names[p] for p in shown_people], ["original", "k components"])
style(fig, height=520, title=slider_title(k_values[0]))
fig.update_layout(margin=dict(l=130, r=30, t=60, b=100))
fig.show()

Figure 5: One photograph each of five people, reconstructed from k principal components. Drag the slider to change k.

At $k=0$ every face is the mean face, and at $k=1$ each is only brighter or darker, which is what $\bfv_1$ predicted. Identity comes back slowly: the faces only start to become recognisable around $k=101$, and even at $k=207$, where the $95\%$ rule is met, they are soft, with Powell’s glasses and hand only faint outlines. With $5749$ different people in the data, a few hundred directions still cannot capture every individual face. If I were to restrict to people that appear at least a few times in the dataset, likely this would get less dramatic.

--- title: "PCA of a collection of faces" subtitle: "Math 123: Mathematical Aspects of Data Analysis — Fall 2026" date: 2026-09-24 number-sections: false jupyter: python3 format: html: toc: true toc-depth: 2 code-line-numbers: true code-tools: true out-width: "100%" column: page execute: echo: true draft: false --- This is a sequel to [the SVD of a collection of faces](../svd-olivetti-faces/svd-olivetti-faces.qmd), run on a harder dataset. The [Labeled Faces in the Wild](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.fetch_lfw_people.html) people dataset is $13233$ news photographs of $5749$ people, most of whom appear only once. ```{python} import base64 import io import numpy as np import plotly.graph_objects as go from matplotlib import colormaps from PIL import Image from plotly.subplots import make_subplots from sklearn.datasets import fetch_lfw_people from sklearn.decomposition import PCA lfw = fetch_lfw_people() images = lfw.images # (13233, 62, 47), brightness in [0, 1] A = lfw.data.astype(float) # (13233, 2914), the same images flattened person = lfw.target # which of the 5749 people is in each image names = [name.split()[-1] for name in lfw.target_names] # surnames, for labels n, d = A.shape H, W = images.shape[1:] counts = np.bincount(person) top = np.argsort(counts)[::-1] # people, most photographed first print(f"{n} images of {len(names)} people, each {H} x {W} pixels") print(f"Data matrix A has shape {n} x {d}") print(f"{(counts == 1).sum()} people appear exactly once") print("Most photographed:", dict(zip([names[p] for p in top[:7]], counts[top[:7]].tolist()))) ``` ```{python} #| code-fold: true #| code-summary: "Figure style and image helpers" # The same muted palette as the Olivetti example, which reads correctly against # both the light and the dark version of this site. INK = "#7c7c78" GRID = "rgba(128, 128, 128, 0.28)" BLUE = "#2a78d6" GAP = 3 def style(fig, height=420, title=None, showlegend=None): """Apply the shared layout to a plot with axes.""" fig.update_layout( template="simple_white", paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)", font=dict(color=INK, size=13), title=dict(text=title, x=0, xanchor="left", font=dict(size=15)), margin=dict(l=70, r=30, t=60 if title else 30, b=55), height=height, showlegend=showlegend, hovermode="closest", ) fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) return fig def colour(tile, cmap="gray", vmin=0.0, vmax=1.0): """Colour one H x W image with a matplotlib colormap, as RGBA bytes.""" scaled = np.clip((tile - vmin) / (vmax - vmin), 0, 1) return (colormaps[cmap](scaled) * 255).astype(np.uint8) def montage(rows): """Tile a grid of coloured H x W images into one array; gaps and None stay transparent.""" n_rows, n_cols = len(rows), max(len(row) for row in rows) out = np.zeros((n_rows * (H + GAP) - GAP, n_cols * (W + GAP) - GAP, 4), np.uint8) for i, row in enumerate(rows): for j, tile in enumerate(row): if tile is not None: top, left = i * (H + GAP), j * (W + GAP) out[top:top + H, left:left + W] = tile return out def png_uri(rgba): """Encode an RGBA array as a PNG data URI, which is far smaller than a heatmap.""" buffer = io.BytesIO() Image.fromarray(rgba, "RGBA").save(buffer, "PNG") return "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode() def tick_centres(count, size): return [j * (size + GAP) + (size - 1) / 2 for j in range(count)] def image_axes(fig, n_rows, n_cols, col_labels=None, row_labels=None): """Label the columns and rows of a montage and hide everything else.""" fig.update_xaxes( tickvals=tick_centres(n_cols, W), ticktext=col_labels or [""] * n_cols, showgrid=False, zeroline=False, showline=False, ticks="", ) fig.update_yaxes( tickvals=tick_centres(n_rows, H), ticktext=row_labels or [""] * n_rows, showgrid=False, zeroline=False, showline=False, ticks="", ) def image_figure(rows, col_labels=None, row_labels=None, title=None, height=None, cmaps=None): """Show a grid of H x W images as one picture; cmaps is one (cmap, vmin, vmax) per row.""" n_rows, n_cols = len(rows), max(len(row) for row in rows) cmaps = cmaps or [("gray", 0.0, 1.0)] * n_rows coloured = [[None if tile is None else colour(tile, *cmap) for tile in row] for row, cmap in zip(rows, cmaps)] fig = go.Figure(go.Image(source=png_uri(montage(coloured)), hoverinfo="skip")) image_axes(fig, n_rows, n_cols, col_labels, row_labels) style(fig, height=height or 170 * n_rows + 90, title=title) fig.update_layout(margin=dict(l=130 if row_labels else 30, r=30, t=60 if title else 20, b=60)) return fig ``` ## The data {#sec-data} Each image is $62\times47$ with brightness between $0$ and $1$, so the [data matrix](../../lectures/lecture-03.qmd#def-data-matrix-review) is $\bfA\in\Rbb^{13233\times2914}$, one flattened photograph per row. Unlike the Olivetti faces, these were not taken in a studio: lighting, background, expression, and camera angle are all uncontrolled. The people are also very unevenly represented: $4069$ of them appear exactly once, while Bush alone appears $530$ times. ```{python} #| label: fig-face-gallery #| fig-cap: "One photograph of each of the seven most photographed people, with the number of photographs of each." #| fig-alt: "A single row of seven small greyscale face photographs, one per person, labelled with surnames and photo counts." #| code-fold: true #| code-summary: "Plot one face each of the most photographed people" fig = image_figure( [[images[person == p][0] for p in top[:7]]], col_labels=[f"{names[p]} ({counts[p]})" for p in top[:7]], title="Every row of A is a picture in R2914", height=300, ) fig.show() ``` ## The mean face {#sec-mean-face} Recall, PCA first subtracts the sample mean $\bar{\bfx}\in\Rbb^{2914}$ from every row, giving the [centered data matrix](../../lectures/lecture-05.qmd#def-centered-data) $\bfX$ whose $i$-th row is $(\bfx_i-\bar{\bfx})^\top$. The mean is itself a face, and a centered row is no longer interpretable as a greyscale image. It is positive where the photograph is brighter than average and negative where it is darker. ```{python} x_bar = A.mean(axis=0) X = A - x_bar print(f"Column means of X are zero up to {np.abs(X.mean(axis=0)).max():.1e}") print(f"Entries of A lie in [{A.min():.2f}, {A.max():.2f}], " f"entries of X in [{X.min():.2f}, {X.max():.2f}]") ``` ```{python} #| label: fig-mean-face #| fig-cap: "The mean face, four photographs (top), and the same photographs with the mean subtracted (bottom). Red is brighter than average and blue is darker, on a scale symmetric about zero." #| fig-alt: "Top row: a blurry average face followed by four photographs. Bottom row: red-and-blue difference images of the four photographs, showing mostly lighting and background differences." #| code-fold: true #| code-summary: "Plot the mean face and some centered rows" shown_mean = [np.flatnonzero(person == p)[2] for p in top[[5, 3, 4, 6]]] centered = [X[i].reshape(H, W) for i in shown_mean] limit = float(np.abs(centered).max()) fig = image_figure( [[x_bar.reshape(H, W)] + [images[i] for i in shown_mean], [None] + centered], col_labels=["mean face"] + [names[person[i]] for i in shown_mean], row_labels=["photograph", "minus the mean"], title="Centering subtracts the mean face from every row", height=430, cmaps=[("gray", 0.0, 1.0), ("RdBu_r", -limit, limit)], ) fig.show() ``` ## PCA of the data matrix {#sec-pca} By [definition](../../lectures/lecture-06.qmd#def-pca), PCA with $k$ components is read off from the SVD $\bfX=\sum_j\sigma_j\bfu_j\bfv_j^\top$ of the centered matrix: the principal directions are $\bfv_1,\ldots,\bfv_k\in\Rbb^{2914}$, and the scores are $\bfZ=\bfX\bfV_k=\bfU_k\Sigma_k\in\Rbb^{13233\times k}$. ### Computing PCA with scikit-learn In the Olivetti example we asked NumPy for every singular vector and then threw most of them away. Here we ask `sklearn.decomposition.PCA` for only the first $k=300$. For a request that small relative to $\min\{n,d\}=2914$, scikit-learn switches to a *randomized* SVD, which never forms the full factorization. I hope to get to this topic later in the course! The fitted object stores the pieces of [definition](../../lectures/lecture-06.qmd#def-pca) under its own names: | `PCA` attribute | In our notation | |---|---| | `mean_` | $\bar{\bfx}$ | | `components_` | $\bfV_k^\top$, one principal direction per row | | `singular_values_` | $\sigma_1,\ldots,\sigma_k$ of the centered $\bfX$ | | `transform(A)` | the scores $\bfZ=\bfX\bfV_k$ | | `explained_variance_` | $\sigma_j^2/(n-1)$, **not** $\lambda_j=\sigma_j^2/n$ | | `explained_variance_ratio_` | $\sigma_j^2/\sum_{i}\sigma_i^2$ | ```{python} k_max = 300 pca = PCA(n_components=k_max, random_state=0).fit(A) V = pca.components_.T # (2914, 300), the principal directions sigma = pca.singular_values_ lam = sigma**2 / n # eigenvalues of C = X^T X / n Z = pca.transform(A) # (13233, 300), the scores print(f"solver used: {pca._fit_svd_solver}") print(f"mean_ = x_bar up to {np.abs(pca.mean_ - x_bar).max():.1e}") print(f"transform(A) = X V_k up to {np.abs(Z - X @ V).max():.1e}") print(f"explained_variance_ = sigma^2/(n-1) up to " f"{np.abs(pca.explained_variance_ - sigma**2 / (n - 1)).max():.1e}") # The scores are uncorrelated, with variances lambda_j (lecture 6). score_cov = Z.T @ Z / n print(f"Z^T Z / n = diag(lambda) up to {np.abs(score_cov - np.diag(lam)).max():.1e}") ``` The randomized solver is approximate, so it is worth checking how well each $\bfv_j$ satisfies the eigenvalue equation $\bfC\bfv_j=\lambda_j\bfv_j$ from [PCA diagonalizes the covariance matrix](../../lectures/lecture-06.qmd#prp-covariance-and-pca). ```{python} C = X.T @ X / n residual = np.linalg.norm(C @ V - V * lam, axis=0) / lam for j in (1, 10, 50, 100, 207, 300): print(f"j = {j:>3}: ||C v_j - lambda_j v_j|| / lambda_j = {residual[j - 1]:.1e}") ``` The leading directions are exact to machine precision and the last few are only rough. That trade is usually fine, since each of the last directions carries well under a tenth of a percent of the variance, as we see next. ### Explained variance ```{python} print("sigma_1 through sigma_6:", sigma[:6].round(1)) print(f"sigma_1 / sigma_2 = {sigma[0] / sigma[1]:.2f}") cosine = V[:, 0] @ x_bar / np.linalg.norm(x_bar) print(f"cos(angle between v_1 and the mean face) = {cosine:.3f}") ``` The [total variance](../../lectures/lecture-06.qmd#def-total-variance) $\operatorname{TV}=\|\bfX\|_F^2/n$ needs no SVD at all, which is why scikit-learn can report the [explained variance ratio](../../lectures/lecture-06.qmd#def-explained-variance) $$ R_k=\frac{\sum_{j\leq k}\sigma_j^2}{\sum_{j}\sigma_j^2}=\frac{\sum_{j\leq k}\lambda_j}{\operatorname{TV}} $$ after computing only $300$ singular values. ```{python} total_variance = np.linalg.norm(X) ** 2 / n R = np.cumsum(pca.explained_variance_ratio_) # R[k - 1] is R_k print(f"TV = ||X||_F^2 / n = {total_variance:.3f}, trace(C) = {np.trace(C):.3f}") print(f"R_k from sklearn matches sum(lambda) / TV up to " f"{np.abs(R - np.cumsum(lam) / total_variance).max():.1e}") thresholds = (0.80, 0.90, 0.95) for t in thresholds: k_t = int(np.searchsorted(R, t)) + 1 print(f"smallest k with R_k >= {t:.0%}: " + (f"{k_t}" if k_t <= k_max else f"more than {k_max}")) print(f"R_{k_max} = {R[-1]:.1%}") ``` ```{python} #| label: fig-explained-variance #| fig-cap: "Left: the share of the total variance along each principal direction, on a logarithmic scale. Right: the explained variance ratio R_k of the first k directions, with the 80%, 90%, and 95% levels dotted." #| fig-alt: "Two line charts against k from one to three hundred. The left curve falls from about twenty-seven percent to about a hundredth of a percent. The right curve rises from twenty-seven percent to about ninety-seven percent, crossing eighty percent at k equals thirty-seven, ninety percent at k equals one hundred one, and ninety-five percent at k equals two hundred seven." #| code-fold: true #| code-summary: "Plot the explained variance" ks = np.arange(1, k_max + 1) fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.1, subplot_titles=["Share of variance along vj", "Explained variance ratio Rk"]) fig.add_trace( go.Scatter(x=ks, y=pca.explained_variance_ratio_, mode="lines", line=dict(color=BLUE, width=2), hovertemplate="j = %{x} %{y:.2%} of the variance<extra></extra>"), row=1, col=1, ) fig.add_trace( go.Scatter(x=ks, y=R, mode="lines", line=dict(color=BLUE, width=2), hovertemplate="k = %{x} Rk = %{y:.1%}<extra></extra>"), row=1, col=2, ) for t in thresholds: fig.add_hline(y=t, line=dict(color=INK, width=1, dash="dot"), row=1, col=2) k_t = int(np.searchsorted(R, t)) + 1 text = f"{t:.0%}: k = {k_t}" if k_t <= k_max else f"{t:.0%}: not reached" fig.add_annotation(x=k_max, y=t, text=text, showarrow=False, xanchor="right", yanchor="bottom", font=dict(color=INK, size=12), row=1, col=2) if k_t <= k_max: fig.add_trace( go.Scatter(x=[k_t], y=[R[k_t - 1]], mode="markers", hoverinfo="skip", marker=dict(color=BLUE, size=9, line=dict(color="white", width=2))), row=1, col=2, ) fig.update_xaxes(title="Index j", row=1, col=1) fig.update_xaxes(title="Number of components k", row=1, col=2) fig.update_yaxes(type="log", tickformat=".1%", row=1, col=1) fig.update_yaxes(tickformat=".0%", range=[0, 1], row=1, col=2) style(fig, height=420, title="How much of the spread do the first k directions keep?", showlegend=False) fig.show() ``` The usual rule of thumb, $R_k\geq95\%$, takes $k=207$ directions, out of $2914$ pixels. Had we asked for only $150$ components we could not have met it, since $R_{150}\approx93\%$. Settling for $90\%$ takes $k=101$. ### What do the principal directions look like? Each $\bfv_j\in\Rbb^{2914}$ reshapes to a $62\times47$ image, often called an *eigenface*. ```{python} #| label: fig-principal-directions #| fig-cap: "The first ten principal directions, each reshaped to 62 × 47. Red is positive and blue is negative, on a scale symmetric about zero." #| fig-alt: "Ten red-and-blue face-shaped patterns. The first is almost uniformly red like a mean face; the next contrast left against right lighting, face against background, and highlight the eyes, mouth, and hairline." #| code-fold: true #| code-summary: "Plot the leading principal directions" tiles = [V[:, j].reshape(H, W) for j in range(10)] limit = float(np.abs(tiles).max()) fig = image_figure( [tiles[:5], tiles[5:]], row_labels=["v1 to v5", "v6 to v10"], title="Principal directions v1 through v10", height=430, cmaps=[("RdBu_r", -limit, limit)] * 2, ) fig.show() ``` ## Low-rank reconstructions {#sec-reconstruct} The [reconstruction](../../lectures/lecture-06.qmd#def-pca) of photograph $i$ from $k$ components is $\widehat{\bfx}_i=\bar{\bfx}+\bfV_k\bfz_i$, where $\bfz_i$ is the $i$-th row of $\bfZ$. The mean is added back on, so $k=0$ already gives every photograph the mean face, rather than the blank image the SVD would give. Storing the reconstructions costs $d$ numbers for $\bar{\bfx}$, $kd$ for $\bfV_k$, and $kn$ for the scores, and the relative error of the centered reconstruction is $\|\bfX-\bfX_k\|_F/\|\bfX\|_F=\sqrt{1-R_k}$. ```{python} def reconstruct(rows, k): """x_bar + V_k z_i for the photographs with the given row indices.""" return x_bar + Z[rows, :k] @ V[:, :k].T print(f"{'k':>4} {'R_k':>7} {'rel. F error':>12} {'storage / A':>11}") for k in (0, 1, 5, 20, 37, 101, 207, 300): R_k = R[k - 1] if k else 0.0 X_k = reconstruct(np.arange(n), k) - x_bar error = np.linalg.norm(X - X_k) / np.linalg.norm(X) storage = (d + k * (n + d)) / (n * d) print(f"{k:>4} {R_k:>7.3f} {error:>12.3f} {storage:>11.2%}") ``` The error column matches $\sqrt{1-R_k}$: at $k=101$ we keep $90\%$ of the variance, but the Frobenius error is still over $30\%$, because variance is measured in squared norms. Move the slider to watch five faces come back into focus, with $R_k$ in the title. ```{python} #| label: fig-reconstruction-slider #| fig-cap: "One photograph each of five people, reconstructed from k principal components. Drag the slider to change k." #| fig-alt: "An interactive grid of five original faces above their reconstructions, with a slider controlling k from zero, where every reconstruction is the mean face, to three hundred." #| code-fold: true #| code-summary: "Build the reconstruction slider" shown_people = top[[1, 3, 0, 4, 2]] shown = np.array([np.flatnonzero(person == p)[1] for p in shown_people]) k_values = [0, 1, 2, 5, 10, 20, 37, 50, 101, 207, 300] def slider_title(k): R_k = R[k - 1] if k else 0.0 return f"Reconstructions from k principal components, Rk = {R_k:.0%}" fig = go.Figure() for k in k_values: rows = [images[shown], reconstruct(shown, k).reshape(-1, H, W)] fig.add_trace( go.Image(source=png_uri(montage([[colour(t) for t in row] for row in rows])), visible=k == k_values[0], hoverinfo="skip") ) fig.update_layout( sliders=[dict( active=0, currentvalue=dict(prefix="k = ", font=dict(color=INK)), pad=dict(t=40), len=0.9, x=0.05, steps=[ dict(method="update", label=str(k), args=[{"visible": [k == other for other in k_values]}, {"title.text": slider_title(k)}]) for k in k_values ], )] ) image_axes(fig, 2, len(shown), [names[p] for p in shown_people], ["original", "k components"]) style(fig, height=520, title=slider_title(k_values[0])) fig.update_layout(margin=dict(l=130, r=30, t=60, b=100)) fig.show() ``` At $k=0$ every face is the mean face, and at $k=1$ each is only brighter or darker, which is what $\bfv_1$ predicted. Identity comes back slowly: the faces only start to become recognisable around $k=101$, and even at $k=207$, where the $95\%$ rule is met, they are soft, with Powell's glasses and hand only faint outlines. With $5749$ different people in the data, a few hundred directions still cannot capture every individual face. If I were to restrict to people that appear at least a few times in the dataset, likely this would get less dramatic.