Cholesky Example — Matrix Decompositions

This example constructs a 4×4 symmetric positive-definite matrix and applies three factorisation routines from the eigensystems library:

• eigensystems.cholesky() — lower-triangular L such that A = L @ Lᵀ

• eigensystems.lu_decomp() — pivoted LU so that A = P @ L @ U

• eigensystems.qr_decomp() — QR so that A = Q @ R

Reconstruction errors confirm the numerical accuracy of each factorisation.

Example Code

"""IMSL Cholesky example: matrix decompositions on a symmetric positive-definite matrix.

Demonstrates:
- Create a 4x4 symmetric positive-definite matrix.
- Compute Cholesky decomposition via :func:`eigensystems.cholesky`.
- Compute LU decomposition via :func:`eigensystems.lu_decomp`.
- Compute QR decomposition via :func:`eigensystems.qr_decomp`.
- Verify reconstructions: L@L.T == A, P@L@U == A, Q@R == A.
- Plot heatmaps of A, L (Cholesky factor), and Q (from QR).
- Save to test_output/example_imsl_cholesky.svg

Outputs:
- Printed reconstruction errors for all three decompositions.
- SVG figure saved to test_output/.
- Console output captured to docs/_generated_text/cholesky_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 cholesky, lu_decomp, qr_decomp


def run_demo_cholesky() -> dict:
    """Run IMSL Cholesky/LU/QR example on a 4x4 symmetric positive-definite matrix.

    Constructs the matrix A = B @ B.T + 4*I to guarantee positive definiteness,
    then validates each factorisation by measuring the Frobenius reconstruction
    error.

    Args:
        None

    Returns:
        dict: Result dictionary with keys ``chol_error`` (float),
            ``lu_error`` (float), ``qr_error`` (float), and
            ``plot_path`` (str).
    """
    B = np.array([
        [4.0, 1.0, 0.5, 0.2],
        [1.0, 3.0, 0.8, 0.3],
        [0.5, 0.8, 2.5, 0.6],
        [0.2, 0.3, 0.6, 2.0],
    ])
    # Ensure positive definiteness
    A = B @ B.T + 4.0 * np.eye(4)

    # --- Cholesky ---
    L_chol = cholesky(A)
    chol_error = float(np.linalg.norm(A - L_chol @ L_chol.T, ord="fro"))

    # --- LU ---
    P_lu, L_lu, U_lu = lu_decomp(A)
    lu_error = float(np.linalg.norm(A - P_lu @ L_lu @ U_lu, ord="fro"))

    # --- QR ---
    Q_qr, R_qr = qr_decomp(A)
    qr_error = float(np.linalg.norm(A - Q_qr @ R_qr, ord="fro"))

    print("\nIMSL Cholesky / LU / QR Example: 4×4 Symmetric Positive-Definite Matrix")
    print("-" * 60)
    print(f"{'Decomposition':<20} {'Verification':>28} {'Frob. Error':>12}")
    print("-" * 60)
    print(f"  Cholesky          A = L @ L.T              {chol_error:>12.2e}")
    print(f"  LU (pivoted)      A = P @ L @ U            {lu_error:>12.2e}")
    print(f"  QR                A = Q @ R                {qr_error:>12.2e}")
    print("-" * 60)

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

    fig, axes = plt.subplots(1, 3, figsize=(13, 4))
    fig.suptitle(
        "IMSL Cholesky / LU / QR: Decomposition Factors of a 4×4 SPD Matrix",
        fontsize=11,
    )

    datasets = [
        ("Original A", A),
        ("Cholesky L\n(A = L @ Lᵀ)", L_chol),
        ("QR factor Q\n(A = Q @ R)", Q_qr),
    ]
    for ax, (title, mat) in zip(axes, datasets):
        im = ax.imshow(mat, cmap="coolwarm", aspect="auto")
        ax.set_title(title, fontsize=10)
        ax.set_xlabel("Column")
        ax.set_ylabel("Row")
        for r in range(mat.shape[0]):
            for c in range(mat.shape[1]):
                ax.text(c, r, f"{mat[r, c]:.2f}", ha="center", va="center",
                        fontsize=7, color="black")
        fig.colorbar(im, ax=ax, shrink=0.8)

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

    return {
        "chol_error": chol_error,
        "lu_error": lu_error,
        "qr_error": qr_error,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_cholesky()

Plot Output

Left to right: original SPD matrix A, lower Cholesky factor L, and orthogonal factor Q from the QR decomposition. Cell values are annotated directly on each heatmap.

Console Output

IMSL Cholesky / LU / QR Example: 4×4 Symmetric Positive-Definite Matrix
------------------------------------------------------------
Decomposition                        Verification  Frob. Error
------------------------------------------------------------
  Cholesky          A = L @ L.T                  4.21e-15
  LU (pivoted)      A = P @ L @ U                9.93e-16
  QR                A = Q @ R                    3.27e-15
------------------------------------------------------------

Plot saved: test_output\example_imsl_cholesky.svg