12. Numerical and Computational Linear Algebra

Every earlier chapter assumed exact arithmetic and ignored cost. Computers have neither: they round every operation to about 16 digits, and a dense $100{,}000 \times 100{,}000$ system takes days to solve by elimination. This chapter studies how linear algebra is actually computed: floating point and pivoting, operation counts and sparse storage for direct solvers, iterative methods that solve huge sparse systems with only matrix-vector products, and the iterative algorithms behind every eigenvalue routine. These ideas explain what np.linalg.solve, scipy.sparse.linalg, and eigh do internally, and when they can be trusted.

<a id="121-direct-methods"></a>

12.1 Direct Methods

This subchapter studies methods that reach the answer in a fixed, finite number of steps (elimination and factorizations), and the three things that decide whether they work in practice: rounding error, pivoting, and cost in time and memory.

Floating Point Arithmetic and Roundoff

Why it matters. Section 1.2 warned that 0.1 + 0.2 != 0.3. In linear algebra, millions of such rounding errors accumulate, and they can be amplified. A solver can return a vector with no correct digits, even though every individual operation was accurate to 16 digits. To know when to trust a result, you need a model of rounding, and you need to separate two questions: is the problem sensitive to small changes in its data (conditioning), and does the algorithm introduce more error than necessary (stability)?

Analogy. Floating-point numbers are like a ruler whose tick marks get farther apart the farther you go from zero. Near 1 the ticks are $2.2 \times 10^{-16}$ apart; near $10^{16}$ they are 2 apart. Every measurement is snapped to the nearest tick, so the error is always small relative to the size of the number, never small in absolute terms.

The IEEE 754 double format. A double stores a number as

$$x = \pm,(1.b_1b_2\ldots b_{52})_2 \times 2^{e}, \qquad -1022 \le e \le 1023,$$

using 64 bits: 1 sign bit, 11 exponent bits, and 52 fraction bits. Key constants (np.finfo(float)):

Quantity	Value	Meaning
machine epsilon $\varepsilon_{\text{mach}}$	$2^{-52} \approx 2.2 \times 10^{-16}$	gap between 1 and the next double
unit roundoff $u$	$2^{-53} \approx 1.1 \times 10^{-16}$	maximum relative rounding error
largest double	$\approx 1.8 \times 10^{308}$	above this: overflow to inf
smallest normal double	$\approx 2.2 \times 10^{-308}$	below this: gradual underflow (subnormals), then 0
precision	about 16 significant decimal digits	$\log_{10}(1/u) \approx 15.95$

Single precision (float32, $u \approx 6 \times 10^{-8}$, about 7 digits) and half precision are used in graphics and machine learning when speed matters more than accuracy.

The standard model of rounding. IEEE 754 requires each basic operation to be computed exactly and then rounded. For $\circ \in {+, -, \times, \div}$ and $\sqrt{\cdot}$,

$$\operatorname{fl}(x \circ y) = (x \circ y)(1 + \delta), \qquad \lvert \delta \rvert \le u.$$

Every rounding-error analysis starts from this one line. Its consequences:

• Associativity and distributivity fail. $(a + b) + c$ and $a + (b + c)$ can differ, so the order of a summation changes its result, and parallel reductions are not bit-reproducible.
• Absorption. Adding a number smaller than half the local tick spacing changes nothing. Near $10^{16}$ the spacing is 2, so $10^{16} + 1$ is a tie, which rounds to the even neighbor $10^{16}$.
• Catastrophic cancellation. Subtracting two nearly equal numbers is exact, but it exposes the errors already present in them. If $a$ and $b$ agree in their first 10 digits, $a - b$ has at most 6 meaningful ones. The cure is algebraic: rewrite the formula to avoid the subtraction.
• Never test floats with ==. Compare with a tolerance relative to the size of the quantities involved.

Worked example: the quadratic formula. The roots of $x^2 + 10^8x + 1 = 0$ are approximately $-10^8$ and $-10^{-8}$. The textbook formula computes the small root as $(-b + \sqrt{b^2 - 4ac})/2a$. Here $\sqrt{b^2 - 4ac}$ agrees with $b$ in almost every digit, so the subtraction cancels and only a few correct digits survive. The algebraically equivalent form $2c/(-b - \sqrt{b^2 - 4ac})$ adds two numbers of the same sign and is accurate to full precision.

Conditioning versus stability. Two different things go wrong in numerical computation, and they must not be confused.

	Conditioning	Stability
Property of	the problem (e.g. $A\mathbf{x} = \mathbf{b}$ with this $A$)	the algorithm (e.g. elimination with partial pivoting)
Question	how much can the exact answer change when the data change slightly?	does the computed answer equal the exact answer for slightly changed data?
Measure	condition number $\kappa(A) = \lVert A \rVert,\lVert A^{-1} \rVert$ (Section 11.3)	backward error
Fixed by	reformulating the problem, more precision	choosing a better algorithm

An algorithm is backward stable if the computed $\hat{\mathbf{x}}$ is the exact solution of a nearby problem: $(A + \delta A)\hat{\mathbf{x}} = \mathbf{b}$ with $\lVert \delta A \rVert / \lVert A \rVert$ a modest multiple of $u$. The backward error is measured by the scaled residual $\lVert \mathbf{b} - A\hat{\mathbf{x}} \rVert / (\lVert A \rVert,\lVert \hat{\mathbf{x}} \rVert)$, which is computable. The forward error, the distance to the true solution, then obeys the rule of thumb

$$\frac{\lVert \hat{\mathbf{x}} - \mathbf{x} \rVert}{\lVert \mathbf{x} \rVert} \lesssim \kappa(A) \times \text{backward error}.$$

So you lose about $\log_{10}\kappa(A)$ of your 16 digits, even with a perfect algorithm. Householder QR and Cholesky are backward stable, and so is elimination with partial pivoting in practice (next topic). Classical Gram-Schmidt (Section 9.3) and forming the normal equations $A^TA$ for least squares (Section 9.5) are examples of avoidable instability.

import numpy as np
from scipy.linalg import hilbert

eps = np.finfo(float).eps
print(eps, eps / 2)                           # 2.220446049250313e-16 1.1102230246251565e-16: epsilon and u
print(1.0 + eps / 2 == 1.0, 1.0 + eps == 1.0) # True False: the gap between 1 and the next double is eps
print(np.spacing(1.0), np.spacing(1e16))      # 2.220446049250313e-16 2.0: spacing is relative, not absolute
print(1e16 + 1.0 - 1e16)                      # 0.0: the 1.0 was absorbed

# Floating-point addition is not associative
print((0.1 + 0.2) + 0.3 == 0.1 + (0.2 + 0.3)) # False

# Catastrophic cancellation: the small root of x^2 + 1e8 x + 1 = 0 is about -1e-8
a, b, c = 1.0, 1e8, 1.0
d = np.sqrt(b * b - 4 * a * c)
naive = (-b + d) / (2 * a)                    # subtracts two nearly equal numbers
stable = (2 * c) / (-b - d)                   # algebraically equal, no cancellation
print(naive, stable)                          # -7.450580596923828e-09 -1e-08: the naive value is 25% off

# Conditioning: a backward-stable solver still loses about log10(kappa) digits
for n in (4, 8, 12):
    H = hilbert(n)                            # H_ij = 1 / (i + j - 1): famously ill conditioned
    x_true = np.ones(n)
    b = H @ x_true
    x = np.linalg.solve(H, b)
    fwd = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)                  # forward error
    res = np.linalg.norm(b - H @ x) / (np.linalg.norm(H, 2) * np.linalg.norm(x))   # backward error
    print(n, f"{np.linalg.cond(H):.1e} {fwd:.1e} {res:.1e}")
# 4 1.6e+04 4.1e-14 0.0e+00
# 8 1.5e+10 6.1e-08 5.7e-17
# 12 1.8e+16 3.2e-01 9.8e-17                 tiny backward error, yet no correct digits: kappa ~ 1/u

Key takeaways

• Doubles carry about 16 significant digits; each operation satisfies $\operatorname{fl}(x \circ y) = (x \circ y)(1 + \delta)$ with $\lvert \delta \rvert \le u \approx 1.1 \times 10^{-16}$.
• Rounding error is relative, floating-point addition is not associative, and subtracting nearly equal numbers exposes earlier errors (catastrophic cancellation).
• Conditioning belongs to the problem, stability to the algorithm; a backward-stable algorithm loses about $\log_{10}\kappa$ digits, and no algorithm can do better in general.
• The scaled residual measures the backward error and is cheap to compute; the forward error can be much larger when $\kappa$ is large.

🧪 Practice

1. Compute $\sum_{k=1}^{10^6} 1/k^2$ in np.float32 summing forward ($k = 1, 2, \dots$) and backward ($k = 10^6, \dots, 1$). Which is closer to $\pi^2/6$, and why?
2. Rewrite $f(x) = \sqrt{x + 1} - \sqrt{x}$ and $g(x) = (1 - \cos x)/x^2$ to avoid cancellation for large $x$ and small $x$, respectively, and compare both versions numerically.
3. For the $10 \times 10$ Hilbert matrix, predict from $\kappa$ how many correct digits np.linalg.solve will deliver, then verify.
4. Interview: A colleague says, "The residual $\lVert \mathbf{b} - A\hat{\mathbf{x}} \rVert$ is $10^{-15}$, so the solution is accurate to 15 digits." When is this wrong, and what would you report instead? Hint: a small residual certifies a small backward error; the forward error can be $\kappa(A)$ times larger; estimate $\kappa$ (for example with np.linalg.cond or a LAPACK condition estimator).

Pivoting Strategies

Why it matters. Section 3.3 showed that a tiny pivot can destroy a solution, and Section 4.4 introduced partial pivoting, which keeps every multiplier at most 1 in absolute value. Small multipliers are not the whole story, though. What actually controls the rounding error of elimination is how large the entries become during elimination. This topic introduces the measure of that growth, the worst case, the alternatives to partial pivoting, and the matrices for which no pivoting is needed at all.

The growth factor. Let $a^{(k)}_{ij}$ be the entries after $k$ elimination steps. The growth factor is

$$\rho = \frac{\max_{i,j,k}\lvert a^{(k)}{ij} \rvert}{\max{i,j}\lvert a_{ij} \rvert}.$$

The backward error analysis of Gaussian elimination (Wilkinson) shows that the computed solution satisfies $(A + \delta A)\hat{\mathbf{x}} = \mathbf{b}$ with

$$\frac{\lVert \delta A \rVert_\infty}{\lVert A \rVert_\infty} \le c,n^3\rho,u$$

for a modest constant $c$, and the factor $n^3$ is very pessimistic in practice. Elimination is backward stable exactly when $\rho$ stays small. Without pivoting, $\rho$ can be unbounded, as the $10^{-20}$ example of Section 3.3 shows.

Partial pivoting and its worst case. With partial pivoting, each step at most doubles the largest entry, because $\lvert a_{ij} - l_{ik}a_{kj} \rvert \le 2\max\lvert a \rvert$ when $\lvert l_{ik} \rvert \le 1$. So $\rho \le 2^{n-1}$, and the bound is attained by Wilkinson's matrix:

$$W_n = \begin{bmatrix} 1 & & & & 1 \ -1 & 1 & & & 1 \ -1 & -1 & 1 & & 1 \ \vdots & & & \ddots & \vdots \ -1 & -1 & \cdots & -1 & 1 \end{bmatrix}.$$

No row exchanges occur (all candidates have absolute value 1), and each step doubles the last column: $1, 2, 4, \dots, 2^{n-1}$. For $n = 60$, $\rho = 2^{59} \approx 5.8 \times 10^{17}$, and the computed solution has no correct digits, even though $\kappa(W_{60}) \approx 27$ and the problem itself is perfectly well conditioned.

Why partial pivoting is used anyway. Matrices like $W_n$ essentially never occur in applications. Experiments with random matrices give growth factors around $n^{1/2}$ to $n^{2/3}$, and decades of practice show the same for real problems. This empirical fact, which is still not fully explained theoretically, is why LAPACK's default solver (getrf) uses partial pivoting: it is cheap ($O(n^2)$ comparisons in total) and stable in practice.

The pivoting menu.

Strategy	Pivot chosen at step $k$	Extra cost	Growth bound	Use
none	$a_{kk}$ as it stands	0	unbounded	SPD and diagonally dominant matrices only
partial	largest $\lvert a_{ik} \rvert$ in column $k$ (row swaps)	$O(n^2)$ total	$2^{n-1}$, small in practice	the default (np.linalg.solve, scipy.linalg.lu)
scaled partial	largest $\lvert a_{ik} \rvert / s_i$, where $s_i$ is the largest entry of row $i$	$O(n^2)$	as partial	rows with very different scales
rook	an entry that is largest in both its row and its column	usually $O(n^2)$	much smaller than $2^{n-1}$	when extra safety is wanted cheaply
complete	largest $\lvert a_{ij} \rvert$ in the whole remaining block (row and column swaps)	$O(n^3)$ comparisons	about $n^{1/2 + \frac14\ln n}$ (slowly growing)	rank-revealing factorizations, research
symmetric (Bunch-Kaufman)	$1 \times 1$ or $2 \times 2$ diagonal pivots, preserving symmetry	$O(n^2)$	bounded	symmetric indefinite $LDL^T$ (scipy.linalg.ldl)
threshold	any $\lvert a_{ik} \rvert \ge \tau\max_i\lvert a_{ik} \rvert$, e.g. $\tau = 0.1$	$O(n^2)$	$(1 + 1/\tau)^{n-1}$	sparse LU: trades stability for less fill-in (next topics)

When no pivoting is needed.

• Symmetric positive definite matrices. Cholesky needs no pivoting, and its growth factor is 1, since every $\lvert l_{ij} \rvert \le \sqrt{a_{ii}}$ (Section 10.3).
• Strictly diagonally dominant matrices (by columns): the pivot is automatically the largest entry in its column, so partial pivoting would never swap.
• Tridiagonal systems that are diagonally dominant, solved by the $O(n)$ Thomas algorithm.

import numpy as np

def lu_growth(A, pivot="partial"):
    """Gaussian elimination; returns the growth factor max|entries during elimination| / max|A|."""
    U = A.astype(float).copy()
    n = len(U)
    biggest = np.abs(U).max()
    for k in range(n - 1):
        if pivot == "partial":                       # largest entry in column k, rows k..n-1
            p = k + np.argmax(np.abs(U[k:, k]))
            U[[k, p]] = U[[p, k]]
        elif pivot == "complete":                    # largest entry in the whole trailing block
            i, j = np.unravel_index(np.argmax(np.abs(U[k:, k:])), U[k:, k:].shape)
            U[[k, k + i]] = U[[k + i, k]]            # row swap
            U[:, [k, k + j]] = U[:, [k + j, k]]      # column swap (reorders the unknowns)
        U[k+1:, k:] -= np.outer(U[k+1:, k] / U[k, k], U[k, k:])
        biggest = max(biggest, np.abs(U).max())
    return biggest / np.abs(A).max()

def wilkinson(n):
    """1 on the diagonal, -1 below it, 1 in the last column: the worst case for partial pivoting."""
    W = np.eye(n) - np.tril(np.ones((n, n)), -1)
    W[:, -1] = 1
    return W

print(wilkinson(4))
# [[ 1.  0.  0.  1.]
#  [-1.  1.  0.  1.]
#  [-1. -1.  1.  1.]
#  [-1. -1. -1.  1.]]
for n in (10, 30, 60):
    W = wilkinson(n)
    print(n, f"{lu_growth(W, 'partial'):.3g}", lu_growth(W, "complete"))
