Part III: Decompositions | Chapter 3

Matrix Decompositions

LU, QR, Cholesky, and Schur: the computational workhorses of numerical linear algebra

Historical Context

Matrix decompositions form the backbone of computational linear algebra. Gaussian elimination (formalized as LU decomposition) dates to ancient Chinese mathematics but was systematized by Gauss for least squares. André-Louis Cholesky developed his decomposition for geodetic survey computations during World War I. Householder (1958) and Givens introduced numerically stable QR algorithms, while Issai Schur's triangularization theorem (1909) provided the theoretical foundation for modern eigenvalue algorithms.

3.1 LU Decomposition

Theorem: LU Factorization

For $A \in \mathbb{R}^{n \times n}$, there exist a permutation matrix $P$, unit lower triangular $L$, and upper triangular $U$ such that $PA = LU$. Without pivoting ($P = I$), the factorization exists when all leading principal minors are nonzero.

The LU decomposition costs $\frac{2}{3}n^3$ flops and enables solving $Ax = b$ in$O(n^2)$ via forward and back substitution. Partial pivoting (row exchanges) ensures numerical stability; complete pivoting additionally permutes columns but is rarely needed in practice.

3.2 QR Decomposition

Theorem: QR Factorization

Every $A \in \mathbb{R}^{m \times n}$ ($m \geq n$) with full column rank has a unique factorization $A = QR$ where $Q \in \mathbb{R}^{m \times n}$ has orthonormal columns and $R \in \mathbb{R}^{n \times n}$ is upper triangular with positive diagonal entries.

Three methods compute the QR decomposition:

Gram-Schmidt: Conceptually simple, $O(mn^2)$, but can lose orthogonality
Householder reflections: $O(2mn^2 - \frac{2}{3}n^3)$, excellent stability, the standard method
Givens rotations: Best for sparse matrices, zeroes one element at a time

3.3 Cholesky Decomposition

Theorem: Cholesky Factorization

If $A$ is symmetric positive definite, then there exists a unique lower triangular $L$with positive diagonal entries such that $A = LL^T$. The computation requires only$\frac{1}{3}n^3$ flops—half of LU.

Cholesky is the method of choice for solving symmetric positive definite systems (covariance matrices, normal equations, finite element stiffness matrices). It is also a convenient test for positive definiteness: the factorization exists if and only if the matrix is SPD.

3.4 Schur Decomposition

Theorem (Schur)

Every $A \in \mathbb{C}^{n \times n}$ can be written as $A = UTU^*$ where $U$ is unitary and $T$ is upper triangular with eigenvalues on the diagonal. For real matrices, the real Schur form uses an orthogonal $Q$ and quasi-upper-triangular $T$ (with $1 \times 1$ and $2 \times 2$ diagonal blocks).

The Schur decomposition is the numerically stable alternative to the Jordan form. The implicit QR algorithm computes it in $O(n^3)$ operations and is the basis of all practical eigenvalue computations. Unlike the Jordan form, the Schur form is stable under perturbation—a crucial property for floating-point arithmetic.

3.5 Selection Guide

Decomposition	Requirements	Cost	Best For
LU	Square, nonsingular	$\frac{2}{3}n^3$	General linear systems
QR	Any shape	$2mn^2$	Least squares, eigenvalues
Cholesky	SPD	$\frac{1}{3}n^3$	SPD systems, sampling
Schur	Square	$O(n^3)$	Eigenvalues, matrix functions
SVD	Any shape	$O(mn^2)$	Rank, pseudoinverse, PCA

Computational Laboratory

This simulation computes LU, QR, Cholesky, and Schur decompositions, verifies each factorization, compares performance, and visualizes the sparsity patterns of the resulting factors.

Matrix Decompositions: LU, QR, Cholesky, Schur

Python

matrix_decompositions.py112 lines

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.linalg import lu, schur
import time

# ============================================================
# Matrix Decompositions: LU, QR, Cholesky, Schur
# ============================================================

np.random.seed(42)

# --- 1. LU Decomposition ---
print("=" * 60)
print("LU DECOMPOSITION (with partial pivoting)")
print("=" * 60)

A = np.array([[2, 1, 1], [4, 3, 3], [8, 7, 9]], dtype=float)
P_lu, L, U = lu(A)
print(f"A:\n{A}")
print(f"\nP (permutation):\n{P_lu}")
print(f"L (lower):\n{np.round(L, 6)}")
print(f"U (upper):\n{np.round(U, 6)}")
print(f"PA = LU? {np.allclose(P_lu @ A, L @ U)}")

