Sup-Galilée – MACS 2

Analyse numérique avancée

Michel KERN 1

2019–2020

1. Inria, Centre de Paris, 75012 Paris, [email protected]

2
Table des matières

1 La méthode du gradient conjugué 6

1.1 Notations – généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Définition de la méthode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 L’algorithme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.1 Développement de l’algorithme . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.2 Mise en oeuvre de l’algorithme . . . . . . . . . . . . . . . . . . . . . . . . 13
1.4 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2 Systèmes linéaires non symétriques : la méthode GMRES 20

2.1 Définition de la méthode GMRES, et premières propriétés . . . . . . . . . . . . . 20
2.2 L’algorithme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2.1 L’algorithme d’Arnoldi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2.2 Résolution du problème de moindres carrés . . . . . . . . . . . . . . . . . 24
2.2.3 Conséquences de l’arithmétique flottante . . . . . . . . . . . . . . . . . . . 25
2.3 Remarques sur la convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3 Préconditionnement 29
3.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.1.1 À gauche, à droite, ou des deux cotés ? . . . . . . . . . . . . . . . . . . . . 30
3.2 Préconditionnement du gradient conjugué . . . . . . . . . . . . . . . . . . . . . . 31
3.3 Exemples de préconditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3.1 Préconditionnement diagonal . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3.2 Préconditionnment par des méthodes itératives classiques . . . . . . . . . 34
3.3.3 Factorisations incomplètes . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.3.4 Préconditionnement polynomial . . . . . . . . . . . . . . . . . . . . . . . . 39

4 Matrices creuses 42
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.2 Stockage par diagonales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.2 Algorithmes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.3 Stockage CSR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3.1 Stockage « coordonnées » . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3.2 Le stockage CSR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3
A Rappels d’algèbre linéaire 50
A.1 Matrices symétriques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
A.2 Normes vectoriellles et matricielles . . . . . . . . . . . . . . . . . . . . . . . . . . 50
A.3 Conditionnement d’un système linéaire . . . . . . . . . . . . . . . . . . . . . . . . 52
A.4 Problèmes de moindres carrés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
A.4.1 Propriétés mathématiques . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
A.4.2 Méthodes numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
A.5 Méthodes itératives classiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4
 Avertissement

Cette édition du cours est encore un travail non–achevé. Il est possible, et malheureusement
probable, qu’il reste des erreurs, des fautes de frappe, des explications peu claire, et même peut-
être des fautes de mathématiques. Par conséquent, lisez ce cours avec attention et un esprit
critique.
Dans l’espoir d’améliorer ce document pour les promotions à venir, je vous sera gré de me
communiquer les erreurs que vous pourriez découvrir, par un message à [email protected]

5
Chapitre 1

La méthode du gradient conjugué

Nous présentons dans ce chapitre la méthode gradient conjugué. Cette méthode, associée
à une méthode de préconditionnement (voir le Chapitre 3 pour une introduction à ce vaste
sujet) s’est imposée depuis maintenant 40 ans comme la méthode itérative la plus utilisée pour
résoudre des systèmes linéaires de tr1es grande taille, comme ceux provenant de la discrétisation
d’équations aux dérivées partielles.
Ce chapitre suit essentiellement [14].

1.1 Notations – généralités

Dans tout ce chapitre A ∈ Rn×n désigne une matrice définie positive. Nous allons étudier
une méthode itérative pour approcher la solution d’un résolution du système linéaire Ax = b (on
notera x∗ la solution exacte (x∗ = A−1 b). La méthode, que nous étudierons aux paragraphes 1.2
et suivants, va générer une suite xk (convergeant vers x∗ , la démonstration en sera donnée au
paragraphe 1.4). Étant donné un itéré xk , nous pouvons lui associer :
le résidu rk = b − Axk . Il s’agit d’une quantité calculable, elle jouera un rôle important
dans l’algorithme ;
l’erreur ek = x∗ − xk . Il s’agit d’une quantité qui n’est pas calculable sans connaissance
préalable de la solution exacte. Son intérêt est essentiellement théorique.
Notons la relation suivante entre ces 2 quantités, que l’on obtient en multipliant la définition de
l’erreur par A, et notant que Ax∗ = b :

(1.1) Aek = rk .

Pour préciser le lien entre le résidu et l’erreur, il est commode d’introduire le conditionnement
de la matrice A (on devrait plus précisément parler du conditionnement du système linéaire).

Définition 1.1. Le conditionnement de la matrice A (pour la norme matricielle k.k) est

(1.2) κ(A) = kAk A−1 .

La définition et les principales propriétés des normes matricielles sont rappelées dans l’an-
nexe A.
Notons les propriétés suivantes du conditionnement.

6
Proposition 1.1. Soit A ∈ Rn×n une matrice inversible.
— Pour λ ∈ R, κ(λA) = κ(A) ;
— κ(A) ≥ 1, et κ2 (A) = 1 si et seulement si A est une matrice orthogonale ;
— Si A est une matrice symétrique et définie positive, en notant 0 < λ1 ≥ · · · ≥ λn ses
valeurs propres, on a
λn
κ2 (A) = .
λ1
Démonstration. La démonstration se trouve dans de nombreux livres d’algèbre linéaire, voir les
références données dans l’annexe A. 

Nous pouvons maintenant donner le lien entre résidu et erreur

Proposition 1.2. Soit xk ∈ Rn et soit rk = b − Axk et ek = x∗ − xk respectivement le résidu

et l’erreur associés. On a
kek k krk k
(1.3) ≤ κ(A) .
ke0 k kr0 k

Démonstration. En inversant la relation entre résidu et erreur, il vient

kek k ≤ A−1 krk k .

De même (mais sans inverser), les quantités initiales sont reliées par

kr0 k ≤ kAk ke0 k .

L’inégalité cherchée découle des deux précédentes. 

Dans le cas important où la matrice A est symétrique et définie positive, la résolution du

système linéaire Ax = b est équivalente à la minimisation d’une fonctionnelle

Proposition 1.3. Un vecteur x∗ est solution du système linéaire Ax = b si et seulement si x∗

minimise la fonctionnelle J : Rn → R définie par
1
(1.4) J(x) = (Ax, x) − (b, x)
2
Démonstration. Une preuve rapide consiste à remarquer que la fonctionnelle J est différentiable,
et que son gradient est donné par
∇J(x) = Ax − b,
et que la fonctionnelle J est strictement convexe (car A est définie positive)
On peut également donner une démonstration élémentaire du résultat. On calcule, pour tout
y ∈ Rn :
1
J(x∗ + y) = J(x) + (Ay, y) + (Ax∗ − b, y)
2
où l’on a utilisé la symétrie de A pour obtenir le dernier terme.
Si x∗ est solution de Ax = b , le dernier terme de cette somme est nul, et comme A est
définie positive, on voit que x∗ réalise le minimum (global, strict) de J ;

7
 1
réciproquement, si x∗ minimise J , la quantité (Ay, y) + (Ax∗ − b, y) doit être positive
2
pour tout y 6= 0. Mais si x∗ n’est pas solution du système linéaire, le choix de y = εr∗ ,
avec r∗ = b − Ax∗ , qui n’est donc pas nul, donne
1
J(x∗ + y) − J(x∗ ) = ε2 (Ar∗ , r∗ ) − ε kr∗ k2
2
qui est négatif pour ε > 0 assez petit.

Lorsque A est définie positive, elle définit une norme, appelée norme de l’énergie selon la
définition
(1.5) kxk2A = (Ax, x), ∀x ∈ Rn
(exercice : vérifier que cela définit bien une norme). La fonctionnelle J est liée à cette norme par
la relation (noter que le dernier terme est constant)
(1.6) kx − x∗ k2A = 2J(x) + (Ax∗ , x∗ ).
qui exprime que minimiser J revient à minimiser l’erreur dans la norme de l’énergie.
Pour vérifier cette égalité, on développe le membre de gauche, en utilisant la symétrie de A
et la définition de x∗ :
kx∗ − xk2A = (Ax, x) − 2(Ax∗ , x) + (Ax∗ , x∗ ) = (Ax − 2b, x) + (Ax∗ , x∗ ).
Remarque 1.1. Lorsque la matrice A provient de la discrétisation d’un problème coercif par la
méthode des éléments finis, la norme de l’énergie est naturelle, et provient du produit scalaire
associé à la forme bilinéaire a. Soit uh et vh deux fonctions de l’espace discret Vh associé à
une triangulation d’un ouvert Ω ⊂ Rd (supposé polygonal ou polyèdral). On note (φi )i∈Nh les
fonctions de la base nodale de Vh (avec Nh = dim Vh ), et on note respectivement x et y les
vecteurs de coordonnées de uh et vh :
Nh
X Nh
X
uh = xi φi , vh = yi φ i .
i=1 i=1

On note enfin A la matrice de masse aij = a(φj , φi ) et on a la relation fondamentale :

a(uh , vh ) = (y, Ax) = (Ax, y),
qui implique entre autres que
kxk2A = a(uh , uh ).
On retrouve le lien bien connu entre le théorème de Lax–Milgram et la minimisation de la
fonctionnelle
1
j(uh ) = a(uh , uh ) − L(uh )
2
où L est la forme linéaire figurant au second membre de la formulation variationnelle.
Lorsque la forme bilinéaire a correspond à un opérateur de type « − div K grad u » dans un
domaine Ω borné (avec une fonction K ∈ L∞ (Ω) vérifiant K(x) ≥ k∗ > 0 pour presque tout x ∈
Ω), la forme bilinéaire s’écrit Z
a(u, u) = K |∇u|2
Ω
et correspond effectivement à une énergie (thermique).

8
1.2 Définition de la méthode
La méthode du gradient conjugué construit une suite d’itérés en cherchant à minimiser l’er-
reur dans la nome de l’énergie k[kA ]e sur une suite de sous-espaces de dimensions croissantes. Par
construction, la méthode trouvera forcément la solution exacte du problème (après n itération),
mais sa réelle utilité vient de la qualité des approximations obtenues après un petit nombre
d’itérations (quelques centaines pour des problèmes avec plusieurs millions d’inconnues). Com-
mençons par définir les sous-espaces en questions.
Définition 1.2. Étant donné un vecteur r0 ∈ Rn , l’espace de Krylov de degré k relatif à A et
r0 ) est
n o
(1.7) Kk (A, r0 ) = Kk (r0 ) = vect r0 , Ar0 , . . . , Ak−1 r0

On notera les propriétés suivantes, dont la démonstration est laissée au lecteur :

— Kk (r0 ) ⊂ Kk+1 (r0 ), et donc la dimension de Kk (r0 ) est une fonction (non strictement
croissante) de k, avec dim Kk (r0 ) ≤ k.
— on a x ∈ Kk (r0 ) ⇐⇒ x = p(A)r0 , pour un polynôme p de degré inférieur ou égal à k − 1.
Notons également que si la suite Kk (r0 ) stagne, c’est-à-dire si Kk (r0 ) ⊂ Kk+1 (r0 ) alors
la solution du système linéaire Ax = b est contenue dans l’espace affine x0 + Kk (r0 ) (avec
r0 = b − Ax0 ).
En effet, dans ce cas, Ak r0 Kk (r0 ), et donc Ak r0 = αr0 + Aq(A)r0 pour un polynôme q de
degré strictement inférieur à k − 1, et un scalaire α.
α ne peut pas être nul, car sinon on aurait Ak r0 = Aq(A)r0 , soit (Ak−1 − q(A))r0 = 0,
puisque A est inversible, et comme q est de degré strictement inférieur à k − 1, A − q(A) n’est
pas nul, et donc r0 = 0, et x0 est solution du problème initial.
On peut donc écrire
1  
b = Ax0 + r0 = Ax0 + A Ak−1 − q(A) r0 ,
α
et nous venons de voir que la matrice Ak−1 −q(A) n’est pas nulle. La solution du système linéaire
1
Ak−1 − q(A) r0 ∈ x0 + Kk (r0 ).


est donc x0 +
α
Par conséquent, lors de l’étude d’un algorithme dont les itérées sont dans l’espace de Krylov
Kk (r0 ), on pourra sans perte de généralité faire l’hypothèse que l’inclusion Kk (r0 ) ⊂ Kk+1 (r0 )
est stricte (sinon, l’algorithme a déjà trouvé la solution du système linéaire).
Nous pouvons maintenant donner la définition de la méthode du gradient conjugué.
Définition 1.3. Étant donné un itéré initial x0 ∈ Rn , l’itéré xk de la méthode du gradient
conjugué réalise le minimum de la fonctionnelle J définie en (1.4) sur l’espace de Krylov Kk (r0 ).
Notons tout de suite que xk est bien défini de façon unique par cette condition. On peut le
voir en appliquant le théorème de projection sur un convexe (ici un sous-espace affine) pour le
produit scalaire associé à la norme de l’énergie (voir l’équation (1.6)). On peut préciser cette
caractérisation. Dans toute la suite, nous noterons rk = b − Axk .
Lemme 1.1. A l’itération k, xk est caractérisé par les conditions suivantes :
i) xk ∈ x0 + Kk (r0 ),

9
 ii) rk ⊥ Kk (r0 ).

Démonstration. En utilisant (1.6), nous voyons que xk est la solution du problème de minimi-
sation
min kxk − x∗ kA ,
x0 +Kk (r0 )

autrement dit, xk est la projection de la solution exacte x∗ sur le sous-espace affine (convexe,
fermé) x0 + Kk (r0 ). Cette projection est caractérisée par les conditions
i) xk ∈ x0 + Kk (r0 ) (identique à i) ci-dessus,
ii) xk − x∗ ⊥A Kk (r0 ) (orthogonalité par rapport au produit scalaire de l’énergie).
Il reste à expliciter la seconde condition :

xk − x∗ ⊥A Kk (r0 ) ⇐⇒ ∀y ∈ Kk (r0 ), (y, A(xk − x∗ )) = 0

⇐⇒ ∀y ∈ Kk (r0 ), (y, (Axk − b)) = 0 ⇐⇒ rk ⊥ Kk (r0 ).

Cette définition de l’algorithme du gradient conjugué reste abstraite. Nous allons maintenant
en donner une version « réalisable ».

1.3 L’algorithme
Le développement de la version plus algorithmique du gradient conjugué utilise diverses
propriétés des itérés que nous établirons dans une série de lemmes. Nous donnerons ensuite des
éléments concernant la mise en oeuvre de la méthode.

1.3.1 Développement de l’algorithme

Lemme 1.2. Pour l < k, on a rl ∈ Kk (r0 ), et (rl , rk ) = 0.
Par conséquent, Kk (r0 ) = vect {r0 , · · · , rk−1 }.

Démonstration. Par définition, rl = b − Axl = r0 + A(x0 − xl ) ∈ Kl (r0 ) + AKl (r0 ) ⊂ Kl+1 (r0 ) ⊂
Kk (r0 ), puisque l + 1 ≤ k.
D’après le lemme 1.1, rk est orthogonal à Kk (r0 ), et d’après le premier point, ce sous espace
contient tous les résidus précédents.
Le premier point montre l’inclusion Kk (r0 ) ⊃ vect {r0 , · · · , rk−1 }, et l’orthogonalité des
résidus montre l’égalité des dimensions. 

Ce résultat n’est évidemment valable que si rk 6= 0, mais si rk = 0, alors par définition xk

est solution du système.
Par ailleurs, l’algorithme ne peut pas stagner : en effet, si xk = xk+1 , alors rk = rk+1 , et le
lemme précédent montre que krk k22 = (rk , rk+1 ) = 0, et donc ce cas ne peut se produire que si
xk est solution du système.
Nous pouvons donc définir un vecteur pk comme la direction du changement entre deux itérés
successifs. Il est commode de poser

(1.8) xk+1 = xk + αk pk , et p0 = r0

10
où αk est un scalaire qui reste à déterminer. Notons qu’à ce stade, pk est défini à une constante
près, et αk pourrait être « absorbé » dans la définition de pk . Nous verrons que cette indétermi-
nation sera résolue naturellement au cours de la détermination de l’algorithme.

