The SVD of a collection of faces

Math 123: Mathematical Aspects of Data Analysis — Fall 2026

Published

September 17, 2026

The Olivetti faces are $400$ greyscale photographs, ten each of $40$ people. We will compress one face on its own first, and then take the SVD of the full $400\times4096$ data matrix.

import base64
import io

import numpy as np
import plotly.graph_objects as go
from matplotlib import colormaps
from PIL import Image
from sklearn.datasets import fetch_olivetti_faces

faces = fetch_olivetti_faces()     # downloads about 4 MB on first use
images = faces.images              # (400, 64, 64), brightness in [0, 1]
A = faces.data.astype(float)       # (400, 4096), the same images flattened
person = faces.target              # which of the 40 people is in each image
m, n = A.shape

print(f"{m} images of {person.max() + 1} people, each 64 x 64 pixels")
print(f"Data matrix A has shape {m} x {n}")
print(f"Brightness runs from {A.min():.2f} to {A.max():.2f}")

400 images of 40 people, each 64 x 64 pixels
Data matrix A has shape 400 x 4096
Brightness runs from 0.00 to 1.00

Figure style and image helpers

# A muted ink colour and a translucent grid read correctly against both the
# light and the dark version of this site, so the figures need no theme switch.
INK = "#7c7c78"
GRID = "rgba(128, 128, 128, 0.28)"
BLUE, ORANGE, AQUA, RED = "#2a78d6", "#eb6834", "#1baf7a", "#e34948"
SIDE = 64


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 montage(rows, gap=4):
    """Tile a grid of 64x64 images into one array, separated by blank (NaN) pixels."""
    n_rows, n_cols = len(rows), max(len(row) for row in rows)
    out = np.full((n_rows * (SIDE + gap) - gap, n_cols * (SIDE + gap) - gap), np.nan)
    for i, row in enumerate(rows):
        for j, tile in enumerate(row):
            top, left = i * (SIDE + gap), j * (SIDE + gap)
            out[top:top + SIDE, left:left + SIDE] = tile
    return out


def png_uri(array, cmap="gray", vmin=0.0, vmax=1.0):
    """Encode an array as a PNG data URI; NaN pixels become transparent.

    A 64x64 heatmap costs thousands of numbers in the page, while the same
    picture as a PNG is a few kilobytes, so every face below is drawn this way.
    """
    scaled = np.clip((array - vmin) / (vmax - vmin), 0, 1)
    rgba = (colormaps[cmap](np.nan_to_num(scaled)) * 255).astype(np.uint8)
    rgba[np.isnan(array), 3] = 0
    buffer = io.BytesIO()
    Image.fromarray(rgba, "RGBA").save(buffer, "PNG")
    return "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode()


def tick_centres(count, gap=4):
    return [j * (SIDE + gap) + (SIDE - 1) / 2 for j in range(count)]


def image_figure(rows, col_labels=None, row_labels=None, title=None, height=None,
                 cmap="gray", vmin=0.0, vmax=1.0, gap=4):
    """Show a grid of 64x64 images as a single picture."""
    n_rows, n_cols = len(rows), max(len(row) for row in rows)
    fig = go.Figure(
        go.Image(source=png_uri(montage(rows, gap), cmap, vmin, vmax), hoverinfo="skip")
    )
    fig.update_xaxes(
        tickvals=tick_centres(n_cols, gap), ticktext=col_labels or [""] * n_cols,
        showgrid=False, zeroline=False, showline=False, ticks="",
    )
    fig.update_yaxes(
        tickvals=tick_centres(n_rows, gap), ticktext=row_labels or [""] * n_rows,
        showgrid=False, zeroline=False, showline=False, ticks="",
    )
    style(fig, height=height or 130 * n_rows + 90, title=title)
    fig.update_layout(margin=dict(l=110 if row_labels else 30, r=30,
                                  t=60 if title else 20, b=50))
    return fig


def rank_k(U, s, Vt, k):
    """The truncated SVD A_k = U_k Sigma_k V_k^T."""
    return (U[:, :k] * s[:k]) @ Vt[:k, :]


def relative_frobenius_tail(s):
    """Entry k is ||A - A_k||_F / ||A||_F, computed from the singular values alone."""
    tail = np.concatenate([np.cumsum(s[::-1] ** 2)[::-1], [0.0]])
    return np.sqrt(tail / tail[0])

The data

Each image is $64\times64$ with brightness between $0$ (black) and $1$ (white), so flattening one gives a vector in $\mathbb{R}^{4096}$. Stacking the $400$ flattened images as rows gives the data matrix \[ A= \begin{pmatrix} - & \mathbf{x}_1^\top & - \\ & \vdots & \\ - & \mathbf{x}_{400}^\top & - \end{pmatrix} \in\mathbb{R}^{400\times4096}. \] The photographs of each person vary in pose, expression, lighting, and glasses, but they share a common framing: every face is roughly centred and the same size.

Plot a gallery of faces

people = [0, 5, 10, 14, 20, 26, 29, 35]
fig = image_figure(
    [[images[person == p][pose] for p in people] for pose in range(4)],
    col_labels=[f"person {p}" for p in people],
    title="Every row of A is a picture in R<sup>4096</sup>",
    height=560,
)
fig.show()

Figure 1: Four of the ten photographs of eight of the forty people. Each column is one person.

One face as a matrix

Before touching the collection, take one photograph on its own as a matrix $F\in\mathbb{R}^{64\times64}$ with compact SVD $F=\sum_{i=1}^r \sigma_i \mathbf{u}_i\mathbf{v}_i^\top$.

F = images[0]
U_f, s_f, Vt_f = np.linalg.svd(F, full_matrices=False)

