Alternative to direct solvers for large, sparse systems where matrix factorization is impractical due to memory or compute limits.
Direct methods factor the matrix. Iterative methods repeatedly improve an approximate solution using cheap operations, usually matrix-vector products. They are preferred when \(\bA\) is large, sparse, or available only as a function that computes \(\bA\bv\).
Direct vs. Iterative Solvers
(Sparse fill-in) Even if \(\bA\) is sparse, its LU or Cholesky factors are typically much denser (fill-in). For a 2D Poisson problem on an \(n \times n\) grid:
Iterative methods exploit sparsity by using only matrix-vector products (\(O(\text{nnz})\)).
(Matrix-free viewpoint) Many scientific codes never assemble \(\bA\) explicitly. They only provide a routine that maps \(\bv\mapsto\bA\bv\). Krylov methods such as CG and GMRES are designed for this setting.
Given a matrix \(\bA\) and a vector \(\bv\), the \(k\)-th Krylov subspace is \[
\begin{align}
\mathcal{K}_k(\bA,\bv)
=
\operatorname{span}\{\bv,\bA\bv,\bA^2\bv,...,\bA^{k-1}\bv\}.
\end{align}
\] It is the smallest subspace that contains \(\bv\) and is built by repeated multiplication by \(\bA\).
Let \(\br^{(0)}=\bb-\bA\bx^{(0)}\). Any iterate of the form \[
\begin{align}
\bx^{(k)}=\bx^{(0)}+p_{k-1}(\bA)\br^{(0)},
\end{align}
\] where \(p_{k-1}\) is a polynomial of degree at most \(k-1\), lies in the affine space \[
\begin{align}
\bx^{(0)}+\mathcal{K}_k(\bA,\br^{(0)}).
\end{align}
\] Thus a matrix-free method can improve \(\bx^{(0)}\) by choosing increasingly good corrections from Krylov subspaces using only matrix-vector products.
Write \(p_{k-1}(t)=c_0+c_1t+\cdots+c_{k-1}t^{k-1}\) and expand \(p_{k-1}(\bA)\br^{(0)}\).
Use the result above to show that this correction belongs to \(\mathcal{K}_k(\bA,\br^{(0)})\).
List \(\mathcal{K}_1\), \(\mathcal{K}_2\), and \(\mathcal{K}_3\) explicitly. What new vector must be computed to move from \(\mathcal{K}_k\) to \(\mathcal{K}_{k+1}\)?
Explain why the result above makes Krylov methods compatible with the matrix-free viewpoint above.
Preview the later chapters: CG chooses the best correction in an energy norm for SPD systems, while GMRES chooses the correction with smallest residual norm.
An iteration of the form , where:
\(\bT = \bI - \bM^{-1}\bA\) is the iteration matrix.
\(\bc = \bM^{-1}\bb\) is the updated right-hand side.
\(\bM\) is the splitting matrix (ideally \(\bM \approx \bA\) and \(\bM^{-1}\) is cheap).
Stationary iterations are fixed-point iterations. Instead of solving \(\bA\bx=\bb\) directly, we rewrite the problem as \(\bx=\bT\bx+\bc\) and repeatedly apply the map. The exact solution is a fixed point of this map.
(Splitting form) If \(\bA=\bM-\mathbf{N}\), then \[
\begin{align}
\bA\bx=\bb
\quad \Longleftrightarrow \quad
\bM\bx=\mathbf{N}\bx+\bb.
\end{align}
\] The iteration is \[
\begin{align}
\bx^{(k+1)}=\bM^{-1}\mathbf{N}\bx^{(k)}+\bM^{-1}\bb.
\end{align}
\] The art is choosing \(\bM\) easy to invert while making the iteration converge quickly.
Convergence and Spectral Radius
\(\rho(\bB) = \max \{|\lambda| : \lambda \in \sigma(\bB)\}\).
(Spectral radius intuition) The eigenvalue of largest magnitude controls the long-run behavior of repeated multiplication. If every eigenvalue of \(\bT\) lies inside the unit disk, powers of \(\bT\) decay and the iteration converges.
For any submultiplicative norm, \(\rho(\bB) = \lim_{k\to\infty} \|\bB^k\|^{1/k}\). For large \(k\), \(\|\bB^k\| \approx \rho(\bB)^k\).
The iteration \(\bx^{(k+1)} = \bT\bx^{(k)} + \bc\) converges for any \(\bx^{(0)}\) iff \(\rho(\bT) < 1\).
Let \(\bx^*\) be the exact solution satisfying \(\bx^* = \bT\bx^* + \bc\). Subtracting this from the iteration formula gives the error evolution: \[
\begin{align}
\bx^{(k+1)} - \bx^* &= (\bT\bx^{(k)} + \bc) - (\bT\bx^* + \bc) \\
\be^{(k+1)} &= \bT \be^{(k)}.
\end{align}
\] Applying this relation recursively from \(k=0\): \[
\begin{align}
\be^{(k)} = \bT^k \be^{(0)}.
\end{align}
\] The error \(\be^{(k)} \to \bzero\) for any \(\be^{(0)}\) iff \(\bT^k \to \bzero\) as \(k \to \infty\), which is guaranteed iff \(\rho(\bT) < 1\) by Gelfand’s formula.
(Iteration count) To reduce error by \(10^{-p}\), the required iterations are \(k \approx \frac{p}{\log_{10}(1/\rho(\bT))}\). As \(\rho(\bT) \to 1^-\), the cost explodes.
Stopping Criteria
Iterative methods need a stopping rule. The true error \(\|\bx^{(k)}-\bx^*\|\) is usually unknown, so algorithms monitor the residual \(\br^{(k)}=\bb-\bA\bx^{(k)}\).
(Common stopping test) Stop when \[
\begin{align}
\frac{\|\br^{(k)}\|}{\|\bb\|} \leq `tol`.
\end{align}
\] This is a relative residual test. It is cheap and scale-aware, but it is not the same as a relative error test unless the problem is well-conditioned.
(Preconditioning preview) Most useful iterative solvers are preconditioned. A preconditioner \(\bM\approx\bA\) should be cheap to solve with and should make the transformed system easier for the iterative method. Good preconditioning changes the problem from “possible but slow” to practical.
For \(\bA = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}\), find \(\rho(\bT)\) for the splitting \(\bM = 3\bI\).
If \(\rho(\bT) = 0.99\), how many iterations are needed to gain 4 digits of precision?
Start \(\bx^{(0)}=\bzero\) and solve \(\bA\bx = (4, 4)^T\) via 3 steps of this iteration.
Fill-in and Sparse Storage
Fill-in Demo: Generate a large sparse tridiagonal matrix. Compare nonzeros in \(\bA\) vs. LU factors.
Timing: Use scipy.sparse.linalg.spsolve vs. dense np.linalg.solve for \(n=2000\). Record memory and time.