Analyse Numer I Que I

Analyse Numérique
Radhia Bessi & Maher Moakher
Version 2020–2021
Table des matières
1 Méthodes directes pour la résolution des systèmes linéaires 1

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Principe des méthodes directes . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 Rappels et notations sur les vecteurs et les matrices . . . . . . . . . 2
1.2.2 Résolution d’un système triangulaire . . . . . . . . . . . . . . . . . 4
1.3 Factorisation LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3.1 Rappel sur la méthode de Gauss . . . . . . . . . . . . . . . . . . . . 5
1.3.2 Factorisation LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Factorisation de Cholesky des matrices symétriques définies positives . . . 12
1.4.1 Rappels sur les matrices symétriques . . . . . . . . . . . . . . . . . 12
1.4.2 Factorisation des matrices symétriques . . . . . . . . . . . . . . . . 14
1.4.3 Factorisation de Cholesky . . . . . . . . . . . . . . . . . . . . . . . 15
2 Méthodes itératives pour la résolution des systèmes linéaires 19

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2 Rappels sur les normes matricielles . . . . . . . . . . . . . . . . . . . . . . 19
2.2.1 Normes matricielles subordonnées . . . . . . . . . . . . . . . . . . . 20
2.3 Méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.3.1 Méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.3.2 Méthodes de relaxation . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.3.3 Méthode de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . 29
2.4 Conditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3 Optimisation sans contraintes 37

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2 Optimisation sur Rn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2.1 Conditions d’optimalité . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.2 Problème d’optimisation quadratique . . . . . . . . . . . . . . . . . 43
3.2.3 Problème de moindre carré . . . . . . . . . . . . . . . . . . . . . . . 43
3.3 Algorithmes de descentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.4 Méthodes des gradients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.4.1 Méthode du gradient à pas fixe : . . . . . . . . . . . . . . . . . . . . 46
3.4.2 Méthodes du gradient à pas optimal . . . . . . . . . . . . . . . . . 46
i
TABLE DES MATIÈRES
3.4.3 Méthodes du gradient conjugué . . . . . . . . . . . . . . . . . . . . 49

3.4.4 Vitesse de convergence . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.4.5 Méthodes des gradients et préconditionnement . . . . . . . . . . . . 53
4 Optimisation avec contraintes linéaires 55

4.1 Problèmes d’optimisations sous contraintes . . . . . . . . . . . . . . . . . . 55
4.1.1 Existence et unicité de minimum . . . . . . . . . . . . . . . . . . . 56
4.1.2 Condition d’optimalité . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1.3 Conditions d’optimalités pour les contraintes d’inégalité linéaires . . 59
4.1.4 Cas de contraintes d’égalités et d’inégalités linéaires . . . . . . . . . 60
4.1.5 Problème quadratique . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2 Quelques algorithmes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2.1 Méthode du gradient projeté . . . . . . . . . . . . . . . . . . . . . . 63
4.2.2 Méthode d’Uzawa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
ii
Chapitre 1
Méthodes directes pour la résolution

des systèmes linéaires
1.1 Introduction
La résolution des systèmes linéaires de grandes tailles est l’un des plus importants
problèmes en analyse numérique. Le but de ce chapitre est d’étudier des méthodes de
résolution numérique d’un linéaire Ax = b, où A est une matrice carrée inversible.
Pour motivation, commençons par citer le problème mécanique classique suivant qui
conduit à la résolution d’un système linéaire.
La déformation x d’une corde élastique, fixée aux bords et soumise à un champ de
force f , peut se traduire par l’équation différentielle suivante
−x00 (t) = f (t),

