Solving Linear Systems

The backbone of scientific computing: solving \(Ax = b\) efficiently, stably, and at scale. Direct methods for dense systems, iterative methods for sparse ones.

Why Linear Systems Are Everywhere

Almost every numerical method ultimately reduces to solving a linear system. Finite element methods, implicit ODE solvers, optimization algorithms, least-squares fitting -- all require solving \(Ax = b\) repeatedly. The cost of solving this system often dominates the total computation time.

Dense: \(O(N^3)\)

Gaussian elimination, LU, Cholesky

Sparse Iterative: \(O(N \cdot k)\)

Jacobi, Gauss-Seidel, CG

Condition Number

\(\kappa(A) = \|A\| \cdot \|A^{-1}\|\)

1. Gaussian Elimination

The fundamental direct method. Transform \(Ax = b\) to upper-triangular form \(Ux = c\)using row operations, then back-substitute. Cost: \(\frac{2}{3}N^3\) flops.

Partial Pivoting

Never divide by small numbers. At each step, swap the current row with the row having the largest absolute value in the pivot column. This prevents catastrophic growth of round-off errors and is essential for numerical stability. Without pivoting, even well-conditioned systems can produce garbage results.

Gaussian Elimination with Partial Pivoting

Python

script.py45 lines

import numpy as np

def gauss_eliminate(A, b):
    """Gaussian elimination with partial pivoting."""
    n = len(b)
    # Augmented matrix
    Ab = np.hstack([A.astype(float), b.reshape(-1, 1).astype(float)])

for k in range(n):
        # Partial pivoting: find max in column k
        max_row = k + np.argmax(np.abs(Ab[k:, k]))
        Ab[[k, max_row]] = Ab[[max_row, k]]

if abs(Ab[k, k]) < 1e-14:
            raise ValueError("Matrix is singular or nearly singular")

# Eliminate below
        for i in range(k + 1, n):
            factor = Ab[i, k] / Ab[k, k]
            Ab[i, k:] -= factor * Ab[k, k:]

# Back substitution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - np.dot(Ab[i, i+1:n], x[i+1:n])) / Ab[i, i]

return x

# Test: solve a 4x4 system
A = np.array([
    [2, 1, -1, 1],
    [4, 5, -3, 5],
    [-2, 5, -2, 6],
    [4, 11, -4, 8]
], dtype=float)
b = np.array([5, 9, 4, 2], dtype=float)

x = gauss_eliminate(A, b)
x_numpy = np.linalg.solve(A, b)

