SVD Example — Low-Rank Approximation

This example builds a 4×6 matrix with low-rank structure and decomposes it using eigensystems.svd(). Rank-1, rank-2, and rank-3 truncated approximations are computed, and reconstruction errors show how quickly the singular values capture the matrix’s energy.

Example Code

"""IMSL SVD example: singular value decomposition and low-rank approximation.

Demonstrates:
- Create a 4x6 image-like matrix (low-rank structure).
- Compute the full SVD using :func:`eigensystems.svd`.
- Build rank-1, rank-2, and rank-3 approximations.
- Plot original matrix and approximations as heatmaps.
- Bar chart of singular values.
- Save to test_output/example_imsl_svd.svg

Outputs:
- Printed table of singular values and reconstruction errors.
- SVG figure saved to test_output/.
- Console output captured to docs/_generated_text/svd_output.txt
"""
from __future__ import annotations

from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

from eigensystems import svd


def run_demo_svd() -> dict:
    """Run IMSL SVD example: low-rank approximation of a 4x6 matrix.

    Constructs a matrix with clear low-rank structure (sum of outer products),
    computes the full SVD via :func:`eigensystems.svd`, and shows how
    successive rank-k truncations reconstruct the original data.

    Args:
        None

    Returns:
        dict: Result dictionary with keys ``singular_values`` (np.ndarray),
            ``errors`` (list[float]), and ``plot_path`` (str).
    """
    rng = np.random.default_rng(42)

    # Build a rank-3 matrix: sum of three outer products
    u_true = np.array([[1, 0, 0], [0, 1, 0], [1, 1, 0], [0, 0, 1]], dtype=float)
    v_true = np.array([
        [2, 1, 0, 1, 0, 0],
        [0, 1, 3, 0, 1, 0],
        [1, 0, 0, 2, 0, 1],
    ], dtype=float)
    A = u_true @ v_true + 0.1 * rng.standard_normal((4, 6))

    U, s, Vh = svd(A)

    print("\nIMSL SVD Example: Low-Rank Approximation")
    print(f"Matrix shape: {A.shape}   (rank ≈ 3)")
    print("-" * 48)
    print(f"{'Index':<8} {'Singular value':>18} {'Energy %':>10}")
    print("-" * 48)
    total_energy = np.sum(s ** 2)
    for i, sv in enumerate(s):
        energy_pct = 100.0 * sv ** 2 / total_energy
        print(f"  σ{i+1:<5} {sv:>18.6f} {energy_pct:>10.2f}%")
    print("-" * 48)

    errors = []
    approx_matrices = {}
    for rank in (1, 2, 3):
        A_k = (U[:, :rank] * s[:rank]) @ Vh[:rank, :]
        err = float(np.linalg.norm(A - A_k, ord="fro") / np.linalg.norm(A, ord="fro"))
        errors.append(err)
        approx_matrices[rank] = A_k
        print(f"  Rank-{rank} approx  Frobenius relative error: {err:.6f}")

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_svd.svg"

    vmin, vmax = A.min(), A.max()
    fig, axes = plt.subplots(2, 3, figsize=(13, 8))
    fig.suptitle("IMSL SVD: Low-Rank Approximations of a 4×6 Matrix", fontsize=12)

    # Row 0: original + rank-1 + rank-2
    for col, (title, mat) in enumerate([
        ("Original A", A),
        ("Rank-1 Approx", approx_matrices[1]),
        ("Rank-2 Approx", approx_matrices[2]),
    ]):
        ax = axes[0, col]
        im = ax.imshow(mat, vmin=vmin, vmax=vmax, cmap="RdYlBu_r", aspect="auto")
        ax.set_title(title, fontsize=10)
        ax.set_xlabel("Column")
        ax.set_ylabel("Row")
        fig.colorbar(im, ax=ax, shrink=0.8)

    # Row 1: rank-3 + error plot + singular values bar chart
    ax = axes[1, 0]
    im = ax.imshow(approx_matrices[3], vmin=vmin, vmax=vmax, cmap="RdYlBu_r", aspect="auto")
    ax.set_title("Rank-3 Approx", fontsize=10)
    ax.set_xlabel("Column")
    ax.set_ylabel("Row")
    fig.colorbar(im, ax=ax, shrink=0.8)

    ax = axes[1, 1]
    ax.bar([1, 2, 3], errors, color="#0369a1", edgecolor="#075985", alpha=0.85)
    ax.set_xticks([1, 2, 3])
    ax.set_xticklabels(["Rank-1", "Rank-2", "Rank-3"])
    ax.set_ylabel("Relative Frobenius error")
    ax.set_title("Reconstruction Error vs Rank", fontsize=10)
    ax.set_ylim(0, max(errors) * 1.2)

    ax = axes[1, 2]
    idx = np.arange(1, len(s) + 1)
    ax.bar(idx, s, color="#be123c", edgecolor="#9f1239", alpha=0.85)
    ax.set_xticks(idx)
    ax.set_xticklabels([f"σ{i}" for i in idx])
    ax.set_ylabel("Singular value")
    ax.set_title("Singular Values", fontsize=10)

    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)
    print(f"\nPlot saved: {plot_path}")

    return {
        "singular_values": s,
        "errors": errors,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_svd()

Plot Output

Top row: original 4×6 matrix and rank-1, rank-2 approximations. Bottom row: rank-3 approximation, relative Frobenius error per rank, and bar chart of all singular values.

Console Output

IMSL SVD Example: Low-Rank Approximation
Matrix shape: (4, 6)   (rank Γëê 3)
------------------------------------------------
Index        Singular value   Energy %
------------------------------------------------
  σ1               5.545273      71.68%
  σ2               3.156083      23.22%
  σ3               1.472810       5.06%
  σ4               0.144370       0.05%
------------------------------------------------
  Rank-1 approx  Frobenius relative error: 0.532195
  Rank-2 approx  Frobenius relative error: 0.225938
  Rank-3 approx  Frobenius relative error: 0.022042

Plot saved: test_output\example_imsl_svd.svg