1.

2 Méthodes directes : Élimination de Gauss et ses variantes

La formule de Cramer nous apporte une façon systématique de résoudre un
système linéaire :
det(Ai (b))
Trouver un vecteur x ∈ Rn , tel que Ax = b ⇐⇒ ∀i ∈ [1, . . . , n] xi = .
det(A)
Toutefois il n’est pas toujours possible appliquer cette approche en pratique. Effec-
tivement, les calculs des déterminants se font avec une relation récursive, et l’effort
de calcul est de l’ordre (n + 1)! flops. Or la dimension typique d’un problème dans
les applications est de l’ordre de n = 102 à n = 107 .

Figure 5 – cf. Livre A. Quarteroni, R. Sacco, F. Saleri “Numerical Mathematics”

Objectif : Construire des algorithmes de resolution avec un coût de calculs

raisonnable.

1.2.1 Algorithmes de descente et remontée

Définition 1.33 (Matrices triangulaires).
(i) On dit qu’une matrice U est triangulaire supérieure si ∀ i > j Ui,j = 0,
(ii) On dit qu’une matrice L est triangulaire inférieur si ∀ i < j Li,j = 0,
C’est à dire :
   
u11 u21 . . . u1n l11 0 ... 0
 0 u22 . . . u2n   l21 l22 ... 0 
U =  .. .. . . . ..  , L =  .. .. ... .. .
   
 . . .   . . . 
0 0 . . . unn ln1 ln2 . . . lnn

Proposition 1.34. Soit L1 , L2 (resp. U1 , U2 ) deux matrices triangulaires infé-

rieures (resp. supérieures), le produit L1 L2 (resp. U1 U2 ) est une matrice triangulaire
inférieure (resp. supérieure).

Démonstration.

17
Proposition 1.35. Soit L (resp. U ) une matrice triangulaire inférieure (resp. supé-
rieure) inversible, son inverse L−1 (resp. U −1 ) est une matrice triangulaire inférieure
(resp. supérieure).

Démonstration.

Dans le cas ou le système est de la forme :

Ly = b ou U x = y

on a des algorithme de résolution très efficaces.

Exemple. Démontrons l’algorithme sur un exemple simple : Soit L ∈ M3×3 (R) une
matrice inversible triangulaire inférieure.

A présent, soient L ∈ Mn×n (R) une matrice inversible triangulaire inférieure et

18
b ∈ Rn un vecteur colonne. La résolution de Ly = b sous la forme

    
l11 0 . . . 0 y1 b1  
 l21 l22 . . . 0   y2   b2  
 .. .. . . ..   ..  =  ..  
    
 . . . .  .   .  
ln1 ln2 . . . lnn yn bn  y

s’écrit sous la forme de l’algorithme de descente :

Remarque. Pour que cet algorithme puisse être appliqué, on doit avoir lii 6= 0 ce
qui revient à dire que la matrice L est inversible.

Algorithme 1 : [de descente] Résoudre Ly = b.

pour i = 1 à n faire
yi ← bi ;
pour j = 1 à i − 1 faire
yi ← yi − lij yj ;
fin
yi ← yi /lii
fin

Pour calculer la complexité de l’algorithme on doit estimer le nombre d’opérations

algébriques à chaque étape : à chaque itération i, on effectue i − 1 sommes, i − 1
produits et une division, soit au total
n n
!
X X
(2(i − 1) + 1) = 2 i − n = n2 opérations.
i=1 i=1

De même, on construit l’algorithme de remontée pour résolution d’un système

triangulaire supérieure donné par :
x
    
u11 u21 . . . u1n x1 y1  
 0 u22 . . . u2n   x2   y2  
 .. .. . . ..   ..  =  ..  
    
 . . . .  .   . 

0 0 . . . unn xn yn  

Soit U ∈ Mn×n (R) une matrice inversible triangulaire supérieure, alors x ∈ Rn , tel
que U x = y, vérifie :

19
 Algorithme 2 : [de remontée] Résoudre U x = y