print("Gaussian Elimination with Partial Pivoting")
print(f"Our solution:   x = {x}")
print(f"NumPy solution: x = {x_numpy}")
print(f"Max difference: {np.max(np.abs(x - x_numpy)):.2e}")
print(f"Residual ||Ax - b||: {np.linalg.norm(A @ x - b):.2e}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

2. LU Decomposition

Factor \(A = LU\) where \(L\) is lower-triangular and \(U\) is upper-triangular. Then solve \(Ax = b\) in two steps: \(Ly = b\) (forward substitution) and \(Ux = y\) (back substitution). The key advantage: factorize once \((O(N^3))\), solve for multiple right-hand sides cheaply \((O(N^2)\) each).

\(PA = LU\)

\(P\) = permutation matrix (from pivoting), \(L\) = lower triangular with 1s on diagonal, \(U\) = upper triangular

LU Decomposition with SciPy

Python

script.py37 lines

import numpy as np
from scipy import linalg

# Create a test matrix
A = np.array([
    [2, 1, 1, 0],
    [4, 3, 3, 1],
    [8, 7, 9, 5],
    [6, 7, 9, 8]
], dtype=float)

# LU decomposition
P, L, U = linalg.lu(A)

print("LU Decomposition: PA = LU")
print(f"\nL (lower triangular):\n{L.round(4)}")
print(f"\nU (upper triangular):\n{U.round(4)}")
print(f"\nVerification ||PA - LU||: {np.linalg.norm(P @ A - L @ U):.2e}")

# Solve for multiple right-hand sides efficiently
b1 = np.array([1, 2, 3, 4], dtype=float)
b2 = np.array([4, 3, 2, 1], dtype=float)

# Factor once
lu, piv = linalg.lu_factor(A)

# Solve for each RHS (cheap O(N^2) each)
x1 = linalg.lu_solve((lu, piv), b1)
x2 = linalg.lu_solve((lu, piv), b2)

print(f"\nSolution for b1: x = {x1.round(6)}")
print(f"Residual: {np.linalg.norm(A @ x1 - b1):.2e}")
print(f"\nSolution for b2: x = {x2.round(6)}")
print(f"Residual: {np.linalg.norm(A @ x2 - b2):.2e}")

print("\nKey insight: O(N^3) factorization done once,")
print("each solve is only O(N^2)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

3. Cholesky Decomposition

If \(A\) is symmetric positive definite (SPD), we can factor \(A = LL^T\) where \(L\) is lower-triangular with positive diagonal. This is twice as fast as LU and more numerically stable. SPD matrices arise naturally in least-squares, covariance matrices, and FEM stiffness matrices.

Cholesky Decomposition

Python

script.py29 lines

import numpy as np
from scipy import linalg

# Create a symmetric positive definite matrix
# A = B^T B guarantees SPD
B = np.random.RandomState(42).randn(4, 4)
A = B.T @ B + 0.1 * np.eye(4)  # add small diagonal for conditioning

print("Cholesky Decomposition: A = LL^T (for SPD matrices)")
print(f"A is symmetric: {np.allclose(A, A.T)}")
print(f"Eigenvalues (all positive = SPD): {np.linalg.eigvalsh(A).round(4)}")

# Cholesky factorization
L = linalg.cholesky(A, lower=True)
print(f"\nL (lower Cholesky factor):\n{L.round(4)}")
print(f"\nVerification ||A - LL^T||: {np.linalg.norm(A - L @ L.T):.2e}")

# Solve Ax = b using Cholesky
b = np.array([1, 2, 3, 4], dtype=float)
x = linalg.cho_solve(linalg.cho_factor(A), b)
print(f"\nSolution: x = {x.round(6)}")
print(f"Residual: {np.linalg.norm(A @ x - b):.2e}")

# Cost comparison
n = 1000
print(f"\nCost for {n}x{n} system:")
print(f"  Gaussian/LU: ~{2*n**3//3:,} flops")
print(f"  Cholesky:    ~{n**3//3:,} flops (2x faster!)")
print(f"  No pivoting needed for SPD matrices")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. Condition Number

How much can the solution change when the input data is slightly perturbed? The condition number quantifies this sensitivity:

\(\kappa(A) = \|A\| \cdot \|A^{-1}\| = \frac{\sigma_{\max}}{\sigma_{\min}}\)

Ratio of largest to smallest singular value. \(\kappa(A) \geq 1\) always.

Rule of thumb: You lose about \(\log_{10} \kappa(A)\)digits of accuracy when solving \(Ax = b\). With double precision (16 digits) and\(\kappa(A) = 10^{10}\), you get only ~6 reliable digits. If \(\kappa(A) \gtrsim 10^{16}\), the matrix is effectively singular in double precision.

Condition Number and Its Effects

Python

script.py39 lines

import numpy as np

def solve_and_analyze(A, b, name=""):
    """Solve Ax=b and analyze conditioning."""
    cond = np.linalg.cond(A)
    x = np.linalg.solve(A, b)
    residual = np.linalg.norm(A @ x - b)

# Perturb b slightly and see how much x changes
    db = np.random.RandomState(0).randn(len(b)) * 1e-10
    x_pert = np.linalg.solve(A, b + db)
    rel_change_x = np.linalg.norm(x_pert - x) / np.linalg.norm(x)
    rel_change_b = np.linalg.norm(db) / np.linalg.norm(b)

print(f"\n{'='*50}")
    print(f"{name}")
    print(f"{'='*50}")
    print(f"Condition number: {cond:.2e}")
    print(f"Digits of accuracy lost: ~{np.log10(cond):.1f}")
    print(f"Residual ||Ax-b||: {residual:.2e}")
    print(f"Relative perturbation in b: {rel_change_b:.2e}")
    print(f"Relative change in x: {rel_change_x:.2e}")
    print(f"Amplification factor: {rel_change_x/rel_change_b:.1f} (vs cond = {cond:.1f})")

# Well-conditioned system
A1 = np.array([[4, 1], [1, 3]], dtype=float)
b1 = np.array([1, 2], dtype=float)
solve_and_analyze(A1, b1, "Well-conditioned (kappa ~ 2)")

# Moderately ill-conditioned
A2 = np.array([[1, 1], [1, 1.0001]], dtype=float)
b2 = np.array([2, 2.0001], dtype=float)
solve_and_analyze(A2, b2, "Ill-conditioned (kappa ~ 40000)")

# Hilbert matrix: notoriously ill-conditioned
n = 8
A3 = np.array([[1.0/(i+j+1) for j in range(n)] for i in range(n)])
b3 = np.ones(n)
solve_and_analyze(A3, b3, f"Hilbert matrix {n}x{n} (extremely ill-conditioned)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Iterative Methods

For large sparse systems (\(N > 10^4\)), direct methods are too expensive. Iterative methods start from an initial guess and improve it, requiring only matrix-vector products (which are cheap for sparse matrices).

Jacobi Method

Split \(A = D + (L + U)\) where \(D\) is the diagonal. Iterate:

\(x^{(k+1)} = D^{-1}(b - (L+U)x^{(k)})\)

Converges if \(A\) is strictly diagonally dominant. Embarrassingly parallel.

Gauss-Seidel

Like Jacobi but uses updated values immediately:

\(x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j<i} a_{ij} x_j^{(k+1)} - \sum_{j>i} a_{ij} x_j^{(k)}\right)\)

Typically converges ~2x faster than Jacobi. Sequential.

Jacobi and Gauss-Seidel Iteration

Python

script.py59 lines

import numpy as np

def jacobi(A, b, x0=None, tol=1e-10, max_iter=1000):
    """Jacobi iterative method."""
    n = len(b)
    x = np.zeros(n) if x0 is None else x0.copy()
    D_inv = 1.0 / np.diag(A)
    R = A - np.diag(np.diag(A))  # off-diagonal
    history = []
    for k in range(max_iter):
        x_new = D_inv * (b - R @ x)
        res = np.linalg.norm(b - A @ x_new)
        history.append(res)
        if res < tol:
            return x_new, history
        x = x_new
    return x, history

def gauss_seidel(A, b, x0=None, tol=1e-10, max_iter=1000):
    """Gauss-Seidel iterative method."""
    n = len(b)
    x = np.zeros(n) if x0 is None else x0.copy()
    history = []
    for k in range(max_iter):
        for i in range(n):
            sigma = np.dot(A[i, :i], x[:i]) + np.dot(A[i, i+1:], x[i+1:])
            x[i] = (b[i] - sigma) / A[i, i]
        res = np.linalg.norm(b - A @ x)
        history.append(res)
        if res < tol:
            return x, history
    return x, history

# Create a diagonally dominant system (ensures convergence)
n = 6
np.random.seed(42)
A = np.random.randn(n, n)
A = A + A.T  # symmetric
A += n * np.eye(n)  # make diagonally dominant
b = np.random.randn(n)

x_exact = np.linalg.solve(A, b)

x_jac, hist_jac = jacobi(A, b)
x_gs, hist_gs = gauss_seidel(A, b)

print("Iterative Methods Comparison")
print(f"System size: {n}x{n}, diagonally dominant")
print(f"\n{'Iter':>4}  {'Jacobi ||r||':>14}  {'Gauss-Seidel ||r||':>18}")
print("-" * 42)
max_show = max(len(hist_jac), len(hist_gs))
for i in range(min(max_show, 20)):
    j = f"{hist_jac[i]:.6e}" if i < len(hist_jac) else "converged"
    g = f"{hist_gs[i]:.6e}" if i < len(hist_gs) else "converged"
    print(f"{i:4d}  {j:>14}  {g:>18}")

print(f"\nJacobi iterations: {len(hist_jac)}")
print(f"Gauss-Seidel iterations: {len(hist_gs)}")
print(f"Gauss-Seidel is typically ~2x faster")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Conjugate Gradient Method

The conjugate gradient (CG) method is the optimal Krylov subspace method for SPD matrices. It minimizes\(\|x - x^*\|_A = \sqrt{(x-x^*)^T A (x-x^*)}\) over successive Krylov subspaces. Converges in at most \(N\) iterations (exact arithmetic), but in practice converges much faster: \(O(\sqrt{\kappa(A)})\) iterations.

Conjugate Gradient for Sparse SPD Systems

Python

script.py34 lines

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

# Create a large sparse SPD system (1D Laplacian)
n = 200
# -1, 2, -1 tridiagonal matrix (arises from finite differences)
A = diags([-1, 2, -1], [-1, 0, 1], shape=(n, n), format='csr').toarray()
b = np.ones(n)

# Track convergence
residuals = []
def callback(xk):
    residuals.append(np.linalg.norm(b - A @ xk))

x_cg, info = cg(A, b, callback=callback, tol=1e-12)

print("Conjugate Gradient: 1D Laplacian (tridiagonal SPD)")
print(f"System size: {n}x{n}")
print(f"Condition number: {np.linalg.cond(A):.1f}")
print(f"CG iterations: {len(residuals)}")
print(f"Expected ~sqrt(kappa): {np.sqrt(np.linalg.cond(A)):.1f}")
print(f"Residual: {np.linalg.norm(b - A @ x_cg):.2e}")

# Show convergence
print(f"\n{'Iter':>4}  {'||residual||':>14}")
print("-" * 22)
for i in range(0, len(residuals), max(1, len(residuals)//15)):
    print(f"{i:4d}  {residuals[i]:14.6e}")
if len(residuals) > 0:
    print(f"{len(residuals):4d}  {residuals[-1]:14.6e}")

print(f"\nCG is ideal for large sparse SPD systems")
print(f"Cost per iteration: one matrix-vector multiply + O(N)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Practical Linear Algebra with NumPy

NumPy/SciPy Linear Algebra Toolkit

Python

script.py47 lines

import numpy as np
from scipy import linalg

np.random.seed(42)
n = 5
A = np.random.randn(n, n)
A = A + A.T + n * np.eye(n)  # SPD
b = np.random.randn(n)

print("=== NumPy/SciPy Linear Algebra Toolkit ===\n")

# 1. Basic solve
x = np.linalg.solve(A, b)
print(f"1. np.linalg.solve: residual = {np.linalg.norm(A@x - b):.2e}")

# 2. Determinant
print(f"2. det(A) = {np.linalg.det(A):.4f}")

# 3. Eigenvalues
evals = np.linalg.eigvalsh(A)  # eigvalsh for symmetric
print(f"3. Eigenvalues: {evals.round(3)}")

# 4. Condition number
print(f"4. Condition number: {np.linalg.cond(A):.4f}")

# 5. SVD
U, s, Vt = np.linalg.svd(A)
print(f"5. Singular values: {s.round(3)}")

# 6. Least squares (overdetermined system)
A_tall = np.random.randn(10, 3)
b_tall = np.random.randn(10)
x_ls, residuals, rank, sv = np.linalg.lstsq(A_tall, b_tall, rcond=None)
print(f"6. Least squares (10x3): x = {x_ls.round(4)}, rank = {rank}")

# 7. Matrix exponential (important for ODEs!)
from scipy.linalg import expm
B = np.array([[0, 1], [-1, 0]], dtype=float)  # rotation generator
expB = expm(B * np.pi/4)
print(f"7. exp(B*pi/4) for B=[0,1;-1,0]:\n{expB.round(4)}")
print(f"   (This is a 45-degree rotation matrix)")

# 8. Rank and null space
C = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float)
print(f"\n8. Rank of [[1,2,3],[4,5,6],[7,8,9]]: {np.linalg.matrix_rank(C)}")
print(f"   Singular values: {np.linalg.svd(C)[1].round(6)}")
print(f"   (Third singular value ~ 0: rank-deficient!)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server