Preconditioned Conjugate Gradients algorithm

Let’s consider each step of the naive CG algorithm with a preconditioner $C$. This is the step where we multiply by the matrix:

$$\begin{array}{rl}{q}^{(k)}& =CA{C}^T{p}^{(k)}\\ {\mu }_k& =\frac{\Vert s^{(k-1)}{\Vert }_2^2}{(p^{(k)}{)}^T{q}^{(k)}}\end{array}$$

This suggests using the following definition:

$${\underline{p}}^{(k)}=C^T{p}^{(k)}.$$

The following recurrence will work

$${\underline{q}}^{(k)}=A{C}^T{p}^{(k)}=A{\underline{p}}^{(k)}.$$

Note that $C$ is missing on the left.

Then we have

$${\mu }_k=\frac{\Vert s^{(k-1)}{\Vert }_2^2}{({\underline{p}}^{(k)}{)}^T{\underline{q}}^{(k)}}$$

with $q^{(k)}=C{\underline{q}}^{(k)}$ from our definition of ${\underline{q}}^{(k)}$.

Solution update step

The solution update step for $y^{(k)}$ is:

$$y^{(k)}=y^{(k-1)}+{\mu }_k{p}^{(k)}$$

Recall that $x^{(k)}=C^T{y}^{(k)}.$ Let’s multiply the equation above by $C^T$

$$C^T{y}^{(k)}=C^T{y}^{(k-1)}+{\mu }_k{C}^T{p}^{(k)}$$

This is just

$$x^{(k)}=x^{(k-1)}+{\mu }_k{\underline{p}}^{(k)}$$

Residual

Let’s focus on the residual update step. The definition of the residual for CG with symmetric preconditioning is

$$s^{(k)}=Cb-CA{C}^T{y}^{(k)}=Cb-CA{x}^{(k)}$$

Recall that

$$r^{(k)}=b-A{x}^{(k)}$$

So:

$$s^{(k)}=C{r}^{(k)}.$$

Let’s see how we can update the residual.

$$s^{(k)}=s^{(k-1)}-{\mu }_{k}C{\underline{q}}^{(k)}$$

Let’s rewrite this in terms of the residual $r^{(k)}$:

$$C{r}^{(k)}=C{r}^{(k-1)}-{\mu }_{k}C{\underline{q}}^{(k)}$$

We multiply by $C^{-1}$ to the left to simplify and get our final expression:

$$r^{(k)}=r^{(k-1)}-{\mu }_k{\underline{q}}^{(k)}$$

Search direction

The search direction update is:

$$p^{(k+1)}=s^{(k)}+{\tau }_k{p}^{(k)}$$

With our definition ${\underline{p}}^{(k)}=C^T{p}^{(k)},$ we get:

$$C^{-T}{\underline{p}}^{(k+1)}=s^{(k)}+{\tau }_k{C}^{-T}{\underline{p}}^{(k)}$$

Multiply to the left by $C^T$:

$${\underline{p}}^{(k+1)}=C^T{s}^{(k)}+{\tau }_k{\underline{p}}^{(k)}$$

Let’s introduce the residual $r^{(k)}$. Recall that:

$$s^{(k)}=C{r}^{(k)}$$

So

$$C^T{s}^{(k)}=C^{T}C{r}^{(k)}=M{r}^{(k)}.$$

We can substitute in the previous equation:

$${\underline{p}}^{(k+1)}=M{r}^{(k)}+{\tau }_k{\underline{p}}^{(k)}$$

In PCG, we use a new temporary vector $z^{(k)}=M{r}^{(k)}$ to reduce the time cost. We finally get:

$${\underline{p}}^{(k+1)}=z^{(k)}+{\tau }_k{\underline{p}}^{(k)}$$

Norm of residual

We are almost done. The last quantity we need is $\Vert s^{(k)}{\Vert }_2.$ We again use the relation for the residual $s^{(k)}=C{r}^{(k)}.$ So

$$\Vert s^{(k)}{\Vert }_2^2=(r^{(k)}{)}^T{C}^{T}C{r}^{(k)}=(r^{(k)}{)}^{T}M{r}^{(k)}$$

As desired, $C$ has disappeared in favor of $M$!

We can even use our temporary work vector $z^{(k)}=M{r}^{(k)}$ to simplify the expression further:

$$\Vert s^{(k)}{\Vert }_2^2=[r^{(k)}{]}^T{z}^{(k)}$$

This is the last step in the Preconditioned CG algorithm!

Summary of PCG

All the steps above are summarized to form the PCG algorithm:

$$\begin{array}{rl}{\underline{q}}^{(k)}& =A{\underline{p}}^{(k)}\\ {\mu }_k& =\frac{[r^{(k-1)}{]}^T{z}^{(k-1)}}{({\underline{p}}^{(k)}{)}^T{\underline{q}}^{(k)}}\\ x^{(k)}& =x^{(k-1)}+{\mu }_k{\underline{p}}^{(k)}\\ r^{(k)}& =r^{(k-1)}-{\mu }_k{\underline{q}}^{(k)}\\ z^{(k)}& =M{r}^{(k)}\\ {\tau }_k& =\frac{[r^{(k)}{]}^T{z}^{(k)}}{[r^{(k-1)}{]}^T{z}^{(k-1)}}\\ {\underline{p}}^{(k+1)}& =z^{(k)}+{\tau }_k{\underline{p}}^{(k)}\end{array}$$

We see that this is only a very small modification of the original CG algorithm!