Demo: Low-Rank Approximation of Images and Video

We will compress a grayscale version of the Ivy and Stripes photograph, then recover the background of a highway video. The road and trees appear in every frame, while cars move through the scene. Putting the frames into one matrix lets the SVD find this shared background.

The notebook uses a photo, sampled traffic frames, and the course’s plotting style. The setup downloads any missing files once. The video plays directly in the notebook and on the course page.

from pathlib import Path
from urllib.request import Request, urlopen

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from PIL import Image, ImageOps
from IPython.display import HTML

course_source = "https://csci145.rtealwitter.com/topics/"
image_path = Path("images/ivy_and_stripes.jpg")
video_path = Path("media/highway_frames.npz")
style_path = Path("../figure_style.mplstyle")
if not style_path.exists():
    style_path = Path("figure_style.mplstyle")
for local, remote in [(image_path, "15_low_rank_approximation/images/ivy_and_stripes.jpg"),
                      (style_path, "figure_style.mplstyle"),
                      (video_path, "15_low_rank_approximation/media/highway_frames.npz")]:
    if not local.exists():
        local.parent.mkdir(parents=True, exist_ok=True)
        request = Request(course_source + remote, headers={"User-Agent": "CSCI145-demo"})
        with urlopen(request, timeout=30) as response:
            local.write_bytes(response.read())
plt.style.use(style_path)

photo = ImageOps.exif_transpose(Image.open(image_path)).convert("L")
photo.thumbnail((384, 384))
A = np.asarray(photo, dtype=float) / 255
m, n = A.shape
r = np.linalg.matrix_rank(A)
U, singular_values, Vt = np.linalg.svd(A, full_matrices=False)
U, singular_values, Vt = U[:, :r], singular_values[:r], Vt[:r, :]
total_energy = np.sum(singular_values ** 2)
print(f"Ivy and Stripes: {m} × {n} pixels, rank {r}.")

Ivy and Stripes: 384 × 288 pixels, rank 288.

Question 1: What Survives at Low Rank?

The grayscale image is a matrix \(\mathbf A\in\mathbb R^{m\times n}\), with \(r=\operatorname{rank}(\mathbf A)\). We keep the first \(k\) singular components to form \(\mathbf A_k\). Predict how many directions it will take to distinguish the two cats and recover the texture of their fur.

ranks = [1, 4, 16, 64, 256]
fig, axes = plt.subplots(2, 3, figsize=(8, 6.5))
for ax, k in zip(axes.flat, ranks):
    approximation = (U[:, :k] * singular_values[:k]) @ Vt[:k, :]
    ax.imshow(approximation, cmap="gray", vmin=0, vmax=1)
    ax.set_title(f"Rank {k}")
    ax.set_axis_off()
axes.flat[-1].imshow(A, cmap="gray", vmin=0, vmax=1)
axes.flat[-1].set_title("Original")
axes.flat[-1].set_axis_off()
fig.tight_layout()
plt.show()

Broad shapes appear before fine detail. The singular values let us measure how much squared error remains at each rank.

Question 2: Does the Spectral Tail Predict Pixel Error?

Choose \(k\) below and compare the error measured from the pixels with the error predicted by the discarded singular values:

\[ \|\mathbf A-\mathbf A_k\|_F =\sqrt{\sum_{i=k+1}^{r}\sigma_i^2}. \]

The cumulative curve shows the fraction \(\sum_{i=1}^k\sigma_i^2/\sum_{i=1}^r\sigma_i^2\) retained. The factorization stores \(k(m+n+1)\) numbers, compared with \(mn\) pixel values.

k = 16  # Change the retained rank, from 1 through r.
approximation = (U[:, :k] * singular_values[:k]) @ Vt[:k, :]
measured_error = np.linalg.norm(A - approximation, ord="fro")
tail_error = np.sqrt(np.sum(singular_values[k:] ** 2))
assert np.isclose(measured_error, tail_error)

