Editing Conjugate gradient method (section)

==The preconditioned conjugate gradient method==
{{See also|Preconditioner}}
In most cases, [[preconditioning]] is necessary to ensure fast convergence of the conjugate gradient method. If <math>\mathbf{M}^{-1}</math> is symmetric positive-definite and <math>\mathbf{M}^{-1}\mathbf{A}</math> has a better condition number than <math>\mathbf{A},</math> a preconditioned conjugate gradient method can be used. It takes the following form:<ref>
{{cite book
| first1 = Richard
| last1 = Barrett
| first2 = Michael
| last2 = Berry
| first3 = Tony F.
| last3 = Chan
| first4 = James
| last4 = Demmel
| first5 = June
| last5 = Donato
| first6 = Jack
| last6 = Dongarra
| first7 = Victor
| last7 = Eijkhout
| first8 = Roldan
| last8 = Pozo
| first9 = Charles
| last9 = Romine
| first10 = Henk
| last10 = van der Vorst
| title = Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods
| edition = 2nd
| language = en
| publisher = SIAM
| page = 13
| url = http://www.netlib.org/templates/templates.pdf
| location = Philadelphia, PA
| access-date = 2020-03-31
}}</ref>

:<math>\mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0</math>
:<math> \textrm{Solve:}\mathbf{M}\mathbf{z}_0 := \mathbf{r}_0</math>
:<math>\mathbf{p}_0 := \mathbf{z}_0</math>
:<math>k := 0 \, </math>
:'''repeat'''
::<math>\alpha_k := \frac{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}{\mathbf{p}_k^\mathsf{T} \mathbf{A p}_k}</math>
::<math>\mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k</math>
::<math>\mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k</math>
::'''if''' '''r'''<sub>''k''+1</sub> is sufficiently small '''then''' exit loop '''end if'''
::<math>\mathrm{Solve}\ \mathbf{M}\mathbf{z}_{k+1} := \mathbf{r}_{k+1}</math>
::<math>\beta_k := \frac{\mathbf{r}_{k+1}^\mathsf{T} \mathbf{z}_{k+1}}{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}</math>
::<math>\mathbf{p}_{k+1} := \mathbf{z}_{k+1} + \beta_k \mathbf{p}_k</math>
::<math>k := k + 1 \, </math>
:'''end repeat'''
:The result is '''x'''<sub>''k''+1</sub>

The above formulation is equivalent to applying the regular conjugate gradient method to the preconditioned system<ref>{{cite book|first1=Gene H.|last1=Golub|first2= Charles F.|last2= Van Loan|title=Matrix Computations|edition=4th|at=sec. 11.5.2|publisher=Johns Hopkins University Press| isbn=978-1-4214-0794-4|date=2013}}</ref>
:<math>\mathbf{E}^{-1}\mathbf{A}(\mathbf{E}^{-1})^\mathsf{T}\mathbf{\hat{x}}=\mathbf{E}^{-1}\mathbf{b}</math>
where 
:<math>\mathbf{EE}^\mathsf{T}=\mathbf{M}, \qquad \mathbf{\hat{x}}=\mathbf{E}^\mathsf{T}\mathbf{x}.</math>

The Cholesky decomposition of the preconditioner must be used to keep the symmetry (and positive definiteness) of the system. However, this decomposition does not need to be computed, and it is sufficient to know <math>\mathbf{M}^{-1}</math>. It can be shown that <math>\mathbf{E}^{-1}\mathbf{A}(\mathbf{E}^{-1})^\mathsf{T}</math> has the same spectrum as <math>\mathbf{M}^{-1}\mathbf{A}</math>.

The preconditioner matrix '''M''' has to be symmetric positive-definite and fixed, i.e., cannot change from iteration to iteration. 
If any of these assumptions on the preconditioner is violated, the behavior of the preconditioned conjugate gradient method may become unpredictable.

An example of a commonly used [[preconditioner]] is the [[incomplete Cholesky factorization]].<ref>{{cite journal |first1=P. |last1=Concus |first2=G. H. |last2=Golub |first3=G. |last3=Meurant |year=1985 |title=Block Preconditioning for the Conjugate Gradient Method |journal=SIAM Journal on Scientific and Statistical Computing |volume=6 |issue=1 |pages=220–252 |doi=10.1137/0906018 |url=https://escholarship.org/uc/item/0j60b61v }}</ref>

=== Using the preconditioner in practice ===
It is important to keep in mind that we don't want to invert the matrix <math>\mathbf{M}</math> explicitly in order to get <math>\mathbf{M}^{-1}</math> for use it in the process, since inverting <math>\mathbf{M}</math> would take more time/computational resources than solving the conjugate gradient algorithm itself. As an example, let's say that we are using a preconditioner coming from incomplete Cholesky factorization. The resulting matrix is the lower triangular matrix <math>\mathbf{L}</math>, and the preconditioner matrix is:

<math>\mathbf{M}=\mathbf{LL}^\mathsf{T}</math>

Then we have to solve:

<math>\mathbf{Mz}=\mathbf{r}</math>

<math>\mathbf{z}=\mathbf{M}^{-1}\mathbf{r}</math>

But:

<math>\mathbf{M}^{-1}=(\mathbf{L}^{-1})^\mathsf{T}\mathbf{L}^{-1}</math>

Then:

<math>\mathbf{z}=(\mathbf{L}^{-1})^\mathsf{T}\mathbf{L}^{-1}\mathbf{r}</math>

Let's take an intermediary vector <math>\mathbf{a}</math>:

<math>\mathbf{a}=\mathbf{L}^{-1}\mathbf{r}</math>

<math>\mathbf{r}=\mathbf{L}\mathbf{a}</math>

Since <math>\mathbf{r}</math> and <math>\mathbf{L}</math> and known, and <math>\mathbf{L}</math> is lower triangular, solving for <math>\mathbf{a}</math> is easy and computationally cheap by using [[Triangular matrix#Forward and back substitution|forward substitution]]. Then, we substitute <math>\mathbf{a}</math> in the original equation:

<math>\mathbf{z}=(\mathbf{L}^{-1})^\mathsf{T}\mathbf{a}</math>

<math>\mathbf{a}=\mathbf{L}^\mathsf{T}\mathbf{z}</math>

Since <math>\mathbf{a}</math> and <math>\mathbf{L}^\mathsf{T}</math> are known, and <math>\mathbf{L}^\mathsf{T}</math> is upper triangular, solving for <math>\mathbf{z}</math> is easy and computationally cheap by using [[Triangular matrix#Forward and back substitution|backward substitution]].

Using this method, there is no need to invert <math>\mathbf{M}</math> or <math>\mathbf{L}</math> explicitly at all, and we still obtain <math>\mathbf{z}</math>.