t ∈ [0, 1]
(1.1)
x(0) = x(1) = 0,
pour f une fonction continue sur [0, 1]. En général, il n’est pas possible d’expliciter la
solution exacte de ce problème. L’idée donc est de chercher une solution approchée xh de
x en prenant une subdivision 0 = t0 ≤ t1 ≤ · · · ≤ tn+1 = 1, avec ti = ih, i = 0, . . . , n + 1
1
et h = . On se limitera à calculer xh (ti ) = xi ' x(ti ), pour i = 0, . . . , n + 1 et par
n+1
interpolation par exemple, on peut avoir xh sur tout l’intervalle [0, 1].
Si on suppose que notre solution x est de classe C 2 , alors on a les deux développement
de Taylor suivants :
h2
x(ti+1 ) = x(ti + h) = x(ti ) + x0 (ti )h + x00 (ti ) + O(h3 ),
2
et
h2
x(ti−1 ) = x(ti − h) = x(ti ) − x0 (ti )h + x00 (ti ) + O(h3 ).
2
Si on fait la somme de deux égalités on obtient
x(ti + h) − 2x(ti ) + x(ti − h)
x00 (ti ) = + O(h).
h2
1
Chapitre 1. Méthodes directes pour la résolution des systèmes linéaires
Si on néglige les termes d’ordre O(h) dans l’expression de la dérivée seconde, le problème
discrétisé pour résoudre notre équation devient
( x − 2x + x
i+1 i i−1
− = f (xi ), i = 1, . . . , n
h2 (1.2)
x0 = xn+1 = 0
qui est équivalent au système linéaire Ah x = bh , où

 
2 −1 0 ··· 0
..     
−1 2 −1 .  x1 f (t1 )

1  .. .. ..  .  . 
Ah = 2 
h  0 . . . 0  x =  ..  et bh =  ..  .
 . xn f (tn )
 ..

−1 2 −1
0 ··· 0 −1 2
La solution approchée est d’autant plus proche de la solution exacte x lorsque n est grand.
On rappele que la solution unique x = (x1 , . . . , xn ) d’un système Ax = b, pour
A = (aij )1≤i,j≤n une matrice inversible et b = (b1 , . . . , bn )T est donnée par les formules de
Cramer suivantes :
det Ai
xi = , i = 1, . . . , n,
det A
où det désigne le déterminant et Ai est la matrice d’ordre n obtenue à partir de A en
remplaçant sa colonne i par le vecteur de second membre b. Donc la résolution d’un
système linéaire d’ordre n, par les formules de Cramer fait intervenir le calcul de n +
1 déterminants dont chaque déterminant nécessite de l’ordre de n! multiplications. La
méthode de Cramer devient trop coûteuse même pour des matrices de tailles assez petites.
D’où l’idée de concevoir des méthodes qui en général donnent la solution à une valeur près,
mais avec un nombre raisonnable d’opérations.
On distinguera deux méthodes numériques pour résoudre un système linéaire. Les
méthodes directes, dont certaines font l’objet de ce chapitre et les méthodes itératives qui
seront développées dans les chapitres suivants.
1.2 Principe des méthodes directes

On appelle méthode directe pour résoudre un système de type Ax = b, une méthode
qui donne x après un nombre fini d’opérations élémentaires. Le principe de ces méthodes
est de se ramener à un système linéaire équivalent, mais qui est plus simple à résoudre.
C’est le cas par exemple où la matrice du système est diagonale, triangulaire ou orthogo-
nale.
1.2.1 Rappels et notations sur les vecteurs et les matrices

Tout au long de ce cours on utilisera les notations suivantes :
2
1.2 Principe des méthodes directes
n
Pour n un entier naturel non  nul fixé,
 on note par R l’espace vectoriel des vecteurs
x1
 .. 
colonnes à coefficients réels x =  . . Le vecteur ligne xT = (x1 , . . . , xn ) désigne la
xn
transposée du vecteur x ∈ Rn .
Le produit scalaire dans Rn sera noté (·, ·) et est défini par
n
X
(x, y) = xi y i ,
i=1
et la norme associée est la norme euclidienne donnée par
n
!1
X 2 1
kxk2 = |xi |2 = (x, x) 2 .
i=1
On notera l’espace vectoriel des matrices à n lignes et p colonnes et à coefficients dans

R par Mn,p (R) et par Mn (R) l’ensemble des matrices carrées d’ordre n à coefficients dans
R.
La transposée d’une matrice A = (aij )1≤i≤n,1≤j≤p ∈ Mn,p (R) par
AT = (aji )1≤i≤p,1≤j≤n ∈ Mp,n (R).
La matrice transposée vérifie
(Ax, y) = (Ax)T y = xT AT y = (x, AT y) pour tout x ∈ Rp et y ∈ Rn .
Définition 1.2.1 Soit A = (aij )1≤i,j≤n ∈ Mn (R). A est dite :

• diagonale si aij = 0 pour tout i 6= j.
• triangulaire inférieure si aij = 0 pour tout i < j.
• triangulaire supérieure si aij = 0 pour tout i > j.
• inversible si son déterminant est non nul ou s’il existe une matrice B ∈ Mn (R)
telle que
AB = BA = In
où In est la matrice unité d’ordre n donnée par
 
1
In = 
 .. .

.
1
La matrice B est appelée inverse de A et sera notée A−1 .

• semblable à une matrice D ∈ Mn (R) s’il existe une matrice inversible P telle que
D = P −1 AP .
3
• diagonalisable si elle est semblable à une matrice diagonale.
On peut vérifier facilement que

Proposition 1.2.1
1. Le produit de deux matrices triangulaires supérieures (repectivement inférieures)
est une matrice triangulaire supérieure (respectivement inférieure).
2. Une matrice A = (aij )1≤,i,j≤n ∈ Mn (R) triangulaire est inversible si et seulement
si aii 6= 0, pour tout i = 1, . . . , n. De plus A−1 est de même type que A et
1
(A−1 )ii = , pour tout i = 1, . . . , n.
aii
1.2.2 Résolution d’un système triangulaire

Soit le système trianglaire supérieur


 u11 x1 + u12 x2 + · · · u1n xn = b1
u22 x2 + · · · u2n xn = b2


... .. . (1.3)


 .
 unn xn = bn
Si U est la matrice triangulaire supérieure suivante :

 
u11 u12 · · · u1n
 u22 · · · u2n 
U = ..  ,
 
. . .
 . 
unn
alors, le système (1.3) est équivalent à U x = b où b = (bi )1≤i≤n . On utilise la méthode de
remontée pour résoudre ce système dont l’algorithme est le suivant :

xn = bn /unn
 Pour k = (n − 1), . . . , 1
 n
!
 X
 xk = b k − uki xi /ukk

 i=k+1
Fin de la boucle sur k
A chaque étape k on a une division, (n − k) multiplication pour 1 ≤ k ≤ n et (n − k)

additions pour k ≤ n − 1. Au total le nombre d’opérations est
n n−1
X X n(n + 1) n(n − 1)
(n − k + 1) + (n − k) = + = n2 .
k=1 k=1
2 2
4
1.3 Factorisation LU
Pour la résolution d’un système triangulaire inférieur Lx = b, on utilise la méthode de

descente dont l’algorithme est le suivant :

b1
 x1 = l11

 Pour k = 2, . . . , n
 k−1
!

 x = b − X 1
k k lki xi
lkk

i=1

Fin de la boucle sur k
Cette méthode nécessite aussi n2 opérations.
1.3.1 Rappel sur la méthode de Gauss
Parmi les méthodes directes classiques pour résoudre un système linéaire Ax = b, pour
A une matrice carrée inversible d’ordre n, on cite la méthode d’élimination de Gauss dont
le principe est d’effectuer uu nombre fini d’opérations algébriques linéaires sur les lignes
de la matrice A et sur b, pour se ramener à un système triangulaire supérieur équivalent.
Exemple 1.3.1 Soit à résoudre le système linéaire d’ordre 4, suivant :



 x1 − 3x2 − x3 = 2
−x1 + x2 + 2x4 = 3


 x2 − x3 = 1
2x1 + x2 − x4 = 0

Sous forme matricielle ce système s’écrit sous la forme Ax = b, pour

   
1 −3 −1 0 2
−1 1 0 2  3
A= et b =  .
0 1 −1 0   1
2 1 0 −1 0
Dans la première étape on élimine la première inconnue x1 des équations 2,3 et 4 en

combinant chacune avec la première équation. Afin d’éliminer il faut d’abord vérifier que
x1 apparait dans la première équation. Si ce n’est pas le cas, il faut permuter l’équation
avec une autre dont le coefficient de x1 est non nul. On choisit comme premier pivot α1
le coefficient de x1 dans la nouvelle première équation qui est appelée ligne de pivot. Dans
notre exemple α1 = 1. Eliminer x1 des autre équations revient à annuler les coefficients
de la première colonne de A en dessous de la diagonale. Ceci revient dans la méthode de
5
Gauss à remplacer chaque ligne Li de A et de b par Li − aαi11 L1 , pour i = 2, 3 et 4. Dans

notre cas
L2 ←− L2 + L1 , L3 ←− L3 , L4 ←− L4 − 2L1 .
Après cette première étape, le système équivalent a comme nouvelle matrice et comme
second membre    
1 −3 −1 0 2
 0 −2 −1 2  
 5 .
 
A(1) =  et b (1)
=
0 1 −1 0  1
0 7 2 −1 −4
Dans la deuxième étape c’est la deuxième ligne qui joue le role de ligne de pivot si x2 est
présent (Sinon, on permute la deuxième équation avec une autre sans toucher la première).
Le coefficient de x2 devient le nouveau pivot α2 qui vaut −2 dans cet exemple. Pour annuler
les coefficients de la deuxième colonne en dessous de la diagonale, on fait les opérations
7
L3 ←− L3 + 12 L2 , L4 ←− L4 + L2 .
2
Le système équivalent a pour matrice et second membre
1 −3 −1 0
   
2
0 −2 −1 2
 et b(2) =  5  .
 
A(2) = 0 0 −3/2 1  7/2 
0 0 −3/2 6 27/2
Enfin, pour éliminer x3 de la quatrième équation, on utilise le pivot α3 = − 23 et on
fait L4 ←− L4 − L3 . La matrice et le second membre de système équivalent sont
   
1 −3 −1 0 2
0 −2 −1 2 
 5 .
 
A(3) = 
0 0 −3/2 1 et b (3)
= 7/2
0 0 0 5 10
Le dernier système est triangulaire supérieur. On
 résout
 par la méthode de remontée ce
1
0
système triangulaire supérieur pour obtenir x = 
−1 qui est aussi la solution de (S).

2
Remarque 1.3.1 Si on note U = A(3) , alors on peut vérifier qu’on a
 
1 0 0 0
−1 1 0 0
A = LU, avec L =   0 −1/2 1 0

2 −7/2 1 1
est la matrice triangulaire inférieure dont les termes diagonaux valent 1 et dont le coeffi-
(k)
a (k)
cient lik de la matrice L sous la diagonale est exactement le terme ik , où aik désigne
αk
les coefficients de A(k) , matrice d’élimination de Gauss à l’étape k.
6
Algorithme de Gauss et nombre d’opérations

Si pour i = 1, . . . , n et p = 1, . . . , n, on désigne par A[i, p : n] les termes de la ligne
i d’une matrice A situés dans l’ordre entre les colonnes p et n, alors l’algorithme de la
méthode d’élimination de Gauss pour résoudre un système Ax = b, s’écrit :

Pour k = 1, . . . , n − 1,
Recherche du pivot et permutation des lignes si necessaire


 Pour i = k + 1, . . . , n
 ìk = aik


  akk
  A[i, k : n] = A[i, k : n] − ìk A[k, k : n]
 
  bi = bi − ìk bk

 Fin de la boucle sur i
Fin de la boucle sur k.
Puis, on résout le système triangulaire obtenu.

A chaque étape k de l’élimination de xk il faut effectuer :
— (n − k) divisions,
— (n − k)2 multiplications,
— (n − k)2 additions.
Le nombre total d’opérations est donc
n−1
X n−1
X n−1
X n−1
X
(n − k) + 2 (n − k)2 = p+2 p2 .
k=1 k=1 p=1 p=1
n n
X n(n + 1) X n(n + 1)(2n + 1)
Si on utilise le fait que p= et que p2 = on tire que le
p=1
2 p=1
6
n(n − 1) n(n − 1)(2n − 1)
nombre d’opérations dans la méthode d’élimination de Gauss est +2
2 6
qui est de l’ordre de 23 n3 lorsque n est très grand.
Remarque 1.3.2 Le choix du pivot peut avoir une grande importance dans le calcul
comme le montrera l’exemple suivant.
Soit à résoudre le système
−p
10 1 x1 1
= ,
1 1 x2 0
1
 
− −p 
dont la solution exacte est  1 −1 10 .
1 − 10−p
7
Après la première étape d’élimination de Gauss, on obtient
10−p

1 x1 1
= .
0 1 − 10p x2 −10p

p p x1 0
Pour p tel que sur une machine 1 − 10 ' −10 on obtient ' qui n’est pas
x2 1
proche de la solution exacte.
Cet exemple montre l’effet des erreurs d’arrondi qui proviennent de la division par des
pivots trop petits. La stratégie de pivot partiel consiste à choisir comme pivot αk vérifiant
(k)
|αk | = max |aik | le coefficient de module maximum de la partie de la colonne k en dessous
k≤i≤n
de la diagonale.
Dans l’exemple (2.4), l’élimination de Gauss était faite sans stratégie de pivot. Avec
la stratégie de pivot partiel, et dans la première étape, on doit choisir le premier pivot
a31 = 2 au lieu de a11 = 1.
1.3.2 Factorisation LU
Dans cette partie on développe la technique de la factorisation LU qui permet aussi
de résoudre un système Ax = b en se ramenant à la résolution d’un système triangulaire
inférieur puis un système triangulaire supérieur.
Proposition 1.3.1 Soit A = (aij )1≤i,j≤n ∈ Mn (R) une matrice carrée d’ordre n telle
que toutes les n sous-matrices A(k) = (aij )1≤i,j≤k , pour k = 1, . . . , n, sont inversibles.
Alors, il existe une matrice triangulaire inférieure L dont les coefficients diagonaux
sont égaux à 1, et une matrice triangulaire supérieure U telle que A = LU . De plus cette
décomposition, dite décomposition ou factorisation LU de A, est unique.
Preuve. Unicité : On suppose que A = L1 U1 = L2 U2 avec Li , i = 1, 2 triangulaire

inférieure vérifiant (Li )kk = 1 pour tout k = 1, . . . , n, Ui , i = 1, 2 triangulaire supérieure.
Alors les matrices U1 et U2 sont inversibles, puisque 0 6= det A = det U1 = det U2 . On a
alors L−1 −1 −1
2 L1 = U2 U1 , ce qui implique que L2 L1 est à la fois une matrice triangulaire
inférieure et triangulaire supérieure, donc L−1
2 L1 est une matrice diagonale. Mais comme
−1 −1
(L1 )kk = (L2 )kk = 1 on en déduit que L2 L1 = In et par suite L2 = L1 et U2 = U1 .
Existence : Notons d’abord que si une telle factorisation existe pour U = (uij )1≤i,j≤n
triangulaire inférieure U = (uij )1≤i,j≤n est triangulaire supérieure, alors
min(i,j)
X
aij = lik ukj , pour tout 1 ≤ i, j ≤ n (1.4)
k=1
On montre l’existence par récurrence sur n, l’ordre de la matrice A à factoriser.
8

a b
Si n = 2, A = vérifiant a 6= 0 et det(A) = ad − bc 6= 0, alors A admet la
c d
factorisation ! !
1 0 a b
A= c ad − bc .
1 0
a a
Supposons que le résultat est vrai jusqu’à l’ordre n − 1. Soit A = (aij ) une matrice
d’ordre n telle que les sous-matrices A(k) , pour k = 1, . . . , n sont inversibles. D’après
l’hypothèse de récurrence la matrice A(n−1) qui est d’ordre n − 1 admet une factorisation
L̃Ũ , où L̃ = (lij )1≤i,j≤n−1 est une matrice triangulaire inférieure d’ordre n − 1 à diagonale
unité et Ũ = (uij )1≤i,j≤n−1 est une matrice triangulaire supérieure d’ordre n − 1 dont
les termes diagonaux sont tous non nuls (uii 6= 0 ∀ i = 1, . . . , n − 1). On détermine les
coefficients 


 u1n = a1n ,

u2n = a2n − l21 u2n ,


..


.




p−1


 X
upn = apn − lpk ukn ,

 k=1
..


.




n−1



 X
unn = ann − lnk ukn ,



k=1
puis les coefficients  an1

 ln1 = ,


 u11

 an2 − ln1 u12
ln2 = ,


u22




 ..
.






 p−1
X
anp −

 lnk ukp
k=1
lnp =
upp



..





 .

 n−2

 X



 an,n−1 − lnk ukn
 k=1
ln,n−1 = .


un−1,n−1
Ensuite, on définit les deux matrices d’ordre n :
   
0 u1n
L=
 ..  et U =  ..  .
L̃ .  Ũ . 
ln1 · · · ln,n−1 1 0 ··· 0 unn
9
On note que L est triangulaire inférieure à diagonale unité, U est triangulaire supérieure
avec coefficients diagonaux non nuls et on peut vérifier facilement que A = LU .
Algorithme de la factorisation LU
De (1.4) et pour p = 1, . . . , n, en lisant dans l’ordre les termes de la p-ième ligne
eu dessus de la diagonale de la matrice à factoriser A, on obtient la p-ième ligne de U ,
puis au fur et on mesure qu’on écrit les termes de la p-ième colonnes de A en dessus de
la diagonale, on déduit la p-ième colonne de L. Ceci peut se traduire par l’algorithme
suivant :

Pour
 p = 1, . . . , n
 Pour j = p, . . . , n
p−1

  X
 upj := apj − lpk ukj
 

 
 k=1

 Fin de la boucle sur j

 

 Pour i = p + 1, . . . , n
p−1
 !
 X
lip := aip − lik ukp /upp
 
 
 
 k=1
 Fin de la boucle sur i
Fin de la boucle sur p
Remarque 1.3.3 De la preuve ou de l’algorithme de la factorisation LU de la matrice A

on tire que toutes les sous matrices A(k) de A admettent aussi une factorisation L(k) U (k)
pour L(k) = (lij )1≤i,j≤k et U (k) = (uij )1≤i,j≤k , pour k = 1, . . . , n. Ainsi la matrice A admet
une factorisation LU si et seulement si toutes les matrices A(k) sont inversibles.
Nombre d’opérations
A chaque étape p, pour calculer upj , j = p, . . . , n il faut effectuer
• p − 1 multiplications,
• p − 1 additions,
donc 2(p − 1) opérations.
Le nombre total d’opérations pour déterminer la p-ième ligne de U est
n
X
2(p − 1) = 2(n − p + 1)(p − 1).
j=p
Ensuite pour calculer le terme lip , i = p + 1, . . . , n il faut
10
• 1 division,
• p − 1 multiplications,
• p − 1 additions.
Donc le calcul de la p-ième colonne de L nécessite
X n
(2(p − 1) + 1) = (n − p)(2p − 1).
i=p+1
Le nombre total d’opérations pour effectuer la factorisation LU de A est

Xn Xn
[2(n − p + 1)(p − 1) + (n − p)(2p − 1)] = 2 (n − p)p + 2.
p=1 p=1
n n
X n(n + 1)(2n + 1) X n(n + 1)
Si on utilise le fait que p2 = et que p= on tire que
p=1
6 p=1
2
le nombre d’opérations pour effectuer la factorisation LU d’une matrice d’ordre n est de
l’ordre de 32 n3 .
Résolution de Ax = b pour A = LU
Le système Ax = b donne LU x = b. Si on pose y = U x, il suffit de résoudre le système
triangulaire inférieur Ly = b puis le système triangulaire supérieur U x = y, i.e.,
(
Résoudre Ly = b,
Résoudre Ax = b ⇔
puis résoudre U x = y.
Puisque la résolution d’un système triangulaire, nécessite n2 opérations, le nombre
d’opérations pour résoudre un système linéaire Ax = b, par la factorisation A = LU est
de l’ordre 32 n3 .
Remarques 1.3.1
1. On utilise la factorisation LU surtout d’une matrice A pour résoudre plusieurs
systèmes linéaires avec la même matrice A mais avec des seconds membres différents.
Si on a par exemple à résoudre p systèmes linéaires avec la même matrice A, le
coût est de l’ordre de 32 n3 + 2pn2 au lieu de p 32 n3 .
2. Puisque
Y n
det A = det L det U = uii ,
i=1
la factorisation LU permet aussi de calculer det A en seulement 23 n3 opérations
au lieu de (n! − 1)(n − 1) opérations si on utilise la formule de Leibniz
X n
Y
det A = (σ) aσ(i),i ,
σ∈Sn i=1
où Sn est l’ensemble des permutations des nombres {1, . . . , n} et (σ) est la si-
gnature de la permutation σ.
11
Matrices bandes
Définition 1.3.1 Une matrice A = (aij )1≤i,j≤n ∈ Mn (C) est dite matrice bande s’il
existe nb ∈ N, tel que aij = 0 pour |i − j| > nb . Autrement dit les coefficients non nuls
de la matrice A sont situés dans une bande de longueur 2nb + 1, centrée sur la diagonale
principale. Pour nb = 1, la matrice est dite tridiagonale.
Proposition 1.3.2 La factorisation A = LU conserve la structure bande.
Preuve. On montre cette proprièté par récurrence sur les étapes de la factorisation de la
matrice A.
La matrice A étant bande, la proprièté est vraie à l’étape 0. Supposons que la pro-
priété est vraie jusqu’à l’étape p − 1 où on a donc lij = uij = 0 pour 1 ≤ i, j ≤ p − 1,
|i − j| > nb . A l’étape p on a
p−1
X
upj := apj − lpk ukj , pour j = p, . . . , n
k=1
et
p−1
X
aip − lik ukp
k=1
lip := , pour i = p + 1, . . . , n.
upp
Soit p tel que |p − j| > nb . Alors apj = 0 puisque A est une matrice bande.
Si p − j > nb > 0 alors upj = 0, puisque U est triangulaire supérieure.
Si p − j < −nb < 0, donc pour 1 ≤ k ≤ p − 1, on a k − j ≤ p − 1 − j nb , alors apj = 0.
Si i − p < −nb < 0 alors lip = 0, puisque L est triangulaire inférieure.
Si i − p > nb < 0, donc pour 1 ≤ k ≤ p − 1 < p, on a i − k ≤ i − p < −nb . D’après
l’hypothèse de récurrence lik = 0. Ainsi lip = 0 et la propriété est vérifiée aussi pour la
matrice L à l’étape k.
1.4 Factorisation de Cholesky des matrices symétriques

définies positives
1.4.1 Rappels sur les matrices symétriques
Définition 1.4.1 Soit A = (aij )1≤i,j≤n ∈ Mn (R). La matrice transposée de A, notée AT
est définie par
AT = (aji )1≤i,j≤n .
L’opérateur transposé vérifie, pour toutes marices A et B de Mn (R),
12
1.4 Factorisation de Cholesky des matrices symétriques définies positives
1. (A + B)T = AT + B T
2. (AB)T = B T AT
3. Si A est inversible, alors AT est inversible et (AT )−1 = (A−1 )T .
4. (Ax, y) = (x, AT y) pout tout x, y ∈ Rn .
Définition 1.4.2 Soit A ∈ Mn (R). La matrice A est dite :

1. symétrique si AT = A.
2. orthogonale si AAT = AT A = In .
Remarque 1.4.1 Si A ∈ Mn (R) est une matrice orthogonale, alors A est inversible et
A−1 = AT puisqu’on a AAT = AT A = In . De plus les vecteurs lignes ainsi que les vecteurs
colonnes forment une base orthonormée de Rn .
Définition 1.4.3 Vecteurs et valeurs propres

Soit A ∈ Mn (R). On rappelle que λ ∈ C est une valeur propre de A si det(A−λIn ) =
0. Un vecteur x ∈ Cn non nul est dit vecteur propre de A associé à la valeur propre λ si
Ax = λx.
Un résultat connu concernant la réduction de matrices symétriques est le suivant :

Proposition 1.4.1 Soit A ∈ Mn (R) une matrice symétrique. Alors, il existe une ma-
trice orthogonale O ∈ Mn (R) telle que la matrice D = OT AO soit diagonale réelle.
Remarques 1.4.1 Si A ∈ Mn (R) symétrique et O une matrice orthogonale telle que

 
λ1
OT AO = D = 
 .. ,

.
λn
alors λ1 , λ2 , . . . , λn sont les valeurs propres de A. De plus si f1 , f2 , . . . , fn sont les vecteurs

colonnes de O, alors fi est un vecteur propre de A associé à la valeur λi pour i = 1, . . . , n
et la famille (f1 , . . . , fn ) forme une base orthonormée des vecteurs propres de A.
Conséquences 1.4.1
1. Toute matrice symétrique est diagonalisable.
2. Les valeurs propres d’une matrice symétrique sont toutes réelles.
Définition 1.4.4 Soit A ∈ Mn (R) une matrice symétrique. La matrice A est dite semi-
définie positive si
(Ax, x) ≥ 0 pour tout x ∈ Rn .
Si de plus A vérifie, (Ax, x) = 0 si et seulement si x = 0, alors A est dite définie positive.
Proposition 1.4.2 Soit A ∈ Mn (R) une matrice symétrique.
13
1. La matrice A est semi-définie positive si et seulement si toutes ses valeurs propres

sont positives.
2. La matrice A est définie positive si et seulement toutes ses valeurs propres sont
strictement positives.
Preuve. Soit A ∈ Mn (R). On suppose que A est semi-définie positive (resp. définie
positive). Soit λ une valeur propre de A et x ∈ Rn (non nul) un vecteur propre associé à
la valeur propre λ. On a Ax = λx, donc
(Ax, x) = λ(x, x) = λkxk22 .
Comme A est semi-définie positive (resp. définie positive), donc (Ax, x) ≥ 0 (resp. > 0)
et par conséquent
λkxk22 ≥ 0 (resp. > 0)
ou λ ≥ 0 (resp. > 0).
Supposons que toutes les valeurs propres de A sont positives (resp. strictement po-
sitives) et montrons que A est semi-définie positive (resp. strictement positives). On sait
qu’il existe une base orthonormée (f1 , . . . , fn ) de vecteurs propres de A. Soit x ∈ Rn ,
(resp. x ∈ Rn \ {0},
X n
x= xi f i .
i=1
Donc
n
X
Ax = xi λ i f i ,
i=1
et n
X
(Ax, x) = λi x2i ≥ 0 (resp. > 0). (1.5)
i=1
Remarque 1.4.2 De (1.5) on tire facilement que si A est une matrice symétrique d’ordre
n et si λ1 , λn sont respectivement la plus petite et la plus grande valeur propre de A, alors
pour tout x ∈ Rn , on a
λ1 kxk22 ≤ (Ax, x) ≤ λn kxk22 . (1.6)
L’inégalité (1.6) sera fort utile pour la suite.
1.4.2 Factorisation des matrices symétriques

Proposition 1.4.3 Soit A = (aij )1≤i,j≤n ∈ Mn (R) une matrice symétrique définie posi-
tive. Alors la matrice A admet une factorisation LU avec les éléments diagonaux de U
sont strictement positifs.
14
Preuve. Soit pour k = 1, . . . , n, A(k) = (aij )1≤i,j≤k ∈ A(k) (R). Montrons que la matrice
symétrique A(k) est inversible. En effet, pour tout x = (x1 , . . . , xk )T ∈ Rk \ {0}, on pose
x̃ = (x1 , . . . , xk , 0, . . . , 0) ∈ Rn . Alors x̃ 6= 0 et on a (A(k) x, x) = (Ax̃, x̃) > 0. La matrice
A(k) est par conséquent définie positive, donc inversible. La matrice A admet donc une
factorisation LU . Montrons que uii > 0. Posons la matrice diagonale
 
u11
D=
 .. .

.
unn
Puisque
n
Y
det A = det L det U = uii 6= 0,
i=1
donc la matrice diagonale D est inversible. Soit

uij
S = D−1 U = ( )1≤i,j≤n ,
uii
alors U = DS et
A = LDS = S T DLT .
La matrice S T est triangulaire inférieure à diagonale unité ((S T )ii = 1, i = 1, . . . , n) et
la matrice DLT est une matrice triangulaire supérieure. D’après l’unicité de factorisation
LU , on a L = S T et U = DLT . Ainsi A = LDLT .
Soit x ∈ Rn \ {0} et y tel que x = LT y, alors y est aussi non nul et on a
(Dx, x) = (DLT y, LT y) = (LDLT y, y) = (Ay, y) > 0,
donc la matrice diagonale symétrique D est définie positive, ses valeurs propres uii , i =
1, . . . , n sont donc strictement positives.
Remarque 1.4.3 (Factorisation de Crout)

De la preuve de la proposition précédente on déduit que si A est une matrice symétrique in-
versible et admettant une factorisation LU , alors A admet la factorisation A = LDLT où
D est une matrice diagonale. La factorisation LDLT est appelée factorisation de Crout.
1.4.3 Factorisation de Cholesky

Proposition 1.4.4 Soit A une matrice symétrique définie positive. Alors il existe une
matrice triangulaire inférieure B à éléments diagonaux strictement positifs telle
que A = BB T . De plus cette factroisation est unique. Une telle factorisation est dite
factorisation de Cholesky de A.
15
Preuve. On sait que A admet une factoriasion LDLT = LU avec (L)ii = 1 et (D)ii =
uii > 0 pour i = 1, . . . , n. Si on pose
√ 
u11
D0 = 
 ... ,

√
unn
alors
A = (LD0 )(D0 LT ) = BB T .
L’unicité découle de celle de la factorisation LU .
Remarque 1.4.4 La matrice B de la factorisation de Cholesky de la matrice symétrique

√
A = LU vérifie B = L diag( uii ) où
√ 
u11
√ ...
diag( uii ) =  .
 
√
unn
Algorithme de la factorisation de Cholesky

Il suffit d’identifier le produit A et BB T pour B = (bij ) une matrice triangulaire
inférieure. On en déduit
min(i,j)
X
aij = bik bjk ∀ 1 ≤ i, j ≤ n.
k=1
On commence par calculer terme à terme la première colonne de la matrice triangulaire B,

puis on passe à la deuxième colonne et ensuite connaissant les (i−1) premières colonnes de
B, on peut déduire les termes de la i-ème colonne se trouvant au dessous de la diagonale
jusqu’on on obtient toute la matrice triangulaire inférieure B. Ceci se traduit aussi par
l’algorithme suivant :

Pour i =v1, . . . , n
 u
Xi−1
 u

 bii := taij − b2ik
 k=1
 

 Pour j = i + 1, . . . , p − 1!
  Xi−1
 bij := aij − bik bkj /bii
 

 
 k=1
 Fin de la boucle sur j
Fin de la boucle sur i
16
3
Remarque 1.4.5 Clairement, le coût de la méthode de Cholesky est de l’ordre de n3
qui est la moitié de la factorisation LU , puisque ici le calcul de la matrice B donne
immédiatement sa matrice transposée B T .
17
18
Chapitre 2
Méthodes itératives pour la

résolution des systèmes linéaires
2.1 Introduction
Pour des systèmes linéaires de grande taille, les méthodes directes de factorisations (de
type LU ou de Cholesky) deviennent couteuses en temps de calcul ou en place mémoire.
L’idée alors de ne plus chercher à résoudre exactement le système linéaire, mais d’appro-
cher sa solution x par une suite de vecteurs (x(k) ) vérifiant
lim kx(k) − xk = 0.
k→+∞
Dans ce chapitre, la suite (x(k) ) est construite à l’aide d’une relation de récurrence simple
de type x(k+1) = F (x(k) ), pour une application affine F (x) = Bx + c.
2.2 Rappels sur les normes matricielles

On rapelle que l’espace vectoriel complexe Cn est muni du produit hermitien défini
pour x = (xi )1≤i≤n et y = (yi )1≤i≤n ∈ Cn , par
n
X
(x, y) = xi ȳi ,
i=1
et la norme euclidienne associé à ce produit scalaire est

n
X 1
kxk2 = ( |xi |2 ) 2 .
i=1
Dans ce chapitre, on note par Mn (K) l’espace vectoriel des matrices carrées d’ordre
n à coefficients dans K où K = R ou C. On rappelle que pour x = (xi )1≤i≤n ∈ Kn on a
les normes vectorielles suivantes :
19
Chapitre 2. Méthodes itératives pour la résolution des systèmes linéaires
n
X 1
• kxkp = ( |xi |p ) p , pour p ∈ [1, +∞[ ; norme de Hölder d’indice p.
i=1
• kxk∞ = max |xi |, norme infinie.
1≤i≤n
2.2.1 Normes matricielles subordonnées

Définition 2.2.1 On appelle norme matricielle sur Cn toute norme k · k sur Mn (C)
vérifiant kABk ≤ kAkkBk, c.-à-d. toute application
k · k : Mn (C) → R+
A 7 → kAk
vérifiant pour tout A, B ∈ Mn (C) et pour tout α ∈ C :

— kAk = 0 ⇔ A = 0n ; matrice nulle d’ordre n,
— kαAk = |αkAk,
— kA + Bk ≤ kAk + kBk,
— kABk ≤ kAkkBk.
Définition 2.2.2 (Normes matricielles subordonnées)

Toute norme vectorielle k · k de Cn définit une norme matricielle de la façon suivante
kAxk
∀A ∈ Mn (C), kAk = sup = sup kAxk = sup kAxk,
x∈Cn , x6=0 kxk kxk≤1 kxk=1
dite norme matricielle subordonnée ou induite (à cette norme vectorielle).
Cette norme matricielle subordonnée vérifie donc

1. kAxk ≤ kAkkxk, pour toute matrice A ∈ Mn (C) et pour tout vecteur x ∈ Cn .
2. kIn k = 1.
On notera k · kp la norme matricielle subordonnée assoiée à la norme vectorielle d’indice
p.
Exemple 2.2.1 La norme (dite de Frobenius) définie, pour A = (aij )1≤i,j≤n ∈ Mn (R),
par
n
X 1 1
kAk = ( a2ij ) 2 = (trace(AT A)) 2 ,
i,j=1
√
n’est pas subordonnée puisque kIn k = n 6= 1.
Définition 2.2.3 (Rayon spectral)

Soient A ∈ Mn (C) et λi ∈ C, 1 ≤ i ≤ n les valeurs propres de A. On rappelle que le
spectre de A, qu’on note Sp(A), est l’ensemble des valeurs propres de A.
On appelle rayon spectral de A, le réel positif, noté ρ(A) = max |λi | qui est le maxi-
1≤i≤n
mum des modules des valeurs propres de A.
20
On a alors :
Proposition 2.2.1 Les normes matricielles subordonnées aux normes vectorielles k · k1 ,
k · k2 et k · k∞ sont
(pdonnées, pour une matrice A = (aij )1≤i,j≤n , par :
p
ρ(AT A) = ρ(AAT ) = kAT k2 si A est une matrice réelle,
— kAk2 = p p
ρ(A∗ A) = ρ(AA∗ ) = kA∗ k2 si A est une matrice complexe,
Xn
— kAk1 = max |aij |,
1≤j≤n
i=1
Xn
— kAk∞ = max |aij |.
1≤i≤n
j=1
Preuve. Soit A = (aij )1≤i,j≤n ∈ Mn (C) une matrice non nulle.

1. Si A ∈ Mn (R),
kAxk22 (Ax, Ax) (AT Ax, x)

kAk22 = sup 2
= sup = sup .
x∈Rn ,x6=0 kxk2 x∈Rn ,x6=0 kxk22 x∈Rn ,x6=0 kxk22
La matrice B = AT A étant symétrique semi définie positive puisque B T = B et

(Bx, x) = (Ax, Ax) = kAxk22 ≥ 0, pour tout x ∈ Rn . D’après (1.6), on a
kAk22 ≤ ρ(AT A).
D’autre part, si x est un vecteur propre de B = AT A associé à la valeur propre

ρ(AT A), alors
kAxk22 = (Bx, x) = ρ(AT A)kxk2 ,
donc
kAxk2 p
= ρ(AT A) ≤ kAk2 ,
kxk2
et finalement p
kAk2 = ρ(AT A).
2. A faire en exercice.
3. Soit x = (xi )1≤i≤n ∈ Cn . On a alors,
n
X n
X n
X
kAxk∞ = max |(Ax)i | = max | aij xj | ≤ max |aij | |xj | ≤ ( max |aij |)kxk∞ ,
1≤i≤n 1≤i≤n 1≤i≤n 1≤i≤n
j=1 j=1 j=1
n
X
et par conséquent kAk∞ ≤ max |aij |. D’autre part, si i0 est tel que
1≤i≤n
j=1
n
X n
X
|ai0 j | = max |aij |.
1≤i≤n
j=1 j=1
21
Comme A 6= 0, alors sa i0 -ème colonne est non nulle et il existe au moins j tel
que ai0 j 6= 0. Posons, pour j = 1, . . . , n,
 āi0 j

si ai0 j 6= 0
xj = |ai0 j | .
0 sinon
Alors kxk∞ = 1 et
n n n
X āi0 j X X
|(Ax)i0 | = | ai 0 j |= |ai0 j | = max |aij |.
j=1
|ai0 j | j=1
1≤i≤n
j=1
n
X
Par suite kAk∞ = max |aij |.
1≤i≤n
j=1
Remarque 2.2.1 Si A ∈ Mn (R) est une matrice symétrique, alors
kAk2 = ρ(A). (2.1)
Car les valeurs propres de A2 sont les carrées des valeurs propres de A lorsque cette
dernière est diagonalisable et donc
ρ(A2 ) = ρ(A)2 .
La relation (2.1) devient une inégalité pour une matrice quelconque et pour les autres
normes. Plus précisément :
Proposition 2.2.2 On a :
1. Pour toute matrice A et pour toute norme matricielle k · k subordonnée ou non
on a
ρ(A) ≤ kAk.
2. Pour tout > 0 et pour toute matrice A ∈ Mn (C), il existe une norme matricielle
subordonnée notée k · kA, telle que
kAkA, ≤ ρ(A) + .
Preuve.
1. Soient A ∈ Mn (C), λ une valeur propore de A telle que |λ| = ρ(A) et x ∈ Cn , x 6=
0 un vecteur propore associé à la valeur propre λ. Alors Ax = λx.
Si y ∈ Cn tel que la matrice xy T ne soit pas nulle, alors Axy T = λxy T et par suite
kAxy T k = ρ(A)kxy T k ≤ kAkkxT k.
Ce qui donne ρ(A) ≤ kAk.
22
2. Soit A ∈ Mn (C), on sait que A est semblable à une matrice triangulaire supérieure
complexe, il existe donc (f1 , . . . , fn ) une base de Cn et une matrice triangulaire
supérieure (λij )1≤i,j≤n telles que
X
Afi = λij fj ,
j≥i
avec, si P = (f1 , . . . , fn ), alors P −1 AP = T .

Pour un réel η ∈]0, 1[, on définit B = (ẽ1 , ẽ2 , . . . , ẽn ) la base de l’espace vectoriel
Cn par,
ẽ1 = f1 , ẽ2 = ηf2 , ...., ẽn = η n−1 fn .
n
X
Si x = αi ẽi ∈ Cn , on définit la norme vectorielle de x par
j=1
kxk := max |αj |.

1≤j≤n
Il est clair que cette norme dépend de A et de η.

Soit > 0, montrons que si η est bien choisi, la norme subordonnée associée vérifie
kAk ≤ ρ(A) + .
Pour tout i ∈ {1, . . . , n},

X X
Aẽi = η i−1 Afi = η i−1 λij fj = η i−1 η 1−j λij ẽj .
1≤i≤j 1≤i≤j
Donc X
Aẽi = η i−j λij ẽj .
1≤i≤j
n
X
Si x = αi ẽi , alors
i=1
n
X X n X
X n
i−j
Ax = αi η λij ẽj = ( η i−j λij αi )ẽj .
i=1 1≤i≤j j=1 i=j
n
" n
#
X X
kAxk = max | η i−j λij αi | ≤ max |λjj kαj | + η |λij kαi |
1≤j≤n 1≤j≤n
i=j i=1
n
X
≤ ρ(A)kxk + η max |λij |kxk.
1≤j≤n
i=1
23
Par conséquent
kAxk
≤ ρ(A) + ηkTA k1 ,
kxk
donc
kAk ≤ ρ(A) + ηkTA k1 .
Pour η ∈]0, 1[ tel que ηkTA k1 ≤ , on a bien
k · kA,ε = kAk ≤ ρ(A) + .
On utilse ce dernier résultat pour montrer la proposition suivante :
Proposition 2.2.3 Soit B une matrice carrée de Mn (C). Les conditions suivantes sont
équivalentes :
i) lim B k x = 0 pour tout vecteur x de Cn ,
k→+∞
ii) ρ(B) < 1,

iii) kBk < 1 pour au moins une norme matricielle subordonnée.
Preuve.
i) ⇒ ii) Soit λ une valeur propre de B telle que ρ(B) = |λ| et soit x un vecteur
propre associé à la valeur propre λ de B. On a Bx = λx, donc B k x = λk x et par
conséquent kB k xk = ρ(B)k kxk → 0 quand k → +∞ si et seulement si ρ(B) < 1.
1 − ρ(B)
ii) ⇒ iii) Pour = , il existe une norme matricielle subordonnée k · k telle
2
que
1 − ρ(B) ρ(B) + 1
kBk ≤ ρ(B) + < ρ(B) + = < 1.
2 2
iii) ⇒ i) On a kB k xk ≤ kBkk kxk → 0 quand k → +∞.
2.3 Méthodes itératives

