Département de mathématiques

Master MASI

Résolution itérative de systèmes linéaires de

grande taille par des méthodes de Krylov
Pr. A. ARCHID 25/09/2018

1
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
1 / 82
Plan

1 Rappels

2 Méthodes de Krylov
Introduction
Méthode d’Arnoldi
Méthode de Lanczos
Lanczos Bi-orthogonal
Méthode de gradient conjugué (GC)
Lien entre Gradient conjugué et Lanczos
Méthode de GMRES

2
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
2 / 82
Plan

1 Rappels

2 Méthodes de Krylov
Introduction
Méthode d’Arnoldi
Méthode de Lanczos
Lanczos Bi-orthogonal
Méthode de gradient conjugué (GC)
Lien entre Gradient conjugué et Lanczos
Méthode de GMRES

2
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
2 / 82
Plan

1 Rappels

2 Méthodes de Krylov
Introduction
Méthode d’Arnoldi
Méthode de Lanczos
Lanczos Bi-orthogonal
Méthode de gradient conjugué (GC)
Lien entre Gradient conjugué et Lanczos
Méthode de GMRES

2
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
2 / 82
Plan

1 Rappels

2 Méthodes de Krylov
Introduction
Méthode d’Arnoldi
Méthode de Lanczos
Lanczos Bi-orthogonal
Méthode de gradient conjugué (GC)
Lien entre Gradient conjugué et Lanczos
Méthode de GMRES

2
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
2 / 82
Plan

1 Rappels

2 Méthodes de Krylov
Introduction
Méthode d’Arnoldi
Méthode de Lanczos
Lanczos Bi-orthogonal
Méthode de gradient conjugué (GC)
Lien entre Gradient conjugué et Lanczos
Méthode de GMRES

2
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
2 / 82
Plan

1 Rappels

2 Méthodes de Krylov
Introduction
Méthode d’Arnoldi
Méthode de Lanczos
Lanczos Bi-orthogonal
Méthode de gradient conjugué (GC)
Lien entre Gradient conjugué et Lanczos
Méthode de GMRES

2
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
2 / 82
Plan

1 Rappels

2 Méthodes de Krylov
Introduction
Méthode d’Arnoldi
Méthode de Lanczos
Lanczos Bi-orthogonal
Méthode de gradient conjugué (GC)
Lien entre Gradient conjugué et Lanczos
Méthode de GMRES

2
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
2 / 82
Plan

1 Rappels

2 Méthodes de Krylov
Introduction
Méthode d’Arnoldi
Méthode de Lanczos
Lanczos Bi-orthogonal
Méthode de gradient conjugué (GC)
Lien entre Gradient conjugué et Lanczos
Méthode de GMRES

2
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
2 / 82
Plan

1 Rappels

2 Méthodes de Krylov
Introduction
Méthode d’Arnoldi
Méthode de Lanczos
Lanczos Bi-orthogonal
Méthode de gradient conjugué (GC)
Lien entre Gradient conjugué et Lanczos
Méthode de GMRES

2
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
2 / 82
 Rappels

Orthonormalisation

Étant donnée une famille de vecteurs {a1 , · · · , an } d’un espace vectoriel de Rn , on

cherche à construire une base orthonormée {v1 , · · · , vn } de Vect{a1 , · · · , an }
comme suit :
Procédé d’orthonormalisation de Gram-Schmidt :
u1 = a1 ; v1 = kuu11 k
u2 = a2 − ha2 , v1 iv1 ; v2 = kuu22 k
.. ..
. .
n−1
P un
un = an − han , vi ivi = an − han , v1 iv1 − ..... − han , vn−1 ivn−1 = ; vn = kun k
i=1

3
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
3 / 82
 Rappels

Matrice creuse

Matrice creuse
X Une matrice creuse est une matrice contenant beaucoup de zéros.
X Une matrice creuse de taille importante est une matrice d’ordre pouvant aller
jusqu’à 108 mais dont la plupart des coefficients sont nuls.
X Certaines matrices creuses de grande dimension ne sont pas manipulables par
les algorithmes classiques.
X Les matrices creuses de très grande taille apparaissent souvent en science ou
en ingénierie (mécanique des fluides, traitement d’image satellite...) pour la
résolution des équations aux dérivées partielles.

4
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
4 / 82
 Rappels

Matrice de Hessenberg

Définition
Une matrice carrée Hm , d’ordre m, est dite de Hessenberg supérieure (resp.
inférieure) si tous ses coefficients situés au dessous de sa première sous-diagonale
sont nuls c-à-d ∀i > j + 1 hij = 0 (resp. ∀j > i + 1 hij = 0)
 
h11 h12 . . . h1m−1 h1m
 h21 h22 . . . h2m−1 h2m 
.. .. 
 
 ..
Hm = (hij ) = 
 0 h 32 . . .  ∈ Mm (K)

 .. . .. . .. .
.. .
.. 
 . 
0 . . . 0 hmm−1 hmm

5
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
5 / 82
 Rappels

Méthode directes - Méthode itératives

Définition
Il existe deux grandes classes de méthodes de résolution de systèmes linéaires :
1 Les méthodes directes (Gauss, LU, ...) qui aboutissent à la solution au bout
d’un nombre fini d’opérations.
2 Les méthodes itératives classiques (Jacobi, Gauss-Siedel, SOR, ...) qui
génèrent, à partir d’une solution de départ x0 , une suite d’estimations
(xk )k∈N qui tend vers la solution du problème, où l’itéré xk se calcul en
fonction des itérés précédents x0 , x1 , · · · , xk−1 .

X Elles nécessitent en général moins d’espace mémoire que les méthodes

directes. (elles ne font appel qu’à des produits matrice-vecteur)
XLeur convergence n’est pas toujours assurée et est en général lente, surtout
pour les systèmes de grande taille.
XElles servent à accélérer la convergence d’une classe importante de méthodes,
appelées méthodes polynomiales.
6
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
6 / 82
 Rappels

Méthode directes - Méthode itératives

Un des principaux inconvénients de ces méthodes, est qu’elles ne sont pas

utilisables lorsque la matrice du système est creuse de grande taille . (Par
exemple,les méthodes basées sur des factorisations matricielles sont trop coûteuses
en temps de calcul et en mémoire).

7
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
7 / 82
 Méthodes de Krylov Introduction

Méthodes de Krylov

Introduction
La résolution de nombreux problèmes que nous rencontrons en physique, en
mécanique, en chimie etc., et d’une façon générale dans les sciences de l’ingénieur,
se ramène à la résolution de grands systèmes linéaires creux Ax = b,(( où A est
une matrice creuse inversible à coefficients réels)), pour lesquelles il est
fondamental d’introduire des méthodes numériques plus fiables et faciles à utiliser,
parmi ces méthodes nous citons les méthodes de type Krylov.

Alexeï Nikolaïevitch Krylov était un mathématicien et ingénieur russe qui a vécu

de 1863 à 1945. Les espaces qu’il a utilisés pour des calculs de valeurs propres
sont associés à son nom, mais les méthodes « dites de Krylov », en particulier
pour les systèmes non symétriques, ont été découvertes beaucoup plus tard.

8
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
8 / 82
 Méthodes de Krylov Introduction

Définition : Espace de Krylov

Soit x0 ∈ Rn un vecteur initial. Le vecteur résidu associé à x0 est défini par
r0 = b − Ax0 . On appelle espace de Krylov d’ordre m associé à la matrice carrée
inversible A ∈ Mn (R) et r0 noté Km (A, r0 ), l’espace vectoriel généré par r0 et ses
m − 1 produits itérés par A.

