Eigenvalues and SVD without tears: compress an image to see them work

LEARN · MATHEMATICS FOR MACHINE LEARNING

An RGB image looks visual to us, but to NumPy it is three matrices: one each for red, green, and blue. Real images contain repeated structure—smooth gradients, correlated colors, recurring edges, and regions that vary together.

Singular value decomposition, or SVD, exposes that structure. Keep only the strongest components, and you can rebuild a recognizable image from far fewer numbers.

In this lesson, you will:

Generate a deterministic RGB image using NumPy only.
Compute an SVD for each color channel.
Reconstruct the image at several ranks.
Compare visual quality, numerical error, and logical storage.
Connect singular values to eigenvalues.
See why LoRA fine-tuning uses the same low-rank idea for model weight updates.

The intuition: an image is a stack of simple patterns

Take one grayscale channel represented by a matrix A with m rows and n columns. SVD factors it into:

A = U · Σ · Vᵀ

The pieces have distinct jobs:

U describes patterns along the image’s rows.
Vᵀ describes patterns along its columns.
Σ contains nonnegative singular values ordered from largest to smallest.

The decomposition can also be written as a sum:

A = σ₁u₁v₁ᵀ + σ₂u₂v₂ᵀ + … + σᵣuᵣvᵣᵀ

Each outer product uᵢvᵢᵀ is a rank-one image pattern. Its singular value σᵢ says how strongly that pattern contributes.

A rank-k approximation keeps only the first k terms:

Aₖ = σ₁u₁v₁ᵀ + … + σₖuₖvₖᵀ

This is not an arbitrary shortcut. Truncated SVD gives the best rank-k approximation under both the Frobenius norm and spectral norm. For a fixed rank, no other rank-k matrix can reproduce A with less error under those measures.

Build a NumPy-only image compressor

We will generate an image instead of loading one, so the project needs no Pillow, OpenCV, or Matplotlib. The script writes binary PPM files, a simple RGB format supported by many image editors and command-line viewers.

NumPy 2.5.1 is a current patch release and supports Python 3.12 through 3.14.

Install it:

python -m pip install numpy==2.5.1

Save the following as svd_image.py:

from pathlib import Path

import numpy as np


OUTPUT_DIR = Path("svd_output")
RANKS = (2, 5, 10, 20, 40)
Factors = list[tuple[np.ndarray, np.ndarray, np.ndarray]]


def save_ppm(path: Path, image: np.ndarray) -> None:
    """Write a height × width × 3 array as a binary PPM image."""
    image_u8 = np.clip(np.rint(image), 0, 255).astype(np.uint8)

    if image_u8.ndim != 3 or image_u8.shape[2] != 3:
        raise ValueError(
            "Expected an RGB array with shape (height, width, 3)."
        )

    height, width, _ = image_u8.shape
    with path.open("wb") as file:
        file.write(f"P6\n{width} {height}\n255\n".encode("ascii"))
        file.write(image_u8.tobytes())


def make_demo_image(
    height: int = 192,
    width: int = 256,
) -> np.ndarray:
    """Create a deterministic image with gradients and sharp shapes."""
    y, x = np.mgrid[0:height, 0:width]
    xn = x / (width - 1)
    yn = y / (height - 1)

    red = 255 * (
        0.15
        + 0.60 * xn
        + 0.18 * np.sin(6 * np.pi * yn) ** 2
    )
    green = 255 * (
        0.10
        + 0.58 * yn
        + 0.25 * np.cos(5 * np.pi * (xn + yn)) ** 2
    )
    blue = 255 * (
        0.48
        + 0.22
        * np.sin(8 * np.pi * xn)
        * np.cos(6 * np.pi * yn)
    )

    image = np.stack((red, green, blue), axis=-1)

    circle = (
        (x - 0.30 * width) ** 2
        + (y - 0.60 * height) ** 2
        < (0.16 * min(height, width)) ** 2
    )
    image[circle] = (245, 215, 55)

    rectangle = (
        (x > 0.62 * width)
        & (x < 0.90 * width)
        & (y > 0.18 * height)
        & (y < 0.42 * height)
    )
    image[rectangle] = (45, 80, 225)

    diagonal = (
        np.abs(y - (0.76 * height - 0.42 * x)) < 2.5
    )
    image[diagonal] = (250, 250, 250)

    return np.clip(image, 0, 255).astype(np.float64)


