Preconditioned steepest descent (source code)

= Preconditioned steepest descent
{title2=$Mz_k=r_k,\quad \alpha_k=r_k^Tz_k/(z_k^TAz_k)$}

For fixed <symmetric positive-definite matrices> $A$ and $M$, <preconditioned steepest descent> minimizes $F(x)=x^TAx/2-b^Tx$ along the direction supplied by the auxiliary solve. Given $x_k$, form the <linear-system residual> $r_k=b-Ax_k$. If $r_k=0$, stop: $x_k=x_*$. Otherwise solve $Mz_k=r_k$ and use
$$
\alpha_k=\frac{r_k^Tz_k}{z_k^TAz_k},\qquad x_{k+1}=x_k+\alpha_kz_k,\qquad r_{k+1}=r_k-\alpha_kAz_k.
$$
This is <exact line search for a positive-definite quadratic> along $z_k$. Indeed,
$$
F(x_k+\alpha z_k)=F(x_k)-\alpha r_k^Tz_k+\frac12\alpha^2z_k^TAz_k,
$$
whose <derivative> vanishes at the displayed positive <step size>. Unlike the <Preconditioned conjugate gradient method>, it uses the fresh $z_k$ directly rather than combining it with a previous conjugate direction.

For <diagonal scaling for a two-dimensional linear system>, start at $x_0=0$. Without <matrix preconditioning>, take $M=I$: $z_0=(1,1)^T$, $\alpha_0=10/507$, and
$$
x_1=\frac1{507}\begin{pmatrix}10\\10\end{pmatrix},\qquad r_1=\frac{165}{169}\begin{pmatrix}1\\-1\end{pmatrix}.
$$
With the <Jacobi preconditioner> $M=\operatorname{diag}(1,100)$, $z_0=(1,1/100)^T$, $\alpha_0=505/507$, and
$$
x_1=\begin{pmatrix}505/507\\101/10140\end{pmatrix},\qquad r_1=\begin{pmatrix}33/16900\\-33/169\end{pmatrix}.
$$
The original <linear-system residual> becomes much smaller after this correction. Quantify progress in the same <energy norm> in both cases: $E_k=F(x_k)-F(x_*)$ and $E_0=2515/4998$. Exact substitution gives
$$
\frac{E_1}{E_0}=q_I=\frac{81675}{85007}\approx0.960803228,\qquad \frac{E_1^{(M)}}{E_0}=q_M=\frac{3267}{8500700}\approx0.000384321291.
$$
In this particular example the entire histories are $E_k/E_0=q_I^k$ and $E_k^{(M)}/E_0=q_M^k$. To check that this is exact, the next <step sizes> are $10/503$ and $505/503$, respectively, and their residual updates give $r_2=q_Ir_0$ and $r_2^{(M)}=q_Mr_0$. Multiplying a <linear-system residual> by a nonzero <scalar> leaves the next <step size> unchanged; every pair of steps repeats the same directions at a smaller scale. The two first-step energy ratios then establish these formulas for all $k$.

\Image[/matrix-preconditioning-steepest-descent-energy-errors.png]
{title=Computed <steepest descent method> energy errors before and after <Jacobi preconditioning>}

The curves come from exact rational updates of the specified <linear system>, converted to decimal values for plotting. The left panel counts iterations. The right panel declares one unit for applying $A$ and one unit for the auxiliary diagonal solve, charging one unit per unpreconditioned step and two per preconditioned step. It excludes setup and <vector> operations, so it compares a declared cost model rather than measured runtime. Scaling is advantageous even with that extra charge.

These curves compare <steepest descent method> iterations. On this two-dimensional <linear system>, both the <Conjugate gradient method> and the <Preconditioned conjugate gradient method> terminate within two steps in exact arithmetic. The plot therefore establishes no CG iteration-count advantage for this example; the benefit of scaling for larger problems depends on their <eigenvalues> and auxiliary-solve costs.