Km (A, r0 ) = Vect r0 , Ar0 , ...., Am−1 r0 = {p(A)r0 /p ∈ Pm−1 }



où Pm−1 est l’ensemble des polynômes de degré au plus m − 1.

Prpriété
La suite des sous-espaces de Krylov (Km )m>0 est une suite croissante au sens
de l’inclusion Km ⊂ Km+1 .

9
Résolution itérative de systèmes linéaires de grande taille par
11des
décembre
méthodes
2018
de Krylov
9 / 82
 Méthodes de Krylov Introduction

Théorème
Si Am r0 ∈ Km (A, r0 ), alors Am+p r0 ∈ Km (A, r0 ) pour tout p > 0.

Preuve : La démonstration se fait par récurrence sur p

m−1
Am+p r0 ∈ Km (A, r0 ) alors Am+p r0 = αk Ak r0
P
Si pour p ≥ 0,
k=0
m−1
X
D’où Am+p+1 r0 = αk Ak+1 r0
k=0
m−2
X
= αk Ak+1 r0 + αm−1 Am r0
k=0
m−2 m−1
!
X X
= αk Ak+1 r0 + αm−1 βk Ak r0
k=0 k=0
m−1
X m−1
X
= αk−1 Ak r0 + αm−1 βk Ak r0
k=1 k=0

puisque ces deux termes appartiennent à Km (A, r0 ), on peut donc noter

m−1
Am+p+1 r0 = λk Ak r0 ∈ Km (A, r0 )
P
10
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
10 / 82
 Méthodes de Krylov Introduction

Théorème
La suite d’espaces de Krylov Km (A, r0 ) est strictement croissante de 1 à mmax ,
puis stagne à partir du rang mmax .

C’est-à-dire : Kq = Kmmax , ∀q ≥ mmax

Preuve : Si m est le plus petit entier pour lequel Am r0 est dépendant des
vecteurs précédents, alors, les vecteurs r0 , Ar0 , ...., Am−1 r0 sont linéairement
indépendants et par conséquent Kk est de dimension k, pour tout 1 6 k < m. En
particulier Km est de dimension m.
De plus, et d’après ce qui précède on a
Am r0 ∈ Km =⇒ Am+p r0 ∈ Km , ∀p > 0
Donc Km+p = Km , pour tout p > 0.
Par suite K1 · · · Km = Km+p , ∀p > 0, ce qui prouve que la suite d’espaces de
Krylov a atteint son point de stagnation en mmax = m, c’est-à-dire,
Kq = Kmmax , ∀q ≥ mmax
. 11
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
11 / 82
 Méthodes de Krylov Introduction

Théorème
La solution du système linéaire Ax = b appartient à l’espace affine x0 + Kmmax .

Les vecteurs r0 , Ar0 , ..., Ammax −1 r0 sont linéairement indépendants, et

max −1
mX
Ammax r0 = αk Ak r0
k=0

P−1
mmax
Si α0 = 0 alors, on a Ammax r0 = αk Ak r0
k=1
donc :
P−1
mmax
Ammax −1 r0 = αk Ak−1 r0
k=1

ce qui est absurde puisque les vecteurs r0 , Ar0 , ..., Ammax −1 r0 sont linéairement
indépendants. Donc α0 est forcément non nul.
Par suite
1 mmax
P−1
mmax
αk k 1 mmax α0 P−1
mmax
αk k
α0 A r0 = α0 A r0 =⇒ α0 A r0 = α0 r0 + α0 A r0
k=0 k=1
12
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
12 / 82
 Méthodes de Krylov Introduction

P−1
mmax
αk k 1 mmax
=⇒ −r0 − α0 A r0 + α0 A r0 = 0
k=1
P−1 αk k
mmax
1 mmax
=⇒ Ax0 − b − α0 A r0 + α0 A r0 = 0
k=1
P−1 αk k
mmax
1 mmax
=⇒ Ax0 − α0 A r0 + α0 A r0 = b
k=1
P−1 −αk k−1
mmax
=⇒ A(x0 + α0 A r0 + α10 Ammax −1 r0 ) = b
k=1
P−1 −αk k−1
mmax
En posant x̃ = x0 + α0 A r0 + α10 Ammax −1 r0
k=1

P−1
mmax
αk k−1 1 mmax −1
et comme α0 A r0 − α0 A r0 ∈ Kmmax (A, r0 )
k=1

Alors x̃ ∈ x0 + Kmmax (A, r0 )

13
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
13 / 82
 Méthodes de Krylov Introduction

Méthodes de projection

Soit à résoudre le système linéaire

Ax = b

où A ∈ Rn×n est une matrice inversible, généralement creuse de grande taille

(d’ordre pouvant aller jusqu’à 108 ), b ∈ Rn le second membre et le vecteur x est
la solution cherchée.

Principe
+ Le principe général des méthodes de projection consiste à projeter le
problème Ax = b sur un sous-espace Km de dimension m inférieure à n,
orthogonalement à un deuxième sous-espace Lm , également de dimension m,
appelé sous-espace des contraintes. Le système linéaire résultant est de petite
taille et donc facile à résoudre.

14
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
14 / 82
 Méthodes de Krylov Introduction

En général, on cherche à extraire une approximation x̃ = x0 + δ de la solution,

dans l’espace x0 + Km vérifiant la condition d’orthogonalité

r˜ = b − Ax̃ ⊥ Lm (Petrov-Galerkin)

Avec x0 donné initial dans Rn , r0 = b − Ax0 est le résidu initial et δ ∈ Km .

Cela peut être reformulé sous la forme générale :


x̃ = x0 + δ δ ∈ Km
< r0 − Aδ, v >= 0 ∀v ∈ Lm

Remarque
+ Lorsque Km = Lm , la solution approchée est exacte sur Km · On dit que la
projection est orthogonale. Dans le cas contraire, lorsque Km 6= Lm , la projection
est dite oblique.

15
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
15 / 82
 Méthodes de Krylov Introduction

Représentation matricielle
Nous allons formuler les conditions précédentes sous forme matricielle. Choisissons
deux bases V = [v1 , v2 , ..., vm ] et W = [w1 , w2 , ..., wm ] respectivement pour Km et
Lm .
La solution approximative peut s’écrire :

x̃ = x0 + Vy y ∈ Rm

