21  Iterative Methods for Linear Systems

Alternative to direct solvers for large, sparse systems where matrix factorization is impractical due to memory or compute limits.

21.1 Direct vs. Iterative Solvers

Tip

The Sparsity Problem: 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:

  • \(\bA\) has \(O(n^2)\) nonzeros.

  • LU factors have \(O(n^3)\) nonzeros.

Iterative methods exploit sparsity by using only matrix-vector products (\(O(\text{nnz})\)).

NoteDefinition: Stationary Iteration

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).

21.2 Convergence and Spectral Radius

NoteDefinition: Spectral Radius

\(\rho(\bB) = \max \{|\lambda| : \lambda \in \sigma(\bB)\}\).

NoteTheorem: Gelfand’s Formula

For any submultiplicative norm, \(\rho(\bB) = \lim_{k\to\infty} \|\bB^k\|^{1/k}\). Insight: For large \(k\), \(\|\bB^k\| \approx \rho(\bB)^k\).

NoteTheorem: Fundamental Convergence Condition

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.

Tip

Practical Cost: 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.

WarningExercise

  1. For \(\bA = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}\), find \(\rho(\bT)\) for the splitting \(\bM = 3\bI\).

  2. If \(\rho(\bT) = 0.99\), how many iterations are needed to gain 4 digits of precision?

  3. Start \(\bx^{(0)}=\bzero\) and solve \(\bA\bx = (4, 4)^T\) via 3 steps of this iteration.

21.3 Fill-in and Sparse Storage

WarningExercise

  1. Fill-in Demo: Generate a large sparse tridiagonal matrix. Compare nonzeros in \(\bA\) vs. LU factors.

  2. Timing: Use scipy.sparse.linalg.spsolve vs. dense np.linalg.solve for \(n=2000\). Record memory and time.