Pour résoudre un système linéaire Ax = b pour A ∈ Mn (R) une matrice inversible,
on utilise dans ce chapitre des méthodes itératives dont le principe est d’écrire A comme
la différence de deux matrices A = M − N , où M est une matrice inversible (M est en
général diagonale, triangulaire ou facile à inverser). Le système Ax = b est équivalent à
x = M −1 (N x + b). On approche donc la solution du système Ax = b par la suite (x(k) )
défine par
x(k+1) = M −1 (N x(k) + b), k > 0,
en partant d’un x(0) donné.
24
Remarques 2.3.1
1. On n’a pas toujours besoin de calculer M −1 , mais il faut savoir calculer la solution
de M x(k+1) = N x(k) + b, pour x(0) donnée.
2. Il est clair que si la suite (x(k) ) converge, elle converge vers la solution unique
x de Ax = b. On dit dans ce cas que la méthode itérative correspondante est
convergente pour la résolution de (S).
3. Si on considère e(k) l’erreur à l’étape k, e(k) = x(k) − x, alors M x(k+1) = b + N x(k)
et M x = b + N x et par conséquent e(k+1) = M −1 N e(k) = · · · = (M −1 )k e(k) . La
matrice M −1 N = B est appelée matrice d’itération de la méthode. La suite
(x(k) ) vérifie donc x(k+1) = Bx(k) + c, pour c = M −1 b.
En général on a :
Proposition 2.3.1 Si B = M −1 N et c = M −1 b, alors la suite (x(k) ) donnée par
x(k+1) = Bx(k) + c, pour x(0) donné, (2.2)
converge vers x la solution unique de x = Bx + c, pour tout choix de x(0) , si et seulement

si la matrice d’itération B vérifie ρ(B) < 1.
Preuve. La suite donnée par x(0) et x(k+1) = Bx(k) + c converge vers x pour tout choix
de x(0) si et seulement si (B k e(0) ) → 0 quand k → +∞ pour tout e(0) = x(0) − x. Ceci
étant vrai si et seulement si ρ(B) < 1.
Remarque 2.3.1 Ainsi, pour montrer que la suite générée par la méthode itérative (2.2)
est convergente, il suffit de vérifier que le rayon spectral de sa matrice d’itération B est
plus petit que 1 ou encore que pour au moins une norme quelconque k · k, (subordonnée
ou non) kBk < 1.
Un premier résultat de convergence d’une méthode itérative concerne les matrices

symétriques définies positives.
Proposition 2.3.2 (Cas d’une matrice symétrique définie positive)

Soit A une matrice symétrique définie positive et M, N deux matrices telles que A =
M − N , avec M inversible. Si la matrice symétrique M T + N est définie positive, alors
ρ(M −1 N ) < 1 et la suite définie par x(k+1) = M −1 N x(k) + M −1 b converge vers x solution
de Ax = b pour tout choix de x(0) .
Preuve. Si A est symétrique, alors la matrice M T + N est toujours symétrique. En effet
(M T + N )T = M + N T = (A + N ) + N T = AT + N T + N = (A + N )T + N = M T + N.
Montrons que si A est définie positive et M T + N est définie positive, alors
ρ(M −1 N ) < 1.
25
Soit λ ∈ Sp(M −1 N ) et soit x ∈ Cn , x = y + iz 6= 0, y, z ∈ Rn tel que M −1 N x = λx. On