print("Rank of F              :", np.linalg.matrix_rank(F))
print("Largest singular values:", s_f[:6].round(2))
print("Smallest singular value:", f"{s_f[-1]:.1e}")

Rank of F              : 63
Largest singular values: [41.33  5.22  2.6   1.94  1.22  0.99]
Smallest singular value: 2.6e-04

A photograph is essentially never exactly low rank, but its singular values fall by an order of magnitude almost immediately.

Plot the singular values of one face

fig = go.Figure(
    go.Scatter(
        x=np.arange(1, len(s_f) + 1), y=s_f, mode="lines+markers",
        line=dict(color=BLUE, width=2), marker=dict(size=6),
        hovertemplate="i = %{x}<br>sigma = %{y:.3g}<extra></extra>",
    )
)
fig.update_xaxes(title="Index i", dtick=8)
fig.update_yaxes(title="Singular value", type="log", exponentformat="power")
style(fig, height=360, title="Singular values of a single face")
fig.show()

Figure 2: Singular values of one 64 × 64 photograph, on a logarithmic scale.

The truncated SVD $F_k=\sum_{i=1}^k\sigma_i\mathbf{u}_i\mathbf{v}_i^\top$ is the best rank-$k$ approximation of $F$ in both norms, and by the spectral and Frobenius Eckart–Young theorems its relative errors are read off from the singular values: \[ \frac{\|F-F_k\|_2}{\|F\|_2} = \frac{\sigma_{k+1}}{\sigma_1}, \qquad \frac{\|F-F_k\|_F}{\|F\|_F} = \sqrt{\frac{\sum_{i>k}\sigma_i^2}{\sum_{i\geq1}\sigma_i^2}} . \]

Plot rank-k approximations of one face

ranks = [1, 2, 5, 10, 20, 30]
single_spectral_err = s_f / s_f[0]                 # entry k is sigma_{k+1} / sigma_1
single_frobenius_err = relative_frobenius_tail(s_f)
fig = image_figure(
    [[F] + [rank_k(U_f, s_f, Vt_f, k) for k in ranks]],
    col_labels=["original"] + [
        f"k = {k}<br>2-norm: {single_spectral_err[k]:.1%}"
        f"<br>F-norm: {single_frobenius_err[k]:.1%}"
        for k in ranks
    ],
    title="Rank-k approximations of one 64 x 64 face",
    height=290,
)
fig.update_layout(margin=dict(b=80))
fig.show()

Figure 3: Rank-k approximations of one photograph, with the relative spectral and Frobenius errors beneath each rank.

The SVD of the data matrix

The rows of $A\in\mathbb{R}^{400\times4096}$ are whole faces, and a column is one pixel position measured across all $400$ photographs. Its SVD \[ A=\sum_{i=1}^r\sigma_i\mathbf{u}_i\mathbf{v}_i^\top, \qquad \mathbf{u}_i\in\mathbb{R}^{400}, \quad \mathbf{v}_i\in\mathbb{R}^{4096}, \] has right singular vectors that are themselves $64\times64$ images, and left singular vectors that assign a number to each of the $400$ photographs.

Computing the economy SVD

NumPy computes the SVD directly with np.linalg.svd. By default it returns the full SVD, whose $U$ is a square $400\times400$ orthogonal matrix and whose $V$ is a square $4096\times4096$ orthogonal matrix. That $V$ alone has almost $17$ million entries, and most of them are wasted. Since $\Sigma\in\mathbb{R}^{400\times4096}$ has only $p=\min\{400,4096\}=400$ diagonal entries, the last $4096-400=3696$ rows of $V^\top$ are multiplied by columns of zeros and never contribute to $A$.

Passing full_matrices=False asks for the economy SVD instead, which keeps only the first $p$ singular vectors on each side: \[ A=U_p\Sigma_pV_p^\top, \qquad U_p\in\mathbb{R}^{400\times400}, \quad \Sigma_p\in\mathbb{R}^{400\times400}, \quad V_p\in\mathbb{R}^{4096\times400}. \] Both $U_p$ and $V_p$ still have orthonormal columns, but now $V_p$ is tall and $V_pV_p^\top\neq I$. If $\operatorname{rank}(A)=r<p$, the last $p-r$ diagonal entries of $\Sigma_p$ are zero, and dropping them gives the compact SVD $U_r\Sigma_rV_r^\top$ from lecture.

U, sigma, Vt = np.linalg.svd(A, full_matrices=False)
rank = np.linalg.matrix_rank(A)

# Singular vectors are only defined up to sign. Fix one choice: make the entry
# of largest magnitude in each v_i positive, and flip u_i to match, so the
# product sigma_i u_i v_i^T is unchanged.
signs = np.sign(Vt[np.arange(len(sigma)), np.abs(Vt).argmax(axis=1)])
U, Vt = U * signs, Vt * signs[:, None]

print(f"U: {U.shape},  sigma: {sigma.shape},  Vt: {Vt.shape}")
print(f"rank(A) = {rank}, smallest singular value = {sigma[-1]:.3f}")
print(f"U^T U = I   up to {np.abs(U.T @ U - np.eye(len(sigma))).max():.1e}")
print(f"V^T V = I   up to {np.abs(Vt @ Vt.T - np.eye(len(sigma))).max():.1e}")
print(f"||A - U diag(sigma) V^T||_F / ||A||_F = "
      f"{np.linalg.norm(A - (U * sigma) @ Vt) / np.linalg.norm(A):.1e}")

U: (400, 400),  sigma: (400,),  Vt: (400, 4096)
rank(A) = 400, smallest singular value = 0.719
U^T U = I   up to 3.1e-15
V^T V = I   up to 2.9e-15
||A - U diag(sigma) V^T||_F / ||A||_F = 8.7e-15

The collection has full rank $r=p=400$: no photograph is exactly a combination of the others, so every singular value is positive and the economy and compact SVDs coincide.

