Numerical Linear Algebra

Numerical Linear Algebra

Numerical linear algebra is the bridge between elegant mathematical theory and practical computation. In theory, we solve $Ax = b$ by computing $x = A^{-1}b$. In practice, computing $A^{-1}$ directly is almost never done — it is slow, numerically unstable, and wasteful. Instead, we use carefully designed algorithms that exploit matrix structure and maintain numerical stability. This module covers the essential numerical concepts every ML practitioner must understand: conditioning, stability, floating-point arithmetic, direct and iterative solvers, and their applications in modern AI systems.

1. Why It Matters

Computational Efficiency

Naive matrix inversion via cofactor expansion is $O(n!)$ — completely impractical. Even the adjugate method is $O(n^4)$. Modern algorithms achieve:

For sparse matrices with $k$ nonzero entries, iterative methods reduce complexity to $O(k)$ per iteration, making problems with millions of variables tractable.

The Cost of Ignoring Numerics

import numpy as np

# Hilbert matrix: H_ij = 1/(i+j-1)
# Theoretical solution: x = [1, 1, ..., 1]
n = 12
H = np.array([[1/(i+j-1) for j in range(n)] for i in range(n)])
b = H @ np.ones(n)

# "Solve" Ax = b
x_numpy = np.linalg.solve(H, b)
print(f"Error (np.linalg.solve): {np.linalg.norm(x_numpy - np.ones(n)):.2e}")
# ~1e-4 (terrible for a 12x12 system!)

# The condition number tells you WHY
print(f"Condition number: {np.linalg.cond(H):.2e}")
# ~10^16 — the matrix is nearly singular in floating point

2. Condition Number

What the Condition Number Tells You

Sensitivity Analysis

For the system $Ax = b$, if $b$ is perturbed by $\delta b$, the relative error in $x$ satisfies:

Similarly, if $A$ is perturbed by $\delta A$:

This is a worst-case bound — the actual error may be much smaller, but you cannot know this without computing $\kappa(A)$.

import numpy as np
import matplotlib.pyplot as plt

# Condition number growth for Hilbert matrices
ns = range(2, 15)
kappas = []
for n in ns:
    H = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])
    kappas.append(np.linalg.cond(H))

for n, k in zip(ns, kappas):
    print(f"H_{n}: κ = {k:.2e}")

# For n >= 12, κ exceeds 10^16 — beyond double precision capability

3. Stability

Stability of Common Algorithms

import numpy as np

def gaussian_elimination_no_pivot(A, b):
    """Unstable Gaussian elimination (no pivoting)."""
    n = len(b)
    Ab = np.hstack([A.astype(float), b.reshape(-1, 1)])
    for col in range(n):
        for row in range(col + 1, n):
            factor = Ab[row, col] / Ab[col, col]
            Ab[row, col:] -= factor * Ab[col, col:]
    # Back substitution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - Ab[i, i+1:n] @ x[i+1:]) / Ab[i, i]
    return x

# Near-singular example
A = np.array([[1e-10, 1], [1, 1]], dtype=float)
b = A @ np.array([1.0, 1.0])

# numpy uses stable LAPACK (LU with partial pivoting)
x_stable = np.linalg.solve(A, b)
print(f"Stable: {x_stable} (error: {np.linalg.norm(x_stable - [1,1]):.2e})")

# No-pivot version may fail
try:
    x_unstable = gaussian_elimination_no_pivot(A, b)
    print(f"Unstable: {x_unstable} (error: {np.linalg.norm(x_unstable - [1,1]):.2e})")
except ZeroDivisionError:
    print("Unstable: Division by zero!")

4. Floating Point Arithmetic

IEEE 754 Standard

Roundoff Error

Every floating-point operation introduces a relative error of at most $\epsilon$. After $n$ operations, errors can accumulate. For a sum of $n$ numbers:

This linear growth in error is why compensated summation (Kahan summation) is important for large reductions.