pour i = n à 1 faire
xi ← y i ;
pour j = i + 1 à n faire
xi ← xi − uij xj ;
fin
xi ← xi /uii
fin

1.2.2 Méthode de Gauss : Réduction

L’idée de la méthode d’élimination de Gauss est de réduire le système linéaire
Ax = b à un système équivalent (c’est à dire à un système qui admet la même
solution que le système initial) de la forme suivante :

U x = b,
e

où U est une matrice triangulaire supérieure et b

e est le vecteur de second membre
modifié. Ensuite ce nouveau système se résout par l’algorithme de remontée.

Elimination Remontee

1
1

0
1

Figure 6 – Idée de l’algorithme de Gauss

Notons le système initial par A(1) x = b(1) . On suppose maintenant a11 6= 0 et on

introduit les multiplicateurs
ai1
mi1 = , i = 2, . . . , n.
a11
La première étape de l’algorithme est alors :
 (1) (1) (1)     (1)   (1) (1) (1)     (1) 
a11 a12 . . . a1n x1 b1 a11 a12 . . . a1n x1 b1
a(1) a(1) (1)  
. . . a2n   x2  b2 
  (1)  0 a(2) . . . a(2)   x2  b(2) 
 21 22 22 2n     2 
 . .. ..   ..  =  .  ⇐⇒  . . .   ..  =  . 
 
 .. . .  .   .  . .
 . .
. .
.   .   .. 
(1) (1) (1) xn (1) (2) (2) xn (2)
an1 an2 . . . ann bn 0 an2 . . . ann bn

où,

On note ce dernier système par A(2) x = b(2) . En procédant ainsi de suite, on

20
construit un système A(k) x = b(k) où
 (1) (1) (1) 
a11 a12 . . . . . . . . . a1n
 0 a(2) (2)
a2n 
22
 . .. .. 
 
 .. . . 
(k)
A = 
(k) (k)

 0 . . . 0 akk . . . akn 
 . .. .. .. 

 .. . . . 
(k) (k)
0 . . . 0 ank . . . ann
(k) (k)
Le coefficient akk est appelé le pivot. Si on suppose ∀ k = 1, . . . , n − 1 akk 6= 0, pour
k = n on obtient A(n) x = b(n) avec
 (1) (1) (1) 
a11 a12 . . . . . . . . . a1n
 0 a(2) (2)
a2n 
22
 . . .
 
 . . . . .
. 

(n)
A =  (k)
 est une matrice triangulaire supérieure.
(k) 
0 a kk . . . a kn 
 . . .. 

 . . . . . 
(n)
0 ann
(k)
Résumé de passage d’une étape k à une étape k + 1 : Soit akk 6= 0, alors en définis-
sant les multiplicateurs :
(k)
aik
mik = (k)
, i = k + 1, . . . , n (2)
akk

on définit
(k+1) (k) (k)
aij = aij − mik akj , i, j = k + 1, . . . , n
(k+1) (k) (k) (3)
bi = bi − mik bk , i = k + 1, . . . , n.

Figure 7 – Voir pour un exemple détaillé d’application de la méthode de Gauss.

Exercice. Calculer la solution H3 x = b où H3 est une matrice de Hilbert d’ordre

3 définie dans exemple 1.1.7 et b = [11/6, 13/12, 47/60]t , en utilisant (2) et (3).

En résumé :

21
 Algorithme 3 : Méthode de Gauss
Données : A une matrice inversible de taille n, b un vecteur colonne de
taille n
Résultat : x un vecteur de taille n tel que Ax = b
// Triangularisation de A et transformation de b
pour k = 1 à n − 1 faire
pour i = k + 1 à n faire
mik = aik /akk
pour j = k à n faire
aij = aij − mik akj
fin
bi = bi − mik bk
fin
fin
// Résolution du système transformé par l’algorithme de remontée
pour i = n à 1 faire
x i = bi
pour j = i + 1 à n faire
xi = xi − aij xj
fin
xi = xi /aii
fin
 
1 2 3
Remarque (Choix du pivot (voir TD)). Considérons A = 2 4 5 une matrice
7 8 9
inversible. L’application
 de procédure d’élimination de Gauss au première étape
