2NM4 Week 3: Linear systems: Iterative

Overview

The lecture compares stationary iterative solvers, Krylov methods, and preconditioners for solving linear systems, explaining when each is useful and how convergence, memory, and sparsity affect performance. It derives Jacobi, Gauss-Seidel, SOR, conjugate gradient, and GMRES at a conceptual level, then introduces incomplete LU, domain decomposition, multigrid, and the Moore-Penrose pseudoinverse for rectangular systems.

Key takeaways

Dense LU has O(n^3) computational cost and can become infeasible in memory, while iterative methods are most valuable for large sparse systems where their convergence can be accelerated.
Jacobi is easy to parallelize because every update uses xₖ; Gauss-Seidel often converges faster and needs less memory because it reuses new values, but its updates depend on sequence.
Conjugate gradient is the preferred low-memory Krylov method for symmetric positive-definite matrices because it preserves A-conjugacy using only a few vectors.
GMRES supports nonsymmetric systems and typically reduces the residual monotonically within a cycle, but its stored vector basis grows with iteration count, making restart and preconditioning important.
Incomplete LU controls fill-in with a drop tolerance; the lecture's example reduced stored nonzeros from 627 to 284 while improving GMRES convergence.
The Moore-Penrose pseudoinverse generalizes matrix inversion: it returns a least-squares solution for inconsistent overdetermined systems and a minimum-norm solution for underdetermined systems.

Chapters

0:00 Why Use Iterative Solvers Instead of Dense LU?

Dense LU typically costs O(n^3), so doubling problem size can multiply time and memory needs by about eight.
Iterative methods often use less memory and can be faster on large problems, but their stability and convergence rates are less predictable.
Sparse matrices and matrix-free methods can make iterative approaches practical when storing or factoring A is too costly.

4:15 Choosing Residual and Solution Tolerances

For Ax = b, define the residual as r = Ax - b and stop when a norm such as the Euclidean or L-infinity norm falls below tolerance.
An absolute residual tolerance can accept a poor answer when the solution scale is very small; a relative tolerance can lose sensitivity near convergence.
A solution-change tolerance can detect stagnation rather than a correct answer, so residual and solution criteria must be interpreted carefully.

8:40 Stationary Methods Update the Solution Element by Element

Stationary methods apply the same update rule at each iteration and revise individual components of the solution vector.
Each component is adjusted to make its equation consistent with the current values of the other components.
The iteration continues until the residual satisfies a chosen tolerance.

11:00 Jacobi Iteration from the LDU Splitting

Write each row of Ax = b in index notation, isolate the diagonal term aᵢᵢxᵢ, and solve for the next value of xᵢ.
The matrix form separates A into its diagonal D and off-diagonal parts, yielding an update based on D inverse.
Jacobi computes every component of xₖ₊₁ from the same old vector xₖ, so it requires separate current and next vectors.

14:40 Jacobi Convergence, Parallelism, and a 2×2 Example

Jacobi is guaranteed to converge when A is diagonally dominant; row permutations can sometimes restore that property.
A 2×2 demonstration takes about 10 to 12 iterations, illustrating why direct methods are preferable for tiny systems.
Because components use only the previous iterate, Jacobi updates can be split across processors or executed in parallel.

18:30 Gauss-Seidel Reuses Newly Updated Components

Gauss-Seidel inserts each newly computed component into the solution immediately and uses it to calculate later components in the same iteration.
In-place updates reduce memory use and often converge faster than Jacobi; diagonal dominance or symmetric positive definiteness guarantees convergence.
The dependence on newly updated values makes straightforward parallelization harder than for Jacobi.

22:59 SOR Controls Gauss-Seidel Step Size with Omega

Successive over-relaxation combines the current estimate with the Gauss-Seidel update using a relaxation parameter ω.
ω = 1 recovers Gauss-Seidel, 0 < ω < 1 damps updates, and ω > 1 amplifies them.
A poorly chosen ω can overshoot or destabilize convergence, so selecting it usually requires empirical or heuristic tuning.

27:42 When Stationary Methods Outrun LU in Benchmarks

Hand-coded Python Jacobi and Gauss-Seidel were compared with optimized SciPy LU on randomly generated diagonally dominant matrices.
LU was substantially faster for small cases, while Gauss-Seidel overtook it around a 3,000-by-3,000 example and used less memory.
The comparison is affected by Python loop overhead and does not account for LU's larger memory requirements.