def factorize_rgb(image: np.ndarray) -> Factors:
    """Compute a reduced SVD independently for each RGB channel."""
    channels = np.moveaxis(image, -1, 0)
    return [
        np.linalg.svd(channel, full_matrices=False)
        for channel in channels
    ]


def reconstruct_rgb(
    factors: Factors,
    rank: int,
) -> np.ndarray:
    """Reconstruct an RGB image from at most `rank` components."""
    if rank < 1:
        raise ValueError("Rank must be at least 1.")

    channels = []

    for u, singular_values, vh in factors:
        r = min(rank, singular_values.size)
        channel = (
            u[:, :r] * singular_values[:r]
        ) @ vh[:r, :]
        channels.append(channel)

    return np.stack(channels, axis=-1)


def quality_metrics(
    original: np.ndarray,
    reconstructed: np.ndarray,
) -> tuple[float, float]:
    """Return RMSE and PSNR for 8-bit image values."""
    rmse = float(
        np.sqrt(np.mean((original - reconstructed) ** 2))
    )
    psnr = (
        float("inf")
        if rmse == 0
        else float(20 * np.log10(255.0 / rmse))
    )
    return rmse, psnr


def retained_energy(
    factors: Factors,
    rank: int,
) -> float:
    """Return the fraction of squared singular values retained."""
    kept = 0.0
    total = 0.0

    for _, singular_values, _ in factors:
        r = min(rank, singular_values.size)
        squared = singular_values**2
        kept += float(np.sum(squared[:r]))
        total += float(np.sum(squared))

    return kept / total


def main() -> None:
    OUTPUT_DIR.mkdir(exist_ok=True)

    image = make_demo_image()
    factors = factorize_rgb(image)
    height, width, channels = image.shape

    save_ppm(OUTPUT_DIR / "original.ppm", image)

    print("rank | RMSE   | PSNR     | energy  | factor ratio")
    print("-----+--------+----------+---------+-------------")

    reconstructions = []

    for rank in RANKS:
        reconstructed = reconstruct_rgb(factors, rank)
        reconstructions.append(reconstructed)

        save_ppm(
            OUTPUT_DIR / f"rank_{rank:02d}.ppm",
            reconstructed,
        )

        rmse, psnr = quality_metrics(image, reconstructed)
        energy = retained_energy(factors, rank)

        original_values = image.size
        factor_values = (
            channels * rank * (height + width + 1)
        )
        ratio = original_values / factor_values

        print(
            f"{rank:>4} | "
            f"{rmse:>6.2f} | "
            f"{psnr:>7.2f} dB | "
            f"{energy:>6.2%} | "
            f"{ratio:>10.2f}x"
        )

    gap = np.full((height, 8, channels), 255.0)
    montage = [image]

    for reconstructed in reconstructions:
        montage.extend((gap, reconstructed))

    save_ppm(
        OUTPUT_DIR / "comparison.ppm",
        np.concatenate(montage, axis=1),
    )

    print(f"\nImages written to {OUTPUT_DIR.resolve()}")
    print("Comparison order: original, ranks 2, 5, 10, 20, 40")


if __name__ == "__main__":
    main()

Run it:

python svd_image.py

The output directory will contain:

original.ppm
rank_02.ppm
rank_05.ppm
rank_10.ppm
rank_20.ppm
rank_40.ppm
comparison.ppm

Open comparison.ppm. From left to right, you will see the original followed by increasingly accurate low-rank reconstructions.

Read the results

A representative run produces:

Rank	RMSE	PSNR	Energy retained	Logical factor ratio
2	29.13	18.84 dB	95.97%	54.73×
5	19.10	22.51 dB	98.27%	21.89×
10	14.44	24.94 dB	99.01%	10.95×
20	8.90	29.15 dB	99.62%	5.47×
40	5.02	34.11 dB	99.88%	2.74×

RMSE and PSNR

RMSE measures the typical pixel-value error:

RMSE = √mean((original − reconstruction)²)

Lower is better. PSNR expresses the error on a logarithmic scale:

PSNR = 20 · log₁₀(255 / RMSE)

Higher is better here, although PSNR is not a perfect model of human perception.

Retained energy

The squared singular values add up to the matrix’s squared Frobenius norm:

Energy retained = sum of first k squared singular values / sum of all squared singular values

Notice the surprise: rank 2 keeps nearly 96% of the numerical energy, yet obvious visual errors remain.

Smooth backgrounds and broad color fields dominate the calculation. A thin white edge may contribute little energy while remaining highly noticeable to your eyes.

Is this real image compression?

The script reports a logical ratio based on stored scalar values.

For an RGB image:

Original values = 3mn
Rank-k factor values = 3k(m + n + 1)

For a 192 × 256 image at rank 10:

Original: 147,456 values
SVD factors: 13,470 values
Logical ratio: about 10.95×

That does not guarantee a 10.95× smaller file. The original pixels are 8-bit integers, while the factors are floating-point values. A practical format also needs quantization, metadata, and entropy coding.

The PPM files contain full reconstructed pixels, so this project demonstrates approximation rather than a competitive codec. Production formats such as JPEG, WebP, and AVIF are engineered specifically for storage and transmission.

Where eigenvalues enter the picture

Eigenvalue decomposition applies directly to square matrices. An image channel A is usually rectangular, but AᵀA and AAᵀ are square and symmetric.

Given:

A = U · Σ · Vᵀ

then:

AᵀA = V · Σ² · Vᵀ

AAᵀ = U · Σ² · Uᵀ

Therefore:

The right singular vectors are eigenvectors of AᵀA.
The left singular vectors are eigenvectors of AAᵀ.
Each nonzero singular value is the square root of a corresponding eigenvalue.

NumPy’s stable SVD documentation describes these shapes, descending singular values, and reduced reconstruction as (u * s) @ vh when full_matrices=False.

You can verify the relationship for the red channel:

import numpy as np

from svd_image import make_demo_image


channel = make_demo_image()[:, :, 0]

singular_values = np.linalg.svd(
    channel,
    compute_uv=False,
)

eigenvalues = np.linalg.eigvalsh(
    channel.T @ channel
)
from_eigenvalues = np.sqrt(
    np.clip(eigenvalues[::-1], 0.0, None)
)[: singular_values.size]

significant = (
    singular_values
    > singular_values[0] * 1e-8
)

print(
    np.allclose(
        singular_values[significant],
        from_eigenvalues[significant],
        rtol=1e-5,
        atol=1e-6,
    )
)

Expected output:

True

Do not compute an SVD in production by explicitly forming AᵀA. That squares the condition number, making small singular values less accurate. Call np.linalg.svd directly.

Why `full_matrices=False` matters

Let p = min(m, n). Reduced SVD returns:

U with shape m × p
Singular values with length p
Vᵀ with shape p × n

A full decomposition may allocate larger square matrices that add no useful information for truncated reconstruction. For low-rank image work, reduced SVD is the natural choice.

Cherry on the cake: this is the idea behind LoRA

Low-Rank Adaptation, or LoRA, applies the same structural idea to neural-network fine-tuning.

Suppose a model contains a frozen weight matrix W. Instead of learning an unrestricted update ΔW with the same dimensions, LoRA represents the update with two narrow matrices:

ΔW = B · A

If their shared inner dimension is r, then ΔW has rank at most r. The model trains the smaller factors rather than every element of a full-sized update.

Current Hugging Face PEFT documentation exposes this rank through the LoRA configuration and describes LoRA as a parameter-efficient fine-tuning method.

The connection is exact at the structural level:

Image SVD represents an image using low-rank components.
LoRA represents a weight update using a low-rank factorization.
In both cases, rank controls compactness versus expressive power.

The procedures are not identical. Truncated SVD analyzes an existing matrix and keeps its strongest singular directions. Standard LoRA learns its factors during training; it does not normally run SVD on the weight update first.

That blurry rank-2 reconstruction is a useful mental model for model adaptation: a few directions can capture a large share of the change, but fine details require more rank.

Sources

Run the script, open comparison.ppm, and choose the lowest rank that still looks acceptable to you. Then carry that visual intuition into the next lesson: dimensionality reduction with PCA.