1 2 3

donne : A = 0 0 −1 (et l’exécution s’arrête si on pense à rajouter la vérification
0 −6 −12
(k)
dans le code). Si akk = 0 on cherche un autre pivot :
(k) (k)
1. Soit en permutant les lignes k et p où |apk | = maxj=k+1,...,n |ajk |.
(k) (k)
2. Soit en permutant les colonnes k et p où |akp | = maxj=k+1,...,n |akj |.
(k) (k)
3. Soit en permutant les deux |apq | = maxi=k+1,...,n maxj=k+1,...,n |aij |.

Considérons maintenant le coût de calculs de la méthode de Gauss. Pour l’étape

d’élimination on a besoin de

2(n − 1)n(n + 1)/3 + n(n − 1)

opérations. Ensuite n2 pour l’algorithme de remontée. Ainsi au total 2n3 /3 + 2n2

opérations. Si on néglige les termes d’ordre inférieure on obtient le coût O(2n3 ).

22
 10 4

Temps d execution[s] 10 2

10 0

10 -2

10 -4

10 -6
10 1 10 2 10 3 10 4
Taille de matrice[n]
Figure 8 – Vérification numérique

Figure 9 – Coût Numérique

1.2.3 Méthode de Gauss : factorisation LU

On peut voir l’algorithme décrit dans la section précédente comme une factori-
sation. Effectivement, en posant
 
1 ... 0 0 ... 0
 .. . . . .. .. .. 
. . . .
0 1 0 0
 
Mk =   = In − mk etk

0 −mk+1,k 1 0
. . .. .. . . .. 
 .. .. . . . .
0 . . . −mn,k 0 ... 1

où mk = (0, ..., 0, mk+1,k , ..., mn,k )t ∈ Rn , on trouve

A(k+1) = Mk A(k) , et donc Mn−1 Mn−2 . . . M1 A = U.

La matrice Mk à diagonale unité admet une matrice inverse définie par

Mk−1 = 2In − Mk = In + mk eTk .

23
On remarque que (mj eTj )(mk eTk ) est une matrice nulle si j 6= k. Par conséquent,
 
1 0 ... ... 0
.. 
 m21 1 .

n−1
 . ... .. 
mk ek =  ..
X
−1 −1 −1
M1 M2 . . . Mn−1 = In +  m32 . .
 . .. ..
k=1
 . . . .

0
mn1 mn2 . . . mn,n−1 1
On définit alors
L = (Mn−1 Mn−2 . . . M1 )−1 = M1−1 M2−1 . . . Mn−1
−1

Et ainsi, A = LU .

Elimination