r ⊥ Lm ), on a :
Pour la condition d’orthogonalité (˜

W T r˜ = 0 ⇒ W T (r0 − AVy ) = 0 ⇒ W T AVy = W T r0

La solution approchée peut donc s’exprimer sous la forme :

x̃ = x0 + V (W T AV )−1 W T r0

16
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
16 / 82
 Méthodes de Krylov Méthode d’Arnoldi

Méthode d’Arnoldi

Principe
La méthode d’Arnoldi est une méthode de projection orthogonale sur un
sous-espace de Krylov, généralement appliquée aux matrices non symétriques.
Cette méthode a été introduite par Arnoldi en 1950, dans le but de réduire une
matrice A ∈ Mn (K) sous forme Hessenberg supérieure de taille m × m si
l’algorithme arrive à la m-ième itération m  n,
 
h11 h12 . . . h1m−1 h1m
 h21 h22 . . . h2m−1 h2m 
.. .. 
 
 ..
Hm = (hij ) =  0 h32
 . .  ∈ Mm (K)
. 
 .. . . . . .
. .
. 
 . . . . . 
0 . . . 0 hmm−1 hmm
et de construire une matrice orthonormale Vm de taille n × m dont les vecteurs
colonnes v1 , ..., vm forment une base orthonormée du sous-espace de Krylov
Km (A, v1 ) = [v1 , Av1 , ..., Am−1 v1 ].

17
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
17 / 82
 Méthodes de Krylov Méthode d’Arnoldi

Algorithme d’Arnoldi classique

r0
v1 ← avec β = kr0 k
β
for j = 1, 2, · · · , m do
X j
w ← Avj − hi,j vi avec hi,j =< Avj , vi >
i=1
hj+1,j ← kw k
if hj+1,j = 0 then
stop
else
w
vj+1 =
hj+1,j
end if
end for

18
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
18 / 82
 Méthodes de Krylov Méthode d’Arnoldi

L’algorithme d’Arnoldi conduit aux résultats suivants :

A · v1 = h1,1 v1 +h2,1 v2
A · v2 = h1,2 v1 +h2,2 v2 + h3,2 v3
.. .. ..
. . .
A · vm = h1,m v1 +h2,m v2 + · · · + hm,m vm + hm+1,m vm+1

Matriciellement, en posant Vm = [v1 , v2 , · · · , vm ], on obtient

 
h1,1 h1,2 h1,3 ··· h1,m
 h2,1 h2,2 h2,3 ··· h2,m 
.. .. 
 
 .. .. >
A · Vm = Vm ·  0 . . . .  + hm+1,m (vm+1 · em ) (1)
 .. . .. .. .. .. 
 . . . . 
0 0 0 hm−1,m hm,m

où em est le m-ième vecteur de la base canonique de Rm .

19
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
19 / 82
 Méthodes de Krylov Méthode d’Arnoldi

Donc
>
AVm = Vm Hm + hm+1,m (vm+1 · em )
= Vm+1 Hm
e

La matrice Hem correspond à la matrice Hm définie ci-dessus, augmentée de la ligne

m + 1, le terme hm+1,m étant le m-ième et seul élément non nul de cette ligne.
 
h11 h12 . . . h1m−1 h1m
 h21 h22 . . . h2m−1 h2m 
. ..
 
 0 h32 . .
 . .. 
. 
Hem = 
 .. . .

.. .. .. .. 
 .
 . . 

 0 . . . 0 hmm−1 hmm 
0 ... 0 0 hm+1m
Par conséquant
VmT AVm = Hm

20
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
20 / 82
 Méthodes de Krylov Méthode d’Arnoldi

Résolution du système linéaire par Arnoldi

Rappel
la solution approchée xm du système linéaire Ax = b appartient à l’espace affine
x0 + Km . Soit Vm une base de Km , alors xm s’écrit sous la forme :

xm = x0 + Vm ym

où ym est un vecteur de dimension m .

Théorème
La solution de l’équation Ax = b est obtenue en déterminant le vecteur ym
solution du système
Hm ym = βe1
(où β est la norme du résidu initial β =k r0 k2 =k b − Ax0 k2 et e1 est le premier
vecteur de la base canonique de Rm .)
La solution approchée xm est donnée par :
−1
xm = x0 + Vm Hm k r0 k2 e1
21
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
21 / 82
 Méthodes de Krylov Méthode d’Arnoldi

Preuve
Soit rm = b − Axm le résidu associé à xm . et r0 = b − Ax0 le résidu correspondant
à x0 .
Comme xm = x0 + Vm ym , alors

rm = b − A(x0 + Vm ym )
= b − Ax0 − AVm ym
= r0 − AVm ym

D’où

VmT rm = VmT r0 − VmT AVm ym

= VmT r0 − Hm ym

D’après la condition de Galerkin on a rm ⊥ Vm , d’où VmT rm = 0

22
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
22 / 82
 Méthodes de Krylov Méthode d’Arnoldi

donc
VmT r0 − Hm ym = 0
⇒ VmT r0 = Hm ym
−1 T
⇒ ym = H m Vm r0 (∗)
Or
v1 = r0 / k r0 k2 et v1 = Vm e1 ⇒ r0 = Vm e1 k r0 k2
Donc (∗) devient :
−1 T
ym = H m Vm Vm e1 k r0 k2
−1
= Hm k r0 k2 e1

Finalement, on obtient
−1
xm = x0 + Vm Hm k r0 k2 e1

23
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
23 / 82
 Méthodes de Krylov Méthode d’Arnoldi

Algorithme d’Arnoldi modifié

r0
v1 ← où β = kr0 k
β
for j = 1, 2, · · · , m do
w ← Avj
for i = 1, 2, · · · , j do
hi,j ←< w , vi > et w ← w − hi,j vi
end for
hj+1,j ← kw k
if hj+1,j = 0 then stop
else
w
vj+1 =
hj+1,j
end if
end for

24
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
24 / 82
 Méthodes de Krylov Méthode d’Arnoldi

Remarque
L’algorithme d’Arnoldi modifié est mathématiquement équivalent à celui d’Arnoldi
classique, par contre il présente une importante stabilité numérique.

25
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
25 / 82
 Méthodes de Krylov Méthode de Lanczos

Méthode de Lanczos

Cornelius Lanczos (1893- 1974), est un mathématicien et physicien hongrois, il

a enseigné les mathématiques en Hongrie, en Allemagne, aux États-Unis et en
Irlande. Ses travaux ont eu un impact profond sur le développement de la science
au XXe siècle, parmi lesquelles nous citons :
- L’algorithme de Lanczos, pour trouver les valeurs propres de matrices
symétriques.
- L’approximation de Lanczos, pour la fonction gamma (fonction complexe).
- La méthode du gradient conjugué pour la résolution itérative des systèmes
d’équations linéaires de grande taille.

Principe
La méthode de Lanczos est l’une des méthodes itératives les plus utilisées pour
la résolution de grands systèmes linéaires creux, son processus peut être vu
comme une simplification du processus d’Arnoldi lorsque la matrice associée au
système est symétrique, il consiste à construire une base orthonormée Vm de
l’espace de Krylov Km (A, v1 ) ayant un vecteur v1 comme vecteur de départ. et Tm
une matrice symétrique tridiagonale de telle sorte que Tm = VmT AVm .
26
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
26 / 82
 Méthodes de Krylov Méthode de Lanczos

Théorème
La matrice Hm calculée par la méthode d’Arnoldi appliquée à une matrice
symétrique A est symétrique tridiagonale.
c’est-à-dire

hij = 0 pour 1 6 i < j − 1
hj+1,j = hj,j+1 pour j = 1, 2, · · · , m − 1

Preuve : D’après l’algorithme d’Arnoldi on a

hij =< Avj , vi >=< vj , Avi > pour i = 1, 2, · · · j

i+1
X
et puisque Avi = hki vk .
k=1
Alors, pour i = 1, · · · , j − 2 on a k varie de 1, · · · , j − 1
Donc hij =< vj , Avi >= 0 (car < vj , vk >= 0 ∀k = 1, 2, · · · j − 1) ça
d’une part.
D’autre part on a ∀j = 1, 2, · · · , m − 1 :

hj+1,j =< Avj , vj+1 >=< vj , Avj+1 >=< Avj+1 , vj >= hj,j+1 .
27
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
27 / 82
 Méthodes de Krylov Méthode de Lanczos

Nous réécrivons donc la matrice de Hessenberg sous la forme :

 
α1 β2 0 ··· 0
 .. .. 
 β2
 α2 β3 . . 

Tm =  0
 .. 
, βi > 0, i = 2, ..., m
 β3 α3 . 0 
. .. .. ..
 ..

. . . βm 
0 ··· 0 βm αm

Remarque
Le processus d’Arnoldi à été donc simplifié car seuls les vecteurs vj et vj−1 sont
nécessaires à la construction de vj+l de telle sorte que

βj+1 vj+l = Avj − βj vj−1 − αj vj

28
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
28 / 82
 Méthodes de Krylov Méthode de Lanczos

Algorithme de Lanczos symétrique

Début
1. Choix d’un vecteur initial x0 et r0 = b − Ax0
2. β =k r0 k2 , v1 = r0 /β, on fixe β1 = 1 et v0 = 0
3. Pour j = 1, 2, .....m Faire :
i. αj = hAvj , vj i
ii. wj = Avj − βj vj−1 − αj vj
iii. βj+1 =k wj k2
iv. si βj+1 = 0 Alors stop
sinon
v. vj+1 = wj /βj+1
vi. Fin si
4. Fin Pour
5. Tm = tridiag (βj+1 , αj , βj+1 ) et Vm = [v1 , v2 ....., vm ]
Fin

29
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
29 / 82
 Méthodes de Krylov Méthode de Lanczos

Avj = βj+1 vj+1 + βj vj−1 + αj vj ∀j = 1, 2, · · · , m

Grace à cette relation de récurrence à trois termes On a :

A · v1 = α1 v1 +β2 v2
A · v2 = β2 v 1 +α2 v2 + β3 v3
.. .. ..
. . .
A · vm = βm vm−1 +αm vm + βm+1 vm+1

Donc
T
AVm = Vm Hm + βm+1 vm+1 em
où em est le m-ième vecteur canonique de l’espace Rm et vm+1 sera le prochain
vecteur colonne de la base de Km .

30
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
30 / 82
 Méthodes de Krylov Méthode de Lanczos

Montrons effectivement que la matrice Vm est orthogonale, c’est-à-dire montrons

que l’assertion suivante est vraie :
(Pk ) vi ⊥ vk , ∀i = 1, 2, . . . , k − 1
On procède par réccurence sur k :

1 1
< v1 , v2 >= < v1 , A v1 − α1 v1 >= (−α1 < v1 , v1 > + < v1 , Av1 >) = 0
β2 β2

(Pk+1 ) vi ⊥ vk+1 , ∀i = 1, 2, . . . , k est-elle vraie ?

Si i = k, on a
1
< vk+1 , vk > = < w k , vk >
βk+1
1
= < A vk − βk vk−1 − αk vk , vk >
βk+1
 
1 
= < A vk , vk > −αk 
βk+1 | {z }
=0
=0
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
31
31 / 82
 Méthodes de Krylov Méthode de Lanczos

Si i = k − 1, on a

1
< vk+1 , vk−1 > = < A vk − βk vk−1 − αk vk , vk−1 >
βk+1
1
= (< A vk , vk−1 > −βk )
βk+1
1
= (< vk , A vk−1 > −βk )
βk+1
1
= (< vk , βk vk + αk−1 vk−1 + βk−1 vk−2 > −βk )
βk+1

1 
= βk < vk , vk > +αk−1 < vk , vk−1 > +βk−1 < vk , vk−2 > −β
βk+1 | {z } | {z }
=0 =0
1
= (βk − βk )
βk+1
= 0.

32
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
32 / 82
 Méthodes de Krylov Méthode de Lanczos

Pour i < k − 1 on a

1
< vk+1 , vi > = < A vk − βk vk−1 − αk vk , vi >
βk+1
1
= < A v k , vi >
βk+1
1
= < vk , A v i >
βk+1
1
= < vk , βi+1 vi+1 + αi vi + βi vi−1 >
βk+1
= 0 car i + 1 < k

Ce qui prouve que (Pk ) est vraie pour ∀i = 1, 2, ..., k − 1

33
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
33 / 82
 Méthodes de Krylov Méthode de Lanczos

On sait que la solution de Ax = b est dans le sous-espace affine x0 + Km (A, r0 ).

Donc la solution approchée xm s’écrit sous la forme :
xm = x0 + Vm ym ( ? ) avec ym ∈ Rm

Or rm = b − Axm
= b − A(x0 + Vm ym )
= r0 − AVm ym

et d’après la condition de Petrov Galerkin ( rm ⊥Vm ), on obtient :

VmT .rm = 0 ⇒ VmT r0 − VmT AVm ym = 0

| {z }
Tm

⇒ Tm ym = VmT r0
⇒ ym = Tm−1 VmT r0

par conséquent (?) devient :

xm = x0 + Vm Tm−1 VmT r0
34
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
34 / 82
 Méthodes de Krylov Méthode de Lanczos

VmT r0 = VmT βv1

= β VmT v1
| {z }
= βe1
=k r0 k e1

On obtient finalement

xm = x0 + Vm Tm−1 k r0 k e1

35
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
35 / 82
 Méthodes de Krylov Méthode de Lanczos

Méthode de Lanczos pour les systèmes linéaires

Début
1. Choix d’un vecteur x0 et r0 = b − Ax0
2. β =k r0 k2 , v1 = r0 /β et β1 v0 = 0
3. Pour j = 1, 2, .....m Faire
i. αj = hAvj , vj i
ii. wj = Avj − βj vj−1 − αj vj
iii. βj+1 =k wj k2
iv. si βj+1 = 0 Alors stop
sinon
v. vj+1 = wj /βj+1
vi. FinSi
4. FinPour
5. Tm = tridiag (βj+1 , αj , βj+1 ) et Vm = [v1 , v2 ....., vm ]
6. Résoudre Tm ym = (βe1 )
7. xm = x0 + Vm ym
Fin

36
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
36 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Lanczos bi-orthogonal
Introduction
La méthode de Lanczos symétrique peut être généralisée au cas non symétrique
appelée méthode de biorthogonalisation de Lanczos, également appelée méthode
de Lanczos bi-orthogonal.

Principe
L’idée de Lanczos bi-orthogonal est d’utiliser simultanément deux procédures de
Lanczos pour construire une paire de bases Vm = {v1 , v2 , ...,
 vm } et
1 si i = j
Wm = {w1 , w2 , ..., wm } bi-orthogonales (i.e viT wj = δij = ) pour
0 si i 6= j
les deux sous-espaces de Krylov :

Km (A, v1 ) = Vect v1 , Av1 , ...., Am−1 v1



et
Km (AT , w1 ) = Vect w1 , AT w1 , ...., (AT )m−1 w1


.
37
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
37 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Principe
L’objectif majeur qui vise à réaliser la procédure de Lanczos bi-orthogonal
consiste à transformer une matrice non symétrique A en une matrice tridiagonale
non symétrique Tm de telle sorte que

WmT AVm = Tm

Avec Vm , Wm ∈ Rn×m qui sont respectivement les bases de Km (A, v1 ) et

Km (AT , w1 ), satisfaisant WmT Vm = Im

38
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
38 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

A l’aide de l’algorithme de Lanczos on a :

AVm = Vm+1 T̃m et AT Wm = Wm+1 R̃m . (∗∗)

 
α1 β2
 δ2
 α2 β3 0 

 .. .. .. 
Où R̃m = 
 . . .  ∈ R(m+1)×m


 δm−1 αm−1 βm 

 0 δm αm 
δm+1
 
α1 δ2
 β2
 α2 δ3 0 

 .. .. .. 
T̃m = 
 . . .  ∈ R(m+1)×m


 βm−1 αm−1 δm 

 0 βm αm 
βm+1

39
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
39 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Avec TmT = Rm , en effet ;

TmT =(WmT AVm )T

=VmT Wm Rm
=Rm

Les deux équations (**) peuvent être reformulées aussi par :


Avj = βj vj−1 + αj vj + δj+1 vj+1
AT wj = δj wj−1 + αj wj + βj+1 wj+1

∀ j = 1, 2, · · · , m., avec β1 = δ1 = 0

40
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
40 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Algorithme Bi-Lanczos

L’algorithme de Lanczos bi-orthogonal ou simplement Bi-Lanczos se déroule

comme suit :

On commence les itérations par deux vecteurs de départ v1 , w1 , vérifiant

v1 w1
< v1 , w1 >= 1, (si < v1 , w1 >= γ 6= 0 alors v1 ←− p ; w1 ←− p ).
 |γ| sign(γ) |γ|
Av1 = α1 v1 + δ2 v2
Comme Alors α1 =< w1 , A · v1 > pour avoir
A> w1 = α1 w1 + β2 w2
< v2 , w1 >= 0 et < v1 , w2 >= 0 (bi-orthogonalité). Puis on pose :

z = A · v1 − α1 v1 = δ2 v2
t = A> · w1 − α1 w1 = β2 w2
p p
Si < z, t >= η, on prend δ2 = |η| et β2 = sign(η) |η| pour avoir
< v2 , w2 >= 1 et ainsi de suite on construit les autres éléments.

41
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
41 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Algorithme Bi-Lanczos

Algorithme
1. Choisir deux vecteurs v1 et w1 tel que < v1 , w1 >= 1
2. Poser β1 = δ1 = 0, w0 = v0 = 0
3. Pour j = 1, ..., m Faire
4. αj =< Avj , wj >
5. ṽj+1 = Avj − αj vj − βj vj−1
6. AT wj − αj wj − δj wj−1
w̃j+1 = p
7. δj+1 = | < ṽj+1 , w̃j+1 > |. Si δj+1 = 0 stop
< ṽj+1 , w̃j+1 >
8. βj+1 =
δj+1
w̃j+1
9. wj+1 =
βj+1
ṽj+1
10. vj+1 =
δj+1
11. Fin

42
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
42 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Théorème
Si l’algorithme ne s’arrête pas avant le rang m, alors les vecteurs (vi )16i6m , et ,
(wi )16i6m , sont biorthogonaux i.e.,

< vj , wi >= δij 1 ≤ i, j ≤ m

43
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
43 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Preuve
Par hypothèse < v1 , w1 >= 1. Supposons que les vecteurs v1 , . . . , vj et w1 , . . . , wj
sont biorthogonaux, et montrons que les vecteurs v1 , . . . , vj+1 et w1 , . . . , wj+1
restent biorthogonaux.
c’est-à-dire, vérifions les assertions suivantes :

 (ṽj+1 , wi ) = 0 pour ∀i ≤ j.
(w̃j+1 , vi ) = 0 pour ∀i ≤ j.
(vj+1 , wj+1 ) = 1


44
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
44 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Si i = j,on a :

< vj , w̃j+1 > =< vj , AT wj − αj wj − δj wj−1 >

= < vj , AT wj > −αj < vj , wj > −δj < vj , wj−1 >
| {z } | {z } | {z }
αj 1 0

= αj − αj
=0

Pour le cas i < j, deux sous cas peuvent être distingués i = j − 1 et i < j − 1
Si i = j − 1

< vj−1 , w̃j+1 > =< vj−1 , AT wj > −αj < vj−1 , wj > −δj < vj−1 , wj−1 >
| {z } | {z }
0 1
=< Avj−1 , wj > −δj
= βj−1 < vj−2 , wj > −αj−1 < vj−1 , wj > +δj < vj , wj > −δj
= δj − δj
=0
45
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
45 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Si i < j − 1

< vi , w̃j+1 > =< vi , AT wj > −αj < vi , wj > −δj < vi , wj−1 >
| {z } | {z }
0 0
=< vi , AT wj >
= δj < vi , wj−1 > +αj < vi , wj > +βj+1 < vi , wj+1 >
| {z } | {z } | {z }
0 0 0
=0

De manière analogue, on peut montrer que < ṽj+1 , wi >= 0 ∀i 6 j.

Pour conclure, il reste à montrer que < vj+1 , wj+1 >= 1

ṽj+1 w̃j+1
< vj+1 , wj+1 > =< , >
βj+1 δj+1
1
= < ṽj+1 , w̃j+1 >
βj+1 δj+1
1
= < ṽj+1 , w̃j+1 >
< ṽj+1 , w̃j+1 >
=1 46
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
46 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Théorème
(vi )16i6m , et , (wi )16i6m , forment respectivement une base de Km (A, v1 ) et
Km (AT , w1 ) , et on a les relations suivantes :
T
AVm = Vm Tm + δm+1 vm+1 em

AT Wm = Wm TmT + βm+1 wm+1 em

T

WmT AVm = Tm

47
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
47 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Résolution de Ax = b

Résolution de Ax = b
Pour résoudre le système linéaire Ax = b, on procède de la manière suivante :
1 En partant d’une donnée initiale x0 , on calcule le résidu correspondant
r0
r0 = b − Ax0 , puis β = kr0 k et v1 =
β
2 On construit w1 vérifiant < w1 , v1 >= 1
3 On applique Lanczos bi-orthogonal pour construire Vm et Wm les bases
bi-orthogonales de Km (A, v1 ) et Km (A> , w1 )
4 On effectue la projection sur Km (A, v1 ) orthogonalement à Km (A> , w1 ), tel
que : (
xm ∈ x0 + Km (A, v1 )
rm ⊥ Km (A> , w1 )

48
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
48 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Résolution de Ax = b
y m ∈ Rm

xm = x0 + Vm ym ;
⇔
rm = r0 − AVm ym ⊥ Wm
D’où

Wm> r (m) = 0 ⇒ Wm> r0 − Tm ym = 0

⇒ βe1 − Tm y (m) = 0 ( à vérifier)
⇒ ym = Tm−1 βe1

où e1 désigne le premier vecteur unitaire dans Rm . Par suite

xm = x0 + Vm Tm−1 βe1

De même pour le calcul de l’approximation du système A> x = b 0

49
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
49 / 82
 Méthodes de Krylov Lanczos Bi-orthogonal

Avantages-Inconvénients
X L’algorithme de Lanczos bi-orthogonal présente les avantages suivants :
- Il s’appuie sur une récurrence à trois termes, donc ce procédé est beaucoup
moins coûteux en place mémoire et en temps de calcul.
-Il permet de résoudre simultanément les systèmes linéaires associés à A et à
sa transposé AT .
X Son inconvénient principal est que la convergence peut être assez irrégulière.
- Un autre inconvénient : il n’est pas facile d’exploiter la forme tridiagonale (
l’algorithme QR, qui est stable et facile à utiliser, augmentera la matrice jusqu’à
la forme de Hessenberg).

50
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
50 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

Méthode de gradient conjugué

Principe
La méthode du gradient conjugué est couramment utilisée pour la résolution de
grands systèmes linéaires creux, dont la matrice associée est symétrique, définie
positive. Il s’agit d’une méthode itérative qui consiste à partir d’un vecteur initial
x0 est à déterminer à chaque étape k un vecteur dk ∈ Kk+1 (A, r0 ) appelé direction
de descente et un scalaire αk permettant de calculer xk+1 à partir de xk tel que :

xk+1 = xk + αk dk

+ Son principe se base donc sur la recherche de directions successives

permettant d’atteindre la solution exacte du système étudié.

Théorème
Les directions de descente di sont A-conjuguées ( ou conjuguées par rapport
à A), c’est-à-dire elles vérifient la condition suivante :

diT Adj = 0 pour tout i 6= j

51
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
51 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

Théorème (Exercice)
Tant que les résidus rk sont non nuls, la méthode du gradient conjugué vérifie :
r0 = b − Ax0 , d0 = r0
1 rk+1 = rk − αk Adk , αk ∈ R
2 Kk+1 (A, r0 )=Vect{d0 , d1 , · · · , dk } = Vect{r0 , r1 , · · · , rk }
3 < rk , ri >= 0 ; < rk , di >= 0 , pour 0 ≤ i ≤ k − 1
4 Adi ∈ Kk+1 (A, r0 ) ∀i ≤ k − 1
5 dk+1 = rk+1 + βk dk , βk ∈ R
6 < Adk , rk > = < Adk , dk >
< rk , rk >
7 αk =
< Adk , dk >
< rk+1 , rk+1 >
8 βk =
< rk , rk >

52
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
52 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

Preuve
1 On sait que xk ∈ x0 + Kk (A, r0 ), ou de manière équivalente

xk+1 ∈ xk + Kk+1 (A, r0 )

Donc, il existe dk ∈ Kk+1 (A, r0 ) et αk ∈ R, vérifiant xk+1 = xk + αk dk ,

entraine b − Axk+1 = b − Axk − αk Adk
Par suite rk+1 = rk − αk Adk

2 Or d0 = r0 , alors Vect{d0 } = Vect{r0 }, donc la propriété est vraie au premier

rang 0.
Supposons que Kk (A, r0 )=Vect{d0 , d1 , · · · , dk−1 } = Vect{r0 , r1 , · · · , rk−1 }
et montrons que l’égalité est encore satisfaite au rang suivant k + 1.

On peut y distinguer deux cas :

1er cas : dim(vect{d0 , d1 , ..., dk−1 }) = dim(vect{d0 , d1 , ..., dk }) et
puisque les directions (di )i∈{0,1,...,k−1} sont linéairement indépendantes,
alors : dk = 0 d’où xk+1 = xk
Donc la convergence est assurée, par conséquent rk+1 = 0.
53
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
53 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

2ème cas : Si dim(Vect{d0 , d1 , ..., dk }) = dim(Vect{d0 , d1 , ..., dk−1 }) + 1 et

comme dk ∈ Kk+1 (A, r0 ) alors on a :

Vect{d0 , d1 , ..., dk } ⊆ Kk+1 (A, r0 )

et puisque dim(Vect{d0 , d1 , ..., dk }) = dim Kk+1 (A, r0 ), alors

Vect{d0 , d1 , ..., dk } = Kk+1 (A, r0 )

D’autre part on a : rk = rk−1 − αk−1 Adk−1 , avec rk−1 ∈ Kk (A, r0 ) et

dk−1 ∈ Kk (A, r0 ) =⇒ Adk−1 ∈ Kk+1 (A, r0 )
Donc rk ∈ Kk+1 (A, r0 )
entraine vect{r0 , r1 , ..., rk } ⊆ Kk+1 (A, r0 ) et comme
dim(vect{r0 , r1 , ..., rk }) = dim(Kk+1 (A, r0 ))
Alors

Vect{r0 , r1 , ..., rk } = Kk+1 (A, r0 )

Par conséquent Kk+1 (A, r0 ) = Vect{r0 , r1 , ..., rk } = Vect{d0 , d1 , ..., dk }

54
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
54 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

3) D’après la condition de Galerkin (rk ⊥ Kk (A, r0 )) , et puisque

Kk (A, r0 ) = Vect{r0 , r1 , ..., rk−1 } = Vect{d0 , d1 , ..., dk−1 } on a

< rk , ri >= 0 ; < rk , di >= 0 ∀i ≤ k − 1

4) ∀i = 0, 1, 2, · · · , k − 1, on a di ∈ Ki+1 (A, r0 )
D’où Adi ∈ Ki+2 (A, r0 ) ⊂ Kk+1 (A, r0 ), cette inclusion est vraie car si i ≤ k − 1,
on a i + 2 ≤ k + 1.
Donc on a bien Adi ∈ Kk+1 (A, r0 ), pour i ≤ k − 1