Lemme 1.3. La direction de descente pk est déterminée à un scalaire près par les conditions :
i) pk ∈ Kk+1 (r0 ) ;
ii) pk ⊥A Kk (r0 ) (c’est-à-dire (y, Apk ) = 0, ∀y ∈ Kk (r0 ).

Démonstration. Comme les espaces de Krylov sont emboîtés, xk et xk+1 sont dans le même
espace affine, donc leur différence est dans l’espace vectoriel Kk+1 (r0 ), et la remarque précédent
le lemme montre que l’on peut supposer αk 6= 0.
Pour montrer que pk ⊥A Kk (r0 ), partons de

(1.9) rk+1 = b − Axx − αk Apk = rk − αk Apk .

Comme rk+1 est orthogonal à Kk+1 (r0 ), et rk est orthogonal à Kk (r0 ), leur différence est bien
orthogonale à Kk (r0 ).
Comme pk est orthogonal à un sous-espace de codimension 1 dans Kk+1 (r0 ), sa direction est
bien déterminé de manière unique. 

On dit que pk est A-conjugué à Kk (r0 ), ce qui explique l’origine du nom de la méthode.
Ce lemme montre que l’on peut déterminer pk (à un facteur multiplicatif près) sans connaître
xk+1 , et (1.8) permet alors de calculer xk+1 . On pourrait croire que la détermination de pk est
coûteuse. Il n’en est rien, et c’est ce qui constitue l’élément le plus étonnant de l’algorithme du
gradient conjugué !
On montre, par un raisonnement similaire à celui du lemme 1.2, que

Kk (r0 ) = vect (p0 , . . . , pk−1 ).

Lemme 1.4. On a

(1.10) pk = rk + βk−1 pk−1 ,

avec
(rk−1 , Ark )
βk−1 = − .
(rk−1 , Apk−1 )
Démonstration. La première condition du lemme 1.3 et la remarque précédente montrent que,
a priori, pk s’écrit sous la forme
k−1
X
pk = rk + γkj pj
j=0

où les coefficients γkj , j = 0, . . . , k − 1 doivent être déterminés. Pour cela, nous utilisons la
deuxième condition du lemme 1.3. Multiplions l’équation précédente par A, et prenons le produit
scalaire avec rl , pour l = 0, . . . , k − 1
k−1
X
(rl , Apk ) = (rl , Ark ) + γkj (rl , Apj ), l = 0, . . . , k − 1
j=0

11
 Pour l = 0, . . . , k −1, on a (rl , Apk ) = 0 d’après le lemme 1.3, puisque rl ∈ Kk (r0 ). De même,
mais pour l = 0, . . . , k − 2, (rl , Ark ) = (rk , Arl ) = 0 car Arl ∈ Kk (ro ) et rk est orthogonal à ce
sous espace.
Pour les termes de la somme, pour l = 0, . . . , k − 2, utilisons la relation (1.9). Pour l = 0, le
seul terme non-nul dans la somme est le premier (avec j = 0) qui vaut γk0 /α0 ((r1 − r0 ), r0 ) =
−γk0 /α0 kr0 k22 . On en déduit que γk0 = 0. On procède de même pour les autres termes, pour
montrer successivement que γkj = 0 pour j = 0, . . . , k − 2.
Il ne reste plus que le terme correspondant à j = k − 1, qui conduit à la condition
0 = (rk−1 , Ark ) + γk,k−1 (rk−1 , Apk−1 )
et donne la valeur annoncée pour βk−1 = γk,k−1 , à condition de montrer que (rk−1 , Apk−1 ) 6= 0.
Pour cela, notons que, comme rk = rk−1 − αk−1 Apk−1 , on a
(rk , rk−1 ) = krk−1 k22 − αk−1 (rk−1 , Apk−1 ).
Le terme de gauche est nul par orthogonalité des résidus successifs. Comme le premier terme
du second membre n’est pas nul (sauf si la méthode a convergé), il en est de même du second
terme. 

Il reste donc à calculer le coefficient αk . Pour cela, nous pouvons utiliser le calcul itératif du
résidu (1.9), et l’orthogonalité de rk+1 et pk (conséquence de la propriété de minimisation vue
au lemme 1.1 et de ce que pk ∈ Kk+1 (ro )) :
(pk , rk+1 ) = (pk , rk ) − αk (pk , Apk )
qui permet de déterminer αk (comme A est définie positive, le dénominateur est strictement
positif).
L’algorithme du gradient conjugué consiste essentiellement à séquencer les équations (1.8), (1.9)
et (1.10). Avant de détailler l’algorithme, notons que les formules pour αk et βk utilisées en pra-
tique ne sont pas celles que nous avons obtenues ci-dessus, mais des formes équivalentes qui ont
de meilleures propriétés de stabilité numérique.
Lemme 1.5. Les quantités αk et βk peuvent se calculer par les formules suivantes
krk k22
(1.11) αk =
(pk , Apk )

krk+1 k22
(1.12) βk =
krk k22
Démonstration. Pour αk il suffit de noter que (pk , rk ) = krk k22 . Pour le voir, il suffit d’utili-
ser (1.10), et l’orthogonalité de rk avec pk−1 ∈ Kk (r0 ).
Pour βk , la relation de récurrence (1.10) donne
(Apk , pk+1 ) = (Apk , rk+1 ) + βk (Apk , pk )
et le membre de gauche est nul puisque les directions successives sont A-conjuguées. On en
déduit
(Apk , rk+1 )
βk = − .
(Apk , pk )
Calculons séparément le dénominateur et le numérateur.

12
 — D’après l’expression de αk obtenue dans la première partie du lemme,

krk k22
(Apk , pk ) = ;
αk
— En utilisant (1.9), il vient (par orthogonalité des résidus successifs)

krk+1 k22 = −αk (Apk , rk+1 ).

On en déduit l’expression (1.12) pour βk .



Nous pouvons maintenant donner une version de l’algorithme est proche de la version im-
plémentable, mais qui reste idéale sur plusieurs points (critère d’arrêt, gestion de la mémoire).

Algorithme 1.1 Algorithme du gradient conjugué

Données: b ∈ Rn , x0 ∈ Rn , A ∈ Rn×n
Résultats: k, xk
r0 = b − Ax0 , p0 = r0 , k = 0
Tant que k ≤ itmax et pas de convergence faire
krk k22
αk =
(pk , Apk )
xk+1 = xk + αk pk
rk+1 = rk − αk Apk
krk+1 k22
βk =
krk k22
pk+1 = rk+1 + βk pk
fin Tant que

1.3.2 Mise en oeuvre de l’algorithme

L’algorithme précédent 1.1 est très simple à mettre en oeuvre : il n’est pas nécessaire de
conserver les itérés successifs, et l’algorithme peut être implémenté avec seulement 4 vecteur. Le
coût d’une itération est de 4 opérations vectorielles (ces opérations sont toutes du type y = y+ax,
où x e t y sont des vecteurs et a est un scalaire), le calcul de 2 produits scalaires et une norme,
et le produit de la matrice A par un vecteur. Dans beaucoup d’applications (notamment dans la
résolution d’équations aux dérivées partielles), la matrice A possède une structure particulière
permettant de calculer ce produit pour un coût très inférieur à n, ce qui serait le cas si A est
stocké comme une matrice pleine générale. Comme nous le verrons au Chapitre 4, quand A est
une matrice creuse, le produit matrice vecteur peut se calculer pour un coût proportionnel aux
nombres d’éléments non-nuls de la matrice. Quand la matrice provient de la discrétisation d’une
équation aux dérivées partielles par différences finies, éléments finis, ou volumes finis, le nombre
d’éléments non-nuls sur une ligne est approximativement constant, ou au moins borné par une
constante. Dans ce cas, le produit matrice – vecteur à un coût linéaire par rapport aux nombre
d’inconnues. Cette propriété se transmet au coût total d’une itération.
Le coût total de l’algorithme dépend naturellement du nombre d’itérations. Cette question
sera traitée au paragraphe 1.4.

13
Remarque 1.2 (Mis en oeuvre sans matrice). En poursuivant la discussion précédente, on voit
qu’il n’est pas nécessaire de la matrice A sous la forme d’un adressage (i, j) → Aij . Il suffit
de disposer d’une fonction (un sous-programme) qui prend en entrée un vecteur v et renvoie
le résultat du produit Av. Cela veut dire que l’on revient à la vision de l’application linéaire
sous-jacente à la matrice. En Matlab, et aussi avec la bibliothèque Scipy de Pyhon, il est possible
d’utiliser le gradient conjugué soit avec une matrice stockée en mémoire, soit avec un opérateur
linéaire (voir des exemples au paragraphe ??.
Cette remarque donne lieu à ce que l’on appelle en anglais un mise en oeuvre matrix free.
Remarque 1.3 (Arrèt de l’algorithme). Nous n’avons pas réellement spécifié la condition d’arrêt
de l’algorithme 3.2, en indiquant simplement « pas de convergence ». Bien entendu, ceci doit être
précisé. Nous pouvons d’ores et déjà noter que l’algorithme 3.2 doit s’arrêter après au plus n
itérations : en effet, le lemme 1.2 montre que les résidus successifs (r0 , . . . , rk−1 forment )une base
orthonormale de Kk (r0 ), et le nombre de vecteurs orthogonaux ne peut excéder la dimension de
l’espace. Cependant ce résultat n’est jamais utile en pratique. On utilise couramment la méthode
du gradient conjugué pour résoudre des systèmes comportant plusieurs millions d’inconnues,
pour lesquels on arrête l’algorithme après quelques centaines d’itérations, voire moins.
Ce fait remarquable justifie a posteriori le choix des espaces de Krylov pour chercher la
solution : ils fournissent très rapidement des approximations précises de la solution. Ce point
sera quantifié au paragraphe 1.4, et l’on parle (improprement) de la convergence de l’algorithme.
Nous pouvons indiquer dès maintenant le choix du critère d’arrêt le plus utilisé en pratique :
— Un critère d’arrêt ne doit utiliser que des quantités calculables. L’erreur n’entre évidem-
ment pas dans ce cadre, et l’on doit utiliser le résidu, les deux quantités étant liées par
l’équation 1.1. On arrête l’algorithme quand le résidu est suffisamment petit ;
— La notion de « petit » ne peut être absolue, et doit relative à la taille des données.
On recommande donc le critère d’arrêt suivant

(1.13) krk ≤ η kbk ,

où η (la tolérance) est une quantité spécifiée par l’utilisateur de la méthode, indiquant quand le
résidu est suffisamment petit (on prend en général une valeur entre 10−12 et 10−6 . Nous n’avons
pas précisé la norme utilisée, qui sera le plus souvent la norme euclidienne.
Une discussion détaillée des critères d’arrêt pour les méthodes itératives (ainsi que de nom-
breuses informations pratiques) se trouve dans le livre [4]. Ce livre explique entre autres pourquoi
il ne faut pas remplacer, dans (1.13), la norme du second membre kbk par celle du résidu initial
kr0 k.
Nous pouvons maintenant donner une forme réalisable de l’algorithme du gradient conjugué,
prenant en compte à la fois la gestion de la mémoire et le critère d’arrêt. Le seul point délicat
est qu’il faut conserver les normes de deux résidus successifs.

1.4 Convergence
Commençons par une observation simple. Par sa définition, la méthode du gradient conjugué
trouve la solution exacte du système Ax = b en au plus n itérations, si n est la dimension du
système. En effet, les sous-espaces de Krylov sont de dimension au plus n. Soit l’algorithme a
trouvé une solution après k < n itérations, soit à l’itération n la minimisation est posée sur tout
l’espace, et l’itéré xn est la solution.

14
Algorithme 1.2 Algorithme du gradient conjugué
Données: b ∈ Rn , x ∈ Rn , A ∈ Rn×n , η > 0
Résultats: k, xk
r = b − Ax, p = r, k = 0, rold = krk22
Tant que k ≤ itmax et rold ≥ η kbk2 faire
roll
α=
(p, Ap)
x = x + αp
r = r − αAp
rnew
rnew = krk22 , β =
rold
p = r + βp
rold = rnew , k = k + 1
fin Tant que

Au moment où elle a été proposée, la méthode du gradient a d’abord été considérée comme
une méthode directe, à cause de cette propriété. Mais ce point de vue n’est plus celui adopté à
l’heure actuelle :
— tout d’abord la propriété concernant la solution exacte n’est plus vraie lorsque la méthode
est implémentée en arithmétique en précision limitée. Dans ce cas, l’algorithme peut ne
pas trouver la solution avant bien plus de n itérations ;
— Mais même si l’on considère la situation idéalisée, ce résultat n’a pas d’intérêt pratique.
La méthode du gradient conjugué est couramment utilisée pour résoudre des systèmes
linéaires à plusieurs centaines de milliers, voire des millions, d’inconnues. Or nous verrons
dans la suie de ce paragraphe que la méthode peut produire des approximations précises
de la solution en quelques centaines d’itérations, parfois encore moins.
Ainsi, le point de vue « moderne » sur la méthode du gradient conjugué est comme une
méthode itérative, produisant une très bonne approximation de la solution exacte du système
en un nombre limité d’itérations. C’est en ce sens que nous parlerons de « convergence », c’est
à dire en majorant la différence entre l’itéré xk et la solution exacte. Nous verrons que cela est
possible en choisissant une norme liée au problème.

La première étape pour établir ces résultats est de réinterpréter la méthode du gradient
conjugué en terme de polynômes. Nous noterons Pk l’espace vectoriel des polynômes de degré
inférieur ou égal à k.

Proposition 1.4. Soit x∗ la solution du système Ax = b, et soi xk l’itéré produit par l’algorithme
du gradient conjugué après k itérations.
i) il existe un polynôme p ∈ Pk , vérifiant p(0) = 1, tel que

(1.14) x∗ − xk = p(A)(x∗ − x0 );

ii) L’itéré xk est caractérisé par la condition

(1.15) kx∗ − xk kA = min kp(A)(x∗ − x0 )kA

p∈Pk
p(0)=1

15
Preuve. i) Puisque xk ∈ x0 + Kk (r0 ), on a, pour certains coefficients γj , j = 0, . . . , k − 1
k−1
X
xk − x∗ = x0 − x∗ + γj Aj r0 .
j=0

En utilisant l’égalité r0 = Ax∗ − Ax0 on obtient (1.14) si l’on définit le polynôme p(z) = 1 −
Pk−1 j+1
j=0 γj z
ii) L’égalité (1.15) est une conséquence immédiate du point précédent et de la définition 1.3.


Pour exploiter ce résultat, nous utiliserons les valeurs propres et vecteurs propres de la
matrice A. Notons qu’il s’agit là d’un outil uniquement théorique, la connaissance des valeurs
propres de A n’est aucunement nécessaire à la mise en oeuvre de l’algorithme, comme on l’a
vu au paragraphe 1.3.2. Comme A est symétrique, elle est diagonalisable avec une matrice de
passage orthogonale, et ses valeurs propres sont réelles (et strictement positives, puisqu’elle est
définie positive).
On a donc l’égalité

(1.16) A = U ΛU T ,

où la matrice (carrée) U est orthogonale (U T U = I), et la matrice Λ est diagonale, Λ =

diag(λ1 , . . . , λn ), avec 0 < λ1 ≤ . . . ≤ λn . On notera σ(A) = {λ1 , . . . , λn } le spectre de A
(l’ensemble de ses valeurs propres).
Une conséquence immédiate est que, pour tout polynôme p, on également

p(A) = U p(Λ)U T , avec p(Λ) = diag(p(λ1 ), . . . , p(λn )).

On notera également A1/2 la matrice

1/2
A1/2 = U Λ1/2 U T avec p(Λ) = diag(λ1 , . . . , λ1/2
n ).

C’est l’unique matrice symétrique et définie positive vérifiant A1/2 · A1/2 = A.

Le résultat suivant exploite ces notions.

Lemme 1.6. Pour tout vecteur x ∈ Rm , et pour tout polynôme p on a :

i) x A
= A1/2 x 2 ,
ii) kp(A)xkA ≤ max |p(A)| kxkA .
λ∈σ(A)

Preuve. Il s’agit dans les deux cas d’un simple calcul.

2
i) x A
= xT Ax = X T A1/2 A1/2 x = (A1/2 x)T (A1/2 x) = A1/2 x 2
;
2 2
ii) D’après le point précédent, p(A)x A
= A1/2 p(A)x 2
, or A1/2 et p(A) commutent, d’où
2
2 2
kxkA ≤ kp(A)k2 A1/2 x ,
2

où kp(A)k2 est la norme matricielle, et la conclusion suit en remarquant que, pour x 2 = 1,

p(A)x 2 = U p(Λ)U T x 2 = p(Λ)y 2 , avec y 2 = U T x 2 = 1 (dans les 2 cas parce que U
est orthogonale), et que la norme euclidienne d’une matrice diagonale est égale à son plus grand
élément.

16
 

Ces résultats préliminaires nous permettent d’énoncer le résultat fondamental suivant.

Théorème 1.1. les itérés de la méthode du gradient conjugué vérifient l’inégalité suivante :

(1.17) kxk − x∗ kA ≤ min max |p(λ)| kx0 − x∗ kA .

p∈Pk λ∈σ(A)
p(0)=1

Preuve. La preuve consiste simplement à expliciter la norme dans la définition 1.3 en utilisant le
lemme 1.6. 

La version simplifiée donnée au corollaire suivant, dont la preuve est immédiate, est souvent
suffisante pour les applications.
Corollaire 1.1. Pour tout polynôme p ∈ Pk , avec p(0) = 1, on a

kxk − x∗ kA ≤ max |p(λ)| kx0 − x∗ kA .

λ∈σ(A)

Pour appliquer le Corollaire 1.1, il suffit de trouver un polynôme adapté à la situation. Une
application particulièrement simple est de retrouver la terminaison en un nombre fini d’itérations
de la méthode du gradient conjugué. On choisit pour cela le polynôme caractéristique de A, qui
est explicite si l’on connaît les valeurs propres.
n  
Y λ
p(λ) = 1− .
λi
i=1

Ce polynôme s’annule en chaque valeur propre, le second membre de la borne du Corollaire 1.1
est donc nul, d’où xn = x∗ .
Dans le cas général, c’est-à-dire sans autre hypothèse sur la répartition des valeurs propres
de la matrice A, on obtient le résultat de convergence suivant :
Proposition 1.5. On note respectivement λ1 et λn la plus petite et la plus grande valeur propre
λn
de A, et on note κ(A) = le conditionnement de A. Les itérés de la méthode du gradient
λ1
conjugué vérifient :
p !k
κ(A) − 1
(1.18) kxk − x∗ kA ≤ 2 p kxk − x∗ kA .
κ(A) + 1