A U
xM ... xM
n-1 1
0
{
-1
L

Exercice. Calculer la décomposition LU pour la matrice A = H3 .

Algorithme 4 : Factorisation LU
Données : A une matrice inversible de taille n, b un vecteur colonne de
taille n
Résultat : x un vecteur de taille n tel que Ax = b
// Calcul de la décomposition LU . On initialise L = In . U est stockée dans
A.
pour k = 1 à n − 1 faire
pour i = k + 1 à n faire
`ik = aik /akk
pour j = k à n faire
aij = aij − `ik akj
fin
fin
fin
// Résolution de Ly = b par descente
pour i = 1 à n faire
yi = bi − i−1
P
k=1 `ik yk
fin
// Résolution de U x = y par remontée
pour i = n à 1 faire 
xi = 1/aii yi − nj=i+1 aij xj
P

fin

Un autre avantage de la factorisation LU (par rapport à la réduction de Gauss)

est qu’elle ne dépend pas du vecteur b. On peut donc la calculer une fois pour A et
ensuite seulement appliquer les algorithmes de descente et remontée pour chaque b.

24
1.2.4 Méthode de Gauss : conditions sur A
Pour s’assurer que la méthode de Gauss atteint l’étape (n−1), il ne suffit pas que
A soit inversible, car cela ne garanti pas la condition nécessaire ∀ k = 1, . . . , n − 1
(k)
akk 6= 0.
Theorem 1.36. Soit A une matrice (pas forcément inversible). Elle admet une
décomposition LU unique où L est à diagonale unité, i.e. ∀ i = 1, . . . , n Lii = 1 et U
est une matrice triangulaire supérieure, si et seulement si toutes les sous matrices
Ak , pour k = 1, . . . , n − 1, sont inversibles, où

Ak ∈ Mk×k (R) et (Ak )ij = Aij , ∀(i, j) ∈ [1, . . . , k]2 .

Démonstration.

25
Exercice. Dire si les matrices suivantes admettent une unique décomposition LU :
     
1 2 0 1 0 1
B= , C= , D= .
1 2 1 0 0 2

En fait, C n’admet aucune décomposition LU alors que D admet une infinité de de

décomposition de la forme :
   
1 0 0 1
D = Lα Uα , avec Lα = , Uα = .
α 1 0 2−α

Proposition 1.37. Si A est une matrice symétrique définie positive, alors ses sous
matrices sont inversibles, et donc elle admet une décomposition LU .
Démonstration.

26
Remarque (Stabilité de la méthode LU). Pour une analyse profond d’erreurs d’ar-
rondi (erreurs machine) vous pouvez consulter A. Quarteroni, R. Sacco, F. Saleri. Ici
on mentionne juste que la stabilité de l’algorithme d’élimination de Gauss dépend
(k)
max |Aij |
de magnitudes des éléments de A , plus précisément de ρ =
(k) i,j,k
max |Aij |
. Mauvaise
i,j
nouvelle : croissance exponentielle ρ ∼ 2n est possible.
Exemple (Wilkinson). Considérons une matrice définie comme suit
 
1 0 0 0 0 0 0 0 0 1
−1 1 0 0 0 0 0 0 0 1
 
−1 −1 1 0 0 0 0 0 0 1
 
−1 −1 −1 1 0 0 0 0 0 1
, si i = j ∨ j = n

 1 
−1 −1

−1 −1 1 0 0 0 0 1
aij = −1 , si i > j , A=
−1 −1
.
−1 −1 −1 1 0 0 0 1
0 , sinon.
  
−1 −1 −1 −1 −1 −1 1 0 0 1
 
−1 −1 −1 −1 −1 −1 −1 1 0 1
 
−1 −1 −1 −1 −1 −1 −1 −1 1 1
−1 −1 −1 −1 −1 −1 −1 −1 −1 1

Alors on a :
, si i = j

, si i = j,  1

 1 
A = LU , lij = −1 , si i > j, uij = 2 i-1 , si j = n
0 , sinon ; , sinon.
 
 0

t
√
Cependant, bonne nouvelle : en pratique, ρ ∼ n.

1.2.5 Méthode Choleski

On va maintenant étudier la méthode de Choleski, qui est une méthode directe
adaptée au cas où A est symétrique définie positive. On rappelle (Théorème 1.28)
que si A est symétrique définie positive elle est en particulier inversible.

27
Proposition 1.38. Si A est une matrice symétrique définie positive, alors elle admet
une décomposition unique de la forme A = HH t où H est une matrice triangulaire
inférieure à diagonale positive.
Démonstration.

28
 Les éléments de H sont définis à la façon suivante :

Exercice. Démontrez les formules ci-dessus.

Algorithme 5 : Méthode de Choleski
pour k = 1 à n faire
hkk = akk
pour i = 1 à k − 1 faire
hkk = hkk − h2ik
fin
√
hkk = hkk
pour i = k + 1 à n faire
hik = aik
pour j = 1 à k − 1 faire
hik = hik − hij hkj
fin
hik = hik /hkk
fin
fin

Remarque (Pivot partiel et Choleski.). Considérons une matrice A symétrique

définie positive. On n’a pas besoin de permutation pour obtenir sa décomposition
HH t de Choleski. Par contre, numériquement, on utilise malgré tout la technique
de pivot partiel pour minimiser les erreurs d’arrondi.
Comparaison Gauss/Choleski : Il est important de souligner que l’algorithme
de Choleski permet à garder la nature symétrique de la matrice A ce qui permet
d’économiser de la mémoire. Le coût de calcul de la résolution d’un système linéaire
par la méthode de Choleski est n3 /3 + O(n2 ), tandis que la méthode de Gauss
demande 2n3 /3 + O(n2 ) opérations (cf TD). Dans le cas d’une matrice symétrique
définie positive, la méthode de Choleski est donc environ deux fois moins chère. Mais
attention dans le contexte actuel des ordinateurs modernes ça ne veut pas dire qu’elle
est deux fois plus rapide (la gestion de la mémoire joue un rôle très important).

1.2.6 Factorisation QR : orthogonalisation de Gram-Schmidt

Idée : Factoriser la matrice A sous la forme A = QR où R est une matrice trian-
gulaire supérieure et Q est une matrice orthogonale (Q−1 = Qt ). La résolution du
problème linéaire se fait alors en appliquant l’algorithme de remontée au vecteur
Qt b.
Remarque. Dans un cas d’une matrice à coefficients complexes, Q est une matrice
unitaire (Q−1 = Q∗ ).
Cette factorisation se construit en utilisant soit l’algorithme d’orthogonalisation
de Gram-Schmidt soit les matrices de transformation : Householder ou Givens (cf
TD).

29
Theorem 1.39 (Gram-Schmidt). Soit un ensemble de k vecteurs (v1 , . . . , vk ) li-
néairement indépendant. On peut construire un ensemble de vecteurs (q1 , . . . , qk )
vérifiant :
(i) Vect(v1 , . . . , vk ) = Vect(q1 , . . . , qk )
(ii) (qi , qj ) = δij

Démonstration.

Ainsi les nouveaux vecteurs qi sont construit à la façon suivante :

i−1
X
e i = vi −
q (vi , qj ) qj , ei / ke
qi = q qi k
j=1

Algorithme 6 : Algorithme de Gram-Schmidt

pour k = 1 à n faire
qi = vi
pour k = 1 à i − 1 faire
qi = qi − (vi , qj )qj
fin
qi = qi /kqi k
fin

30
Remarque. Attention, le calcul de la complexité pour cet algorithme ne peut pas
être déterminée directement. En effet, les produits scalaires ne sont pas des opéra-
tions élémentaires

Theorem 1.40. Soit A une matrice inversible. Elle admet une décomposition QR
où Q est une matrice orthogonale et R est une matrice triangulaire supérieure.

Démonstration.

31
Remarque. Numériquement, on évite l’utilisation de cette approche car l’algo-
rithme d’orthogonalisation de G.S. est instable et conduit à une matrice Q non
parfaitement orthogonale.
Exemple. Les ordinateurs ne font pas les opérations sur R “correctement”. Calculons
décomposition QR de matrice de Hilbert H10 ∈ M10×10 en utilisant l’algorithme 6,
et faisons une simple vérification :
1.0000 0.0000 −0.0000 0.0000 −0.0000 0.0000 −0.0000 −0.0000 −0.0000 −0.0000
 
 0.0000 1.0000 −0.0000 0.0000 −0.0000 0.0000 −0.0000 −0.0000 −0.0000 −0.0000
−0.0000 −0.0000 1.0000 0.0000 −0.0000 0.0000 −0.0000 −0.0000 −0.0000 −0.0000
 
 0.0000 0.0000 0.0000 1.0000 −0.0000 0.0000 −0.0000 −0.0000 −0.0000 −0.0000
 
−0.0000 −0.0000 −0.0000 −0.0000 1.0000 0.0000 −0.0008 −0.0007 −0.0007 −0.0006
 
t
Q Q=
 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 −0.0540 −0.0430 −0.0360 −0.0289

−0.0000 −0.0000 −0.0000 −0.0000 −0.0008 −0.0540 1.0000 0.9999 0.9998 0.9996 
 
−0.0000 −0.0000 −0.0000 −0.0000 −0.0007 −0.0430 0.9999 1.0000 1.0000 0.9999 
 
−0.0000 −0.0000 −0.0000 −0.0000 −0.0007 −0.0360 0.9998 1.0000 1.0000 1.0000 
−0.0000 −0.0000 −0.0000 −0.0000 −0.0006 −0.0289 0.9996 0.9999 1.0000 1.0000

Il faut donc proposer un moyen stable de calculer la décomposition QR.

1.2.7 Méthode de Householder

Soit v ∈ Rn . On définit la matrice élémentaire de Householder H(v) comme suit

v vt
H(v) = I − 2
kvk2

Interprétation géométrique :

Proposition 1.41. La matrice H(v) est symétrique et orthogonale et H(v±kvkek )v =

∓kvkek

Démonstration.

32
Exercice. De la même manière que dans la preuve précédente, montrez que si
v = b − a 6= 0 avec kbk = kak, alors H(v)a = b.

Remarque. On déduit de la proposition précédente qu’on peut transformer un vec-

teur x quelconque en un vecteur n’ayant qu’une coordonnée non nulle en multipliant
par la matrice de Householder qui convient, via le choix du bon vecteur unitaire de
la base canonique

Exemple. Considérons x = [1, 1, 1, 1]t .

On remarque également que les matrices H(u) (pour des vecteurs u bien choisis)
nous permettent de réduire la matrice A à une matrice triangulaire inférieure. Cette
idée forme une base pour la méthode de Householder :

33
  
a11 a12 . . . a2n
 a21 a22 . . . a2n 
..  × H (a1 + ka1 k e1 ) où a1 = (a11 , . . . , an1 )
(1) t
A =  .. .. ..
 
 . . . . 
an1 an2 . . . ann
 (2) (2) (2) 
a11 a12 . . . a1n
 0 a(2) (2)
. . . a2n   t
22 (2) (2)
→ A(2) =  . .  × H (ã2 + kã2 k e2 ) où ã2 = 0, a22 , . . . , an2
(2)
.. ..
 
 .. . . .. 
(2) (2)
0 an2 . . . ann
 (3) (3) (3) (3)

a11 a12 a13 . . . a1n
(3) (3) (3) 
 0 a22 a23 . . . a2n 

 (3) (3) 
→ A(3)  0
= 0 a33 . . . a3n  ...
 .. .. .. .. . 
. .. 

 . . .
(3)
0 ... 0 . . . ann
A l’étape k on construit une matrice A(k) comme
A(k) = H (k−1) H (k−2) . . . H (1) A
 
(k) (k) (k) (k) (k)
avec∀ k H (k) = H ãk + kãk kek , ãk = (0 . . . , akk , . . . ank )t et
 (k) (k) (k) (k)

a11 a12 . . . a1k . . . a1n
(k) (k) (k) 
 0 a22 . . . a2k . . . a2n 

 .. .. .. .. 
(k)
 0 0 . . . . 
A = . . .
 
 . . .
.
(k)
0 akk . . . akn 
(k)

 .. .. .. .. . . .. 
 
 . . . . . . 
(k) (k)
0 . . . 0 ank . . . ann
On déduit alors :
A(n) = H (n−1) H (n−2) . . . H (1) A = R(triangulaire supérieure !)
et

t
t
t
Q = (H (n−1) H (n−2) . . . H (1) )−1 = H (1) . . . H (n−2) H (n−1) = H (1) . . . H (n−2) H (n−1)
Remarque. Attention ! ! La décomposition QR conduit a une matrice triangulaire
supérieure différente de la matrice U de la décomposition LU , R 6= U .
 3
La complexité de la méthode de Householder s’estime à O 4n3 (donc plus élevée
que pour l’algorithme de la décomposition LU ). Mais une des avantages importants
de la méthode de Householder est sa stabilité numérique.
Cette méthode s’applique également pour le cas des systèmes surdéterminés où
A ∈ Mn×p , n ≥ p est une matrice rectangulaire, ainsi que pour des problèmes de
minimisation écrits :
Trouver minp kAx − bk
x∈R

où A ∈ Mn×p , (n 6= p) (cf TD).

34