# 10 512 2.0                                 2^9
# 30 5.37e+08 2.0                            2^29
# 60 5.76e+17 2.0                            2^59: partial pivoting fails, complete pivoting does not

# The damage is real: LAPACK (partial pivoting) on W with n = 60
n = 60
W = wilkinson(n)
x_true = np.random.default_rng(0).standard_normal(n)
x = np.linalg.solve(W, W @ x_true)
print(f"{np.linalg.cond(W):.1f}", f"{np.linalg.norm(x - x_true) / np.linalg.norm(x_true):.1e}")
# 26.8 2.5e+00                               well-conditioned problem, 250% error

# Random matrices: growth under partial pivoting stays small
rng = np.random.default_rng(1)
print([round(float(lu_growth(rng.standard_normal((n, n)))), 1) for n in (50, 100, 200)])   # [3.0, 7.8, 6.9]

Key takeaways

• The growth factor $\rho$, the largest entry produced during elimination relative to the largest entry of $A$, controls the backward error of Gaussian elimination.
• Partial pivoting bounds the multipliers by 1 and $\rho$ by $2^{n-1}$; the bound is attained by Wilkinson's matrix, but in practice $\rho$ is small.
• Complete and rook pivoting have far better worst cases at a higher cost; Bunch-Kaufman pivoting preserves symmetry for indefinite matrices; threshold pivoting trades some stability for sparsity.
• SPD and diagonally dominant matrices need no pivoting at all.

🧪 Practice

1. Carry out elimination with partial pivoting on $W_4$ by hand and verify that the last column becomes $(1, 2, 4, 8)$.
2. Show that if $\lvert l_{ik} \rvert \le 1$ for all multipliers, then each elimination step at most doubles $\max_{i,j}\lvert a_{ij} \rvert$, and conclude that $\rho \le 2^{n-1}$.
3. Explain why a strictly column diagonally dominant matrix never triggers a row exchange under partial pivoting. (Show that diagonal dominance is preserved by one elimination step.)
4. Interview: If complete pivoting has a much better worst case, why does every standard library use partial pivoting by default? Hint: compare the cost of searching the whole trailing block at every step ($O(n^3)$ comparisons, poor memory access) with the observed behavior of partial pivoting on real and random matrices; mention when you would switch (rank-revealing needs, or a detected large growth factor).

Operation Counts and Complexity

Why it matters. Before running an algorithm on a large matrix, you should be able to predict whether it will take a millisecond or a week, and whether the matrix even fits in memory. Counting floating-point operations (flops) gives that prediction, explains rules such as "never form $A^{-1}$" and "factor once, solve many times", and shows why $10^5$ unknowns are routine for sparse problems but out of reach for dense elimination.

Counting flops. One flop is one floating-point addition, subtraction, multiplication, or division. Counts are reported by their leading term, because lower-order terms are negligible for large $n$. For example, elimination performs, at step $k$, about $(n - k)^2$ multiply-subtract pairs on the trailing block, so the total is

$$\sum_{k=1}^{n-1} 2(n - k)^2 \approx 2\int_0^n t^2,dt = \tfrac23n^3.$$

Reference table ($n \times n$ unless stated):

Operation	Flops (leading term)	Notes
dot product, axpy ($\mathbf{y} \leftarrow a\mathbf{x} + \mathbf{y}$)	$2n$	BLAS level 1
matrix-vector product ($m \times n$)	$2mn$	BLAS level 2
triangular solve	$n^2$	forward or back substitution
matrix-matrix product	$2n^3$	BLAS level 3; Strassen-type methods: $O(n^{2.81})$
LU with partial pivoting	$\tfrac23n^3$	then $2n^2$ per right-hand side
Cholesky	$\tfrac13n^3$	SPD only (Section 10.3)
Householder QR ($m \times n$)	$2mn^2 - \tfrac23n^3$	$\tfrac43n^3$ when square (Section 9.3)
inverse $A^{-1}$	$2n^3$	three times LU, and less accurate
determinant	$\tfrac23n^3$	product of the pivots of LU
symmetric eigenvalues (values / with vectors)	$\approx \tfrac43n^3$ / $\approx 9n^3$	tridiagonalization plus QR (Section 12.3)
general eigenvalues (values / with vectors)	$\approx 10n^3$ / $\approx 25n^3$	Hessenberg reduction plus QR
banded LU, bandwidth $p$	$\approx 2np^2$	tridiagonal ($p = 1$): $\approx 8n$
sparse matrix-vector product	$2,\operatorname{nnz}(A)$	nnz = number of stored nonzeros

Consequences.

• Factor once, solve many. Solving with $k$ right-hand sides costs $\tfrac23n^3 + 2kn^2$ with one LU, versus $k \cdot \tfrac23n^3$ if each solve refactors. In SciPy: lu_factor once, then lu_solve; in NumPy: pass all right-hand sides as columns of one matrix.
• Never form the inverse to solve a system. inv(A) @ b costs three times as much as solve(A, b) and has a larger backward error. For an ill-conditioned matrix the residual can be $10^{10}$ times worse, as the example shows. If you need $A^{-1}B$, solve $AX = B$.
• Exploit structure. Using Cholesky instead of LU halves the work, and so does eigh instead of eig (even better). Bandedness turns $n^3$ into $np^2$.
• Cubic growth is unforgiving. Doubling $n$ multiplies the time by 8. At $10^{10}$ flops per second, dense LU takes 0.07 s for $n = 1000$, about 70 s for $n = 10^4$, and about 19 hours for $n = 10^5$, which would also need 80 GB just to store $A$.

Flops are not the whole story: memory traffic. Modern processors can do arithmetic far faster than they can fetch data from main memory. What matters is arithmetic intensity, the number of flops per number moved:

BLAS level	Example	Data	Flops	Flops per datum
1	axpy	$3n$	$2n$	$O(1)$: memory-bound
2	matrix-vector	$n^2$	$2n^2$	$O(1)$: memory-bound
3	matrix-matrix	$3n^2$	$2n^3$	$O(n)$: can run at peak speed

LAPACK reorganizes LU, Cholesky, and QR into blocked algorithms whose inner work is matrix-matrix products. That is why a dense factorization on a laptop runs at tens of billions of flops per second, while a naive triple loop in Python runs about a million times slower. Sparse iterative methods, built on matrix-vector products, are memory-bound, so their speed is measured in bytes per second rather than flops.

import time
import numpy as np
from scipy.linalg import lu_factor, lu_solve

rng = np.random.default_rng(0)

def seconds(f, repeat=3):
    best = np.inf
    for _ in range(repeat):
        t0 = time.perf_counter(); f(); best = min(best, time.perf_counter() - t0)
    return best

# Cubic scaling: doubling n multiplies the cost of LU by about 8 (timings are machine-dependent)
for n in (1000, 2000):
    A = rng.standard_normal((n, n))
    t = seconds(lambda: lu_factor(A))
    print(n, f"{t:.3f} s", f"{(2 / 3) * n**3 / t / 1e9:.0f} GFLOP/s")
# 1000 0.013 s 53 GFLOP/s                    blocked LAPACK code runs near the processor's peak
# 2000 0.093 s 57 GFLOP/s

# Factor once, solve many: the O(n^3) step is paid once, each new right-hand side costs O(n^2)
n = 2000
A = rng.standard_normal((n, n))
B = rng.standard_normal((n, 100))
t_refactor = seconds(lambda: [np.linalg.solve(A, B[:, j]) for j in range(100)], repeat=1)
lu_piv = lu_factor(A)
t_reuse = seconds(lambda: lu_solve(lu_piv, B), repeat=1)
print(f"{t_refactor / t_reuse:.0f}x faster when the factorization is reused")   # several hundred times

# solve beats inv on an ill-conditioned matrix (condition number 1e12): smaller residual, smaller error
n = 500
U, _ = np.linalg.qr(rng.standard_normal((n, n)))
V, _ = np.linalg.qr(rng.standard_normal((n, n)))
A = U @ np.diag(np.logspace(0, -12, n)) @ V.T       # singular values from 1 down to 1e-12
x_true = rng.standard_normal(n)
b = A @ x_true
for name, x in [("solve", np.linalg.solve(A, b)), ("inv", np.linalg.inv(A) @ b)]:
    print(name, f"{np.linalg.norm(A @ x - b) / np.linalg.norm(b):.1e}",
          f"{np.linalg.norm(x - x_true) / np.linalg.norm(x_true):.1e}")
# solve 1.4e-15 1.2e-05                      relative residual, relative error
# inv 2.7e-05 2.7e-04                        residual 10^10 times worse, error 20 times worse

Key takeaways

• Count flops by the leading term: matrix-vector $2n^2$, LU $\tfrac23n^3$, Cholesky $\tfrac13n^3$, QR $\tfrac43n^3$, inverse $2n^3$, matrix product $2n^3$.
• Factor once and reuse the factorization for every right-hand side; never form $A^{-1}$ to solve a system.
• Structure (symmetry, definiteness, bandedness, sparsity) reduces cost by large factors.
• Speed depends on memory traffic as well as flops: matrix-matrix work (BLAS 3) runs near peak, matrix-vector work (BLAS 2) is memory-bound.

🧪 Practice

1. Derive the leading term $n^3/3$ for Cholesky by counting the work in the column-by-column formulas of Section 10.3.
2. A simulation solves $A\mathbf{x}_t = \mathbf{b}_t$ for 10,000 time steps with the same $2000 \times 2000$ matrix $A$. Estimate the total flops with and without reusing the LU factorization.
3. Show that computing $A^{-1}$ via LU costs about $2n^3$ flops, and explain why exploiting the zeros in the identity's columns matters for reaching that count.
4. Interview: Your service computes np.linalg.inv(X.T @ X) @ X.T @ y for a $10^6 \times 500$ design matrix, and it is slow and occasionally inaccurate. Estimate its cost and propose a faster, more accurate replacement. Hint: forming $X^TX$ costs $2mn^2 \approx 5 \times 10^{11}$ flops and squares the condition number; compare Cholesky of $X^TX$ with solve instead of inv, and QR or lstsq on $X$ directly (Section 9.5) when accuracy matters.

Sparse Matrix Storage and Solvers

Why it matters. The matrices of discretized PDEs, circuits, social networks, and recommender systems have millions of rows but only a handful of nonzeros per row. Storing them densely is impossible: a $10^6 \times 10^6$ double matrix needs 8 TB. Storing only the nonzeros makes memory and matrix-vector products proportional to the number of nonzeros, $\operatorname{nnz}(A)$. Direct solvers can also exploit sparsity, but they must fight fill-in: elimination creates new nonzeros, and the order of the unknowns decides how many.

Storage formats. For the matrix

$$A = \begin{bmatrix} 4 & 0 & 0 & 1 \ 0 & 3 & 0 & 0 \ 2 & 0 & 5 & 0 \ 0 & 0 & 6 & 7 \end{bmatrix}:$$

COO (coordinate): one (row, col, value) triple per nonzero; easy to build, then convert
    row   = [0, 0, 1, 2, 2, 3, 3]
    col   = [0, 3, 1, 0, 2, 2, 3]
    value = [4, 1, 3, 2, 5, 6, 7]

CSR (compressed sparse row): rows stored one after another
    data    = [4, 1 | 3 | 2, 5 | 6, 7]      values, row by row
    indices = [0, 3 | 1 | 0, 2 | 2, 3]      column of each value
    indptr  = [0, 2, 3, 5, 7]               row i is data[indptr[i] : indptr[i+1]]