32:15 Jacobi and Gauss-Seidel Trade Parallelism for Faster Updates

Jacobi computes a whole new vector from the previous vector, while Gauss-Seidel uses updated values as soon as they are available.
That in-place dependency helps Gauss-Seidel converge faster but limits splitting one iteration across CPUs or GPU threads.
The methods have different hardware advantages, and neither is universally superior.

38:21 Krylov Subspaces Search Along Constructed Directions

Krylov methods build a basis of search directions instead of updating each solution component independently.
An n-dimensional system could require up to n independent directions in exact arithmetic, although practical convergence may occur earlier.
Basis storage and roundoff accumulation are key numerical costs, making preconditioning important for difficult systems.

41:12 Conjugate Gradient Minimizes a Quadratic for SPD Matrices

For symmetric positive-definite A, solving Ax = b is equivalent to finding the minimum of ½xᵀAx - bᵀx.
Symmetry lets A act as the Hessian of a scalar quadratic, while positive definiteness makes that quadratic a bowl with a unique minimum.
The minimum occurs where the gradient is zero, which recovers the original linear system.

45:25 Steepest Descent Finds the Line-Minimizing Step Length

An iteration can be written as xₖ₊₁ = xₖ + αₖpₖ, where pₖ is a direction and αₖ is its step length.
For a quadratic objective, differentiating along a chosen direction gives the step length that minimizes the objective on that line.
Steepest descent chooses the negative gradient, related here to the residual, but can zigzag and take many iterations.

50:13 A-Conjugate Directions Prevent CG from Backtracking

Conjugate gradient selects directions that are orthogonal under the A-weighted inner product, rather than ordinary Euclidean orthogonality.
A-conjugate steps avoid undoing earlier progress; the illustrated 2D SPD example reaches its solution in two iterations.
The demo also checks that transformed step directions have dot products near machine zero, around 10⁻¹⁶.

55:28 CG Residuals, Stopping Criteria, and Roundoff

SciPy's conjugate-gradient call accepts A, b, an initial guess, and absolute and relative tolerances.
Residuals can fall sharply after gradual progress, and the solution components need not move monotonically toward their final values.
A seven-dimensional example takes eight iterations, showing that finite-precision roundoff can require more than n iterations.

59:00 GMRES Restarting and Matrix-Free Products

GMRES builds and stores a history of orthogonal search vectors, so memory use grows with the number of iterations.
A restart limit, often around 30 in common practice, discards older vectors to control memory and accumulated roundoff, at the cost of forgetting prior directions.
Krylov methods need matrix-vector products such as Av rather than the entries of A themselves, enabling matrix-free computation when Av is available directly.

1:03:00 GMRES Handles Nonsymmetric Systems

Conjugate gradient requires a symmetric positive-definite matrix; nonsymmetric systems generally call for GMRES.
GMRES builds ordinary orthogonal directions and tracks previous directions to keep later ones orthogonal.
SciPy's optimized GMRES does not necessarily expose each internal solution vector, but its residual history can still show convergence.

1:08:38 CG and GMRES Differ in Residuals, Memory, and Robustness

GMRES is designed to minimize the residual over its search space, so its residual typically decreases monotonically within a cycle.
CG may show a temporarily increasing Euclidean residual because it minimizes an energy norm instead.
GMRES can require a large vector history, while CG maintains conjugacy using only a few vectors, making CG far more memory-efficient for SPD systems.

1:11:57 Why CG Needs Few Vectors and GMRES Needs a Restart

CG's recurrence preserves conjugacy by using the recent search directions, avoiding storage of the full iteration history.
GMRES must orthogonalize each new direction against prior directions, increasing both memory use and sensitivity to roundoff.
Restarting GMRES reduces storage but can slow progress by discarding useful directions.

1:20:39 Preconditioning Improves Conditioning Before Iteration

Iterative convergence is sensitive to the condition number: poorly conditioned systems can converge extremely slowly or fail to converge usefully.
A preconditioner P approximates A so that applying P inverse produces an equivalent system closer to the identity.
Left, right, and two-sided preconditioning place P inverse differently; two-sided forms can preserve symmetry needed by CG.

1:26:26 Jacobi Preconditioning Scales the Diagonal

The Jacobi preconditioner uses the diagonal of A, whose inverse is inexpensive to apply and requires little storage.
A 2D Poisson heat-transfer example produces a banded sparse matrix connecting each grid point to nearby horizontal and vertical points.
The diagonal preconditioner improves GMRES convergence in the example, but the gain is modest because it mainly rescales the system.

