Méthode du gradient biconjugué

1. Choisir ${\displaystyle x_{0}}$, ${\displaystyle y_{0}}$, un préconditionneur régulier ${\displaystyle M}$ (on utilise fréquemment ${\displaystyle M^{-1}=1}$) et ${\displaystyle c}$;
2. ${\displaystyle r_{0}\leftarrow b-Ax_{0},s_{0}\leftarrow c-A^{*}y_{0}\,}$;
3. ${\displaystyle d_{0}\leftarrow M^{-1}r_{0},f_{0}\leftarrow M^{-*}s_{0}\,}$;
4. for ${\displaystyle k=0,1,\dots \,}$ do
5. ${\displaystyle \alpha _{k}\leftarrow {s_{k}^{*}M^{-1}r_{k} \over f_{k}^{*}Ad_{k}}\,}$;
6. ${\displaystyle x_{k+1}\leftarrow x_{k}+\alpha _{k}d_{k}\,}$ ${\displaystyle \left(y_{k+1}\leftarrow y_{k}+{\overline {\alpha _{k}}}f_{k}\right)\,}$;
7. ${\displaystyle r_{k+1}\leftarrow r_{k}-\alpha _{k}Ad_{k}\,}$, ${\displaystyle s_{k+1}\leftarrow s_{k}-{\overline {\alpha _{k}}}A^{*}f_{k}\,}$ (${\displaystyle r_{k}=b-Ax_{k}}$ et ${\displaystyle s_{k}=c-A^{*}y_{k}}$ sont le résidus);
8. ${\displaystyle \beta _{k}\leftarrow {s_{k+1}^{*}M^{-1}r_{k+1} \over s_{k}^{*}M^{-1}r_{k}}\,}$;
9. ${\displaystyle d_{k+1}\leftarrow M^{-1}r_{k+1}+\beta _{k}d_{k}\,}$, ${\displaystyle f_{k+1}\leftarrow M^{-*}s_{k+1}+{\overline {\beta _{k}}}f_{k}\,}$.

La méthode est numériquement instable, mais on y remédie par la méthode stabilisée du gradient biconjugué (en), et elle reste très importante du point de vue théorique : on définit l'itération par ${\displaystyle x_{k}:=x_{j}+P_{k}A^{-1}\left(b-Ax_{j}\right)}$ et ${\displaystyle y_{k}:=y_{j}+A^{-*}P_{k}^{*}\left(c-A^{*}y_{j}\right)}$ (${\displaystyle j<k}$) en utilisant les projections suivantes :

${\displaystyle P_{k}:=\mathbf {u_{k}} \left(\mathbf {v_{k}} ^{*}A\mathbf {u_{k}} \right)^{-1}\mathbf {v_{k}} ^{*}A}$,

Avec ${\displaystyle \mathbf {u_{k}} =\left(u_{0},u_{1},\dots u_{k-1}\right)}$ et ${\displaystyle \mathbf {v_{k}} =\left(v_{0},v_{1},\dots v_{k-1}\right)}$. On peut itérer les projections elles-mêmes, comme

${\displaystyle P_{k+1}=P_{k}+\left(1-P_{k}\right)u_{k}\otimes {v_{k}^{*}A\left(1-P_{k}\right) \over v_{k}^{*}A\left(1-P_{k}\right)u_{k}}}$.

Les nouvelles directions de descente ${\displaystyle d_{k}:=\left(1-P_{k}\right)u_{k}}$ et ${\displaystyle f_{k}:=\left(A\left(1-P_{k}\right)A^{-1}\right)^{*}v_{k}}$ sont alors orthogonales aux résidus : ${\displaystyle v_{i}^{*}r_{k}=f_{i}^{*}r_{k}=0}$ et ${\displaystyle s_{k}^{*}u_{j}=s_{k}^{*}d_{j}=0}$, qui satisfont aux mêmes ${\displaystyle r_{k}=A\left(1-P_{k}\right)A^{-1}r_{j}}$ et ${\displaystyle s_{k}=\left(1-P_{k}\right)^{*}s_{j}}$(${\displaystyle i,j<k}$).

La méthode du gradient biconjugué propose alors le choix suivant :

${\displaystyle u_{k}:=M^{-1}r_{k}}$ et ${\displaystyle v_{k}:=M^{-*}s_{k}}$.

Ce choix particulier permet alors d'éviter une évaluation directe de ${\displaystyle P_{k}}$ et ${\displaystyle A^{-1}}$, et donc augmenter la vitesse d'exécution de l'algorithme.

• Si ${\displaystyle A=A^{*}}$ est auto-adjointe, ${\displaystyle y_{0}=x_{0}}$ et ${\displaystyle c=b}$, donc ${\displaystyle r_{k}=s_{k}}$, ${\displaystyle d_{k}=f_{k}}$, et la méthode du gradient conjugué produit la même suite ${\displaystyle x_{k}=y_{k}}$.

• En dimensions finies ${\displaystyle x_{n}=A^{-1}b}$, au plus tard quand ${\displaystyle P_{n}=1}$ : La méthode du gradient biconjugué rend la solution exacte après avoir parcouru tout l'espace et est donc une méthode directe.

• La suite produite par l'algorithme est biorthogonale (en) : ${\displaystyle f_{i}^{*}Ad_{j}=0}$ et ${\displaystyle s_{i}^{*}M^{-1}r_{j}=0}$ pour ${\displaystyle i\neq j}$.

• SI ${\displaystyle p_{j'}}$ est un polynôme avec ${\displaystyle deg\left(p_{j'}\right)+j<k}$, alors ${\displaystyle s_{k}^{*}p_{j'}\left(M^{-1}A\right)u_{j}=0}$. L'algorithme est donc composé de projections sur des sous-espaces de Krylov ;

• SI ${\displaystyle p_{i'}}$ est un polynôme avec ${\displaystyle i+deg\left(p_{i'}\right)<k}$, alors ${\displaystyle v_{i}^{*}p_{i'}\left(AM^{-1}\right)r_{k}=0}$.

• (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Biconjugate gradient method » (voir la liste des auteurs).