Part IV: Advanced Topics | Chapter 3

Numerical Linear Algebra

Algorithms, stability, and the art of computing with matrices

Historical Context

Numerical linear algebra was born from the computing revolution. Alan Turing analyzed the stability of Gaussian elimination (1948), James Wilkinson developed backward error analysis at the National Physical Laboratory, and Gene Golub pioneered computational SVD algorithms. The LINPACK (1979) and LAPACK (1992) libraries codified best practices, while iterative methods (Lanczos 1950, Hestenes-Stiefel conjugate gradients 1952) enabled solving systems too large for direct methods. Today, numerical linear algebra underpins scientific computing, machine learning, and every application of matrix computation.

4.1 Floating-Point Arithmetic

IEEE 754 double precision uses 64 bits: 1 sign, 11 exponent, 52 mantissa, giving machine epsilon $\epsilon_{\text{mach}} \approx 2.2 \times 10^{-16}$. Every floating-point operation introduces relative error bounded by $\epsilon_{\text{mach}}$: $\text{fl}(x \oplus y) = (x \oplus y)(1 + \delta)$with $|\delta| \leq \epsilon_{\text{mach}}$.

Catastrophic cancellation occurs when subtracting nearly equal numbers, amplifying relative error. Algorithms must be designed to avoid such operations—this is the art of numerical stability.

4.2 Condition Numbers and Stability

Condition Number

The condition number $\kappa(A) = \|A\|\|A^{-1}\| = \sigma_{\max}/\sigma_{\min}$ bounds the worst-case amplification of input perturbations:$\frac{\|\delta x\|}{\|x\|} \leq \kappa(A) \frac{\|\delta b\|}{\|b\|}$.

An algorithm is backward stable if the computed solution is the exact solution of a slightly perturbed problem. Gaussian elimination with partial pivoting, QR factorization via Householder, and the SVD algorithm are all backward stable.

4.3 Iterative Methods

For large sparse systems, direct methods are too expensive. Iterative methods produce a sequence $x^{(k)} \to x^*$:

Jacobi: $x^{(k+1)} = D^{-1}(b - (L+U)x^{(k)})$, converges for diagonally dominant matrices
Gauss-Seidel: Uses updated components immediately, converges faster
Conjugate gradients: Optimal for SPD systems, converges in $\leq n$ steps exactly
GMRES: Generalized minimum residual for non-symmetric systems

4.4 Eigenvalue Algorithms

Power iteration: Finds dominant eigenvalue, convergence rate $|\lambda_2/\lambda_1|^k$
Inverse iteration: Finds eigenvalue nearest a target $\mu$
QR algorithm: Computes all eigenvalues via $A_k = R_kQ_k$ iterations, $O(n^3)$
Lanczos: For large sparse symmetric matrices, builds a tridiagonal approximation

4.5 Sparse Matrix Techniques

Sparse matrices store only nonzero entries in formats like CSR (compressed sparse row) and CSC (compressed sparse column). For a matrix with $\text{nnz}$ nonzeros, storage is $O(\text{nnz})$ vs $O(n^2)$ for dense. Sparse direct solvers use fill-reducing orderings (Cuthill-McKee, nested dissection) to minimize fill-in during factorization.

Computational Laboratory

This simulation demonstrates floating-point pitfalls, condition number analysis, iterative solver convergence, the QR eigenvalue algorithm, and sparse vs dense performance comparison.

Numerical Linear Algebra: Algorithms and Stability

Python

numerical_la.py175 lines

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

# ============================================================
# Numerical Linear Algebra
# ============================================================
np.random.seed(42)

# --- 1. Floating-point arithmetic ---
print("=" * 60)
print("FLOATING-POINT ARITHMETIC")
print("=" * 60)

eps = np.finfo(float).eps
print(f"Machine epsilon: {eps}")
print(f"1.0 + eps/2 == 1.0? {1.0 + eps/2 == 1.0}")
print(f"1.0 + eps   == 1.0? {1.0 + eps == 1.0}")

# Catastrophic cancellation
a = 1e15
b = 1e15 + 1
print(f"\nCatastrophic cancellation:")
print(f"(a+1) - a where a = 10¹⁵: {(a + 1) - a} (should be 1)")
print(f"(a+0.5) - a where a = 10¹⁵: {(a + 0.5) - a} (should be 0.5)")

# --- 2. Condition numbers ---
print("\n" + "=" * 60)
print("CONDITION NUMBERS")
print("=" * 60)

# Hilbert matrix (notoriously ill-conditioned)
for n in [3, 5, 8, 10, 12]:
    H = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])
    cond = np.linalg.cond(H)
    x_true = np.ones(n)
    b = H @ x_true
    x_computed = np.linalg.solve(H, b)
    rel_error = np.linalg.norm(x_computed - x_true) / np.linalg.norm(x_true)
    print(f"Hilbert {n}×{n}: cond = {cond:.2e}, relative error = {rel_error:.2e}")

# --- 3. Iterative methods ---
print("\n" + "=" * 60)
print("ITERATIVE METHODS")
print("=" * 60)