Catastrophic Cancellation

import numpy as np

# Classic catastrophic cancellation
a = 1.0
b = 1.0 + 1e-15
print(f"Direct subtraction: {(b - a) - 1e-15}")  # May not be exactly 0

# Catastrophic cancellation in quadratic formula
# For b >> sqrt(b^2 - 4ac), the standard formula loses precision
a_coeff, b_coeff, c_coeff = 1, 1e8, 1
sqrt_disc = np.sqrt(b_coeff**2 - 4*a_coeff*c_coeff)

x1_standard = (-b_coeff + sqrt_disc) / (2 * a_coeff)
x2_standard = (-b_coeff - sqrt_disc) / (2 * a_coeff)
print(f"Standard: x1 = {x1_standard:.6e}, x2 = {x2_standard:.6e}")

# Stable alternative: use the other root to avoid cancellation
x2_stable = c_coeff / (a_coeff * x1_standard)  # = c/(a*x1)
x1_stable = (-b_coeff - np.sign(b_coeff) * sqrt_disc) / (2 * a_coeff)
print(f"Stable:   x1 = {x1_stable:.6e}, x2 = {x2_stable:.6e}")

Kahan Summation

def kahan_sum(values):
    """Kahan compensated summation — nearly exact for large sums."""
    s = 0.0
    c = 0.0  # compensation for lost low-order bits
    for x in values:
        y = x - c
        t = s + y
        c = (t - s) - y  # algebraically zero, but captures rounding error
        s = t
    return s

# Demonstration: sum 1e7 values of 1.0
n = 10_000_000
values = [1.0] * n
print(f"Naive sum:   {sum(values) - n}")       # May show rounding error
print(f"Kahan sum:   {kahan_sum(values) - n}")  # Exact (or nearly so)

# More dramatic: sum decreasing values
vals = [1e16, 1.0, -1e16]
print(f"Naive:   {sum(vals)}")          # 0.0 (correct by luck)
vals2 = [1e16, 1.0, -1e16, 1.0, -1e16]
print(f"Naive:   {sum(vals2)}")         # 0.0 (lost the 1.0!)
print(f"Kahan:   {kahan_sum(vals2)}")   # 1.0 (correct)

5. Direct Methods

Gaussian Elimination

LU Decomposition

import numpy as np
from scipy.linalg import lu, solve_triangular

A = np.array([[2, 1, 1], [4, 3, 3], [8, 7, 9]], dtype=float)
b = np.array([1, 2, 3], dtype=float)

# LU decomposition via scipy
P, L, U = lu(A)

# Verify: PA = LU
print("PA = LU check:", np.allclose(P @ A, L @ U))

# Solve using LU factors
y = solve_triangular(L, P @ b, lower=True)
x = solve_triangular(U, y)
print(f"Solution: {x}")
print(f"Verify:   {A @ x}")  # Should equal b

# numpy's solve uses LAPACK's LU implementation internally
x_numpy = np.linalg.solve(A, b)
print(f"numpy:    {x_numpy}")

Cholesky Decomposition

import numpy as np
from scipy.linalg import cho_factor, cho_solve

# Create SPD matrix
A = np.array([[4, 2, 1], [2, 5, 3], [1, 3, 6]], dtype=float)
b = np.array([1, 2, 3], dtype=float)

# Cholesky decomposition
L = np.linalg.cholesky(A)
print("A = LL^T check:", np.allclose(L @ L.T, A))

# Solve using Cholesky
c, low = cho_factor(A)
x = cho_solve((c, low), b)
print(f"Solution: {x}")
print(f"Verify:   {A @ x}")

6. Iterative Methods

Iterative methods are essential for large-scale problems where direct methods are too expensive (e.g., finite element systems with millions of unknowns, or sparse systems arising in graph algorithms).

Jacobi Iteration

Gauss-Seidel Iteration

Conjugate Gradient

