Méthode du gradient conjugué

En analyse numérique, la méthode du gradient conjugué est un algorithme pour résoudre des systèmes d'équations linéaires dont la matrice est symétrique définie positive. Cette méthode, imaginée en 1950 simultanément par Cornelius Lanczos, Eduard Stiefel et Magnus Hestenes^[1], est une méthode itérative qui converge en un nombre fini d'itérations (au plus égal à la dimension du système linéaire). Toutefois, son grand intérêt pratique du point de vue du temps de calcul vient de ce qu’une initialisation astucieuse (dite « préconditionnement ») permet d'aboutir en seulement quelques passages à une estimation très proche de la solution exacte : c'est pourquoi, en pratique, on se borne à un nombre d'itérations bien inférieur au nombre d'inconnues.

La méthode du gradient biconjugué fournit une généralisation pour les matrices non symétriques.

Principe[modifier | modifier le code]

L'objectif est de minimiser la fonction $f:x\mapsto {\frac {1}{2}}(\mathbf {A} x,x)-(b,x)$ où $A$ est une matrice carrée symétrique définie positive de taille n.

Le calcul montre qu'une solution du problème est la solution du système $\mathbf {A} x=b$ : en effet, on a $\nabla f\left(x\right)=\mathbf {A} x-b$ .

Intuitivement, la fonction $f$ peut donc être vue comme une primitive (littéralement un potentiel scalaire) du résidu $\mathbf {A} x-b$ . En annulant le gradient de $f$ , on obtient le vecteur $x$ qui minimise l'erreur.

La méthode du gradient conjugué vue comme une méthode directe[modifier | modifier le code]

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

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

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

\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 $A x = b$ dans cette base :

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

Les coefficients sont donnés par

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

{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

\forall i\neq k,p_{i},p_{k}

sont conjugués deux à deux)

\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 $A x = b$ : trouver une suite de $n$ directions conjuguées, et calculer les coefficients $α k$ .

La méthode du gradient conjugué vue comme une méthode itérative[modifier | modifier le code]

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 0$ , qu'on pourra supposer nul (sinon, il faut considérer le système $A z = b - A x 0$ ). L'algorithme va consister, partant de $x 0$ , à 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 :

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 1$ comme l'opposé du gradient de $f$ à $x = x 0$ . Le gradient vaut $A x 0 - 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 :

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

{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

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

le pas

α k +1

est déterminé de manière à minimiser

g(\alpha )=f({x}_{k}+\alpha {p}_{k+1})

:

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

{\frac {\mathrm {d} g}{\mathrm {d} \alpha }}(\alpha _{k+1})=0

et comme

A

est définie positive,

{p}_{k+1}^{\mathrm {T} }\mathbf {A} {p}_{k+1}>0

,

donc :

\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}}}

Algorithme[modifier | modifier le code]

Pour amorcer la récurrence, il faut partir d’une estimation initiale $x 0$ du vecteur $x$ recherché ; et le nombre d'itérations N nécessaire pour que $\|x_{N}-x\|<\varepsilon$ (où ε est un nombre positif arbitrairement proche de zéro) dépend du $x 0$ 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 0$ 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].

Algorithme itératif en pseudo-code[modifier | modifier le code]

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

{\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}}

Convergence[modifier | modifier le code]

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

\|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ù ${\sqrt {\kappa (\mathbf {A} )}}$ désigne le conditionnement de la matrice et $\|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.

Solveur[modifier | modifier le code]

(en) M1CG1 - A solver of symmetric linear systems by conjugate gradient iterations, using/building a BFGS/ℓ-BFGS preconditioner. Écrit en Fortran-77. Le solveur a l'intérêt d'offrir la possibilité de construire un préconditionneur BFGS ou ℓ-BFGS (en), qui pourra être utile pour la résolution d'un système linéaire avec une matrice proche et un second membre différent.

Notes et références[modifier | modifier le code]

Notes[modifier | modifier le code]

↑ Magnus Hestenes et Eduard Stiefel, « Methods of Conjugate Gradients for Solving Linear Systems », Journal of Research of the National Bureau of Standards, vol. 49, n^o 6,‎ 1952 (lire en ligne [PDF])
↑ Selon Dianne O'Leary (cf. bibliographie), l'article de J. Meijerink et H. van der Vorst, « An iterative solution method for linear systems of which the coefficient matrix is symmetric a M-matrix », Mathematics of Computation, n^o 31,‎ 1977, p. 148-162 marque une étape décisive dans l'idée de préconditionner un système linéaire avant de lui appliquer l'algorithme. Cet article pionnier proposait le pré-conditionnement par Décomposition LU incomplète. Suivirent entre autres le pré-conditionnement par SOR interrompu (M. DeLong et J. Ortega, « SOR as a preconditionner », Applied Numerical Mathematics, n^o 18,‎ 1995, p. 431-440), et par Méthode de Gauss-Seidel interrompue (Y. Saad et Gene Golub (dir.), Parallel preconditionners for general sparse matrices, Recent Advances in Iterative Methods, Springer Verlag, 1994, 165-199 p.). On trouvera un aperçu des différentes techniques dans, entre autres : A. Bruaset, A survey of preconditionned iterative methods, Longman Scientific & Technical, coll. « Pitman Research Notes in Mathematics », 1995 ; J. Erhel et K. Burrage, On the performance of various adaptive preconditioned GMRES strategies, INRIA/IRISA, 1997, etc.

Bibliographie[modifier | modifier le code]

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 »