1:30:07 Incomplete LU Trades Fill-In for a Useful Preconditioner

Full LU can create fill-in inside the original matrix bandwidth, increasing storage even when the input matrix is sparse.
Incomplete LU drops selected small entries using a drop tolerance, retaining sparse L and U factors that are cheaper to store and apply.
In the demonstrated matrix, ILU reduced nonzeros from 627 to 284; changing the drop tolerance altered fill-in and convergence.

1:34:12 Domain Decomposition, Multigrid, and Block Preconditioning

Additive Schwarz divides a geometry into overlapping subdomains, with overlap regions enabling communication between neighboring solves.
Multigrid moves between fine and coarse meshes: coarse levels address long-wavelength error, while fine levels resolve local detail.
Block preconditioning applies physics-specific strategies, such as fluid and elasticity solvers, while accounting for their coupling through block structure.

1:40:02 Solver Selection: LU, GMRES, CG, and Schwarz

Use LU when the matrix fits in memory; it is robust and usually the first choice for small or moderate systems.
Use CG for symmetric positive-definite matrices and GMRES as a general option when LU is impractical.
For large problems, established physics-specific preconditioners, additive Schwarz, or geometric multigrid can improve scalability; the lecture cites roughly 10,000 degrees of freedom per Schwarz subproblem as a rule of thumb.

1:48:53 Moore-Penrose Pseudoinverse and Least-Squares Systems

For an overdetermined m-by-n system with m > n, an exact solution may not exist, so the goal becomes minimizing the residual.
Differentiating the squared residual gives the normal equations AᵀAx = Aᵀb and the full-rank formula A⁺ = (AᵀA)⁻¹Aᵀ.
Although the normal equations explain the least-squares solution, AᵀA can be poorly conditioned; singular value decomposition is a more stable way to compute the pseudoinverse.

1:55:09 Overdetermined Fits and Minimum-Norm Underdetermined Solutions

For inconsistent overdetermined systems, the pseudoinverse finds the least-squares solution that balances residuals across equations, as illustrated by three lines in two variables.
When equations are linearly dependent, the pseudoinverse can still select a defined least-squares solution among redundant constraints.
For underdetermined systems with more unknowns than equations, it returns a minimum-norm solution, the solution closest to the origin among those that fit.

2:01:00 Pseudoinverse Generalizes the Inverse and Multigrid Recap

For a full-rank square matrix, the pseudoinverse reduces to the ordinary inverse; for rectangular systems it extends solving to least-squares or minimum-norm answers.
Multigrid can use coarse meshes to correct long-wavelength error and fine meshes for high-frequency detail, with V, W, and F cycle patterns.
The pseudoinverse is useful for curve fitting, but its output is a best-fit or minimum-norm answer rather than necessarily an exact solution.

Transcript excerpt

0:00 So, um, last time we spoke about direct methods for solving linear systems. Um, if you remember one of the problems with that was that um, it scaled really badly. So, it went up typically with dense LU solves, it would take order n cubed. Right? So, you double it, your time and memory cost roughly uh, increases by a factor of eight. Um, there were some mitigations about like you could put in sparsity, you could do uh, some of the other things. Um, that would all work okay if you activated the sparse solver and then typically you could bring it down more in lines of like n uh, order n versus order n squared. Um, so that was a big benefit, but it's still like not fantastic, right? So, whereas direct solvers were very robust and in infinite precision they would be able to find the exact answer in a fixed number of steps.

0:48 Um, iterative methods are a lot less so. Typically they they can be Typically they are much less memory intensive. Um, every once in a while they'll also be faster, but um, they there are limitations in terms of their stability. Uh, generally speaking direct solvers are much more stable than iterative solvers. Okay. Uh, we don't actually need to execute that one. So, uh, I just said the first part there. And the second one there too. Yeah. And um, even even if the iterative solver does converge, sometimes the rate of convergence is is problematic and it just like you're you're just kind of left staring at a screen that's slowly taking off the iterations without making actual progress. So, um, some of the iterative solvers typically they have they they have to have in

1:35 order to be useful, they have like a

This is the start of the transcript. Save this video free for the whole transcript, searchable and time-coded.

Related notes

Summary, takeaways, and chapters were generated by AI from the video's transcript and may contain errors. The video belongs to its creator, Numerical Methods for Engineering.