Méthode du gradient conjugué

Illustration de la 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, 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

L'objectif est de minimiser la fonction $f:x\mapsto \frac{1}{2}(𝐀x,x)-(b,x)$ où Modèle:Math 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 $𝐀x=b$ : en effet, on a $\mathrm{\nabla }f\left(x\right)=𝐀x-b$.

Intuitivement, la fonction Modèle:Mvar peut donc être vue comme une primitive (littéralement un potentiel scalaire) du résidu $𝐀x-b$. En annulant le gradient de Modèle:Mvar, on obtient le vecteur Modèle:Mvar qui minimise l'erreur.

La méthode du gradient conjugué vue comme une méthode directe

On rappelle que deux vecteurs non nuls Modèle:Mvar et Modèle:Mvar sont conjugués par rapport à Modèle:Math si

$u^{\mathrm{T}}𝐀v=0.$

Sachant que Modèle:Math est symétrique définie positive, on en déduit un produit scalaire

$⟨u,v{⟩}_{𝐀}:=⟨𝐀u,v⟩=⟨u,{𝐀}^{\mathrm{T}}v⟩=⟨u,𝐀v⟩=u^{\mathrm{T}}𝐀v.$

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

La conjugaison est une relation symétrique : si Modèle:Mvar est conjugué à Modèle:Mvar pour Modèle:Math, alors Modèle:Mvar est conjugué à Modèle:Mvar.

Supposons que Modèle:Formule est une suite de Modèle:Mvar directions conjuguées deux à deux. Alors les Modèle:Formule forment une base de Rⁿ, ainsi la solution Modèle:Math de Modèle:Formule dans cette base :

$x_{*}=\sum _{i=1}^n{\alpha }_i{p}_i$

Les coefficients sont donnés par

$b=𝐀x_{*}=\sum _{i=1}^n{\alpha }_i𝐀p_i.$

$p_k^{\mathrm{T}}b=p_k^{\mathrm{T}}𝐀x_{*}=\sum _{i=1}^n{\alpha }_i{p}_k^{\mathrm{T}}𝐀p_i={\alpha }_k{p}_k^{\mathrm{T}}𝐀p_k.$ (car $\mathrm{\forall }i\ne k,p_i,p_k$ sont conjugués deux à deux)

${\alpha }_k=\frac{p_k^{\mathrm{T}}b}{p_k^{\mathrm{T}}𝐀p_k}=\frac{⟨p_k,b⟩}{⟨p_k,p_k{⟩}_{𝐀}}=\frac{⟨p_k,b⟩}{\Vert p_k{\Vert }_{𝐀}^2}.$

On a ainsi l'idée directrice de la méthode pour résoudre le système Modèle:Formule : trouver une suite de Modèle:Mvar directions conjuguées, et calculer les coefficients Modèle:Mvar.

La méthode du gradient conjugué vue comme une méthode itérative

En choisissant correctement les directions conjuguées Modèle:Mvar, il n'est pas nécessaire de toutes les déterminer pour obtenir une bonne approximation de la solution Modèle:Math. 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 Modèle:Math, qu'on pourra supposer nul (sinon, il faut considérer le système Modèle:Formule). L'algorithme va consister, partant de Modèle:Math, à se « rapprocher » de la solution Modèle:Math inconnue, ce qui suppose la définition d'une métrique. Cette métrique vient du fait que la solution Modèle:Math est l'unique minimiseur de la forme quadratique :

$f(𝐱)=\frac{1}{2}{x}^{\mathrm{T}}𝐀x-x^{\mathrm{T}}b,x\in {ℝ}^n.$

Ainsi, si Modèle:Formule diminue après une itération, alors on s'approche de Modèle:Math.

Ceci suggère donc de prendre la première direction Modèle:Math comme l'opposé du gradient de Modèle:Mvar à Modèle:Formule. Le gradient vaut Modèle:Formule, 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 Modèle:Mvar le résidu à la kModèle:E itération :

$r_k=b-𝐀x_k.$

Notons que Modèle:Mvar est l'opposé du gradient de Modèle:Mvar en Modèle:Formule, ainsi, l'algorithme du gradient indique d'évoluer dans la direction Modèle:Mvar. On rappelle que les directions Modèle:Mvar 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=r_k-\sum _{i\le k-1}\frac{p_i^{\mathrm{T}}𝐀r_k}{p_i^{\mathrm{T}}𝐀p_i}{p}_i$

Suivant cette direction, le point suivant est donné par

