2.3.5 Conjugate Gradient Method

Task: Solve $\mbox{${\bf A}$} \cdot \vec{x}=\vec{b}$. This may be interpreted as a minimization problem:

Let

$$f(\vec{x}) \equiv \frac{1}{2} \, \left\vert \mbox{${\bf A}$} \cdot \vec{x} - \vec{b} \right\vert^{2} \, ,$$

Then the $N$-vector $\vec{x}$ which minimizes $f(\vec{x})$ (with the minimum value $f=0$) is the desired solution.



Minimization of a quadratic function of $N$ variables:

Figure: Conjugate gradients: $\vec{g}_{0}=-\vec{\nabla} f(P_{0})$ ... direction of steepest descent at $P_{0}$, etc.; $\vec{h}_{1}$ ... direction of the gradient conjugate to $\vec{g}_{0}$.

Warm-up: Cauchy's Steepest Descent Method

Assume $f$ to depend on two variables $\vec{x}=(x_{1},x_{2})$ only. The lines of equal elevation of a quadratic function are ellipses that may be very elongated (see the figure).

• Start from $P_{0}$ with position vector $\vec{x}_{0}$ and follow the local gradient
  

  $$\vec{g}_{0} = - \nabla f(P_{0})$$ (2.2)

  
  
  This direction will not lead directly to the extremal point of $f$.
• Find the lowest point $P_{1}$ along this path. Determine once more the local gradient $\vec{g}_{1}= - \nabla f(P_{1})$; by construction it must be perpendicular to the previous direction $\vec{g}_{0}$.
• Iterate this procedure to arrive, after many mutually orthogonal bends, at the bottom of the channel (or top of the hill).

Now for the serious work: Conjugate Gradients

From point $P_{1}$ take the direction $\vec{h}_{1}$ instead of $\vec{g}_{1}$.

How to find $\vec{h}_{1}$?

Require that along $\vec{h}_{1}$ the change of the gradient of $f$ has no component parallel to $\vec{g}_{0}$. (In contrast: along $\vec{g}_{1}$, the gradient of $f$ has initially no $\vec{g}_{0}$-component). As a consequence, a gradient in the direction of $\vec{g}_{0}$ will not develop immediately (on quadratic surfaces never!).

$\Longrightarrow$This requirement leads to the prescription given in the Figure.

Figure 2.2: The CG method

For dimension $N=2$ this is all there is; $\vec{x} = \vec{x}_{2}$ is the solution vector. For systems of higher dimension one has to go on from $\vec{x}_{2}$ in the direction

$$\mbox{$\vec{h}_{2}$} = \mbox{$\vec{g}_{2}$} - \frac{(\mbox{$... ...\vert \mbox{${\bf A}$} \cdot \vec{h}_{1} \right\vert^{2}} \, .$$

until the next ``low point'' is reached at

$$\vec{x}_{3} = \vec{x}_{2} + \lambda_{3} \, \mbox{$\vec{h}_{2}$} \, ,$$

with

$$\lambda_{3} = \frac{\left\vert \mbox{$\vec{g}_{2}$} \cdot \m... ...ox{${\bf A}$} \cdot \mbox{$\vec{h}_{2}$} \right\vert^{2}} \, .$$

Next we determine $\vec{x}_{4}$ and so forth, until we arrive at the solution $\vec{x} = \vec{x}_{N}$.

Advantage of CG: No matrix inversion!

But: Several multiplication $\mbox{${\bf A}$} \cdot \vec{x}$; therefore most efficient for sparse matrices.

EXAMPLE: To see the mechanics of the method, assume

$$\mbox{${\bf A}$}=\mbox{$\left( \begin{array}{cc}3&1\\ \vspac... ...in{array}{r} 1.2 \\ \vspace{-9 pt}\\ 0.2 \end{array} \right)$}$$

The gradient vector at $\vec{x}_{0}$ is

$$\mbox{$\vec{g}_{0}$} = - \mbox{${\bf A}$}^{T} \cdot \left[ \... ...in{array}{r} 4.8 \\ \vspace{-9 pt}\\ 5.6 \end{array} \right)$}$$

and

$$\lambda_{1} = \frac{\left\vert \mbox{$\vec{g}_{0}$} \right\v... ...${\bf A}$} \cdot \mbox{$\vec{g}_{0}$} \right\vert^{2}} = 0.038$$

such that

$$\vec{x}_{1}=\vec{x}_{0}+\lambda_{1}\vec{g}_{0}=\mbox{$\left(... ...ray}{r} 1.017 \\ \vspace{-9 pt}\\ -0.014 \end{array} \right)$}$$

Similarly,

$$\mbox{$\vec{g}_{1}$}=\mbox{$\left( \begin{array}{r} -0.063 \... ... 0.0532 \end{array} \right)$} \, , \;\; \lambda_{2}=0.262 \, ,$$

and thus

$$\vec{x}_{2}= \mbox{$\left( \begin{array}{r} 1.000 \\ \vspace{-9 pt}\\ 0.000 \end{array} \right)$}$$

which is the desired solution vector in this 2-dimensional case. If $\mbox{${\bf A}$}$ were a $N \times N$ matrix, $N$ such steps would be necessary.