Iterative Methods Comparison

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import cg, gmres

# Create a large SPD tridiagonal system
n = 10000
A = diags([-1, 2, -1], [-1, 0, 1], shape=(n, n), format='csr')
b = np.ones(n)

# Conjugate Gradient
x_cg, info = cg(A, b, tol=1e-10, maxiter=1000)
print(f"CG: {np.linalg.norm(A @ x_cg - b):.2e} residual, info={info}")

# Jacobi iteration (manual)
def jacobi(A, b, x0=None, tol=1e-10, maxiter=1000):
    x = x0 if x0 is not None else np.zeros_like(b)
    D = A.diagonal()
    for k in range(maxiter):
        r = b - A @ x
        if np.linalg.norm(r) < tol:
            return x, k
        x = x + r / D  # Jacobi step: x_new = D^{-1}(b - (A-D)x)
    return x, maxiter

x_jac, iters = jacobi(A, b)
print(f"Jacobi: {np.linalg.norm(A @ x_jac - b):.2e} residual, {iters} iterations")

7. Matrix Norms for Numerical Analysis

8. Least Squares Problems

When $Ax = b$ has no exact solution (overdetermined system, $m > n$), we seek $x$ that minimizes $\|Ax - b\|_2$.

Normal Equations

QR Decomposition Approach

SVD Approach

import numpy as np

# Overdetermined system: 5 equations, 3 unknowns
np.random.seed(42)
A = np.random.randn(5, 3)
b = np.random.randn(5)

# Method 1: Normal equations (UNSTABLE)
x_normal = np.linalg.solve(A.T @ A, A.T @ b)
print(f"Normal equations: residual = {np.linalg.norm(A @ x_normal - b):.2e}")
print(f"  κ(A^T A) = {np.linalg.cond(A.T @ A):.2e}")

# Method 2: QR decomposition (STABLE)
x_qr, res_qr, rank_qr, sv_qr = np.linalg.lstsq(A, b, rcond=None)
print(f"QR (lstsq):       residual = {np.linalg.norm(A @ x_qr - b):.2e}")

# Method 3: SVD (MOST STABLE)
U, s, Vt = np.linalg.svd(A, full_matrices=False)
x_svd = Vt.T @ np.diag(1/s) @ U.T @ b
print(f"SVD:              residual = {np.linalg.norm(A @ x_svd - b):.2e}")

# Rank-deficient example
A_rankdef = np.array([[1, 2], [2, 4], [3, 6]], dtype=float)  # rank 1
b_rankdef = np.array([1, 2, 3.1], dtype=float)
x_rd = np.linalg.lstsq(A_rankdef, b_rankdef, rcond=None)[0]
print(f"\nRank-deficient least squares: x = {x_rd}")
print(f"  Residual: {np.linalg.norm(A_rankdef @ x_rd - b_rankdef):.2e}")

9. Sparse Matrix Algorithms

Why Sparse?

Sparse Storage Formats

Preconditioning

import numpy as np
from scipy.sparse import diags, random
from scipy.sparse.linalg import cg, spilu, LinearOperator

# Create a difficult SPD system (convection-diffusion)
n = 500
A = diags([-1, 4, -1], [-1, 0, 1], shape=(n, n), format='csr')
A = A + 0.1 * random(n, n, density=0.01, format='csr')  # Add some non-symmetric perturbation
A = A @ A.T  # Make SPD
b = np.random.randn(n)

# Without preconditioning
x_no_precon, info1 = cg(A, b, tol=1e-10, maxiter=1000)
print(f"No preconditioner: {np.linalg.norm(A @ x_no_precon - b):.2e}")

# With Jacobi (diagonal) preconditioner
M_inv = 1.0 / A.diagonal()
M = LinearOperator(A.shape, matvec=lambda x: M_inv * x)
x_jac, info2 = cg(A, b, tol=1e-10, maxiter=1000, M=M)
print(f"Jacobi precon:     {np.linalg.norm(A @ x_jac - b):.2e}")