$x_{k+1}=x_k+{\alpha }_k{p}_k$

où le pas ${\alpha }_k$ est déterminé de manière à minimiser $g(\alpha )=f(x_k+\alpha p_k)=\frac{1}{2}{\alpha }^2{p}_k^{\mathrm{T}}𝐀p_k+\alpha p_k^{\mathrm{T}}(𝐀x_k-b)+constante$;

le minimum de Modèle:Mvar est atteint pour $\frac{\mathrm{d}g}{\mathrm{d}\alpha }({\alpha }_k)=0$ et comme Modèle:Math est définie positive, $p_k^{\mathrm{T}}𝐀p_k>0$ ,

on a donc : ${\alpha }_k=\frac{p_k^{\mathrm{T}}(b-𝐀x_k)}{p_k^{\mathrm{T}}𝐀p_k}=\frac{p_k^{\mathrm{T}}{r}_k}{p_k^{\mathrm{T}}𝐀p_k}$

Une analyse plus détaillée de cette algorithme (un raisonnement par récurrence) montre que ${𝐫}_i$ est orthogonal à ${𝐫}_j$, i.e. ${𝐫}_i^{𝖳}{𝐫}_j=0$ pour $i\ne j$ (voir ci-après) et que ${𝐩}_i$ est $𝐀$-orthogonal à ${𝐩}_j$ , i.e. ${𝐩}_i^{𝖳}𝐀{𝐩}_j=0$ , pour $i\ne j$.

Algorithme

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

Algorithme itératif en pseudo-code

L'algorithme ci-dessous résout Modèle:Math, où Modèle:Math est une matrice réelle, symétrique, et définie positive. Le vecteur d'entrée Modèle:Math peut être une approximation de la solution initiale ou 0. Cette algorithme est issu de la méthode itérative exacte présentée dans le paragraphe précédent, les valeurs des coefficients ${\alpha }_k$ semblent différentes mais en utilisant les relations de l'algorithme ci-dessous, en particulier : ${𝐩}_k={𝐫}_k+{\beta }_{k-1}{𝐩}_{k-1}$, et le fait que les résidus ${𝐫}_k$ et ${𝐫}_{k-1}$ soient orthogonaux , on peut montrer par récurrence que l'on a bien : ${\alpha }_k=\frac{p_k^{\mathrm{T}}{r}_k}{p_k^{\mathrm{T}}𝐀p_k}=\frac{r_k^{\mathrm{T}}{r}_k}{p_k^{\mathrm{T}}𝐀p_k}$. Ci-dessous le coefficient ${\alpha }_k$ est choisi pour que les résidus ${𝐫}_{k+1}$ et ${𝐫}_k$ soient orthogonaux et non pas minimiser g comme dans le paragraphe précédent, mais en fait les deux approches reviennent à la même formule pour ${\alpha }_k$ et le coefficient ${\beta }_k$ est choisi pour que ${𝐩}_{k+1}$ soit A-conjugué de ${𝐩}_k$.

$\begin{array}{cc}& {𝐫}_0:=𝐛-{𝐀𝐱}_0\\ & {𝐩}_0:={𝐫}_0\\ & k:=0\\ & \text{répéter}\\ & {\alpha }_k:=\frac{{𝐫}_k^{𝖳}{𝐫}_k}{{𝐩}_k^{𝖳}{𝐀𝐩}_k}\\ & {𝐱}_{k+1}:={𝐱}_k+{\alpha }_k{𝐩}_k\\ & {𝐫}_{k+1}:={𝐫}_k-{\alpha }_k{𝐀𝐩}_k\\ & \text{si }{r}_{k+1}\text{ est suffisamment petit, alors on sort de la boucle}\\ & {\beta }_k:=\frac{{𝐫}_{k+1}^{𝖳}{𝐫}_{k+1}}{{𝐫}_k^{𝖳}{𝐫}_k}\\ & {𝐩}_{k+1}:={𝐫}_{k+1}+{\beta }_k{𝐩}_k\\ & k:=k+1\\ & \text{fin de répéter}\\ & \text{Le résultat est }{𝐱}_{k+1}\end{array}$

Cet algorithme du gradient conjugué est celui qui est très souvent utilisé. Il est intéressant de constater que le coefficient ${\beta }_k$ de cet algorithme a la même expression que celui qui est utilisé dans la méthode 'Fletcher-Reeves' du gradient conjugué non-linéaire.