CSC (compressed sparse column): the same idea by columns (the transpose's CSR)

Format	Strength	Typical use
COO	simple construction, duplicates are summed on conversion	assembling finite element matrices
CSR	fast row access and matrix-vector products	iterative solvers, A @ x
CSC	fast column access	sparse direct solvers (splu, spsolve)
banded / DIA	stores a few diagonals	tridiagonal and banded systems
LIL / DOK	cheap insertion of single entries	incremental construction (then convert)

CSR needs $2\operatorname{nnz} + n + 1$ numbers, and $A\mathbf{x}$ costs $2\operatorname{nnz}$ flops.

Fill-in. Eliminating a variable connects all of its neighbors. In graph terms, the vertices adjacent to the eliminated vertex become a clique. New nonzeros appear where $A$ had zeros, and they cost memory and time. The ordering of the unknowns changes the amount of fill dramatically. The classic example is the arrow matrix:

  dense row/column FIRST                         dense row/column LAST
  A                L (Cholesky)                  A                L (Cholesky)
  * * * * * *      *                             * . . . . *      *
  * * . . . .      * *                           . * . . . *      . *
  * . * . . .      * * *                         . . * . . *      . . *
  * . . * . .      * * * *                       . . . * . *      . . . *
  * . . . * .      * * * * *                     . . . . * *      . . . . *
  * . . . . *      * * * * * *                   * * * * * *      * * * * * *
  eliminating x1 couples everything:             eliminating x1..x5 touches only the last row:
  L is completely full                           no fill at all

Finding the ordering with the least fill is NP-hard, but good heuristics work well:

Ordering	Idea	Good for
reverse Cuthill-McKee (RCM)	breadth-first search to shrink the bandwidth	banded and profile solvers
(approximate) minimum degree, AMD / COLAMD	eliminate the vertex with fewest neighbors first	general sparse matrices; default in many solvers
nested dissection	split the graph with a small separator, order the separator last, recurse	2D/3D meshes: fill $O(n\log n)$ in 2D instead of $O(n^{3/2})$ for banded

Sparse direct solvers. Libraries such as SuperLU (scipy.sparse.linalg.splu, spsolve), CHOLMOD, UMFPACK, MUMPS, and PARDISO work in four phases: (1) choose a fill-reducing ordering; (2) symbolic factorization, which predicts the nonzero structure of the factors; (3) numeric factorization, using threshold pivoting so that stability does not destroy the chosen ordering; (4) triangular solves. They are robust and give accurate answers. Their weakness is memory: for 3D meshes, fill grows quickly, and iterative methods (Section 12.2) take over somewhere between $10^5$ and $10^7$ unknowns.

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve, splu

A = np.array([[4., 0., 0., 1.],
              [0., 3., 0., 0.],
              [2., 0., 5., 0.],
              [0., 0., 6., 7.]])
S = sp.csr_matrix(A)
print(S.data)       # [4. 1. 3. 2. 5. 6. 7.]: nonzero values, row by row
print(S.indices)    # [0 3 1 0 2 2 3]: column index of each value
print(S.indptr)     # [0 2 3 5 7]: row i occupies data[indptr[i]:indptr[i+1]]

def csr_matvec(data, indices, indptr, x):
    y = np.zeros(len(indptr) - 1)
    for i in range(len(y)):
        for k in range(indptr[i], indptr[i + 1]):   # only the stored entries of row i
            y[i] += data[k] * x[indices[k]]
    return y
x = np.array([1., 2., 3., 4.])
print(csr_matvec(S.data, S.indices, S.indptr, x), A @ x)   # [ 8.  6. 17. 46.] [ 8.  6. 17. 46.]

# 2D Poisson equation on an m x m grid: n = m^2 unknowns, about 5 nonzeros per row
m = 300
T = sp.diags([-1., 2., -1.], [-1, 0, 1], shape=(m, m))
P = (sp.kron(sp.identity(m), T) + sp.kron(T, sp.identity(m))).tocsc()
n = P.shape[0]
print(n, P.nnz, f"{P.nnz * 12 / 1e6:.1f} MB sparse vs {n * n * 8 / 1e9:.0f} GB dense")
# 90000 448800 5.4 MB sparse vs 65 GB dense
b = np.ones(n)
u = spsolve(P, b)
print(f"{np.linalg.norm(P @ u - b) / np.linalg.norm(b):.1e}")   # 2.1e-12

# Fill-in depends on the ordering of the unknowns
for order in ("NATURAL", "COLAMD"):
    lu = splu(P, permc_spec=order)
    print(order, lu.L.nnz + lu.U.nnz)
# NATURAL 54000598                           natural (row-by-row grid) ordering
# COLAMD 8902568                             fill-reducing ordering: 6 times fewer nonzeros

# Arrow matrices: dense first row/column fills everything; dense last row/column fills nothing
k = 8
arrow = 4 * np.eye(k); arrow[0, :] = arrow[:, 0] = 1; arrow[0, 0] = k
for name, M in (("arrow first", arrow), ("arrow last", arrow[::-1, ::-1])):
    L = np.linalg.cholesky(M)
    print(name, int((np.abs(L) > 1e-14).sum()), "nonzeros in L, out of", k * (k + 1) // 2)
# arrow first 36 nonzeros in L, out of 36
# arrow last 15 nonzeros in L, out of 36

Key takeaways

• Sparse formats store only nonzeros: COO for assembly, CSR for row access and $A\mathbf{x}$, CSC for direct solvers; memory and $A\mathbf{x}$ cost $O(\operatorname{nnz})$.
• Elimination creates fill-in; in graph terms, eliminating a vertex turns its neighbors into a clique.
• The ordering of the unknowns controls fill; AMD/COLAMD, RCM, and nested dissection are the standard heuristics.
• Sparse direct solvers run ordering, symbolic factorization, numeric factorization with threshold pivoting, and triangular solves; they are robust but memory-hungry for 3D problems.

🧪 Practice

1. Write the CSR arrays (data, indices, indptr) of $\begin{bmatrix} 0 & 2 & 0 \ 1 & 0 & 3 \ 0 & 0 & 4 \end{bmatrix}$ by hand, and compute $A\mathbf{x}$ for $\mathbf{x} = (1, 1, 1)$ using only those arrays.
2. Show that the Cholesky factor of a tridiagonal SPD matrix is bidiagonal, so there is no fill-in, and deduce an $O(n)$ algorithm.
3. Draw the graph of the $5 \times 5$ arrow matrix with the dense row first. Eliminate vertex 1 and show that the remaining vertices form a clique. Then explain why ordering the hub last avoids this.
4. Interview: A finite element model with $2 \times 10^6$ unknowns in 3D runs out of memory in spsolve, although the matrix itself needs only 1 GB. Explain why, and list your options. Hint: 3D fill-in grows much faster than nnz even with nested dissection; options are a better ordering or an out-of-core solver, or switching to a preconditioned iterative method (CG or GMRES with multigrid or incomplete factorization, Section 12.2).

<a id="122-iterative-methods-for-linear-systems"></a>

12.2 Iterative Methods for Linear Systems

This subchapter solves $A\mathbf{x} = \mathbf{b}$ by improving an approximate solution step by step. It touches $A$ only through cheap operations, matrix-vector products or sweeps over its rows, so the work is proportional to $\operatorname{nnz}(A)$ and no fill-in is created.

Jacobi Method

Why it matters. Direct methods pay $O(n^3)$ for dense matrices, or suffer fill-in for sparse ones, before producing any answer. Iterative methods start from a guess and refine it. Each step costs about one sweep over the nonzeros, memory stays at a few vectors, and you can stop as soon as the answer is good enough. Jacobi is the simplest such method. It is rarely the best solver on its own, but its structure explains every method that follows, and it survives in practice as a parallel smoother inside multigrid.

Analogy. Imagine $n$ people, each responsible for one unknown. Each equation says how person $i$'s value depends on the others. In a Jacobi round, everyone looks at the values announced in the previous round, solves their own equation for their own value, and then everyone announces simultaneously. Repeat until nobody changes much.

The method. Solve equation $i$ for $x_i$, using the old values of all other unknowns:

$$x_i^{(k+1)} = \frac{1}{a_{ii}}\Big(b_i - \sum_{j \ne i}a_{ij}x_j^{(k)}\Big), \qquad i = 1, \dots, n.$$

Matrix form: splitting. Write $A = D + L + U$, with $D$ diagonal, $L$ strictly lower triangular, and $U$ strictly upper triangular. These are pieces of $A$, not the LU factors. Then $A\mathbf{x} = \mathbf{b}$ becomes $D\mathbf{x} = \mathbf{b} - (L + U)\mathbf{x}$, and Jacobi iterates

$$\mathbf{x}^{(k+1)} = D^{-1}\big(\mathbf{b} - (L + U)\mathbf{x}^{(k)}\big) = G_J\mathbf{x}^{(k)} + D^{-1}\mathbf{b}, \qquad G_J = -D^{-1}(L + U).$$

Every iterative method of this subchapter has the form $\mathbf{x}^{(k+1)} = G\mathbf{x}^{(k)} + \mathbf{c}$. The iteration matrix $G$ decides convergence (see Convergence Criteria).

When it converges. If $A$ is strictly diagonally dominant, $\lvert a_{ii} \rvert > \sum_{j \ne i}\lvert a_{ij} \rvert$ for every row, then $\lVert G_J \rVert_\infty = \max_i\sum_{j \ne i}\lvert a_{ij} \rvert / \lvert a_{ii} \rvert < 1$. So the error shrinks by at least that factor in every sweep, for any starting vector. Without a dominant diagonal, Jacobi can diverge.

Worked example. Solve

$$\begin{aligned} 10x_1 - x_2 + 2x_3 &= 6 \ -x_1 + 11x_2 - x_3 &= 22 \ 2x_1 - x_2 + 10x_3 &= -10 \end{aligned} \qquad (\text{exact solution } (1, 2, -1)).$$

Starting from $\mathbf{x}^{(0)} = \mathbf{0}$:

$k$	$x_1 = (6 + x_2 - 2x_3)/10$	$x_2 = (22 + x_1 + x_3)/11$	$x_3 = (-10 - 2x_1 + x_2)/10$
1	$0.6$	$2.0$	$-1.0$
2	$1.0$	$1.9636$	$-0.92$
3	$0.9804$	$2.0073$	$-1.0036$

Each row of the table uses only the row above it.

Properties.

Property	Jacobi
Work per sweep	one pass over the nonzeros, $O(\operatorname{nnz})$
Memory	two vectors (old and new)
Parallelism	perfect: all $n$ updates are independent
Convergence guaranteed for	strictly diagonally dominant matrices; SPD $A$ only if $2D - A$ is also SPD
Speed	slow; for the Poisson matrix, $O(n^2)$ sweeps in 1D
Variant	weighted (damped) Jacobi $\mathbf{x} \leftarrow \mathbf{x} + \omega D^{-1}(\mathbf{b} - A\mathbf{x})$, a standard multigrid smoother

import numpy as np

A = np.array([[10., -1., 2.],
              [-1., 11., -1.],
              [2., -1., 10.]])                # strictly diagonally dominant
b = np.array([6., 22., -10.])                 # chosen so that x = (1, 2, -1)
x_exact = np.linalg.solve(A, b)

def jacobi(A, b, tol=1e-10, max_iter=500):
    D = np.diag(A)
    R = A - np.diag(D)                        # off-diagonal part L + U
    x = np.zeros_like(b)
    for k in range(1, max_iter + 1):
        x = (b - R @ x) / D                   # every component uses only the OLD x: parallel
        if np.linalg.norm(b - A @ x) <= tol * np.linalg.norm(b):
            return x, k
    return x, max_iter

x, iters = jacobi(A, b)
print(np.round(x, 10), iters, np.allclose(x, x_exact))   # [ 1.  2. -1.] 17 True

# The first three sweeps, as in the table
x = np.zeros(3)
for k in range(3):
    x = (b - (A - np.diag(np.diag(A))) @ x) / np.diag(A)
    print(k + 1, np.round(x, 4))
# 1 [ 0.6  2.  -1. ]
# 2 [ 1.      1.9636 -0.92  ]
# 3 [ 0.9804  2.0073 -1.0036]

# Jacobi can diverge when the diagonal does not dominate
B = np.array([[1., 2.], [3., 1.]])
G = -np.diag(1 / np.diag(B)) @ (B - np.diag(np.diag(B)))   # Jacobi iteration matrix
print(round(np.max(np.abs(np.linalg.eigvals(G))), 4),      # 2.4495: spectral radius sqrt(6) > 1
      f"{np.linalg.norm(jacobi(B, np.array([3., 4.]), max_iter=50)[0]):.1e}")   # 4.0e+19: diverged

Key takeaways

• Jacobi solves each equation for its own unknown using the previous iterate: $\mathbf{x}^{(k+1)} = D^{-1}(\mathbf{b} - (L + U)\mathbf{x}^{(k)})$.
• It is a splitting method with iteration matrix $G_J = -D^{-1}(L + U)$; strict diagonal dominance guarantees convergence.
• All updates are independent, so Jacobi parallelizes perfectly, but it converges slowly.
• Damped Jacobi is an effective smoother inside multigrid, which is where Jacobi lives on in practice.

🧪 Practice

1. Perform two Jacobi sweeps by hand for $4x_1 - x_2 = 3$, $-x_1 + 4x_2 = 3$ starting from $\mathbf{0}$, and compare with the exact solution $(1, 1)$.
2. Show that $\lVert G_J \rVert_\infty < 1$ for every strictly row diagonally dominant matrix, and conclude that the error satisfies $\lVert \mathbf{e}^{(k)} \rVert_\infty \le \lVert G_J \rVert_\infty^k\lVert \mathbf{e}^{(0)} \rVert_\infty$.
3. Reorder the equations of the divergent example $x_1 + 2x_2 = 3$, $3x_1 + x_2 = 4$ so that Jacobi converges, and explain why reordering helps.
4. Interview: Why would anyone use Jacobi on a GPU, when Gauss-Seidel converges about twice as fast on a CPU? Hint: compare data dependencies within a sweep; Jacobi's $n$ updates are independent and map onto thousands of threads, while Gauss-Seidel's are sequential unless the unknowns are recolored.

Gauss-Seidel Method

Why it matters. Jacobi ignores information it already has: when it updates $x_2$, a new value of $x_1$ is already available, but Jacobi uses the old one. Gauss-Seidel uses each new value immediately. This costs nothing extra, roughly halves the number of sweeps for typical problems, and needs only one vector of storage. For symmetric positive definite matrices it has a clean interpretation as exact minimization along one coordinate at a time.

Analogy. In the room of $n$ people, Gauss-Seidel lets people speak in turn. Each person updates using the latest values already announced in this round, so good news spreads within a sweep instead of waiting for the next one.

The method.

$$x_i^{(k+1)} = \frac{1}{a_{ii}}\Big(b_i - \sum_{j < i}a_{ij}x_j^{(k+1)} - \sum_{j > i}a_{ij}x_j^{(k)}\Big), \qquad i = 1, \dots, n \text{ in order}.$$

In matrix form, $(D + L)\mathbf{x}^{(k+1)} = \mathbf{b} - U\mathbf{x}^{(k)}$, so

$$G_{GS} = -(D + L)^{-1}U.$$

Applying $(D + L)^{-1}$ is just forward substitution, which the sweep performs implicitly.

Worked example. For the system of the Jacobi topic, starting from $\mathbf{0}$:

• $x_1 = (6 + 0 - 0)/10 = 0.6$;
• $x_2 = (22 + \mathbf{0.6} + 0)/11 = 2.0545$, using the new $x_1$;
• $x_3 = (-10 - 2(\mathbf{0.6}) + \mathbf{2.0545})/10 = -0.9145$, using the new $x_1$ and $x_2$.

Gauss-Seidel reaches a relative residual of $10^{-10}$ in 8 sweeps, against 17 for Jacobi.

Gauss-Seidel as coordinate descent. For SPD $A$, solving $A\mathbf{x} = \mathbf{b}$ is equivalent to minimizing $\varphi(\mathbf{x}) = \tfrac12\mathbf{x}^TA\mathbf{x} - \mathbf{b}^T\mathbf{x}$ (Section 10.3). Setting $\partial\varphi/\partial x_i = (A\mathbf{x})_i - b_i = 0$ and solving for $x_i$ gives exactly the Gauss-Seidel update. So each update minimizes $\varphi$ exactly along one coordinate axis, $\varphi$ never increases, and the iteration converges for every SPD matrix.

Convergence.

Matrix class	Gauss-Seidel
strictly diagonally dominant	converges
symmetric positive definite	always converges (coordinate descent)
"consistently ordered" (e.g. tridiagonal, Poisson with natural ordering)	$\rho(G_{GS}) = \rho(G_J)^2$: one GS sweep is worth two Jacobi sweeps

The price: sequential updates. Within a sweep, $x_i$ depends on $x_1, \dots, x_{i-1}$, so the sweep cannot be parallelized directly. The standard fix is multicoloring. For the grid Poisson matrix, color the unknowns like a checkerboard (red-black ordering). Red unknowns depend only on black ones and vice versa, so each color is updated in parallel. Symmetric Gauss-Seidel (a forward sweep followed by a backward sweep) gives a symmetric iteration, which matters when it is used as a preconditioner for CG.

import numpy as np

A = np.array([[10., -1., 2.],
              [-1., 11., -1.],
              [2., -1., 10.]])
b = np.array([6., 22., -10.])

def gauss_seidel(A, b, tol=1e-10, max_iter=500):
    n = len(b)
    x = np.zeros_like(b)
    for k in range(1, max_iter + 1):
        for i in range(n):                    # x[:i] already holds NEW values from this sweep
            x[i] = (b[i] - A[i, :i] @ x[:i] - A[i, i+1:] @ x[i+1:]) / A[i, i]
        if np.linalg.norm(b - A @ x) <= tol * np.linalg.norm(b):
            return x, k
    return x, max_iter

x, iters = gauss_seidel(A, b)
print(np.round(x, 10), iters)                 # [ 1.  2. -1.] 8   (Jacobi needed 17)

# 1D Poisson matrix: Gauss-Seidel needs about half as many sweeps as Jacobi
def jacobi(A, b, tol=1e-8, max_iter=100_000):
    D = np.diag(A); R = A - np.diag(D); x = np.zeros_like(b)
    for k in range(1, max_iter + 1):
        x = (b - R @ x) / D
        if np.linalg.norm(b - A @ x) <= tol * np.linalg.norm(b):
            return x, k
    return x, max_iter

n = 30
T = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
f = np.ones(n)
print(jacobi(T, f)[1], gauss_seidel(T, f, tol=1e-8, max_iter=100_000)[1])   # 3564 1783

# Spectral radii of the iteration matrices: rho_GS = rho_J^2 for this matrix
D, L, U = np.diag(np.diag(T)), np.tril(T, -1), np.triu(T, 1)
rho = lambda G: np.max(np.abs(np.linalg.eigvals(G)))
rho_J = rho(-np.linalg.solve(D, L + U))
rho_GS = rho(-np.linalg.solve(D + L, U))
print(round(rho_J, 6), round(rho_GS, 6), round(rho_J**2, 6))   # 0.994869 0.989765 0.989765

Both counts are huge for only 30 unknowns. That is the motivation for SOR and for conjugate gradients.

Key takeaways

• Gauss-Seidel uses each new component immediately: $(D + L)\mathbf{x}^{(k+1)} = \mathbf{b} - U\mathbf{x}^{(k)}$.
• For SPD matrices it is exact coordinate descent on $\tfrac12\mathbf{x}^TA\mathbf{x} - \mathbf{b}^T\mathbf{x}$, so it always converges.
• For consistently ordered matrices $\rho(G_{GS}) = \rho(G_J)^2$: about half the sweeps of Jacobi.
• Sweeps are sequential; multicoloring (red-black) restores parallelism.

🧪 Practice

1. Perform two Gauss-Seidel sweeps for $4x_1 - x_2 = 3$, $-x_1 + 4x_2 = 3$ from $\mathbf{0}$ and compare the error with that of two Jacobi sweeps.
2. Verify that the Gauss-Seidel update for $x_i$ is the exact minimizer of $\varphi(\mathbf{x}) = \tfrac12\mathbf{x}^TA\mathbf{x} - \mathbf{b}^T\mathbf{x}$ over $x_i$ with the other components fixed.
3. For the $2 \times 2$ matrix $\begin{bmatrix} 2 & -1 \ -1 & 2 \end{bmatrix}$, compute $G_J$ and $G_{GS}$ explicitly, find their spectral radii, and confirm $\rho(G_{GS}) = \rho(G_J)^2$.
4. Interview: You need to parallelize Gauss-Seidel for a 2D five-point stencil. How do you do it, and does the result still converge at the Gauss-Seidel rate? Hint: red-black (checkerboard) ordering makes each half-sweep embarrassingly parallel; it is Gauss-Seidel for a reordered matrix, which is still consistently ordered, so the rate is preserved.

Successive Over-Relaxation

Why it matters. On grid problems, Gauss-Seidel moves in the right direction but takes steps that are too short: the error decays smoothly and slowly. Successive over-relaxation (SOR) lengthens every Gauss-Seidel step by a factor $\omega$. With the right $\omega$, the iteration count drops from $O(n^2)$ to $O(n)$ for 1D Poisson, a dramatic gain for a one-line change. SOR was the workhorse of scientific computing before conjugate gradients and multigrid, and its symmetric version remains a common preconditioner.

Analogy. If every correction you make undershoots the target by the same amount, deliberately overshoot. Push each correction a bit further than Gauss-Seidel suggests.

The method. Compute the Gauss-Seidel value $\tilde x_i$, then move past it:

$$x_i^{(k+1)} = (1 - \omega),x_i^{(k)} + \omega,\tilde x_i, \qquad \tilde x_i = \frac{1}{a_{ii}}\Big(b_i - \sum_{j < i}a_{ij}x_j^{(k+1)} - \sum_{j > i}a_{ij}x_j^{(k)}\Big).$$

Here $\omega = 1$ is Gauss-Seidel, $\omega > 1$ is over-relaxation, and $\omega < 1$ is under-relaxation, which is occasionally useful for stabilizing a divergent iteration. In matrix form,

$$(D + \omega L),\mathbf{x}^{(k+1)} = \omega\mathbf{b} - \big(\omega U + (\omega - 1)D\big)\mathbf{x}^{(k)}.$$

Choosing $\omega$.

• Kahan's theorem: $\rho(G_\omega) \ge \lvert \omega - 1 \rvert$ for every matrix, so SOR can converge only for $0 < \omega < 2$.
• Ostrowski-Reich theorem: if $A$ is SPD, SOR converges for every $0 < \omega < 2$.
• Young's theorem: for consistently ordered matrices whose Jacobi iteration matrix has real eigenvalues with $\rho_J = \rho(G_J) < 1$, the optimal parameter and the resulting rate are

$$\omega^* = \frac{2}{1 + \sqrt{1 - \rho_J^2}}, \qquad \rho(G_{\omega^}) = \omega^ - 1.$$

How big the gain is. For the 1D Poisson matrix with mesh width $h = 1/(n + 1)$, $\rho_J = \cos\pi h$, and:

Method	Spectral radius	Sweeps to gain one digit
Jacobi	$\cos\pi h \approx 1 - \tfrac{\pi^2h^2}{2}$	$O(h^{-2}) = O(n^2)$
Gauss-Seidel	$\cos^2\pi h \approx 1 - \pi^2h^2$	half of Jacobi, still $O(n^2)$
SOR with $\omega^*$	$\approx 1 - 2\pi h$	$O(h^{-1}) = O(n)$

The rate is sensitive to $\omega$: overestimating $\omega^*$ costs less than underestimating it, so when in doubt, err slightly high.

Symmetric SOR (SSOR). A forward SOR sweep followed by a backward one gives a symmetric iteration matrix. As a standalone solver SSOR is rarely faster than SOR, but it is a good preconditioner for CG, because the preconditioner it defines is SPD.

import numpy as np

def sor(A, b, omega, tol=1e-8, max_iter=100_000):
    n = len(b)
    x = np.zeros_like(b)
    for k in range(1, max_iter + 1):
        for i in range(n):
            gs = (b[i] - A[i, :i] @ x[:i] - A[i, i+1:] @ x[i+1:]) / A[i, i]   # Gauss-Seidel value
            x[i] = (1 - omega) * x[i] + omega * gs                            # over-relax toward it
        if np.linalg.norm(b - A @ x) <= tol * np.linalg.norm(b):
            return x, k
    return x, max_iter

n = 30
T = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
f = np.ones(n)

rho_J = np.cos(np.pi / (n + 1))               # Jacobi spectral radius for this matrix (known formula)
omega_opt = 2 / (1 + np.sqrt(1 - rho_J**2))
print(round(omega_opt, 4))                    # 1.8163

for omega in (0.5, 1.0, 1.5, 1.7, 1.8, omega_opt, 1.95):
    print(round(omega, 4), sor(T, f, omega)[1])
# 0.5 5359                                   under-relaxation: slower than Gauss-Seidel
# 1.0 1783                                   Gauss-Seidel
# 1.5 585
# 1.7 293
# 1.8 152
# 1.8163 115                                 optimal: 15 times fewer sweeps than Gauss-Seidel
# 1.95 399                                   too large: slower again

Key takeaways

• SOR blends the old value with the Gauss-Seidel value: $x_i \leftarrow (1 - \omega)x_i + \omega\tilde x_i$; $\omega = 1$ is Gauss-Seidel.
• Convergence requires $0 < \omega < 2$ (Kahan), which is also sufficient for SPD matrices (Ostrowski-Reich).
• For consistently ordered matrices, $\omega^* = 2/(1 + \sqrt{1 - \rho_J^2})$ and the rate is $\omega^* - 1$; for Poisson this turns $O(n^2)$ sweeps into $O(n)$.
• The optimal $\omega$ depends on $\rho_J$, which is rarely known exactly; SSOR is mainly used as a preconditioner.

🧪 Practice

1. For $\rho_J = 0.9$, $0.99$, and $0.999$, compute $\omega^$ and $\rho(G_{\omega^})$, and estimate how many sweeps SOR and Gauss-Seidel need to reduce the error by $10^{-6}$.
2. Prove Kahan's bound: $\det G_\omega = (1 - \omega)^n$, so the product of the eigenvalues forces $\rho(G_\omega) \ge \lvert 1 - \omega \rvert$. (Use the matrix form and the fact that $D + \omega L$ is triangular.)
3. Run SOR on the $3 \times 3$ system of the Jacobi topic for $\omega = 0.8, 1.0, 1.1, 1.3$. Is over-relaxation always beneficial for a strongly diagonally dominant matrix? Explain.
4. Interview: You inherit an SOR solver for a heat-transfer code with $\omega = 1.9$ hard-coded. After a mesh refinement it suddenly needs ten times more iterations. What happened, and how would you choose $\omega$ robustly? Hint: $\omega^*$ depends on $\rho_J$, which approaches 1 as the mesh is refined; estimate $\rho_J$ from the observed convergence of a few Gauss-Seidel or Jacobi sweeps, or switch to CG with a preconditioner, whose behavior does not hinge on one parameter.

Convergence Criteria

Why it matters. An iterative method has two questions to answer that a direct method does not: will it converge at all, and when should it stop? The first is decided by the spectral radius of the iteration matrix. The second needs a stopping test that is cheap to compute and actually controls the error you care about. Getting either wrong produces solvers that silently diverge, stop too early with inaccurate answers, or waste hours iterating below the level of meaningful accuracy.

The error recursion. For any stationary iteration $\mathbf{x}^{(k+1)} = G\mathbf{x}^{(k)} + \mathbf{c}$ with fixed point $\mathbf{x}^* = G\mathbf{x}^* + \mathbf{c}$ (the solution of $A\mathbf{x} = \mathbf{b}$), subtracting the two equations gives

$$\mathbf{e}^{(k+1)} = G,\mathbf{e}^{(k)}, \qquad \mathbf{e}^{(k)} = G^k\mathbf{e}^{(0)}, \qquad \mathbf{e}^{(k)} = \mathbf{x}^{(k)} - \mathbf{x}^*.$$

Convergence of the method is therefore the question of whether $G^k \to 0$.

Theorem (Convergence of stationary iterations). The iteration converges for every starting vector if and only if the spectral radius

$$\rho(G) = \max{\lvert \lambda \rvert : \lambda \text{ an eigenvalue of } G}$$

satisfies $\rho(G) < 1$. Moreover, $\lVert \mathbf{e}^{(k)} \rVert$ decreases asymptotically like $\rho(G)^k$.

Proof sketch. If $G\mathbf{v} = \lambda\mathbf{v}$ with $\lvert \lambda \rvert \ge 1$, starting with $\mathbf{e}^{(0)} = \mathbf{v}$ gives $\mathbf{e}^{(k)} = \lambda^k\mathbf{v}$, which does not tend to zero. (If $\mathbf{v}$ is complex, use its real or imaginary part.) Conversely, if $\rho(G) < 1$, the Jordan form (Section 8.3) shows that each entry of $G^k$ is a polynomial in $k$ times $\lambda^k$, which tends to zero. $\blacksquare$

Rates and predictions. The asymptotic convergence rate is $-\log_{10}\rho(G)$ correct digits per iteration. To reduce the error by a factor $\tau$ takes about

$$k \approx \frac{\ln\tau}{\ln\rho(G)}$$

iterations. For Jacobi on the 1D Poisson matrix with $n = 30$, $\rho = 0.994869$ predicts 3581 sweeps to reduce the error by $10^{-8}$, which is exactly what the code observes.

Norms versus spectral radius. Any norm with $\lVert G \rVert < 1$ proves convergence, since $\lVert \mathbf{e}^{(k)} \rVert \le \lVert G \rVert^k\lVert \mathbf{e}^{(0)} \rVert$. That is how diagonal dominance arguments work. But $\rho(G) \le \lVert G \rVert$, and the gap can be large for non-normal $G$ (Section 10.1). Then the error can grow for many steps before the asymptotic decay sets in. For $G = \begin{bmatrix} 0.5 & 10 \ 0 & 0.5 \end{bmatrix}$, $\rho = 0.5$ but $\lVert G \rVert_2 \approx 10$, and the error grows tenfold before it decays.

Sufficient conditions at a glance.

Matrix class	Jacobi	Gauss-Seidel	SOR ($0 < \omega < 2$)
strictly diagonally dominant	converges	converges	converges for $0 < \omega \le 1$
symmetric positive definite	converges iff $2D - A$ is also SPD	always converges	always converges
general	check $\rho(G_J)$	check $\rho(G_{GS})$	$0 < \omega < 2$ is necessary, not sufficient

Stopping criteria. The error is unknown, so practical tests use computable quantities:

Test	Formula	Comments
relative residual	$\lVert \mathbf{b} - A\mathbf{x}^{(k)} \rVert \le \text{tol}\cdot\lVert \mathbf{b} \rVert$	the standard choice; scale-invariant in $\mathbf{b}$
normwise backward error	$\dfrac{\lVert \mathbf{r}^{(k)} \rVert}{\lVert A \rVert\lVert \mathbf{x}^{(k)} \rVert + \lVert \mathbf{b} \rVert} \le \text{tol}$	meaningful down to $u$ (Section 12.1)
increment	$\lVert \mathbf{x}^{(k+1)} - \mathbf{x}^{(k)} \rVert \le \text{tol}\cdot\lVert \mathbf{x}^{(k+1)} \rVert$	dangerous for slow methods: the error can be about $\tfrac{\rho}{1 - \rho}$ times the increment
safeguards	maximum iteration count; detect stagnation or growth	always include them

The residual is not the error. From $\mathbf{e} = -A^{-1}\mathbf{r}$,

$$\frac{\lVert \mathbf{x}^{(k)} - \mathbf{x}^* \rVert}{\lVert \mathbf{x}^* \rVert} \le \kappa(A),\frac{\lVert \mathbf{r}^{(k)} \rVert}{\lVert \mathbf{b} \rVert}.$$

A residual of $10^{-10}$ guarantees about 10 correct digits only if $\kappa(A) \approx 1$. A tolerance should also not be tighter than the accuracy the problem deserves: iterating below the discretization error of a PDE, or below $\kappa(A),u$, wastes time without improving the answer.

import numpy as np

rho = lambda G: np.max(np.abs(np.linalg.eigvals(G)))

# Predicting the iteration count from the spectral radius
n = 30
T = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
D = np.diag(np.diag(T))
G = -np.linalg.solve(D, T - D)                # Jacobi iteration matrix
c = np.linalg.solve(D, np.ones(n))
x_star = np.linalg.solve(T, np.ones(n))
x, k = np.zeros(n), 0
while np.linalg.norm(x - x_star) > 1e-8 * np.linalg.norm(x_star):
    x, k = G @ x + c, k + 1
print(round(rho(G), 6), k, round(np.log(1e-8) / np.log(rho(G))))   # 0.994869 3581 3581: actual = predicted

# rho(G) < 1 decides convergence, but a large ||G|| can cause transient growth first
G2 = np.array([[0.5, 10.], [0., 0.5]])
e = np.array([0., 1.])
norms = []
for k in range(40):
    norms.append(np.linalg.norm(e)); e = G2 @ e
print(rho(G2), round(np.linalg.norm(G2, 2), 2))           # 0.5 10.02
print([round(float(v), 2) for v in norms[:6]], f"{norms[-1]:.1e}")
# [1.0, 10.01, 10.0, 7.5, 5.0, 3.13] 1.4e-09   grows tenfold, then decays like 0.5^k

# A small residual does not guarantee a small error when A is ill conditioned
A = np.array([[1., 1.], [1., 1. + 1e-10]])
b = np.array([2., 2. + 1e-10])                # exact solution (1, 1)
x_bad = np.array([2., 0.])                    # far from (1, 1)
r = b - A @ x_bad
print(f"{np.linalg.norm(r) / np.linalg.norm(b):.1e}", np.linalg.norm(x_bad - 1) / np.linalg.norm([1., 1.]),
      f"{np.linalg.cond(A):.1e}")
# 3.5e-11 1.0 4.0e+10                        residual 3.5e-11, error 100%, kappa 4e10

Key takeaways

• A stationary iteration $\mathbf{x}^{(k+1)} = G\mathbf{x}^{(k)} + \mathbf{c}$ converges for every start if and only if $\rho(G) < 1$; the error decays asymptotically like $\rho(G)^k$.
• About $\ln\tau/\ln\rho$ iterations reduce the error by a factor $\tau$; as $\rho \to 1$, this blows up.
• $\lVert G \rVert < 1$ in some norm is a sufficient condition; non-normal $G$ can show transient growth even when $\rho(G) < 1$.
• Stop on the relative residual or backward error, with an iteration cap; remember that the error can be up to $\kappa(A)$ times the relative residual.

🧪 Practice

1. An iteration has $\rho(G) = 0.95$. How many iterations reduce the error by $10^{-6}$? What if $\rho = 0.995$?
2. Show that if $\lVert G \rVert < 1$ in some induced norm, then $\lVert \mathbf{x}^{(k)} - \mathbf{x}^* \rVert \le \frac{\lVert G \rVert}{1 - \lVert G \rVert}\lVert \mathbf{x}^{(k)} - \mathbf{x}^{(k-1)} \rVert$. Why does this make the increment test misleading when $\lVert G \rVert$ is close to 1?
3. Construct a $2 \times 2$ matrix $G$ with $\rho(G) < 1$ but $\lVert G \rVert_\infty > 1$, and compute $\lVert G^k \rVert_\infty$ for $k = 1, \dots, 10$.
4. Interview: A solver reports "converged: relative residual $10^{-12}$" but the physical answer is clearly wrong. List the possible explanations, in order of what you would check. Hint: ill conditioning ($\kappa$ times the residual), a wrong residual definition (preconditioned versus true residual), a bug in assembling $A$ or $\mathbf{b}$, or a model error; compute the true residual $\mathbf{b} - A\mathbf{x}$ yourself and estimate $\kappa$.

Conjugate Gradient Method

Why it matters. For symmetric positive definite systems, the conjugate gradient method (CG) is the method of choice when the matrix is large and sparse. It needs one matrix-vector product per iteration and four vectors of memory. It has no parameter to tune. It converges in $O(\sqrt{\kappa})$ iterations instead of the $O(\kappa)$ of simpler methods, and in exact arithmetic it would even terminate in at most $n$ steps. Hestenes and Stiefel introduced it in 1952. Together with preconditioning, it solves systems with billions of unknowns in structural analysis, reservoir simulation, and machine learning.

From linear system to minimization. For SPD $A$, solving $A\mathbf{x} = \mathbf{b}$ is the same as minimizing the convex quadratic (Section 10.3)

$$\varphi(\mathbf{x}) = \tfrac12\mathbf{x}^TA\mathbf{x} - \mathbf{b}^T\mathbf{x}, \qquad \nabla\varphi(\mathbf{x}) = A\mathbf{x} - \mathbf{b} = -\mathbf{r}.$$

The residual is the direction of steepest descent.

Steepest descent and its failure. Move along $\mathbf{r}_k$ with the exactly optimal step $\alpha_k = \mathbf{r}_k^T\mathbf{r}_k / \mathbf{r}_k^TA\mathbf{r}_k$. On a round bowl ($\kappa = 1$) this finishes in one step. On a long, narrow valley ($\kappa \gg 1$), consecutive steps are perpendicular and the iterates zigzag across the valley. The error shrinks only by the factor $(\kappa - 1)/(\kappa + 1)$ per step:

   level curves of phi for an ill-conditioned A (a long, narrow valley)

   steepest descent: each step is orthogonal        conjugate gradient: each new direction is
   to the previous one, so it zigzags               A-conjugate to all previous ones (n = 2: two steps)

     x0                                               x0
       \    /\    /\    /\                              \
        \  /  \  /  \  /  \   ...  -> x*                 \
         \/    \/    \/    \/                             '-----------------------> x*

The idea of CG: conjugate directions. Two directions are $A$-conjugate if $\mathbf{p}_i^TA\mathbf{p}_j = 0$: they are orthogonal in the inner product $\langle \mathbf{u}, \mathbf{v} \rangle_A = \mathbf{u}^TA\mathbf{v}$ (Section 9.1). Minimizing $\varphi$ along $n$ mutually conjugate directions, one after another, reaches the exact minimizer. The minimization along a new conjugate direction never spoils the minimization along the previous ones. CG builds these directions on the fly: each new direction is the current residual, corrected by a multiple of the previous direction so as to make it $A$-conjugate.

The algorithm. Start with $\mathbf{x}_0$, $\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0$, $\mathbf{p}_0 = \mathbf{r}_0$. For $k = 0, 1, 2, \dots$:

$$\begin{aligned} \alpha_k &= \frac{\mathbf{r}_k^T\mathbf{r}_k}{\mathbf{p}_k^TA\mathbf{p}k} &&\text{(exact line search)} \ \mathbf{x}{k+1} &= \mathbf{x}_k + \alpha_k\mathbf{p}k, \qquad \mathbf{r}{k+1} = \mathbf{r}k - \alpha_kA\mathbf{p}k \ \beta_k &= \frac{\mathbf{r}{k+1}^T\mathbf{r}{k+1}}{\mathbf{r}k^T\mathbf{r}k}, \qquad \mathbf{p}{k+1} = \mathbf{r}{k+1} + \beta_k\mathbf{p}_k &&\text{(new conjugate direction)} \end{aligned}$$

Each iteration costs one product $A\mathbf{p}_k$, two dot products, and three vector updates.

What CG guarantees (in exact arithmetic):

• The residuals are mutually orthogonal, and the directions are mutually $A$-conjugate.
• $\mathbf{x}_k$ is the best approximation in the Krylov subspace (Section 12.3): it minimizes the $A$-norm of the error $\lVert \mathbf{x} - \mathbf{x}^* \rVert_A = \sqrt{(\mathbf{x} - \mathbf{x}^)^TA(\mathbf{x} - \mathbf{x}^)}$ over $\mathbf{x}_0 + \operatorname{span}{\mathbf{r}_0, A\mathbf{r}_0, \dots, A^{k-1}\mathbf{r}_0}$.
• It terminates with the exact solution in at most $n$ steps, and in at most $d$ steps if $A$ has only $d$ distinct eigenvalues.
• Convergence bound, with $\kappa = \lambda_{\max}/\lambda_{\min}$:

$$\frac{\lVert \mathbf{x}_k - \mathbf{x}^* \rVert_A}{\lVert \mathbf{x}_0 - \mathbf{x}^* \rVert_A} \le 2\left(\frac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\right)^k.$$

So the iteration count grows like $\sqrt\kappa$, compared with $\kappa$ for steepest descent. For the 2D Poisson matrix, $\kappa = O(h^{-2})$, so CG needs $O(h^{-1})$ iterations. The bound is often pessimistic, because CG also profits from clustered eigenvalues.

In floating point, the residuals and directions gradually lose orthogonality and conjugacy (the same effect as in Lanczos, Section 12.3). Finite termination is lost, so CG is used as an iterative method with a stopping test, and badly conditioned problems can need more than $n$ iterations.

Beyond SPD. CG requires $A$ symmetric positive definite. Other Krylov methods cover the rest:

Method	Matrices	Notes
CG	SPD	short recurrence, minimizes the $A$-norm of the error
MINRES	symmetric indefinite	short recurrence, minimizes the residual
GMRES	general	minimizes the residual; memory grows with $k$, so it is restarted
BiCGSTAB	general	short recurrences, no optimality, sometimes erratic
LSQR / CGLS	least squares $\min\lVert A\mathbf{x} - \mathbf{b} \rVert$	CG on the normal equations without forming $A^TA$

import numpy as np
import scipy.sparse as sp

def cg(A, b, tol=1e-8, max_iter=10_000):
    x = np.zeros_like(b)
    r = b.copy()                              # residual b - A x, also minus the gradient of phi
    p = r.copy()                              # first search direction: steepest descent
    rr = r @ r
    for k in range(1, max_iter + 1):
        Ap = A @ p                            # the only use of A: one matrix-vector product
        alpha = rr / (p @ Ap)                 # exact line search along p
        x += alpha * p
        r -= alpha * Ap
        rr_new = r @ r
        if np.sqrt(rr_new) <= tol * np.linalg.norm(b):
            return x, k
        p = r + (rr_new / rr) * p             # new direction, A-conjugate to all previous ones
        rr = rr_new
    return x, max_iter

def steepest_descent(A, b, tol=1e-8, max_iter=200_000):
    x = np.zeros_like(b); r = b.copy()
    for k in range(1, max_iter + 1):
        alpha = (r @ r) / (r @ (A @ r))
        x += alpha * r
        r = b - A @ x
        if np.linalg.norm(r) <= tol * np.linalg.norm(b):
            return x, k
    return x, max_iter

# 2D Poisson matrix on a 50 x 50 grid (n = 2500, SPD)
m = 50
T = sp.diags([-1., 2., -1.], [-1, 0, 1], shape=(m, m))
A = (sp.kron(sp.identity(m), T) + sp.kron(T, sp.identity(m))).tocsr()
b = np.ones(m * m)
x, k_cg = cg(A, b)
print(k_cg, steepest_descent(A, b)[1], f"{np.linalg.norm(A @ x - b) / np.linalg.norm(b):.1e}")
# 93 9711 8.4e-09                            CG: 93 iterations; steepest descent: 9711

# Theory: error reduction per step is about (sqrt(kappa) - 1) / (sqrt(kappa) + 1)
lam = np.linalg.eigvalsh(A.toarray())
kappa = lam[-1] / lam[0]
print(round(kappa), round(np.log(2 / 1e-8) / -np.log((np.sqrt(kappa) - 1) / (np.sqrt(kappa) + 1))))
# 1053 310                                   kappa, and the iteration count the bound predicts (pessimistic)

# Few distinct eigenvalues: CG finishes in that many steps, whatever n is
rng = np.random.default_rng(0)
Q, _ = np.linalg.qr(rng.standard_normal((500, 500)))
B = Q @ np.diag(rng.choice([1., 4., 9., 16.], 500)) @ Q.T    # 500 x 500, only 4 distinct eigenvalues
print(cg(B, rng.standard_normal(500), tol=1e-12)[1])          # 4

Key takeaways

• For SPD $A$, CG minimizes $\tfrac12\mathbf{x}^TA\mathbf{x} - \mathbf{b}^T\mathbf{x}$ along $A$-conjugate directions built from the residuals.
• Each iteration needs one matrix-vector product and $O(n)$ extra work; memory is four vectors; there are no parameters.
• $\mathbf{x}_k$ is optimal in the $A$-norm over a Krylov subspace; the iteration count scales like $\sqrt\kappa$, and clustered eigenvalues help further.
• Exact termination in $n$ steps is lost in floating point; CG is used as an iterative method.
• For non-SPD problems use MINRES, GMRES, BiCGSTAB, or LSQR.

🧪 Practice

1. Run CG by hand on $\begin{bmatrix} 4 & 1 \ 1 & 3 \end{bmatrix}\mathbf{x} = \begin{bmatrix} 1 \ 2 \end{bmatrix}$ from $\mathbf{x}_0 = \mathbf{0}$. Verify that it finishes in two steps and that $\mathbf{p}_0^TA\mathbf{p}_1 = 0$.
2. Show that if $\mathbf{p}0, \dots, \mathbf{p}{n-1}$ are nonzero and mutually $A$-conjugate, they are linearly independent, and that $\mathbf{x}^* = \sum_k \frac{\mathbf{p}_k^T\mathbf{b}}{\mathbf{p}_k^TA\mathbf{p}_k}\mathbf{p}_k$.
3. Using the convergence bound, estimate the number of CG iterations needed for a $10^{-8}$ error reduction when $\kappa = 10^2$, $10^4$, and $10^6$, and compare with steepest descent.
4. Interview: Your colleague applies CG to a symmetric matrix and it diverges, and to a nonsymmetric one and it stagnates. What went wrong in each case, and what would you use instead? Hint: CG needs positive definiteness (a negative eigenvalue makes $\mathbf{p}^TA\mathbf{p}$ possibly zero or negative, breaking the line search) and symmetry (otherwise conjugacy and the short recurrence fail); MINRES for symmetric indefinite, GMRES or BiCGSTAB for nonsymmetric.

Preconditioning

Why it matters. Krylov methods converge at a rate set by the spectrum of $A$: roughly $\sqrt\kappa$ iterations for CG, and similar behavior for GMRES. Real problems from fine meshes, multiscale physics, or badly scaled variables have huge condition numbers, and unpreconditioned iterations can need thousands of steps or never converge. A preconditioner transforms the system into an equivalent one with a better spectrum. In practice, the choice of preconditioner matters more than the choice of Krylov method.

Analogy. Steepest descent and CG struggle in a long, narrow valley. A preconditioner is a change of coordinates that stretches the valley back into a round bowl. A perfect preconditioner makes the bowl exactly round, but finding it is as hard as solving the problem. A good one makes it round enough, cheaply.

The idea. Choose a matrix $M \approx A$ such that systems $M\mathbf{z} = \mathbf{r}$ are cheap to solve. Then solve an equivalent system:

Variant	System solved	Remarks
left	$M^{-1}A\mathbf{x} = M^{-1}\mathbf{b}$	GMRES then minimizes the preconditioned residual
right	$AM^{-1}\mathbf{y} = \mathbf{b}$, $\mathbf{x} = M^{-1}\mathbf{y}$	true residual is minimized
split, $M = LL^T$	$L^{-1}AL^{-T}\mathbf{y} = L^{-1}\mathbf{b}$, $\mathbf{x} = L^{-T}\mathbf{y}$	keeps symmetry: the basis of preconditioned CG

Two goals pull against each other. $M = A$ gives convergence in one iteration, but applying $M^{-1}$ is the original problem. $M = I$ is free but does nothing. A good preconditioner makes the eigenvalues of $M^{-1}A$ clustered (ideally near 1), and costs $O(\operatorname{nnz})$ to build and apply.

Preconditioned CG (PCG). For SPD $A$ and SPD $M$, CG on the split system can be rewritten so that it only needs products with $A$ and solves with $M$, never the factor $L$. The changes to plain CG are one extra line, $\mathbf{z}_k = M^{-1}\mathbf{r}_k$, and $\mathbf{r}^T\mathbf{z}$ in place of $\mathbf{r}^T\mathbf{r}$. The convergence bound becomes the CG bound with $\kappa(M^{-1}A)$ in place of $\kappa(A)$. The preconditioner must be SPD, which rules out plain Gauss-Seidel or SOR but allows symmetric Gauss-Seidel and SSOR.

The standard preconditioners.

Preconditioner	$M$	Cost to apply	Strength
Jacobi (diagonal)	$\operatorname{diag}(A)$	$n$ divisions	fixes bad scaling; trivially parallel
block Jacobi	diagonal blocks of $A$	small block solves	variables grouped by node or subdomain
SSOR	$(D + \omega L)D^{-1}(D + \omega U)$ (up to scaling)	two triangular sweeps	cheap, symmetric
incomplete Cholesky / ILU	$\tilde L\tilde U \approx A$, dropping fill-in by position (IC(0), ILU(0)) or size (ILUT)	two sparse triangular solves	general purpose; quality tunable via drop tolerance
multigrid (geometric or algebraic)	one V-cycle	$O(n)$	$\kappa(M^{-1}A) = O(1)$ for elliptic PDEs: iteration counts independent of mesh size
domain decomposition	local subdomain solves plus a coarse problem	parallel local solves	scales on large parallel machines

The ILU drop tolerance shows the trade-off clearly. A smaller tolerance means more fill, so more memory and setup time, but fewer iterations. The best total time usually lies in between.

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import gmres, spilu, LinearOperator

def pcg(A, b, M_solve, tol=1e-8, max_iter=10_000):
    """Preconditioned CG: M_solve(r) applies M^{-1}, where M is SPD and M ~ A."""
    x = np.zeros_like(b)
    r = b.copy()
    z = M_solve(r)                            # preconditioned residual
    p = z.copy()
    rz = r @ z
    for k in range(1, max_iter + 1):
        Ap = A @ p
        alpha = rz / (p @ Ap)
        x += alpha * p
        r -= alpha * Ap
        if np.linalg.norm(r) <= tol * np.linalg.norm(b):
            return x, k
        z = M_solve(r)
        rz_new = r @ z
        p = z + (rz_new / rz) * p
        rz = rz_new
    return x, max_iter

# A badly scaled SPD matrix: 1D Poisson with unknowns measured in wildly different units
n = 200
T = sp.diags([-1., 2., -1.], [-1, 0, 1], shape=(n, n))
S = sp.diags(np.logspace(0, 3, n))            # scales from 1 to 1000
A = (S @ T @ S).tocsr()
b = A @ np.ones(n)
d = A.diagonal()
print(pcg(A, b, lambda r: r)[1],              # no preconditioner (M = I)
      pcg(A, b, lambda r: r / d)[1])          # Jacobi preconditioner (M = diag(A))
# 6253 200                                   without: far more than n steps (rounding destroys conjugacy)

# A nonsymmetric convection-diffusion matrix on a 60 x 60 grid: GMRES with and without ILU
m = 60
I = sp.identity(m)
T1 = sp.diags([-1., 2., -1.], [-1, 0, 1], shape=(m, m))
C1 = sp.diags([-1., 1.], [-1, 0], shape=(m, m))         # upwind first derivative
K = (sp.kron(I, T1) + sp.kron(T1, I) + 0.5 * sp.kron(I, C1)).tocsc()
f = np.ones(m * m)

def gmres_iterations(M=None):
    count = []
    x, info = gmres(K, f, M=M, rtol=1e-8, restart=50, maxiter=1000,
                    callback=lambda res: count.append(res), callback_type="pr_norm")
    return len(count), bool(np.linalg.norm(K @ x - f) / np.linalg.norm(f) < 1e-7)

print("none", gmres_iterations())
for drop in (1e-1, 1e-2, 1e-4):
    ilu = spilu(K, drop_tol=drop, fill_factor=20)        # incomplete LU: drops small fill-in entries
    M = LinearOperator(K.shape, ilu.solve)               # applying M^{-1} = two sparse triangular solves
    print(drop, ilu.L.nnz + ilu.U.nnz, gmres_iterations(M))
print("nnz(K) =", K.nnz)
# none (401, True)                           GMRES iterations, converged?
# 0.1 39820 (35, True)                       ILU factors with about 2x nnz(K): 11x fewer iterations
# 0.01 77597 (12, True)
# 0.0001 144054 (4, True)                    nearly a full LU: few iterations, much more memory
# nnz(K) = 17760

Key takeaways

• A preconditioner $M \approx A$, cheap to invert, replaces the system by one whose spectrum (ideally $M^{-1}A \approx I$) gives fast Krylov convergence.
• PCG needs $M$ SPD; it changes plain CG by one preconditioner solve per iteration, and its rate depends on $\kappa(M^{-1}A)$.
• Jacobi fixes scaling, ILU/IC are general-purpose, multigrid gives mesh-independent iteration counts for elliptic PDEs.
• Stronger preconditioners cost more memory and setup; the best choice minimizes total time, not iteration count.

🧪 Practice

1. Show that for $A = D_sTD_s$, with $D_s$ a positive diagonal scaling and $T$ SPD, the Jacobi-preconditioned matrix $D^{-1}A$ is similar to $\tilde D^{-1}T$, where $\tilde D = \operatorname{diag}(T)$, so the bad scaling disappears.
2. Derive PCG from CG applied to $L^{-1}AL^{-T}\mathbf{y} = L^{-1}\mathbf{b}$ with $M = LL^T$, and check that $L$ itself never appears, only $M^{-1}$.
3. For the convection-diffusion example, measure the total time (factorization plus solve) for several ILU drop tolerances and find the best one. Is it the one with the fewest iterations?
4. Interview: A CG solver for a structural mechanics model needs 3000 iterations, and the count grows every time the mesh is refined. What would you try, in which order, and why? Hint: start with diagonal scaling (cheap), then incomplete Cholesky; for mesh-independent counts on elliptic problems, algebraic multigrid; check symmetry and positive definiteness of the preconditioner, and compare total time rather than iterations.

<a id="123-eigenvalue-algorithms"></a>

12.3 Eigenvalue Algorithms

Eigenvalues are roots of a polynomial of degree $n$, so for $n \ge 5$ no finite formula computes them (Section 8.1), and every eigenvalue algorithm is iterative. This subchapter builds them up from repeated multiplication by $A$ (the power method), through shifts and inversion, to the QR algorithm behind eig and eigh, and to the Krylov methods behind scipy.sparse.linalg.eigs and eigsh.

Power Method

Why it matters. Many applications need only the dominant eigenvalue, the one of largest modulus, and its eigenvector. Examples are the PageRank vector (Section 8.4), the long-run behavior of $\mathbf{x}_{k+1} = A\mathbf{x}k$ (Section 8.4), the spectral norm $\lVert A \rVert_2 = \sqrt{\lambda{\max}(A^TA)}$, and the first principal component. The power method finds it using nothing but matrix-vector products. It works for matrices too large to store densely, and even for matrices available only as a function that computes $A\mathbf{x}$. It is also the conceptual seed of every other algorithm in this subchapter.

Analogy. Repeatedly applying $A$ is like repeatedly shouting into a canyon with many echo paths. Each echo is amplified by a different factor, and after a few rounds only the most amplified path, the dominant eigenvector, is still audible.

The method. Start with a vector $\mathbf{x}_0$ and repeat

$$\mathbf{y} = A\mathbf{x}k, \qquad \mathbf{x}{k+1} = \frac{\mathbf{y}}{\lVert \mathbf{y} \rVert}, \qquad \mu_{k+1} = \mathbf{x}{k+1}^TA\mathbf{x}{k+1}.$$

The normalization prevents overflow or underflow. $\mu_k$, the Rayleigh quotient (Section 10.3), estimates the eigenvalue. For a nonsymmetric matrix, the conjugate transpose $\mathbf{x}^*$ replaces $\mathbf{x}^T$.

Why it converges. Suppose $A$ is diagonalizable with eigenvalues ordered $\lvert \lambda_1 \rvert > \lvert \lambda_2 \rvert \ge \cdots \ge \lvert \lambda_n \rvert$ and eigenvectors $\mathbf{v}_i$, and write $\mathbf{x}_0 = c_1\mathbf{v}_1 + \cdots + c_n\mathbf{v}_n$ with $c_1 \ne 0$. Then

$$A^k\mathbf{x}_0 = \lambda_1^k\Big(c_1\mathbf{v}1 + \sum{i \ge 2}c_i\Big(\frac{\lambda_i}{\lambda_1}\Big)^k\mathbf{v}_i\Big).$$

Every ratio $\lvert \lambda_i/\lambda_1 \rvert < 1$, so after normalization $\mathbf{x}_k$ approaches $\pm\mathbf{v}_1$ (or a unimodular multiple of it in the complex case). The error decays like $\lvert \lambda_2/\lambda_1 \rvert^k$. For a symmetric matrix, the Rayleigh quotient has quadratic accuracy (Section 10.3), so the eigenvalue error decays like $\lvert \lambda_2/\lambda_1 \rvert^{2k}$.

Practical facts.

Issue	What happens	Remedy
ratio $\lvert \lambda_2/\lambda_1 \rvert$ close to 1	very slow convergence	shifts or inverse iteration (next topic); Krylov methods
$\lvert \lambda_1 \rvert = \lvert \lambda_2 \rvert$ (e.g. $\pm\lambda$, or a complex pair of a real matrix)	the iterates oscillate and never settle	use a Krylov method or the QR algorithm
$c_1 = 0$ in the starting vector	in exact arithmetic, convergence to $\mathbf{v}_2$	random start; rounding errors usually reintroduce $\mathbf{v}_1$
need $\lambda_2$ as well	the power method finds only $\lambda_1$	deflation: for symmetric $A$, iterate with $A - \lambda_1\mathbf{v}_1\mathbf{v}_1^T$, or iterate on a block of vectors (next topics)
cost	one matrix-vector product per step	ideal for sparse or matrix-free $A$

For PageRank, the second eigenvalue of the Google matrix is at most the damping factor $0.85$ in modulus, so about 100 iterations reach seven digits ($0.85^{100} \approx 10^{-7}$) regardless of the size of the web graph.

import numpy as np

rng = np.random.default_rng(0)
Q, _ = np.linalg.qr(rng.standard_normal((5, 5)))
A = Q @ np.diag([10., 9., 4., 2., 1.]) @ Q.T  # symmetric, known eigenvalues, ratio 9/10

def power_method(A, x0, iters):
    x = x0 / np.linalg.norm(x0)
    history = []
    for _ in range(iters):
        y = A @ x
        x = y / np.linalg.norm(y)             # normalize to avoid overflow or underflow
        history.append(x @ A @ x)             # Rayleigh quotient estimate of lambda_1
    return x, history

x, hist = power_method(A, rng.standard_normal(5), 100)
q1 = Q[:, 0]
for k in (10, 20, 40, 80):
    print(k, f"{abs(hist[k - 1] - 10):.1e}")  # error in the eigenvalue estimate
# 10 6.5e-02
# 20 8.3e-03
# 40 1.2e-04
# 80 2.7e-08                                 40 more steps gain a factor (0.9^2)^40 ~ 2e-4
print(f"{1 - abs(x @ q1):.1e}", f"{0.9**200:.1e}")   # 2.0e-10 7.1e-10: same order as (9/10)^(2k)

# Faster when lambda_2 / lambda_1 is small
B = Q @ np.diag([10., 2., 1.5, 1., 0.5]) @ Q.T
_, hist_b = power_method(B, rng.standard_normal(5), 20)
print(f"{abs(hist_b[-1] - 10):.1e}")          # 1.8e-15: ratio 0.2, machine precision in 20 steps

# Failure: equal magnitudes (eigenvalues 1 and -1) -> the iterates oscillate forever
R = np.array([[0., 1.], [1., 0.]])
x = np.array([1., 0.3])
for _ in range(4):
    x = R @ x / np.linalg.norm(R @ x)
    print(np.round(x, 4))
# [0.2873 0.9578]
# [0.9578 0.2873]
# [0.2873 0.9578]
# [0.9578 0.2873]

Key takeaways

• The power method repeats $\mathbf{x} \leftarrow A\mathbf{x}/\lVert A\mathbf{x} \rVert$ and converges to the dominant eigenvector when $\lvert \lambda_1 \rvert > \lvert \lambda_2 \rvert$.
• The error decays like $\lvert \lambda_2/\lambda_1 \rvert^k$; for symmetric matrices the Rayleigh quotient eigenvalue estimate converges twice as fast in digits.
• It needs only matrix-vector products, so it suits huge sparse or matrix-free problems such as PageRank.
• It fails or crawls when the dominant eigenvalue is not well separated in modulus.

🧪 Practice

1. Run three steps of the power method by hand on $\begin{bmatrix} 2 & 1 \ 1 & 2 \end{bmatrix}$ from $\mathbf{x}_0 = (1, 0)$ and compare the Rayleigh quotients with $\lambda_1 = 3$.
2. A matrix has eigenvalues $5, 4.9, 1$. Estimate how many power iterations are needed to get the dominant eigenvector to 6 digits. How much does the shift $A - 2.45I$ help, and why can no shift help much here?
3. Explain how the power method applied to $A^TA$ estimates $\lVert A \rVert_2$ without forming $A^TA$, and what it costs per iteration.
4. Interview: Why does PageRank converge in about 100 iterations for a graph with billions of nodes, and what would happen without the damping factor? Hint: the damping factor bounds $\lvert \lambda_2 \rvert \le 0.85$ independently of the graph; without it, a disconnected or periodic graph can have $\lvert \lambda_2 \rvert = 1$, so the iteration never converges or the limit is not unique.

Inverse and Shifted Power Iteration

Why it matters. The power method finds only the eigenvalue of largest modulus, and only at the rate $\lvert \lambda_2/\lambda_1 \rvert$. Often you need a specific eigenvalue instead: the smallest (the fundamental vibration frequency), the one nearest a given value (a resonance near the operating frequency), or the eigenvector for an eigenvalue that is already known. Applying the power method to $(A - \sigma I)^{-1}$ turns any eigenvalue into the dominant one, and updating the shift on the fly gives one of the fastest converging iterations in numerical analysis.

Inverse iteration. If $A\mathbf{v} = \lambda\mathbf{v}$, then $(A - \sigma I)^{-1}\mathbf{v} = \frac{1}{\lambda - \sigma}\mathbf{v}$. The eigenvalue of $A$ closest to the shift $\sigma$ becomes the dominant eigenvalue of $(A - \sigma I)^{-1}$. The power method on $(A - \sigma I)^{-1}$ is

$$\text{solve } (A - \sigma I)\mathbf{y} = \mathbf{x}k, \qquad \mathbf{x}{k+1} = \mathbf{y}/\lVert \mathbf{y} \rVert,$$

and it converges to the eigenvector of the eigenvalue $\lambda_J$ nearest $\sigma$, at the rate

$$\frac{\lvert \lambda_J - \sigma \rvert}{\lvert \lambda_K - \sigma \rvert}, \qquad \lambda_K = \text{the second-closest eigenvalue to } \sigma.$$

With $\sigma = 0$ this finds the smallest eigenvalue in modulus. A shift close to $\lambda_J$ makes the ratio tiny and the convergence very fast.

Never invert: factor once. Each step is a linear solve with the same matrix $A - \sigma I$. Factor it once (LU, Cholesky, or a sparse factorization, Section 12.1) and each step costs only two triangular solves.

"But the system is nearly singular." When $\sigma \approx \lambda_J$, the matrix $A - \sigma I$ is badly conditioned, and the computed $\mathbf{y}$ has a large error. That error lies almost entirely in the direction of $\mathbf{v}_J$, which is the direction we want, and it disappears in the normalization. Inverse iteration with a nearly exact eigenvalue as the shift typically delivers the eigenvector in one or two steps. This is how LAPACK computes eigenvectors of selected eigenvalues once they are known.

Rayleigh quotient iteration (RQI). Use the best available eigenvalue estimate as the shift at every step:

$$\sigma_k = R_A(\mathbf{x}_k) = \frac{\mathbf{x}_k^TA\mathbf{x}_k}{\mathbf{x}_k^T\mathbf{x}_k}, \qquad \text{solve } (A - \sigma_kI)\mathbf{y} = \mathbf{x}k, \qquad \mathbf{x}{k+1} = \mathbf{y}/\lVert \mathbf{y} \rVert.$$

For symmetric $A$, an eigenvector error $\varepsilon$ gives a shift error $O(\varepsilon^2)$ (Section 10.3). The next step then reduces the eigenvector error by a factor proportional to the shift error, so $\varepsilon_{k+1} = O(\varepsilon_k^3)$: cubic convergence. The number of correct digits roughly triples per step. For nonsymmetric matrices the convergence is quadratic.

Method	Shift	Cost per step	Convergence
power method	none	one product $A\mathbf{x}$	linear, ratio $\lvert \lambda_2/\lambda_1 \rvert$
inverse iteration	fixed $\sigma$	two triangular solves (after one factorization)	linear, ratio $\lvert \lambda_J - \sigma \rvert / \lvert \lambda_K - \sigma \rvert$
Rayleigh quotient iteration	$\sigma_k = R_A(\mathbf{x}_k)$	a new factorization every step	cubic (symmetric), quadratic (general)

The price of RQI is a new factorization per step, and less control: it converges to some eigenvalue near the starting estimate, not necessarily the one you wanted. A common strategy is a few steps of inverse iteration with a fixed shift to approach the target, then RQI to finish.

import numpy as np
from scipy.linalg import lu_factor, lu_solve

rng = np.random.default_rng(0)
Q, _ = np.linalg.qr(rng.standard_normal((6, 6)))
lams = np.array([1., 2., 3.1, 3.5, 7., 10.])
A = Q @ np.diag(lams) @ Q.T

def shifted_inverse_iteration(A, sigma, x0, iters):
    lu = lu_factor(A - sigma * np.eye(len(A)))   # factor once: each step is then two triangular solves
    x = x0 / np.linalg.norm(x0)
    for _ in range(iters):
        y = lu_solve(lu, x)                       # y = (A - sigma I)^{-1} x
        x = y / np.linalg.norm(y)
    return x @ A @ x, x

x0 = rng.standard_normal(6)
for sigma in (0.0, 3.0, 3.4, 6.0):
    lam, _ = shifted_inverse_iteration(A, sigma, x0, 10)
    print(sigma, round(lam, 10))              # converges to the eigenvalue nearest the shift
# 0.0 1.0000014501                           ratio 1/2: slow
# 3.0 3.1                                    ratio 0.1/0.5: fast
# 3.4 3.4999999998                           ratio 0.1/0.3
# 6.0 6.9999999556                           ratio 1/2.5

# Rayleigh quotient iteration: update the shift every step
def rqi(A, x0, iters):
    x = x0 / np.linalg.norm(x0)
    sigma = x @ A @ x
    for k in range(iters):
        try:
            y = np.linalg.solve(A - sigma * np.eye(len(A)), x)
        except np.linalg.LinAlgError:         # shift hit an eigenvalue exactly: converged
            break
        x = y / np.linalg.norm(y)
        sigma = x @ A @ x
        res = np.linalg.norm(A @ x - sigma * x)
        print(k + 1, f"{sigma:.15f}", f"{res:.1e}")

rqi(A, Q[:, 2] + 0.1 * rng.standard_normal(6), 4)   # start near the eigenvector for 3.1
# 1 3.105563640462558 1.8e-01
# 2 3.100000020655925 1.9e-04
# 3 3.100000000000001 4.0e-12                residual 1.9e-4 -> 4e-12: about (1.9e-4)^3, cubic
# 4 3.100000000000001 8.5e-16

Key takeaways

• Inverse iteration is the power method applied to $(A - \sigma I)^{-1}$; it converges to the eigenvalue nearest $\sigma$, at the rate of the ratio of the nearest and second-nearest distances to $\sigma$.
• Factor $A - \sigma I$ once; near-singularity is harmless, because the error points along the wanted eigenvector.
• Rayleigh quotient iteration updates the shift with $R_A(\mathbf{x}_k)$ and converges cubically for symmetric matrices, at the cost of a new factorization per step.
• Inverse iteration is the standard way to get eigenvectors of eigenvalues that are already known.

🧪 Practice

1. For eigenvalues $1, 2, 3.1, 3.5, 7, 10$, compute the convergence ratio of inverse iteration for the shifts $\sigma = 0, 3, 3.4, 6$, and check the results of the code.
2. Show that $(A - \sigma I)^{-1}$ has the same eigenvectors as $A$, with eigenvalues $1/(\lambda_i - \sigma)$. What happens if $\sigma$ is exactly an eigenvalue, and how should code handle it?
3. Starting from an eigenvector error $\varepsilon_0 = 10^{-1}$, how many RQI steps reach machine precision if the error cubes each step? How many steps of inverse iteration with ratio $0.1$?
4. Interview: You need the 5 eigenvalues of a $10^6 \times 10^6$ sparse symmetric stiffness matrix nearest a target frequency. How would you do it, and what is the expensive part? Hint: shift-invert: factor $A - \sigma I$ once with a sparse direct solver and apply Lanczos (or block inverse iteration) to $(A - \sigma I)^{-1}$; the eigenvalues nearest $\sigma$ become the largest; the factorization dominates the cost (scipy.sparse.linalg.eigsh(A, k=5, sigma=sigma)).

QR Algorithm

Why it matters. When you call np.linalg.eig or eigh on a dense matrix, all $n$ eigenvalues are computed by the QR algorithm, developed independently by Francis and by Kublanovskaya around 1961. It is backward stable, it costs $O(n^3)$, comparable to a few matrix multiplications, and it computes the Schur form of Section 8.3. Even numpy.roots is an eigenvalue computation: it finds polynomial roots as eigenvalues of the companion matrix (Section 8.1).

The basic iteration. Set $A_0 = A$ and repeat: factor $A_k = Q_kR_k$ (Section 9.3), then multiply in reverse order,

$$A_{k+1} = R_kQ_k = Q_k^TA_kQ_k.$$

Each $A_{k+1}$ is orthogonally similar to $A_k$, so all the $A_k$ have the same eigenvalues, and no conditioning is lost along the way. Under mild conditions (eigenvalues of distinct modulus), $A_k$ converges to an upper triangular matrix, the Schur form, with the eigenvalues on the diagonal in decreasing order of modulus. For symmetric $A$, the limit is diagonal. The subdiagonal entry $(i + 1, i)$ decays like $\lvert \lambda_{i+1}/\lambda_i \rvert^k$.

Why it works: simultaneous power iteration. Let $\hat Q_k = Q_1Q_2\cdots Q_k$ and $\hat R_k = R_k\cdots R_1$. By induction, $A^k = \hat Q_k\hat R_k$, so $\hat Q_k$ is the Q factor of $A^k$. Its first column is the power method applied to $\mathbf{e}_1$. Its first $j$ columns span $A^k\operatorname{span}(\mathbf{e}_1, \dots, \mathbf{e}_j)$, which is the power method applied to a whole $j$-dimensional subspace, re-orthonormalized at every step. These subspaces converge to the dominant invariant subspaces, which is exactly the Schur basis. The QR algorithm runs the power method on all subspaces at once, expressed in the coordinates of the current basis.

Making it fast: the four ingredients. The basic iteration costs $O(n^3)$ per step and converges only linearly. Practical implementations change four things:

Ingredient	What it does	Effect
Hessenberg reduction	once, by Householder reflections (Section 9.3): $A = Q_0HQ_0^T$, $H$ zero below the first subdiagonal (tridiagonal if $A$ is symmetric)	each QR step preserves the Hessenberg form and costs $O(n^2)$ (Givens rotations), or $O(n)$ for tridiagonal
shifts	factor $A_k - \mu_kI = Q_kR_k$, set $A_{k+1} = R_kQ_k + \mu_kI$, with $\mu_k$ near an eigenvalue (Wilkinson shift: the eigenvalue of the trailing $2 \times 2$ block closest to its last entry)	the last subdiagonal entry converges quadratically, or cubically for symmetric matrices (the shifted step is closely related to RQI)
deflation	once a subdiagonal entry is negligible, split the matrix and continue on the smaller blocks	the problem shrinks by one eigenvalue every few iterations
implicit double shifts (Francis)	for a real matrix with complex eigenvalues, apply a conjugate pair of shifts in real arithmetic by "chasing a bulge" down the Hessenberg matrix	real matrices stay real; the result is the real Schur form

  full A           Householder           Hessenberg H          shifted QR steps           deflate
                   (once, O(n^3))                              (O(n^2) each)

  * * * * *                              * * * * *             * * * * *                  * * * *  | *
  * * * * *                              * * * * *             * * * * *                  * * * *  | *
  * * * * *        ---------->           0 * * * *  -------->  0 * * * *     ------->     0 * * *  | *
  * * * * *                              0 0 * * *             0 0 * * *                  0 0 * *  | *
  * * * * *                              0 0 0 * *             0 0 0 e l                  ---------+---
                                                               e = negligible             0 0 0 0  | l
                                                                                          continue on the 4x4 block

With all four ingredients, about two iterations per eigenvalue suffice. The total cost is about $10n^3$ flops for all eigenvalues of a general matrix, and about $\tfrac43n^3$ for a symmetric one, most of it in the initial reduction. Eigenvectors come from accumulating the orthogonal transformations, or from inverse iteration with the computed eigenvalues. Symmetric problems also have alternative algorithms (divide and conquer, MRRR, Jacobi's method) that LAPACK uses for speed or accuracy.

import numpy as np
from scipy.linalg import hessenberg

A = np.array([[4., 1., 2., 0.],
              [1., 3., 0., 1.],
              [2., 0., 2., 1.],
              [0., 1., 1., 1.]])                # symmetric
print(np.round(np.linalg.eigvalsh(A)[::-1], 6)) # [ 5.672824  3.        1.648061 -0.320884]

# Unshifted QR: A_{k+1} = R_k Q_k = Q_k^T A_k Q_k is similar to A_k
Ak = A.copy()
for k in range(1, 61):
    Qk, Rk = np.linalg.qr(Ak)
    Ak = Rk @ Qk
    if k in (1, 10, 30, 60):
        print(k, np.round(np.diag(Ak), 6), f"{np.abs(np.tril(Ak, -1)).max():.1e}")
# 1 [ 5.47619   3.089744  1.677968 -0.243902] 6.2e-01
# 10 [ 5.672823  3.000001  1.648061 -0.320884] 1.7e-03
# 30 [ 5.672824  3.        1.648061 -0.320884] 4.9e-09    linear: slowest ratio 1.648 / 3 = 0.55
# 60 [ 5.672824  3.        1.648061 -0.320884] 2.5e-17

# Step 1 of practical QR: reduce to Hessenberg form (tridiagonal for symmetric A) once, by orthogonal similarity
H = hessenberg(A)
print(np.round(H, 4) + 0.0)
# [[ 4.     -2.2361  0.      0.    ]
#  [-2.2361  2.2    -1.4     0.    ]
#  [ 0.     -1.4     1.3918 -0.867 ]
#  [ 0.      0.     -0.867   2.4082]]

# Shifted QR with deflation (Wilkinson shift), for a symmetric matrix
def qr_eigenvalues(A, tol=1e-14):
    H = hessenberg(A)
    eigs, iters = [], 0
    while len(H) > 1:
        n = len(H)
        a, b, c = H[n-2, n-2], H[n-1, n-2], H[n-1, n-1]
        d = (a - c) / 2                       # Wilkinson shift: eigenvalue of the trailing 2x2 block nearest c
        mu = c - np.sign(d if d != 0 else 1.0) * b**2 / (abs(d) + np.hypot(d, b))
        Q, R = np.linalg.qr(H - mu * np.eye(n))
        H = R @ Q + mu * np.eye(n)
        iters += 1
        if abs(H[n-1, n-2]) < tol * (abs(H[n-1, n-1]) + abs(H[n-2, n-2])):
            eigs.append(H[n-1, n-1])          # converged: deflate, continue on the leading block
            H = H[:n-1, :n-1]
    eigs.append(H[0, 0])
    return np.sort(eigs)[::-1], iters

vals, iters = qr_eigenvalues(A)
print(np.round(vals, 6), iters)               # [ 5.672824  3.        1.648061 -0.320884] 7: 7 steps instead of 60

# numpy.roots is an eigenvalue computation: roots of x^3 - 6x^2 + 11x - 6 via the companion matrix
C = np.array([[6., -11., 6.], [1., 0., 0.], [0., 1., 0.]])
print(np.round(np.sort(np.linalg.eigvals(C).real), 10), np.round(np.sort(np.roots([1, -6, 11, -6])), 10))
# [1. 2. 3.] [1. 2. 3.]

Key takeaways

• The QR algorithm iterates $A_k = Q_kR_k$, $A_{k+1} = R_kQ_k$: a sequence of orthogonal similarities that converges to the Schur form (diagonal for symmetric matrices).
• It is simultaneous power iteration on all subspaces at once; unshifted, the subdiagonal entries decay like $\lvert \lambda_{i+1}/\lambda_i \rvert^k$.
• Practical versions reduce to Hessenberg (or tridiagonal) form first, use Wilkinson or Francis double shifts, and deflate converged eigenvalues.
• The total cost is about $10n^3$ flops (general) or $\tfrac43n^3$ (symmetric, values only), and the algorithm is backward stable; it is what eig, eigh, and roots call.

🧪 Practice

1. Perform one unshifted QR step by hand on $\begin{bmatrix} 2 & 1 \ 1 & 2 \end{bmatrix}$, and check that the result has the same trace and determinant and a smaller off-diagonal entry.
2. Prove that if $H$ is upper Hessenberg and $H = QR$, then $RQ$ is also upper Hessenberg. (Show that $Q$ is Hessenberg, using Givens rotations or $Q = HR^{-1}$ for invertible $H$.)
3. Show by induction that $A^k = (Q_1\cdots Q_k)(R_k\cdots R_1)$, and explain why this links the QR algorithm to the power method.
4. Interview: Why does every eigenvalue library first reduce the matrix to Hessenberg or tridiagonal form instead of running QR steps on the full matrix? Hint: compare the cost of one QR step on a full matrix ($O(n^3)$) with one on a Hessenberg matrix ($O(n^2)$) or a tridiagonal one ($O(n)$); the reduction costs $O(n^3)$ once, the iteration needs $O(n)$ steps, and the form is preserved by every step.

Krylov Subspaces and Arnoldi Iteration

Why it matters. The QR algorithm needs the dense matrix and $O(n^3)$ work, impossible for $n = 10^6$. Large applications (vibrations of a finite element model, stability of a fluid flow, the top principal components of a huge data matrix) need only a few eigenvalues, and can afford only matrix-vector products. The power method uses those products wastefully: it keeps only the latest vector $A^k\mathbf{b}$. Krylov methods keep all of them and extract the best possible approximations from the space they span. The same subspaces power GMRES and CG for linear systems (Section 12.2).

Krylov subspaces. For a matrix $A$ and a starting vector $\mathbf{b}$, the $k$-th Krylov subspace is

$$\mathcal{K}_k(A, \mathbf{b}) = \operatorname{span}{\mathbf{b}, A\mathbf{b}, A^2\mathbf{b}, \dots, A^{k-1}\mathbf{b}} = {p(A)\mathbf{b} : \deg p < k}.$$

Every vector in it is a polynomial in $A$ applied to $\mathbf{b}$. The power method uses the single polynomial $p(t) = t^{k-1}$. A Krylov method chooses the best polynomial of degree $< k$, and well-chosen polynomials, such as Chebyshev polynomials, amplify the extreme eigenvalues far more strongly than $t^{k-1}$. That is why Krylov methods converge much faster than the power method with the same number of matrix-vector products. Once the dimension of $\mathcal{K}_k$ stops growing, the Krylov subspace is invariant under $A$ (Section 8.3); it is then the smallest invariant subspace containing $\mathbf{b}$.

Why not use the basis $\mathbf{b}, A\mathbf{b}, A^2\mathbf{b}, \dots$ directly? These vectors all converge toward the dominant eigenvector, so they become nearly parallel. With $k = 20$, the condition number of this basis is around $10^{23}$, and it is numerically useless.

Arnoldi iteration. Build an orthonormal basis $\mathbf{q}_1, \dots, \mathbf{q}_k$ of $\mathcal{K}_k$ one vector at a time. Multiply the newest basis vector by $A$, orthogonalize the result against all previous basis vectors (with modified Gram-Schmidt, Section 9.3), and normalize:

$$\mathbf{q}1 = \mathbf{b}/\lVert \mathbf{b} \rVert; \qquad \text{for } j = 1, \dots, k: \quad \mathbf{v} = A\mathbf{q}j, \quad h{ij} = \mathbf{q}i^T\mathbf{v}, \quad \mathbf{v} \leftarrow \mathbf{v} - \sum{i \le j}h{ij}\mathbf{q}i, \quad h{j+1,j} = \lVert \mathbf{v} \rVert, \quad \mathbf{q}{j+1} = \mathbf{v}/h{j+1,j}.$$

The coefficients form a $(k + 1) \times k$ upper Hessenberg matrix $\tilde H_k$, and the Arnoldi relation holds:

$$AQ_k = Q_{k+1}\tilde H_k, \qquad Q_k^TAQ_k = H_k \ \ (\text{the top } k \times k \text{ block}).$$

$H_k$ is the matrix of $A$ compressed onto the Krylov subspace, the orthogonal projection of $A$ onto $\mathcal{K}k$. If $h{j+1,j} = 0$, the Krylov subspace is invariant and the eigenvalues of $H_j$ are exact eigenvalues of $A$, a "lucky breakdown".

Rayleigh-Ritz: eigenvalue estimates. The eigenvalues of the small matrix $H_k$ are the Ritz values. They approximate eigenvalues of $A$, the extreme and well-separated ones first. If $H_k\mathbf{y} = \theta\mathbf{y}$, then the Ritz vector $Q_k\mathbf{y}$ has residual norm $\lVert AQ_k\mathbf{y} - \theta Q_k\mathbf{y} \rVert = h_{k+1,k}\lvert y_k \rvert$, which is a free convergence test.

Costs and restarts. Step $j$ costs one matrix-vector product plus $O(nj)$ for orthogonalization, and all $k$ basis vectors must be stored. To limit memory, practical codes restart: they compress the basis to the most promising Ritz vectors and continue. ARPACK, behind scipy.sparse.linalg.eigs, uses implicitly restarted Arnoldi; Krylov-Schur is a newer variant. The which argument selects the target: "LM" (largest magnitude), "SR" (smallest real part), and so on.

Arnoldi for linear systems. GMRES (Section 12.2) takes $\mathbf{x}_k \in \mathcal{K}_k(A, \mathbf{b})$ minimizing $\lVert \mathbf{b} - A\mathbf{x}_k \rVert$. With the Arnoldi relation, this becomes a $(k + 1) \times k$ least-squares problem with $\tilde H_k$, solved with Givens rotations (Section 9.3).

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

def arnoldi(A, b, k):
    """k steps of Arnoldi: A Q_k = Q_{k+1} H_k, with Q orthonormal and H (k+1) x k upper Hessenberg."""
    n = len(b)
    Q = np.zeros((n, k + 1))
    H = np.zeros((k + 1, k))
    Q[:, 0] = b / np.linalg.norm(b)
    for j in range(k):
        v = A @ Q[:, j]                        # one matrix-vector product per step
        for i in range(j + 1):                 # modified Gram-Schmidt against all previous vectors
            H[i, j] = Q[:, i] @ v
            v -= H[i, j] * Q[:, i]
        H[j + 1, j] = np.linalg.norm(v)
        if H[j + 1, j] < 1e-12:                # invariant subspace found: Ritz values are exact
            return Q[:, :j + 1], H[:j + 1, :j + 1]
        Q[:, j + 1] = v / H[j + 1, j]
    return Q, H

rng = np.random.default_rng(0)
n = 400
# Nonsymmetric matrix: three outlying eigenvalues 12, 10, 8 and 397 others spread over [-1, 1]
lam = np.concatenate([[12., 10., 8.], rng.uniform(-1, 1, n - 3)])
X = rng.standard_normal((n, n))
A = X @ np.diag(lam) @ np.linalg.inv(X)          # similar to diag(lam), far from normal

b = rng.standard_normal(n)
for k in (5, 10, 20, 40):
    Q, H = arnoldi(A, b, k)
    ritz = np.sort(np.linalg.eigvals(H[:k, :k]).real)[::-1][:3]   # Ritz values: eigenvalues of the small H_k
    print(k, np.round(ritz, 8))
# 5 [11.94594238  8.49510055  0.87374802]
# 10 [11.99999995 10.00000029  7.99999863]  10 products: 7-8 digits for all three outliers
# 20 [12. 10.  8.]
# 40 [12. 10.  8.]

Q, H = arnoldi(A, b, 20)
print(np.allclose(A @ Q[:, :20], Q @ H), np.allclose(Q.T @ Q, np.eye(21)))   # True True: Arnoldi relation, orthonormality

# Monomial Krylov basis b, Ab, A^2 b, ... becomes numerically dependent very quickly
K = np.column_stack([np.linalg.matrix_power(A, j) @ b for j in range(20)])
print(f"{np.linalg.cond(K):.1e}", f"{np.linalg.cond(Q):.1e}")   # 1.8e+23 1.0e+00

# Library version: implicitly restarted Arnoldi (ARPACK)
print(np.round(np.sort(eigs(A, k=3, which="LM", return_eigenvectors=False).real)[::-1], 8))   # [12. 10.  8.]

Key takeaways

• The Krylov subspace $\mathcal{K}_k(A, \mathbf{b}) = \operatorname{span}{\mathbf{b}, A\mathbf{b}, \dots, A^{k-1}\mathbf{b}}$ contains every $p(A)\mathbf{b}$ with $\deg p < k$; Krylov methods pick the best such polynomial, the power method only $t^{k-1}$.
• Arnoldi builds an orthonormal basis with modified Gram-Schmidt and produces $AQ_k = Q_{k+1}\tilde H_k$ with $H_k$ upper Hessenberg.
• Ritz values (eigenvalues of $H_k$) approximate the extreme eigenvalues of $A$ after far fewer matrix-vector products than the power method; the residual $h_{k+1,k}\lvert y_k \rvert$ tests convergence for free.
• Memory and orthogonalization grow with $k$, so practical codes restart (ARPACK's eigs); the same machinery gives GMRES.

🧪 Practice

1. For $A = \begin{bmatrix} 2 & 1 & 0 \ 0 & 2 & 1 \ 0 & 0 & 3 \end{bmatrix}$, explain why Arnoldi started from $\mathbf{b} = \mathbf{e}_1$ stops after one step. Then run two steps by hand from $\mathbf{b} = \mathbf{e}_3$, write down $H_2$, and compare its Ritz values with the eigenvalues $2, 2, 3$.
2. Prove that if $A$ has only $d$ distinct eigenvalues and is diagonalizable, then $\dim\mathcal{K}_k(A, \mathbf{b}) \le d$ for every $k$, so Arnoldi breaks down after at most $d$ steps with exact eigenvalues.
3. Derive the Ritz residual formula $\lVert AQ_k\mathbf{y} - \theta Q_k\mathbf{y} \rVert = h_{k+1,k}\lvert y_k \rvert$ from the Arnoldi relation.
4. Interview: A fluid-dynamics code needs the 10 eigenvalues with largest real part of a $5 \times 10^6$ sparse Jacobian to check stability. Why not the QR algorithm, why not the power method, and what would you use? Hint: QR needs $O(n^2)$ memory and $O(n^3)$ time; the power method finds only the largest-magnitude eigenvalue, and slowly; use implicitly restarted Arnoldi (eigs(J, k=10, which="LR")), possibly with shift-invert or a Cayley transform to target the right half-plane.

Lanczos Method

Why it matters. Most large eigenproblems in practice are symmetric: vibration modes, quantum Hamiltonians, graph Laplacians, covariance and kernel matrices. For symmetric matrices, Arnoldi simplifies dramatically. Each new vector needs to be orthogonalized only against the previous two, the compressed matrix is tridiagonal, and the cost per step no longer grows with $k$. This is the Lanczos method, behind scipy.sparse.linalg.eigsh, and it is closely related to CG. Its one serious flaw, the loss of orthogonality in floating point, explains the spurious duplicate eigenvalues mentioned in Section 9.3.

From Arnoldi to Lanczos. If $A$ is symmetric, then $H_k = Q_k^TAQ_k$ is symmetric. A symmetric upper Hessenberg matrix is tridiagonal:

$$T_k = \begin{bmatrix} \alpha_1 & \beta_1 & & \ \beta_1 & \alpha_2 & \ddots & \ & \ddots & \ddots & \beta_{k-1} \ & & \beta_{k-1} & \alpha_k \end{bmatrix}.$$

So $h_{ij} = 0$ for $i < j - 1$: $A\mathbf{q}_j$ is automatically orthogonal to $\mathbf{q}1, \dots, \mathbf{q}{j-2}$. The Arnoldi step collapses to the three-term recurrence

$$\beta_j\mathbf{q}_{j+1} = A\mathbf{q}j - \alpha_j\mathbf{q}j - \beta{j-1}\mathbf{q}{j-1}, \qquad \alpha_j = \mathbf{q}_j^TA\mathbf{q}_j, \qquad \beta_j = \lVert A\mathbf{q}j - \alpha_j\mathbf{q}j - \beta{j-1}\mathbf{q}{j-1} \rVert.$$

	Arnoldi	Lanczos
Matrices	general	symmetric (Hermitian)
Compressed matrix	upper Hessenberg $H_k$	tridiagonal $T_k$
Orthogonalization per step	against all $j$ previous vectors, $O(nj)$	against the previous two, $O(n)$
Vectors needed to continue	all $k$	the last two (all $k$ only to form Ritz vectors)
Ritz values	eigenvalues of $H_k$ (complex in general)	eigenvalues of $T_k$ (real, $O(k^2)$ or less to compute)
Linear-system cousin	GMRES	CG and MINRES

Convergence. The eigenvalues of $T_k$ are Rayleigh-Ritz approximations from $\mathcal{K}k$, so by Courant-Fischer (Section 10.3) the largest Ritz value increases toward $\lambda{\max}$ and the smallest decreases toward $\lambda_{\min}$. The extreme eigenvalues converge first, and fastest when they are well separated from the rest of the spectrum (Kaniel-Paige theory, based on Chebyshev polynomials). Interior eigenvalues need shift-invert: Lanczos on $(A - \sigma I)^{-1}$ (eigsh(A, sigma=...)).

Lanczos and CG. CG applied to $A\mathbf{x} = \mathbf{b}$ is the Lanczos process started from $\mathbf{r}_0$ in disguise. The CG coefficients $\alpha_k, \beta_k$ determine the entries of $T_k$, and the CG iterate is the Krylov approximation obtained by solving with $T_k$. That is why CG converges fast exactly when the Ritz values capture the spectrum quickly.

Loss of orthogonality and ghost eigenvalues. In exact arithmetic, the three-term recurrence keeps all the $\mathbf{q}_j$ orthonormal. In floating point it does not. Paige's analysis (1970s) showed that orthogonality is lost precisely when a Ritz value converges: the new vectors regain components along the converged eigenvector, and Lanczos "finds" that eigenvalue again. The result is spurious duplicate copies of converged eigenvalues, ghosts, in $T_k$. The eigenvalues themselves remain accurate; it is their multiplicities that become wrong.

Remedy	Cost	Idea
full reorthogonalization	$O(nk)$ per step, store all $\mathbf{q}_j$	Gram-Schmidt against all previous vectors (back to Arnoldi's cost)
selective / partial reorthogonalization	much less	reorthogonalize only against converged Ritz vectors, or only when an estimate of the lost orthogonality exceeds $\sqrt u$
implicit restarts (ARPACK eigsh)	bounded memory	keep $k$ small, restart with the best Ritz vectors, with reorthogonalization inside
accept and filter ghosts (Cullum-Willoughby)	none	no reorthogonalization; identify spurious copies afterwards

import numpy as np
from scipy.linalg import eigh_tridiagonal
from scipy.sparse.linalg import eigsh

def lanczos(A, b, k, reorth=False):
    """k steps of Lanczos for symmetric A: returns the tridiagonal T_k (alpha on, beta off the diagonal)."""
    n = len(b)
    Q = np.zeros((n, k + 1))
    alpha, beta = np.zeros(k), np.zeros(k)
    Q[:, 0] = b / np.linalg.norm(b)
    for j in range(k):
        v = A @ Q[:, j]
        alpha[j] = Q[:, j] @ v
        v -= alpha[j] * Q[:, j] + (beta[j - 1] * Q[:, j - 1] if j > 0 else 0)   # three-term recurrence
        if reorth:                                   # full reorthogonalization (costly, but robust)
            v -= Q[:, :j + 1] @ (Q[:, :j + 1].T @ v)
        beta[j] = np.linalg.norm(v)
        Q[:, j + 1] = v / beta[j]
    return alpha, beta[:-1], Q[:, :k]

n = 1000
lam = np.concatenate([np.linspace(0, 1, n - 1), [2.0]])    # one well-separated top eigenvalue
A = np.diag(lam)                                           # diagonal: eigenvalues known exactly
b = np.ones(n)

for k in (20, 60, 120):
    for reorth in (False, True):
        alpha, beta, Q = lanczos(A, b, k, reorth)
        ritz = eigh_tridiagonal(alpha, beta, eigvals_only=True)
        copies = int(np.sum(np.abs(ritz - 2.0) < 1e-8))    # Ritz values that have converged to 2
        loss = np.abs(Q.T @ Q - np.eye(k)).max()           # departure from orthonormality
        print(k, reorth, copies, f"{loss:.1e}")
# 20 False 1 5.4e-05                         the top eigenvalue has converged; orthogonality starts to go
# 20 True 1 6.7e-16
# 60 False 3 6.9e-01                         orthogonality lost: two ghost copies of the eigenvalue 2
# 60 True 1 8.9e-16
# 120 False 6 6.9e-01                        five ghosts: a new copy appears every ~20 steps
# 120 True 1 8.9e-16

# Library version: implicitly restarted Lanczos (ARPACK)
print(np.round(eigsh(A, k=3, which="LA", return_eigenvectors=False), 8))   # [0.998998 1.       2.      ]

Key takeaways

• For symmetric $A$, Arnoldi becomes Lanczos: a three-term recurrence that produces a tridiagonal $T_k = Q_k^TAQ_k$ at $O(n)$ extra cost per step.
• Ritz values (eigenvalues of $T_k$) converge first to the extreme, well-separated eigenvalues; shift-invert targets interior ones.
• CG is Lanczos applied to a linear system.
• In floating point, orthogonality is lost as Ritz values converge, which creates ghost copies of converged eigenvalues; reorthogonalization (full, selective, or via implicit restarts as in eigsh) prevents them.

🧪 Practice

1. Show that if $A$ is symmetric, the Arnoldi coefficients satisfy $h_{ij} = 0$ for $i < j - 1$ and $h_{j-1,j} = h_{j,j-1}$, so that $H_k$ is symmetric tridiagonal.
2. Run two Lanczos steps by hand for $A = \operatorname{diag}(1, 2, 3)$ and $\mathbf{b} = (1, 1, 1)$. Compute $T_2$ and its eigenvalues, and compare them with the extreme eigenvalues $1$ and $3$.
3. Using the code, plot (or print) the smallest singular value of $Q_k$ against $k$ without reorthogonalization, and check that orthogonality is lost right after the top Ritz value converges.
4. Interview: You inherit an Arnoldi or Lanczos implementation whose computed eigenvalues show spurious duplicate copies after a few hundred iterations (Section 9.3). Explain the cause, how you would confirm it, and two ways to fix it. Hint: Paige: orthogonality is lost exactly as Ritz values converge; check $\lVert Q_k^TQ_k - I \rVert$; fix with full or selective reorthogonalization (using modified Gram-Schmidt or a second pass), or switch to an implicitly restarted library routine such as ARPACK's eigsh.

<a id="13-advanced-topics-and-applications"></a>