55
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
55 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)


d0 = r0
5) Montrons que
dk+1 = rk+1 + βk dk
Comme Vect{d0 , d1 , · · · , dk } = Vect{r0 , r1 , · · · , rk }= Kk+1 (A, r0 ) , alors
d0 ∈ Vect{r0 }, on peut donc prendre d0 = r0 .
k+1
P
Puisque rk+1 ∈ Kk+2 (A, r0 )= Vect{d0 , d1 , · · · , dk+1 }, alors rk+1 = σj dj
j=1

Soit i ∈ {0, 1, 2, · · · , k − 1}, on a :

k+1
P
< rk+1 , Adi > = σj < dj , Adi >
j=1
= σi < di , Adi > car < dj , Adi >= 0 si i 6= j
< rk+1 , Adi >
Ainsi σi = = 0 pour i ≤ k − 1 car Adi ∈ Kk+1 (A, r0 ), et
< di , Adi >
rk+1 ⊥ Kk+1 (A, r0 ), avec < di , Adi > > 0 pour A symétrique et définie
positive.

56
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
56 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

On peut donc écrire : rk+1 = σk dk + σk+1 dk+1

Si σk+1 = 0 alors rk+1 = σk dk , et puisque dk ∈ Kk+1 (A, r0 ) alors
rk+1 ∈ Kk+1 (A, r0 ) ce qui est absurde car rk+1 ⊥ Kk+1 (A, r0 ), on a donc
forcément σk+1 6= 0, prenons σk+1 = 1 et σk = βk
finalement on obtient dk+1 = rk+1 − βk dk

