Mathematics.

numerical linear algebra

Iterative Solvers for Linear Systems

Numerical Analysis30 minDifficulty7 out of 10

You should know: jacobi and gauss seidel methods

Overview

Simple iterative methods like Jacobi and Gauss-Seidel reduce error at every step, but they do so unevenly: they are very effective at damping high-frequency (rapidly oscillating) components of the error, yet painfully slow at removing smooth, low-frequency error components, so residual reduction stalls after a promising start on real, large-scale problems such as PDE discretizations. Multigrid methods exploit this observation directly: after a few 'smoothing' sweeps of Jacobi or Gauss-Seidel wipe out the high-frequency error, the remaining smooth error is transferred ('restricted') to a coarser grid, where it looks relatively higher-frequency and can again be smoothed away cheaply; the correction found on the coarse grid is then transferred back ('prolongated') to the fine grid. Applied recursively across a whole hierarchy of grids — the essence of the V-cycle and W-cycle schedules — multigrid can solve many elliptic PDE-derived linear systems in a number of operations that scales linearly with the number of unknowns, in stark contrast to the O(n³) of direct elimination or the slow, ever-larger iteration counts required by Jacobi/Gauss-Seidel alone on fine grids.

Intuition

Imagine trying to smooth out a lumpy bedsheet by only ever tugging at points a few centimeters apart: you can quickly flatten small wrinkles right under your hands, but a long, gentle wave running the whole length of the sheet barely responds no matter how many times you tug — your hands are simply too close together to 'see' that large-scale wave. Multigrid solves this by periodically stepping back and working with a coarser grid of hands spaced further apart: what looked like a slow, smooth wave up close now looks like a wrinkle those widely-spaced hands can fix quickly. Repeating this zoom-out, fix, zoom-back-in process across several grid resolutions lets every spatial scale of error get flattened efficiently by whichever grid resolution is best suited to seeing it.

Formal Definition

Definition

A two-grid multigrid cycle alternates smoothing on the fine grid with an exact (or near-exact) correction on a coarser grid, using restriction and prolongation operators to move information between grids:

x(k+1)=S(x(k)),S=a Jacobi or Gauss-Seidel smoothing sweepx^{(k+1)} = S(x^{(k)}), \qquad S = \text{a Jacobi or Gauss-Seidel smoothing sweep}
Pre-smoothing step on the fine grid
r(k)=bAx(k),r2h=Rrh(restriction to a coarser grid with spacing 2h)r^{(k)} = b - Ax^{(k)}, \qquad r_{2h} = R\, r_h \quad \text{(restriction to a coarser grid with spacing } 2h\text{)}
Residual computation and restriction to the coarse grid
A2he2h=r2h  (solved/approximated on the coarse grid),xhxh+Pe2hA_{2h}e_{2h} = r_{2h} \;\text{(solved/approximated on the coarse grid)}, \qquad x_h \leftarrow x_h + P\, e_{2h}
Coarse-grid correction and prolongation (P) back to the fine grid

Worked Examples

  1. Iteration 1 gives x⁽¹⁾=(1.2,1.3,1.4); compute the residual b − Ax⁽¹⁾ component-wise.

    r(1)=(12,13,14)A(1.2,1.3,1.4)T    r(1)26.836r^{(1)} = (12,13,14) - A(1.2,1.3,1.4)^T \;\Rightarrow\; \|r^{(1)}\|_2 \approx 6.836
  2. Continue the Jacobi sweeps (as computed in Jacobi/Gauss-Seidel methods): x⁽²⁾=(0.93,0.92,0.90), x⁽³⁾=(1.018,1.024,1.03).

    r(2)21.883,r(3)20.544\|r^{(2)}\|_2 \approx 1.883, \qquad \|r^{(3)}\|_2 \approx 0.544
  3. Continuing to x⁽⁴⁾=(0.9946,0.9934,0.9916) and x⁽⁵⁾=(1.0015,1.0019,1.0024).

    r(4)20.154,r(5)20.0439\|r^{(4)}\|_2 \approx 0.154, \qquad \|r^{(5)}\|_2 \approx 0.0439

Answer: Residual norms 6.836 → 1.883 → 0.544 → 0.154 → 0.0439: each sweep shrinks the residual by roughly a constant factor of about 0.28-0.29, the hallmark of a simple iterative smoother — steady but geometrically slow, which is exactly the behavior multigrid accelerates by handling the stubborn smooth error components on a coarser grid instead of grinding them out with many more fine-grid sweeps.

Practice Problems

Difficulty 5/10

Using the residual sequence 6.836, 1.883, 0.544, 0.154, 0.0439 from the worked example, estimate the approximate ratio between consecutive residual norms (r₂/r₁).

Difficulty 6/10

If plain Jacobi/Gauss-Seidel smoothing scales as O(n²) total operations to converge on an n-unknown discretized PDE, and multigrid scales as O(n), what is the operation-count ratio at n = 10,000?

Difficulty 8/10

A computational fluid dynamics (CFD) solver needs to solve a pressure-Poisson equation on a 3D grid with tens of millions of cells at every timestep of a simulation. Explain why multigrid, rather than plain Gauss-Seidel iteration alone, is the standard choice in production CFD codes.

Quiz

Simple iterative smoothers like Jacobi and Gauss-Seidel are effective at reducing:
In a multigrid cycle, the residual is transferred to a coarser grid using the:
The main computational advantage of multigrid over plain Jacobi/Gauss-Seidel on large discretized PDE systems is that multigrid:

Summary

  • Jacobi/Gauss-Seidel smoothers reduce high-frequency error quickly but stall on smooth, low-frequency error, e.g. residual norms 6.836 → 1.883 → 0.544 → 0.154 → 0.0439 shrink by a roughly constant ~0.28 factor per sweep rather than accelerating.
  • Multigrid restricts the smooth residual to a coarser grid (where it looks higher-frequency), corrects it there cheaply, then prolongates the correction back to the fine grid.
  • Applied recursively across a grid hierarchy (V-cycles/W-cycles), multigrid achieves near-O(n) total cost versus O(n²) or worse for plain smoothing alone as resolution grows, e.g. a ~1000x speedup at n=1000 unknowns.
  • This makes multigrid the standard choice for solving the large elliptic-PDE-derived linear systems that arise every timestep in production CFD and other large-scale simulation codes.

References

  1. BookBriggs, W. L., Henson, V. E., & McCormick, S. F. A Multigrid Tutorial, 2nd ed. SIAM, 2000.