Iterative Methods for Large Linear Systems
Iterative methods build successive approximations to solutions of large sparse linear systems. Stationary iterations converge from every start when their iteration matrix has spectral radius below 1. Suitable preconditioning can reduce the work needed.
What a learner can do afterwards
- Run a splitting step, calculate a residual to check a guess, and state the iteration's convergence condition
- Say why the conjugate gradient method needs a symmetric positive definite matrix
- Explain what a preconditioner changes and why that helps
1 · Read
You have seen how LU factors can serve repeated solves. For a large sparse matrix, factorisation may require too much storage or work, but it can still be suitable. An iterative method builds successive guesses. Check a guess x with its residual r = b - Ax: multiply by A, then subtract from b. For an ill-conditioned matrix, a small residual need not mean a small solution error.
A splitting writes A = B - C, where B is invertible and cheap to solve. Repeat B x_next = C x_old + b; the iteration matrix is T = B^-1 C. Jacobi takes B as the diagonal of A, so those diagonal entries must be nonzero. It computes every new coordinate from old values. Gauss-Seidel instead uses each fresh coordinate as soon as it is available.
For 5x + y = 11 and x + 4y = 9, Jacobi uses x_next = (11 - y_old)/5 and y_next = (9 - x_old)/4. From (0,0), the first guess is (2.2,2.25). The next guess is (1.75,1.7). Both second-step updates use the first guess: y_next uses x_old = 2.2, not the fresh x value 1.75.
A stationary iteration fixes T and c: x_next = T x_old + c. If the solution x_star is a fixed point, subtract its equation: e = x - x_star obeys e_next = T e_old. An eigenvalue lambda satisfies T v = lambda v for nonzero v; lambda and v may be complex. The spectral radius rho(T) is the largest eigenvalue magnitude. Convergence to x_star from every start holds exactly when rho(T) < 1. Error need not decrease at each step. Starting at x_star stays there even if rho(T) >= 1.
Strict row diagonal dominance means that, in each row, the absolute diagonal entry exceeds the sum of the absolute off-diagonal entries. It guarantees convergence of Jacobi and Gauss-Seidel from every start. It is sufficient, not necessary: without it, Jacobi may still converge. In the worked system, 5 > 1 and 4 > 1, so the guarantee applies.
A real matrix is SPD (symmetric positive definite) if A = A^T and v^T A v > 0 for every nonzero v. For symmetric A, this is equivalent to all eigenvalues being positive. A diagonal matrix is symmetric; its diagonal entries are its eigenvalues. One nonzero v with v^T A v < 0 rules out positive definiteness; positive diagonal entries alone do not prove it. The span of vectors is the set of all their linear combinations. The conjugate gradient method is abbreviated CG.
A = diag(3,6) is symmetric. Its quadratic form is positive for every nonzero (a,b), so it is SPD.
For symmetric A, f(x) = (1/2)x^T A x - b^T x has gradient Ax - b, the vector of first partial derivatives. SPD makes f strictly convex: strictly below each chord between its endpoints. Its unique minimum solves Ax = b. Standard CG needs both assumptions. From start x0, its kth iterate minimizes f over nested affine sets x0 + K_k in exact arithmetic. For k >= 1, K_k is the span of r0, Ar0, ..., A^(k-1)r0, where r0 = b - Ax0. CG need not be faster than another solver.
Choose an invertible preconditioner M whose systems are cheap to solve. Left preconditioning gives M^-1 A x = M^-1 b. Transform both sides to preserve the exact solution. Apply M^-1 by solving with M, rather than forming an inverse. For M = diag(2,5) and b = (8,15), solve M c = b to get c = (4,3). A diagonal solve divides each component by its diagonal entry.
For CG, use a symmetric positive definite M through preconditioned CG. Do not apply ordinary CG blindly to M^-1 A: it need not be symmetric in the usual inner product. A smaller ratio of largest to smallest preconditioned eigenvalue can improve the convergence bound; clustering can also help. Neither improvement nor a faster total solve is automatic. Solving with M also has a cost.
Jacobi uses old values. Stationary convergence from every start needs rho(T) < 1. Standard CG needs SPD assumptions. Judge preconditioning by its effect and cost.
The lesson, practice questions and flashcards were written for this dataset.
The readings, practice questions and flashcards were written for this dataset. The linked video Introduction to eigenvalues and eigenvectors is listed as a Khan Academy resource. Its reuse terms have not been verified.
More to read
Comparing the whole solve
A sparse matrix has mostly zero entries. Those zeros can save storage and arithmetic, but elimination can create new nonzero entries. A factorisation may therefore need much more storage than the original matrix.
Iteration offers another approach. Keep a guess, calculate an update and repeat. A cheap update is useful only if enough updates reach the required accuracy. A method that diverges is not rescued by a low cost per step.
Compare the full calculation. A direct factorisation can be reused for several right-hand sides. An iterative solver may avoid that factorisation but require many matrix-vector products and solves with a preconditioner. Storage limits, convergence and requested accuracy all affect the choice.
Both approaches use numerical arithmetic in a computer. Calling a method direct does not mean its computed answer is free of rounding error.
The value used in the next equation
Consider 3x + y = 8 and x + 2y = 5, starting at (1,1). The rearranged equations are x = (8 - y)/3 and y = (5 - x)/2.
Jacobi keeps the old pair while calculating both new coordinates. It gives x_new = (8 - 1)/3 = 7/3 and y_new = (5 - 1)/2 = 2. The new pair is (7/3,2).
Gauss-Seidel, with x updated first, also begins with x_new = 7/3. Its next calculation uses that fresh value: y_new = (5 - 7/3)/2 = 4/3. Its new pair is (7/3,4/3).
The methods used the same equations and starting guess. Their results differ because the y calculation used different x values. Name the method and update order before checking an iteration.
Finding the solution through a minimum
For A = diag(2,5) and b = (4,5), solving the two diagonal equations gives x = (2,1). The energy is f(x1,x2) = x1^2 + (5/2)x2^2 - 4x1 - 5x2.
Writing this as f = (x1 - 2)^2 + (5/2)(x2 - 1)^2 - 13/2 shows its minimum. Both squared terms are nonnegative and vanish together only at (2,1). The unique energy minimum is therefore the solution of Ax = b.
For any real symmetric positive definite A, use f(x) = (1/2)x^T A x - b^T x. Symmetry gives gradient Ax - b. Positive definiteness makes f strictly convex, so the equation Ax - b = 0 identifies its unique minimum.
The CG search domain depends on the starting guess x0. Define the initial residual r0 = b - Ax0. For k >= 1, let K_k be the span of r0, Ar0, ..., A^(k-1)r0. In exact arithmetic, the kth CG iterate minimizes f over x0 + K_k: all points x0 + v with v in K_k. These affine sets are nested. They are linear subspaces when x0 = 0, but need not be linear for other starts.
The standard CG guarantee requires both symmetry and positive definiteness. Merely writing down a quadratic expression does not establish that guarantee for an arbitrary matrix.
A small residual can hide an error
Take A = diag(1,0.001) and b = (1,0.001). The exact solution is (1,1). Suppose an iterative solver stops at (1,0).
At this guess, Ax = (1,0). The residual b - Ax is therefore (0,0.001). Its components are small, but the second coordinate of the guess differs from the solution by 1.
For an invertible matrix, x_star - x_guess = A^-1 r. Here A^-1 multiplies the second residual component by 1000. The small residual and larger solution error are consistent with this relation.
A residual is useful because it can be calculated without knowing x_star. Its interpretation still depends on scaling and conditioning. A sensitive system can amplify the remaining residual; that is different from an algorithm introducing excessive rounding error.
Following the error instead of the guess
The scalar update z_next = -0.25 z_old + 5 has fixed point z_star = 4. Substitution checks it: -0.25 times 4 + 5 = 4.
Subtract this fixed-point equation from the update. If e = z - 4, then e_next = -0.25 e_old. Each update changes the sign of a nonzero error and divides its magnitude by four. Every starting error tends to zero.
For a matrix iteration, the corresponding rule is e_next = T e_old. The largest eigenvalue magnitude is the spectral radius rho(T). Convergence to the fixed-point solution from every start holds exactly when rho(T) < 1. Unlike the scalar example, a general matrix can temporarily increase the error norm before it eventually decreases.
Starting exactly at a fixed point leaves the guess there, even for an iteration that fails from other starts. This is why the phrase from every start matters.
Explore more
2 · Watch
Take it off screen
Where it sits
Learn first
Optional connections
This opens up
Nothing builds on it yet.
14 questions wait behind this lesson, each with its answer explained. Every answer feeds the sky: stars light as they are learned, and dim when it is time to come back.