6) Montrons maintenant que < Adk , rk > = < Adk , dk >

< Adk , rk > = < Adk , dk + βk−1 dk−1 >
= < Adk , dk > +βk−1 < Adk , dk−1 >
= < Adk , dk >
Sachant que les dk sont A-conjugués, c’est-à-dire diT Adk = 0 si k 6= i.

57
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
57 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

7) On sait que les rk sont orthogonaux entre eux deux à deux, donc
< rk+1 , rk >= 0
Comme rk+1 = rk − αk Adk alors

< rk − αk Adk , rk >=< rk , rk > −αk < Adk , rk >= 0

< rk , rk >
⇐⇒ αk =
< Adk , rk >

< rk , rk >
Ainsi αk =
< Adk , dk >

58
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
58 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

< rk+1 , rk+1 >

8) Montrons que βk = −
 < rk , rk >
rk+1 = rk − αk Adk
Or et puisque < Adk , dk+1 >= 0, alors on a :
dk+1 = rk+1 − βk dk
< rk+1 , rk+1 > = < rk − αk Adk , rk+1 >
= < rk , rk+1 > −αk < Adk , rk+1 >
= −αk < Adk , rk+1 >
= −αk < Adk , dk+1 + βk dk >
= −βk αk < Adk , dk >
< rk , rk >
= −βk < Adk , dk >
< Adk , dk >
< rk+1 , rk+1 >
Par conséquent βk = −
< rk , rk >

