Méthode du gradient conjugué

On rappelle que deux vecteurs non nuls u et v sont conjugués par rapport à A si

${\displaystyle u^{\mathrm {T} }\mathbf {A} v=0.}$

Sachant que A est symétrique définie positive, on en déduit un produit scalaire

${\displaystyle \langle u,v\rangle _{\mathbf {A} }:=\langle \mathbf {A} {u},{v}\rangle =\langle {u},\mathbf {A} ^{\mathrm {T} }{v}\rangle =\langle {u},\mathbf {A} {v}\rangle ={u}^{\mathrm {T} }\mathbf {A} {v}.}$

Ainsi, deux vecteurs sont conjugués s'ils sont orthogonaux pour ce produit scalaire.

La conjugaison est une relation symétrique : si u est conjugué à v pour A, alors v est conjugué à u.

Supposons que {p_k} est une suite de n directions conjuguées deux à deux. Alors les {p_k} forment une base de Rⁿ, ainsi la solution x_* de Ax= b dans cette base :

${\displaystyle x_{*}=\sum _{i=1}^{n}\alpha _{i}p_{i}}$

Les coefficients sont donnés par

${\displaystyle {b}=\mathbf {A} {x}_{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} {p}_{i}.}$

${\displaystyle {p}_{k}^{\mathrm {T} }{b}={p}_{k}^{\mathrm {T} }\mathbf {A} {x}_{*}=\sum _{i=1}^{n}\alpha _{i}{p}_{k}^{\mathrm {T} }\mathbf {A} {p}_{i}=\alpha _{k}{p}_{k}^{\mathrm {T} }\mathbf {A} {p}_{k}.}$ (car ${\displaystyle \forall i\neq k,p_{i},p_{k}}$ sont conjugués deux à deux)

${\displaystyle \alpha _{k}={\frac {{p}_{k}^{\mathrm {T} }{b}}{{p}_{k}^{\mathrm {T} }\mathbf {A} {p}_{k}}}={\frac {\langle {p}_{k},{b}\rangle }{\,\,\,\langle {p}_{k},{p}_{k}\rangle _{\mathbf {A} }}}={\frac {\langle {p}_{k},{b}\rangle }{\,\,\,\|{p}_{k}\|_{\mathbf {A} }^{2}}}.}$

On a ainsi l'idée directrice de la méthode pour résoudre le système Ax=b : trouver une suite de n directions conjuguées, et calculer les coefficients α_k.

En choisissant correctement les directions conjuguées p_k, il n'est pas nécessaire de toutes les déterminer pour obtenir une bonne approximation de la solution x_*. Il est ainsi possible de considérer la méthode du gradient conjugué comme une méthode itérative. Ce choix permet ainsi de considérer la résolution de systèmes de très grande taille, où le calcul de l'ensemble des directions aurait été très long.

On considère ainsi un premier vecteur x₀, qu'on pourra supposer nul (sinon, il faut considérer le système Az=b − Ax₀). L'algorithme va consister, partant de x₀, à se « rapprocher » de la solution x_* inconnue, ce qui suppose la définition d'une métrique. Cette métrique vient du fait que la solution x_* est l'unique minimiseur de la forme quadratique :

${\displaystyle f(\mathbf {x} )={\frac {1}{2}}x^{\mathrm {T} }\mathbf {A} x-x^{\mathrm {T} }b,\quad x\in \mathbb {R} ^{n}.}$

Ainsi, si f(x) diminue après une itération, alors on s'approche de x_*.

Ceci suggère donc de prendre la première direction p₁ comme l'opposé du gradient de f à x=x₀. Le gradient vaut Ax₀ – b= – b, d'après notre première hypothèse. Les vecteurs suivants de la base seront ainsi conjugués au gradient, d'où le nom « méthode du gradient conjugué ».

Soit r_k le résidu à la k^e itération :

${\displaystyle r_{k}=b-\mathbf {A} x_{k}.\,}$

Notons que r_k est l'opposé du gradient de f en x=x_k, ainsi, l'algorithme du gradient indique d'évoluer dans la direction r_k. On rappelle que les directions p_k sont conjuguées deux à deux. On veut aussi que la direction suivante soit construite à partir du résidu courant et des directions précédemment construites, ce qui est une hypothèse raisonnable en pratique.

La contrainte de conjugaison est une contrainte d'orthonormalité, aussi le problème partage des similitudes avec le procédé de Gram-Schmidt.

On a ainsi

${\displaystyle {p}_{k+1}={r}_{k}-\sum _{i\leq k}{\frac {{p}_{i}^{\mathrm {T} }\mathbf {A} {r}_{k}}{{p}_{i}^{\mathrm {T} }\mathbf {A} {p}_{i}}}{p}_{i}}$

