Direct methods such as Gaussian elimination and Cholesky factorization compute exact solutions up to machine roundoff in a finite number of arithmetic operations. However, for the very large linear systems arising from the discretization of three-dimensional partial differential equations or network models, direct factorizations become computationally prohibitive. Direct factorizations suffer from fill-in, in which sparse matrices produce dense triangular factors that exceed available computer memory.
Iterative methods provide an essential alternative. Instead of transforming the matrix through elimination, an iterative solver generates a sequence of improving approximations \(\bx^{(0)}, \bx^{(1)}, \bx^{(2)}, ...\) that converges toward the true solution \(\bx^* = \bA^{-1}\bb\). The fundamental computational primitive of an iterative method is the matrix-vector multiplication \(\bv \mapsto \bA\bv\). Because sparse matrix-vector products require only \(O(\text{nnz})\) operations and zero additional memory, iterative methods scale efficiently to millions of degrees of freedom.
Matrix Splittings and Stationary Iterations
The classical approach to designing an iterative solver is matrix splitting. We decompose the coefficient matrix \(\bA\) into the difference of two matrices, \[
\begin{align}
\bA = \bM - \mathbf{N},
\end{align}
\] where the splitting matrix \(\bM\) is nonsingular and easily invertible. The linear system \(\bA\bx = \bb\) is then algebraically equivalent to \(\bM\bx = \mathbf{N}\bx + \bb\).
Given a splitting \(\bA = \bM - \mathbf{N}\), the corresponding stationary iteration computes successive approximations via the fixed-point recurrence \[
\begin{align}
\bM\bx^{(k+1)} = \mathbf{N}\bx^{(k)} + \bb \iff \bx^{(k+1)} = \bT\bx^{(k)} + \bc,
\end{align}
\] where \(\bT = \bM^{-1}\mathbf{N} = \bI - \bM^{-1}\bA\) is the iteration matrix and \(\bc = \bM^{-1}\bb\) is the modified right-hand side vector.
(Fixed-point interpretation) A stationary iteration is a linear fixed-point mapping. The exact solution \(\bx^*\) satisfies \(\bx^* = \bT\bx^* + \bc\). The effectiveness of the method depends entirely on the choice of \(\bM\): the matrix \(\bM\) should be as close to \(\bA\) as possible to minimize the magnitude of \(\bT\), while remaining computationally cheap to invert at each step.
Convergence Theory and the Spectral Radius
The convergence of a linear stationary iteration is determined entirely by the spectral properties of the iteration matrix \(\bT\).
Let \(\bT \in \fR^{n \times n}\) with spectrum \(\sigma(\bT) \subset \fC\). The spectral radius of \(\bT\) is defined as the maximum modulus among all its eigenvalues: \[
\begin{align}
\rho(\bT) = \max \{ |\lambda| : \lambda \in \sigma(\bT) \}.
\end{align}
\]
For any matrix norm \(\|\cdot\|\) induced by a vector norm on \(\fC^n\) and any matrix \(\bT \in \fR^{n \times n}\), \[
\begin{align}
\rho(\bT) = \lim_{k \to \infty} \|\bT^k\|^{1/k}.
\end{align}
\] In particular, for any \(\varepsilon > 0\), there exists an induced matrix norm \(\|\cdot\|_*\) such that \(\|\bT\|_* \leq \rho(\bT) + \varepsilon\).
The stationary iteration \(\bx^{(k+1)} = \bT\bx^{(k)} + \bc\) converges to the unique solution \(\bx^* = \bA^{-1}\bb\) for every initial guess \(\bx^{(0)} \in \fR^n\) if and only if \[
\begin{align}
\rho(\bT) < 1.
\end{align}
\]
Let \(\be^{(k)} = \bx^{(k)} - \bx^*\) denote the error vector at step \(k\). Subtracting the fixed-point identity \(\bx^* = \bT\bx^* + \bc\) from the iteration formula \(\bx^{(k+1)} = \bT\bx^{(k)} + \bc\) yields the homogeneous linear error recurrence \[
\begin{align}
\be^{(k+1)} = \bT \be^{(k)}.
\end{align}
\] Applying this recurrence recursively from the initial error \(\be^{(0)}\) gives \[
\begin{align}
\be^{(k)} = \bT^k \be^{(0)}.
\end{align}
\] The error \(\be^{(k)}\) converges to the zero vector for every choice of \(\be^{(0)}\) if and only if the matrix sequence \(\bT^k \to \bzero\) as \(k \to \infty\). By the Jordan canonical form, every Jordan block \(J_i = \lambda_i \bI + N_i\) satisfies \(J_i^k \to \bzero\) if and only if \(|\lambda_i| < 1\). Therefore, \(\lim_{k \to \infty} \bT^k = \bzero\) if and only if every eigenvalue satisfies \(|\lambda_i| < 1\), which is precisely the condition \(\rho(\bT) < 1\).
(Asymptotic convergence rate and iteration count) The quantity \(-\log_{10} \rho(\bT)\) represents the asymptotic rate of convergence, measuring the number of decimal digits of accuracy gained per iteration. To reduce the initial error by a relative tolerance factor \(\varepsilon\), the required number of iterations satisfies \[
\begin{align}
k \approx \frac{\ln(1/\varepsilon)}{-\ln \rho(\bT)}.
\end{align}
\] If \(\rho(\bT) = 0.99\), we have \(-\ln(0.99) \approx 0.01005\), which requires approximately \(230\) iterations to reduce the error by a factor of ten. As \(\rho(\bT) \to 1\), convergence slows drastically.
Krylov Subspaces and Matrix-Free Solvers
While stationary methods rely on fixed matrix splittings, modern iterative methods project the problem onto expanding polynomial subspaces generated by the initial residual.
Let \(\bA \in \fR^{n \times n}\) and let \(\br^{(0)} \in \fR^n\) be a nonzero vector. The \(k\)-th Krylov subspace generated by \(\bA\) and \(\br^{(0)}\) is the linear subspace \[
\begin{align}
\mathcal{K}_k(\bA, \br^{(0)}) = \operatorname{span} \left\{ \br^{(0)}, \bA\br^{(0)}, \bA^2\br^{(0)}, ..., \bA^{k-1}\br^{(0)} \right\}.
\end{align}
\]
Let \(\bx^{(0)}\) be an initial guess with residual \(\br^{(0)} = \bb - \bA\bx^{(0)}\). Any iterate \(\bx^{(k)}\) belonging to the affine Krylov space \(\bx^{(0)} + \mathcal{K}_k(\bA, \br^{(0)})\) can be represented as \[
\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\). The corresponding residual vector satisfies \[
\begin{align}
\br^{(k)} = \bb - \bA\bx^{(k)} = q_k(\bA)\br^{(0)},
\end{align}
\] where \(q_k(t) = 1 - t p_{k-1}(t)\) is a polynomial of degree at most \(k\) satisfying the normalization constraint \(q_k(0) = 1\).
(Matrix-free computation) Krylov subspace methods never require explicit access to the individual entries of \(\bA\). They interact with the linear operator exclusively by computing matrix-vector products \(\bv \mapsto \bA\bv\). Consequently, linear systems can be solved efficiently even when the matrix is defined implicitly as a differential operator or a numerical simulation routine.
Because the true error \(\be^{(k)} = \bx^{(k)} - \bx^*\) is unobservable, termination is determined using the computable residual vector \(\br^{(k)} = \bb - \bA\bx^{(k)}\). Common stopping criteria include:
Relative residual test: \(\frac{\|\br^{(k)}\|_2}{\|\bb\|_2} \leq `tol`\).
Preconditioned relative residual test: \(\frac{\|\bM^{-1}\br^{(k)}\|_2}{\|\bM^{-1}\bb\|_2} \leq `tol`\).
Step difference test: \(\frac{\|\bx^{(k+1)} - \bx^{(k)}\|_2}{\|\bx^{(k+1)}\|_2} \leq `tol`\).
Consider the matrix \(\bA = \begin{pmatrix} 4 & 1 \\ 2 & 3 \end{pmatrix}\) and the splitting \(\bM = 4\bI\). Compute the iteration matrix \(\bT = \bI - \bM^{-1}\bA\), find its eigenvalues, and evaluate its spectral radius \(\rho(\bT)\).
For the splitting above, determine the number of iterations required to reduce the initial error by a factor of \(10^{-6}\).
Prove that if \(\|\bT\| < 1\) for some induced matrix norm, then the error satisfies \(\|\be^{(k)}\| \leq \frac{\|\bT\|}{1 - \|\bT\|} \|\bx^{(k)} - \bx^{(k-1)}\|\).
Let \(\bA \in \fR^{n \times n}\) and let \(\br^{(0)} \in \fR^n\). Prove that the sequence of dimensions \(\dim \mathcal{K}_k(\bA, \br^{(0)})\) strictly increases by one at each step until it reaches an invariant subspace of \(\bA\), after which the dimension remains constant.