59
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
59 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

GC-Algorithme
Initialisation :
x0 donné ;
k = 0;
r0 = b − Ax0 ;
d0 = r0
While rk 6= 0 ;
r T rk
αk = Tk ;
rk Adk
xk+1 = xk + αk dk
rk+1 = rk − αk Adk
r T rk+1
βk = − k+1T ;
rk rk
dk+1 = rk+1 − βk dk
k =k +1
end while

60
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
60 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

Convergence de GC

Proposition
Tant que le résidu rk 6= 0, les coefficients αk et βk existent et sont uniques.

Preuve
T
rk+1 rk+1
+ Si rk 6= 0, alors βk = − est bien défini et est unique.
rkT rk
+ Si rk 6= 0, alors dk 6= 0, et puisque la matrice A est définie positive alors
r T rk
dkT Adk 6= 0, ce qui prouve que αk = Tk existe et est unique.
dk Adk

61
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
61 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

Convergence de GC
Proposition
La méthode de gradientr conjugué à une convergence strictement monotone, tel
1
que kek+1 kA ≤ 1 − kek kA < kek kA
cond2 (A)
où [Link] désigne la norme associée à la matrice symétrique définie positive A.

Preuve
On a rk = b − Axk = Ax − Axk = Aek , avec ek = x − xk
D’où T
rk+1 A−1 rk+1 = ek+1
T
Aek+1 ⇐⇒ kek+1 k2A = krk+1 k2A−1
Or rk+1 = rk − αk Adk , alors