Suivant cette direction, le point suivant est donné par

${\displaystyle {x}_{k+1}={x}_{k}+\alpha _{k+1}{p}_{k+1}}$

le pas α_k+1 est déterminé de manière à minimiser ${\displaystyle g(\alpha )=f({x}_{k}+\alpha {p}_{k+1})}$:

${\displaystyle g(\alpha )={\frac {1}{2}}\alpha ^{2}{p}_{k+1}^{\mathrm {T} }\mathbf {A} {p}_{k+1}+\alpha {p}_{k+1}^{\mathrm {T} }(\mathbf {A} {x}_{k}-b)+constante}$

le minimum de g est atteint pour ${\displaystyle {\frac {\mathrm {d} g}{\mathrm {d} \alpha }}(\alpha _{k+1})=0}$ et comme A est définie positive, ${\displaystyle {p}_{k+1}^{\mathrm {T} }\mathbf {A} {p}_{k+1}>0}$ ,

donc :

${\displaystyle \alpha _{k+1}={\frac {{p}_{k+1}^{\mathrm {T} }(b-\mathbf {A} {x}_{k})}{{p}_{k+1}^{\mathrm {T} }\mathbf {A} {p}_{k+1}}}={\frac {{p}_{k+1}^{\mathrm {T} }{r}_{k}}{{p}_{k+1}^{\mathrm {T} }\mathbf {A} {p}_{k+1}}}}$

Pour amorcer la récurrence, il faut partir d’une estimation initiale x₀ du vecteur x recherché ; et le nombre d'itérations N nécessaire pour que ${\displaystyle \|x_{N}-x\|<\varepsilon }$ (où ε est un nombre positif arbitrairement proche de zéro) dépend du x₀ choisi. Malheureusement, les méthodes de « préconditionnement » à la fois sûres et générales (c'est-à-dire efficaces pour toutes sortes de matrices symétriques positives) pour former un x₀ correct sont aussi elles-mêmes coûteuses en temps de calcul. En pratique, l'intuition physique, guidée par la nature physique du problème à résoudre, suggère parfois une initialisation efficace : ces idées ont donné lieu depuis plus de trente ans à une littérature spécialisée abondante[2].

L'algorithme ci-dessous résout Ax = b, où A est une matrice réelle, symétrique, et définie positive. Le vecteur d'entrée x₀ peut être une approximation de la solution initiale ou 0.

${\displaystyle {\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\hbox{répéter}}\\&\qquad \alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}\\&\qquad \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}\\&\qquad \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}\\&\qquad {\hbox{si }}r_{k+1}{\hbox{ est suffisamment petit, alors on sort de la boucle}}\\&\qquad \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}\\&\qquad \mathbf {p} _{k+1}:=\mathbf {r} _{k+1}+\beta _{k}\mathbf {p} _{k}\\&\qquad k:=k+1\\&{\hbox{fin de répéter}}\\&{\hbox{Le résultat est }}\mathbf {x} _{k+1}\end{aligned}}}$

On peut montrer le résultat suivant sur la convergence de l'algorithme :

${\displaystyle \|x^{*}-x_{k}\|_{\mathbf {A} }\leqslant 2\left({\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\|x^{*}-x_{0}\|_{\mathbf {A} },}$

où ${\displaystyle {\sqrt {\kappa (\mathbf {A} )}}}$ désigne le conditionnement de la matrice et ${\displaystyle \|z\|_{\mathbf {A} }={\sqrt {z^{T}\mathbf {A} z}}.}$

La méthode du gradient conjugué a donc une convergence superlinéaire, qui peut être mise à mal par un mauvais conditionnement de la matrice. Elle reste toutefois meilleure que les algorithmes à direction de plus forte pente.

• Philippe Ciarlet, Introduction à l’analyse numérique matricielle et à l’optimisation, Masson, coll. « Math. Appl. pour la Maîtrise », 1985 (réimpr. 2001) (ISBN 2-225-68893-1)
• Dianne P. O'Leary (dir.), Linear and nonlinear Conjugate gradient-related Methods, AMS-SIAM, 1996, « Conjugate gradient and related KMP algorithms : the Beginnings »

• Direction de descente
• Algorithme à directions de descente
• Algorithme de Lanczos

• (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Conjugate gradient method » (voir la liste des auteurs).
• (fr) La méthode du gradient conjugué
• (en) Le gradient conjugué sans peine
• (fr) Introduction à la méthode du gradient et du gradient conjugué