# With incomplete Cholesky
ilu = spilu(A.tocsc(), drop_tol=1e-4)
M_ilu = LinearOperator(A.shape, matvec=ilu.solve)
x_ilu, info3 = cg(A, b, tol=1e-10, maxiter=1000, M=M_ilu)
print(f"ILU precon:        {np.linalg.norm(A @ x_ilu - b):.2e}")

10. Python Implementation

NumPy Linear Algebra

import numpy as np

# === Core Operations ===

# Solve Ax = b
A = np.array([[3, 1], [1, 2]], dtype=float)
b = np.array([9, 8], dtype=float)
x = np.linalg.solve(A, b)
print(f"Solve: x = {x}")

# Determinant and inverse (AVOID inverting in production!)
print(f"det(A) = {np.linalg.det(A):.4f}")
print(f"inv(A) = \n{np.linalg.inv(A)}")

# Matrix factorizations
Q, R = np.linalg.qr(A)
print(f"\nQR: Q = \n{Q}\nR = \n{R}")

L = np.linalg.cholesky(A @ A.T + np.eye(2))  # Ensure SPD
print(f"Cholesky: L = \n{L}")

U, s, Vt = np.linalg.svd(A)
print(f"\nSVD: singular values = {s}")

# Eigenvalues
eigvals, eigvecs = np.linalg.eig(A)
print(f"Eigenvalues: {eigvals}")

# Matrix norms
print(f"\n||A||_F = {np.linalg.norm(A, 'fro'):.4f}")
print(f"||A||_2 = {np.linalg.norm(A, 2):.4f}")
print(f"||A||_1 = {np.linalg.norm(A, 1):.4f}")
print(f"||A||_inf = {np.linalg.norm(A, np.inf):.4f}")

# Condition number
print(f"κ(A) = {np.linalg.cond(A):.4f}")
print(f"κ(A^T A) = {np.linalg.cond(A.T @ A):.4f}")  # κ^2

SciPy Sparse Linear Algebra

from scipy import sparse
from scipy.sparse import linalg as spla
import numpy as np

# Create sparse matrix (5-point Laplacian)
n = 100
A = sparse.diags([-1, 4, -1], [-1, 0, 1], shape=(n, n), format='csr')
b = np.random.randn(n)

# Solve sparse system
x, info = spla.cg(A, b, tol=1e-10)
print(f"CG solve: residual = {np.linalg.norm(A @ x - b):.2e}")

# Sparse eigenvalues (largest/smallest)
vals, vecs = spla.eigs(A, k=5, which='LM')
print(f"Largest eigenvalues: {np.sort(vals.real)}")

vals_sm, _ = spla.eigs(A, k=3, which='SM')
print(f"Smallest eigenvalues: {np.sort(vals_sm.real)}")

# Sparse SVD (useful for large-scale PCA)
U, s, Vt = spla.svds(A, k=5)
print(f"Top 5 singular values: {s}")

Practical Patterns

11. Applications in AI/ML

Training Stability

Numerical Precision in Deep Learning

import numpy as np

# Simulating mixed-precision issues
def forward_pass_fp32(W, x):
    """Full precision forward pass."""
    return W @ x

def forward_pass_fp16(W, x):
    """Simulated mixed-precision forward pass."""
    W_16 = W.astype(np.float16)
    x_16 = x.astype(np.float16)
    result = W_16 @ x_16
    return result.astype(np.float32)

np.random.seed(42)
W = np.random.randn(100, 100) * 0.1
x = np.random.randn(100)

y_fp32 = forward_pass_fp32(W, x)
y_fp16 = forward_pass_fp16(W, x)

