Nous allons étudier une méthode directe de résolution de système linéaire: la méthode de Cholesky. L’objectif est de construire pour une matrice A symétrique définie positive une matrice triangulaire L telle que A soit égale au produit de la matrice triangulaire L et de la transposée de cette dernière.

Description du problème

Considérons le système suivant:

\[Ax=b,\]

où $A$ est une matrice de taille $n\times n$ symétrique définie positive ($A^\top =A$ et $x^\top Ax>0$, pour tout vecteur $x\in \mathbf{R}^n$ non nul). Soit $x_\star$ la solution exacte de ce système.

Directions conjuguées

Comme la matrice $A$ est symétrique définie positive, on peut définir le produit scalaire suivant sur $\mathbf{R}^n$: $\langle u,v \rangle_A = u^\top A v.$ Deux éléments $u,v\in \mathbf{R}^n$ sont dit $A$-conjugués si:

\[u^\top A v = 0\]

La méthode du gradient conjugué consiste à  construire une suite $(p_k)_{k\in\mathbf{N}^*}$ de $n$ vecteurs $A$-conjugués. Dès lors, la suite $p_1,p_2,\ldots,p_n$ forme une base de $\mathbf{R}^n$. La solution exacte $x_\star$ peut donc se décomposer comme suit:

\[x_\star = \alpha_1 p_1 + \cdots + \alpha_n p_n\]

où $\alpha_k = \dfrac{p_k^\top b}{p_k^\top A p_k},k=1,\ldots,n.$

Construction des directions conjuguées

La solution exacte $x_\star$ peut-àªtre également vu comme l’unique minimisant de la fonctionnelle

\[J(x) = \frac12 x^\top A x - b^\top x, \quad x\in\mathbf{R}^n\]

On a donc clairement

\(\nabla J(x) = A x - b, \quad x\in\mathbf{R}^n\) d’où \(\nabla J(x_\star)=0_{\mathbf{R}^n}\)

On définit le résidu du système d’équation comme suit

\[r_k = b - Ax_k=-\nabla J(x_k)\]

$r_k$ représente donc la direction du gradient de la fonctionnelle $J$ en $x_k$ (à  un signe près).

La nouvelle direction de descente $p_{k+1}$ suit donc celle du résidu modulo, sa $A$-conjugaison avec $p_k$, on a alors:

\[p_{k+1} = r_k - \frac{p_k^\top A r_k}{p_k^\top A p_k} p_k\]

C’est le choix du coefficient $\dfrac{p_k^\top A r_k}{p_k^\top A p_k}$ qui assure la $A$-conjugaison des directions $p_k$. Pour vous en assurez calculer $(Ap_{k+1},p_k)$, cette quantité est nulle !

Algorithme du gradient conjugué

Calcul du résidu $r_0=b-Ax_0$ pour un vecteur $x_0$ quelconque. On fixe $p_0=r_0$.

Pour $k=1,2,\ldots$

\[q_k=A p_k\]

If $k=1$

\[p_1=r_{0}\]

Else

\[\beta_{k-1}= \frac{r_{k-1}^\top r_{k-1}}{r_{k-2}^\top r_{k-2}}\] \[p_k=r_{k-1} + \beta_{k-1}p_{k-1}\]

EndIf

\[\alpha_k=\frac{r_{k-1}^\top r_{k-1}}{p_k^\top q_k}\] \[x_k=x_{k-1}+\alpha_k p_k\] \[r_k=r_{k-1}-\alpha_k q_k\]

FinPour