a M −1 N x = λx, donc
N x = λM x. (2.3)
Par suite
(N x, x) = λ(M x, x).
Comme A = M − N , donc A = M (I − M −1 N ) et par conséquent
Ax = M (x − M −1 N x) = M (x − λx) = (1 − λ)M x.
D’où
(Ax, x) = (Ay, y) + (Az, z) = (1 − λ)(M x, x) > 0,
et donc λ 6= 1. De plus on a
(N x, x) = λ(M x, x) et (M T y, y) = (y, M y) = (M y, y)
impliquant
((M T + N )x, x) = ((M T + N )y, y) + ((M T + N )z, z) = (1 + λ)(M x, x)
1+λ 1+λ
= (Ax, x) = [(Ay, y) + (Az, z)].
1−λ 1−λ
Comme les matrices A et M T + N sont définies positives, alors nécessairement le nombre
1−λ
complexe λ est tel que > 0. Ceci est équivalent à λ ∈ R et −1 < λ < 1.
1+λ
Dans les trois méthodes itératives classiques qu’on va étudier, on utilise la décomposition
suivante de A = (aij )1≤i,j≤n ∈ Mn (R) :
A = D − E − F,
avec
— La matrice diagonale  
a11
D=
 .. .

.
ann
 
0 0 ... ... 0
 −a21 0 ... ··· 0
 .. 
— La matrice triangulaire inférieure E =  −a31 −a32 0 .  qui représente
 
 . ... ...
 ..


−an1 −an2 . . . −ann−1 0
l’opposé de la partie en dessous de la diagonale.
0 −a12 −a13 · · · −a1n
 
 ..
. 0 −a23 ... −a2n 

— La matrice triangulaire supérieure F = 
 .. .. ..  qui
 . . .  
0 0 −an−1n 
0 ... ... ... 0
représente l’opposé de la partie en dessus de la diagonale.
26
2.3.1 Méthode de Jacobi

Dans la méthode de Jacobi on choisit M = D et donc N = E + F = D − A où on
suppose que la matrice diagonale D est inversible (aii 6= 0 pour tout i = 1, . . . , n). La
matrice d’itération est J = D−1 (E + F ) et la suite x(k) dans la méthode de Jacobi est
donnée par :
x(0) ∈ Rn donnée,

Pour k = 0, 1, 2, . . .
Dx(k+1) = (E + F )x(k) + b
Les composantes de x(k+1) sont données en fonction de celles de x(k) par

" #
(k+1) 1 X (k)
xi = bi − aij xj , i = 1, . . . , n.
aii i6=j
(k+1)
Il est bien de noter que les composantes xi de x(k+1) sont calculées les uns indépendement
des autres ce qui peux s’effectuer en parallèle.
On montrera la convergence de cette méthode pour les matrices à diagonale stricte-
ment dominante :
Définition 2.3.1 Une matrice A = (aij )1≤i,j≤n ∈ Mn (C) est dite à diagonale strictement
dominante si elle vérifie
X
|aii | > |aij | ∀ i ∈ {1, . . . , n}.
j6=i
Proposition 2.3.3 Soit A une matrice à diagonale strictement dominante. Alors A est
inversible et la méthode de Jacobi converge vers la solution du système Ax = b.
Preuve. Montrons que A est inversible. Soit x = (xi )1≤i≤n ∈ Cn tel que Ax = 0. Alors,
pour tout i = 1, . . . , n,
Xn
aij xj = 0.
j=1
Ou encore n
X
aii xi = aij xj .
j6=i
Par conséquent
n
X
|aii kxi | ≤ |aij kxj |.
j6=i
27
Si x 6= 0, alors kxk∞ 6= 0 et il existe i0 ∈ {1, . . . , n} tel que
|xi0 | = max |xi | = kxk∞ .

1≤i≤n
Ainsi
X |xj | X
|ai0 i0 | ≤ |ai0 j | ≤ |ai0 j | < 1.
j6=i0
|xi0 | j6=i
0
Ce qui est une contradiction, donc x = 0 et A ne peut être que inversible.

Pour la méthode de Jacobi on a
− aa11 − aa1n
 
0 12
... 11
 − a21 0 ... − aa2n 
 a2 22 
J = . ..
 ..

. 
an1 ann−1
− ann . . . − ann 0
qui vérifie
X |aij |
kJk∞ = max < 1.
1≤i≤n
j6=i
|aii |
Comme
ρ(J) < kJk∞ < 1,
et donc la méthode de Jacobi converge.
2.3.2 Méthodes de relaxation

Dans la méthode de relaxation, on introduit un paramètre réel non nul ω et on prend
1 1−ω
M= D − E, N = M − A = D + F.
ω ω
La matrice d’itération est donc
Lω = ( ω1 D − E)−1 ( 1−ω
ω
D + F)
= ω(I − ωD−1 E)D−1 D ω1 ((1 − ω)In + ωD−1 F )
= (I − ωD−1 E)((1 − ω)In + ωD−1 F ).
— Si ω < 1, on parle de sous relaxation.

— Si ω > 1, on parle de sur relaxation.
— Si ω = 1 la méthode est dite de Gauss-Seidel.
Indépendamment de la matrice A, on a le résultat de divergence suivant :
Proposition 2.3.4 La méthode de relaxation diverge pour tout ω ∈ R\]0, 2[.
28
Preuve. On a
det(Lω ) = det(In − ωD−1 E)−1 det((1 − ω)In + ωD−1 F ).
La matrice (In − ωD−1 E) est triangulaire inférieure dont les coefficients diagonaux sont
égaux à 1, donc (In − ωD−1 E)−1 est aussi triangulaire inférieure à coefficients diagonaux
sont aussi 1. La matrice (1 − ω)In + ωD−1 F est triangulaire supérieure ayant 1 − ω comme
termes diagonaux. Ainsi det(Lω ) = (1 − ω)n . Comme |1 − ω|n = | det(Lω )| ≤ ρ(Lω )n , donc
|1 − ω| ≤ ρ(Lω ) et si ω ∈ R \ [0, 2], alors |1 − ω| > 1. Par conséquent ρ(Lω ) > 1 et la
méthode de relaxation diverge.
Ainsi, la méthode de relaxation n’est intéressante que pour 0 < ω < 2.
Proposition 2.3.5 Soit A ∈ Mn (R) une matrice symétrique définie positive. Alors la
méthode de relaxation converge si et seulement ω ∈]0, 2[.
Preuve. Il suffit de montrer que M T + N est définie positive.
1 1−ω 1 1−ω 2−ω

M T + N = ( D − E)T + D+F = D−F + D+F = D.
ω ω ω ω ω
Comme A est définie positive, les coefficients diagonaux de A qui sont les éléments de la
diagonale de la matrice D, sont strictement positifs, puisque chaque terme aii = (Aei , ei ),
pour i = 1, . . . , n, pour ei est le i−ième vecteur de la base canonique de Rn . La matrice
M T + N est donc définie positive si et seulement si 2−ω ω
> 0. Ce qui est équivalent à
ω ∈]0, 2[.
2.3.3 Méthode de Gauss-Seidel

C’est la méthode de relaxation pour ω = 1. On choisit donc pour cette méthode la
matrice triangulaire inférieure M = D −E et donc N = F qui est une matrice triangulaire
supérieure à diagonale nulle. On suppose que M est inversible (aii 6= 0 pour tout i =
1, . . . , n). Dans ce cas la matrice d’itération est
L1 = (D − E)−1 F,
et la méthode de Gauss-Seidel s’écrit :
x(0) ∈ Rn donnée,

Pour k = 0, 1, 2, . . .
(D − E)x(k+1) = F x(k) + b
29
(k+1)
La i-ème composante de x(k+1) est donnée en fonction des composantes xj ,j i de x(k) par
" #
(k+1) 1 X (k+1)
X (k)
xi = bi − aij xj − aij xj , i = 1, . . . , n.
aii ji
(k+1)
Il faut noter que le calcul des composantes xi de x(k+1) doit s’effectuer dans l’ordre
i = 1, . . . , n.
On sait déja que si A est une matrice symétrique définie positive, alors la méthode
de Gauss-Seidel est convergente (ω = 1 ∈]0, 2[). On montre aussi que
Proposition 2.3.6 Si A ∈ Mn (R) une matrice à diagonale strictement dominante,

alors, la méthode de Gauss-Seidel pour la résolution de système Ax = b est convergente.
Preuve. Montrons que si A est à diagonale strictement dominante, alors ρ(L1 ) < 1.
Soit x = (xi )1≤i≤n un vecteur propre de L1 pour une valeur propre λ de L1 =
(D − E)−1 F . Alors L1 x = λx et donc λ(D − E)x = F x.
Par conséquent
X X
λ(aii xi + aij xj ) = aij xj ∀i ∈ {1, . . . , n}.
ji
Pour i ∈ {1, . . . , n}, si xi 6= 0, alors

X xj X xj
λ(aii + aij )= aij ∀i ∈ {1, . . . , n}.
ji
x i
Ainsi
X |xj | X |xj |
|aij | |aij |
j>i
|xi | j>i
|xi |
|λ| ≤ X xj ≤ .
|aii + aij |
X |xj |
xi |aii | − |aij |
ji0
|xi0 |
|λ| ≤ < 1.
X |xj |
|ai0 i0 | − |ai0 j |
j<i0
|xi0 |
Par suite ρ(L1 ) < 1.
30
Vitesse de convergence
En pratique, il faut tenir compte de la vitesse de convergence. Autrement dit, entre
deux méthodes itératives convergentes, on choisit celle qui donne la solution plus rapide-
ment que l’autre. Un critère de mesure de cette vitesse est l’évolution de l’erreur kx(k) −xk
qui dépend du rayon spectral de la matrice d’itération. En effet, si on a une suite conver-
gente (x(k) ) définie par x(k+1) = Bx(k) + c, pour x(0) donné, sa limite est x vérifiant
x = Bx + c, et on a
x(k) − x = B(x(k−1) − x) = B k (x(0) − x)
et par conséquent
kx(k) − xk2 = kB k (x(0) − x)k2 ≤ kBkk2 kx(0) − xk2 .
En particulier, si B est symétrique, alors, kBk2 = ρ(B) et il vient
kx(k) − xk2 ≤ ρ(B)k kx(0) − xk2 .
Ainsi, la suite (x(k) ) converge plus vite vers x d’autant que ρ(B) est plus petit. Ceci reste
encore vrai pour une norme matricielle quelconque et pour une matrice B quelconque et
on énonce le résultat de comparaison de deux méthodes itératives suivant :
Proposition 2.3.7 Pour la recherche d’une solution x d’un système Ax = b, une méthode
itérative de matrice d’itération B est dite plus rapide qu’une autre de matrice d’itération
B̃, si ρ(B) ≤ ρ(B̃).
En général, on ne peut rien dire sur la convergence d’une de deux méthodes de Jacobi
et de Gauss-seidel connaissant la nature de convergence de l’autre. Cependant, pour une
matrice tridiagonale on a :
Proposition 2.3.8 : Cas d’une matrice tridiagonale

Soit A une matrice tridiagonale. Alors les rayons spectraux de Jacobi et Gauss-Seidel
sont liés par la relation
ρ(L1 ) = (ρ(J))2 .
Donc les deux méthodes convergent ou divergent simultanèment. Lorsqu’elles convergent,
la méthode de Gauss-Seidel est plus rapide.
Preuve. A = D − E − F , avec
   
a11 a12 a11
... ... ...
a
   
A =  21 . , D = ,
  
.. ...   ...
 an−1n   
ann−1 ann ann
31
−a12
   
0 0
. .. ..
−a21 . . 0 . .
   
E=  et F =  .
 
.. .. .. ..
 . .   . . −an−1n 
−ann−1 0 0
Montrons que pour tout α, µ ∈ C∗ , det(µD − E − F ) = det(µD − αE − α−1 F ). On pose
 
µa11 a12
.. ..
a
 . . 
Aµ := µD − E − F =  21 .

.. .. 
 . an−1n 
ann−1 µann
et
µa11 α−1 a12
 
... ...
αa
 
Aα,µ := µD − αE − α−1 F =  21 .

.. .. −1
 . . α an−1n 
αann−1 µann
Soit la matrice diagonale  
α
 α2 
P = .
 
 ... 
αn
Alors P est inversible et
α−1 µa11 α−2 a12

 
 −1 .. ..
 α a21 . . 
P −1 Aα,µ = ,

. .. ...
 α−n an−1n 
α1−n ann−1 α−n µann
 
µa11 a12
.. ..
a
 . . 
P −1 Aα,µ P =  21 .  = Aµ .

.. ...
 an−1n 
ann−1 µann
Donc
det(P −1 Aα,µ P ) = det Aµ .
Montrons que si P1 le polynôme caractéristique de L1 et PJ le polynôme caractéristique
de J, alors
n 1
P1 (λ) = λ 2 PJ (λ 2 ).
32
2.4 Conditionnement
En effet, sachant que D − E est triangulaire supérieure de diagonale celle de la matrice

inversible D, donc elle est aussi inversible et on a
P1 (λ) = det((D − E)−1 F − λIn )

= det(D − E)−1 det(F − λD − λE)
= det D−1 det(−λE + F − λD)
n 1 1 1
= λ 2 det D−1 det(−λ 2 D + λ 2 E + λ− 2 F )
n 1
= λ 2 det D−1 det(−λ− 2 D + E + F )
n 1
= λ 2 det(D−1 (E + F ) − λ 2 In )
n
= λ 2 PJ (λ1/2 ).
En conclusion, λ ∈ Sp(J) si et seulement si λ2 ∈ Sp(L1 ). Ainsi ρ(L1 ) = (ρ(J))2 .
Critère ou test d’arrêt

En général, dans les méthodes itératives pour la détermination de x̄ solution d’un
problème, on construit une suite (x(k) ) qui converge vers x̄. En pratique, et si on n’exige pas
un critère d’arret, le nombre d’itérations pour calculer les termes de cette suite pourrait
être infini. Le calcul donc devrait être interrompu dès qu’on obtient une solution approchée
à ε prés, vérifiant par exemple kAx(k) − bk ≤ ε si x̄ est aussi solution de Ax = b, ou
kx(k+1) −x(k) k ≤ ε dans le cas général, pour une tolérence ε donnée et une norme vectorielle
k · k choisie. Il existe aussi d’autres critères d’arrêt comme
kAx(k) − bk
≤ε
kbk
ou
kx(k) − x(k−1) k
≤ε
kx(k) k
pour les algorithmes de résolution des systèmes linéaires.
2.4 Conditionnement
La notion du conditionnent d’une matrice peut servir à établir des majorations des
erreurs d’arrondi dues au erreurs sur les données. En pratique, par exemple pour un
système linéaire Ax = b, les données A et b ne sont données qu’à une erreur près. Cette
perturbation des données peut engendrer une grande erreur de la solution. C’est le cas de
l’exemple suivant :
Soit le système linéaire Ax = b, pour

1 2 3
A= et b = .
2 4 + 10−6 6 + 10−6
33

1
Ce système admet une solution unique x = .
1
3
Si on change légèrement b par b̃, avec b̃ = , la solution du système Ax =
6 − 10−6
5
b̃ est x̃ = . Ainsi une petite perturbation de la donnée b implique une grande
−1
modification de la solution.
En général, si x est tel que Ax = b et x + δx vérifiant A(x + δx) = b + δb, alors Aδx = δb
ou encore δx = A−1 b. Par conséquent
kδxk ≤ kA−1 kkδbk, (2.4)
pour k · k une norme matricielle subordonnée.

D’autre part kbk = kAxk ≤ kAkkxk et donc
1 kAk
≤ . (2.5)
kxk kbk
De (2.4) et (2.5), on déduit la majoration d’erreur suivante :
kδxk kδbk
≤ kAkkA−1 k . (2.6)
kxk kbk
On remarque que pour kδxk soit assez petit pour une petite perturbation de b, il suffit
que le terme kAkkA−1 k soit aussi assez petit,
Définition 2.4.1 On appelle conditionnement d’une matrice inversible A, pour une norme
matricielle subordonnée k · k, le réel positif, qu’on note Cond(A)
Cond(A) = kAkkA−1 k.
Proposition 2.4.1 Pour toute matrice inversible réelle A on a :

1. Cond(A) = Cond(A−1 ).
2. Cond(αA) = Cond(A), pour tout réel α non nul.
3. Cond(A) ≥ 1.
4. Cond(In ) = 1.
max |λi |
1≤i≤n
5. Cond2 (A) = , si A est une matrice symétrique inversible de valeurs
min |λi |
1≤i≤n
propres (λi )1≤i≤n .
En plus de la majoration (2.6), on a en général
Proposition 2.4.2 Soient A et δA deux matrices d’ordre n, A inversible et soient b et

δb deux vecteurs de Rn , b 6= 0.
34
2.4 Conditionnement
1. Si Ax = b et A(x + δx) = b + δb, alors on a
kδxk kδbk
≤ Cond(A) .
kxk kbk
2. Si Ax = b et (A + δA)(x + δx) = b, alors
kδxk kδAk
≤ Cond(A)
kx + δxk kAk
Preuve.
1. Déja montré
2. Si Ax = b et (A + δA)(x + δx) = b, alors δA(x + δx) = −Aδx et donc δx =
−A−1 δA(x + δx) et on obtient ainsi
kδAk
kδxk ≤ kA−1 kkδAkkx + δxk = Cond(A) kx + δxk.
kAk
1
Il suffit de multiplier par kx+δxk
pour conclure le résultat désiré.
Préconditionnement d’un système

Au lieu de résoudre un système Ax = b, on multiplie à gauche par une matrice
inversible P et résoudre le système P Ax = P b, où la matrice P est choisie pour que P A
soit bien conditionnée. Ce dernier système s’appelle système préconditionné.
35
36
Chapitre 3
Optimisation sans contraintes
3.1 Introduction
On rappelle que tout minimum ou maximum x d’une fonction réelle dérivable f :
R → R satisfait f 0 (x) = 0. Ce résultat s’applique aussi pour les fonctions J : Rn → R
dérivables.
Ce chapitre vise dans sa première partie à rappeler quelques notions élémentaires
concernant les problèmes d’optimisation non contraints sur Rn . Ensuite, et dans la dernière
partie de ce chapitre, des méthodes de descente de type gradient seront développées pour
calculer numériquement un minimum x̄ d’une fonction J sur Rn .
Comme application on considère le problème de minimisation d’une fonction quadra-
tique dont la résolution est équivalente à la résolution d’un système linéaire Ax = b, où
A est une matrice symétrique définie positive.
3.2 Optimisation sur Rn

Dans de nombreux domaines d’applications d’analyse numérique, on est amené à
minimiser une fonction J et chercher x̄ vérifiant
J(x̄) ≤ J(x), ∀ x ∈ Rn , (3.1)
où la fonction dite objectif J est une fonction numérique de variable un vecteur x de Rn .
C’est un problème de minimisation qu’on peut noter aussi
min J(x) (P ).
x∈Rn
Pour maximiser une fonction J il suffit de minimiser (−J), un problème d’optimisation

peut donc s’écrire toujours sous la forme de (P ).
Tout vecteur x̄ solution de (3.1) est appelé solution optimale de (P ). On dit aussi
que x̄ présente (ou tout simplement) un minimum de J sur Rn .
Du théorique au numérique, pour la résolution du problème (P ) ou de (3.1) on s’in-
teresse aux questions suivantes :
37
Chapitre 3. Optimisation sans contraintes
1. Existe-il une solution et est elle unique ?

2. Comment caractériser cette solution ?
3. Comment approcher d’une manière efficace cette solution ?
Commençons par développer la première étape concernant l’existence puis l’unicité
d’une solution d’un problème d’optimisation sans contraintes.
Existence et unicité d’un minimum

On rappelle tout d’abord qu’une suite de vecteurs (x(k) ) est dite bornée de Rn s’il
existe R > 0 tel que
kx(k) k ≤ R, pour tout entier k,
i.e., la suite (x(k) ) est contenue dans une boule de centre 0 et de rayon R > 0.
On rappelle aussi que de toute suite bornée de Rn , on peut en extraire une sous suite
convergente.
Proposition 3.2.1
On suppose que J est continue et qu’elle vérifie
lim J(x) = +∞. (3.2)

kxk→+∞
Alors le problème de minimisation (P ) admet au moins une solution.
Preuve. Soit
α = infn J(x).
x∈R
Montrons que α est atteint, i.e. il existe x̄ ∈ Rn tel que J(x̄) = α.

D’après la propriété de la borne inférieure, il existe une suite (x(k) )k de K
lim J(x(k) ) = α.
p→∞
Alors la suite (x(k) ) est bornée. En effet, si la suite (x(k) ) est non bornée, elle vérifie, pour
une sous suite, (notée aussi (x(k) )), lim ||x(k) || = +∞ et comme lim J(x(k) ) = +∞,
k→+∞ kxk→+∞
donc α = +∞, ce qui est une contradiction.
Comme (x(k) ) est bornée, on peut en extraire une sous suite (x(k) ) convergente vers
x̄.
Comme J est continue, alors
J(x) = lim J(x(k) ) = lim J(x(k) ) = inf J(x).

k→+∞ k→+∞ x∈K
D’où l’existence d’au moins d’une solution pour le problème (P ).
Remarques 3.2.1 1. Si J vérifie (3.2), on dit qu’elle infinie à l’infinie ou qu’elle

est coercive.
38
2. Une fonction J continue et minorée sur Rn , alors elle admet une borne inf, mais
pas de minimum si elle n’est pas infinie à l’infini, voici un contre exemple : La
fonction x 7→ exp(x) est positive, mais elle n’a pas de minimum sur R.
Remarque 3.2.1 Ce résultat donne l’existence d’un minimum mais il ne donne pas l’uni-
cité. On peut avoir des résultats d’unicité si on utilise la convexité.
Définition 3.2.1 Soit J une fonction définie sur Rn à valeurs dans R.

1. On dit que J est convexe sur Rn si :
∀ (x, y) ∈ Rn × Rn , ∀ t ∈ [0, 1], J((1 − t)x) + ty) ≤ (1 − t)J(x) + tJ(y).
2. On dit que J est strictement convexe sur Rn si :
∀ (x, y) ∈ Rn × Rn / x 6= y, ∀ t ∈]0, 1[, J((1 − t)x) + ty) < (1 − t)J(x) + tJ(y).
Remarques 3.2.2 1. Graphiquement, une fonction est convexe si sa courbe se trouve