indices = np.arange(1, r + 1)
cumulative = np.cumsum(singular_values ** 2) / total_energy
fig, axes = plt.subplots(1, 2, figsize=(8, 3))
axes[0].plot(indices, singular_values)
axes[0].set_yscale("log")
axes[0].set_xlabel("Singular-Value Index $i$")
axes[0].set_ylabel(r"Singular Value $\sigma_i$")
axes[1].plot(indices, cumulative)
axes[1].set_xlabel("Rank $k$")
axes[1].set_ylabel("Fraction of Energy Kept")
axes[1].set_ylim(0, 1.02)
for ax in axes:
    ax.axvline(k, color="#333333", linestyle="--", linewidth=1)
fig.tight_layout()
plt.show()
print(f"Measured error: {measured_error:.6f}; spectral-tail prediction: {tail_error:.6f}")
print(f"Rank {k}: {cumulative[k-1]:.1%} energy retained, {k*(m+n+1)/(m*n):.1%} of the original number of values.")

Measured error: 23.190532; spectral-tail prediction: 23.190532
Rank 16: 97.8% energy retained, 9.7% of the original number of values.

The two errors agree. Increasing \(k\) reduces reconstruction error and increases storage, so a high-rank factorization can cost more than storing the pixels directly.

Question 3: Can One Singular Vector Reveal the Background?

A fixed camera looks down on cars passing along a highway. We want to recover the road behind the cars using the structure shared across frames. The clip below contains 96 grayscale frames sampled over 24 seconds of A40 highway footage.

Let \(T\) be the number of frames, each with height \(h\) and width \(w\). Flatten each frame into one row of \(\mathbf B\in\mathbb R^{T\times(hw)}\). In its SVD \(\mathbf B=\mathbf U\boldsymbol\Sigma\mathbf V^\top\), the leading right singular vector \(\mathbf v_1\in\mathbb R^{hw}\) has one entry per pixel, and the leading left singular vector \(\mathbf u_1\in\mathbb R^T\) has one entry per frame. The rank-one approximation is

\[ \mathbf B_1=\sigma_1\mathbf u_1\mathbf v_1^\top. \]

Every reconstructed frame is a multiple of the same image: \(\mathbf v_1\) reshaped to \(h\times w\). Its brightness in frame \(t\) is set by \(a_t=\sigma_1(\mathbf u_1)_t\), for \(t=1,\ldots,T\). Because the road and trees recur throughout the clip, this shared image should look mostly like the background.

We keep the frames uncentered to preserve the background. Subtracting each pixel’s average across frames would remove the persistent scene before we compute the SVD.

with np.load(video_path) as video:
    frames = video["frames"].astype(float) / 255
    fps = float(video["fps"])
T, h, w = frames.shape
B = frames.reshape(T, h * w)  # One row per frame, one column per pixel.
U_video, s_video, Vt_video = np.linalg.svd(B, full_matrices=False)

u1, v1 = U_video[:, 0].copy(), Vt_video[0, :].copy()
# Either simultaneous sign choice gives the same rank-one matrix.
if v1.mean() < 0:
    u1, v1 = -u1, -v1
coefficients = s_video[0] * u1
B1 = np.outer(coefficients, v1)
background_frames = B1.reshape(T, h, w)

# Rescale the unit vector to the average brightness of the reconstruction.
background = (coefficients.mean() * v1).reshape(h, w)
assert np.isclose(np.linalg.norm(v1), 1)
assert np.allclose(B1, (U_video[:, :1] * s_video[:1]) @ Vt_video[:1, :])
assert np.isclose(np.linalg.norm(B - B1, ord="fro") ** 2, np.sum(s_video[1:] ** 2))
print(f"Video matrix: {T} frames × {h*w:,} pixels ({h} × {w} per frame).")
print(f"Leading right singular vector: {v1.size:,} pixel values.")

fig, axes = plt.subplots(1, 2, figsize=(8, 3))
axes[0].imshow(frames[0], cmap="gray", vmin=0, vmax=1)
axes[0].set_title("One Video Frame")
axes[1].imshow(background, cmap="gray", vmin=0, vmax=1)
axes[1].set_title(r"Background From $\mathbf{v}_1$")
for ax in axes:
    ax.set_axis_off()