kek+1 k2A = (rk − αk Adk )T A−1 (rk − αk Adk )

= rkT A−1 rk − 2αk rkT dk + αk2 dkT Adk
= kek k2A − 2αk rkT dk + αk2 dkT Adk

rkT rk
Comme αk = , et rkT dk = rkT rk (à vérifier)
dkT Adk 62
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
62 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

Alors kek+1 k2A = kek k2A − αk krk k22 . Nous allons minorer k rk k22 , et αk
k ek k2A
k ek k2A = rkT A−1 rk 6k A−1 k2 k rk k22 , d’où k rk k22 >
k A−1 k2

rkT rk k rk k22
Il est clair que αk = = , on a
dkT Adk dkT Adk

dkT Adk = (rk − βk−1 dk−1 )T Adk = rkT Adk

= rkT A(rk − βk−1 dk−1 )
= rkT Ark − βk−1 rkT Adk−1
= rkT Ark − βk−1 (dk + βk−1 dk−1 )T Adk−1
2
= rkT Ark − βk−1 T
dk−1 Adk−1
2
Comme −βk−1 T
dk−1 Adk−1 6 0, alors dkT Adk 6 rkT Ark 6k A k2 k rk k22 , entraine

1
αk >
k A k2 63
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
63 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

Ainsi
k ek k2A k ek k2A
αk k rk k22 > =
k A−1 k2 k A k2 cond2 (A)
On en déduit donc
k ek k2A
k ek+1 k2A 6k ek+1 k2A −
cond2 (A)
finalement, on obtient
s
1
kek+1 kA ≤ 1− kek kA < kek kA
cond2 (A)

Remarque ? ? ?
+ La vitesse de convergence du GC dépendrdu conditionnement de la matrice A.
1
+ Si le cond2 (A) est grand, le coefficient 1 − devient proche de 1
cond2 (A)
et la convergence risque d’être lente.

64
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
64 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

Théorème
La méthode de gradient conjugué converge en un nombre fini d’itérations au
plus égal à la dimension n du système linéaire.
Autrement dit : pour tout choix de x0 , si r0 = b − Ax0 6= 0, alors

∃k ∈ {0, 1, · · · , n} tel que Ax k = b

Preuve
Supposons par l’absurde que l’algorithme ne converge pas en moins de n
itérations, alors r0 6= 0, r1 6= 0, ..., rn 6= 0, donc si rn 6= 0 le vecteur dn 6= 0, on a
donc une famille {d0 , d1 , ..., dn } A-orthogonale de vecteurs non nuls, et donc une
famille libre de n + 1 éléments de Rn ce qui est absurde, car il ne peut exister dans
un espace de dimension n plus de n vecteurs linéairement indépendants. Ainsi
l’algorithme converge vers x ∗ en au plus n itérations.

65
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
65 / 82
 Méthodes de Krylov Méthode de gradient conjugué (GC)

Exercice
T
Appliquer la méthode de gradient conjugué, avec 
donnée initiale x0 =
 (0, 0, 0)
2 −1 0
pour résoudre le système linéaire Ax = b où A =  −1 2 −1  ;
0 −1 2
 
−1
b= 2 
−1

     
−1 −3 1
3 1 1
r0 = d0 =  2  ; α0 = ; x1 =  6  . r1 =  1  ;
10 10 5
−1   −3 1
9
−1 1  5
β0 = ; d1 = 12  ; α1 =
50 50 3
  9 
0 0
x2 =  1  ; r2 =  0  La méthode a donc convergé à la 2ème itération.
0 0
66
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
66 / 82
 Méthodes de Krylov Lien entre Gradient conjugué et Lanczos

Lien entre Gradient conjugué et Lanczos

Dans ce qui suit, nous allons montrer que l’algorithme de CG n’est rien d’autre
qu’une variante de l’algorithme de Lanczos.

La solution approchée obtenue par la méthode de Lanczos sur Km est donnée par

xm = x0 + Vm ym
ym = Tm−1 βe1
 
α1 β2
 β2 α2 β3 
 
 .. .. .. 
 . . . 
Avec Tm =  
 .. .. 

 . . βm−1 

 βm−1 αm−1 βm 
βm αm

67
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
67 / 82
 Méthodes de Krylov Lien entre Gradient conjugué et Lanczos

Soit Tm = Lm Um la décomposition LU de Tm où Lm est bidiagonale inférieure à

diagonale
 unité, et Um est bidiagonale supérieure.
  
1 η1 β2
 λ2 1   η 2 β3 
   

 . .
.. .. 
 
 . .
.. .. 

Tm =   ×  
 . .. . .. 
  . .. . .. 

   
 λm−1 1   ηm−1 βm 
λm 1 ηm
| {z } | {z }
Lm Um
βm
Avec λm = , et ηm = αm − λm βm , en effet ;
ηm−1

68
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
68 / 82
 Méthodes de Krylov Lien entre Gradient conjugué et Lanczos

On a Tm (m, m − 1) = βm d’après la forme de T

m 
0
 .. 
 . 

 
et T (m, m − 1) = 0 0 ··· 0 λm 1  0  = λm ηm−1 d’après
 
 βm−1 
 
 ηm−1 
0
la décomposition Lm Um de Tm .
βm
D’où βm = λm ηm−1 ce qui implique que λm = .
ηm−1

69
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
69 / 82
 Méthodes de Krylov Lien entre Gradient conjugué et Lanczos

Or Tm (m, m) = αm d’après la forme de Tm  

0
  0
..
 
 
.


et Tm (m, m) = 0 0 · · · 0 λm 1   = λm βm + ηm d’après la

 0 
 
 βm 
ηm
décomposition Lm Um de Tm . Donc αm = λm βm + ηm
Par conséquent ηm = αm − λm βm .

Cette factorisation nous permet de réécrire l’approximation de la façon suivante,

−1 −1
xm = x0 + Vm Um Lm (βe1 ).
−1
Pm = Vm Um
Posons −1
zm = Lm (βe1 )

Ainsi xm = x0 + Pm zm
.
70
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
70 / 82
 Méthodes de Krylov Lien entre Gradient conjugué et Lanczos

Notons par pm la dérnière colonne de Pm , et montrons que

−1
pm = η m [vm − βm pm−1 ], avec vm = Vm (:, m)
−1
On a Pm = Vm Um , donc Vm = Pm Um  
0
 0 
..
 
 
 . 
d’où vm = Vm (:, m) = Pm Um (:, m) = [p1 , p2 , · · · · · · , pm−1 , pm ] 
 ..
.


 . 

 βm 
ηm
1
Donc, vm = βm pm−1 + ηm pm . Par suite pm = [vm − βm pm−1 ] .
ηm

Avec βm est un coefficient calculé par l’algorithme de Lanczos, tandis que ηm

résulte de la factorisation LU de Tm

71
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
71 / 82
 Méthodes de Krylov Lien entre Gradient conjugué et Lanczos

   
zm (1 : m − 1) zm−1
Posons zm = =
ξm ξm