en dessus de ses cordes.
2. Une fonction J est dite concave si (−J) est convexe.
Exemples 3.2.1
1. L’application x 7→ kxk22 est strictement convexe sur Rn .
2. Toute forme linéaire J(x) = (b, x) = bT x où b ∈ Rn est convexe mais non stricte-
ment convexe.
Proposition 3.2.2 Si la fonction J : Rn −→ R est convexe, alors l’ensemble des solu-

tions du problème (P ) est convexe.
Si de plus J est strictement convexe et elle admet un minimum sur Rn , il est unique.
Preuve. Soit S l’ensemble des solutions optimales. Montrons que S est convexe. Soit
(x1 , x2 ) ∈ S × S, si t ∈ [0, 1], et comme J est convexe donc :
J((1 − t)x1 + tx2 ) ≤ (1 − t)J(x1 ) + tJ(x2 )

≤ (1 − t)J(x) + tJ(x) = J(x), ∀ x ∈ Rn
On en déduit que (1 − t)x1 + tx2 ∈ S, pour tout t ∈ [0, 1]. D’où la convexité de S.
Si de plus J est strictement convexe, on suppose qu’il existe deux solutions x1 et x2
x1 + x2
de (P ) tel que x1 6= x2 . Alors ∈ S et en vertu de la stricte convexité de J il vient
2
que,
x1 + x2
J( ) < 21 J(x1 ) + 12 J(x2 ) = J(x1 ).
2
Contradiction, ainsi x1 = x2 .
39
3.2.1 Conditions d’optimalité

Fonction différentiable
Définition 3.2.2 On dit que J : Rn → R est différentiable au point x ∈ Rn , s’il existe
p ∈ Rn et une fonction : Rn → R tels que
J(x + h) = J(x) + (p, h) + khk(h) et lim (h) = 0.

h→0
Le vecteur p s’il existe est appelé la différentielle (ou la dérivée) de J en x et il est noté
J 0 (x).
On dit que J est différentiable sur Rn si elle est différentiable en tout point de Rn .
Remarque 3.2.2 Si on note (e1 , . . . , en ) la base canonique de Rn , on peut vérifier que

J : Rn → R est différentiable en x, de plus
∂J
 
 ∂x1 (x) 
J 0 (x) = ∇J(x) = 
 .. 
(3.3)
 . 

 ∂J 
(x)
∂xn
∂J J(x + tei ) − J(x)

qui est le gradient de J en x où (x) = lim , i = 1, . . . , n sont les
∂xi t→0 t
dérivées partielles de J.
Exercice 3.2.1 : On pourra montrer en exercice que :

1. Toute forme linéaire J(x) = (b, x) = bT .x où b ∈ Rn , est différentiable et on a
J 0 (x).h = (b, h) = bT h,
donc
J 0 (x) = ∇J(x) = b.
2. La fonctionnelle J : x 7→ J(x) = 21 a(x, x), où a est une forme bilinéaire symétrique
continue sur Rn , alors
J 0 (x).h = a(x, h).
Donc
J 0 (x) = a(x, .).
1
En particulier si A est une matrice symétrique, alors J : x → (Ax, x) est
2
différentiable et on a :
J 0 (x) = ∇J(x) = Ax. (3.4)
40
Fonction différentiable convexe

On peut caractériser la convexité et la convexité stricte des fonctions différentiables.
C’est l’objet de deux propositions suivantes :
Proposition 3.2.3 Soit J : Rn → R une fonction différentiable, alors les propriétés
suivantes sont équivalentes
i) J est convexe sur Rn
ii) pour tout (x, y) ∈ Rn × Rn on a :

J(y) ≥ J(x) + (∇J(x), y − x)
iii) ∇J est un opérateur monotone i.e : ∀ (x, y) ∈ Rn × Rn on a :
(∇J(x) − ∇J(y), x − y) ≥ 0.
Preuve.
i) ⇒ ii) Supposons que J est convexe, alors pour tout t ∈]0, 1] et pour tout (x, y) ∈
K × K on a :
J(tx + (1 − t)y) ≤ (1 − t)J(x) + tJ(y).
Ou encore
J(x + t(y − x)) − J(x)
J(y) ≥ J(x) +
t
En faisant tendre t vers 0 on obtient
J(y) ≥ J(x) + (∇J(x), y − x).
ii) ⇒ iii) Si x, y ∈ Rn , alors
J(y) ≥ J(x) + (∇J(x), y − x) ≥ J(y) + (∇J(y), x − y) + (∇J(x), y − x).
Par conséquent
(∇J(y) − ∇J(x), y − x) ≥ 0.
iii) ⇒ i) Pour montrer que J est convexe il suffit de montrer que, pour tout x, y
dans Rn , la fonction g définie sur [0, 1] par
g(t) = (1 − t)J(x) + tJ(y) − J((1 − t)x + ty)
est positive. Il est clair que g est dérivable sur [0, 1] et que
g 0 (t) = −J(x) + J(y) − (∇J(x + t(y − x)), y − x).
Pour t1 6= t2 on a :
(g 0 (t1 )−g 0 (t2 ))(t1 −t2 ) = (∇J(x+t2 (y −x))−∇J(x+t1 (y −x)), y −x)(t2 −t1 ) ≤ 0.
Donc g 0 est décroissante. Comme g est continue sur [0, 1] et dérivable sur ]0, 1[
vérifiant g(0) = g(1) = 0, d’après le théorème de Rolle, il existe c ∈]0, 1[ tel que
g 0 (c) = 0, ce qui entraı̂ne que g est positive pour tout t ∈ [0, 1].
41
On peut montrer d’une façon analogue que
Proposition 3.2.4 Si J : Rn → R une fonction différentiable, alors les propriétés sui-

vantes sont équivalentes
i) J est strictement convexe sur Rn .
ii) pour tout (x, y) ∈ Rn × Rn tels que x 6= y, on a :
J(y) > J(x) + (∇J(x), y − x) (3.5)
iii) ∇J est un opérateur strictement monotone i.e : ∀ (x, y) ∈ Rn × Rn tels que

x 6= y, on a :
(∇J(x) − ∇J(y), x − y) > 0
Exemple 3.2.1 Si A est une matrice symétrique, la fonction J : Rn → R; x 7→ 21 (Ax, x)−

(b, x) est strictement convexe sur Rn si et seulement si
(∇J(x) − ∇J(y), x − y) = (A(x − y), x − y) > 0 pour tout x, y ∈ Rn , x 6= y.
On rappelle que cette dernière propriété sur A est vérifiée si et seulement si A est définie
positive.
Condition nécessaire et suffisante d’optimalité

En général pour déterminer un minimum d’une fonction J on cherche parmi les points
critiques qui vérifient la condition nécessaire d’optimalité, qui devient aussi suffisante
lorsque J est convexe, suivante.
Proposition 3.2.5 Soit J : Rn → R une fonction différentiable en x̄

1. Si x̄ réalise un minimum de J sur Rn . alors,
∇J(x̄) = 0. (3.6)
2. Si la fonction J est convexe et si ∇J(x̄) = 0, alors x̄ réalise un minimum de J

sur Rn .
La condition (3.6) est connue sous le nom de l’équation d’Euler qui est sans la convexité
de J est en général une condition nécessaire mais non suffisante d’optimalité.
La fonction J : x 7→ x3 vérifie J 0 (0) = 0 mais 0 n’est ni um minimum ni un maximum
de J sur R.
42
3.2.2 Problème d’optimisation quadratique

Si la fonction objectif est de type quadratique, c.-à.-d, si elle est de la forme
J(x) = 12 (Ax, x) − (b, x),
pour A ∈ Mn (R) est symétrique et b ∈ Rn , alors le problème de minimisation (P ) est dit

quadratique et on a :
Proposition 3.2.6 Si A est symétrique définie positive, alors le problème quadratique

(P ) admet une solution unique x̄ vérifiant
Ax̄ = b.
Preuve. La fonction J est continue, infinie à l’infini, puisqu’elle vérifie (1.6) puis d’après
l’inégalité de Cauchy-Schawrz,
λ1 λ1
J(x) ≥ kxk22 − (b, x) ≥ kxk22 − kbk2 kxk2 → ∞ .
2 2 kxk2 →+∞
De plus J est strictement convexe, donc J admet un minimum unique x̄ sur Rn vérifiant
∇J(x̄) = Ax̄ − b = 0.
3.2.3 Problème de moindre carré

Pour A ∈ Mnp (R) et b ∈ Rp , avec n > p, le système Ax = b ne possède pas, en
général, de solution. On peut se contenter alors à chercher x qui minimise le carré de la
norme de Ax − b et donc ainsi formuler le problème de moindres carrés suivant :
min 1 kAx − bk22 . (3.7)

x∈Rn 2
Remarque 3.2.3 Pour chercher la solution de (3.7) et afin de pouvoir appliquer l’équation
d’Euler, on minimise la norme au carré de Ax−b au lieu de kAx−bk2 puisque l’application
norme n’est jamais différentiable en 0. Clairement, les deux problèmes sont équivalents.
Proposition 3.2.7 Le problème (3.7) admet au moins une solution. De plus toute solu-
tion x̄ de (3.7) est aussi solution de système
AT Ax = AT b.
Le problème (3.7) admet une solution unique si et seulement si la matrice A est injective.
Preuve.
La fonction objectif ici est
J(x) = 12 (AT Ax, x) − (AT b, x) + 12 kbk22 .
43
La matrice AT A est toujours symétrique semi définie positive puisque (AT Ax, x) =
kAxk22 ≥ 0 pour tout vecteur x de Rn . La fonction quadratique J est donc convexe.
Par conséquent, x̄ est solution de (3.7) si et seulement
∇J(x̄) = AT Ax̄ − AT b = 0.
Utilisant le résultat classique d’algèbre linéaire suivant : Im(AT ) = Im(AT A), puisque les
deux espaces sont de même dimension et Im(AT A) ⊂ Im(AT ), on déduit que le système
AT Ax = AT b admet au moins une solution, d’où l’existence d’au moins d’un minimum
de J sur Rn . L’unicité de minimum est assurée si et seulement si le système AT Ax = AT b
admet une solution unique, pour tout vecteur b, donc si et seulement si la matrice AT A
est bijective. Ceci étant vérifié si et seulement A est injective.
3.3 Algorithmes de descentes