# Solve Ax = b via LU
b = np.array([1, 2, 3], dtype=float)
y = np.linalg.solve(np.tril(L), P_lu @ b)
x = np.linalg.solve(np.triu(U), y)
print(f"\nSolve Ax = [1,2,3]: x = {np.round(x, 6)}")
print(f"Verify: Ax = {np.round(A @ x, 6)}")

# --- 2. QR Decomposition ---
print("\n" + "=" * 60)
print("QR DECOMPOSITION")
print("=" * 60)

B = np.array([[1, 1, 0], [1, 0, 1], [0, 1, 1]], dtype=float)
Q, R = np.linalg.qr(B)
print(f"B:\n{B}")
print(f"Q (orthogonal):\n{np.round(Q, 6)}")
print(f"R (upper tri):\n{np.round(R, 6)}")
print(f"Q^TQ = I? {np.allclose(Q.T @ Q, np.eye(3))}")
print(f"QR = B? {np.allclose(Q @ R, B)}")

# --- 3. Cholesky Decomposition ---
print("\n" + "=" * 60)
print("CHOLESKY DECOMPOSITION")
print("=" * 60)

M = np.array([[4, 2, 1], [2, 5, 3], [1, 3, 6]], dtype=float)
L_chol = np.linalg.cholesky(M)
print(f"Positive definite M:\n{M}")
print(f"Eigenvalues: {np.round(np.linalg.eigvalsh(M), 4)} (all > 0)")
print(f"\nCholesky L:\n{np.round(L_chol, 6)}")
print(f"LL^T = M? {np.allclose(L_chol @ L_chol.T, M)}")

# --- 4. Schur Decomposition ---
print("\n" + "=" * 60)
print("SCHUR DECOMPOSITION")
print("=" * 60)

C = np.array([[1, 2, 0], [-1, 3, 1], [0, 1, 2]], dtype=float)
T_schur, Z = schur(C)
print(f"C:\n{C}")
print(f"T (quasi-upper triangular):\n{np.round(T_schur, 6)}")
print(f"Z^TZ = I? {np.allclose(Z.T @ Z, np.eye(3))}")
print(f"ZTZ^T = C? {np.allclose(Z @ T_schur @ Z.T, C)}")
print(f"Eigenvalues from T diagonal: {np.round(np.diag(T_schur), 4)}")
print(f"Eigenvalues from eig: {np.round(np.sort(np.linalg.eigvals(C).real), 4)}")

# --- 5. Performance comparison ---
print("\n" + "=" * 60)
print("DECOMPOSITION TIMING (100×100)")
print("=" * 60)

n = 100
A_large = np.random.randn(n, n)
A_spd = A_large.T @ A_large + n * np.eye(n)

for name, func, mat in [("LU", lambda m: lu(m), A_large),
                          ("QR", lambda m: np.linalg.qr(m), A_large),
                          ("Cholesky", lambda m: np.linalg.cholesky(m), A_spd),
                          ("Schur", lambda m: schur(m), A_large),
                          ("SVD", lambda m: np.linalg.svd(m), A_large)]:
    t0 = time.perf_counter()
    for _ in range(50):
        func(mat)
    elapsed = (time.perf_counter() - t0) / 50 * 1000
    print(f"  {name:10s}: {elapsed:.3f} ms")

# --- 6. Visualization ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Plot 1: Sparsity patterns
decomps = {"L (LU)": L, "Q (QR)": Q, "L (Chol)": L_chol}
for idx, (name, mat) in enumerate(decomps.items()):
    ax = axes[idx]
    im = ax.imshow(np.abs(mat), cmap="Blues", vmin=0)
    ax.set_title(name, fontsize=12, fontweight="bold")
    for i in range(mat.shape[0]):
        for j in range(mat.shape[1]):
            ax.text(j, i, f"{mat[i,j]:.2f}", ha="center", va="center",
                    fontsize=8, color="white" if abs(mat[i,j]) > 0.5 else "gray")
    ax.set_facecolor("#0f172a")

plt.tight_layout()
plt.savefig("output.png", dpi=150, bbox_inches="tight", facecolor="#0f172a")
plt.show()
print("\nVisualization saved.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Singular Value Decomposition Least Squares →