fig.tight_layout()
plt.show()

Video matrix: 96 frames × 36,864 pixels (144 × 256 per frame).
Leading right singular vector: 36,864 pixel values.

The right panel displays \(\mathbf v_1\) at the average reconstructed brightness, since a unit vector has a much smaller scale than the original pixel intensities. The road and trees remain, while the moving vehicles largely disappear.

Question 4: What Changes When We Remove the Background?

Subtract the rank-one reconstruction from the video matrix. The residual \(\mathbf R=\mathbf B-\mathbf B_1\) records where each frame differs from the shared scene. Predict where it will be brightest, then play the video to follow those differences.

The right panel shows the absolute residual with brightness multiplied by four and clipped to the display range, so smaller differences are visible. All three panels show the same frame at the same time.

def show_video(sequences, titles, fps=12):
    """Play synchronized frame sequences; the SVD is computed in the cells above."""
    import re

    fig, axes = plt.subplots(1, len(sequences), figsize=(9, 2.3), dpi=80)
    pictures = []
    for ax, sequence, title in zip(axes, sequences, titles):
        pictures.append(ax.imshow(sequence[0], cmap="gray", vmin=0, vmax=1))
        ax.set_title(title)
        ax.set_axis_off()
    fig.tight_layout()

    def update(t):
        for picture, sequence in zip(pictures, sequences):
            picture.set_data(sequence[t])
        return pictures

    movie = FuncAnimation(fig, update, frames=len(sequences[0]), interval=1000/fps)
    # JPEG keeps the player compact; SVD calculations use the original pixels.
    with plt.rc_context({"animation.frame_format": "jpeg"}):
        html = movie.to_jshtml(default_mode="loop")
    # Collapse JavaScript line continuations so Quarto can render longer clips.
    html = html.replace("\\\n", "")
    plt.close(fig)
    # Use built-in text symbols so playback needs no external icon fonts or codec.
    html = re.sub(r'<link[^>]*font-awesome[^>]*>', '', html)
    icons = {"minus": "−", "fast-backward": "|◀", "step-backward": "◀|",
             "play fa-flip-horizontal": "◀", "pause": "Ⅱ", "play": "▶",
             "step-forward": "|▶", "fast-forward": "▶|", "plus": "+"}
    for icon, symbol in icons.items():
        html = html.replace(f'<i class="fa fa-{icon}"></i>', symbol)
    html = html.replace('class="animation"', 'class="animation low-rank-video"')
    html += """<style>
    .low-rank-video { display: block; max-width: 100%; margin: 1em auto; }
    .low-rank-video img { display: block; max-width: 100%; height: auto; margin: auto; }
    .low-rank-video input.anim-slider { width: 100%; max-width: 374px; }
    .low-rank-video .anim-buttons { display: flex; flex-wrap: wrap; justify-content: center; gap: 4px; }
    .low-rank-video .anim-buttons button { background: white; color: #333; border: 1px solid #ccc; border-radius: 3px; }
    </style>"""
    return HTML("\n".join(line.rstrip() for line in html.splitlines()))

residual_frames = (B - B1).reshape(T, h, w)
show_video([frames, background_frames, np.clip(4 * np.abs(residual_frames), 0, 1)],
           ["Original", "Rank-One Background", "Difference (×4)"], fps=fps)

Once Loop Reflect

The moving cars stand out in the residual. Shadows and small camera movements also contribute, and vehicles that occupy the same pixels for many frames can leave traces in the background. Rank one minimizes squared reconstruction error, so it has no rule for deciding which objects count as background. The separation works here because the fixed scene is repeated much more consistently than the cars.

PCA in the embeddings reading uses centered observations to describe variation around their mean. Here the persistent mean image is part of the signal we want, which is why we applied the SVD to the uncentered frames.

Footage by Kathinka Engels / Sounds of Changes, Motorway A40: on a bridge above the traffic, licensed under CC BY 3.0. We sampled, resized, and converted the clip to grayscale, without audio. Source and preparation details.