Les méthodes numériques de descentes en général sont des méthodes itératives dont
le but est de déterminer x̄ réalisant le minimum d’une fonction J sur Rn , en utilisant
des directions de déplacement permettant de se rapprocher le plus possible de x̄. Dans ce
chapitre on va utiliser des méthodes numériques de descente du gradient pour un problème
de minimisation sans contrainte quelconque dont la convergence sera étudiée seulement
pour le problème quadratique.
1
J(x̄) = minn (Ax, x) − (b, x), (3.8)
x∈R 2
où A est une matrice symétrique définie positive et b un vecteur de Rn .

Ces méthodes s’appliquent alors aussi pour résoudre numériquement un système
linéaire Ax = b, pour A une matrice symétrique définie positive.
Définition 3.3.1 Soit J une fonction de Rn à valeurs dans R. Soit x ∈ Rn . On dit que
d ∈ Rn , avec d 6= 0, est une direction de descente de J en x s’il existe α0 > 0 tel que
J(x + αd) ≤ J(x), ∀ α ∈ [0, α0 ]
Ainsi, une méthode de descente pour la recherche de x̄ solution de :
min J(x)
x∈Rn
consiste à construire une suite (x(k) )k∈IN de la manière suivante :
44
3.4 Méthodes des gradients
Initialisation : x(0) ∈ Rn ,
Pour k ≥ 0 :
On cherche la direction de descente d(k) au point x(k) et on détermine le pas αk ,
Ensuite, on fait la mise à jour : x(k+1) = x(k) + αk d(k) .
Proposition 3.3.1 Soient J : Rn → R une fonction différentiable, x, d ∈ Rn avec d 6= 0.

Si d est une direction de descente en x, alors :
(∇J(x), d) ≤ 0.
Preuve. Soit la fonction ϕ de R dans R définie par : ϕ(α) = J(x + αd). Alors ϕ est
dérivable
ϕ0 (α) = (∇J(x + αd), d)
On considère d ∈ Rn une direction de descente au point x alors par définition, il existe
α0 > 0 tel que :
J(x + αd) ≤ J(x), ∀ α ∈ [0, α0 ].
Comme d est une direction de descente, on peut écrire, pour tout α ∈ [0, α0 ],
ϕ(α) ≤ ϕ(0),
et donc pour tout α ∈]0, α0 ]

ϕ(α) − ϕ(0)
≤ 0.
α−0
En passant à la limite lorsque α tend vers 0, on déduit que ϕ0 (0) ≤ 0, ce qui signifie que
(∇J(x), d) ≤ 0.

Dans ces méthodes les directions de descentes s’expriment en fonction du gradient de
la fonction à minimiser.
Remarque 3.4.1 Ces méthodes itératives de descente ne sont pas aussi finies. Il faut
donc un critère d’arrêt qui permet d’arréter les itérations dès que l’itére x(k) s’approche
de la solution x̄ du problème à minimiser. Comme test d’arrêt classique pour ce type
d’algorithmes consiste à s’arrêter à l’itération k si
k∇J(x(k) )k ≤ ε ou si kx(k+1) − x(k) k ≤ ε,
pour ε une précision donnée.
45
3.4.1 Méthode du gradient à pas fixe :

La méthode du gradient à pas fixe α, s’écrit pour α > 0, comme suit :
1. On choisit x(0) ∈ Rn
2. Pour k ≥ 0, on calcule :
d(k) = −∇J(x(k) ) = b − Ax(k) ,
x(k+1) = x(k) + αd(k)
Proposition 3.4.1 : Convergence de la méthode du gradient à pas fixe pour un

