Fundamentals of Multigrids Methods

Some good references: McCormick's ``Multilevel Methods for PDE's'' (SIAM) and ``Multigrid Tutorial'' by W. Briggs (SIAM). These are very accessible and inexpensive SIAM books. Here, we'll follow Iserles rather closely.

A good software and information source for multigrid

Multigrid methods are nested techniques for the iterative solution of the linear algebraic problem

$$A {\bf x} = {\bf b}$$

and is a current active research area. Since they are fast, they are used not only in the solution of PDE's but also in applications related to image processing, filtering, etc.

Suppose we want to solve the 5-point FD approximation of the Poisson equation using Gauss-Seidel. As we have discussed in 475A, the rate of convergence of the iterative solution $\rho (A)$, which in this case is approximately $1-\pi^2/m^2$ for the $m\times m$ matrix $A$. This is an asymptotic result. In reality, we'd see that $\ln\vert\vert r^k\vert\vert$, where ${\bf r}^k=A{\bf x}^k-{\bf b}$, will drop as shown in Figure 38

Figure 38: The logarithm of the norm of the residual as a function of the iteration count $k$ for Gauss-Seidel on the Poisson problem

we see a severe drop in the first few iterates, followed by the linear rate predicted by the asymptotic result. This is true for any $m$!!

Why? Because the Gauss Seidel acts as a ``smoother,'' altervating high wave numbers faster than low wave numbers. Understanding why provides a technique for accelerating iterative schemes:

Subtract the 5-point equations

$$u_{j-1,l}+u_{j,l-1}+u_{j+1,l}+u_{j,l+1}-4u_{j,l}=\Delta x^2 f_{j,l}\quad j,l=1,2\ldots m$$ (180)

from Gauss Seidel Scheme:

$$u^{k+1}_{j-1, l}+u^{k+1}_{j,l-1}+u^k_{j+1,l}+u^k_{j,l+1}-4u^{k+1}_{j,l}= \Delta x^2 f_{j,l}\quad j,l=1,2,\ldots m$$

$$\mbox {to obtain }$$ (181)

$$\varepsilon ^{k+1}_{j-1,l}+\varepsilon ^{k+1}_{j,l-1}+\varepsilon... ...+1,l}+\varepsilon ^k_{j,l+1} -4\varepsilon ^{k+1}_{j,l}=0 \quad j,l=1,2\ldots m$$ (182)

where $\varepsilon ^k_{j,l}\equiv u^k_{j,l}-u_{j,l}$ is the error after $k$ iterations at the $(j,l)$ grid point.

Since we're assuming Dirichlet boundary conditions $u_{j,l}$ and $u^k_{j,l}$ are identical at the boundary $\Rightarrow \varepsilon ^k_{j,l}=0$ there.

Let

$$p^k(\theta,\psi)=\sum^m_{j=1}\sum^m_{l=1}\varepsilon ^k_{j,l}e^{i(j\theta +l\psi)}, \qquad 0\le \theta\psi\le 2\pi$$

be the 2-D Fourier transform of the sequence $$\{\varepsilon ^k_{j,l}\}^m_{j,l=1}$$.

and denote the Euclidean Norm

$$\vert\vert\vert g\vert\vert\vert=\left[\frac{1}{4\pi^2}\int^{\pi}... ...i}\int^{\pi}_{-\pi}\vert g(\theta,\psi)\vert^2 d\theta d\psi\right]^{\frac12}.$$

By Parseval's Theorem, can show that

$$\vert\vert\vert p^k\vert\vert\vert^2=\vert\vert\varepsilon ^k\vert\vert^2$$

where $$\vert\vert y\vert\vert=\left(\sum^m_{j=1}\sum^{m}_{l=1}\vert y_{j,l}\vert^2\right)^{\frac12}$$.

We wish to establish the rate of decay of residuals. Multiply (183) by $e^{i(j\theta+l\psi)}$ and sum over $j,l=1,2\ldots m$

Since

