Chapitre 1

Introduction
Dans ce cours, nous nous intéressons à la résolution de grands systèmes linéaires :
Ax = b

A est une matrice donnée, elle peut être issue de la discrétisation d'équations aux dérivées partielles. Le
vecteur b appelé second membre est aussi donné. Le vecteur x est l'inconnue qu'on cherche à calculer. On
suppose ici que la matrice A est toujours inversible, et donc que ce système possède une unique solution.
On a deux grandes classes de matrices
 Matrices pleines (souvent issues d'une résolution par équations intégrales). Ces matrices sont coû-
teuses à stocker, on dispose souvent d'un produit matrice vecteur, et la matrice n'est pas stockée.
 Matrices creuses (souvent issues d'une discrétisation par éléments nis ou volumes nis) de grande
taille.
Pour résoudre ces systèmes linéaires, on a deux classes de méthodes :
 Méthodes directes : par exemple méthodes de Cholesky ou LU. Ce sont des méthodes qui ont été
vues dans les cours précedents. Pour les matrices creuses, on obtient un phénomène de remplissage :
un zéro initial de la matrice A peut devenir un élément non-nul de L (ou U ). En conséquence, les
facteurs L et U sont beaucoup plus volumineux que la matrice initiale. Pour les matrices pleines,
ces méthodes imposent de stocker la matrice. Globalement, ce sont des méthodes qui sont coûteuses
en stockage mémoire.
 Méthodes itératives : par exemple gradient conjugué, Gmres. Ce sont des méthodes qui construisent
une suite x qui converge vers la solution x. Ce sont des méthodes peu coûteuses en espace mémoire
car on n'est pas obligé de stocker la matrice (on utilise une fonction calculant le produit matrice
n

vecteur), et le nombre de vecteurs utilisés pour calculer la suite x est réduit. Elles sont aussi souvent
plus facilement parallélisables que les méthodes directes. Le problème majeur de ces méthodes est
n

que la convergence peut être parfois assez lente.

L'objet de ce cours est d'étudier les méthodes itératives qui sont vraiment intéressantes pour résoudre de
grands systèmes linéaires, notamment pour des simulation numériques 3-D. Pour pallier au problème de
convergence, il est courant de préconditionner le système linéaire, ce qui consiste à résoudre :
M −1 Ax = M −1 b

au lieu de Ax = b, avec M qui est la matrice de préconditionnement. En choisissant bien le précondi-

−1

tionneur M , la convergence peut être beaucoup plus rapide.

−1

Une méthode itérative consiste à calculer une suite de vecteurs x telle que
k

lim xk = x
k→+∞

chaque vecteur x est construit à partir des vecteurs précédents x , x , · · · . L'itéré initial x est souvent
pris égal à 0. Pour tester la convergence de la méthode, on ne peut pas calculer l'erreur ||x − x|| car on
k+1 k k−1 0

ne connait pas la solution x. Un mauvais choix serait de considérer ||x − x || car certaines méthodes
k

peuvent stagner ou présenter des plateaux de convergence, donc deux itérés successifs peuvent être assez
k+1 k

proches. En pratique, on calcule le résidu noté r k

rk = b − Axk

et on vérie que la norme de ce résidu est assez petite. On se donne un critère d'arrêt du type
||rk ||
≤ε
||b||
avec ε appelé la tolérance. On notera qu'on divise par la norme de b pour avoir un résidu relatif, ce qui
permet d'avoir une tolérance ε qui ne dépend pas de b. Dans les applications, il arrive souvent d'avoir des
seconds membres b dont la norme est très petite ou très grande. Une valeur typique de ε est 10 , on verra
−6

qu'il est important de prendre ε assez petit pour la plupart des matrices.
On remarquera que
rk = b − Axk = Ax − Axk = A(x − xk )

1
2 CHAPITRE 1. INTRODUCTION

et que la norme
||rk ||2 =< rk , rk >=< A(x − xk ), A(x − xk ) >=< AT A(x − xk ), x − xk >= ||x − xk ||2AT A

La norme du résidu correspond à la norme induite par A A de l'erreur x−x . Comme A A est symétrique
T T

dénie positive, il s'agit bien d'une norme.

k
Chapitre 2

Méthodes de relaxation
Ce sont des méthodes qui s'écrivent toutes la forme suivante
M xk+1 = N xk + b
avec des matrices M et N choisies telles que
A=M −N
et avec M facile à inverser. On remarquera que si la suite (x ) converge vers x, alors on a
k

M x = N x + b ⇒ (M − N )x = b ⇒ Ax = b
Donc lorsque la suite converge, elle converge vers la solution attendue. Les méthodes de relaxation sont
basées sur le découpage suivant de A :

. . . −F 
 
... 
A=  D 
−E
D représente la diagonale de A, −E la partie triangulaire inférieure et −F la partie triangulaire supérieure.

2.1 Méthode de Jacobi

La méthode de Jacobi consiste à écrire le système Ax = b et à indicer la diagonale avec x alors que
les termes extra-diagonaux sont indicés avec x , on obtient le système suivant :
k+1
k

b1 − a1,2 xk2 − a1,3 xk3 − · · · − a1,n xkn

a1,1 xk+1
1 + a1,2 xk2 + a1,3 xk3 + · · · + a1,n xkn = b1 ⇒ xk+1
1 =
a1,1
Xn
bi − ai,j xkj
j=1,j̸=i
ai,1 xk1 + · · · + ai,i−1 xki−1 + ai,i xk+1
i + ai,i+1 xki+1 + · · · + ai,n xkn = b1 ⇒ xk+1
i =
ai,i
La formule générale est donc n
X
bi − ai,j xkj
j=1,j̸=i
xk+1
i =
ai,i
ce qui revient à prendre
M = D, N =E+F

2.2 Méthode de Gauss-Seidel

La méthode de Gauss-Seidel est similaire à celle de Jacobi, car on prend les nouvelles valeurs de x (au
lieu de garder les anciennes valeurs pour Jacobi). On obtient le système suivant
i

b1 − a1,2 xk2 − a1,3 xk3 − · · · − a1,n xkn

a1,1 xk+1
1 + a1,2 xk2 + a1,3 xk3 + · · · a1,n xkn = b1 ⇒ xk+1
1 =
a1,1
i−1
X n
X
k+1
bi − ai,j xj − ai,j xkj
j=1 j=i+1
ai,1 xk+1
1 + · · · + ai,i−1 xk+1 k+1
i−1 + ai,i xi + ai,i+1 xki+1 + · · · ai,n xkn = b1 ⇒ xk+1
i =
ai,i
La formule générale est donc
i−1
X n
X
bi − ai,j xk+1
j − ai,j xkj
j=1 j=i+1
xk+1
i =
ai,i

3
4 CHAPITRE 2. MÉTHODES DE RELAXATION

ce qui revient à prendre

M = D − E, N = F
La méthode de Gauss-Seidel rétrograde consiste à partir de i=n et de remonter jusqu'à i = 1, ce qui revient
à prendre
M = D − F, N =E

2.3 Méthode de relaxation

La relaxation se base sur l'idée que comme x est plus proche de la solution, on peut essayer d'aller
plus loin vers cette nouvelle valeur en posant :
k+1

xk+1
i = (1 − ω)xki + ω x̃ki
où x̃ est la valeur calculée par la méthode de Gauss-Seidel. ω est un paramètre de relaxation. Si ω > 1,
k

on parle de sur-relaxation. On obtient donc la formule suivante

i

i−1
X n
X
bi − ai,j xk+1
j − ai,j xkj
j=1 j=i+1
xk+1
i = (1 − ω)xki + ω
ai,i
Matriciellement, on a donc
xk+1 = (1 − ω)xk + ωD−1 b + ωD−1 Exk+1 + ωD−1 F xk
ce qui donne
(D − ωE)xk+1 = ((1 − ω)D + ωF ) xk + ωb
on divise par ω pour obtenir    
D (1 − ω)
− E xk+1 = D + F xk + b
ω ω
ce qui prouve que
D (1 − ω)
M= − E, N= D+F
On peut aussi utiliser la méthode de relaxation rétrograde basée sur la méthode de Gauss-Seidel rétrograde,
ω ω

ce qui revient à prendre

D (1 − ω)
M= − F, N= D+E
ω ω

2.4 Bilan et implémentation

Les méthodes de relaxation s'écrivent toutes avec la relation de récurrence

M xk+1 = N xk + b
avec le choix suivant pour les matrices M et N :
Méthode M N
Jacobi D E+F
Gauss-Seidel directe D-E F
Gauss-Seidel rétrograde D-F E
Relaxation directe ω1 D - E 1 −ω ω D + F
Relaxation rétrograde ω1 D - F 1 −ω ω D + E
On remarque que −1 −1 −1
M N = M (M − A) = I − M A
En pratique pour implémenter ces méthodes, on pourra utiliser donc l'expression suivante :
xk+1 = xk + M −1 (b − Axk )
ce qui permet de calculer dans un premier temps le résidu
rk = b − Axk
puis de mettre à jour la solution −1
x k+1=x +M r k k

Bien sûr, on ne calcule pas l'inverse de M , on résout le système M z = r . Ce système est facile à résoudre
car M est diagonale ou triangulaire. L'algorithme géneral sécrit donc
k
2.5. CONVERGENCE DES MÉTHODES DE RELAXATION 5
r = b − Ax0
k=0
et
while ||r||/||b|| > ε k ≤ Nmax do
Résoudre M z = rk
Remplacer par
x x+z
Calculer r = b − Ax
end while
Nmax est le nombre maximal d'itérations, c'est utilisé pour éviter une boucle innie. Il peut arriver que
le résidu stagne, il est donc inutile de faire tourner l'algorithme pour rien.
2.5 Convergence des méthodes de relaxation

On rappelle que le rayon spectral d'une matrice A est noté ρ(A) et vaut
Sp
ρ(A) = max
λi ∈ (A)
|λi |

Il est également plus petit que toute norme matricielle :

ρ(A) ≤ ||A||

On a le résultat de convergence suivant :

Theorème 2.5.1 Une méthode de relaxation converge si et seulement si ρ(M −1 N ) < 1
Preuve 1 On note la matrice d'itération B
B = M −1 N

On a
xk+1 = Bxk + M −1 b
En passant à la limite, la solution x vérie

x = Bx + M −1 b

on retranche les deux équations pour obtenir

xk+1 − x = B(xk − x)

En itérant, on obtient
xk − x = B k (x0 − x)
Or la matrice B k converge vers 0 si et seulement si ρ(B) < 1, on obtient donc le résultat annoncé.

Pour la méthode de Jacobi, on a un résultat pour les matrices à diagonalement strictement dominantes.
Ce sont des matrices telles que X
|ai,i | > |ai,j |, ∀i ∈ [1, n]
j̸=i

Theorème 2.5.2 La méthode de Jacobi converge pour des matrices à diagonale strictement dominante
Preuve 2 On utilise la norme innie de la matrice B = M −1 N
n
X
||B||∞ = max |bi,j |
i
j=1

ai,j
Or bi,i = 0 et bi,j = − , on en déduit donc que
ai,i
X |ai,j |
||B||∞ = max
i |ai,i |
j̸=i

Or la matrice est à diagonale srictement dominante, on a donc

X |ai,j |
< 1, ∀i ∈ [1, n]
|ai,i |
j̸=i

On a donc ||B||∞ < 1, ce qui implique que ρ(B) < 1. La méthode de Jacobi converge.
6 CHAPITRE 2. MÉTHODES DE RELAXATION

Pour la méthode de relaxation, la convergence est plus large. On a un résultat pour les matrices symétriques
dénies positives (sdp). Ce sont des matrices qui vérient
< Ax, x̄ > > 0, ∀x ̸= 0
On note ici < x, y > la quantité suivante
n
X
< x, y >= x i yi

qui correspond au produit scalaire lorsque x et y sont deux vecteurs réels. Dans le cas complexe, le produit
i=1

scalaire est < x̄, y >.

Theorème 2.5.3 La méthode de relaxation directe converge pour les matrices A symétriques dénies
positives et pour ω ∈]0, 2[
Preuve 3 On considère une matrice sdp A et un paramètre ω ∈]0, 2[. On montre d'abord que la matrice
M + M T − A est symétrique dénie positive. On a
D D
M= − E, MT = −F
ω ω
d'où
D 2−ω
M + MT − A = 2 − (E + F ) − (D − E − F ) = D
ω ω
On a donc M +M T −A symétrique dénie positive car D est une diagonale strictement positive et ω ∈]0, 2[.
Montrons maintenant que si M + M T − A est symétrique dénie positive, alors la méthode de relaxation
converge. Soit λ ∈ C une valeur propre de M −1 N et x ∈ Cn le vecteur propre associé. On a donc
M −1 N x = λx ⇔ N x = λM x ⇔ (M − A)x = λM x ⇔ (1 − λ)M x = Ax
En faisant le produit scalaire avec x, on obtient donc
< x̄, M x >=
1
1−λ
< x̄, Ax > (2.1)
On prend le conjugué de cette expression pour obtenir :
1
< x, M̄ x̄ >= < x, Āx̄ >
1 − λ̄
On utilise ici la symétrie de A et que A et M sont des matrices réelles pour obtenir
< M T x, x̄ >=
1
1 − λ̄
< Ax, x̄ > (2.2)
On somme les deux équations (2.1) et (2.2), ce qui donne 
T 1 1
< x̄, (M + M − A)x >= + − 1 < x̄, Ax >
1 − λ 1 − λ̄
On remarque que
1 − |λ|2
 
1 1 1 − λ + 1 − λ̄ − (1 − λ)(1 − λ̄) 1 − λλ̄
+ −1 = = =
1 − λ 1 − λ̄ (1 − λ)(1 − λ̄) |1 − λ|2 |1 − λ|2
On a donc
1 − |λ|2
< x̄, (M + M T − A)x >= < Ax, x̄ >
|1 − λ|2
Or < x̄, (M + M T − A)x > et < Ax, x̄ > sont strictement positifs car A et M + M T − A sont sdp. On en
déduit donc que
1 − |λ|2 > 0 ⇔ |λ| < 1
On a donc montré que ρ(M −1 N ) < 1, la méthode de relaxation converge. La méthode de Gauss-Seidel
converge aussi puisque elle revient à prendre ω = 1.

2.6 Choix du paramètre de relaxation

On a obtenu la convergence pour ω ∈]0, 2[, mais est-ce qu'il y a une valeur optimale pour ω ? Pour
certaines matrices, on peut caluler la valeur optimale de ω. Plus précisément si on note
B = D−1 (E + F )
la matrice d'itération de Jacobi, et que la matrice A vérie la propriété suivante
Les valeurs propres de B(α) = αL + Uα ne dépendent pas de α ∈ C
où L et U sont respectivement la partie triangulaire inférieure et supérieure de B, on obtient la valeur
optimale suivante de ω
ωopt =
2
p
1 + 1 − ρ(B) 2

Cette formule peut être utilisée par exemple pour les matrices issues d'une discrétisation par diérences
nies. On remarquera que la valeur optimale de ω tend vers 2 lorsque le pas de maillage tend vers 0.
Chapitre 3

Méthodes de gradient
Dans ce chapitre, on considére des matrices réelles symétriques. <, > correspond donc au produit
scalaire. Les méthodes de gradient cherchent à minimiser la fonctionnelle
1
f (x) = ⟨Ax, x⟩ − < b, x >
2
On calcule le gradient de f en calculant
1 1 1 1
f (x + h) = < Ax, x > + < Ax, h > + < Ah, x > + < Ah, h > − < b, h > − < b, x >
2 2 2 2
avec h ∈ R . En développant, on obtient donc
n

 
1 T 1
f (x + h) = f (x) + (A + A )x − b, h + < Ah, h >
2 2
Le gradient correspond au terme d'ordre un, on a donc
1
∇f (x) = (A + AT )x − b
2
On voit ici qu'il est essentiel que la matrice A soit symétrique pour qu'un point critique de f corresponde
à la solution x̃ de Ax = b. De plus si A est dénie positive, on obtient que
< Ah, h > ≥ 0 ⇒ f (x̃ + h) ≥ f (x̃), ∀h ∈ Rn

La solution de Ax = b correspond donc au minimum de la fonctionnelle f quand la matrice A est sdp.

Dénition 3.0.1 On appelle méthode de descente, une méthode itérative dénie par la relation de récur-
rence
xk+1 = xk + αk dk
où dk ∈ Rn est la direction de descente et αk ∈ R le pas de descente.
Les méthodes vont calculer α pour que x soit plus proche de la solution que x (par exemple
en cherchant α qui minimise f (x )). On remarquera que les méthodes de descente ne sont pas des
k k+1 k

méthodes linéaires, car le pas α va dépendre de manière non-linéaire de x . Pour certaines méthodes
k k+1

comme le gradient conjugué, la direction de descente d dépendra aussi de manière non-linéaire de x .

k k
k k

3.1 Méthode du gradient à pas optimal

Pour cette méthode, on choisit

dk = rk = −∇f (xk ) = b − Axk

la direction de descente est choisie opposée au gradient de f , ce qui semble assez naturel comme choix,
puisque le gradient correspond à la direction de plus forte croissance de f . La direction d est égale au
résidu r qu'on avait vu en introduction.
k

Le pas α est alors choisi pour minimiser f (x ) qui s'écrit

k
k k+1

1
f (xk+1 ) = < A(xk + αk rk ), xk + αk rk > − < b, xk + αk rk >
2
On obtient donc
1 α2
f (xk+1 ) = f (xk ) + αk (< Ark , xk > + < Axk , rk >) − αk < b, rk > + k < Ark , rk >
2 2
soit en utilisant la symétrie de A
αk2
f (xk+1 ) = f (xk ) + αk < Axk − b, rk > + < Ark , rk >
2

7
8 CHAPITRE 3. MÉTHODES DE GRADIENT

soit αk2
f (xk+1 ) = f (xk ) − αk < rk , rk > + < Ark , rk >
2
On remarque que f (x ) est un polynôme de degré 2 avec un coecient en positif car < Ar , r αk2 >> 0
(avec A matrice sdp). Le minimum est donc atteint pour
k+1 k k

< rk , rk >
αk =
< Ark , rk >
On obtient ainsi la méthode du gradient à pas optimal
r = b − Ax0
k=0
et
while ||r||/||b|| > ε k ≤ Nmax
do
Calculer
z = Ar
Calculer
α =< r, r > / < z, r >
Mettre à jour la solution :
x = x + αr
Mettre à jour le résidu
r = r − αz
k =k+1
end while
On a un produit matrice-vecteur par itération. L'algorithme ne peut être utilisé que pour des matrices
sdp. Pour une matrice sdp, on a le résultat de convergence suivant
 k
κ−1
||xk − x||A ≤ ||x0 − x||A
κ+1
ce qui montre que le gradient à pas optimal converge. Ce résultat sera démontré en TD. Le coecient κ
est le conditionnement en norme 2 de A déni par
κ = cond (A) = ||A|| ||A ||
2 2
−1
2

Pour une matrice A symétrique dénite positive, on remarquera que

||A|| = λmax (A) 2

ce qui donne
λmax (A)
λmin (A)
κ=

3.2 Méthode du résidu minimum

La méthode du résidu minimum utilise la même direction de descente que le gradient à pas optimal
dk = rk = b − Axk

En revanche, elle minimise le résidu ||r k+1 || (au lieu de f (x )), ce qui donne l'expression suivante pour
k+1
αk
< Ark , rk >
αk =
< Ark , Ark >
La démonstration de cette formule est faite en TD. On notera que la méthode du résidu converge vers la
solution si la matrice A + A est symétrique dénie positive, la symétrie de A n'est pas nécessaire.
T

3.3 Méthode des directions conjuguées

L'idée de la méthode est de choisir des directions d orthogonales. On choisit alors α pour rendre
l'erreur orthogonale à la direction de descente. Par cette technique, on veut converger en au plus n itéra-
k k

tions, puisqu'au bout de n itérations les directions de descente auront formé une base de R . Si on exige n

l'orthogonalité usuelle, on impose

< ek+1 , dk >= 0
avec e k l'erreur. On impose aussi l'orthogonalité des directions
= xk − x

< d , d >= 0, si i ̸= j
i j

On remarquera qu'aux itérés suivants, on aura l'erreur orthogonale aux directions précédentes car
xk+2 = xk+1 + αk+1 dk+1

on obtient en retranchant la solution x et en faisant le produit scalaire avec d : k

< ek+2 , dk >=< ek+1 , dk > +αk < dk+1 , dk >= 0

3.4. MÉTHODE DU GRADIENT CONJUGUÉ 9
Le problème de cette approche réside dans le calcul du pas α qui vient de la relation
k

0 =< ek+1 , dk >=< ek , dk > +αk < dk , dk >

On en déduit que pour avoir l'orthogonalité, il faut choisir

< ek , dk >
αk = −
< dk , dk >

ce qui est pas pratiquable, car on aurait besoin de connaître l'erreur (et donc la solution) pour calculer
α !
En revanche, si on choisit une A-orthogonalité, c'est à dire
k

< Ad , d >= 0 si i ̸= j
i j

le calcul de α devient
k
0 =< Aek+1 , dk >=< Aek , dk > +αk < Adk , dk >
soit (en utilisant que Ae k = −rk )
< rk , dk >
αk =
< Adk , dk >
Cette fois l'expression est exploitable. Pour naliser la méthode des directions conjuguées, il faut construire
une base de direction d qui soit A-orthogonale. On peut utiliser un procédé de Gram-Schmidt par exemple.
k

3.4 Méthode du gradient conjugué

La méthode du gradient conjugué consiste à construire une base A-orthogonale pour les directions
notées p (au lieu de d ). De plus, on impose que les résidus soient orthogonaux (pour le calcul de α ). Ces
deux conditions susent à déterminer entièrement la méthode. On rappelle que la solution vérie donc
k k k

xk+1 = xk + αk pk

Les directions de descente sont calculées comme

k−1
X
pk = rk − βj pj
j=0

avec
p0 = r0
Pour déterminer les coecients α et β , on impose
k j

< rk+1 , rk >= 0, < Apk , pj >= 0 ∀j ∈ [0, k − 1]

Theorème 3.4.1 Les coecients du gradient conjugué valent

< rk , rk > < rk , rk >
αk = , βk−1 = −
< Apk , pk > < rk−1 , rk−1 >

Les autres coecient βj sont nuls :

βj = 0, ∀j ∈ [0, k − 2]
De plus, on obtient l'orthogonalité des résidus

< rk+1 , rj >= 0, ∀j ∈ [0, k]

Preuve 4 On part de la formule de récurrence des résidus

rk+1 = rk − αk Apk

en faisant le produit scalaire avec rk , on obtient

0 =< rk+1 , rk >=< rk , rk > −αk < Apk , rk >

L'égalité est vériée si on prend

< rk , rk >
αk =
< Apk , rk >
Or on a
k−1
X
rk = pk + βj pj
j=0
10 CHAPITRE 3. MÉTHODES DE GRADIENT

et les directions pj forment une base A-orthogonale, d'où

< Apk , rk >=< Apk , pk >

On en déduit donc que

< rk , rk >
αk =
< Apk , pk >
On a aussi
j−1
X
rj = pj + βm pm
m=0
Pour j < k, on obtient en faisant le produit scalaire avec Apk que
j−1
X
< Apk , rj >=< Apk , pj > + βm < Apk , pm >= 0
m=0

On calcule
< rk+1 , rj >=< rk , rj > −αk < Apk , rj >= 0
car on a supposé les rj orthogonaux jusqu'au rang k (hypothèse de récurrence). On a donc obtenu que les
résidus sont orthogonaux jusqu'au rang k + 1 :
< rk+1 , rj >= 0, ∀j ∈ [0, k]

Traitons maintenant le calcul de βj . On part de la relation

k−1
X
pk = rk − βj pj
j=0

pour obtenir
< Apk , pj >=< Ark , pj > −βj < Apj , pj >
Pour obtenir que les directions soient A-orthogonales, il sut de prendre
< Ark , pj >
βj =
< Apj , pj >

Or, on a
rj+1 − rj = −αj Apj
On en déduit donc que
rj − rj+1
Apj =
αj
On a donc en utilisant l'orthogonalité des résidus :
si j ≤ k − 2


rj − rj+1
  0,
< Ark , pj >= rk , = < rk , rk >
αj  − , si j = k − 1
αk−1

Les coecients βj sont nuls de j = 0 à k − 2, seul βk−1 est non-nul et vaut

< rk , rk > < rk , rk > < Apk−1 , pk−1 >
βk−1 = − =−
< Apk−1 , pk−1 > αk−1 < Apk−1 , pk−1 > < rk−1 , rk−1 >
D'où l'expression nale
< rk , rk >
βk−1 = −
< rk−1 , rk−1 >

On remarque que comme les coecients β sont nuls jusqu'au rang k − 2, les directions p se calculent
avec une récurrence courte
j k

pk = rk + γk pk−1
en ayant posé γ k = −βk−1 . On obtient ainsi l'algorithme du gradient conjugué
x = x0
r = b − Ax0
k=0
p=r
ρ prev= < r, r >
while ||r||/||b|| > εet k ≤ Nmax do
Produit matrice-vecteur : q = Ap
δ = < p, q >
α = ρprev /δ
3.4. MÉTHODE DU GRADIENT CONJUGUÉ 11
x = x + αp
r = r − αq
ρ = < r, r >
γ = ρ/ρ prev
p = γp + r
ρprev =ρ
k =k+1
end while
Par construction, comme les résidus forment une base orthogonale, la méthode du gradient conjugué
convergera en au plus n itérations. Cet algorithme fonctionne pour toute matrice symétrique. Il n'est pas
nécessaire que la matrice soit en plus dénie positive. Néanmoins, si la matrice est symétrique dénie
positive, on est sûr de ne pas rencontrer de division par zéro. De plus on a un résidu monotone, on a le
résultat de convergence suivant (quand A est sdp) :
√ k
κ−1
||xk − x̃||A ≤ 2 √ ||x0 − x̃||A
κ+1

avec x̃ la solution du système Ax = b. On a priori une convergence beaucoup plus rapide que le gradient
à pas optimal, car on a κ au lieu de κ.
√
12 CHAPITRE 3. MÉTHODES DE GRADIENT
Chapitre 4

Préconditionnement
On a vu que le conditionnement de la matrice intervenait dans la vitesse de convergence des méthodes
telles que le gradient conjugué. Le conditionnement joue aussi un rôle sur la précision de la solution. A
cause des erreurs numériques, on ne calcule pas la vraie solution x qui vérie
Ax = b

mais une approximation qu'on note x + ∆x qui vérie le système approché

A(x + ∆x) = b + ∆b

∆x est l'erreur commise sur la solution alors que ∆b correspond au résidu. A partir de ces deux équations,
il est possible de démontrer que
||∆x|| ||∆b||
≤ κ(A)
||x|| ||b||
où κ(A) est le conditionnement de A. La démonstration de cette inégalité sera faite en TD. On voit donc
que même si on obtient un résidu petit, comme il est multiplié par le conditionnement de A, on peut
obtenir une erreur importante sur la solution. On voit la double diculté quand on a une matrice mal
conditionnée (i.e. telle que κ(A) est grand) :
 Convergence plus lente
 Nécessité de prendre une tolérance ε assez petite pour obtenir une solution able
4.1 Préconditionnement à gauche et à droite

Le préconditionnement à gauche consiste à multiplier à gauche par une matrice de préconditionnement

qu'on note M pour obtenir le système :
−1

M −1 A = M −1 b

Le préconditionnement sera bon si le conditionnement de M A est bien plus faible que le conditionnement
−1

de A.
Le préconditionnement à droite change d'inconnue en posant
y = Mx

On obtient alors (
AM −1 y = b
Ax = b ⇔
x = M −1 y
En général, on modie l'algorithme sans préconditionneur pour incorporer le préconditionneur. Dans
l'algorithme avec préconditionneur, on élimine les calculs superus et on utilise le vrai résidu
r = b − Ax

pour la condition d'arrêt de l'algorithme itératif. Par exemple si on part de l'algorithme du résidu minimum
sans préconditionnement :
r = b − Ax0
k=0
et Nmax
while ||r||/||b|| > ε k ≤ do
Calculerw = Ar
Calculerα =< r, w > / < w, w >
Mettre à jour la solution :
x = x + αr
Mettre à jour le résidu
r = r − αz
k =k+1
end while
Il devient avec préconditionnement à gauche
13
14 CHAPITRE 4. PRÉCONDITIONNEMENT

r = b − Ax0
Résoudre Mq = r
k=0
et Nmax
while ||r||/||b|| > ε k ≤ do
Calculerw = Aq
Résoudre Mz = w
Calculerα =< q, z > / < z, z >
Mettre à jour la solution :
x = x + αz
Mettre à jour le vrai résidu
r = r − αw
Mettre à jour le résidu
q = q − αz
k =k+1
end while
On remarque ici qu'on a besoin de mettre à jour le résidu de l'algorithme et le vrai résidu r = b − Ax qui
est utilisé pour tester la convergence de l'algorithme. Dans cet algorithme, on a un produit matrice-vecteur
et un appel au préconditionneur à chaque itération, ce qui est optimal. Dans les TDs, le cas du résidu
minimum avec préconditionnement à droite sera traité.
Il est à noter que les valeurs propres de M A sont les mêmes que celles de AM , par conséquent
−1 −1

la vitesse de convergence devrait être similaire pour les deux cas. Lorsqu'on utilise des algorithmes qui
exigent la symétrie de la matrice comme le gradient conjugué, on ne peut pas l'utiliser directement pour
résoudre M Ax = M b car la matrice M A (ou AM ) n'est en général pas symétrique même si A
−1 −1 −1 −1

et M le sont. Une parade est d'utiliser le produit scalaire < x, y > =< M x, y > pour lequel la matrice
M A est symétrique. On obtient ainsi l'algorithme du gradient conjugué préconditionné qui sera établi
M
−1

dans les TDs.

4.2 Préconditionnement par des méthodes de relaxation

Les préconditionnements sont souvent basés sur une résolution approchée du système linéaire Ax =
b. Les préconditionneurs basés sur les méthodes de relaxation consistent à eectuer une itération de
l'algorithme de relaxation en partant de x = 0.
Par exemple si on veut préconditionner par la méthode de Jacobi, le calcul de
0

z = M −1 r

revient à eectuer une itération de Jacobi sur le système Az = r :

z = z0 + D−1 (r − Az0 )

avec z 0 =0, ce qui donne

z = D−1 r
cela correspond à un préconditionnement diagonal. Appliquer le préconditionneur de Jacobi revient à
diviser le vecteur r par la diagonale de A (terme à terme).
Pour une méthode de relaxation, on a donc
z = M r ⇔ Résoudre
 
−1 D
−E z =r
ω
On a un système triangulaire à résoudre.
Pour des méthodes comme le gradient conjugué, on a besoin que la matrice M soit symétrique, on
peut utiliser par exemple la méthode de Gauss-Seidel symétrique abordée en TD pour obtenir
 Résoudre le système triangulaire (D-E) q = r


z=M r⇔ −1
Multiplier par la diagonale p = D q

Résoudre le système triangulaire (D-F) z = p
4.3 Préconditionnement par la factorisation incomplète

L'idée de la factorisation incomplète est de calculer les facteurs L et U que sur les éléments non-nuls
de A. On interdit ainsi tout remplissage, mais c'est une approximation. Si on note
P = {(i, j) tel que a ̸= 0} i,j

L'algorithme de factorisation incomplète s'écrit

M =A
for i = 1, n − 1 do
et
for j = i + 1, n (j, i) ∈ P do
mj,i
mj,i =
et
mi,i
for k = i + 1, n (j, k) ∈ P do
4.4. PRÉCONDITIONNEMENT PAR DES MÉTHODES MULTIGRILLES 15
mj,k = mj,k − mj,i mi,k
end for
end for
end for
La matrice M contient ici les facteurs incomplets L (partie triangulaire inférieure) et U (partie triangulaire
supérieure).
L'application du préconditionneur revient á résoudre le système avec la factorisation LU ainsi calculée :

Résoudre Ly = r
Mz = r ⇔
Résoudre U x = y
De la même manière, on peut calculer la factorisation de Cholesky incomplète (pour une matrice A
sdp) :
M =A
for i = 1,v
n do
u
u i−1
X
mi,i = tmi,i − m2 i,k

for j = i + 1, n et (j, i) ∈ P do
k=1

i−1
X
aj,i − mi,k mj,k
k=1
mj,i =
mi,i
end for
end for
Pour obtenir, une meilleure convergence, on peut élargir le prol P , on obtient ainsi la méthode ILU(1)
en prenant
P = {(i, j) tel que ∃k tel que a ̸= 0 et a ̸= 0}
i,k k,j

et la méthode ILU(k), où k est le nombre maximal d'intermédiaires entre i et j. Une autre approche
est de ne pas conserver les valeurs petites de L, on obtient ainsi la méthode ILUT (Incomplete LU with
Threshold) paramétrée par un seuil ε.
4.4 Préconditionnement par des méthodes multigrilles

L'idée des méthodes multigrille est de calculer la solution sur une grille plus grossière an de construire
une solution approchée. Lorsqu'on décompose la solution sur les modes, la résolution sur une grille grossière
permet d'atténuer l'erreur pour les modes basses fréquences, mais les modes hautes fréquences sont peu
ou pas atténués. Il est donc nécessaire d'adjoindre à ces algorithmes des étapes de lissage qui permettront
de réduire l'erreur haute fréquence. Si on considère par exemple la résolution d'un laplacien
−∆u = f

en 1-D, sur le segment [0, 1]. Une discrétisation par diérences nies donne la matrice
 
2 −1 0

... ... ...

 −1 2 −1 
1  
A∆x =
 
∆x2 
 

 −1 2 −1 
0 −1 2

avec A une matrice de taille n et 1

∆x =
n+1
on obtient le système linéaire
A∆x U = F
A∆x/2 désigne la matrice sur la grille ne et A la matrice sur la grille grossière. La résolution de ce
système par une méthode bigrille s'écrit comme suit :
∆x

x = x0
r = b − Ax0
k=0
while ||r||/||b|| > ε k ≤ et do Nmax
Pré-lissage
x = S1 (x, b)
Calcul du résidur = b − A∆x/2 x
Restriction sur grille grossière
r̃ = Rr
Résolution de A∆x ỹ = r̃
Prolongement sur grille ne
y = P ỹ
16 CHAPITRE 4. PRÉCONDITIONNEMENT

Mise à jour solution x = x + y

Post-lissage x = S (x, b)
Mise à jour résidu r = b − A
2
∆x/2 x
k =k+1
end while
où S et S sont des lisseurs, R un opérateur de restriction, P un opérateur de prolongement. Pour
l'opérateur de prolongement, on prend une interpolation linéaire :
1 2

(P u)2i = ui
ui + ui+1
(P u)2i+1 =
2
avec la convention u 0 = un+1 = 0 (condition de Dirichlet). Matriciellement, on a donc
 
1 0 ···
1  2 0 ··· 
P =  
2 1 1 ··· 
0 2 0 ···

Pour l'opérateur de restriction, on choisit 1

R = PT
2
ce qui revient à prendre v2i−1 + 2v2i + v2i+1
(Rv)i =
4
L'utilisation de la transposée de P permet d'obtenir un préconditionneur symétrique (si on choisit aussi
S = S ), ce qui est important si on utilise un algorithme comme le gradient conjugué. Pour le prélisseur,
T

on pourra par exemple considérer la méthode de Jacobi pondérée

2 1

S1 (x, b) = x + ωD−1 (b − A∆x/2 x)

où ω est le paramètre de relaxation. Si on n'a pas de nécessité d'un préconditionneur symétrique, on

pourra prendre S (x, b) = x (pas de post-lissage). Si on regarde l'action du lisseur, on trouve que les
valeurs propres de I − ωD A sont égales à
2
−1

 
2 kπ
λk = 1 − 2ω sin , k = 1..2n + 1
2(2n + 2)
et le mode associé vaut 
j kπ

vj = sin , j = 1 · · · 2n + 1
2n + 2
Les modes basses fréquences correspondent à j petit ou proche de 2n + 1, alors que les modes hautes
fréquences correspondent à j proche de n. Si on prend ω = 2 , la valeur propre vaut
1

 
2 kπ
λk = cos
2(2n + 2)

Les modes basses fréquences sont peu atténués (cos proche de 1), alors que les modes hautes fréquences
2

sont fortement atténuées. C'est ce qu'on attend d'un lisseur.

L'algorithme bigrille est une méthode linéaire qui peut s'écrire
M xk+1 = N xk + b

même si en pratique, on n'implémente pas l'algorithme ainsi. D'un point de vue théorique, il est quand
même intéressant d'utiliser cette forme pour étudier la convergence. Lorsqu'on n'a pas de post-lissage
(S = I ), on obtient
2
M −1 N = (I − P A−1
∆x RA∆x/2 )J

où J est la matrice d'itération du lisseur (I −ωD pour Jacobi pondéré). Lorsqu'on calcule le rayon
−1 A

spectral ρ(M 1N ), on trouve

∆x/2
−

ρ(M N ) = , pour ω =
1 −1 1
2 2
ce qui est résultat remarquable, car on obtient une vitesse de convergence indépendante du pas de maillage.
La valeur optimale est obtenue pour ω = 3 : 2

pour ω = 23
1
ρ(M −1 N ) = ,
3
L'algorithme multigrille consiste à remplacer la résolution directe par un appel récursif à l'algorithme
multigrille sur la grille grossière :
4.4. PRÉCONDITIONNEMENT PAR DES MÉTHODES MULTIGRILLES 17
Multigrille(∆x, b, x0 , γ )
x = x0
Pré-lissage x = S (x, b)
Calcul du résidu r = b − A x
1

Restriction sur grille grossière r̃ = Rr

∆x/2

if ∆x assez grand then

Résolution de A ỹ = r̃
∆x
else
ỹ = 0
for g = 1, γ do
ỹ = Multigrille (2 ∆x, r̃, ỹ, γ)
end for
end if
Prolongement sur grille ne y = P ỹ
Mise à jour solution x = x + y
Post-lissage x = S (x, b)
Retourner x
2

La méthode multigrille s'écrira alors

x = x0
r = b − Ax0
k=0
while ||r||/||b|| > ε et k ≤ Nmax do
Multigrille( , b, x, γ)
∆x/2
Mise à jour résidu r = b − A x ∆x/2
k =k+1
end while
Le paramètre γ règle le nombre d'itérations à faire pour chaque niveau. En général, on choisit γ = 1
(V-cycle) ou γ = 2 (W-cycle). On peut aussi prendre un γ variable, égal à 2 pour la première itération
(g = 1) et à 1 (g = 2) pour la seconde (F-cycle).
Pour obtenir une convergence encore plus rapide, la méthode multigrille est utilisée comme précondi-
tionneur à des algorithmes itératifs (comme le gradient conjugué). Appliquer le préconditionneur revient
à eectuer une itération multigrille (en partant de z = 0) :
0

Résoudre M z = r ⇔ z = Multigrille(∆ /2, r, 0, γ)

x
18 CHAPITRE 4. PRÉCONDITIONNEMENT
Chapitre 5

Factorisation QR et équations normales

Dans cette section, on s'intéresse à la factorisation QR d'une matrice qui est une étape du GMRes.
Une matrice A peut s'écrire sous la forme
A = QR

avec Q matrice orthogonale et R une matrice triangulaire.

5.1 Transformation de Householder

La transformation de Householder est une rotation qui permet de transformer un vecteur X en un

vecteur Y de même norme. Elle équivaut à une réexion par rapport au plan médian entre les vecteurs X
et Y . On note X −Y
U=
||X − Y ||
la normale à ce plan. La matrice de Householder H est alors dénie comme
H = I − 2U T U

Par construction, on a
HX = Y

H est une matrice orthogonale.

Remarque 5.1.1 Pour appliquer la transformation de Householder à un vecteur V , il est plus ecace de
ne pas calculer la matrice H, mais de stocker uniquement le vecteur U , et de faire le calcul

HV = V − 2(U · V )U

De la même manière, la multiplication à droite de H par une matrice A se fait en calculant HCi pour
chaque colonne Ci de A. La multiplication à gauche revient à calculer HLi pour chaque ligne Li .

Exemple 5.1.1 Calculer la matrice de Householder qui transforme X = (3, 4, 0) en Y = (0, 0, 5). Utiliser
l'algorithme rapide pour calculer HV pour V = (5, 10, −5).

On prend X = (3, 4, 0) , Y = (0, 0, 5) , on a

X −Y 1
U= = √ (3, 4, −5)
||X − Y || 50

On en déduit que
   
32 −24 30 0.64 −0.48 0.6
1 
H = I − 2U T U = −24 18 40  =  −0.48 0.36 0.8 
50
30 40 0 0.6 0.8 0

Pour calculer HV avec V = (5, 10, −5), on calcule

80
U ·V = √
50

On en déduit que      
5 3 −4.6
160  4  =  −2.8 
HV =  10  −
50
−5 −5 11

19
20 CHAPITRE 5. FACTORISATION QR ET ÉQUATIONS NORMALES

5.2 Factorisation QR en utilisant les transformations de Householder

On va multiplier successivement par des matrices de Householder pour former une matrice triangulaire :
Hn Hn−1 · · · H1 A = R

avec qui sera donc triangulaire, et H la transformation de Householder qui transforme la colonne
R
en Y = (α, 0, · · · , 0). On a donc
i
X = (R(i, i), R(i + 1, i), · · · , R(m, i))

Q = H1 H2 · · · Hn

L'algorithme de factorisation QR s'écrit comme suit :

Q = I; R = A
for i = 1, n do
;
X = R(i :, i) α = ±||X|| ; Y = (α, 0, · · · , 0)
Calculer la transfo de Householder H qui transforme X en Y
Appliquer H à la sous-matrice R(i :,i :)
Remplacer Q( :,i :) par Q( :,i :) H
end for
Pour une plus grande stabilité, on pourra choisir
α = −signe(X )||X|| 1

5.3 Méthode de Givens

C'est une méthode alternative pour calculer la factorisation QR d'une matrice. On élimine chaque
coecient R en appliquant une rotation de la forme
j,i

 
c 0 ··· −s 0 · · ·

.. . . . . . . ..
 0 1
 0 ··· 0 

Q1 = 
 

 
 s 0 ··· c 0 
0 ··· ··· ··· 1

à la matrice R. c et s sont choisis pour éliminer le coecient R : j,i

ri,i rj,i
c= q , s = −q
2 + r2
ri,i 2 + r2
ri,i
j,i j,i

En multipliant par les rotations successives Q , Q , · · · Q , on obtient ainsi une matrice triangulaire R.
On a donc l'égalité
1 2 N

QN QN −1 · · · Q2 Q1 A = R

On en déduit que
Q = (QN QN −1 · · · Q2 Q1 )T = QT1 QT2 · · · QTN −1 QTN

Pour former la matrice Q, il faut donc multiplier par la transposée de Q à droite. On obtient ainsi
l'algorithme suivant pour la factorisation QR avec la méthode de Givens
i

Q = I; R = A
for i = 1, n do
for j = i + 1, m do
if R(j, i)
p̸= 0 then
D = R(i, i)2 + R(j, i)2
c = R(i, i)/D; s = −R(j, i)/D
L1 = R(i, :); L2 = R(j, :)
R(i, :) = c ∗ L1 − s ∗ L2
R(j, :) = s ∗ L1 + c ∗ L2
C1 = Q(:, i); C2 = Q(:, j)
Q(:, i) = c ∗ C1 − s ∗ C2
Q(:, j) = s ∗ C1 + c ∗ C2
end if
end for
end for
L'avantage de cet algorithme est que le traitement d'une matrice de Hessenberg est immédiat.
5.4. PROBLÈME AUX MOINDRES CARRÉS 21
5.4 Problème aux moindres carrés

Une application de la factorisation QR se rencontre lorsqu'on cherche à minimiser

m
X n
X
||Ax − b||22 = (bj − ai,j xj )2
j=1 i=1

lorsque la matrice A est rectangulaire. Ce problème est connu comme problème aux moindres carrés. On
a n inconnues (les composantes x ) et m équations. On prend en général m > n, de telle sorte qu'on ait
plus d'équations que d'inconnues, on a donc en général
j

min ||Ax − b||2 ̸= 0

x∈Rn

On utilise la factorisation A = QR pour obtenir

||Ax − b|| = ||QRx − b||2 = ||QT (QRx − b)||2 = ||Rx − QT b||2

On peut écrire la matrice R et le vecteur Q b sous la forme

T

   
T T z1
R= , Q b=
0 z2

où T est une matrice triangulaire carrée, z les n premières composantes de Q b. On a donc

1
T

||Rx − QT b||22 = ||T x − z1 ||22 + ||z2 ||2

Si la matrice T est inversible, on peut annuler le premier terme en prenant

x = T −1 z1

et on obtient alors
min ||Ax − b||2 = ||z2 ||
x∈Rn

5.5 Équations normales

Pour résoudre un système linéaire non-symétrique

Ax = b

une solution est de multiplier par A pour obtenir

T

AT Ax = AT b

On a alors un système symétrique. Ce procédé peut être aussi utilisé pour la minimisation aux moindres
carrés de la section précédente. On peut utiliser la méthode du gradient conjugué pour caluler la solution
de ce système symétrique : c'est la méthode CGNE (Conjugate Gradient Normal Equation). Il est à noter
qu'on a
Cond (A A) = Cond (A)
2
T
2
2

car
||A||22 = ρ(AT A)
Par conséquent, la matrice A A est particulièrement mal conditionnée. La méthode CGNE converge mais
T

beaucoup plus lentement que d'autres méthodes abordés dans le chapitre suivant.
22 CHAPITRE 5. FACTORISATION QR ET ÉQUATIONS NORMALES
Chapitre 6

Méthodes de Krylov
Dans cette section, on s'intéresse aux méthodes de Krylov utilisées pour résoudre des systèmes linéaires
non-symétriques. Ces méthodes sont basées sur les espaces de Krylov
K = Vect(r , Ar , A r , · · · , A
m 0 0
2
0 r ) m−1
0

En général, on cherche à projeter la solution dans cet espace de Krylov.

6.1 Base orthogonale

An de projeter dans l'espace de Krylov K , une approche consiste à construire une base orthonormée
de cet espace par le procédé de Gram-Schmidt dont la version initiale s'écrit :
m
(vi )
v1 = r0 /||r0 ||
for j = 1, m do
w = Avj
for i = 1, m do
hi,j =< w, vi >
end for
w=w− m
P
i=1 hi,j vi
vj+1 = w/||w||
end for
On note V la matrice dont les colonnes sont les vecteurs orthonormés v . Cette matrice est donc une
matrice orthogonale. L'algorithme de Gram-Schmidt initial est modié an de limiter les erreurs d'arrondi
m i

(pour que la base soit plus orthogonale numériquement) :

v1 = r0 /||r0 ||
for j = 1, m do
w = Avj
for i = 1, m do
hi,j =< w, vi >
w = w − hi,j vi
end for
hj+1,j = ||w||
vj+1 = w/hj+1,j
end for
Dans la matrice H formée des coecients h , on a mis la norme de w dans la sous-diagonale (coecient
). On obtient ainsi la relation suivante :
m i,j
hj+1,j

j
X
Avj = hj+1,j vj+1 + hi,j vj
i=1

Matriciellement, on obtient donc

AVm = Vm+1 Hm
On notera que H est une matrice de Hessenberg.
m

6.2 Méthode Gmres

La méthode Gmres (Generalized Minimum Residual) consiste à minimiser le résidu ||b − Ax ||, en
cherchant x de la forme :
k
k
xk = x0 + Vm y
avec y ∈ R . Par construction, on a bien V
m
my ∈ Km . De plus,
b − Axk = b − Ax0 − AVm y = r0 − AVm y

23
24 CHAPITRE 6. MÉTHODES DE KRYLOV

On utilise la relation r 0 = βv1 , avec β = ||r ||, pour obtenir

0

b − Axk = βv1 − AVm y = Vm+1 (βe1 − Hm y)

Comme la matrice V m+1 est orthogonale, on obtient que

||b − Axk || = ||βe1 − Hm y||

Ici, pour miniser la norme de βe − H y, on va utiliser la factorisation QR de la matrice H , qui est

une matrice de petite taille. Le vecteur y est donc le vecteur qui minimise ||βe − H y||. On obtient ainsi
1 m m

l'algorithme suivant pour le Gmres

1 m

r0 = b − Ax0
β = ||r0 ||
v1 = r0 /β
Hm = 0
for j=1, m do
w = Avj
for i=1, j do
hi,j =< w, vi >
w = w − hi,j vi
end for
hj+1,j = ||w||
vj+1 = w/hj+1,j
end for
Calculer y qui minimise ||βe − H y||
Former x = x + V y
1 m
0 m

On va utiliser la méthode de Givens an de factoriser H progressivement. Comme la matrice est de

Hessenberg, le second membre Q βe ne sera modié que sur les lignes j et j+1 à l'itération j. le résidu
m
T

vaudra alors |(Q βe ) |, ce qui nous permettra d'arrêter l'algorithme en testant cette valeur. On obtient
1
T

ainsi l'algorithme du Gmres avec condition d'arrêt :

1 j+1

r0 = b − Ax0
β = ||r0 ||
v1 = r0 /||r0 ||
Hm = 0
j=1
z = (β, 0, · · · , 0)
et
while j ≤ m |z(j)|/||b|| > ε do
w = Avj
for i=1, j do
hi,j =< w, vi >
w = w − hi,j vi
end for
hj+1,j = ||w||
vj+1 = w/hj+1,j
for k=1,j do
L1 = hk,j , L2 = hk+1,j
hk,j = ck L1 − sk L2
hk+1,j = sk L1 + ck L2
end qfor
α = h2j,j + h2j+1,j
cj = hj,j /α, sj = −hj+1,j /α
hj,j = α; hj+1,j = 0
zj+1 = zj sj
zj = zj cj
j =j+1
end while
Résoudre le système triangulaire H(1 : j − 1, 1 : j − 1)y = z(1 : j − 1)
Former x = x + X y v
j−1
0 i i
i=1

Pour limiter l'espace mémoire requis par Gmres, l'algorithme est répété plusieurs fois, avec une valeur xe
de m.
6.3. ALGORITHME DU GRADIENT BICONJUGUÉ 25
6.3 Algorithme du gradient biconjugué

L'algorithme du gradient biconjugué utilise la biorthogonalité pour construire des directions p et p̃

A-biorthogonales (biorthogonales vis-à-vis de A) et des résidus r et r˜ biorthogonaux :
m m
m m

< p̃ , Ap >= 0, si i ̸= j
i j

< r̃ , r >= 0, si i ̸= j
i j

On résout quelque part deux systèmes linéaires Ax = b et A x̃ = b̃ simultanément. Comme on n'a pas
T

besoin de la solution x̃, celle-ci n'est pas calculée. Similairement au gradient conjugué, on a
xm+1 = xm + αm pm

et pour les résidus

rm+1 = rm − αm Apm
r̃m+1 = r̃m − αm AT p̃m
Les directions de descente vérient la relation
m−1
X
pm = rm + βj pj
j=0

m−1
X
p̃m = r̃m + βj p̃j
j=0

Les coecients β sont choisis pour assurer la A-biorthogonalité des p :

j m

< r̃m , Apj >

< p̃m , Apj >= 0 ⇒ βj = −
< p̃j , Apj >
Or on a
rm+1 − rm = −αm Apm
on en déduit que  
rj − rj+1
< r̃m , Apj >= r̃m ,
αj
Comme on a la biorthogonalité des résidus (i.e. < r̃ , r >= 0 si i ̸= j), on obtient que
i j

β = 0, si j < m − 1
j

On a bien une récurrence courte pour calculer les vecteurs p et p̃ . On obtient que
m m

< r̃m , rm >

βm−1 =
αm−1 < p̃m−1 , Apm−1 >
Il nous reste à déterminer α pour écrire l'algorithme nal. La valeur de α est obtenue en écrivant la
biorthogonalité des résidus
m m

< rm+1 , r̃m >= 0 =< rm − αm Apm , r̃m >

On obtient donc < rm , r̃m >
αm =
< Apm , r̃m >
On remarque que la relation < r̃ m+1 , rm >= 0 fournit l'expression
< rm , r̃m >
αm =
< AT p̃m , rm >
Ces deux expressions sont en fait égales car on a
< AT p̃m , rm >=< Apm , r̃m >

cela peut être démontré par récurrence car on a

r0 = r̃0 = p0 = p̃0

et
< AT p̃m , rm >= βm−1 < AT p̃m−1 , rm > + < AT r̃m , rm >
< Apm , r̃m >= βm−1 < Apm−1 , r̃m > + < Arm , r̃m >
En injectant l'expression de α m−1 dans la formule pour β , on obtient m−1

< r̃m , rm >

βm−1 =
< r̃m−1 , rm−1 >
26 CHAPITRE 6. MÉTHODES DE KRYLOV

Pour α , on peut utiliser la A-biorthogonalité des p pour remarquer

m i

< Apm , r̃m >=< Apm , p̃m − βm−1 p̃m−1 >=< Apm , p̃m >

et obtenir < rm , r̃m >

αm =
< Apm , p̃m >
On en déduit l'algorithme BiCg (Biorthogonal Conjugate Gradient) :
r = b − Ax0
r̃ = r
p = r; p̃ = r
k=0
et
while ||r||/||b|| > ε k ≤ N max do
ρ0 =< r, r̃ >
q = Ap
q̃ = AT p̃
α =< r, r̃ > / < q, p̃ >
x = x + αp
r = r − αq
r̃ = r̃ − αq̃
ρ =< r, r̃ >
β = ρ/ρ0
p = r + βp
p̃ = r̃ + β p̃
k =k+1
end while
Il est à noter que cet algorithme coïncide avec celui du gradient conjugué lorsque A est symétrique.