n = 50
# Diagonally dominant matrix
A = np.random.randn(n, n) * 0.1
np.fill_diagonal(A, np.abs(A).sum(axis=1) + 1)
b = np.random.randn(n)
x_exact = np.linalg.solve(A, b)

# Jacobi iteration
D_inv = np.diag(1.0 / np.diag(A))
R = A - np.diag(np.diag(A))
x = np.zeros(n)
jacobi_errors = []
for k in range(100):
    x = D_inv @ (b - R @ x)
    jacobi_errors.append(np.linalg.norm(x - x_exact) / np.linalg.norm(x_exact))

# Gauss-Seidel
x = np.zeros(n)
gs_errors = []
for k in range(100):
    for i in range(n):
        x[i] = (b[i] - A[i, :i] @ x[:i] - A[i, i+1:] @ x[i+1:]) / A[i, i]
    gs_errors.append(np.linalg.norm(x - x_exact) / np.linalg.norm(x_exact))

print(f"Jacobi:       converged to {jacobi_errors[-1]:.2e} in 100 iters")
print(f"Gauss-Seidel: converged to {gs_errors[-1]:.2e} in 100 iters")

# --- 4. Power iteration ---
print("\n" + "=" * 60)
print("EIGENVALUE ALGORITHMS")
print("=" * 60)

M = np.random.randn(20, 20)
M = M + M.T  # symmetric
true_evals = np.sort(np.abs(np.linalg.eigvalsh(M)))[::-1]

# Power iteration
x = np.random.randn(20)
power_evals = []
for k in range(30):
    x = M @ x
    x = x / np.linalg.norm(x)
    rayleigh = x @ M @ x
    power_evals.append(rayleigh)

print(f"Dominant eigenvalue (true): {true_evals[0]:.6f}")
print(f"Power iteration (30 iters): {power_evals[-1]:.6f}")

# QR algorithm (simplified)
A_qr = M.copy()
qr_evals_hist = []
for k in range(50):
    Q, R = np.linalg.qr(A_qr)
    A_qr = R @ Q
    qr_evals_hist.append(np.sort(np.diag(A_qr))[::-1][:3])

print(f"QR algorithm top 3 eigenvalues: {np.round(qr_evals_hist[-1], 4)}")
print(f"True top 3: {np.round(np.sort(np.linalg.eigvalsh(M))[::-1][:3], 4)}")

# --- 5. Sparse vs dense timing ---
print("\n" + "=" * 60)
print("SPARSE vs DENSE (timing)")
print("=" * 60)

from scipy import sparse
from scipy.sparse.linalg import spsolve

for n in [100, 500, 1000]:
    # Tridiagonal (sparse)
    diags = [np.ones(n-1), -2*np.ones(n), np.ones(n-1)]
    A_sp = sparse.diags(diags, [-1, 0, 1], format='csc')
    A_dense = A_sp.toarray()
    b = np.random.randn(n)

t0 = time.perf_counter()
    for _ in range(10): x_sp = spsolve(A_sp, b)
    t_sparse = (time.perf_counter() - t0) / 10

t0 = time.perf_counter()
    for _ in range(10): x_dn = np.linalg.solve(A_dense, b)
    t_dense = (time.perf_counter() - t0) / 10

print(f"n={n:4d}: sparse={t_sparse*1000:.3f}ms, dense={t_dense*1000:.3f}ms, ratio={t_dense/t_sparse:.1f}x")

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

# Plot 1: Iterative convergence
ax = axes[0]
ax.semilogy(jacobi_errors[:50], color="#3b82f6", linewidth=2, label="Jacobi")
ax.semilogy(gs_errors[:50], color="#0ea5e9", linewidth=2, label="Gauss-Seidel")
ax.set_xlabel("Iteration", color="gray")
ax.set_ylabel("Relative Error", color="gray")
ax.set_title("Iterative Solver Convergence", fontsize=12, fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_facecolor("#0f172a")

# Plot 2: Condition number growth
ax = axes[1]
ns = list(range(2, 15))
conds = [np.linalg.cond(np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])) for n in ns]
ax.semilogy(ns, conds, "o-", color="#3b82f6", markersize=5)
ax.set_xlabel("n", color="gray")
ax.set_ylabel("κ(H_n)", color="gray")
ax.set_title("Hilbert Matrix Condition", fontsize=12, fontweight="bold")
ax.grid(True, alpha=0.3)
ax.set_facecolor("#0f172a")

# Plot 3: QR algorithm convergence
ax = axes[2]
qr_hist = np.array(qr_evals_hist)
for i, color in enumerate(["#3b82f6", "#0ea5e9", "#06b6d4"]):
    ax.plot(qr_hist[:, i], color=color, linewidth=2, label=f"λ_{i+1}")
    true_val = np.sort(np.linalg.eigvalsh(M))[::-1][i]
    ax.axhline(y=true_val, color=color, linestyle="--", alpha=0.5)
ax.set_xlabel("QR Iteration", color="gray")
ax.set_title("QR Algorithm Convergence", fontsize=12, fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
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

← Intro to Functional Analysis Applications in Data Science →