How fast do the singular values decay?

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

mean_face = A.mean(axis=0)
cosine = Vt[0] @ mean_face / np.linalg.norm(mean_face)
print(f"cos(angle between v_1 and the average face) = {cosine:.5f}")

sigma_1 through sigma_6: [717.2  67.5  50.5  40.2  36.3  31.6]
sigma_1 / sigma_2 = 10.6
cos(angle between v_1 and the average face) = 0.99999

The first singular value is more than ten times the second. Every row of $A$ is a nonnegative image of a centred face, so all $400$ rows point in roughly the same direction in $\mathbb{R}^{4096}$. The right singular vector $\mathbf{v}_1$ is that shared direction, and it is essentially the average face rescaled to unit length.

Plot the singular values of the data matrix

fig = go.Figure(
    go.Scatter(
        x=np.arange(1, rank + 1), y=sigma, mode="lines",
        line=dict(color=BLUE, width=2),
        hovertemplate="i = %{x}<br>sigma = %{y:.3g}<extra></extra>",
    )
)
fig.add_annotation(x=1, y=np.log10(sigma[0]), text="  sigma<sub>1</sub>: the average face",
                   showarrow=False, xanchor="left", font=dict(color=BLUE, size=13))
fig.update_xaxes(title="Index i", dtick=50)
fig.update_yaxes(title="Singular value", type="log", exponentformat="power")
style(fig, height=380, title="Singular values of the 400 x 4096 face matrix")
fig.show()

Figure 4: All 400 singular values of the face data matrix, on a logarithmic scale. The first stands far above the rest.

Eckart–Young, numerically

The Eckart–Young theorem makes two promises about $A_k=U_k\Sigma_kV_k^\top$: the error is exactly the discarded singular values, \[ \|A-A_k\|_2=\sigma_{k+1}, \qquad \|A-A_k\|_F=\Big(\sum_{i>k}\sigma_i^2\Big)^{1/2}, \] and no other matrix of rank at most $k$ does better. The first promise we can check directly. For the second, we cannot try every rank-$k$ matrix, but we can try a natural competitor: pick $k$ of the photographs at random and replace every face by its closest point in their span. That matrix $B$ also has rank at most $k$.

rng = np.random.default_rng(0)


def random_faces_approximation(k):
    """Project every row of A onto the span of k randomly chosen rows."""
    Q, _ = np.linalg.qr(A[rng.choice(m, size=k, replace=False)].T)   # (4096, k)
    return (A @ Q) @ Q.T


frob_tail = relative_frobenius_tail(sigma) * np.linalg.norm(A)
print(f"{'k':>4}  {'||A-A_k||_2':>12}  {'sigma_k+1':>10}  {'||A-A_k||_F':>12}  "
      f"{'tail':>8}  {'||A-B||_2':>10}  {'||A-B||_F':>10}")
for k in (1, 5, 20, 50, 100):
    Ak = rank_k(U, sigma, Vt, k)
    B = random_faces_approximation(k)
    print(f"{k:>4}  {np.linalg.norm(A - Ak, 2):>12.3f}  {sigma[k]:>10.3f}  "
          f"{np.linalg.norm(A - Ak):>12.3f}  {frob_tail[k]:>8.3f}  "
          f"{np.linalg.norm(A - B, 2):>10.3f}  {np.linalg.norm(A - B):>10.3f}")

   k   ||A-A_k||_2   sigma_k+1   ||A-A_k||_F      tail   ||A-B||_2   ||A-B||_F
   1        67.454      67.454       157.019   157.019     199.348     252.229
   5        31.587      31.587       120.906   120.906      62.516     151.120
  20        15.007      15.007        86.803    86.803      29.401     113.630
  50         7.903       7.903        63.308    63.308      20.913      90.649
 100         4.935       4.935        45.293    45.293      13.146      68.086

The first two pairs of columns match to every printed digit, and the random-faces competitor is worse in both norms at every $k$.

Plot error against rank for both norms

from plotly.subplots import make_subplots

k_grid = [1, 2, 3, 5, 8, 12, 20, 30, 50, 75, 100, 150, 200, 300]
competitor = [A - random_faces_approximation(k) for k in k_grid]
curves = {
    "spectral": (sigma[k_grid] / sigma[0],
                 [np.linalg.norm(E, 2) / sigma[0] for E in competitor]),
    "Frobenius": (frob_tail[k_grid] / frob_tail[0],
                  [np.linalg.norm(E) / frob_tail[0] for E in competitor]),
}

fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.1,
                    subplot_titles=["||A - B||<sub>2</sub> / ||A||<sub>2</sub>",
                                    "||A - B||<sub>F</sub> / ||A||<sub>F</sub>"])
for col, (norm, (svd_curve, other_curve)) in enumerate(curves.items(), start=1):
    for label, curve, colour in (("truncated SVD A_k", svd_curve, BLUE),
                                 ("k random photographs", other_curve, ORANGE)):
        fig.add_trace(
            go.Scatter(
                x=k_grid, y=curve, mode="lines+markers", name=label,
                line=dict(color=colour, width=2), marker=dict(size=6),
                showlegend=col == 1,
                hovertemplate=f"{label}<br>k = %{{x}}<br>{norm} error %{{y:.1%}}<extra></extra>",
            ),
            row=1, col=col,
        )
fig.update_xaxes(title="Rank k", type="log")
fig.update_yaxes(tickformat=".0%", rangemode="tozero")
style(fig, height=420, title="The truncated SVD beats a natural competitor in both norms")
fig.update_layout(legend=dict(orientation="h", y=-0.22, x=0))
fig.update_layout(margin=dict(b=90))
fig.show()

