Preconditioned Steepest Descent Method (PSDM)

When I was implementing an iterative solver using Steepest Descent Method (SDM), I tried to search on the web how to implement a preconditioned version of the algorithm.

However, I only found posts in some forums asking the same question.

Therefore, I started reading this famous document on Conjugate Gradient (CG) – An Introduction to the Conjugate Gradient Method Without the Agonizing Pain by Jonathan Richard Shewchuk (http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf).

Then, I derived PSDM by myself.

The original SDM looks like this:

<br />
x = x0;<br />
r = b - Ax;<br />
For i = 1, 2, ... until convergence<br />
 q = Ar;<br />
 alpha = <r , r> / </r><r , q>;<br />
 x = x + alpha * r;<br />
 r = r - alpha * q;<br />
End<br />
</r>

Based on a symmetric positive definite preconditioner M, M = EE[sup:20nmn0md]T[/sup:20nmn0md] by Cholesky factorization the preconditioned coefficient matrix should be E[sup:20nmn0md]-1[/sup:20nmn0md]AE[sup:20nmn0md]-T[/sup:20nmn0md] because it is also symmetric positive definite the original system Ax = b is transformed to E[sup:20nmn0md]-1[/sup:20nmn0md]AE[sup:20nmn0md]-T[/sup:20nmn0md]y = E[sup:20nmn0md]-1[/sup:20nmn0md]b where y = E[sup:20nmn0md]T[/sup:20nmn0md]x

Thus, given z = M[sup:20nmn0md]-1[/sup:20nmn0md]r, we can update the symbols in our original algorithm: A -> A’ = E[sup:20nmn0md]-1[/sup:20nmn0md]AE[sup:20nmn0md]-T[/sup:20nmn0md] b -> b’ = E[sup:20nmn0md]-1[/sup:20nmn0md]b r -> r’ = E[sup:20nmn0md]-1[/sup:20nmn0md]r q -> q’ = A’r’ = E[sup:20nmn0md]-1[/sup:20nmn0md]Az -> r’[sup:20nmn0md]T[/sup:20nmn0md]r’ = r[sup:20nmn0md]T[/sup:20nmn0md]z -> r’[sup:20nmn0md]T[/sup:20nmn0md]q’ = (E[sup:20nmn0md]-1[/sup:20nmn0md]r)[sup:20nmn0md]T[/sup:20nmn0md] E[sup:20nmn0md]-1[/sup:20nmn0md]Az = r[sup:20nmn0md]T[/sup:20nmn0md]E[sup:20nmn0md]-T[/sup:20nmn0md]E[sup:20nmn0md]-1[/sup:20nmn0md]Az = z[sup:20nmn0md]T[/sup:20nmn0md]Az

The update for x becomes x = x + alpha * r -> x’ = x’ + alpha’ * r’ -> E[sup:20nmn0md]T[/sup:20nmn0md]x = E[sup:20nmn0md]T[/sup:20nmn0md]x + alpha’ * E[sup:20nmn0md]-1[/sup:20nmn0md]r multiply on both sides by E[sup:20nmn0md]-T[/sup:20nmn0md], becomes x = x + alpha’ z (since E[sup:20nmn0md]-T[/sup:20nmn0md]E[sup:20nmn0md]-1[/sup:20nmn0md] = M[sup:20nmn0md]-1[/sup:20nmn0md])

Similarly, the update for r becomes r = r – alpha * q -> r’ = r’ – alpha’ * q’ -> E[sup:20nmn0md]-1[/sup:20nmn0md]r = E[sup:20nmn0md]-1[/sup:20nmn0md]r – alpha’ * E[sup:20nmn0md]-1[/sup:20nmn0md]Az multiply on both sides by E, becomes r = r – alpha’ * q’ (since Az = q’)

And, PSDM looks like this:

<br />
x = x0;<br />
r = b - Ax;<br />
For i = 1, 2, ... until convergence<br />
 solve Mz = r;<br />
 q = Az;<br />
 alpha = <r , z> / <z , q>;<br />
 x = x + alpha * z;<br />
 r = r - alpha * q;<br />
End<br />
</z></r>

I hope this would help others who are also implementing the same method. Correct me if I miss anything.