problème quadratique
Si J(x) = 21 (Ax, x) − (b, x), où A est une matrice symétrique définie positive, alors,
2
la méthode du gradient à pas fixe converge si et seulement si α ∈]0, [.
ρ(A)
Preuve. On a, pour k ≥ 0,
x(k+1) = x(k) − α(Ax(k) − b) = (I − αA)x(k) + αb.
La matrice d’itération est Bα = I − αA dont le spectre est
Sp(Bα ) = {1 − αλ; λ ∈ Sp(A)} .
Si on note λ1 et λn la plus petite et la plus grande valeur propre de A, alors
1 − αλn ≤ 1 − αλ ≤ 1 − αλ1 , ∀λ ∈ Sp(A).
Ainsi,
ρ(Bα ) = max(|1 − αλ1 |, |1 − αλn |)
Par conséquent, la méthode du gradient à pas fixe converge si et seulement ρ(Bα ) < 1.
Donc si et seulement si 0 < α < λ21 et 0 < α < λ21 . Cette dernière condition est bien
équivalente à
2
0<α< .
ρ(A)
3.4.2 Méthodes du gradient à pas optimal

Méthodes de descente à pas optimal
Une méthode de descente en général de type x(k+1) = x(k) + αk d(k) , pour x(0) donné et
pour une direction d(k) connue, est dite à pas optimal si l’on choisit le pas αk de manière
46
à minimiser la fonction α 7→ J(x(k) + αk d(k) ). Calculons le pas αk pour un problème

quadratique. Ceci revient à chercher α ∈ R vérifie, pour x, d ∈ Rn tel que d 6= 0 :
∀ r ∈ R, J(x + αd) ≤ J(x + rd) (Q).
Lemme 3.4.1 Si la matrice A est symétrique définie positive, alors (Q) admet une unique
solution donnée par
(Ax − b, d)
α=− .
(Ad, d)
Preuve. On considère la fonction f : R → R définie par :
f (r) = J(x + rd)
On a
f (r) = J(x + rd)

1
= (A(x + rd), x + rd) − (b, x + rd)
2
1
= ((Ax, x) + r(Ax, d) + r(Ad, x) + r2 (Ad, d)) − (b, x) − r(b, x)
2
1 2 1
= r (Ad, d) + r(Ax − b, d) + (Ax, x) − (b, x).
2 2
Puisque (Ad, d) > 0 (car A est définie positive et d 6= 0), on en déduit que la fonction
polynôme de degré 2 f admet un minimum global au point r = α donné par :
(Ax − b, d)
α=− .
(Ad, d)
Remarque 3.4.2 Pour un problème quadratique de type (3.8), une méthode de descente
à pas optimal αk donne une relation d’orthogonalité entre la direction de descente d(k) et
g (k+1) = ∇J(x(k+1) ) = Ax(k+1) − b = −d(k+1) .
En effet
(g (k+1) , d(k) ) = (Ax(k+1) − b, d(k) ) = (Ax(k) − b, d(k) ) + αk (Ad(k) , d(k) ) = 0. (3.9)
47
Algorithme du gradient à pas optimal

La méthode du gradient à pas optimal est une méthode de descente à pas optimal
dont la direction de descente est donnée par le gradient.
Le principe de cette méthode s’écrit :
Initialisation : x(0) ∈ Rn donné.

Pour k ≥ 0,
On calcule d(k) = −∇J(x(k) ),
On choisit αk ≥ 0 tel que : J(x(k) + αk d(k) ) ≤ J(x(k) + αd(k) ), pour tout α ≥ 0.
On pose x(k+1) = x(k) + αk d(k) .
Ainsi, l’algorithme du gradient à pas optimal pour un problème quadratique s’écrit :
1. Initialisation : x(0) ∈ Rn donné,

2. Pour k ≥ 0, on calcule :
 (k)
d = b − Ax(k) ,
(k) 2
 αk = k d k 2 ,

 hAd(k) , d(k) i
(k+1)
x = x(k) + αk d(k) .
Proposition 3.4.2 Si A est une matrice symétrique définie positive, alors la méthode du
gradient à pas optimal converge.
Preuve. Par construction la suite (J(x(k )) est décroissante et on sait qu’elle est minorée
par J(x̄), donc elle est convergente. Par conséquent
lim J(x(k+1) ) − J(x(k) ) = 0.
k→+∞
Or
αk2
J(x(k+1) ) − J(x(k) ) = (Ad(k) , d(k) ) + αk (Ax(k) − b, d(k) ).
2
Sachant que
k d(k) k22
αk = et d(k) = −g (k) = −(Ax(k) − b)
hAd(k) , d(k) i
donc
kg (k) k42
J(x(k+1) ) − J(x(k) ) = − 21 .
(Ag (k) , g (k) )
48
De (1.6) on déduit
1 (k) 2
kg k2 ≤ J(x(k+1) ) − J(x(k) ) → 0 ,
λn k→+∞
où λn est la plus grande valeur propre de A.

Par conséquent
kx(k) − x̄k2 = kA−1 (A(x(k) − b)k2 ≤ kA−1 k2 kg (k) k2 → 0 .

k→+∞
3.4.3 Méthodes du gradient conjugué

Cette méthode s’applique directement pour chercher la solution unique x̄ du problème
min J(x),
x∈Rn
pour
J(x) = 12 (Ax, x) − (b, x) (P )
où A est une matrice symétrique définie positive.
Dans les algorithmes de deux méthodes de descentes précédentes, les termes de la
suite (x(k)) ) s’approchent de plus en plus vers la solution cherchée x̄ sans généralement
jamais l’atteindre. D’où l’idée de la méthode du gradient conjugué où, comme on le verra,
la suite (x(k) ) générée par l’algorithme de cette méthode est stationnaire, et elle vaut x̄
en moins de n itérations.
Définition 3.4.1 Soit A est une matrice symétrique définie positive de Mn (R).
1. Deux vecteurs non nuls d(1) , d(2) de Rn sont dits A−conjugués si (Ad(1) , d(2) ) = 0
2. Une famille (d(1) , d(2) , . . . , d(p) ) de p vecteurs non nuls de Rn est dite A−conjuguée
si
(Ad(i) , d(j) ) = 0 ∀ i, j ∈ {1, . . . , p}, i 6= j.
Remarque 3.4.3 On peut facilement vérifier que

1. Toute famille A−conjuguée de p vecteurs non nuls de Rn est libre.
2. Une famille A−conjuguée de n vecteurs non nuls de Rn est une base de Rn .
L’idée de la méthode du gradient conjugué, qui est aussi une méthode de descente à pas
optimal, est de construire la direction de descente en fonction du gradient et de la direction
précédente de descente.
On part de x(0) donné et d(0) = −∇J(x(0) ) = −Ax(0) + b = g (0) , Connaissant x(k) ,
d(k−1) et g (k) , et si g (k) 6= 0, on choisit la direction de descente
d(k) = −g (k) + βk−1 d(k−1) (3.10)
49
et
x(k+1) = x(k) + αk d(k) (3.11)
de manière à minimiser la fonction à deux variables
f : (α, β) 7→ J(x(k) + α(g (k) + βd(k−1) ) = J(x(k) + αd(k) ).
Donc le pas ici αk est optimal et on a la relation d’orthogonalité (voir lemme (3.4.1))
suivante
(d(k−1) , g (k) ) = 0. (3.12)
Si g (k) 6= 0 d’après (3.10), on aura aussi d(k) 6= 0 et par conséquent
(Ax(k) − b, d(k) ) (g (k) , −g (k) + βk−1 d(k) ) kg (k) k22

αk = = − = . (3.13)
(Ad(k) , d(k) ) (Ad(k) , d(k) ) (Ad(k) , d(k) )
Il suffit donc de déterminer la direction de descente d(k) . Ceci revient à chercher βk−1 .
Sachant que cette dernière doit vérifier
∂f
(αk , βk−1 ) = αk (∇J(x(k) +αk d(k) ), d(k−1) ) = α(∇J(x(k) +αk (g (k) +βk−1 d(k−1) )), d(k−1) ) = 0.
∂β
(3.14)
et que ∇J(x) = Ax − b, on déduit d’abord puis que
(Ad(k−1) , g (k) )
βk−1 = ,
(Ad(k−1) , d(k−1) )
puis que les directions d(k) et d(k−1) sont A−conjugués,
(Ad(k) , d(k−1) ) = 0. (3.15)
g (k) − g (k−1)
En remplaçant Ad(k−1) par , on déduit que
αk−1
kg (k) k22 kg (k−1) k22

(Ad(k−1) , g (k) ) = et que (Ad(k−1) , d(k−1) ) = .
αk−1 αk−1
On tire alors que

kg (k) k22
βk−1 = .
kg (k−1) k22
D’où l’algorithme du gradient conjugué :
50
1. On choisit x(0) ∈ Rn . On prend g (0) = Ax(0) − b = −d(0) .

2. Pour k ≥ 0, on calcule
kg (k) k22

α k = (k) , d(k) )
, x(k+1) = x(k) + αk d(k) ,
 (k+1) (Ad (k+1)

 g
 = Ax − b = g (k) + αk Ad(k) ,
 kg (k+1) k22
βk = , d(k+1) = −g (k+1) + βk d(k) .
kg (k) k22
On montre alors :
Proposition 3.4.3 Si A est symétrique définie positive, alors l’algorithme du gradient
conjugué converge en au plus n itérations vers le minimum de J(x) := 21 (Ax, x) − (b, x)
sur Rn .
Preuve.
Pour k = 0, si g (0) = 0, alors Ax(0) = b et x̄ = x(0) . Sinon,
d(0) = −g (0) = −Ax(0) + b = −∇J(x(0) ).
(g (0) , d(0) )
α0 = − , x(1) = x(0) + α0 d(0) , g (1) = Ax(1) − b
(Ad(0) , d(0) )
vérifiant
(g (1) , g (0) ) = (g (1) , d(0) ) = (d(1) , Ad(0) ) = 0.
Supposons maintenant, pour 1 ≤ k < n, l’hypothèse de récurrence
(g (k) , g (j) ) = 0 pour 1 ≤ j < k,

(g (k) , d(j) ) = 0 pour 1 ≤ j < k, (3.16)
(k) (j)
(d , Ad ) = 0 pour 1 ≤ j < k,
Si g (k) = 0, alors x̄ = x(k) et l’algorithme converge à l’itération k. Si g (k) 6= 0, on construit

x(k+1) , g (k+1) et d(k+1) par l’algorithme du gradient conjugué.
Vérifions l’hypothèse de récurrence à l’ordre (k + 1)
(g (k+1) , g (j) ) = 0 pour 1 ≤ j < k + 1,

(k+1) (j)
(g ,d ) = 0 pour 1 ≤ j < k + 1, (3.17)
(d(k+1) , Ad(j) ) = 0 pour 1 ≤ j < k + 1,
D’après (3.12), et par construction on a
(g (k+1) , d(k) ) = 0,
51
et pour j < k, on a
(g (k+1) , d(j) ) = (g (k+1) , d(j) ) − (g (k) , d(j) ) = (g (k+1) − g (k) , d(j) ) = αk (Ad(k) , d(j) ) = 0
=0
f
On a
g (k) = −d(k) + βk−1 d(k−1)
donc, pour j ≤ k, on a
(g (k+1) , g (j) ) = (g (k+1) , d(j) ) − βj (g (k+1) , d(j−1) ) = 0.
Or, pour j ≤ k, g (j+1) = g (j) + αj Ad(j) donc
αj (g (k+1) , Ad(j) ) = (g (k+1) , g (j+1) ) − (g (k+1) , g (j) ) = 0.
Sachant que g (j) 6= 0, donc αj 6= 0 et par conséquent (g (k+1) , Ad(j) ) = 0.

Comme
d(k+1) = −g (k+1) + βk d(k) ,
alors, pour j ≤ k, on a
(d(k+1) , Ad(j) ) = (g (k+1) , Ad(j) ) + βk (d(k) , Ad(j) ) = 0
et la récurrence est vérifiée.

Si g (k) 6= 0, pour k = 0, . . . , n−1, alors (d(0) , d(1) , . . . , d(n−1) ) est une famille A−conjuguée
de n vecteurs, donc base de Rn et le vecteur g (n) sera orthogonal à d(0) , d(1) , . . . , d(n−1) .
Nécessairement g (n) = 0 et par suite x(n) = x̄. L’algorithme ainsi converge à l’itération
exactement n.
Remarque 3.4.4 Les directions de descentes sont telles que la base (d(0) , d(1) , . . . , d(n−1) )
est obtenue par le procédé d’orthogonalité de Gram-Shmidt, adapté au produit scalaire
(Ax, y), à partir de la base (−g (0) , . . . , −g (n−1)) ).
3.4.4 Vitesse de convergence

En général une suite (x(k) ) qui converge vers x̄, vérifiant x(k) 6= x̄ pour tout entier k,
est dite d’ordre de convergence au moins p ∈ R∗+ si
kx(k+1) − x̄k = O(kx(k) − x̄kp ).
En particulier, s’il existe γ ∈]0, 1[ tel que
kx(k+1) − x̄k
lim = γ,
k→+∞ kx(k) − x̄kp
52
alors la convergence est dite exactement d’ordre p. La convergence est dite linéaire si
p = 1, quadratique si p = 2.
Les méthodes de descente du gradient à pas fixe et à pas optimal sont à convergence au
moins linéaire et on pourra montrer que
cond2 (A) − 1 (k)

kx(k+1) − x̄k ≤ kx − x̄k, (3.18)
cond2 (A) + 1
où k.k désigne la norme euclidienne pour la méthode du gradient à pas fixe et la norme
1
définie par kxk = kxkA = (Ax, x) 2 pour la méthode du gradient à pas optimal.
Pour la méthode du gradient conjugué, et tant que g (k) = Ax(k) − b 6= 0, on peut montrer
la majoration suivante :
p
(k+1) cond2 (A) − 1 k (0)
kx − x̄kA ≤ 2( p ) kx − x̄kA . (3.19)
cond2 (A) + 1
3.4.5 Méthodes des gradients et préconditionnement

Des majorations (3.18) et (3.19), on constate que l’erreur x(k) − x̄, obtenue au cours
des méthodes des gradients, diminue d’autant plus rapidement que cond2 (A) est petit.
On rappelle aussi que les méthodes des gradients sont aussi des méthodes numériques
pour résoudre un système linéaire Ax = b, pour A ∈ Mn (R) symétrique définie positive.
Ainsi, la matrice du système préconditionné doit aussi être symétrique définie positive. Le
préconditionnement de ce système consiste alors à appliquer les méthodes des gradients
pour le problème de minimisation
min 1 (Ãy, y) − (b̃, y) (Pc ),

y∈Rn 2
où P est une matrice inversible telle que la matrice Ã = P AP T soit bien conditionnée et
pour b̃ = P b. La nouvelle matrice Ã est aussi symétrique définie positive et la solution x̄
du système Ax = b est donnée par x̄ = P T ȳ, pour ȳ solution optimale de (Pc ).
53
54
Chapitre 4
Optimisation avec contraintes

linéaires
Lorsque un minimum x̄ d’une fonction J : Rn → R doit satisfaire certaines conditions,

on parle de problème de minimisation contraint ou avec contraintes. Dans ce cas on cherche
à calculer le minimum d’une fonction J non pas sur Rn , mais sur un ensemble fermé K
de Rn et notre problème sera de type
min J(x) (P ).
x∈K
Ce chapitre comporte un rappel des résultats d’existence et d’unicité d’une solution

du problème (P ), des conditions d’optimalités puis il donne quelques algorithmes pour la
résolution ce problème.
4.1 Problèmes d’optimisations sous contraintes

Un exemple simple de ce type de problèmes est de déterminer dans le plan le point
M (x, y) de R2 le plus proche possible à un point donné (x0 , y0 ) et qui appartient à une
droite donnée d’équation ax + by = c, donc chercher le couple (x, y) solution de
min (x − x0 )2 + (y − y0 )2 .
ax+by=c
Ici l’ensemble K est la droite d’équation ax + by + c = 0.

On appelle problème de minimisation sous contraintes d’égalités linéaires si
K := {x ∈ Rn tel que Cx = d} ,
où C est une matrice de Mm,n (R) et d ∈ Rm .
Pour les problèmes avec contraintes d’inégalités linéaires si
K := {x ∈ Rn tel que Bx ≤ c} ,
où B est une matrice de Mp,n (R) et c ∈ Rp . On parle de problème de minimisation linéaire
si les contraintes sont d’égalités ou d’inégalités linéaires.
55
Chapitre 4. Optimisation avec contraintes linéaires
4.1.1 Existence et unicité de minimum

Définition 4.1.1 Une partie K de Rn est dite convexe si :
tx + (1 − t)y ∈ K ∀ (x, y) ∈ K × K, ∀ t ∈ [0, 1].
Autrement dit K est une partie convexe si elle contient tout segment d’extrémités deux
quelconques de ses points.
Remarques 4.1.1
1. La notion de convexité et de stricte convexité d’une fonction J sur une partie
convexe K, ainsi que leurs caratérisation lorsque J est différentiable, sont iden-
tiques à celles données dans la définition (3.2.1) et dans les propositions (3.2.3)
et (3.2.4) qui s’appliquent pour tout x, y ∈ K au lieu de Rn .
2. Lorsque K est convexe et la fonction à minimiser J est convexe, le problème de
minimisation (P ) est dit convexe.
Comme pour le problème sans contraintes, on a :

Proposition 4.1.1
Soient K un fermé non vide de Rn et J : K → R continue. On suppose que J est infinie
à l’infini ou que K est borné, alors le problème de minimisation (P ) admet au moins
unesolution. Si de plus K est convexe et J est strictement convexe sur K, ce minimum
est unique.
Preuve. La démonstration est similaire au cas sans contraintes, où on prend pour l’exis-
tence une suite minimisante (x(k) ) qui est aussi bornée. Comme K est fermé, donc x̄ la
limite d’une sous suite de (x(k) ) appartient aussi à K. Lorsque le problème est stricte-
ment convexe et il admet un minimum, et comme dans le cas non contraint, l’unicité est
immédiate.
Projection sur un convexe

Soit K une partie non vide de Rn et x un vecteur de Rn qui n’appartient pas à K.
On cherche à définir la distance de ce vecteur x à K, c.à.d, le réel qui réalise la distance
minimale de x à tout les points de K. Cette distance est elle finie ? et existe-il x̄ ∈ K
qui réalise cette distante ? La réponse à ces deux questions n’est pas évidente à priori.
Cependant et lorsque K est convexe et fermé on a :
Proposition 4.1.2 Soit K une partie convexe fermée non vide de Rn et x un point de
Rn . Alors il existe un unique point de K, noté PK (x) tel que :
(
PK (x) ∈ K
(4.1)
kx − PK (x)k2 ≤ ky − xk2 , ∀ y ∈ K
PK (x) est appelé projection du point x sur le convexe fermé K.
56
Preuve. On pose, pour y ∈ K

1 1
J(y) = ky − xk22 = (y − x, y − x).
2 2
La fonction J est strictement convexe et différentiable avec ∇J(y) = y − x. D’après le
théorème d’existence et d’unicité, le problème admet une unique solution.
4.1.2 Condition d’optimalité

Soit x̄ une solution optimale de problème simple suivant :
min ax,
x∈[−1,1]
pour a ∈ R. Si a ≥ 0, alors la fonction J : x 7→ ax atteint son minimum en x̄ = −1

vérifiant donc J 0 (x̄) ≤ 0. Si a < 0, alors x̄ = 1 et dans ce cas J 0 (x̄) ≥ 0. Dans les deux
cas, x̄ vérifie
J 0 (x̄)(x − x̄) ≥ 0, ∀ x ∈ [−1, 1]. (4.2)
Dans la proposition qui suit, on généralise la condition d’optimalité (4.2) pour un ensemble
de contraintes convexe K quelconque.
Proposition 4.1.3 Soit J : Rn → R une fonction différentiable en x̄ et soit K un convexe
fermé de Rn . Si x̄ réalise un minimum de J sur K, alors
(∇J(x̄), x − x̄) ≥ 0 ∀ x ∈ K. (4.3)
Si de plus la fonction J est convexe sur K, alors x̄ réalise un minimum de J sur K, si et
seulement si
(∇J(x̄), x − x̄) ≥ 0 ∀ x ∈ K. (4.4)
L’inéquation (4.3) et connue sous le nom d’inéquation d’Euler.
Preuve.
1. Soit x ∈ K, on a alors pour 0 < t < 1, x̄ + t(x − x̄) = tx + (1 − t)x̄ ∈ K, puisque
K est convexe. Comme x̄ est un minimum de J sur K, alors
J(x̄ + t(x − x̄)) − J(x̄) ≥ 0.
Par conséquent
J(x̄ + t(x − x̄) − J(x̄)
lim+ = (∇J(x̄), x − x̄) ≥ 0.
t→0 t
2. On suppose que J est convexe sur K et soit x̄ ∈ K vérifiant l’inéquation d’Euler.
De la convexité de J on déduit que
J(x) ≥ J(x̄) + (∇J(x̄), x − x̄)
≥ J(x̄), ∀ x ∈ K
Par suite x̄ est une solution de (P ).
57
Exercice 4.1.1 Montrer que la fonction projection sur le convexe fermé K est caractérisée
par
(x − PK (x), y − PK (x)) ≤ 0 ∀ y ∈ K. (4.5)
et qu’elle vérifie
kPK (x) − PK (y)k2 ≤ kx − yk2 ∀ x, y ∈ Rn . (4.6)
Conditions d’optimalités pour les contraintes d’égalité linéaires

On suppose que J est différentiable en x̄ minimum de J sur K. Commençons par le
cas de contraintes d’égalités où on considère le problème
min J(x), (4.7)

Cx=d
pour C ∈ Mm,n (R) et d ∈ Rm .
Proposition 4.1.4 On suppose que J est différentiable en x̄. Si le vecteur x̄ est solution
de (4.7), il existe alors λ̄ ∈ Rm tel que
∇J(x̄) + C T λ̄ = 0. (4.8)
Le multiplicateur λ̄ est unique si la matrice C est de rang m. (rg (C) = m).
Preuve. Il est clair que K est un convexe fermé de Rn . Comme J est différentiable en x̄
qui est un minimum de J sur K, d’après la condition d’optimalité (4.1.3), x̄ vérfie
(∇J(x̄), x − x̄) ≥ 0 ∀ x ∈ K
Or si x ∈ K, alors Cx = d = C x̄ et par conséquent C(x − x̄) = 0. Ainsi x − x̄ ∈ ker C où

ker C est le noyau de C.
Inversement, si y ∈ ker C, alors y = x − x̄ pour x ∈ K. On a alors
(∇J(x̄), y) ≥ 0 ∀ y ∈ ker C.
Comme ker C est un sous espace vectoriel, on en déduit que
(∇J(x̄), y) = 0 ∀ y ∈ ker C.
Ceci est équivalent à ∇J(x̄) ∈ (ker C)⊥ . Or, on sait que (ker C)⊥ = ImC T , donc il existe
λ̄ ∈ Rm tel que
∇J(x̄) + C T λ̄ = 0.
L’unicité de ce multiplicateur λ̄ est évidente si C est surjective (de rang m) puisque cette
dernière propriété est équivalente à ker C T = 0.
58
Remarque 4.1.1 Si C = (cij ) ∈ Mm,n (R) et d = (d1 , . . . , dm )T ∈ Rm , la contrainte

Cx = d implique que x est solution d’un système linéaire de m équations et à n inconnues,
s’écrivant donc de la forme

 g1 (x) := c11 x1 + · · · + c1n xn − d1 = 0

.. .. ..
 . . .
g (x) := c x + · · · + c x − d
m m1 1 mn n m = 0
où Ci = (ci1 , . . . , cin )T , i = 1, . . . , m représentent les lignes de la matrice C. Dans ce cas

Xm
et comme C T λ̄ = λi Ci , alors, la condition d’optimalité s’écrit
i=1
m
X m
X
∇J(x̄) + λi Ci = ∇J(x̄) + λi ∇gi (x̄) = 0.
i=1 i=1
4.1.3 Conditions d’optimalités pour les contraintes d’inégalité

linéaires
Soit le problème
min J(x), (4.9)
Bx≤c
Pour B ∈ Mnp (R) et c ∈ Rp .

On adopte ici les notations suivantes :
∀ µ = (µ1 , . . . , µp ) ∈ Rp , µ ≥ 0 si et seulement si µj ≥ 0, ∀j = 1, . . . , p.
Proposition 4.1.5 Si x̄ réalise un minimum de J sur K = {x ∈ Rn , tel que Bx = c et

J est différentiable en x̄, il existe alors µ̄ ∈ Rp tel que
µ̄ ≥ 0, ∇J(x̄) + B T µ̄ = 0 et µ̄j (B x̄ − c)j = 0, pour tout 1 ≤ j ≤ p. (4.10)
Ce vecteur µ̄ est appelé aussi multiplicateur de Lagrange et il est unique si la matrice

B est de rang p (rg(B) = p).
Preuve. On se limitera dans la démonstration pour une matrice B de rang p où pour
tout vecteur y ∈ Rp il existe x ∈ Rn tel que y = Bx.
Soit d = B x̄. Alors d ≤ c et tout vecteur x vérifiant Bx = d est un vecteur de K.
Clairement, x̄ est aussi solution du problème à contraintes égalités linéaires
min J(x). (4.11)

Bx=d
Il existe donc µ̄ ∈ Rp tel que

∇J(x̄) + B T µ̄ = 0.
59
Montrons que µ̄j ≥ 0 et que µ̄j (B x̄ − d)j = 0, pour tout 1 ≤ j ≤ p.

Pour ε > 0 on a dj − ε ≤ dj ≤ cj et pour j ∈ {1, . . . , p}, si x est tel que Bx = d − εej ,
pour ej j-ième vecteur de la base canonique de Rp , alors x est dans K.
D’après l’inéquation d’Euler, on a
(∇J(x̄), x − x̄) ≤ 0.
En remplaçant ∇J(x̄) par −B T µ̄, on obtient
−(µ̄, Bx − B x̄) = εµ̄j ≥ 0.
Par conséquent µ̄j ≥ 0.

Soit maintenant j ∈ {1, . . . , p} tel que (B x̄ − c)j < 0. Alors si ε > 0 suffisamment petit,
tel que dj + ε ≤ cj et si y solution de By = d + εej , alors y ∈ K et on a
(µ̄, By − B x̄) = −εµ̄j ≥ 0.
Le réel positif µj est aussi négatif, donc nul.
4.1.4 Cas de contraintes d’égalités et d’inégalités linéaires

Soit le problème
min J(x)
x∈K
où
K := {x ∈ Rn tel que Cx = d et Bx ≤ c},
pour
C ∈ Mm,n (R), d ∈ Rm , B ∈ Mp,n (R), c ∈ Rp .
Proposition 4.1.6 Si x̄ réalise un minimum de J sur K et J est différentiable en x̄, il

existe (λ̄, µ̄)T ∈ Rm × Rp ,tel que
∇J(x̄) + C T λ̄ + B T µ̄ = 0 et µ̄I ≥ 0, µj (B x̄ − c)j = 0, j = 1, . . . , n.
Preuve. Il suffit d’appliquer le cas avec contraintes d’inégalités linéaires en remarquant

que la contrainte d’égalité Cx = d est équivalente à Cx ≤ d et −Cx ≤ −d.
4.1.5 Problème quadratique

Problème quadratique sans contraintes
Soit A une matrice symétrique. On considère le problème
1
minn (Ax, x) − (b, x) (4.12)
x∈R 2
60
Proposition 4.1.7 Si A est symétrique semi définie positive, alors le problème (4.12)
admet une solution x̄ si et seulement si
Ax̄ = b.
Preuve.
La fonction J est continue, infinie à l’infini, puisqu’elle vérifie, d’après (1.6) puis l’inégalité
de Cauchy-Schawrz,
λ1 λ1
J(x) ≥ kxk22 − (b, x) ≥ kxk22 − kbk2 kxk2 →0 .
2 2 kxk2 →+∞
De plus J est strictement convexe, donc J admet un minimum unique x̄ sur Rn vérifiant
∇J(x̄) = Ax̄ − b = 0.
Problème quadratique avec contraintes d’égalité

Soit le problème
1
min (Ax, x) − (b, x), (4.13)
2 Cx=d
pour A ∈ Mn (R) une matrice symétrique, C ∈ Mm,n (R), et b ∈ Rn et d ∈ Rm .
Proposition 4.1.8 Si A est symétrique définie positive, alors le problème (4.13) admet
une solution unique x̄ caractérisée par :

m A CT x̄ b
∃λ̄ ∈ R tel que = . (4.14)
C 0 λ̄ d
Preuve.
Il s’agit d’un problème de minimisation dont la fonction objectif est infinie à l’infini,
strictement convexe et différentiable, donc admettant un minimum unique x̄. Comme il
s’agit d’un problème de minimisation avec contraintes d’égalités linéaires et d’après (4.8),
alors il existe un multiplicateur de Lagrange λ̄ ∈ Rm tel que
∇J(x̄) + C T λ̄ = Ax̄ − b + C T λ̄ = 0.
Le problème (4.13) est un problème convexe, donc x̄ solution de (4.13) si et seulement si

il existe λ̄ ∈ Rp tel que
Ax̄ + C T λ̄ = b
C x̄ = d
vérifiant donc
A CT x̄ b
= .
C 0 λ̄ d
61
Remarque 4.1.2 Si A est symétrique semi définie positive, alors le problème (4.13) peut
ne pas admettre de solution comme il peut avoir une ou plusieurs solutions. Comme il
s’agit d’un problème convexe, alors x̄ est une solution de (4.13) si et seulement si elle
vérifie (4.14).
Exemple 4.1.1 Soit le problème

min J(x) (4.15)
x∈K
pour
J(x) = 12 x21 + 12 x22 + 12 x23 − x1 x2 − x1 − 2x2 − x3
et
K = {(x1 , x2 , x3 )T ∈ Rn tel que x1 + x2 − x3 = 1, 2x3 + x2 = 2}.
Il s’agit d’un problème quadratique de type (4.13), pour

   
1 −1 0 1
A=  −1 1 0  et b =  2  .
0 0 1 1
La matrice
 A
ici est semi-définie
  positive, (ses valeurs propres sont 0, 1 et 2). Si on note
1 0
c1 =  1  et c2 =  1  les deux vecteurs représentant les deux lignes de la matrice
−1 2
C, alors
C T λ̄ = λ̄1 c1 + λ̄2 c2 .
 
x1
un vecteur x̄ =  x2  est solution de ce problème si et seulement si il existe λ̄ =
x3
T 2
(λ1 , λ2 ) ∈ R tel que


 x1 − x2 + λ1 =1
 −x1 + x2 + λ1 + λ2 = 2


x3 − λ1 + 2λ2 = 1
x1 + x2 − x3 =2




x2 + 2x3 =1

   
1 1
1
Ce système admet comme solution x̄ =  1  et λ̄ = . Donc x̄ =  1  est
1
0 0
la solution cherchée.
62
4.2 Quelques algorithmes
Problème quadratique avec contraintes d’inégalité

On considère le problème
1
min (Ax, x) − (b, x), (4.16)
Bx≤c 2
pour A ∈ Mn (R) une matrice symétrique, B ∈ Mp,n (R), et b ∈ Rn et c ∈ Rp .
Proposition 4.1.9 Si A est symétrique définie positive, alors le problème (4.16) admet
une solution unique x̄ si et seulement si

p A BT x̄ b
∃ µ̄ ∈ R µ̄ ≥ 0/ = et µ̄j (B x̄ − c)j = 0, pour tout 1 ≤ j ≤ p.
B 0 µ̄ c
Preuve.
Similaire au cas avec contraintes d’égalités où on utilise ici (4.1.5).
Remarque 4.1.3 On peut vérifier facilement que ce multiplicateur λ̄ est solution de
maxp L(x̄, λ)
λ∈R
λ≥0
où
L(x, λ) = 12 (Ax, x) − (b, x) + (Bx − c, λ)
est appelé le Lagrangien du problème quadratique (4.16).
Exercice 4.1.2 Soit

K = Rp+ = {x ∈ Rp , x ≥ 0} .
Montrer que PK (x) = max(0, x) = (max(xi , 0))1≤i≤p

Dans cette section on donnera deux méthodes numériques pour calculer la solution
d’un problème d’optimisation avec contraintes.
On se limitera dans ce chapitre à étudier la méthode du gradient projeté pour un
problème quadratique pour K quelconque et la méthode d’Uzawa pour un problème qua-
dratique avec contraintes d’inégalité.
4.2.1 Méthode du gradient projeté

On rappelle que dans le cas du problème sans contraintes, la méthode du gradient à
pas fixe α s’écrit, pour x(0) donné,
x(k+1) = x(k) − α∇J(x(k) ).
63
Lorsque qu’on minimise sur une partie K supposée fermée convexe, on n’est pas sur qu’à
chaque étape k, l’itéré x(k) reste dans K même si on part d’un x(0) admissible. Pour cette
raison, on choisit
x(k+1) = PK (x(k) − α∇J(x(k) )), (x(0) ∈ K,
si Pk est calculable, où Pk désigne la projection sur le fermé convexe K. L’algorithme du
gradient projeté s’écrit :
1. On choisit x(0) ∈ K,
(k)
d = −∇J(x(k) ),
x(k+1) = PK (x(k) + αd(k) )
Proposition 4.2.1 Soit A une matrice symétrique définie positive et K est un convexe
2
fermé non vide de Rn . Si 0 < α < , alors la suite générée par la méthode du gradient
ρ(A)
projeté converge vers x̄ réalisant le minimum de J(x) = 12 (Ax, x) − (b, x) sur K.
Preuve.
Commençons par montrer que si x̄ est la solution unique du problème (P ), alors, pour
tout α > 0, on a
x̄ = PK (x̄ − α∇J(x̄)).
En effet, comme x̄ est solution du problème de minimisation convexe (P), alors d’après
l’inéquation d’Euler, on a pour tout x ∈ K,
(∇J(x̄), x − x̄) ≥ 0,
donc, si α > 0, pour tout x ∈ K,
(x̄ − (x̄ − α∇J(x̄), x − x̄) = α(∇J(x̄), x − x̄) ≥ 0.
De la propriété caractérisant la projection sur un convexe (voir (4.5)) vient
x̄ = PK (x̄ − α∇J(x̄)).
Montrons que (x(k) ) converge vers x̄. On a pour k ≥ 0,
kx(k+1) − x̄k2 = kPK (x(k) − α∇J(x(k) )) − PK (x̄ − α∇J(x̄))k2

≤ k(x(k) − α∇J(x(k) )) − (x̄ − α∇J(x̄)k2
.
= k(I − αA)(x(k) − x̄)k2
≤ kI − αAk2 kx(k) − x̄k2
64
Par récurrence, on obtient alors

kx(k) − x̄k2 ≤ kI − αAkk2 kx(0) − x̄k2 = ρ(I − αA)k kx(0) − x̄k2 ,
car I − αA est symétrique.
2
Si 0 < ρ < ρ(A) , alors ρ(I − αA) < 1 et donc la suite (x(k) ) converge vers x̄.
Remarque 4.2.1 La méthode du gradient projeté devient intéressante si la projection
sur le convexe K est facile à calculer. C’est le cas par exemple d’une projection sur Rp+
ou sur un pavé de Rp . En général, le calcul de la projection sur K est délicat. D’où la
difficulté de la mise en oeuvre de cette méthode malgré son apparente simplicité.
4.2.2 Méthode d’Uzawa

On rappelle que pour une matrice A symétrique définie positive, si x̄ est le minimum
de
min 1 (Ax, x) − (b, x),
Bx≤c 2
alors il existe λ̄ ∈ Rp+ , tel que

Ax̄ + B T λ̄ = b et λj (B x̄ − c)j = 0, j = 1, . . . , p
et que ce multiplicateur λ̄ est solution de
max L(x̄, λ),
λ∈Rp+
où
L(x, λ) = 12 (Ax, x) − (b, x) + (Bx − c, λ).
Comme x 7→ L(x, λ̄) est convexe et elle vérifie
∂L
(x̄, λ̄) = Ax̄ − b + B T λ̄ = 0,
∂x
donc x̄ est aussi solution de
min L(x, λ̄).
x∈Rn
Ainsi ∀x ∈ Rn , ∀λ ∈ Rp+ , on a
L(x̄, λ) ≤ L(x̄, λ̄) ≤ L(x, λ̄).
On dit que (x̄, λ̄) est un point selle de L.
Le principe de la méthode d’Uzawa consiste à calculer numériquement ce point selle.
On se donne (x(0) , λ(0) ), et si à l’itération k on connaı̂t (x(k) , λ(k) ), on commence par
calculer x(k+1) solution de minn L(x, λ(k) ) sur Rn solution donc de système Ax = b−B T λ(k) ,
x∈R
puis calculer
∂L (k+1) (k)
λ(k+1) = PRp+ (λ(k) + α (x , λ )) = max(0, λ(k) + α(Bx(k) − c)),
∂λ
itéré obtenu par la méthode du gradient projeté, de pas α, appliquée au problème dual
max L(x(k) , λ).
λ∈Λ
65
Algorithme d’Uzawa pour un problème quadratique avec contraintes d’inégalités

Si (P ) est le problème
min 1 (Ax, x) − (b, x),
Bx≤c 2
alors l’algorithme de la méthode d’Uzawa de pas α s’écrit :
1. On choisit (x(0) , λ(0) ) ∈ Rn × Rp+ et α > 0.

(k+1)
x solution de Ax = b − B T λ(k)
(k+1)
λ = max(0, λ(k) + α(Bx(k) − c))
On peut utiliser une des méthodes numériques directes ou itératives pour déterminer
(k+1)
x
Proposition 4.2.2 Soit A une matrice symétrique définie positive et soit λ1 sa plus petite
2λ1
valeur propre. Alors si 0 < α < , la suite (x(k) ) générée par l’algorithme d’Uzawa
kBk22
converge vers x̄ solution du problème min 21 (Ax, x) − (b, x).
Bx≤c
Preuve. On a
Ax̄ − b + B T λ̄ = 0
et
Ax(k+1) − b + B T λ(k) = 0.
Donc
A(x(k+1) − x̄) = −B T (λ(k) − λ̄).
Si on applique le produit scalaire avec x(k+1) − x̄, on obtient
(A(x(k+1) − x̄), x(k+1) − x̄) = −(B T (λ(k) − λ̄), x(k+1) − x̄) = −(λ(k) − λ̄, B(x(k+1) − x̄)).
Ainsi
λ1 kx(k+1) − x̄k22 + (λ(k) − λ̄, B(x(k+1) − x̄)) ≤ 0. (4.17)
Or
kλ(k+1) −λ̄k2 = kPRp+ (λ(k) +α(Bx(k+1) −c)−PRp+ (λ̄+α(B x̄−c)k2 ≤ kλ(k) −λ̄+αB(x(k+1) −x̄)k2 .
kλ(k+1) − λ̄k22 ≤ kλ(k) − λ̄k22 + 2α(λ(k) − λ̄, B(x(k+1) − x̄)) + α2 kB(x(k+1) − x̄)k22 .
D’après (4.17), on déduit que
kλ(k+1) − λ̄k22 ≤ kλ(k) − λ̄k22 − 2αλ1 kx(k+1) − x̄k22 + α2 kBk22 kx(k+1) − x̄k22 (4.18)
66
2λ1
Si on choisit α tel que 0 < α < , alors β = α(2λ1 − αkBk22 ) > 0 et
kBk22
kλ(k+1) − λ̄k22 ≤ kλ(k) − λ̄k22 − βkx(k+1) − x̄k22 . (4.19)
Ainsi la suite (kλ(k) − λ̄k22 ) est décroissante, donc convergente et par conséquent (kx(k+1) −
x̄k22 ) converge vers 0 puisque
βkx(k+1) − x̄k22 ≤ kλ(k) − λ̄k22 − kλ(k+1) − λ̄k22 .
Exercice 4.2.1 Ecrire la méthode d’Uzawa pour un problème quadratique avec contraintes
d’égalités.
67
68
Références
[1] R. Bessi Optimisation quadratique, ENIT 2010.

[2] P.G. Ciarlet : Introduction à l’analyse numérique matricielle et à l’optimisation.
Edition Masson, 1982.
[3] P.G. Ciarlet, B. Miara, J.M. Thomas : Exercices d’analyse numérique matricielle
et d’optimisation avec solutions. Edition Masson, 1991.
[4] P. Ciarlet, H. Zidani : Optimisation quadratique, cours AO101, ENSTA.
[5] H. El Fekih, K. Griba, T. Hadhri : Cours d’analyse numérique matricielle. ENIT
1996.
69

Analyse Numer I Que I

Transféré par

Droits d'auteur :

Formats disponibles

Analyse Numer I Que I

Transféré par

Informations du document

Description originale:

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Analyse Numer I Que I

Transféré par

Droits d'auteur :

Formats disponibles

Analyse Numérique

Radhia Bessi & Maher Moakher

1 Méthodes directes pour la résolution des systèmes linéaires 1

2 Méthodes itératives pour la résolution des systèmes linéaires 19

3 Optimisation sans contraintes 37

3.4.3 Méthodes du gradient conjugué . . . . . . . . . . . . . . . . . . . . 49

4 Optimisation avec contraintes linéaires 55

Méthodes directes pour la résolution

qui est équivalent au système linéaire Ah x = bh , où

1.2 Principe des méthodes directes

1.2.1 Rappels et notations sur les vecteurs et les matrices

et la norme associée est la norme euclidienne donnée par

On notera l’espace vectoriel des matrices à n lignes et p colonnes et à coefficients dans

AT = (aji )1≤i≤p,1≤j≤n ∈ Mp,n (R).

La matrice transposée vérifie

(Ax, y) = (Ax)T y = xT AT y = (x, AT y) pour tout x ∈ Rp et y ∈ Rn .

Définition 1.2.1 Soit A = (aij )1≤i,j≤n ∈ Mn (R). A est dite :

La matrice B est appelée inverse de A et sera notée A−1 .

• diagonalisable si elle est semblable à une matrice diagonale.

On peut vérifier facilement que

1.2.2 Résolution d’un système triangulaire

Si U est la matrice triangulaire supérieure suivante :

A chaque étape k on a une division, (n − k) multiplication pour 1 ≤ k ≤ n et (n − k)

Pour la résolution d’un système triangulaire inférieur Lx = b, on utilise la méthode de

Cette méthode nécessite aussi n2 opérations.

Exemple 1.3.1 Soit à résoudre le système linéaire d’ordre 4, suivant :

Sous forme matricielle ce système s’écrit sous la forme Ax = b, pour

Dans la première étape on élimine la première inconnue x1 des équations 2,3 et 4 en

Gauss à remplacer chaque ligne Li de A et de b par Li − aαi11 L1 , pour i = 2, 3 et 4. Dans

Algorithme de Gauss et nombre d’opérations

Puis, on résout le système triangulaire obtenu.

Après la première étape d’élimination de Gauss, on obtient

Preuve. Unicité : On suppose que A = L1 U1 = L2 U2 avec Li , i = 1, 2 triangulaire

On montre l’existence par récurrence sur n, l’ordre de la matrice A à factoriser.

puis les coefficients  an1

Remarque 1.3.3 De la preuve ou de l’algorithme de la factorisation LU de la matrice A

Ensuite pour calculer le terme lip , i = p + 1, . . . , n il faut

Le nombre total d’opérations pour effectuer la factorisation LU de A est

Proposition 1.3.2 La factorisation A = LU conserve la structure bande.

1.4 Factorisation de Cholesky des matrices symétriques

L’opérateur transposé vérifie, pour toutes marices A et B de Mn (R),

Définition 1.4.2 Soit A ∈ Mn (R). La matrice A est dite :

Définition 1.4.3 Vecteurs et valeurs propres

Un résultat connu concernant la réduction de matrices symétriques est le suivant :

Remarques 1.4.1 Si A ∈ Mn (R) symétrique et O une matrice orthogonale telle que

alors λ1 , λ2 , . . . , λn sont les valeurs propres de A. De plus si f1 , f2 , . . . , fn sont les vecteurs

Proposition 1.4.2 Soit A ∈ Mn (R) une matrice symétrique.

1. La matrice A est semi-définie positive si et seulement si toutes ses valeurs propres

(Ax, x) = λ(x, x) = λkxk22 .

L’inégalité (1.6) sera fort utile pour la suite.

1.4.2 Factorisation des matrices symétriques

donc la matrice diagonale D est inversible. Soit

(Dx, x) = (DLT y, LT y) = (LDLT y, y) = (Ay, y) > 0,

Remarque 1.4.3 (Factorisation de Crout)

1.4.3 Factorisation de Cholesky

Remarque 1.4.4 La matrice B de la factorisation de Cholesky de la matrice symétrique

Algorithme de la factorisation de Cholesky

On commence par calculer terme à terme la première colonne de la matrice triangulaire B,

k · kA,ε = kAk ≤ ρ(A) + .