print(f"FP32 result norm: {np.linalg.norm(y_fp32):.6f}")
print(f"FP16 result norm: {np.linalg.norm(y_fp16):.6f}")
print(f"Relative error:   {np.linalg.norm(y_fp32 - y_fp16) / np.linalg.norm(y_fp32):.2e}")

# Condition number of a typical layer's weight matrix
W_layer = np.random.randn(512, 512) * (2 / 512)**0.5  # Xavier init
print(f"\nLayer condition number: {np.linalg.cond(W_layer):.2f}")
# Usually ~10-100 with good initialization, much worse after training

Key Applications

12. Common Mistakes

import numpy as np

# BAD: Using inv(A) @ b
A = np.array([[1, 2], [3, 4]], dtype=float)
b = np.array([5, 6], dtype=float)

# Wrong
x_bad = np.linalg.inv(A) @ b  # Slower, less stable

# Correct
x_good = np.linalg.solve(A, b)  # Direct, stable

print(f"inv @ b: {x_bad}")
print(f"solve:  {x_good}")
print(f"Equal:  {np.allclose(x_bad, x_good)}")

# BAD: Comparing floats with ==
a = 0.1 + 0.2
print(f"\n0.1 + 0.2 == 0.3? {a == 0.3}")          # False!
print(f"0.1 + 0.2 ≈ 0.3?  {np.isclose(a, 0.3)}")  # True

13. Interview Questions

14. Practice Problems

import numpy as np

for n in [5, 10, 15, 20]:
    H = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])
    kappa = np.linalg.cond(H)
    digits_lost = int(np.log10(kappa))
    print(f"n={n:2d}: κ = {kappa:.2e}, ~{digits_lost} digits lost, "
          f"{max(0, 16-digits_lost)} digits remaining")

import numpy as np

def backward_error(A, b, x_approx):
    """Compute the backward error of an approximate solution."""
    residual = b - A @ x_approx
    return np.linalg.norm(residual) / (np.linalg.norm(A) * np.linalg.norm(x_approx))

# Ill-conditioned system
np.random.seed(42)
n = 100
H = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])
b = H @ np.ones(n)  # True solution is x = [1,...,1]

x_solve = np.linalg.solve(H, b)
x_perturbed = x_solve + np.random.randn(n) * 1e-10  # Slightly wrong solution

print(f"np.linalg.solve backward error: {backward_error(H, b, x_solve):.2e}")
print(f"Perturbed solution backward error: {backward_error(H, b, x_perturbed):.2e}")
print(f"Solve error: {np.linalg.norm(x_solve - np.ones(n)):.2e}")
print(f"Perturbed error: {np.linalg.norm(x_perturbed - np.ones(n)):.2e}")

import numpy as np

A = np.array([[4, 1], [1, 3]], dtype=float)
D = np.diag(np.diag(A))
L = np.tril(A, -1)
U = np.triu(A, 1)

T_J = -np.linalg.solve(D, L + U)
T_GS = -np.linalg.solve(D + L, U)

print(f"ρ(T_Jacobi) = {max(abs(np.linalg.eigvals(T_J))):.4f}")
print(f"ρ(T_GS)     = {max(abs(np.linalg.eigvals(T_GS))):.4f}")
print(f"Ratio:       {max(abs(np.linalg.eigvals(T_J))) / max(abs(np.linalg.eigvals(T_GS))):.2f}x")

15. Quick Reference

16. Cross-References

• Linear Algebra Basics — Vector spaces, subspaces, and the four fundamental subspaces
• Linear Algebra SVD — SVD, eigendecomposition, QR, Cholesky (theory)
• Linear Algebra Eigenvalues — Eigenvalue computation and the QR algorithm
• Linear Algebra Applications — PCA, SVD embeddings, and matrix factorizations in ML
• Optimization Convex — Gradient descent, second-order methods, and Hessian computation
• Probability Distributions — Multivariate Gaussians require positive definite covariance matrices (Cholesky)
• Calculus Differential Equations — Backpropagation, batch normalization, and training dynamics