$$\sum^m_{j=1}\sum^{m}_{k=-1}\varepsilon _{j-1,l}e^{i(j\theta+l\psi... ...1}\varepsilon ^{k+1}_{j,l}e^{i[(j+1)\theta +l\psi]}=e^{i\theta}p^k(\theta,\psi)$$

$$-e^{i(m+1)\theta}\sum^m_{l=1}\varepsilon _{m,l} e^{i l \theta}$$

and applying similar algebra to other term in (183) we obtain

$$(4-e^{i\theta}-e^{i\psi})p^{k+1}(\theta,\psi) = (e^{-i\theta}+e^ {-i\psi})p^k(\theta,\psi)$$

$$- \left\{e^{i(m+1)\theta}\sum^{m}_{l=1}\varepsilon ^{k+1}_{m,l}e^{il\psi} +e^{i(m+1)\psi} \sum^u_{j=1}\varepsilon ^{k+1}_{jm} e^{ij\theta}\right.$$

$$+ \left.\sum^m_{l=1}\varepsilon ^k_{1,l}e^{il\psi}+ \sum^m_{j=1}\varepsilon ^k_{j-1}e^{ij\theta}\right\}$$

Now $(4-e^{i\theta}-e^{i\psi})p^{k+1}(\theta,\psi)\approx (e^{-i\theta}+e^ {-i\psi})p^k(\theta,\psi)$ in general but the term in curly brackets would have disappeared if we were considering periodic boundary conditions. For Dirichlet, the justification is more subtle. We could check a posteriori that it is indeed ok.

Define the local attenuation factor as

$$\rho^k(\theta,\psi)=\Big\vert\;\frac{p^{k+1} (\theta,\psi)}{p^k(\theta,\psi)}$$

for $\theta\vert\le \pi$ and $\vert\psi\vert\le\pi$,

then $\rho^k(\theta,\psi)\approx \overline\rho(\theta,\psi)=\Big\vert \frac{e^{-i\the... ...{i\theta}-e^{i\psi}} \Big\vert\; \qquad \vert\theta\vert, \vert\psi\vert \le\pi$.

$\bar \rho$ is independent of $k$.

In HW10 you will confirm graphically that when $\tilde \rho$ is restricted to the set $\Pi_0\equiv\left\{(\theta,\psi):\frac{1}{2}\pi\le\max\left\{\vert\theta\vert, \vert\psi\vert\right\}\le\pi\right\}$ the function $\tilde \rho$ peaks at $$\frac12$$, whereas $\tilde \rho$ over $\Pi\equiv \left\{(\theta,\psi):\vert\theta\vert\le \pi, \vert\psi\vert\le \pi\right\}$ peaks at 1. In fact, you will show that

$$\max_{(\theta,\psi)\in\Pi}\tilde \rho (\theta,\psi)=\tilde \rho \left(\frac{\pi}{2}, \tan^{-1}\left(\frac{3}{4}\right)\right)=\frac12$$

$\therefore$ if we disregard non-oscillatory wave numbers, the amplitude of the error is halved in each iteration!

$\therefore$ G-S attenuates highly oscillatory components much faster thus contribution of these vanishes quickly after a few iterations.

In the context of a continuum, all wavenumbers are supported. On a discretization of the continuum, a lattice or grid what is ``highly oscillatory depends on the grid spacing. In fact, for a specific grid there are oscillations that are not resolvable. Since what is meant by ``high oscillations'' is with respect to each grid realization, that is, high oscillations are those components with wavelengths that are comparable to the grid size, and the G-S iteration attenuates rapidly ``high oscillations'' we could conceive of an algorithm in which we go back and forth between coarse and fine meshes and iterate at each level ONLY while attenuation rates are high, we can get a fast rate of convergence for the residual by constructing a nesting sequence between coarse and fine meshes.

Suppose that we coarsen a grid by taking out every second point, the outcome being a grid on $[0,1]\times [0,1]$ but with $\Delta x$ replaced by $2\Delta x$. The range of former high frequencies $\Pi_0$ is no longer visible on the coarse grid. The new grid will have its own range of high frequencies on which G-S performs well:

$$\Pi_1=\left\{(\theta,\psi):\frac14\pi\le\max\left\{\vert\theta\vert,\vert\psi\vert\right\} \le\frac12\pi\right\}$$

We could coarsen again and again and form a hierarchy of grids embedded into each other, whose (grid-specific) high frequencies correspond, as far as the fine grid is concerned, to the sets

$$\Pi_s=\left\{(\theta,\psi):2^{-s-1}\pi\le\max\left\{\vert\theta\vert \psi\vert\right\}\le 2^{-s}\pi\right\}$$

$$s=1,2,3\cdots \log_2(m+1)$$

The sets $\Pi_s$ nest inside each other, as shown in Figure 39

Figure 39: The Nested subspaces

The multigrid technique takes advantage of this fact, traveling up and down the grid hierarchy, using G-S iterations to dampen the locally highly oscillatory components of the error (see Figure 40).

Figure 40: Schematic representation of the grid hierarchy

We need to describe how each coarsening or refinement step is performed, as well as to specify the exact strategy of how to start, when to coarsen, when to refine and when to terminate the whole process:

Refinement and Coarsening: Consider just 2 grids, find and coarse. Suppose we wish to solve

$$A_f {\bf x}_f={\bf v}_f$$$$\qquad \mbox {on a fine grid}$$

perform a few G-S iteratives (smoothing out high frequencies), then

$${\bf r}_f\equiv A_f{\bf x}_f-{\bf v}_f$$$$\mbox { is the residual.}$$

To go to coarse grid, use a ``restriction matrix'' $R$

$${\bf r}_c=R{\bf r}_f$$

remember, at this point ${\bf r}_f$ is constructed of low frequency components (relative to fine grid) $\therefore$ makes sense to go on smoothing ${\bf r}_c$ on coarser grid: let ${\bf v}_c\equiv -{\bf r}_c$ and solve

$$A_c{\bf x}_c=-{\bf r}_c$$ (183)

$A_c$ is the coarse matrix, i.e. A restricted to coarse grid.

To ascend, suppose ${\bf x}_c$ is approximate solution to (184) after a few iterations we translate ${\bf x}_c$ into the fine grid using a ``Prolongation matrix'' $P$?

$${\bf y}_f=P {\bf x}_c$$ (184)

and update the old value of ${\bf x}_f$:

$${\bf x}^{\mbox {new}}_f ={\bf x}^{\mbox {old}}_f+{\bf y}_f$$ (185)

Evaluating the residual ${\bf r}^{\mbox {new}}_f$, under the assumption that ${\bf x_c}$ is exact solution of (184). Since

$${\bf r}^{\mbox {new}}_f=A_f {\bf x}^{\mbox{new}}_f{\bf -v}_f=A_f({\bf x}^ {\mbox {old}}_f+{\bf y_f})-{\bf v}_f$$

then using (185) and (186):

$${\bf r}^{\mbox {new}}_{\bf f}={\bf r}^{\mbox {old}}_{{\bf f}}+A_f{\bf y}_f ={\bf r}_{f}^{\mbox {old}}+A_f P {\bf x}_c$$

$$\therefore {\bf r}^{\mbox {new}}_{f}={\bf r}^{\mbox {old}}_f-A_f P A^{-1}_c {\bf r}_c.$$

$$\mbox {Since }{\bf r}_c=R {\bf r}_f$$

$${\bf r}^{\mbox {new}}_f=({\bf I}-A_f P A^{-1}_c R){\bf r}^{\mbox{old}} _{\bf f}$$

$\therefore$ the sole contribution to the new residual comes from replacing a fine grid by a coarser one. similar reasoning is valid even if ${\bf x}_c$ is anapproximate solution of course grid problem provided high freq's have been smoothed out.

Now we need to specify the Restriction and Prolongation Matrices:

A popular $\left\{\begin{array}{l} \mbox {Restriction and Prolongation matrix comes from \lq\lq Full Weighting''}:\\ \end{array}\right.$

It leads to $R=\displaystyle \frac14P^T$, which is very convenient.

let

$${\bf w}_f = P{\bf w}c$$

$${\bf w}_c = \left(w^c_{j,l}\right)^{m}_{j,l=1}$$

$${\bf w}_f = \left(w^f_{j,l}\right)^{2m+1}_{j,l=1}$$

Full Weighting:

$$w^c_{j,l}=\frac14 w^{f}_{2j,2l}+\frac18\left(w^f_{2j-1,2l}+w^f_{2j,2l-1^t} +w^{f}_{2j+1,2l}+w^f_{2j,2l+1}\right)$$

$$+\frac{1}{16}\left(w^f_{2j-1,2l-1}+w^f_{2j+1,2l-1} +w^{f}_{2j-1,2l+1}+w^f _{2j+1,2l+1}\right)\quad j,l=1,2,\ldots m$$

Prolongation: use linear interpolation:

$$w^f_{2j-1, 2l-1}=w^c_{j,l}\quad j,l=1, 2, \ldots m$$

$$w^f_{2j-2, 2l}=\frac12\left(w^c_{j,l}+w^c_{j,l+1}\right)\quad j=1,2, \ldots m-1; l-1,2\ldots m$$

$$w^f_{2j,2l-1}=\frac12(w^c_{j,l}+w^c_{j+1,l})\quad j=1,2\cdots m; l=1,2 \cdots m-1$$

$$w^f_{2j-2l}=\frac14\left(w^c_{j,l}+w^c_{j,l+1}+w^c_{j+1,l}+w^c_{j+1,l+1} \right)\quad j,l=1,2\ldots m-1$$

$wf=0$ at boundary $\ldots$ recall we're dealing with residuals.

$\Box$

ALGORITHM V-CYCLE (One of many, the simplest and most popular) The algorithm is shown schematically in Figure 41.

Figure 41: The V-Cycle algorithm

Start and end on finest grid. To start, stipulate inital guess ${\bf v_f} ={\bf b}$ (original right-hand-side of system). Iterate Gauss-Seidel $n_r$ times.

Evaluate ${\bf r}_f$

Restrict it to coarser grid

Perform $n_r$ Gauss-Seidel iterations

Evaluate ${\bf r}_c$

Restrict on even coarser grid

and repeat process till we reach coarsest grid, with just one single interior point, which we can solve for exactly.

At this stage we've damped out the high frequencies of error, relative to each grid resolution. $\therefore$ damped influence of error components over entire range of wave numbers supported by the finest grid, except for small error introduced by restriction.

Now we go up all the way to the finest grid. At each step we Prolong,

Update residual on new grid,

and Perform no G-S iterations

to eliminate errors (corresponding to high oscillations for each grid resolution ) that might have been introduced by post prolongations.

We're back to the starting point, completed the $V$-Cycle.

Now we check for convergence by measuring size of residual vector.

If residual below some specified tolerance, we quit. Otherwise, repeat $V$-Cycle.

There is one problem with $V$-Cycle algorithm $\ldots$ it could be made faster if we start with a really good initial guess. The ``Full Multigrid'' method usually combines this with the pattern illustrated in Figure 42.

Figure 42: An improved V-cycle algorithm

See references for full details.

$\Box$

$V$-Cycle Computational Cost for Poisson (5-point formula):

let $\gamma$ be the cost of a single G-S iteration on finest grid. Note that a single coarsening decreases operation count by 4 $\therefore$ the cost of $V$-cycle is

$$\left(1+\frac14+\frac{1}{4^2}+\frac{1}{4^3}\cdots\right)(n_r+n_p)\gamma \approx \frac43(n_r+n_p)\gamma$$ (186)

where $n_r$ and $n_p$ are the $\char93$ of iterates in restriction and prolongation phases. $\therefore$ $V$-cycle is linear in the $\char93$ of grid points on finest grid.

Ex) For $m=63$, using $V$-cycle with $n_r=np=1$ for residual $\sim 10^{-5}\rightarrow 8^{th} V$-cycle for same residual: using

$$\begin{array}{ll} SOR & \rightarrow 243 \mbox { iterations}... ...s}\\ \mbox {Conjugate Gradient} & \rightarrow 179 \end{array}$$

Note: Cost of restriction and prolongation is not included in (187)

$\Box$