22  Jacobi and Gauss-Seidel Methods

Standard stationary iterations based on diagonal and triangular splittings.

22.1 Jacobi Iteration

Decomposes \(\bA\) into diagonal \(\bD\) and remainder \(\bR\): \(\bA = \bD + \bR\).

NoteDefinition: Jacobi Scheme

\[ \begin{align} \bx^{(k+1)} = \bD^{-1}(\bb - \bR\bx^{(k)}) \end{align} \]

  • Component form: \(x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j \neq i} a_{ij}x_j^{(k)}\right)\).

  • Concurrency: All components update independently using values from the previous step.

TipRemark

(Jacobi performance) Jacobi iteration is inherently parallel, as all components can be updated simultaneously. However, it often exhibits slow convergence and requires double buffering of the solution vector.

22.2 Gauss-Seidel Iteration

Splits \(\bA = \bL_* + \bU_*\), where \(\bL_*\) includes the diagonal.

NoteDefinition: Gauss-Seidel Scheme

\[ \begin{align} \bL_* \bx^{(k+1)} = \bb - \bU_* \bx^{(k)} \end{align} \]

  • Component form: \(x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j < i} a_{ij}x_j^{(k+1)} - \sum_{j > i} a_{ij}x_j^{(k)}\right)\).

  • Update rule: Use updated values immediately within the same sweep.

TipRemark

(Gauss-Seidel performance) Sequential dependencies restrict efficient parallelization, but Gauss-Seidel typically converges twice as fast as Jacobi and allows for in-place updates.

22.3 Convergence Conditions

NoteTheorem: Stationary Iteration Matrices

For \(\bA=\bD+\bL+\bU\), where \(\bL\) and \(\bU\) are the strict lower and upper triangular parts, \[ \begin{align} \bT_J = -\bD^{-1}(\bL+\bU) \end{align} \] is the Jacobi iteration matrix, and \[ \begin{align} \bT_{GS}=-(\bD+\bL)^{-1}\bU \end{align} \] is the Gauss-Seidel iteration matrix. By the result above, each method converges for every initial guess exactly when the corresponding spectral radius is less than \(1\).

WarningExercise

  1. Rewrite the result above in the form \(\bx^{(k+1)}=\bT_J\bx^{(k)}+\bc_J\).

  2. Rewrite the result above in the form \(\bx^{(k+1)}=\bT_{GS}\bx^{(k)}+\bc_{GS}\).

  3. Apply the result above to prove the convergence criterion in the result above.

  4. Compute \(\bT_J\) and \(\bT_{GS}\) for \(\bA=\begin{pmatrix}2&1\\1&2\end{pmatrix}\).

NoteTheorem: Diagonal Dominance

If \(\bA\) is strictly diagonally dominant (\(|a_{ii}| > \sum_{j \neq i} |a_{ij}|\)), both Jacobi and Gauss-Seidel converge for any \(\bx^{(0)}\).

WarningExercise

  1. For Jacobi, use strict diagonal dominance to show that the infinity-norm row sums of \(\bT_J\) are less than \(1\).

  2. Conclude \(\rho(\bT_J)\leq \|\bT_J\|_\infty<1\).

  3. Repeat the argument for Gauss-Seidel by isolating the largest component of the error equation.

NoteTheorem: SPD Convergence

If \(\bA\) is SPD, Gauss-Seidel is guaranteed to converge. Jacobi may still diverge.

WarningExercise

  1. For \(\bA = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\), compute one step of Jacobi and Gauss-Seidel from \(\bzero\).

  2. Compare \(\rho(\bT_J)\) and \(\rho(\bT_{GS})\). Does \(\rho_{GS} \approx \rho_J^2\) hold?

22.4 Successive Over-Relaxation (SOR)

Accelerates Gauss-Seidel by overshooting the update.

NoteDefinition: SOR Update

\[ \begin{align} x_i^{(k+1)} = (1-\omega)x_i^{(k)} + \omega \left[\text{GS update}\right] \end{align} \]

NoteTheorem: SOR Convergence Window

For SPD systems, SOR converges for relaxation parameters \[ \begin{align} 0<\omega<2. \end{align} \] For model Poisson problems, well-chosen \(\omega>1\) can reduce the iteration count substantially compared with Gauss-Seidel.

TipRemark

(SOR acceleration) For the 1D Poisson problem on a grid with spacing \(h\):

  • Jacobi and Gauss-Seidel require \(O(n^2)\) iterations.

  • Optimal SOR requires only \(O(n)\) iterations.

SOR effectively improves the convergence rate from \(O(h^2)\) to \(O(h)\).

WarningExercise

  1. Optimal Tuning: For \(\rho_J = 0.95\), calculate the optimal \(\omega^* = 2/(1 + \sqrt{1-\rho_J^2})\).

  2. Convergence Sweep: Implement SOR and plot iterations to convergence vs. \(\omega \in [1, 2)\). Identify the ``cliff’’ at \(\omega^*\).

  3. Check numerically that choices \(\omega\leq 0\) or \(\omega\geq 2\) fail on an SPD Poisson matrix, as predicted by the result above.

22.5 Stopping Criteria

NoteDefinition: Iterative Stopping Tests

Common stopping tests include

  1. absolute residual: \(\|\br^{(k)}\| < \varepsilon\),

  2. relative residual: \(\|\br^{(k)}\|/\|\bb\| < \varepsilon\),

  3. increment: \(\|\bx^{(k+1)}-\bx^{(k)}\|<\varepsilon\).

Residual tests are tied to the equation; increment tests are cheap but can stop prematurely if convergence is slow.

WarningExercise

  1. Why is a small residual not enough to guarantee a small error \(\be^{(k)}\)? Use the result above.

  2. Implement a robust solver that combines relative residual with a max\_iter budget.