et puisque Lm zm = βe1 , alors d’après la structure de Lm , on trouve

ξm = −λm ξm−1 , et ξ1 = β =k r0 k2

Avec ξi étant la i-ième composante du vecteur zm .

Finalement, la suite récurrente peut encore

 s’exprimer de la façon suivante :
  zm−1
xm = x0 + Pm−1 pm
ξm
= x0 + Pm−1 zm−1 +ξm pm
| {z }
= xm−1 + ξm pm

72
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
72 / 82
 Méthodes de Krylov Lien entre Gradient conjugué et Lanczos

Il s’ensuit donc que l’approximation xm peut être mise à jour à chaque étape par
la relation
xm = xm−1 + ξm pm

Donc la relation obtenue est analogue à celle de la méthode de gradient

conjugué, et les directions de descentes pi sont A-conjuguées,
i-e.(Api , pj) = 0∀i 6= j. En effet ;
T −T T −1
Pm APm = Um Vm AVm Um
−T −1
= Um Tm Um
−T
= Um Lm

qui est une matrice triangulaire inférieure et symétrique, donc diagonale.

73
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
73 / 82
 Méthodes de Krylov Lien entre Gradient conjugué et Lanczos

Cela nous a permis décrire l’algorithme suivant, que l’on appelle la version directe
de l’algorithme de Lanczos (ou le D-Lanczos) pour les systèmes linéaires.

Algorithme de D-Lanczos
[Link] de r0 = b − Ax0 , ξ1 := β := ||r0 ||, et v1 := r0 /β
[Link] pose λ1 = β1 = 0, p0 = 0
[Link] m = 1, 2, . . .,until convergence Do
[Link] de w := Avm − βm vm−1 et αm = (w , vm )
[Link] m > 1 alors on calcule λm = ηβm−1 m
et ξm = −λm ξm−1
6.ηm = αm − λm βm
−1
[Link] = ηm (vm − βm pm−1 )
[Link] = xm−1 + ξm pm
[Link] xm converge alors Stop
10.w := w − αm vm
11.βm+1 = ||w ||2 , vm+1 = w /βm+1
[Link]

74
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
74 / 82
 Méthodes de Krylov Méthode de GMRES

Méthode de GMRES

Principe
La méthode GMRES (Generalized Minimum Residual) a été formulée par Saad
et Schultz en 1986 comme une méthode de projection oblique sur le sous-espace
de Krylov Km = {r0 , Ar0 , · · · , Am−1 r0 }, avec Lm = AKm pour résoudre les
systèmes linéaires non symétriques Ax = b. Le m-ième itéré xm de la méthode
GMRES est obtenu en minimisant le résidu r = b − Ax sur l’ensemble des x
appartenant à l’espace affine x0 + Km . On écrit donc :

kb − Axm k = min kb − Axk

x∈x0 +Km

75
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
75 / 82
 Méthodes de Krylov Méthode de GMRES

L’algorithme de GMRES repose sur le processus d’Arnoldi pour construire une

base orthonormale Vm = [v1 , v2 , ..., vm ] et la matrice de Hessenberg
Hm = VmT AVm de taille m. Nous obtenons ainsi la matrice Hem , construite à partir
de Hm en lui ajoutant une ligne supplémentaire dont le seul élément non nul est
hm+1,m , à la position (m + 1, m). Les vecteurs vi et la matrice Hem satisfont la
relation suivante : AVm = Vm+1 Hem .

76
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
76 / 82
 Méthodes de Krylov Méthode de GMRES

Principe de GMRES
La minimisation des résidus successivement acquis pour chaque itéré est la base
de cette méthode :

kb − Axm k = min kb − Axk

x∈x0 +Km

= minm kb − Ax0 − AVm y k

y ∈R
 
= minm kr0 − Vm+1 Hem y k AVm = Vm+1 Hem (Arnoldi)
y ∈R

= minm kβVm+1 e1 − Vm+1 Hem y k (β = kr0 k)

y ∈R

= minm kVm+1 (βe1 − Hem y )k (Vm+1 est orthogonale)

y ∈R

= minm kβe1 − Hem y k

y ∈R

La méthode revient donc à minimiser J (y ) = kβe1 − Hem y k

77
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
77 / 82
 Méthodes de Krylov Méthode de GMRES

Ainsi, ym ∈ Rm est la solution par les moindres carrés du système J (y ) de

dimension (m + 1) × m ( m  n) :

ym = arg minm J (y )
y ∈R

Remarque : L’inconvénient principal de GMRES est l’absence de récurrence.

Remarque
Pour résoudre ce problème nous allons utiliser la factorisation QR basée soit sur
les transformations de Householder, qui est très stable numériquement, soit sur les
rotations de Givens.

78
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
78 / 82
 Méthodes de Krylov Méthode de GMRES

QR factorization based on Givens rotations

On transforme la matrice de Hessenberg Hem en une matrice triangulaire
supérieure en utilisant une succession de rotations de Givens ;
La matrice de rotation se présente sous la forme :
 
1

 ··· 


 1 

 ci s i  , avec ci2 + si2 = 1

Ωi = 

 −s i c i



 1 

 ··· 
1

Remarque : les matrices Ωi sont orthogonales.

79
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
79 / 82
 Méthodes de Krylov Méthode de GMRES

Pour éliminer le coefficient hi+1,i , on multiplie la matrice Hem par Ωi , où :

hi,i hi+1,i
ci = q et si = q
2 + h2
hi,i 2 + h2
hi,i
i+1,i i+1,i

Le produit Ωi Hem conduit à une matrice dont le coefficient hi+1,i a été éliminé.
En répétant ce procédé, on construit une matrice triangulaire supérieure
 
h1,1 h1,2 · · · h1,m−1 h1,m
 0 h2,2 h2,m 
 .. ..
 
. . 
Ωm Ωm−1 ...Ω1 Hem =  .
 . .  = Rm
 e
 0 hm−1,m−1 hm−1,m 
 
 0 ··· 0 hm,m 
0 ··· 0 0

80
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
80 / 82
 Méthodes de Krylov Méthode de GMRES

À l’aide des rotations de Givens, on peut donc définir la décomposition suivante :

Hem = Qm Rem

Avec
Qm = (Ωm Ωm−1 · · · Ω1 )T


bem = Q T (βe1 )m

Notre problème de minimisation consiste donc à résoudre par remontée le système

triangulaire :
Rem ym = bem
Autrement dit : la solution ym de minm kβe1 − Hem y k est la solution de l’équation
y ∈R
Rem ym = bem

81
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
81 / 82
 Méthodes de Krylov Méthode de GMRES

Algorithme de GMRES
r0
1. Initialisation : x0 donné, calculer r0 = b − Ax0 , β =k r0 k2 et v1 =
β
2. Construction de la base d’Arnoldi par Gram-schmidt :
i. Pour j = 1,2,....,m
ii. wj = Avj
iii. Pour i=1,2,.....,j
iv. hi,j = viT wj
P j
v. wj = wj − hi,j vi
i=1
vi. Calculer hj+1,j =k wj k2 et vj+1 = wj hj+1,j
vii. Fin Pour
3. Etablir la solution approchée :
Calculer ym qui minimise kβe1 − H̃m y k en utilisant la factorisation QR de H̃m .
Poser xm = x0 + Vm ym et rm = b − Axm
4. Redémarrage :
Si k rm k2 ≤  (tolérance) stop, sinon
r0
Poser x0 := xm , r0 := rm , β =k r0 k2 , et v1 = et retourner à l’étape 2
β
82
Résolution itérative de systèmes linéaires de grande taille par
11 décembre
des méthodes
2018de Krylov
82 / 82