(partielle). Le maximum sur l’ensemble discret des valeurs propres est difficile à déterminer.
Sans autre information que la valeur des valeurs propres extrêmes, il est naturel de chercher le
minimum sur l’intervalle [λ1 , λn ], ce qui a naturellement pour effet d’augmenter le minimum du
second membre de l’inégalité (1.17) en

kxk − x∗ kA ≤ min max |p(λ)| kx0 − x∗ kA .

p∈Pk λ∈[λ1 ,λn ]
p(0)=1

La détermination du min-max est un problème d’approximation polynomiale, dont la solution est

donnée par le polynôme de Tchebychev adapté à l’intervalle en question. Cette partie, technique,
de la démonstration sera omise. Le lecteur pourra trouver les détails dans les références [1, para.
9.5.5.] ou [15]. 

17
 On retrouve ici le conditionnement, non pas comme indicateur de la sensibilité de la matrice
aux perturbations, mais comme indicateur de la vitesse de convergence. Plus le conditionnement
sera grand, plus le facteur gouvernant la convergence sera proche de 1, et plus cette convergence
sera lente. On notera que c’est la racine carrée du conditionnement qui intervient, ce qui constitue
un progrès par rapport à la méthode de relaxation (SSOR), où la vitesse de convergence dépend
du conditionnement.

Exemple 1.1 (Le laplacien sur un carré).

à titre d’illustration, prenons l’exemple de la discrétisation de l’équation de Poisson avec des
conditions de Dirichlet homogène sur le carré unit par une méthode de différences finies. On
résout l’équation aux dérivées partielles