Figure 5: Relative error of the truncated SVD and of projecting onto k random photographs, in the spectral norm (left) and the Frobenius norm (right). The SVD curves come from the singular values alone.

Relative to $\|A\|_2=\sigma_1$ the spectral error collapses to about $10\%$ after a single term, because $\sigma_1$ is so dominant. The Frobenius error falls much more slowly: no single missing piece is large, but there are many of them.

What do the right singular vectors look like?

Each $\mathbf{v}_i\in\mathbb{R}^{4096}$ reshapes to a $64\times64$ image.

Plot the leading right singular vectors

tiles = [Vt[i].reshape(SIDE, SIDE) for i in range(12)]
limit = float(np.abs(tiles).max())
fig = image_figure(
    [tiles[:6], tiles[6:]],
    col_labels=[""] * 6,
    title="Right singular vectors v<sub>1</sub> through v<sub>12</sub>",
    height=400, cmap="RdBu_r", vmin=-limit, vmax=limit,
)
for row in range(2):
    for col in range(6):
        fig.add_annotation(
            x=tick_centres(6)[col], y=row * (SIDE + 4) - 5,
            text=f"v<sub>{row * 6 + col + 1}</sub>", showarrow=False,
            font=dict(color=INK, size=13),
        )
fig.show()

Figure 6: The first twelve right singular vectors, each reshaped to 64 × 64. Red is positive and blue is negative, on a scale symmetric about zero.

Low-rank reconstructions

Storing $A_k$ costs $k(400+4096+1)$ numbers, against $400\times4096$ for $A$ itself. The plot below sets that cost against both relative errors from Eckart–Young, $\|A-A_k\|_2/\|A\|_2=\sigma_{k+1}/\sigma_1$ and $\|A-A_k\|_F/\|A\|_F$, at every rank shown in the figures that follow.

Plot error and storage against rank

k_values = [1, 2, 3, 5, 10, 20, 50, 100, 200, 300, 400]   # every rank shown below
ks = np.array(k_values)

spectral_error = np.append(sigma, 0.0)[ks] / sigma[0]    # sigma_{k+1} / sigma_1
frobenius_error = frob_tail[ks] / frob_tail[0]
storage = ks * (m + n + 1) / (m * n)

fig = make_subplots(specs=[[{"secondary_y": True}]])
for label, curve, colour, secondary in (
    ("relative spectral error", spectral_error, AQUA, False),
    ("relative Frobenius error", frobenius_error, BLUE, False),
    ("storage (fraction of A)", storage, ORANGE, True),
):
    fig.add_trace(
        go.Scatter(
            x=ks, y=curve, mode="lines+markers", name=label,
            line=dict(color=colour, width=2, dash="dash" if secondary else "solid"),
            marker=dict(size=7),
            hovertemplate=f"{label}<br>k = %{{x}}<br>%{{y:.1%}}<extra></extra>",
        ),
        secondary_y=secondary,
    )
fig.add_hline(y=1, line=dict(color=ORANGE, width=1, dash="dot"), secondary_y=True)

fig.update_xaxes(title="Rank k", type="log", tickvals=k_values)
fig.update_yaxes(title="Relative error", tickformat=".0%", rangemode="tozero",
                 secondary_y=False)
fig.update_yaxes(title="Storage of A<sub>k</sub> / storage of A", tickformat=".0%",
                 rangemode="tozero", showgrid=False, title_font=dict(color=ORANGE),
                 tickfont=dict(color=ORANGE), secondary_y=True)
style(fig, height=440, title="What each rank costs, and what it buys", showlegend=True)
fig.update_layout(legend=dict(orientation="h", y=-0.2, x=0),
                  margin=dict(l=70, r=80, t=60, b=100))
fig.show()

Figure 7: Relative spectral and Frobenius error of A_k (left axis) and the storage cost of A_k as a fraction of the raw pixels (right axis), at every rank shown below. The dotted line marks where A_k costs as much as A.

Move the slider to watch six faces come back into focus.

Build the reconstruction slider

shown_people = [6, 10, 14, 26, 29, 35]
shown = np.array([np.flatnonzero(person == p)[0] for p in shown_people])
gap = 4

fig = go.Figure()
for k in k_values:
    reconstruction = rank_k(U[shown], sigma, Vt, k).reshape(-1, SIDE, SIDE)
    fig.add_trace(
        go.Image(source=png_uri(montage([images[shown], reconstruction], gap)),
                 visible=k == k_values[0], hoverinfo="skip")
    )

fig.update_layout(
    sliders=[dict(
        active=0, currentvalue=dict(prefix="Rank 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]}])
            for k in k_values
        ],
    )]
)
fig.update_xaxes(tickvals=tick_centres(len(shown), gap),
                 ticktext=[f"person {p}" for p in shown_people],
                 showgrid=False, zeroline=False, showline=False, ticks="")
fig.update_yaxes(tickvals=tick_centres(2, gap), ticktext=["original", "rank k"],
                 showgrid=False, zeroline=False, showline=False, ticks="")
style(fig, height=460, title="Rows of the truncated SVD A<sub>k</sub>")
fig.update_layout(margin=dict(l=110, r=30, t=60, b=100))
fig.show()

Figure 8: One photograph each of six people, reconstructed from rank k truncations of the data matrix. Drag the slider to change k.

At $k=1$ every face is the same average face, only brighter or darker, which is what $\mathbf{v}_1$ predicted. By $k=20\sim 50$ the people are recognisable, and by $k=100$ the reconstructions are hard to tell from the originals.