On peut également remarquer que $x_1$ est déduit de $x_0$ en utilisant la méthode du gradient et que prendre ${\beta }_k=0$, revient à appliquer la méthode du gradient et peut donc être utilisé pour réinitialiser un calcul de gradient conjugué en cours. Réinitialiser un calcul en cours peut ralentir la convergence, mais peu également en augmenter la stabilité en particulier en réduisant les erreurs dus à l'accumulation d'imprécisions numériques (round-off errors).

Les formules ${𝐱}_{k+1}:={𝐱}_k+{\alpha }_k{𝐩}_k$ et ${𝐫}_k:=𝐛-{𝐀𝐱}_k$, qui sont exactes, impliquent que les deux formules ${𝐫}_{k+1}:={𝐫}_k-{\alpha }_k{𝐀𝐩}_k$ et ${𝐫}_{k+1}:=𝐛-{𝐀𝐱}_{k+1}$ sont mathématiquement équivalentes. La première formule ${𝐫}_{k+1}:={𝐫}_k-{\alpha }_k{𝐀𝐩}_k$ est utilisé dans l'algorithme pour éviter une multiplication supplémentaire par $𝐀$ car le vector ${𝐀𝐩}_k$ est déjà calculé pour évaluer ${\alpha }_k$. La deuxième formule ${𝐫}_{k+1}:=𝐛-{𝐀𝐱}_{k+1}$ peut par contre être plus précise car elle réduit l'accumulation des imprécisions numérique et elle est donc parfois recommandée.

Exemple numérique

Considérons le système linéaire Ax = b suivant :

$𝐀𝐱=\left[\begin{array}{cc}4& 1\\ 1& 3\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}1\\ 2\end{array}\right],$

Deux itérations de l'algorithme du gradient conjugué vont être réalisées pas à pas à en partant du vecteur initial

${𝐱}_0=\left[\begin{array}{c}2\\ 1\end{array}\right]$

Il est rappelé, que s'il n'y avait pas d'imprécision numérique dans les calculs le présent exemple est résolu en seulement n = 2 itérations, puisque ici la matrice est de dimension 2x2 et donc que x₂ devrait aux erreurs numériques près retrouver la solution exacte.

Solution

La solution exacte de ce système linéaire est obtenu en inversant simplement la matrice A :

${𝐀}^{-1}=\left[\begin{array}{cc}3& -1\\ -1& 4\end{array}\right]/11$ et $𝐱={𝐀}^{-1}\left[\begin{array}{c}1\\ \\ 2\end{array}\right]=\left[\begin{array}{c}\frac{1}{11}\\ \\ \frac{7}{11}\end{array}\right]\approx \left[\begin{array}{c}0.0909\\ \\ 0.6364\end{array}\right]$

Les valeurs initiales de l'algorithme itératifs sont :

${𝐫}_0=𝐛-𝐀{𝐱}_0=\left[\begin{array}{c}1\\ 2\end{array}\right]-\left[\begin{array}{cc}4& 1\\ 1& 3\end{array}\right]\left[\begin{array}{c}2\\ 1\end{array}\right]=\left[\begin{array}{c}-8\\ -3\end{array}\right]={𝐩}_0.$

${\alpha }_0=\frac{{𝐫}_0^{𝖳}{𝐫}_0}{{𝐩}_0^{𝖳}{𝐀𝐩}_0}=\frac{\left[\begin{array}{cc}-8& -3\end{array}\right]\left[\begin{array}{c}-8\\ -3\end{array}\right]}{\left[\begin{array}{cc}-8& -3\end{array}\right]\left[\begin{array}{cc}4& 1\\ 1& 3\end{array}\right]\left[\begin{array}{c}-8\\ -3\end{array}\right]}=\frac{73}{331}\approx 0.2205$

Lors de la première itération la valeur x₁ est une meilleur approximation de la solution exacte:

${𝐱}_1={𝐱}_0+{\alpha }_0{𝐩}_0=\left[\begin{array}{c}2\\ 1\end{array}\right]+\frac{73}{331}\left[\begin{array}{c}-8\\ -3\end{array}\right]\approx \left[\begin{array}{c}0.2356\\ 0.3384\end{array}\right].$

et le résidu associé r₁ est égal à :

${𝐫}_1={𝐫}_0-{\alpha }_0𝐀{𝐩}_0=\left[\begin{array}{c}-8\\ -3\end{array}\right]-\frac{73}{331}\left[\begin{array}{cc}4& 1\\ 1& 3\end{array}\right]\left[\begin{array}{c}-8\\ -3\end{array}\right]\approx \left[\begin{array}{c}-0.2810\\ 0.7492\end{array}\right].$