−∆u = f dans Ω =]0, 1[×]0, 1[, u = 0 sur ∂Ω.

La discrétisation par une méthode de différences finies sur une grille régulière de pas h (on
supposera que h = 1/(N + 1) où N est un entier) conduit à un système linéaire Ax = b où la
2 2
matrice A ∈ RN ×N est
 
4 −1 −1
 
−1 4 −1 −1 
 
−1 4 −1 −1
 
 
 
−1 4 −1 −1
 
 
 

 −1 4 −1 

. . . . 
 .. .. .. .. 
 
 .. .. .. .. .. 

 . . . . . 

1 
A = 2
 −1 −1 4 1 −1 
h  
.. .. .. .. .. 

 . . . . . 

 .. .. .. .. 

 . . . .

 

 −1 4 −1 

 

 −1 −1 4 −1 

−1 −1 4 −1
 
 
 
−1 −1 4 −1
 

−1 −1 4

2
et où le vecteur b ∈ RN contient des approximations des valeurs de la fonction f aux sommets
de la grille.
On sait que le conditionnement de la matrice A vaut
 
2 nπ
cos
2(n + 1) 4
κ(A) =   ≈ 2 2,
π π h
sin2
2(n + 1)

18
ce qui conduit à un taux de convergence de
p
κ(A) − 1
ρ= p ≈= 1 − πh + O(h2 ).
κ(A) + 1

On en déduit que le nombre d’itérations nécessaires pour réduire la norme du résidu initial d’un
facteur donné est (approximativement) proportionnel au nombre de points de discrétisation sur
un côté du carré.
Ce résultat sera confirmé numériquement au paragraphe ??.

Cet exemple est représentatif de ce que l’on peut observer lorsqu’on résout une équation
aux dérivées partielles (de type elliptique) par la méthode du gradient conjugué : pour atteindre
un même niveau de réduction de la norme du résidu, le nombre d’itérations augmente lorsque
l’on raffine le maillage. Comme le coût d’une itération augmente également (à peu près linéai-
rement par rapport au nombre d’inconnues), le coût total de la méthode augmente plus vite
que le nombre d’inconnues. On dit que la méthode du gradient conjugué n’est pas optimale. Il
existe des méthodes optimales, la plus connue étant la méthode multigrille (voir la section ??).
Néanmoins, la méthode du gradient conjugué continue d’être utilisée, en partie pour sa simpli-
cité de mise en oeuvre, mais également pour sa plus grande robustesse en comparaison avec la
méthode multigrille. En réalité, les deux méthodes sont utilisées ensemble : on peut soit voir le
gradient conjugué comme un accélérateur de la méthode multigrille, soit le multigrille comme
un préconditionenur (voir la section ??).
Remarque 1.4 (Influence de la répartition des valeurs propres).
Pour diminuer le nombre d’itérations de la méthode, on peut chercher à formuler, et à
résoudre, un système équivalent au système original, mais dont le conditionnement est plus
petit (ou plus généralement, dont les valeurs propres ont une répartition plus favorable pour la
convergence). On parle de préconditionnement du système linéaire, un sujet que nous abordons
au chapitre 3.

19
Chapitre 2

Systèmes linéaires non symétriques : la

méthode GMRES

Dans ce chapitre nous considérons le cas d’un système linéaire

(2.1) Ax = b, A ∈ Rn×n , b ∈ Rn ,
où la matrice A n’est pas nécessairement inversible. Notez que A sera toujours supposée in-
versible. Dans ce cas, la Proposition 1.3 n’est plus valable, la résolution du système n’est plus
équivalente à la minimisation d’une fonctionnelle d’énergie.
Nous suivons également l’approche de Kelley [14, Chap. 3]
Il existe de nombreuses situations conduisant à la résolution d’un système non–symétrique.
Citons :
— la discrétisation d’une équation d’advection–diffusion
−∆u + div(~qu) = f dans Ω, u = 0 sur ∂Ω,
où Ω est un ouvert de Rd (d = 1, 2 ou 3), et ~q est un champ de vecteurs donné ;
— la résolution des équations de Navier–Stokes, compressibles ou incompressibles ;
— certains problèmes d’élasticité non–linéaire ;
— de manière générale la résolution de problèmes non–linéaires par la méthode de Newton
passe par la résolution, à chaque itération, d’un système linéaire dont la matrice est
presque toujours non–symétrique.

2.1 Définition de la méthode GMRES, et premières propriétés

L’idée de la méthode GMRES (Résidu Minimum Généralisé, en anglais « Generalized Mi-
nimimum RESidual », inventée en 1986 par Y. Saad et M. H. Schultz [19]) est de minimiser
la norme euclidienne du résidu sur les espaces de Krylov correspondant au résidu initial. De
manière précise :
Définition 2.1. Étant donné un itéré initial x0 ∈ Rn , l’itéré xk de la méthode GMRES réalise
le minimum de
(2.2) kb − Axk22
sur l’espace affine de Krylov x0 + Kk (r0 ).

20
 L’espace de Krylov Kk (r0 ) a été défini au chapitre 1. Cette définition ne dépend pas de la
symétrie de A, et nous pouvons donc toujours écrire
k−1
X
xk = x0 + γj Aj r0 ,
j=1

pour des coefficients γ1 , . . . , γk−1 . On en déduit

k−1
X k−1
X
j
b − Axk = b − Ax0 − A γj A r0 = r0 − γj−1 Aj r0 .
j=1 j=1

Il existe donc un polynôme p ∈ Pk , avec p(0) = 1 tel que, si xk ∈ x0 + Kk (r0 ), rk = b − Axk

s’écrit sous la forme
(2.3) rk = p(A)r0 .
Une conséquence de cette relation est que, comme le gradient conjugué, GMRES trouve
toujours la solution du système après un nombre fini d’itérations.
Proposition 2.1. Après au plus n itérations, l’itéré calculé par la méthode GMRES est égal à
la solution exacte du système (2.1).
Démonstration. On choisit pour p le polynôme caractéristique de A, qui est de degré n. Plus
précisément, pour tenir compte de la normalisation p(0) = 1, on prend
p(λ) = det(A − λI)/ det(A).
Notons que det(A) 6= 0 puisque nous avons supposé A inversible.
Le théorème de Cayley–Hamilton ([17, Thm 7.1.1], [11, par. 4.2, Thm 3] ou [21, Thm 3.3])
nous indique que p(A) = 0, et d’après les remarques précédant la proposition, cela implique
rn = 0. 

2.2 L’algorithme
2.2.1 L’algorithme d’Arnoldi
Supposons connue une base de l’espace de Krylov Kk (r0 ), et notons Vk ∈ Rk×n = [v1k , . . . , vkk ]
la matrice dont les colonnes sont les vecteurs de cette base. Un choix possible (mais nous verrons
que ce n’est pas le meilleur) est bien entendu de prendre vk = Ak−1 r0 . La solution approchée à
l’étape k de la méthode GMRES s’écrit donc sous la forme
k
X
xk = x0 + yj vjk = x0 + Vk yk , yk = (y1k , . . . , ykk )T ∈ Rk ,
j=1

pour un vecteur yk qui reste à déterminer (l’itéré xk ne dépend pas du choix de la base Vk , mais
sa représentation yk en dépend). Le problème de minimisation qui définit xk devient donc
min kr0 = AVk yk22 .
y∈Rk

Cette expression sera valable pour tout choix de la matrice Vk . Mais le choix naturel de Vk
indiqué précédemment pose des difficultés du point de vue de la mise en oeuvre numérique :

21
 — L’espace de Krylov Kk (r0 ) est un bon espace pour approcher la solution, mais la base
de puissances de Vk = [r0 , Ar0 , . . . , Ak−1 r0 conduit à une matrice vk mal conditionnée,
et ceci d’autant plus que k augmente. Ce constat est lié à la méthode de la puissance
pour calculer les valeurs propres (voir par exemple [16, Chap. 10]), qui montre que les
vecteurs Ak r0 deviennent parallèles au vecteur propre associé à la valeur propre de plus
grand module (si celle-ci est simple, et moyennant une hypothèse raisonnable sur r0 ) ;
— de manière générale, si la matrice Vk n’est pas orthogonale, le problème de moindres
carrés précédent sera difficile et coûteux à résoudre.
Pour ces raisons, une première étape consiste à chercher une base orthonormée de l’espace
de Krylov Kk (r0 ). Cela peut être réalisé par la méthode d’Arnoldi qui consiste à appliquer la
méthode d’orthogonalisation de Gram–Schmidt à la suite des sous–espaces Kj (r0 ), pour j variant
de 1 à k. Dans toute la suite, il sera commode de noter β = kr0 k2 , que l’on peut supposer non–nul.
À chaque étape, on suppose que l’on dispose d’une base orthonormée (v1 , . . . , vk ) de l’espace
de Krylov Kk (r0 ), et on doit donc définir le vecteur vk+1 de sorte que (v1 , . . . , vk+1 ) soit une base
orthonormée de Kk+1 (r0 ). Pour cela, il suffit d’orthogonaliser le nouveau vecteur Avk contre les
éléments précédents de la base. Pour cela, on calcule
k
X
wk+1 = Avk − (Avk , vj )vj ,
j=1
w
k+1
et on définit vk+1 = kwk+1 k2 .
Dans sa version la plus simple, qui sera améliorée plus loin, cela conduit à l’algorithme 2.1.
Par construction, les vecteurs produits par cet algorithme sont orthogonaux, tant que hk+1,k 6= 0.

Algorithme 2.1 Algorithme d’Arnoldi, version « naïve »

Données: k ∈ N, base orthonormée (v1 , . . . , vk ) de Kk (r0 ).
Résultats: Base orthonormée (v1 , . . . , vk+1 ) de Kk+1 (r0 )
v1 = r0 /β
pour j = 1, . . . , k faire
hjk = (Avk , vj )
fin pour
wk+1 = Avk
pour j = 1, . . . , k faire
wk+1 = wk+1 − hjk vj
fin pour
hk+1,k = kwk+1 k2
si hk+1,k = 0 alors stop
vk+1 = wk+1 /hk+1,k

Mais si cette dernière éventualité se produit, alors l’algorithme a trouvé la solution du système
linéaire.
Lemme 2.1. Dans l’algorithme 2.1, si hk+1,k = 0 alors x? ∈ Kk (r0 ).
Pk
Démonstration. Par hypothèse, Avk = j=1 hjk vj , et donc Avk ∈ Kk (r0 ), soit AKk (r0 ) ⊂
Kk (r0 ). Les colonnes de Vk forment une base orthonormée de Kk (r0 ), donc

AVk = Vk H̃k

22
où la matrice H̃k ∈ Rk×k est inversible (puisque A l’est). On en déduit

rk = b − Axk = r0 − A (xk − x0 ) = r0 − Vk H̃k y,

| {z }
∈Kk (r0 )

pour y ∈ Rk .
Par défiintion, r0 = βv1 = kr0 k2 Vk e1 (où e1 est le premier vecteur de la base canonique de
Rk ), et par conséquent

krk k2 = βVk e1 − Vk H̃y = kβe1 − Hk yk2 ,

2

puisque Vk est orthogonale. Comme H̃k est inversible, on peut choisir y de manière à ce que
rk = 0, et xk ∈ x0 + Kk (r0 ) est alors la solution de (2.1). 

Nous pouvons donc faire dans la suite l’hypothèse que l’algorithme d’Arnoldi ne s’arrête pas,
et que nous disposons, à chaque itération k d’une
 base orthonormée
 (v1 , . . . , vk ) de l’espace de
Krylov Kk (r0 ). Autrement dit la matrice Vk = v1 , . . . , vk est orthogonale (on prendra garde
que c’est une matrice rectangulaire)
VkT Vk = I.
L’algorithme 2.1 permet également de définir une matrice Hk ∈ Rk+1×k par

(2.4) (Hk )ji = (Avj , vi ), j = 1, . . . i.

La matrice Hk est formée par colonnes dans l’algorithme. Elle a une structure particulière,
appelée forme de Hessenberg : elle a une seule diagonale non nulle sous la diagonale principale.
Nous verrons que cette structure est la clé pour une mise en oeuvre efficace de GMRES.
L’algorithme d’Arnoldi peut s’écrire
k
X
hk+1,k vk+1 = Avk − hjk vj .
j=1

Il est possible, et utile, de réécrire l’ensemble de ces relations sous forme matricielle :

(2.5) AVk = Vk+1 Hk .

Illustrons la relation 2.5 la matrice Hk dans le cas k = 4

   
h11 h12 h13 h14
   
h h22 h23 h24 
   
i  21
 
 h i h 
 v1 v2 v3 v4 = v1 v2 v3 v4
   
 A v5  h32 h33 h34 
   
   
   h43 h44 
   
h54

L’équation 2.5, l’orthogonalité de Vk (et l’orthogonalité de vk+ à Vk ) permettent d’écrire

(2.6) Hk = VkT AVk .

23
 Nous pouvons maintenant simplifier le problème de moindres carrés définissant xk . Notons
tout d’abord que nous pouvons écrire
r0 = βVk+1 e1 ,
où e1 est le premier vecteur de la base canonique de Rk+1 , et β = kr0 k2 . Le résidu rk peut donc
se mettre sous la forme
rk = b − Axk = r0 − A(xk − x0 ) = r0 − AVk yk ,
pour un vecteur yk ∈ Rk , puisque les colonnes de Vk forment une base de Kk (r0 ). La relation 2.5
permet d’écrire
rk = Vk+1 (βe1 − Hk yk ),
et l’orthogonalité de Vk+1 montre alors que le problème de minimisation qui doit être résolu à
l’itération k de GMRES est :
(2.7) min kHk y − βe1 k2 ,
y∈Rk

qui est un problème de moindres carrés posé en dimension k. Une fois ce problème résolu, en
notant sa solution yk , l’itéré xk de GMRES s’en déduit par
xk = x0 + Vk yk .
On voit ici une différence importante avec le cas du chapitre précédent : pour calculer xk on
a besoin de tous les vecteurs de la base de l’espace de Krylov qui ont déjà été calculés. Ces
vecteurs sont aussi nécessaires pour l’algorithme d’Arnoldi : le coût d’une itération de GMRES
augmente avec k.

2.2.2 Résolution du problème de moindres carrés

Nous donnerons dans ce paragraphe quelques indications sur la résolution de (2.7). Signalons
tout de suite que, puisque la matrice Hk est rectangulaire, le « système » Hh y = βe1 a une
équation de plus que d’inconnues.
L’appendice A.4 contient des rappels sur la résolution numérique des problèmes de moindres
carrés. En appliquant ces résultats au cas qui nous intéresse (le rôle de la matrice A dans
l’appendice est ici joué par la matrice Hk , qui est de dimension (k+1 ) × k, nous obtenons le
résultat suivant :
Théorème 2.1. Il existe une matrice orthogonale Qk ∈ R(k+1)×(k+1) et une matrice triangulaire
supérieure inversible Rk ∈ Rk×k telles que
 
Rk
(2.8) Hk = Qk  
0
La solution du problème de moindres carrés (2.7) est alors obtenue par :
(2.9) Rk yk = βzk ,
 
zk
avec QTk e1 =   (zk ∈ Rk , ρk ∈ R), et le résidu est égal à
ρk
kβe1 − Hk yk k2 = β |ρk | .

24
 Le coût de cette factorisation est (au premier ordre) O(k 3 ). Par ailleurs, il est facile de voir,
en utilisant l’orthogonalité de Qk que R vérifie

HkT Hk = RkT RK ,

autrement dit Rk est le facteur de Cholesky de la matrice des équations normales, mais il a été
obtenu sans former cette matrice. Le coût de la résolution du système triangulaire 2.9 est en
O(k 2 ), donc très inférieur à la factorisation.
En remontant les différentes définitions, nous voyons que

kb − Axk k2 = kβe1 − Hk yk k2 = β |ρk | ,

ce qui veut dire que le résidu à l’itération k est disponible sans qu’il y ait besoin de calculer xk !

2.2.3 Conséquences de l’arithmétique flottante

L’algorithme d’Arnoldi présenté en 2.1 suppose que les calculs sont réalisés exactement.
Si l’on programme cet algorithme sur un ordinateur, en arithmétique flottante (celle utilisée
par la majorité des langages, comme C, C++, Matlab ou Python), avec une précision relative
d’environ 15 chiffres décimaux, on constate que les vecteurs produits ne sont pas orthogonaux. Ce
phénomène bien connu est l’instabilité de la méthode Gram–Schmidt. Heureusement, il existe un
remède simple, qui passe par une modification de l’algorithme : on orthogonalise immédiatement
le vecteur Avk contre chacun des vecteurs de la base (v1 , . . . , vk ). L’algorithme appelé « Gram–
Schmidt modifié » est présenté dans l’algorithme 2.2. On peut montrer que les deux version

Algorithme 2.2 Algorithme d’Arnoldi avec Gram–Schmidt modifié, version numériquement

stable
Données: k ∈ N, base orthonormée (v1 , . . . , vk ) de Kk (r0 ).
Résultats: Base orthonormée (v1 , . . . , vk+1 ) de Kk+1 (r0 )
v1 = r0 /β
wk+1 = Avk
pour j = 1, . . . , k faire
hjk = (wk+1 , vj )
wk+1 = wk+1 − hjk vj
fin pour
hk+1,k = kwk+1 k2
si hk+1,k = 0 alors stop
vk+1 = wk+1 /hk+1,k

sont mathématiquement équivalentes, mais ont des comportements différents en arithmétique

flottante.
En incorporant un critère d’arrêt, on peut maintenant écrire une version complète de la
méthode GMRES
Contrairement au gradient conjugué, le coût d’une itération de l’algorithme n’est pas constant,
ni en calcul ni en stockage, puisque chaque itération doit enrichir l’espace de Krylov, et que la
base orthonormée (v1 , . . . , vk ) doit être construite explicitement. Ainsi, à l’itération k de la
méthode, on doit calculer
— 1 produit matrice–vecteur (le coût dépend de la matrice) ;

25
Algorithme 2.3 GMRES, avec Gram-Schmidt modifié, et résolution du problème de moindres
carrés par QR
Données: Itéré initial x0 ∈ Rn , tolérance η ∈ R, nombre maximum d’itérations kmax .
Résultats: Solution approchée du système xk
r0 = b − Ax0 , β = kr0 k2 , v1 = rβ0 , ρ = β, k = 0
Tant que ρ > η kbk2 et k < kmax faire
k =k+1
// Arnoldi (Gram–Schmidt modifié), calcul de vk+1
// Mise à jour de Vk+1 et de Hk
wk+1 = Avk
pour j = 1, . . . , k faire
hjk = (wk+1 , vj )
wk+1 = wk+1 − hjk vj
fin pour
hk+1,k = kwk+1 k2
si hk+1,k = 0 alors stop
vk+1 = wk+1 /hk+1,k
// Résolution du problème de moindres carrés miny∈Rk kβe1 − Hk yk2
Factoriser Hk = Qk Rk
e1 = (1, 0, . . . , 0)T ∈ Rk+1 , Qk e1 = (zkT , ρk )T
Résoudre Rk yk = βzk
ρ = βρk
fin Tant que
// Calcul de la solution à convergence
xk = x0 + Vk yk

— k produits scalaires de taille n

— k combinaisons linéaires de taille n
— la résolution du problème de moindres carrés, de coût O(k 3 ).

Remarque 2.1. La dernière étape (résolution du problème de moindres carrés) peut être amélio-
rée. On utilise la structure du problème (à chaque itération, on ajoute une colonne à la matrice
Hk ) pour mettre à jour itérativement la factorisation QR, ce qui ramène le coût de cette étape
à O(k). On utilise pour cela des rotations de Givens, voir les livres de Kelley [14, Par. 3.5] ou
Saad [20, Par. 6.5.3] pour des détails sur cette étape technique.

En ce qui concerne le stockage, k itérations de GMRES demandent le stockage de k vecteurs

de taille n (ainsi qu’une matrice de taille k × k dans notre implémentation simplifiée).
On voit donc que la principale limitation de GMRES est son coût, qui augmente plus que
linéairement avec le nombre d’itérations. Une manière de contrôler le coût, en particulier du
stockage, est de redémarrer la méthode après un nombre maximum d’itérations. Cette méthode,
nommée GMRES(m), où m est le paramètre de redémarrage est également décrite dans les livres
de Kelley [14, Alg. 3.5.2] ou Saad [20, Alg. 6.11].

26
2.3 Remarques sur la convergence
La convergence de GMRES est plus délicate à étudier que pour le gradient conjugué. La
principale raison est qu’il n’est pas toujours possible de diagonaliser une matrice quelconque, et
que même si c’est le cas, les valeurs propres ne jouent pas le même rôle que dans le cas défini
positif. Nous nous contenterons donc des résultats les plus simples et d’indiquer des limitations
connues.
Rappelons d’abord que, d’après la proposition 2.1, l’algorithme trouve toujours la solution
exacte du système (2.1) en au plus n itérations (ce résultat est valable en arithmétique exacte).
Comme pour le gradient conjugué, ce résultat n’a pas d’intérêt pratique, et ceci d’autant plus que
les itérations de GMRES sont coûteuses. Remarquons également que la quantité krk k2 décroît
au cours de l’algorithme (on minimise cette quantité sur des espaces de dimension croissante),
mais que cette décroissance n’est pas nécessairement stricte. Il existe des exemples pour lesquels
le résidu stagne (reste constant) jusqu’à la dernière itération. C’est la cas pour la matrice et le
second membre représentés ci–dessous pour n = 4 (mais l’exemple est valable pour n quelconque)
   
0 0 0 1 1
   
1 0 0 0 0
   
A= 
, b =  .
  
0 1 0 0 0
   
0 0 1 0 0

dont la solution exacte est x = (0, 0, 0, 1)T . Si l’algorithme démarre avec x0 = 0, le résidu reste
nul jusqu’à la n-ème itération, où la méthode trouve la solution.
Comme pour le gradient conjugué, mais avec des conséquences qui seront moins spectacu-
laires, nous pouvons exploiter le lien entre l’espace de Krylov et les polynômes pour obtenir le
résultat suivant, dont la démonstration vient simplement de (2.3) :
Théorème 2.2. Soit xk ∈ Kk (ro ) l’approximation obtenue après k itérations de GMRES. Pour
tout polynôme pk ∈ Pk , avec pk (0) = 1, on a
krk k2 ≤ kpk (A)r0 k2
et on en déduit la majoration (plus faible, mais indépendante du résidu initial)
krk k2
(2.10) ≤ kpk (A)k2 .
kr0 k2
Dans le cas où A est diagonalisable (rappelons que ce n’est pas vrai pour toute matrice),
nous obtenons le résultat suivant :
Théorème 2.3. Supposons que la matrice A soit diagonalisable
A = V ΛV −1
où Λ = diag(λ1 , . . . , λn ) est la matrice avec les valeurs propres de A sur la diagonale, et les
colonnes de V sont formées par les vecteurs propres de A. Alors pour tout polynôme p ∈ Pk ,
avec p(0) = 1, on a la majoration
krk k2
(2.11) ≤ κ2 (V ) max |p(z)| .
kr0 k2 z∈{λ1 ,...,λn }

27
Démonstration. Soit p un polynôme de degré k, avec p(0) = 1. On a

p(A) = V p(Λ)V −1 ,

(rappelons que les valeurs propres λ1 , . . . , λn , ainsi que les vecteurs propres, peuvent être com-
plexes). On peut majorer la norme de p(A) par

kp(A)k2 ≤ kV k2 V 1 2
kp(Λ)k2 ,

comme Λ est diagonale, on a encore kp(Λ)k2 = maxz∈{λ1 ,...,λn } |p(z)|. Par ailleurs on reconnaît
kV k2 V 1 2 = κ2 (V ), le conditionnement de la matrice des vecteurs propres. 

Remarque 2.2. Ce résultat appelle les commentaires suivants, qui réduisent son application en
pratique :
— Au moins pour les matrices diagonalisables, on peut donc se ramener à un problème de
minimisation posé sur les polynômes. Ce problème est toutefois plus difficile à résoudre
que celui rencontré pour les matrices symétriques et définies positives, parce qu’il est
posé sur un sous-ensemble a priori quelconque du plan complexe. On peut se ramener à
une partie convexe, en pratique des résultats existent si le spectre de A peut être inclus
dans une ellipse qui ne contient pas l’origine (cette dernière hypothèse est restrictive),
voir Saad[20, Par. 6.11.4] pour ces résultats.
— La quantité κ2 (V ) qui apparaît dans la borne (2.11) est difficile à estimer, et il n’existe
pas de borne générale a priori (rappelons que dans le cas des matrices symétriques et
définies positives, la matrice V peut être choisie orthogonale, ce qui conduit à κ2 (V ) = 1.
En général, même si l’on sait minimiser la quantité maxz∈{λ1 ,...,λn } |p(z)| (et le point
précédent montre les limites de cette question), la borne obtenue est multipliée par le
conditionnement des vecteurs propres, et ne donne donc qu’une indication très faible sur
la convergence réelle de GMRES.
Pour ces raisons, on a cherché d’autres outils pour étudier la convergence de GMRES. On pourra
consulter Greenbaum [12, Par. 3.2] pour une discussion des résultats connus.
Pour finir, signalons le résultat suivant, du à Greenbaum, Pták et Strakoš [13] : étant donné
une suite (finie) décroissante η1 ≥ · · · ≥ ηn , il existe une matrice A et un second membre b
tels que l’application de GMRES au système Ax = b pour résidus les nombres ηn . De plus, les
valeurs propres de la matrice peuvent être choisis à l’avance. Autrement dit, en toute généralité,
la connaissance des valeurs propres ne donne aucune indication sur la convergence.

28
Chapitre 3

Préconditionnement

Ce chapitre donne quelques indications sur les les principales techniques de préconditionne-
ment. Le terme « préconditionnement » vient de la réalisation que (au moins pour le gradient
conjugué, c’est moins vrai pour GMRES), la convergence de la méthode est gouvernée par le
conditionnement du système linéaire (voir en particulier la proposition 1.5). L’idée du précondi-
tionnement est alors de remplacer le système à résoudre par un autre, qui possède bien entendu
la même solution, mais dont la matrice soit « mieux conditionnée ». Ce terme est ici mis entre
guillemets car, comme nous l’avions signalé au chapitre 1, la convergence du gradient conjugué
dépend en réalité de la répartition de l’ensemble des valeurs propres, pas seulement du rapport
des valeurs propres extrêmes.
Après quelques remarques générales sur ce qu’est une méthode de préconditionnement, nous
verrons le cas particulier du gradient conjugué. Puis nous donnerons quelques exemples de mé-
thodes générales. Ce chapitre n’est qu’une brève introduction à un vaste sujet. Une référence
pour approfondir est le livre de Y. Saad [20], qui consacre deux chapitres au sujet, et deux autres
pour étudier en détail des méthodes qui s’appliquent particulièrement bien aux systèmes issus
de la discrétisation des équations aux dérivées partielles. Le livre (accessible en ligne) [4] donne
des indications pratiques sur la mise en oeuvre des méthodes. Chen [7] est entièrement consacré
à ce sujet, tandis que d’autres livres y consacrent une partie importante de leur contenu. On
peut citer ceux de O. Axelsson [3], Axelsson et Barker [2], A. Greenbaum[12], G. Meurant[18]
et H. van der Vorst [23]. Enfin, des articles de revue font régulièrement le point sur les progrès
récents, citons par exemple celui, un peu plus ancien, de Benzi[5] et celui, récent, de Wathen[24]

3.1 Généralités
Étant donné un système linéaire

Ax = b, A ∈ Rn , b ∈ Rn

le but d’une méthode de préconditionnement (on parle de manière quelque peu impropre de
« préconditionneur ») est de remplacer ce système par un autre (qui a la même solution) est
qui soit « plus facile » à résoudre. Il est difficile, et peu utile, de définir précisément ce que
l’on entend par « plus facile ». En pratique, on cherche à diminuer le nombre d’itérations de la
méthode originale, mais nous verrons que le préconditionnement a en général un sur–coût, et
que le but sera finalement de réduire la durée totale de résolution.

29
 Pour préciser quelque peu, une méthode de préconditionnement (à gauche, nous verrons qu’il
y a d’autres possibilités) est donnée par une matrice M inversible, la matrice de précondition-
nement, et le système préconditionné est

(3.1) M −1 Ax = M −1 b.

Comme toujours, il est important de noter que l’écriture M −1 est symbolique. Il n’est jamais
nécessaire de former explicitement l’inverse de M . Nous verrons qu’en général, la méthode se
traduit par la nécessité de résoudre, à chaque itération de la méthode itérative, un système
linéaire M z = r. Le but est alors que
i) Le nombre d’itérations de la méthode avec préconditionnement soit plus faible que le
nombre d’itérations de la méthode originale ;
ii) le sur–coût de chaque itération, induit par la résolution de l’étape de préconditionnement,
reste faible.
Il faut ajouter ce que l’on appelle le coût de mise en place du préconditionnement, c’est-à-dire
le coût du calcul de la matrice M à partir de A.
Il est habituel de dire que la matrice de préconditionnement est M , alors que on doit agir
avec M −1 . La raison est que c’est M qui doit « ressembler » à A.
Deux exemples extrêmes aident à comprendre le compromis à réaliser :
— Choisir M = A−1 . Dans ce cas, le nombre d’itérations est réduit à 1, mais il est bien
évident que l’on n’a rien gagné, puisque l’étape de préconditionnement est aussi difficile
que le problème original ;
— Choisir M = I. On ne fait rien, ce qui ne coûte rien, mais ne change pas le nombre
d’itérations.
Si un compromis « optimal » existe, il doit se situer entre ces deux extrêmes. Mais notons que ce
compromis ne peut être absolu, et dépendra de la nature du problème : difficulté de la résolution
sans préconditionneur, coût d’une itération (pour le gradient conjugué, c’est essentiellement le
coût de la multiplication de la matrice A par un vecteur, pour GMRES, le coût d’une itération
n’est pas constant, il est encore plus important de minimiser la taille de l’espace de Krylov). Les
méthodes de préconditionnement ont donc souvent une composante heuristique. Il est parfois
possible de justifier certains choix, mais ce sera rarement le cas dans le cadre de ce cours.
Pour conclure ces généralités, notons que nous avons insisté sur la présentation algébrique,
mais qu’une méthode de préconditionnement peut souvent provenir d’une méthode de résolu-
tion approchée d’un problème voisin, l’interprétation algébrique est alors simplement un cadre
général.

3.1.1 À gauche, à droite, ou des deux cotés ?

Il y a plusieurs façons de préconditionner un système linéaire. Nous avons introduit précé-
demment le préconditionnement à gauche, qui consiste à remplacer le système original par

(3.2) M −1 Ax = M −1 b.

Cela veut dire qu’à chaque étape d’une méthode de Krylov, on doit remplacer le produit r = Av
par la succession de
r = Av, suivi de la résolution M z = r.

30
Selon les méthodes et les situations, cette vision en deux étapes est plus ou moins explicite.
Avant la résolution, on doit calculer le nouveau second membre M −1 b, mais souvent, cela revient
simplement à préconditionner le résidu initial M −1 (b − Ax0 ).
On peut également considérer le système préconditionné par la droite :

(3.3) AM −1 y = b, avec M x = y,

on doit donc cette fois résoudre une étape de préconditionnement à la fin des itérations pour
retrouver la solution du problème original. Nous verrons que, au moins dans le cas de GMRES,
la méthode avec préconditionnement peut être mise en oeuvre sans faire référence à la variable
auxiliaire y.
Enfin, les deux idées précédents peuvent être combinées. On prend deux matrices M1 et M2
et on résout le système

(3.4) M1−1 AM2−1 y = M1−1 , M2 x = y.

L’approche bilatérale est parfois utile quand la matrice de préconditionnement est obtenue sous
la forme d’un produit :
M = LU,
où L est U peuvent être (mais ne sont pas obligatoirement) des matrices triangulaires, voir le
paragraphe 3.3.3. Dans le cas du gradient conjugué, cette approche est nécessaire pour maintenir
la symétrie du système préconditionné.
Quelles sont les similarités, et les différences, entre le préconditionnement à gauche et à
droite ? Notons tout d’abord que les matrices M −1 A et AM −1 sont semblables, et ont donc les
même valeurs propres (dans les deux cas ce sont les valeurs propres du problème généralisé Ax =
λM x). Si la convergence est gouvernée par les valeurs propres, les deux approches donneront
une convergence similaire.
Dans le cas d’un système non–symétrique résolu par GMRES, on démontre (voir Saad [20,
Prop. 9.1], que la méthode avec préconditionnement à droite minimise encore le résidu de la
méthode originale, alors que la méthode avec préconditionnement à gauche minimise la norme
du résidu préconditionné. Les deux quantités sont souvent proches, sauf dans le cas où la matrice
de préconditionnement M est mal conditionnée.
L’algorithme GMRES avec préconditionnement est presque identique à l’algorithme sans
préconditionnement. La seule différence est de remplacer dans l’algorithme 2.3, la ligne

wk+1 = Avk

par
wk+1 = M −1 Avk pour un préconditionnement à gauche,
ou par
wk+1 = AM −1 vk pour un préconditionnement à droite.

3.2 Préconditionnement du gradient conjugué

Le cas du gradient conjugué demande une attention particulière, et nous verrons que la mise
en oeuvre de la méthode est particulièrement simple.

31
 Rappelons que le gradient conjugué ne peut fonctionner qu’avec une matrice symétrique et
définie positive. A première vue, ni le préconditionnement à droite, ni à gauche ne sont possibles.
Nous nous placerons donc dans le cas d’un préconditionnement bilatéral, où la matrice M est
décomposée sous la forme M = LLT (la matrice L sera souvent obtenue par une factorisation
incomplète, voir paragraphe 3.2, mais cette hypothèse n’est pas nécessaire ici). Si la matrice L
est inversible (ce qui l’hypothèse minimale), alors M sera définie positive, et inversement, toute
matrice définie positive peut se mettre sous la forme M = LLT . Le système préconditionné est
alors

(3.5) L−1 AL−T x̂ = L−1 b, LT x = x̂,

et nous noterons avec des chapeaux toutes les quantités transformées :

 = L−1 AL−T , b̂ = L−1 b.

La matrice  est symétrique et définie positive. Mais nous allons voir que l’introduction de ces
quantités n’est utile qu’au niveau théorique, et que l’algorithme du gradient conjugué précondi-
tionné par M peut se réécrire entièrement avec les données originales, et que seule la matrice M
est nécessaire.
Pour voir cela, définissons les quantités liées au système préconditionné :

x̂k = LT xk , p̂k = LT pk ,

Pour le résidu, on a

r̂k = b̂ − Âx̂k = L−1 b − L−1 AL−T LT xk = L−1 rk ,

et il sera utile d’introduire le vecteur zk défini par

r̂k = LT zk , de sorte que L−1 rk = LT zk , soit M zr = rk .

Nous calculons maintenant le scalaire

(r̂k r̂k , ) (LT zk , LT zk ) (zk , rk )
α̂k = = = ,
(Âp̂k , p̂k ) (L AL−T LT pk , LT pk )
−1 (Apk , pk )

de sorte que l’on a

x̂k+1 =x̂k + α̂k p̂k =⇒ xk+1 =xk + α̂k pk

(3.6)
r̂k+1 =r̂k − α̂k Ap̂k =⇒ rk+1 =rk − α̂k Apk .

De même,
(r̂k+1 , r̂k+1 ) (zk+1 , rk+1 )
β̂k = = ,
(r̂k , r̂k ) (zk , rk )
et

(3.7) p̂k+1 = r̂k + β̂k p̂k =⇒ pk+1 = L−T r̂k + β̂k pk = zk + β̂k pk .

La conclusion de ces calculs est bien que l’algorithme du gradient conjugué préconditionné
peut être mis en oeuvre en ajoutant simplement une étape supplémentaire à chaque itération de
l’algorithme, pour calculer le vecteur zk , appelé le résidu préconditionné. L’algorithme obtenu

32
Algorithme 3.1 Algorithme du gradient conjugué
Données: b ∈ Rn , x0 ∈ Rn , A ∈ Rn×n
Résultats: k, xk
r0 = b − Ax0 , z0 = M −1 r0 , p0 = r0 , k = 0
Tant que k ≤ itmax et pas de convergence faire
(rk , zk )
αk =
(pk , Apk )
xk+1 = xk + αk pk
rk+1 = rk − αk Apk
zk+1 = m−1 rk+1
(rk+1 , zk+1 )
βk =
(rk , zk )
pk+1 = zk+1 + βk pk
fin Tant que

est le suivant, à comparer avec l’algorithme (nous avons noté sans chapeau les deux scalaires αk
et βk ) :
Il n’est pas inutile d’insister sur le fait que, comme toujours, la notation zk+1 = M −1 rk+1
est un raccourci mathématique, et que l’implémentation revient à résoudre, à chaque itération
un système
M zk+1 = rk+1
pour le résidu préconditionné. La matrice de préconditionnement aura été choisie de manière à ce
que ce système soit « facile à résoudre ». Nous présentons plusieurs méthodes au paragraphe 3.3.

3.3 Exemples de préconditionnement

Nous donnons dans ce paragraphe quelques exemples de préconditionnement, parmi les mé-
thodes les plus simples. Nous ne pourrons ici encore, et plus qu’ailleurs qu’effleurer la surface
d’un domaine très vaste, et faisant encore l’objet de recherches actuelles.
Avant de décrire quelques méthodes, notons que l’on peut grossièrement classer les méthodes
de préconditionnement en deux catégories (la réalité est comme toujours plus subtile) :
— des méthodes généralistes, faisant peu d’hypothèses sur la matrice, ou l’origine du pro-
blème ;
— des méthodes plus spécialisées, s’appliquant à une catégorie restreinte de problèmes.
Dans ce cours, nous nous bornerons à quelques exemples de la première catégorie (utilisation
de méthodes itératives, factorisation incomplètes, préconditionnement polynomial). Mais il est
assez clair que plus on utilise une information spécifique au problème à résoudre, plus on a de
chances de concevoir un préconditionneur efficace. On consultera les références données au début
de ce chapitre pour découvrir certaines des méthodes qui ont été proposées.

3.3.1 Préconditionnement diagonal

Cette méthode n’est certainement pas la plus efficace, nous la mentionnons en premier pour
sa simplicité, et parce qu’elle est très généralement applicable.
Si A est une matrice inversible, notons D sa diagonale principale. Nous faisons l’hypothèse
que D est inversible (ou de manière équivalente que Dii 6= 0, f oralli. Cette hypothèse est en

33
particulier vérifiée si A est définie positive (pourquoi ?), et elle le sera très souvent pour des
matrices provenant de la discrétisation de modèles d’EDP. Un cas important où elle n’est pas
vérifiée es celui des matrices de forme « point selle » :
 
P Q
A= 
R 0
qui interviennent dans les problèmes d’optimisation avec contraintes, et également dans les
problèmes de Stokes ou Navier–Stokes (on a le plus souvent P = RT ), voir ??.
Le préconditionnement diagonal consiste dans le choix M = D. L’étape de mise en place
initiale est ici triviale, et l’étape de préconditionnement consiste simplement à effectuer une
boucle sur toutes les inconnues : résoudre M z = r est équivalent à
ri
(3.8) zi = , ∀i = 1, . . . , n.
Dii
Si la matrice A est symétrique et définie positive, ses éléments diagonaux sont tous strictement
positifs, ce qui fait que nous sommes dans la situation de la Section 3.2.
L’intérêt de ce préconditionneur est qu’il conduit à une mise à l’échelle des coefficients de la
diagonale de A qui peut parfois améliorer le comportement des algorithmes. Mais cette remarque
n’est pas toujours vérifiée (exemple du problème de Poisson, la diagonale est égale à 4 fois
l’identité). On peut montrer que dans certains cas, la diagonale de A est la matrice diagonale
qui réduit le plus le conditionnement parmi toutes les matrices diagonales. Mais cette réduction
peut être très modeste (voir Meurant [18, Sec. 8.2] ou Greenbaum [12, Sec. 10.5] pour plus de
détails).

3.3.2 Préconditionnment par des méthodes itératives classiques

Nous allons voir que les méthodes itératives dites « classiques », comme Jacobi, Gauss–Seidel
ou SSOR (voir Section A.5) conduisent à des préconditionnement.
Nous décomposons la matrice A sous la forme
A=M −N
où M sera supposée inversible. On obtient une classe de méthodes itératives en résolvant le
système équivalent x = (M −1 N )x+M −1 b par une méthode de point fixe, qui conduit à l’itération
(3.9) xk+1 = (M −1 N )xk + M −1 b, x0 donné.
Pour mettre en évidence le lien avec les méthodes de préconditionnement, nous éliminons la
matrice N en remarquant que N = M − A, et donc l’itération (3.9) est équivalente à
(3.10) xk+1 = xk + M −1 (b − Axk ),
et l’on reconnaît une itération de point fixe pour résoudre le système préconditionné (3.1). Plutôt
que d’utiliser l’itération simple (dont la convergence est en général lente), on résout le système
préconditionné par une méthode de Krylov (gradient conjugué ou GMRES), en utilisant la
matrice M comme matrice de préconditionnement. À chaque itération de la méthode de Krylov,
on effectue une seule itération de la méthode classique, en résolvant le système M z = rk .
Les méthodes classiques correspondent à des décompositions particulières. Suivant la tradi-
tion, nous notons A = D − E − F , avec D diagonale, E la sous–diagonales stricte de A, et F la
sur–diagonale (stricte) de A. Si A est symétrique, on a F = E T .

34
 Jacobi On prend M = D, on retrouve le préconditionnement diagonal du paragraphe pré-
cédent ;
Gauss-Seidel On prend M = D − E. Le système M z = r est facile à résoudre, puisque
la matrice M est ici triangulaire inférieure. Notons que si la matrice A est symétrique
et définie positive, on doit utiliser la version symétrique de la méthode de Gauss–Seidel.
Elle correspond au choix ω = 1 du point suivant ;
1
SOR Dans le cas où la matrice A n’est pas symétrique, on prend M = (D − ωE), M reste
ω
triangulaire inférieure. Dans le cas symétrique, on doit utiliser la version symétrique de
SOR (SSOR), et on a alors :

M = (D − ωE)D−1 (D − ωE T ).

(nous avons supprimé le facteur ω(2 − ω) puisque le préconditionneur n’est pas sensible
à un changement de facteur d’échelle). On constate que la matrice M est un produit
de matrices triangulaires (et diagonale), donc toujours facile à inverser. En détail, la
résolution de M z = r pour le préconditionnement SSOR s’effectue de la manière suivante :

résoudre (D − E) u = r
(3.11) v = Du
T
résoudre (D − E) z = v

Il existe peu de résultats théoriques concernant ce type de préconditionnement. Dans le

livre[12, Sec. 10.4], A. Greenbaum montre comment généraliser un résultat de comparaison des
rayons spectraux pour comparer des conditionnements.

3.3.3 Factorisations incomplètes

Les méthodes de factorisation incomplète restent les plus utilisées parmi les méthodes géné-
ralistes, c’est-à-dire applicable uniquement avec la connaissance de la matrice.
Le point de départ de ces méthodes est la constatation que la factorisation LU (ou de
Cholesky) d’une matrice creuse contient en général beaucoup plus d’éléments non–nuls que la
matrice de départ. On peut alors chercher à calculer une factorisation approchée, dont les facteurs
L et U contiennent moins d’éléments non–nuls que les « vrais » facteurs. Il existe plusieurs
variantes de la méthode, avec de degrés de sophistication croissants. Nous nous contenterons
comme d’habitude de présenter la variante la plus simple, appelée ILU(0).

Factorisation incomplète ILU(0)

Dans cette méthode, les matrices L et U (comme d’habitude L est triangulaire inférieure à
diagonale unité, et U est triangulaire supérieure) remplissent les conditions suivantes :
(
lij = uij = 0 si aij = 0
(3.12)
(LU )ij = aij si aij 6= 0.

Notons bien la différence avec une factorisation complète : la matrice A n’est pas égale au
produit LU . Il existe une matrice R (avec Rij = 0 si aij 6= 0) telle que

A = LU + R.

35
Par contre, on contrôle a priori la mémoire utilisée par la méthode.
On peut obtenir les facteurs incomplets de A par une modification simple de l’algorithme de
factorisation de Gauss (sans pivotage). Dans l’algorithme 3.2, les seules modifications sont les
deux tests « si aij 6= 0 ».

Algorithme 3.2 Factorisation incomplète ILU(0)

Données: A ∈ Rn×n
Résultats: Facteurs L et U de la factorisation incomplète, vérfiiant (3.12).
// La factorisation est « sur place » : les matrices L et U écrasent A
pour k = 1, . . . , n − 1 faire
pour i = k + 1, . . . , n faire
si aij 6= 0 alors aik = aik /akk
pour j = k + 1, . . . , n faire
si aij 6= 0 alors aij = aij − aik akj
fin pour
fin pour
fin pour

Il n’est pas évident que l’algorithme 3.2 permette toujours de calculer les matrices L et U ,
autrement dit que les conditions (3.12) déterminent ces matrices. Et en fait, ce n’est pas toujours
possible, et l’algorithme 3.2 n’est pas un algorithme sens strict, il ne se termine pas toujours.
Nous présentons donc une condition suffisante relativement simple pour que l’algorithme 3.2
se termine. On pourra consulter [20, Sec. 10.3], [23, Chap. 13] pour plus de détails.

Définition 3.1. Une matrice (inversible A) est une M-matrice si et seulement si elle vérifie les
conditions suivantes :
i) aij ≤ 0, ∀i 6= j ;
ii) tous les éléments de A−1 sont positifs.

Remarque 3.1. La condition ii) ci–dessus est équivalente à

Si x est la solution de Ax = b, et si bi ≥ 0, ∀i, alors xi ≥ 0, ∀i.

Il s’agit d’une forme discrète du principe du maximum.

Exemple 3.1. — Si A est une matrice symétrique et définie positive, dont tous les éléments
hors–diagonaux sont négatifs, alors A est une M-matrice.
— Si A est une matrice avec aii > 0, ∀i et aij ≤ 0, ∀i 6= j, telle que
X
(3.13) |ai | > |aij | ,
i6=j

alors A est une M-matrice.

La condition (3.13) exprime que A est à diagonale strictement dominante

On démontre alors le résultat suivant (voir [20, Thm. 10.2] ou[23, Thm. 13.3] pour la dé-
monstration) :

36
Théorème 3.1. Soit A ∈ Rn×n une M-matrice (inversible). L’algorithme 3.2 calcule des ma-
trices L (triangulaire inférieure, à diagonale unité) et U (triangulaire supérieure) vérifiant les
conditions (3.12).
Dans le cas où A est symétrique et définie positive, il est naturel de chercher une matrice L
(triangulaire inférieure, avec des éléments diagonaux positifs) telle que
(
lij = 0 si aij = 0
(3.14) T
(LL )ij = aij si aij 6= 0.

La factorisation obtenue s’appelle une factorisation de Cholesky incomplète, et existe si A est

une M-matrice. L’algorithme 3.2 peut être adapté à cette situation.
Nous verrons que, dans certaines situations, il peut être plus commode de modifier la facto-
risation de sorte que l’on ait

A = (L + D)D−1 (L + D)T + R

où les éléments diagonaux de L sont cette fois égaux à 1 (D est toujours la diagonale de A).

Un exemple : factorisation incomplète pour le problème modèle

Nous présentons dans ce paragraphe un exemple de calcul de factorisation incomplète. La
motivation vient du problème modèle, la discrétisation par différences finies de l’équation de
Poisson sur un carré, avec un maillage régulier de pas h. Pour simplifier, nous supposons que le
schéma utilisé conduit à une matrice pentadiagonale, comme à l’exemple 1.1, ou la figure 3.1.

 
  
    
     
   
 
 

    

    
    
 
 
A=
      

      

    
    
   
 
 

     

   
  

Figure 3.1 – Un exemple de matrice penta–diagonale

Nous décomposons la matrice A sous la forme

(3.15) A = L + D + LT

où D est la diagonale de A et L sa partie strictement triangulaire inférieure. Pour le problème

modèle, tout les éléments diagonaux sont égaux à 4, et tous les éléments hors–diagonaux sont
0. L n’est donc pas la même matrice que celle de (3.12)

37
égaux à −1. Cependant, la valeur de ces éléments ne joue aucun rôle, seule la structure de la
matrice A est importante, c’est-à-dire que nous considérons une matrice A décomposée comme
en (3.15), où la matrice L n’a que deux diagonales non–nulles. Ces diagonales sont situées à une
distance de 1 et de n de la diagonale principale (n est le nombre de points de grilles, sans compter
les bords, la matrice A est de taille n2 × n2 . Nous ferons l’hypothèse que A est symétrique et
définie positive (c’est le cas pour le problème modèle), et nous chercherons donc une factorisation
de Cholesky incomplète, notée IC(0).
Pour cet exemple, il est commode de modifier légèrement la forme de la factorisation incom-
plète. Nous cherchons la matrice de préconditionnement M sous la forme factorisée
M = (LM + ∆)∆−1 (LM + ∆)T ,
ou ∆ est diagonale et LM strictement triangulaire inférieure.
La condition que M n’aie d’éléments non–nuls qu’aux mêmes places que A force LM a n’avoir
que deux diagonales non nulles, également aux distances 1 et n de la diagonale principale. Donc
sur une ligne 1 ≤ i ≤ n2 donnée, les seuls éléments non nuls de LM sont ceux situés en (i, i − 1)
et (i, i − n). Nous adoptons pur la suite de cet exemple la convention qu’un indice nul ou négatif
n’est pas pris en compte.
En effectuant le produit des matrices définissant M nous obtenons
M = LM + LTM + ∆ + LM ∆−1 LTM .
Les calculs se simplifient si les choses sont prises dans le bon ordre. Nous cherchons tout
d’abord la structure des différents termes, sans nous soucier de calculer les valeurs.
— C’est simple pour les 3 premiers termes : le premier et le second n’ont, par construction,
que deux diagonales non nulles, aux mêmes places que L. Et le troisième est diagonal ;
— L’élément générique du dernier terme (appelons H cette matrice) est :
n 2
X
(3.16) hij = (LM )ik (LM )jk ∆−1
kk .
k=1

Mais du fait de la structure de LM , k ne peut prendre qu’un nombre très restreint de

valeurs : on doit avoir k = i − 1 ou k = i − n, et symétriquement, k = j − 1 ou k = j − n.
Pour que ces conditions soient compatibles, il faut que, soit j = i, soit j = i − n + 1 (avec
j ≤ i). Par conséquent, les termes non nuls de H sont, soit sur la diagonale principale,
soit sur les diagonales situés à une distance de n − 1 de la diagonale principale.
Nous voyons maintenant comment les différents termes du produit définissant M contribuent.
Les deux diagonales à la distance de n − 1 de la diagonale principale n’existent pas dans A, et ne
vont donc pas contribuer au préconditionneur. C’est la deuxième condition (3.12) qui l’impose.
Ces deux termes constituent les éléments de la matrice R définie par A = M + R.
Ensuite, les seuls éléments venant se placer sur les diagonales 1 et n sont ceux provenant de
LM . Par conséquent (toujours à cause de la deuxième condition (3.12), nous devons avoir
LM = L,
et ceci sans aucun calcul. C’est ce qui explique la forme choisie pour la factorisation.
Enfin, il nous reste à déterminer la matrice diagonale ∆. Pour cela nous exploiterons l’égalité
des éléments diagonaux : Mii = Aii , pour 1 ≤ i ≤ n2 . D’après (3.16), les éléments diagonaux de
H sont donnés par
2 2
hii = li,i−1 /∆i−1,i−1 + li,i−n /∆i−n,i−n ,

38
et on doit donc avoir (en notant simplement ∆i l’élément diagonal de la ligne i.
2 2
(3.17) ai,i = ∆i + li,i−1 /∆i−1 + li,i−n /∆i−n, ,
ce qui donne une relation de récurrence pour ∆i,i (l’initialisation se fait en prenant en compte
que les éléments d’indice négatifs ou nuls sont absents). En raisonnant sur la grille (avec une
numérotation par ligne de gauche à droite et de bas en haut), on voit que l’élément de ∆
correspondant à un sommet se calcule en fonction de quantités associées aux points directement
en bas, et à gauche.
On pourra trouver dans [20, Sec. 10.3.4] une approche plus générale, qui se base directement
sur la grille de différences finies.

Généralisations
La méthode de base que nous avons exposée peut être généralisée dans différentes directions :
— On peut autoriser plus d’éléments non–nuls dans L et U que dans A. Pour cela, on
se fixe un sous–ensemble S ⊂ {(i, j), i 6= j, 1 ≤ i, j ≤ n} de positions où l’on impose
lij = uij = 0 (on exclut toujours la diagonale). Le cas particulier vu précédemment est
S = {(i, j) aij = 0}. Dans l’algorithme 3.2 il suffit alors de remplacer la condition aij 6= 0
par (i, j) 6∈ S. Le théorème 3.1 s’étend également, sous les mêmes hypothèses ;
— On peut également décider au cours du calcul quels éléments seront conservés, par exemple
sur un critère de taille (factorisation avec seuil) ;
— on peut sommer tous les éléments que l’on ne conserve pas, et ajouter la somme à la
diagonale de L ou de U . On parle de factorisation modifiée ;
— si la factorisation échoue (lorsque la matrice A n’est pas une M–matrice), on peut ajouter
à A un multiple de l’identité, en espérant que la factorisation aboutisse. La justification
est que de toutes façons la méthode repose sur des heuristiques.
Pour ces différents aspects, on pourra consulter des références plus spécialisées, par ordre
de difficulté approximativement croissante : [20, Chap. 10], [23, Chap. 13], [18, Chap. 8], [3,
Chap. 7].

3.3.4 Préconditionnement polynomial

Une autre classe assez générale de préconditionneurs est le préconditionnement polynomial.
On cherche ici une matrice de préconditionnement M de la forme
M = p(A),
où p est un polynôme. Le système préconditionné est donc
p(A)Ax = p(A)b
et l’on remarque que comme p(A) et A commutent, il est indifférent que le préconditionnement
soit appliqué à droite ou à gauche.
Une première idée est d’utiliser la série de Neumann. Rappelons le résultat suivant :
Lemme 3.1. Soit G une matrice vérifiant kGk < 1 pour une norme subordonnée (cf. Lemme A.1).
On a l’égalité
∞
X
(3.18) (I − G)−1 = Gk = I + G + G2 + . . . + Gk + . . .
k=0

39
où la série est absolument convergente dans L(Rn ).
Il est alors naturel de tronquer la série, et de proposer la somme Gk = I + G + . . . + Gk
comme approximation de l’inverse de G.
Pour appliquer cette idée à une matrice quelconque A, il est nécessaire d’introduire un facteur
de mise à l’échelle ω. La remarque simple que (on désigne toujours par D la diagonale de A)
ωA = D − (D − ωA)
montre que (prendre l’inverse des deux membres pour vérifier cette égalité) :
−1 −1
(ωA)−1 = I − (I − ωD−1 A)

D
En appliquant le Lemme 3.1 à la matrice
G = I − ωD−1 A,
on définit donc M −1 par
(3.19) M −1 = (I + G + . . . + Gk )D−1 .
Puisque par définition D−1 A = 1/ω(I − G), on peut vérifier que
1
M −1 A = (I − Gk+1 ).
ω
L’évaluation de Gk+1 peut être numériquement délicate, la formule du binôme est à éviter !
Nous supposons maintenant que A est symétrique et définie positive, et nous noterons 0 <
λ1 ≤ λ2 ≤ . . . ≤ λn ses valeurs propres. Les valeurs propres de la matrice p(A)A sont donc
les nombres λi p(λi ) mais l’ordre n’est pas nécessairement préservé. L’utilisation de la série de
Neumann revient à privilégier l’information spectrale autour de 0. Il peut être préférable de
considérer le spectre de p(A)A dans son ensemble, et de chercher un polynôme p qui rende le
spectre aussi proche de 1 que possible. Notons q(λ) = λp(λ), et remarquons qu’on doit avoir
q(0) = 0. Comme lorsque nous avons étudié la convergence du gradient conjugué (Section 1.4),
il n’est pas possible de travailler sur le spectre discret, nous ferons donc l’hypothèse que nous
connaissons un intervalle [a, b], avec a > 0, contenant toutes les valeurs propres de A (cette
hypothèse n’est pas particulièrement réaliste). On cherche donc à minimiser l’erreur k1 − qk, sur
l’ensemble des polynômes de degré fixé (mais pas connu a priori), tels que q(0) = 0. Les deux
choix les plus courants sont :
le polynôme du minimax on prend
k1 − qk∞ = max |1 − q(λ)| .
[a,b]

On montre que la solution est donnée encore une fois par un polynôme de Tchebychev
(voir [20, Sec. 12.3] ou [18, Sec. 8.18] pour les détails ;
le polynôme des moindres carrés on prend
Z b
k1 − qk22 = |1 − q(λ)|2 dλ.
a
Il existe une description explicite du polynôme optimal, mais il peut aussi être intéressant
de résoudre le problème de moindres carrés numériquement. Nous renvoyons également
à [20, Sec. 12.3] ou [18, Sec. 8.18] pour les détails.

40
 L’avantage du préconditionnement polynômial est qu’il se prête facilement à une mise en
oeuvre parallèle, puisque tout ce qui est demandé est de multiplier la matrice A par un vec-
teur. Ceci est en contraste avec les méthodes de factorisation incomplètes, qui demandent la
résolution de systèmes triangulaires, une étape qui est intrinsèquement séquentielle. Le principal
inconvénient de ces méthodes est qu’elles reposent sur une estimation des valeurs propres de
A, et qu’une telle estimation est difficile à obtenir. Des méthodes adaptatives (où l’estimation
des valeurs propres est améliorée au cours des itérations) existent, mais leur implémentation est
complexe.
Pour conclure, rappelons comment on peut évaluer le résidu préconditionné. Pour calculer
i = 0k αi λi , il ne faut surtout pas calculer les puissances de A, il est
P
z = p(A)r, où p(λ) =
recommandé d’utiliser la méthode Horner, comme dans l’algorithme 3.3 (cf. [18, 8.18.4]).
Pk ir
Algorithme 3.3 Algorithme de Horner pour le calcul de z = i=0 αi A
bk = αk r
pour i = k − 1, . . . , 0 faire
bi = αi r + Abi+1
fin pour
z = b0

41
Chapitre 4

Matrices creuses

Nous donnons dans ce chapitre quelques indications sur les structures de données permettant
de stocker efficacement des matrices creuses, et les algorithmes pour exploiter ces structures de
données.

4.1 Introduction
Il n’existe pas de définition mathématique d’une matrice creuse, mais on peut donner une
« définition » opérationnelle :

Définition 4.1. Une matrice creuse est une matrice dont un très grande partie des éléments sont
nuls, et pour laquelle une représentation informatique qui ne stocke que les éléments non–nuls
conduit à un gain en performance, en mémoire ou en temps de calcul.

Cette définition reste imprécise sur « la très grande partie », et sur la manière d’exploiter
l’information. Il est clair que le stockage seul ne suffit pas, et que l’exploitation pour utiliser la
matrice dans un code de calcul est critique.
De manière pragmatique, on sait que certaines applications conduisent à des matrices creuses,
et on sait stocker ces matrices de manière à pouvoir réaliser les opérations les plus communes.
On doit donc se poser non–seulement la question de comment stocker de manière économique ces
matrices, mais également celle des algorithmes permettant de réaliser les opérations nécessaires.
L’une des applications rencontrée le plus fréquemment est la discrétisation d’équations aux
dérivées partielles par des méthodes comme les différences finies, les éléments finis ou les vo-
lumes finis. Dans les trois cas, on obtient une matrice creuse car la discrétisation est locale.
Les inconnues en un point ne sont reliées qu’à un petit nombre d’inconnues, en général géomé-
triquement proches. Dans le cas des différences finies, ces inconnues correspondent aux voisins
immédiats dans les directions d’espace (parfois à des voisins à deux niveaux). Dans le cas des
éléments finis, les inconnues en un point géométrique ne sont reliées qu’à un petit nombre de
celles coreepondant à des points partageant un élément avec le point considéré.
Notons toutefois que cette proximité géométrique ne se traduit pas obligatoirement par une
proximité des inconnues discrètes, voir par exemple la matrice obtenue par la dsicrétisation du
Laplacien par le schéma de différences finies « à 5 points ». Avec la numérotation habituelle (par
ligne), si les voisins horizontaux sont numérotés consécutivement, les voisins verticaux sont à
une distance qui augmente avec le nombre de points de discrétisation. Ce qui reste vrai est que

42
le nombre d’inconnues par ligne de la matrice est borné indépendamment du nombre de points
de discrétisation.
Enfin, notons que c’est la combinaison de l’application et de la méthode de discrétisation qui
conduit à une matrice creuse. Les méthodes intégrales ou les méthodes spectrales conduisent à
des matrices pleines, mais de dimension plus faible.
Une autre application, avec laquelle l’auteur est moins familier, provient des matrices repré-
sentant des graphes, que l’on rencontre maintenant de plus en plus souventdans les applications
étudiant différentes formes de réseaux (internet, réseaux sociaux, ...).
Nous montrons sur la figure 4.1 deux exemples de matrices creuses. Elles proviennent toutes
deux de la collection de matrices creuses « SuiteSparse Matrix Collection 1 » [9]. L’Image de
gauche est la matrice Boeing/bcsstk36. Elle provient de la discrétisation par éléments finis
d’un absorbeur de chocs automobile. La matrice est de taille 23052, et a 1143140 éléments non
nuls, ce qui représente 0,2%, une économie d’un facteur 500.
L’image de droite est la matrice SNAP/soc-sign-bitcoin-alpha. Elle provient de la des-
cription du réseau de confiance sur une plateforme d’échange de bitcoinsSNAP/soc-sign-bitcoin-
alpha . Elle est de taille 3783, avec 24186 éléments non–nuls, soit seulement 0,17%. Elle est
proportionnellement plus creuse que la précédente, ce qui ne saute pas aux yeux !

Figure 4.1 – Deux exemples de matrice creuses

Ce chapitre sera une brève introduction au sujet. Il existe un grand nombre de modes de
stockage différents, chacun possédant des variantes. Nous ne chercherons pas à être exhaustif, et
ne présenterons que deux modes de stockage, avec les algorithmes associés :
le stockage par diagonales qui est bien adapté aux matrices provenant de discrétisation
d’EDP par différences finies sur des grilles régulières. Si ce mode est moins utilisé, il reste
toujours utile, et permet une introduction simple des algorithmes ;
le stockage dit « CSR » ou Compressed Sparse Rows, qui signifie stockage par lignes
creuses compressées. Ce mode permet de stocker des matrices creuses générales.
On trouvera une descrption plus détaillée des modes de stockage et des algorithmes associés
dans le livre de Y. Saad[20, Chap. 3]

1. https://sparse.tamu.edu/

43
4.2 Stockage par diagonales

Ce mode de stockage est adapté aux matrices dont les éléments sont disposés le long de
diagonales parallèles à la diagonale principale. Comme signalé plus haut, c’est en particulier le
cas des matrices obtenues par la discrétisation par différences finies d’une équation aux dérivées
partielles sur une grille régulière. Dans le cas le plus simple, en dimension 1, on trouve l’exemple
bien connu de la matrice tridiagonale (h est le pas de discrétisation).

 
2 −1
 
−1 2 −1
 

1  .. .. ..

(4.1) A= 2 . . .
 
h 


 

 −1 2 −1
−1 2

En dimension 2, la figure 3.1 montre une exemple d’une matrice de ce type.

4.2.1 Description

Pour décrire le stockage par diagonales, on considère une matrice A ∈ Rn×n , dont les élé-
ments sont disposés le long de Nd dagonales, dont les positions sont connues. Nous suivrons
la convention de Matlab qui consiste à repérer les diagonales par leur distance à la diagonale
principale. Par convention, celle–ci est la diagonale 0, les diagonales supérieures ont des dis-
tances positives et les diagonales inférieures des distances négatives. Par exemple, la matrice de
l’équation (4.1) a trois diagonales, en position -1, 0 et 1.
Le stocakge diagonal représente ce type de matrice par une structure consistant en
— un tableau à deux dimensions ADIAG, de taille n × Nd , dont chaque colonne correspond à
une diagonale ;
— un tableau IOFF (pour offsets, décalages) qui indique les numéros des diagonales.
Plus précisément, pour 1 ≤ i ≤ n et 1 ≤ k ≤ Nd , ADIAG(i, k). contient l’élément Ai,i+IOFF(k) .

Remarque 4.1. Il est commode d’utiliser un tableau rectangulaire ADIAG, mais on voit que seule
la diagonale principale est de taille n. Les autres diagonales sont incomplètes, soit « au début »
pour les sous–diagonales, soit « à la fin » pour les sur–diagonales. Les positions correspondantes
dans ADIAG, qui correspondraient soit à des indices négatifs, soit à des indices plus grands que
n ne seont pas utilisés (et on ne doit pas y accéder).

Exemple 4.1.
Nous prenons l’exemple de la matrice du Laplacien sur une grille3 × 3 (dont le seul intéret est
de pouvoir représenter toute la matrice !), qui donc une matrice 9 × 9. Au facteur 1/16 près, la

44
matrice est  
4 −1 −1
 
−1 4 −1 −1
 

 
−1 −1
 
 4 
 
−1 −1 −1
 
4 
 
A= −1 −1 4 −1 −1
 

 
−1 −1 −1
 
 4
 
−1 −1
 
 4 
 
−1 −1 4 −1
 

 
−1 −1 4

Les diagonales sont en positions Ioffs = (−3, −1, 0, 1, 3), et le tableau Adiag contient
 
∗ ∗ 4 −1
−1
 
 ∗ −1 4 −1 −1
 
 
 ∗ −1 4 0 −1
 
 
−1 0 4 −1 −1
 
 
Adiag = −1 −1 4 −1 −1
 
 
−1 −1 4 0 −1
 
 
−1 0 4 −1 ∗
 
 
−1 −1 4 −1 ∗
 
 
−1 −1 4 ∗ ∗

4.2.2 Algorithmes
L’algorithme le plus important pour exploiter le stockage creux est celui qui réalise le produit
de la matrice par un vecteur donné. Nous avons vu aux chapitres 1 et 2 que cette opération permet
d’implémenter les algorithmes du gradient conjugué et de GMRES. On note communément cette
opération matvec. Nous noterons x le vecteur à multplier, et y = Ax le résultat
Dans le cas du stockage par diagonales, il est naturel de boucler sur les diagonales. Comme
on l’a vu plus haut, ADIAG(i, k) contient l’élément Ai,i+d avec d = IOFF(k), et cet élément doit
être multiplié par xi+d pour s’ajouter au résultat yi .
Pour écrire l’algorithme, nous devons prendre en compte la remarque 4.1, qui dit que i et
i + d doivent être compris 1 et n.
On voit que, pour les matrices qui se prêtent à ce mode de stockage, l’algorithme 4.1 est
efficace. Le nombre d’opérations (additions et multiplications) est (presque) égal au nombre
d’éléments non–nuls : chaque passage dans la boucle « pour i » multiplie le vecteur x par une
diagonale (en ne tenant pas compte des éléments auxquels on n’accède pas). De plus, le sur–coût
du au mode de stockage est négligeable, puisqu’il suffit de calculer d pour chaque diagonal, et
que Nd est faible en pratique.
Les algorithmes pour résoudre des systèmes triangulaires dans le stockage par diagonales,
qui sont nécéssaires dans beaucoup de préconditionnement, par exemple pour utiliser une fac-

45
Algorithme 4.1 Produit matrice–vecteur pour le stockage par diagonales
y=0
pour k = 1 : Nd faire
d = IOFF(k)
i1 = max(1, 1 − d); i2 = min(n, n − d)
pour i = i1 : i2 faire
y(i) = y(i) + ADIAG(i, k)x(i + d)
fin pour
fin pour

torisation LU incomplète, sont plus difficiles à mettre en oeuvre. Il n’est pas facile, dans ce cas,
d’exploiter la structure par diagonales, et il est nécessaire de procéder par ligne ou par colonne.

4.3 Stockage CSR

Ce format de stockage est celui qui est le plus utilisé pour des matrices creuses générales, telles
que celles produites dans une méthode d’éĺéments finis. Il est utilisé, avec sa proche vaiante CSC,
qui propose un stockage par colonnes, dans un grand nombre de logiciels ou de bibliothèques
(c’est le stockage creux de Matlab, et Python utilise CSC).
Nous prendrons comme exemple la matrice suivante qui est naturellement trop petite pour
que le stockage creux soit utile, mais qui permet une description explicite des structures de
données.
 
A11 0 0 A14 0
 
A21 A22 0 A24 0 
 
 
(4.2) A = A31 0 A33 A34 A35 
 
 
 
 0 0 A43 A44 0 
 
0 0 0 0 A55

4.3.1 Stockage « coordonnées »

Avant de décrire le stockage CSR, nous discutons rapidement du stockage coordonnées, qui est
le plus naturel, mais ne permet pas d’implémenter facilement les algorithmes comme le produit
matrice vecteur.
Son principe est simple : l’élément Aij de la matrice est représenté par ses « coordonnées »
i et j, aisni que la valeur Aij .
Nous noterons NNZ le nombre d’éléments non–nuls de la matrice. Cette quantité importante
est, dans les cas que nous considérons, très inférieure au nombre d’éléments total de A, qui est
n2 . Pour le stockage coordonnées, nous aurons besoin de trois tableaux de taille NNZ, qui seront
notés IA (pour les indices de ligne), JA (pour les indices de colonne), et AA (pour les valeurs).
Pour 1 ≤ k ≤ NNZ, AA(k) contient la valeur de l’élément placé à la ligne i = IA(k) et la
colonne j = JA(k). La figure 4.2 illustre cela pour notre matrice exemple.
Comme nous le voyons sur la figure 4.2, les éléments peuvent être dans n’importe quel ordre.
Mais il est souvent commode de les trier, d’abord par ligne, puis par colonne croissante, comme

46
 AA A55 A35 A33 A24 A11 A14 A44 A21 A31 A22 A34 A43

IA 5 3 3 2 1 1 4 2 3 2 3 4

JA 5 5 3 4 1 4 4 1 1 2 4 3

Figure 4.2 – Stockage « ccordonnées » de la matrice (4.2), éléments non triés

illustré sur la figure 4.3.

AA A11 A14 A21 A22 A24 A31 A33 A34 A35 A43 A44 A55

IA 1 1 2 2 2 3 3 3 3 4 4 5

JA 1 4 1 2 4 1 3 4 5 3 4 5

Figure 4.3 – Stockage « ccordonnées » de la matrice (4.2), éléments triés

Comme nous l’avons déjà noté, ce format n’est pas très maniable. Sa simplicité fait qu’il
est souvent utilisé comme format d’entrée pour un logiciel avsnt d’être converti en interne en
un format plus maniable, comme le format CSR. C’est le cas de Matlab : la commande sparse
prend une description de la matrice au format coordonnées.

4.3.2 Le stockage CSR

Le stockage CSR se déduit par une variation simple du stockage coordonnée. Il conserve
les deux tableaux AA et JA mais remplace le tableau des indices de ligne IA par un tableau de
pointeurs vers le premier élément de chaque ligne. Comme nous le verrons plus loin, il permet de
mettre en oeuvre de manière efficace la plupart des algorithmes utiliśe par les méhodes itératives.
Pour l’exemple de la matrice A de refeq :matfull, les trois tableaux sont

AA A11 A14 A21 A22 A24 A31 A33 A34 A35 A43 A44 A55

JA 1 4 1 2 4 1 3 4 5 3 4 5

IA 1 3 6 10 12 13

Figure 4.4 – Stockage CSR de la matrice (4.2)

Comme nous l’avons vu plus haut, les tableaux AA et JA sont les mêmes que sur la figure 4.3.
Pour interpreter le tableau IA, notons que
— la première ligne va de l’élément stocké en IA(1) = 1 dans AA à l’élément stocké en
IA(2) − 1 = 2. La ligne 1 de A contient bien deux éléments, stockés dans les colonnes 1
et 4 ;
— la ligne 3 de A est entre les éléments IA(3) = 3 et IA(4) − 1 = 9 de IA, soit 4 éléments,
rangés aux places 6, 7, 8 et 9. Les colonnes correspondantes sont 1,3, 4 et 5.

47
 De manière générale : les éléments de la ligne i sont stockés dans les tableaux JA et AA aux
positions situées entre IA(i) et IA(i+1)−1. De plus, si j = JA(k), pour IA(i) ≤ k ≤ IA(i+1)−1,
alors AA(k) contient Aij .
Comme son nom l’indique, le stockage CSR est orienté par lignes. Cela aura bien entendu une
inluence sur la manière d’implémenter les différents algorithmes. Nous verrons qu’ils ont pour la
plupart une structure commune, formée de deux boucles imbriquées. La première parcourt les
lignes de la matrice, la seconde parcourt les éléments non–nuls ds la ligne en question.
Nous commencerons par l’algorithme qui est à la fois le plus simple et sans doute le plus
important : le produit d’une matrice creuse A stockée en format CSR par un vecteur x, et nous
noterons y = Ax le vecteur résultat. Comme l’algorithme naturel procède en prenant le produit
scalaire d’une ligne de A par x, la mise en oeuvre suit de près la description mathématique.

Algorithme 4.2 Produit matrice–vecteur pour le stockage CSR

y=0
     pour i = 1, . . . , n faire


 
    k1 = IA(i), k2 = IA(i + 1) − 1
   
    
   pour k = k1, . . . , k2 faire
=
  
j = JA(k); y(i) = y(i) + AA(k)x(j)
   
    
     
fin pour
fin pour

Nous pouvons également écrire cet algorithme à un plus haut niveau, en remplaçant la boucle
interne par une somme. Ces deux versions sont mathématiquement équivalentes, et il n’y a pas
de raison objective de préféreer l’une à l’autre. Le niveau de détail supérieur dans la première
version peut être vu comme un avantage ou un inconvénient.

Algorithme 4.3 Produit matrice vecteur en stockage CSR — version alternative

pour i = 1, . . . , n faire
k1 = IA(i), k2 = IA(i + 1) − 1
k2
X
y(i) = AA(k) x(JA(k))
k=k1
fin pour

Dans certains cas, il peut être nécessaire de calculer le produit AT x de la transposée de A

par un vecteur, naturellement sans stocker explicitement AT . Puisque A est stocké par lignes, il
faut trouver un algorithme qui accede aux colonnes de AT . Pour cela, nous devons interprérter
le produit d’une matrice par un vecteur d’une manièere différente.
Soit dinc B ∈ Rn une matrice, x ∈ Rn un vecteur. Plus tard, B sera AT , mais il est plus
clair de généraliser la situation. En revenant à la définition du produit Bx comme l’application
(dans une base fixée) de l’application linéaire représentée par la matrice B au vecteur x, nous
voyons que le produit y = Bx est une combinaison linéaire des colonnes de B par les coefficients
de x. En notant B = [b1 , . . . , bn ] les colonnes de B, nous avons donc :
n
X
y= xj bj ,
j=1

48
ce qui conduit à un algorithme par colonnes pour le produit
pour j = 1, . . . , n faire
pour i = 1, . . . , n faire
yi = yi + Bij xj
fin pour
fin pour
Dans le cas où la matrice B est AT , on obtient, en échageant le rôle des indices i et j :
pour i = 1, . . . , n faire
pour j = 1, . . . , n faire
yj = yj + Aij xi
fin pour
fin pour
et cette fois, on accède aux éléments de A par ligne.
En adaptant cet algorithme au stockage CSR, nous obtenons finalement

Algorithme 4.4 Produit de la transpoée d’une matrice par un vecteur en stockage CSR

y=0
     pour i = 1, . . . , n faire
          k1 = IA(i), k2 = IA(i + 1) − 1




   
   pour k = k1, . . . , k2 faire
  =  j = JA(k); y(j) = y(j) + AA(k)x(i)
    

         
  
 fin pour
fin pour

49
Annexe A

Rappels d’algèbre linéaire

Ce chapitre contient des résultats utiles d’algèbre linéaire, rappelés sans démonstration. Pour
plus de détails, on consultera les références suivantes [1, 8, 10, 15].

A.1 Matrices symétriques

Définition A.1. Une matrice A ∈ Rn×n est symétrique si et seulement si :
A = AT , c’est à dire : aij = aji , ∀(i, j) ∈ [1, n]
ou de manière équivalente :
(Ax, y) = (x, Ay), ∀(x, y) ∈ Rn ,
où le produit scalaire est défini par
(x, y) = y T x, ∀(x, y) ∈ Rn .
Rappelons que toutes les valeurs propres d’une matrice symétrique sont réelles, que les vec-
teurs propres correspondant à des valeurs propres distinctes sont orthogonaux, et qu’une matrice
symétrique est diagonalisable dans une base orthonormée.
Définition A.2. Une matrice symétrique A ∈ Rn×n est définie positive si et seulement si elle
vérifie l’une des conditions (équivalentes) suivantes :
i) (Ax, x) > 0, ∀x ∈ Rn , x 6= 0 ;
ii) toutes les valeurs propres de A sont strictement positives.
iii) Il existe une matrice L ∈ Rn inversible et triangulaire inférieure telle que A = LLT
(factorisation de Cholesky).

A.2 Normes vectoriellles et matricielles

Nous rappelons dans ce paragraphe les principales propriétés des normes matricielles, essen-
tiellement sans définition (voir les références générales citées précédemment).
Les normes permettent de « mesurer » la taille des éléments d’un espace vectoriel. Nous
verrons au paragraphe A.3 qu’elles ont un rôle essentiel à jouer pour quantifier la sensibilité
de la solution d’un système linéaire à des perturbations sur les données (second membre ou
matrice).

50
Définition A.3. Une norme sur un espace vectoriel E est une application de E dans R+ , notée
k.k telle que :
— ∀(x, y) ∈ E 2 , kx + yk ≤ kxk + kyk ;
— ∀x ∈ E, ∀λ ≥ 0, kλxk = |λ| kxk ;
— kxk = 0 ⇒ x = 0.

Exemple A.1.
Quand l’espace vectoriel E est de dimension finie (E = Rn ), les exemples les plus courants sont
P 1/2
n 2
— La norme euclidienne : kxk2 = |x
i=1 i | ;
— La norme du max : kxk∞ = maxi=1,...,n |xi | ;
— Plus généralement, la norme-p est définie pour p ∈ [1, +∞[ par kxkp = ( ni=1 |xi |p )1/p .
P
La norme 2 en est un cas particulier.

Définition A.4. Étant donné une norme vectorielle, la quantité

kAxk
(A.1) kAk = sup = sup kAxk
x6=0 kxk kxk=1

définit une norme sur l’espace vectoriel des matrices, dite norme subordonnée à la norme vecto-
rielle.

Une norme vectorielle subordonnée vérifie l’inégalité

(A.2) kABk ≤ kAk kBk ,

pour deux matrices pour lesquelles le produit est défini, et en particulier

kAxk ≤ kAk kxk .

Dans le cas des normes p pour p = 1, p = 2 et p = ∞, on dispose d’expressions explicites

pour les normes matricielles correspondantes :
Pour p = 1, kAk1 = max1≤j≤n m
P
i=1 |aij | ;
√
Pour p = 2, kAk2 = max1≤i≤n λi , où λi est une valeur propre de AT A ;
Pour p = ∞, kAk∞ = max1≤i≤m nj=1 |aij |.
P

Terminons en signalant qu’il existe des normes sur l’espace des matrices qui ne sont pas des
normes subordonnées. L’exemple le plus courant est la norme de Frobenius :
 1/2
Xm X
n
kAkF =  |aij |2  .
i=1 j=1

Comme une norme subordonnée vérifie kIk = 1 (pourquoi ?), et que kIkF = n, la norme de
Frobenius ne peut-être une norme subordonnée. Par contre, elle vérifie l’inégalité (A.2).
Le Lemme suivant sera utilisé au paragrpahe A.3 pour déterminer la taille d’une perturbation
pouvant laisser une matrice inversible.

51
Lemme A.1. Soit k.k une norme subordonnée, et A une matrice vérifiant kAk < 1. La matrice
I + A est inversible. De plus,
1
(A.3) (I + A)−1 ≤ .
1 − kAk

Preuve. Raisonnons par l’absurde : si I + A est singulière, soit x un vecteur non nul dans le noyau de
cette matrice. On a alors
(I + A)x = 0 ⇒ Ax = −x
soit en passant aux normes
kxk = kAxk ≤ kAk kxk < kxk,
ce qui est évidemment une contradiction.
Pour prouver (A.3), notons B = (I + A)−1 . Partant de B(I − A) = I, il vient B = I + BA, puis en
prenant les normes kBk ≤ 1 + kBk kAk, ce qui donne (A.3). 

A.3 Conditionnement d’un système linéaire

Le conditionnement d’un système linéaire (ou plus précisément d’une matrice) est un nombre
qui mesure comment la solution de ce système varie si les données (le second membre, ou les
coefficients de la matrice) sont perturbés. Cette mesure est une propriété intrinsèque du système,
et ne dépend pas de la méthode de résolution.
Il sera commode de séparer l’effet des perturbations sur le secopnd membre et sur la matrice.

Définition A.5. Étant donné une matrice inversible A ∈ Rn×n , et une norme matricielle
subordonnée k.k, le conditionnement de A, noté κ(A) est le nombre défini par (pour la norme
considérée) :

(A.4) κ(A) = kAk A−1

Théorème A.1. Soit A ∈ Rn×n une matrice inversible, et b ∈ Rn supposé non-nul. Nous
considérons le système linéaire Ax = b, et des perturbations δb ∈ Rn et δA ∈ Rn×n .
— Pour le système perturbé Ax̂ = b̂, avec b̂ = b + δb, et x̂ = x + δx, l’erreur relative sur x
vérifie :

kδxk kδbk
(A.5) ≤ κ(A) .
kxk kbk

— Pour le système perturbé Âx̂ = b, avec  = A + δA et x̂ = x + δx, où δA vérifie

A−1 δA < 1, l’erreur relative sur x vérifie

kδxk κ(A) kδAk

(A.6) ≤ .
kxk kδAk kAk
1 − κ(A)
kAk

Dans le cas où la norme est une norme p, nous noterons κp le conditionnement.

Proposition A.1. Propriétés du conditionnement. Soit A ∈ Rn une matrice inversible.

— Pour toute norme subordonnée κ(A) ≥ 1 ;
— κ(λA) = κ(A), ∀λ ∈ R∗ ;

52
 — Pour la norme euclidienne, κ2 (A) = 1 si A est (proportionnelle à) une matrice orthogo-
nale.
λn
— Si A est une matrice symétrique et définie positive, κ2 (A) = où λn (resp. λ1 ) est la
λ1
plus grande (resp. petite) valeur propre de A.

On dit qu’une matrice est mal conditionnée si son conditionnement est « très grand ». Dans
ce cas, la solution d’un système linéaire est très sensible à des perturbations sur ses données.

A.4 Problèmes de moindres carrés

Nous rappelons dans ce paragraphe les principaux résultats concernant les problèmes de
moindres carrés linéaires. Ces problèmes interviennent dans des domaines très divers chaque
fois que l’on cherche à estimer des paramètres à partir de données expérimentales, dans le cas
où le modèle est linéaire. L’utilisation d’un nombre de mesures plus grand que le nombre de
paramètres permet de réduire l’influence des erreurs expérimentales.
Nous considérons donc une matrice A ∈ Rm×n , avec m ≥ n (et en général m 6= n), et un
vecteur b̂ ∈ Rm . Nous cherchons un vecteur x ∈ Rn solution du système linéaire

Ax = b.

Il est connu que, dans le cas ou m > n, ce système n’a en général pas de solution, et que si
une solution existe (si b ∈ Im A), celle-ci ne sera en général pas unique, mais définie à l’addition
d’un élément de Ker A près. On va alors chercher un vecteur x tel que l’équation soit vérifiée
« au mieux », en un sens qu’il faut naturellement préciser. Le choix le plus commun, dont l’une
des justifications est son lien avec le minimum de vraisemblance est de minimiser la norme
euclidienne du résidu. On résout donc le problème

(A.7) min kAx − bk22

x∈Rn

Comme nous ne ferons qu’effleurer le sujet, donnons quelques références qui permettent de
l’approfondir . En français, citons les livre [15], ou [1]. En anglais, le classique est [10], une
référence encyclopédique sur les problèmes de moindres carrés est [6]. Un autre livre, moins
complets, mais plus abordable, est [22].

A.4.1 Propriétés mathématiques

Nous allons voir que le problème (A.7) est équivalent à une équation linéaire, mais pour un
opérateur différent de A.

Théorème A.2. Soit A ∈ Rm×n et soit b ∈ Rm . Un élément x̂ ∈ Rn est une solution de (A.7)
si et seulement si

(A.8) AT Ax̂ = AT b

Démonstration. On peut obtenir facilement (A.8) en calculant le gradient de la fonctionnelle

1
x → kAx − zk2 . Nous laissons cette méthode en exercice au lecteur. Nous présenterons plutôt
2
l’argument élémentaire suivant, emprunté à Björck [6].

53
 Soit x vérifiant (A.8). On a pour tout y ∈ Rn :

b − Ay = b − Ax + A(x − y).

L’équation normale (A.8) implique que les deux termes de la somme sont orthogonaux (le résidu
b − Ax est orthogonal a l’image de A). Le théorème de Pythagore implique :

kb − Ayk22 = kb − Axk22 + kA(x − y)k22 ≥ kb − Axk22 .

x est donc bien solution de (A.7)

Réciproquement, soit x tel que AT (b − Ax) = w 6= 0. Choisissons y = x + εw, avec ε > 0.
On a alors :

kb − Ayk22 = (b − Ay, b − Ay) = kb − Axk22 − 2ε(b − Ax, w) + kAwk22 < kb − Axk22

si ε est suffisamment petit. x n’est donc pas solution de (A.7). 

Remarque A.1. L’équation normale (A.8) se réécrit :

AT (Ax̂ − b) = 0,

ce qui exprime simplement que le résidu b − Ax̂ est dans le noyau de AT , c’est-à-dire orthogonal
à l’image de A (voir le théorème A.3). Ceci conduit à l’illustration géométrique bien connue :

b − Ax

Ax

Im A

Figure A.1 – Illustration géométrique des moindres carrés

La solution du problème de moindres carrés est telle que Ax est la projection de b sur l’image
de A.

Notons que nous n’avons pour l’instant évoqué ni l’existence, ni l’unicité, pour les solutions
de (A.7) (ou de (A.8)). Nous avons simplement montré l’équivalence des deux problèmes. Nous
allons préciser ces points dans le théorème suivant.

Théorème A.3. Le problème (A.7) admet toujours une solution. De plus, cette solution est
unique si, et seulement si, la matrice A est injective (ou ce qui est équivalent de rang maximal)

54
Démonstration. Rappelons les résultats généraux d’algèbre linéaire :

Ker AT A = Ker A
(A.9)
Im AT A = Im AT

Un sens est évident à chaque fois. Pour l’autre sens, nous avons pour la première égalité :

AT Ax = 0 ⇒ (AT Ax, x) = 0 ⇒ (Ax, Ax)Rm = kAxk2 = 0 ⇒ Ax = 0.

La deuxième s’en déduit sans autre calcul car :

(Im AT ) = (Ker A)⊥ = (Ker AT A)⊥ = Im AT A.

Nous obtenons donc l’existence d’une solution, puisque le second membre de (A.8) est élément
de Im AT = Im AT A.
Enfin, A et AT A sont injectifs en même temps, ce qui donne le résultat, puisque AT A est
carrée donc inversible si elle est injective. 

Sous l’hypothèse que A est de rang maximal, nous pouvons donner un résultat plus précis :

Proposition A.2. Sous l’hypothèse que A est de rang n, la matrice des équations normales
AT A est définie positive.

Démonstration. Il est facile de voir que AT A est semi-définie positive (c’est en fait vrai indé-
pendamment du rang de A) :

(AT Ax, x) = (Ax, Ax) = kxk22 ≥ 0.

De plus, quand A est de rang n, nous avons vu (au lemme A.3) que la matrice (carrée) AT A est
injective, donc inversible. Elle est donc définie positive. 

A.4.2 Méthodes numériques

Équations normales et factorisation de Cholesky
Comme nous l’avons vu au paragraphe A.4.1 (équation (A.8)), la solution d’un problème de
moindres carrés se ramène, en théorie du moins, à la résolution d’un système linéaire pour la
matrice AT A (dite matrice des équation normales) :

(A.10) AT Ax = AT z.

Nous avons vu également (proposition A.2) que, sous l’hypothèse que A est rang n, la ma-
trice AT A est définie positive. Cette matrice est de taille n par n, et les équations normales
représentent une « compression » d’information, puisque n ≤ m . Le système (A.10) peut donc
(toujours en théorie) être résolu par la factorisation de Cholesky. Nous rappelons cette méthode
bien connue :

Proposition A.3. Soit C une matrice symétrique et définie positive. Il existe une unique matrice
R triangulaire supérieure, à éléments diagonaux strictement positifs, telle que

(A.11) C = RT R

55
Preuve. Elle se trouve dans les références citées ci-dessus. Notons que la preuve est constructive : elle
fournit un algorithme pour calculer la matrice R à partir de C. 

Il suffit ensuite de résoudre les deux systèmes linéaires


 R T y = AT z
(A.12)
 Rx = y

pour obtenir la solution. Le premier a une matrice triangulaire inférieure, le second une matrice
triangulaire supérieure.
Le nombre d’opérations flottantes nécessaires à la formation de la matrice des équations
normales est n(n + 1)m + mn (en tirant parti de la symétrie). Le nombre d’opérations de
la factorisation de Choleski est n3 /3 opérations flottantes (additions et multiplications), plus
n divisions et n extractions de racines carrées. La solution des équations triangulaires (A.12)
prend n2 opérations, et est donc négligeable.
Le coût de la méthode des équations normales est donc, en ne gardant que les termes prin-
cipaux : mn2 + 31 n3 .
En dépit de sa simplicité, et de son coût raisonnable (nous verrons que les autres méthodes
de résolution des problèmes de moindres carrés sont plus chères), la méthodes des équations
normales n’est pas recommandée, pour des raisons de stabilité. Elle cumule deux inconvénients :
i) Le simple fait de former la matrice AT A peut faire perdre de l’information sur les petits
coefficients de la matrice A, ce qui peut avoir des conséquences désastreuses, comme le
montre l’exemple A.2 ci-dessous.
ii) De plus, le conditionnement de la matrice AT A est le carré de celui de A (puisque les
valeurs propres de AT A sont les carrés des valeurs singulières de A).

Exemple A.2 (d’après Björck [6]).

Soit le problème minx∈R3 kAx − zk, avec
 
1 1 1
 
 ε 0 0
   T
A=  ,
 z = 1, 0, 0, 0 .
0 ε 0
 
0 0 ε

Cet exemple peut correspondre à un problème où la somme x1 +x2 +x3 est connue avec beaucoup
plus de précision de que les composants de x individuellement. La calcul exact donne
 
1 + ε2 1 1
   T 1  T
T
A= 1 1 + ε2 1  , A z = 1, 1, 1 , x =
 
1, 1, 1
  3 + ε2
1 1 1+ε 2

Supposons que ε = 10−4 et que le calcul soit effectué avec 8 chiffres significatifs. Alors 1 + ε2 =
1.000000001 sera arrondi à 1, et la matrice AT A calculée sera singulière. Tout se passe comme
si les trois dernières lignes de A étaient nulles, et n’avaient contribué aucune information.

56
 Les difficultés illustrées par cet exemple peuvent être évitées, ou tout du moins les limites re-
poussées, en passant en double précision. Dans ces conditions, les équations normales deviennent
une méthode acceptable. Toutefois, la méthode présentée au paragraphe suivant lui est préférable
comme méthode de résolution générale.

La factorisation QR
La méthode « moderne » pour résoudre les problèmes de moindres carrés est basée sur une
factorisation dite QR de la matrice A, où Q est une matrice orthogonale, et R est triangulaire
supérieure. La réduction de A à une forme triangulaire s’effectue par une suite de multiplica-
tions par des matrices orthogonales élémentaires, qui peuvent être soit des symétries (méthode
de Householder), soit des rotations (méthode Givens). Nous ne détaillerons pas cette étape, ren-
voyant encore une fois aux références citées plus haut, et nous nous contenterons de donner le
résultat de la factorisation, et de montrer comment elle est utilisée pour résoudre le problème
de moindres carrés.

Théorème A.4. Soit A ∈ Rm×n une matrice de rang n. Il existe une matrice orthogonale
Q ∈ Rm×m et une matrice triangulaire supérieure inversible telles que
 
R
(A.13) A = Q 
0

On peut montrer (voir [22]) que le nombre d’opérations flottantes est (au premier ordre)
2
2mn2 − n3 . Le calcul de la factorisation QR est donc plus cher que la formation, et la fac-
3
torisation, de la matrice des équations normales. Le sur-coût n’est toutefois que d’un facteur 2
(l’ordre de complexité est le même). De plus, (voir les références citées au début de ce chapitre),
les propriétés numériques de la méthode QR sont supérieures (dues au fait que l’on travaille
avec des matrices orthogonales), et en font la méthode de choix pour une implémentation ro-
buste. En pratique, les bibliothèques numériques fournissent toutes une version de la résolution
de problèmes de moindres carrés basées sur la factorisation QR. Dans les programmes interactifs
comme Scilab où Matlab, la commande x= A\z résolvent le problème de moindres carrés quand
A est rectangulaire, et sont basées sur la factorisation QR de A.

La décomposition QR possède une importante interprétation en terme d’orthogonalisa-

tion. Pour préciser cela, partitionnons la matrice Q en Q = (Q1 Q2 ), avec Q1 ∈ Rm×n ,
Q1 ∈ R(m−n)×n . Les colonnes de Q1 forment une base orthogonale de l’image de A, celle de
Q2 une base orthonormale de Im A⊥ , et l’on a la factorisation réduite

(A.14) A = Q1 R.

Cette factorisation montre que, pour k ≤ n, les k premières colonnes de Q forment une base
orthonormale du sous-espace engendré par les k premières colonnes de A. De plus, comme R est
inversible, ces sous–espace sont égaux pour chaque k.
La factorisation QR construit donc une base orthonormale pour tous les sous espaces engen-
drés par les colonnes de A. On trouve le résultat obtenu habituellement par l’orthogonalisation
de Gram–Schmidt. Il y a toutefois deux différences notables :

57
 — l’algorithme de Gram–Schmidt est notoirement instable : les vecteurs calculés par cette
méthode ne sont pas numériquement orthogonaux. L’algorithme QR calcule le même
résultat de façon numériquement stable (voir ci-dessous). Il est possible d’utiliser une
variante (l’algorithme de Gram–Schmidt modifié), qui possède de meilleures propriétés
de stabilité, mais celles de l’algorithme QR sont encore meilleures.
— L’algorithme de Gram–Schmidt ne calcule que la matrice que nous avons notée Q1 .
Comme nous l’avons vu, l’algorithme QR fournit en plus la matrice Q2 . Selon les ap-
plications, cette différence peut ou on être importante.
Remarque A.2. La relation (A.14) montre que R n’est autre que le facteur de Cholesky de AT A.
En effet, AT A = RT QT1 Q1 R = RT R, puisque Q1 est orthogonale. Cette remarque montre alors
que R est déterminée de façon unique par A (si A est de rang n), et donc Q1 = AR−1 l’est
également. Par contre, Q2 n’est déterminée que par la condition que (Q1 Q2 ) soit orthogonale,
et n’est donc pas unique.
Une fois la factorisation QR obtenue, il est facile de résoudre le problème de moindres
carrés (A.7).
 
R
Théorème A.5. Soir A = Q   la décomposition QR de la matrice A. La solution du
0
problème de moindres carrés est donnée par la résolution du système linéaire

(A.15) Rx = z1 ,
 
z1
avec QT z =  .
z2

Preuve. Du fait que la matrice Q est orthogonale, nous avons

  2   2
2 R R
kAx − zk2 = Q   x − z =   x − QT z
0 0
2 2
 
z1
Posons QT z =   comme indiqué dans l’énoncé du théorème, il vient
z2

2 2 2
kAx − zk2 = kRx − z1 k2 + kz2 k2 .

Le deuxième terme est constant (indépendant de x), et le premier est minimisé en prenant x = R−1 z1 ,
2
ce qui est possible puisque R est inversible. Pour conclure, notons que le terme kze k2 représente la valeur
du résidu à la solution. 

A.5 Méthodes itératives classiques

58
Bibliographie

[1] Grégoire Allaire and Mahmoud Sidi Kaber. Algèbre linéaire numérique : Cours et exercices.
Ellipses, 2002.
[2] O. Axelsson and V. A. Barker. Finite element solution of boundary value problems, vo-
lume 35 of Classics in Applied Mathematics. Society for Industrial and Applied Mathe-
matics (SIAM), Philadelphia, PA, 2001. Theory and computation, Reprint of the 1984
original.
[3] Owe Axelsson. Iterative solution methods. Cambridge University Press, Cambridge, 1994.
[4] Richard Barrett, Michael Berry, Tony F. Chan, and et al. Templates for the solution of
linear systems : building blocks for iterative methods. Society for Industrial and Applied
Mathematics (SIAM), Philadelphia, PA, 1994.
[5] Michele Benzi. Preconditioning techniques for large linear systems : a survey. J. Comput.
Phys., 182(2) :418–477, 2002.
[6] Åke Björck. Numerical methods for least squares problems. Society for Industrial and
Applied Mathematics (SIAM), Philadelphia, PA, 1996.
[7] Ke Chen. Matrix preconditioning techniques and applications, volume 19 of Cambridge
Monographs on Applied and Computational Mathematics. Cambridge University Press,
Cambridge, 2005.
[8] Philippe G. Ciarlet. Introduction à l’analyse numérique matricielle et à l’optimisation.
Collection Mathématiques Appliquées pour la Maîtrise. Masson, Paris, 1982.
[9] T. Davis and Y. Hu. The University of Florida Sparse Matrix Collection. ACM Transactions
on Mathematical Software, 38(1) :1–25, 2011.
[10] Gene H. Golub and Charles F. Van Loan. Matrix computations. Johns Hopkins Studies in
the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, 4ème edition,
2013.
[11] Xavier Gourdon. Algèbre. Ellipses, 2009.
[12] Anne Greenbaum. Iterative methods for solving linear systems, volume 17 of Frontiers in
Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadel-
phia, PA, 1997.
[13] Anne Greenbaum, Vlastimil Pták, and Zdeněk Strakoš. Any nonincreasing convergence
curve is possible for GMRES. SIAM J. Matrix Anal. Appl., 17(3) :465–469, 1996.
[14] C. T. Kelley. Iterative methods for linear and nonlinear equations, volume 16 of Fron-
tiers in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM),
Philadelphia, PA, 1995.

59
[15] P. Lascaux and R. Théodor. Analyse numérique matricielle appliquée à l’art de l’ingénieur.
Tome 1. Masson, Paris, 2 edition, 1993. Méthodes directes.
[16] P. Lascaux and R. Théodor. Analyse numérique matricielle appliquée à l’art de l’ingénieur.
Tome 2. Masson, Paris, 2 edition, 1994. Méthodes itératives.
[17] Roger Mansuy and Rached Mneimné. Algèbre linéaire. Rédution des endomorphismes. Vui-
bert, 2012.
[18] Gérard Meurant. Computer solution of large linear systems, volume 28 of Studies in ma-
thematics and its applicaitons. Elsevier North-Holland, 2 edition, 1999.
[19] Youcef Saad and Martin H. Schultz. GMRES : a generalized minimal residual algorithm for
solving nonsymmetric linear systems. SIAM J. Sci. Statist. Comput., 7(3) :856–869, 1986.
[20] Yousef Saad. Iterative methods for sparse linear systems. Society for Industrial and Applied
Mathematics, Philadelphia, PA, 2 edition, 2003.
[21] Denis Serre. Matrices, volume 216 of Graduate Texts in Mathematics. Springer, New York,
second edition, 2010. Theory and applications.
[22] Lloyd N. Trefethen and David Bau, III. Numerical linear algebra. Society for Industrial
and Applied Mathematics (SIAM), Philadelphia, PA, 1997.
[23] Henk A. van der Vorst. Iterative Krylov methods for large linear systems, volume 13 of
Cambridge Monographs on Applied and Computational Mathematics. Cambridge University
Press, Cambridge, 2009. Reprint of the 2003 original.
[24] A. J. Wathen. Preconditioning. Acta Numer., 24 :329–376, 2015.

60