--- title: "The SVD of a collection of faces" subtitle: "Math 123: Mathematical Aspects of Data Analysis — Fall 2026" date: 2026-09-17 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 --- The [Olivetti faces](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.fetch_olivetti_faces.html) are $400$ greyscale photographs, ten each of $40$ people. We will compress one face on its own first, and then take the SVD of the full $400\times4096$ data matrix. ```{python} import base64 import io import numpy as np import plotly.graph_objects as go from matplotlib import colormaps from PIL import Image from sklearn.datasets import fetch_olivetti_faces faces = fetch_olivetti_faces() # downloads about 4 MB on first use images = faces.images # (400, 64, 64), brightness in [0, 1] A = faces.data.astype(float) # (400, 4096), the same images flattened person = faces.target # which of the 40 people is in each image m, n = A.shape print(f"{m} images of {person.max() + 1} people, each 64 x 64 pixels") print(f"Data matrix A has shape {m} x {n}") print(f"Brightness runs from {A.min():.2f} to {A.max():.2f}") ``` ```{python} #| code-fold: true #| code-summary: "Figure style and image helpers" # A muted ink colour and a translucent grid read correctly against both the # light and the dark version of this site, so the figures need no theme switch. INK = "#7c7c78" GRID = "rgba(128, 128, 128, 0.28)" BLUE, ORANGE, AQUA, RED = "#2a78d6", "#eb6834", "#1baf7a", "#e34948" SIDE = 64 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 montage(rows, gap=4): """Tile a grid of 64x64 images into one array, separated by blank (NaN) pixels.""" n_rows, n_cols = len(rows), max(len(row) for row in rows) out = np.full((n_rows * (SIDE + gap) - gap, n_cols * (SIDE + gap) - gap), np.nan) for i, row in enumerate(rows): for j, tile in enumerate(row): top, left = i * (SIDE + gap), j * (SIDE + gap) out[top:top + SIDE, left:left + SIDE] = tile return out def png_uri(array, cmap="gray", vmin=0.0, vmax=1.0): """Encode an array as a PNG data URI; NaN pixels become transparent. A 64x64 heatmap costs thousands of numbers in the page, while the same picture as a PNG is a few kilobytes, so every face below is drawn this way. """ scaled = np.clip((array - vmin) / (vmax - vmin), 0, 1) rgba = (colormaps[cmap](np.nan_to_num(scaled)) * 255).astype(np.uint8) rgba[np.isnan(array), 3] = 0 buffer = io.BytesIO() Image.fromarray(rgba, "RGBA").save(buffer, "PNG") return "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode() def tick_centres(count, gap=4): return [j * (SIDE + gap) + (SIDE - 1) / 2 for j in range(count)] def image_figure(rows, col_labels=None, row_labels=None, title=None, height=None, cmap="gray", vmin=0.0, vmax=1.0, gap=4): """Show a grid of 64x64 images as a single picture.""" n_rows, n_cols = len(rows), max(len(row) for row in rows) fig = go.Figure( go.Image(source=png_uri(montage(rows, gap), cmap, vmin, vmax), hoverinfo="skip") ) fig.update_xaxes( tickvals=tick_centres(n_cols, gap), ticktext=col_labels or [""] * n_cols, showgrid=False, zeroline=False, showline=False, ticks="", ) fig.update_yaxes( tickvals=tick_centres(n_rows, gap), ticktext=row_labels or [""] * n_rows, showgrid=False, zeroline=False, showline=False, ticks="", ) style(fig, height=height or 130 * n_rows + 90, title=title) fig.update_layout(margin=dict(l=110 if row_labels else 30, r=30, t=60 if title else 20, b=50)) return fig def rank_k(U, s, Vt, k): """The truncated SVD A_k = U_k Sigma_k V_k^T.""" return (U[:, :k] * s[:k]) @ Vt[:k, :] def relative_frobenius_tail(s): """Entry k is ||A - A_k||_F / ||A||_F, computed from the singular values alone.""" tail = np.concatenate([np.cumsum(s[::-1] ** 2)[::-1], [0.0]]) return np.sqrt(tail / tail[0]) ``` ## The data {#sec-data} Each image is $64\times64$ with brightness between $0$ (black) and $1$ (white), so flattening one gives a vector in $\mathbb{R}^{4096}$. Stacking the $400$ flattened images as rows gives the [data matrix](../../lectures/lecture-03.qmd#def-data-matrix-review) $$ A= \begin{pmatrix} - & \mathbf{x}_1^\top & - \\ & \vdots & \\ - & \mathbf{x}_{400}^\top & - \end{pmatrix} \in\mathbb{R}^{400\times4096}. $$ The photographs of each person vary in pose, expression, lighting, and glasses, but they share a common framing: every face is roughly centred and the same size. ```{python} #| label: fig-face-gallery #| fig-cap: "Four of the ten photographs of eight of the forty people. Each column is one person." #| fig-alt: "A four-row gallery of greyscale face photographs, eight people across, showing variation in pose and expression within each column." #| code-fold: true #| code-summary: "Plot a gallery of faces" people = [0, 5, 10, 14, 20, 26, 29, 35] fig = image_figure( [[images[person == p][pose] for p in people] for pose in range(4)], col_labels=[f"person {p}" for p in people], title="Every row of A is a picture in R4096", height=560, ) fig.show() ``` ## One face as a matrix {#sec-one-face} Before touching the collection, take one photograph on its own as a matrix $F\in\mathbb{R}^{64\times64}$ with compact SVD $F=\sum_{i=1}^r \sigma_i \mathbf{u}_i\mathbf{v}_i^\top$. ```{python} F = images[0] U_f, s_f, Vt_f = np.linalg.svd(F, full_matrices=False) print("Rank of F :", np.linalg.matrix_rank(F)) print("Largest singular values:", s_f[:6].round(2)) print("Smallest singular value:", f"{s_f[-1]:.1e}") ``` A photograph is essentially never exactly low rank, but its singular values fall by an order of magnitude almost immediately. ```{python} #| label: fig-single-spectrum #| fig-cap: "Singular values of one 64 × 64 photograph, on a logarithmic scale." #| fig-alt: "A line chart of sixty-four singular values on a log scale, dropping from about forty to about five within two indices and then decaying slowly toward one thousandth." #| code-fold: true #| code-summary: "Plot the singular values of one face" fig = go.Figure( go.Scatter( x=np.arange(1, len(s_f) + 1), y=s_f, mode="lines+markers", line=dict(color=BLUE, width=2), marker=dict(size=6), hovertemplate="i = %{x} sigma = %{y:.3g}<extra></extra>", ) ) fig.update_xaxes(title="Index i", dtick=8) fig.update_yaxes(title="Singular value", type="log", exponentformat="power") style(fig, height=360, title="Singular values of a single face") fig.show() ``` The [truncated SVD](../../lectures/lecture-04.qmd#def-truncated-svd) $F_k=\sum_{i=1}^k\sigma_i\mathbf{u}_i\mathbf{v}_i^\top$ is the best rank-$k$ approximation of $F$ in both norms, and by the [spectral](../../lectures/lecture-04.qmd#thm-eckart-young) and [Frobenius](../../lectures/lecture-04.qmd#thm-eckart-young-frobenius) Eckart–Young theorems its relative errors are read off from the singular values: $$ \frac{\|F-F_k\|_2}{\|F\|_2} = \frac{\sigma_{k+1}}{\sigma_1}, \qquad \frac{\|F-F_k\|_F}{\|F\|_F} = \sqrt{\frac{\sum_{i>k}\sigma_i^2}{\sum_{i\geq1}\sigma_i^2}} . $$ ```{python} #| label: fig-single-reconstruction #| fig-cap: "Rank-k approximations of one photograph, with the relative spectral and Frobenius errors beneath each rank." #| fig-alt: "A row of seven face images: the original followed by its rank one, two, five, ten, twenty, and thirty approximations, from a blocky striped pattern to a nearly exact face." #| code-fold: true #| code-summary: "Plot rank-k approximations of one face" ranks = [1, 2, 5, 10, 20, 30] single_spectral_err = s_f / s_f[0] # entry k is sigma_{k+1} / sigma_1 single_frobenius_err = relative_frobenius_tail(s_f) fig = image_figure( [[F] + [rank_k(U_f, s_f, Vt_f, k) for k in ranks]], col_labels=["original"] + [ f"k = {k} 2-norm: {single_spectral_err[k]:.1%}" f" F-norm: {single_frobenius_err[k]:.1%}" for k in ranks ], title="Rank-k approximations of one 64 x 64 face", height=290, ) fig.update_layout(margin=dict(b=80)) fig.show() ``` ## The SVD of the data matrix {#sec-data-matrix} The rows of $A\in\mathbb{R}^{400\times4096}$ are whole faces, and a column is one pixel position measured across all $400$ photographs. Its SVD $$ A=\sum_{i=1}^r\sigma_i\mathbf{u}_i\mathbf{v}_i^\top, \qquad \mathbf{u}_i\in\mathbb{R}^{400}, \quad \mathbf{v}_i\in\mathbb{R}^{4096}, $$ has right singular vectors that are themselves $64\times64$ images, and left singular vectors that assign a number to each of the $400$ photographs. ### Computing the economy SVD NumPy computes the SVD directly with `np.linalg.svd`. By default it returns the [full SVD](../../lectures/lecture-04.qmd#thm-full-svd), whose $U$ is a square $400\times400$ orthogonal matrix and whose $V$ is a square $4096\times4096$ orthogonal matrix. That $V$ alone has almost $17$ million entries, and most of them are wasted. Since $\Sigma\in\mathbb{R}^{400\times4096}$ has only $p=\min\{400,4096\}=400$ diagonal entries, the last $4096-400=3696$ rows of $V^\top$ are multiplied by columns of zeros and never contribute to $A$. Passing `full_matrices=False` asks for the *economy* SVD instead, which keeps only the first $p$ singular vectors on each side: $$ A=U_p\Sigma_pV_p^\top, \qquad U_p\in\mathbb{R}^{400\times400}, \quad \Sigma_p\in\mathbb{R}^{400\times400}, \quad V_p\in\mathbb{R}^{4096\times400}. $$ Both $U_p$ and $V_p$ still have orthonormal columns, but now $V_p$ is tall and $V_pV_p^\top\neq I$. If $\operatorname{rank}(A)=r<p$, the last $p-r$ diagonal entries of $\Sigma_p$ are zero, and dropping them gives the [compact SVD](../../lectures/lecture-04.qmd#thm-svd) $U_r\Sigma_rV_r^\top$ from lecture. ```{python} U, sigma, Vt = np.linalg.svd(A, full_matrices=False) rank = np.linalg.matrix_rank(A) # Singular vectors are only defined up to sign. Fix one choice: make the entry # of largest magnitude in each v_i positive, and flip u_i to match, so the # product sigma_i u_i v_i^T is unchanged. signs = np.sign(Vt[np.arange(len(sigma)), np.abs(Vt).argmax(axis=1)]) U, Vt = U * signs, Vt * signs[:, None] print(f"U: {U.shape}, sigma: {sigma.shape}, Vt: {Vt.shape}") print(f"rank(A) = {rank}, smallest singular value = {sigma[-1]:.3f}") print(f"U^T U = I up to {np.abs(U.T @ U - np.eye(len(sigma))).max():.1e}") print(f"V^T V = I up to {np.abs(Vt @ Vt.T - np.eye(len(sigma))).max():.1e}") print(f"||A - U diag(sigma) V^T||_F / ||A||_F = " f"{np.linalg.norm(A - (U * sigma) @ Vt) / np.linalg.norm(A):.1e}") ``` The collection has full rank $r=p=400$: no photograph is exactly a combination of the others, so every singular value is positive and the economy and compact SVDs coincide. ### How fast do the singular values decay? ```{python} print("sigma_1 through sigma_6:", sigma[:6].round(1)) print(f"sigma_1 / sigma_2 = {sigma[0] / sigma[1]:.1f}") mean_face = A.mean(axis=0) cosine = Vt[0] @ mean_face / np.linalg.norm(mean_face) print(f"cos(angle between v_1 and the average face) = {cosine:.5f}") ``` The first singular value is more than ten times the second. Every row of $A$ is a nonnegative image of a centred face, so all $400$ rows point in roughly the same direction in $\mathbb{R}^{4096}$. The right singular vector $\mathbf{v}_1$ *is* that shared direction, and it is essentially the average face rescaled to unit length. ```{python} #| label: fig-spectrum #| fig-cap: "All 400 singular values of the face data matrix, on a logarithmic scale. The first stands far above the rest." #| fig-alt: "A line chart of four hundred singular values on a log scale, starting at about seven hundred, dropping to about seventy at the second index, then decaying smoothly to below one." #| code-fold: true #| code-summary: "Plot the singular values of the data matrix" fig = go.Figure( go.Scatter( x=np.arange(1, rank + 1), y=sigma, mode="lines", line=dict(color=BLUE, width=2), hovertemplate="i = %{x} sigma = %{y:.3g}<extra></extra>", ) ) fig.add_annotation(x=1, y=np.log10(sigma[0]), text=" sigma1: the average face", showarrow=False, xanchor="left", font=dict(color=BLUE, size=13)) fig.update_xaxes(title="Index i", dtick=50) fig.update_yaxes(title="Singular value", type="log", exponentformat="power") style(fig, height=380, title="Singular values of the 400 x 4096 face matrix") fig.show() ``` ### Eckart–Young, numerically The [Eckart–Young theorem](../../lectures/lecture-04.qmd#thm-eckart-young) makes two promises about $A_k=U_k\Sigma_kV_k^\top$: the error is exactly the discarded singular values, $$ \|A-A_k\|_2=\sigma_{k+1}, \qquad \|A-A_k\|_F=\Big(\sum_{i>k}\sigma_i^2\Big)^{1/2}, $$ and no other matrix of rank at most $k$ does better. The first promise we can check directly. For the second, we cannot try every rank-$k$ matrix, but we can try a natural competitor: pick $k$ of the photographs at random and replace every face by its closest point in their span. That matrix $B$ also has rank at most $k$. ```{python} rng = np.random.default_rng(0) def random_faces_approximation(k): """Project every row of A onto the span of k randomly chosen rows.""" Q, _ = np.linalg.qr(A[rng.choice(m, size=k, replace=False)].T) # (4096, k) return (A @ Q) @ Q.T frob_tail = relative_frobenius_tail(sigma) * np.linalg.norm(A) print(f"{'k':>4} {'||A-A_k||_2':>12} {'sigma_k+1':>10} {'||A-A_k||_F':>12} " f"{'tail':>8} {'||A-B||_2':>10} {'||A-B||_F':>10}") for k in (1, 5, 20, 50, 100): Ak = rank_k(U, sigma, Vt, k) B = random_faces_approximation(k) print(f"{k:>4} {np.linalg.norm(A - Ak, 2):>12.3f} {sigma[k]:>10.3f} " f"{np.linalg.norm(A - Ak):>12.3f} {frob_tail[k]:>8.3f} " f"{np.linalg.norm(A - B, 2):>10.3f} {np.linalg.norm(A - B):>10.3f}") ``` The first two pairs of columns match to every printed digit, and the random-faces competitor is worse in both norms at every $k$. ```{python} #| label: fig-eckart-young #| fig-cap: "Relative error of the truncated SVD and of projecting onto k random photographs, in the spectral norm (left) and the Frobenius norm (right). The SVD curves come from the singular values alone." #| fig-alt: "Two line charts of relative error against rank k. In each, the truncated SVD curve lies below the random-photographs curve at every rank." #| code-fold: true #| code-summary: "Plot error against rank for both norms" from plotly.subplots import make_subplots k_grid = [1, 2, 3, 5, 8, 12, 20, 30, 50, 75, 100, 150, 200, 300] competitor = [A - random_faces_approximation(k) for k in k_grid] curves = { "spectral": (sigma[k_grid] / sigma[0], [np.linalg.norm(E, 2) / sigma[0] for E in competitor]), "Frobenius": (frob_tail[k_grid] / frob_tail[0], [np.linalg.norm(E) / frob_tail[0] for E in competitor]), } fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.1, subplot_titles=["||A - B||2 / ||A||2", "||A - B||F / ||A||F"]) for col, (norm, (svd_curve, other_curve)) in enumerate(curves.items(), start=1): for label, curve, colour in (("truncated SVD A_k", svd_curve, BLUE), ("k random photographs", other_curve, ORANGE)): fig.add_trace( go.Scatter( x=k_grid, y=curve, mode="lines+markers", name=label, line=dict(color=colour, width=2), marker=dict(size=6), showlegend=col == 1, hovertemplate=f"{label} k = %{{x}} {norm} error %{{y:.1%}}<extra></extra>", ), row=1, col=col, ) fig.update_xaxes(title="Rank k", type="log") fig.update_yaxes(tickformat=".0%", rangemode="tozero") style(fig, height=420, title="The truncated SVD beats a natural competitor in both norms") fig.update_layout(legend=dict(orientation="h", y=-0.22, x=0)) fig.update_layout(margin=dict(b=90)) fig.show() ``` Relative to $\|A\|_2=\sigma_1$ the spectral error collapses to about $10\%$ after a single term, because $\sigma_1$ is so dominant. The Frobenius error falls much more slowly: no single missing piece is large, but there are many of them. ### What do the right singular vectors look like? Each $\mathbf{v}_i\in\mathbb{R}^{4096}$ reshapes to a $64\times64$ image. ```{python} #| label: fig-singular-vectors #| fig-cap: "The first twelve right singular vectors, each reshaped to 64 × 64. Red is positive and blue is negative, on a scale symmetric about zero." #| fig-alt: "Twelve red-and-blue face-shaped patterns. The first is uniformly red and resembles an average face; later ones contrast left against right, top against bottom, and highlight the eyes, nose, and hairline." #| code-fold: true #| code-summary: "Plot the leading right singular vectors" tiles = [Vt[i].reshape(SIDE, SIDE) for i in range(12)] limit = float(np.abs(tiles).max()) fig = image_figure( [tiles[:6], tiles[6:]], col_labels=[""] * 6, title="Right singular vectors v1 through v12", height=400, cmap="RdBu_r", vmin=-limit, vmax=limit, ) for row in range(2): for col in range(6): fig.add_annotation( x=tick_centres(6)[col], y=row * (SIDE + 4) - 5, text=f"v{row * 6 + col + 1}", showarrow=False, font=dict(color=INK, size=13), ) fig.show() ``` ## Low-rank reconstructions {#sec-reconstruct} Storing $A_k$ costs $k(400+4096+1)$ numbers, against $400\times4096$ for $A$ itself. The plot below sets that cost against both relative errors from Eckart–Young, $\|A-A_k\|_2/\|A\|_2=\sigma_{k+1}/\sigma_1$ and $\|A-A_k\|_F/\|A\|_F$, at every rank shown in the figures that follow. ```{python} #| label: fig-error-storage #| fig-cap: "Relative spectral and Frobenius error of A_k (left axis) and the storage cost of A_k as a fraction of the raw pixels (right axis), at every rank shown below. The dotted line marks where A_k costs as much as A." #| fig-alt: "A line chart against rank k on a log scale. Two error curves fall from about ten and twenty percent at rank one to zero at rank four hundred, while a storage curve rises from almost nothing to about one hundred ten percent, crossing a dotted one hundred percent line just before rank four hundred." #| code-fold: true #| code-summary: "Plot error and storage against rank" k_values = [1, 2, 3, 5, 10, 20, 50, 100, 200, 300, 400] # every rank shown below ks = np.array(k_values) spectral_error = np.append(sigma, 0.0)[ks] / sigma[0] # sigma_{k+1} / sigma_1 frobenius_error = frob_tail[ks] / frob_tail[0] storage = ks * (m + n + 1) / (m * n) fig = make_subplots(specs=[[{"secondary_y": True}]]) for label, curve, colour, secondary in ( ("relative spectral error", spectral_error, AQUA, False), ("relative Frobenius error", frobenius_error, BLUE, False), ("storage (fraction of A)", storage, ORANGE, True), ): fig.add_trace( go.Scatter( x=ks, y=curve, mode="lines+markers", name=label, line=dict(color=colour, width=2, dash="dash" if secondary else "solid"), marker=dict(size=7), hovertemplate=f"{label} k = %{{x}} %{{y:.1%}}<extra></extra>", ), secondary_y=secondary, ) fig.add_hline(y=1, line=dict(color=ORANGE, width=1, dash="dot"), secondary_y=True) fig.update_xaxes(title="Rank k", type="log", tickvals=k_values) fig.update_yaxes(title="Relative error", tickformat=".0%", rangemode="tozero", secondary_y=False) fig.update_yaxes(title="Storage of Ak / storage of A", tickformat=".0%", rangemode="tozero", showgrid=False, title_font=dict(color=ORANGE), tickfont=dict(color=ORANGE), secondary_y=True) style(fig, height=440, title="What each rank costs, and what it buys", showlegend=True) fig.update_layout(legend=dict(orientation="h", y=-0.2, x=0), margin=dict(l=70, r=80, t=60, b=100)) fig.show() ``` Move the slider to watch six faces come back into focus. ```{python} #| label: fig-reconstruction-slider #| fig-cap: "One photograph each of six people, reconstructed from rank k truncations of the data matrix. Drag the slider to change k." #| fig-alt: "An interactive grid of six original faces above their reconstructions, with a slider controlling the rank from a single blurred average face at k equals one to exact faces at k equals four hundred." #| code-fold: true #| code-summary: "Build the reconstruction slider" shown_people = [6, 10, 14, 26, 29, 35] shown = np.array([np.flatnonzero(person == p)[0] for p in shown_people]) gap = 4 fig = go.Figure() for k in k_values: reconstruction = rank_k(U[shown], sigma, Vt, k).reshape(-1, SIDE, SIDE) fig.add_trace( go.Image(source=png_uri(montage([images[shown], reconstruction], gap)), visible=k == k_values[0], hoverinfo="skip") ) fig.update_layout( sliders=[dict( active=0, currentvalue=dict(prefix="Rank 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]}]) for k in k_values ], )] ) fig.update_xaxes(tickvals=tick_centres(len(shown), gap), ticktext=[f"person {p}" for p in shown_people], showgrid=False, zeroline=False, showline=False, ticks="") fig.update_yaxes(tickvals=tick_centres(2, gap), ticktext=["original", "rank k"], showgrid=False, zeroline=False, showline=False, ticks="") style(fig, height=460, title="Rows of the truncated SVD Ak") fig.update_layout(margin=dict(l=110, r=30, t=60, b=100)) fig.show() ``` At $k=1$ every face is the same average face, only brighter or darker, which is what $\mathbf{v}_1$ predicted. By $k=20\sim 50$ the people are recognisable, and by $k=100$ the reconstructions are hard to tell from the originals.