Ce résidu est encore grand et donc il est nécessaire de poursuivre la deuxième itération et donc de calculer les paramètres qui vont être utile pour calculer x₂: (i) d'abord Modèle:Math puis (ii) la nouvelle direction de recherche p₁ et ensuite (iii) Modèle:Math :

${\beta }_0=\frac{{𝐫}_1^{𝖳}{𝐫}_1}{{𝐫}_0^{𝖳}{𝐫}_0}\approx \frac{\left[\begin{array}{cc}-0.2810& 0.7492\end{array}\right]\left[\begin{array}{c}-0.2810\\ 0.7492\end{array}\right]}{\left[\begin{array}{cc}-8& -3\end{array}\right]\left[\begin{array}{c}-8\\ -3\end{array}\right]}=0.0088.$

${𝐩}_1={𝐫}_1+{\beta }_0{𝐩}_0\approx \left[\begin{array}{c}-0.2810\\ 0.7492\end{array}\right]+0.0088\left[\begin{array}{c}-8\\ -3\end{array}\right]=\left[\begin{array}{c}-0.3511\\ 0.7229\end{array}\right].$

${\alpha }_1=\frac{{𝐫}_1^{𝖳}{𝐫}_1}{{𝐩}_1^{𝖳}{𝐀𝐩}_1}\approx \frac{\left[\begin{array}{cc}-0.2810& 0.7492\end{array}\right]\left[\begin{array}{c}-0.2810\\ 0.7492\end{array}\right]}{\left[\begin{array}{cc}-0.3511& 0.7229\end{array}\right]\left[\begin{array}{cc}4& 1\\ 1& 3\end{array}\right]\left[\begin{array}{c}-0.3511\\ 0.7229\end{array}\right]}=0.4122.$

${𝐱}_2={𝐱}_1+{\alpha }_1{𝐩}_1\approx \left[\begin{array}{c}0.2356\\ 0.3384\end{array}\right]+0.4122\left[\begin{array}{c}-0.3511\\ 0.7229\end{array}\right]=\left[\begin{array}{c}0.0909\\ 0.6364\end{array}\right].$

x₂ est une bien meilleur approximation de la solution que x₁ et x₀ et cela est numériquement déterminé en évaluant le résidu r₂.

${𝐫}_2={𝐫}_1-{\alpha }_1𝐀{𝐩}_1=\left[\begin{array}{c}-0.2810\\ 0.7492\end{array}\right]-0.4122\left[\begin{array}{cc}4& 1\\ 1& 3\end{array}\right]\left[\begin{array}{c}-0.3511\\ 0.7229\end{array}\right]\approx \left[\begin{array}{c}-0.00009\\ -0.00001\end{array}\right].$

Il est également intéressant de vérifier que : ${𝐫}_1^{𝖳}{𝐫}_1={𝐩}_1^{𝖳}{𝐫}_1$ :

${𝐫}_1^{𝖳}{𝐫}_1=\left[\begin{array}{cc}-0.2810& 0.7492\end{array}\right]\left[\begin{array}{c}-0.2810\\ 0.7492\end{array}\right]=0.6402$ et

${𝐩}_1^{𝖳}{𝐫}_1=\left[\begin{array}{cc}-0.3511& 0.7229\end{array}\right]\left[\begin{array}{c}-0.2810\\ 0.7492\end{array}\right]=0.6402$

Convergence

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

$\Vert x^{*}-x_k{\Vert }_{𝐀}⩽2{\left(\frac{\sqrt{\kappa (𝐀)}-1}{\sqrt{\kappa (𝐀)}+1}\right)}^k\Vert x^{*}-x_0{\Vert }_{𝐀},$

où $\sqrt{\kappa (𝐀)}$ désigne le conditionnement de la matrice et $\Vert z{\Vert }_{𝐀}=\sqrt{z^T𝐀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

• Modèle: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 Modèle:Lien, 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

Bibliographie

• Modèle:Ouvrage
• Modèle:Ouvrage
• Modèle:Ouvrage

Voir aussi

Articles connexes

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

Liens externes

• Modèle:Fr La méthode du gradient conjugué
• Modèle:En Le gradient conjugué sans peine
• Modèle:Fr Introduction à la méthode du gradient et du gradient conjugué