Numi 4

Chapitre IV
Systèmes d’Equations Linéaires

“It’s all linear algebra.” (fréquente proclamation de Gene Golub)
Considérons un système d’équations linéaires (aij , bj donnés)

a11 x1 + a12 x2 + . . . + a1n xn = b1
a21 x1 + a22 x2 + . . . + a2n xn = b2
.. .. .. .. (0.1)
. . . .
an1 x1 + an2 x2 + . . . + ann xn = bn
et cherchons sa solution x1 , . . . , xn . Très souvent, il est commode d’utiliser la notation matricielle
Ax = b. (0.2)
La méthode de Wanner. Prenons l’exemple

x1 + 2x2 − 3x3 = 5 (1)
3x1 + x2 − 2x3 = 1 (2)
2x1 − x2 + x3 = 2 (3)
On soustrait (2) de (1), puis (3) de (2), et (1) de (3). Cela donne
−2x1 + x2 − x3 = 4 (4)
x1 + 2x2 − 3x3 = −1 (5)
x1 − 3x2 + 4x3 = −3 (6)
Pour éliminer x1 , on soustrait (5) et (6) pour trouver
5x2 − 7x3 = 2 ⇒ x2 = 6, x3 = 4 et par (4) x1 = −1.
Tous les trois équations (4), (5) et (6) sont satisfaits. Mais il y a une mauvaise surprise si on essaie
ces “solutions” dans (1), (2) ou (3). Nous comprenons que des éliminations à l’aveugle ne sont pas
sans dangers !...
La règle de Cramer. (Une autre invention genevoise de 1750, cette fois sérieuse...) On apprend
en algèbre linéaire que, si det(aij ) 6= 0, les solutions sont données par
det(aij ...bi ...aij )
xk =
det(aij )
(où la k-ème colonne des aij est remplacée par les bi ) et pour det(aij ) on connait une formule
avec une somme sur le groupe de permutation des indices 1, 2, ..., n, i.e., avec n! de termes (ce qui
76 Systèmes d’Equations Linéaires
donne la même chose que de développer récursivement les déterminants par une ligne). Or si, par
exemple, n = 20, nous avons n! = 2.4 · 1018 , et même sur l’ordinateur le plus rapide du monde
(disons 109 opérations par seconde), ça prendrait 2 · 109 secondes = 3 · 107 minutes = 5 · 105 heures
= 3 · 103 jours = 10 années de calcul.
Bibliographie sur ce chapitre
Å. Björck (1996): Numerical Methods for Least Squares Problems. SIAM. [MA 65/387]
P.G. Ciarlet (1982): Introduction à l’analyse numérique matricielle et à l’optimisation, Masson.
J.J. Dongarra, C.B. Moler, J.R. Bunch & G.W. Stewart (1979): LINPACK Users’ Guide. SIAM.
D.K. Faddeev & V.N. Faddeeva (1963): Computational Methods of Linear Algebra. Freeman &
Co. [MA 65/271]
G.H. Golub & C.F. Van Loan (1989): Matrix Computations. Second edition. John Hopkins Univ.
Press. [MA 65/214]
N.J. Higham (1996): Accuracy and Stability of Numerical Algorithms. SIAM. [MA 65/379]
A.S. Householder (1964): The Theory of Matrices in Numerical Analysis. Blaisdell Publ. Comp.
[MA 65/262]
G.W. Stewart (1973): Introduction to Matrix Computations. Academic Press.
L.N. Trefethen & D. Bau (1997): Numerical Linear Algebra. SIAM. [MA 65/388]
J.H. Wilkinson (1969): Rundungsfehler. Springer-Verlag.
J.H. Wilkinson & C. Reinsch (1971): Handbook for Automatic Computation, Volume II, Linear
Algebra. Springer-Verlag.
IV.1 Elimination de Gauss

L’élimination dite “de Gauss” a été pratiquée pendant des siècles sans grand tam-tam, notamment
par Newton (voir chapitre II.1) et par Lagrange (en 1781 dans ses calculs astronomiques ; Oeuvres
V, p. 125-490). Toutefois, Gauss ayant le souci de prouver l’existence des solutions pour son
principium nostrum des moindres carrés (voir notices historiques du cours d’Algèbre Linéaire,
p. 17) décrit l’algorithme explicitement :
Si a11 6= 0, on peut éliminer la variable x1 dans les équations 2 à n à l’aide de l’équation 1, c.-à-d.,
on calcule
ai1
ℓi1 = pour i = 2, . . . , n (1.1)
a11
et on remplace la ligne i par
ligne i − ℓi1 ∗ ligne 1.
De cette manière, on obtient le système équivalent
(1) (1) (1) (1)
a11 x1 + a12 x2 + . . . + a1n xn = b1
(1) (1) (1)
a22 x2 + . . . + a2n xn = b2
.. .. .. (1.2)
. . .
(1)
an2 x2 + . . . + ann xn = b(1)
(1)
n
où
(1) (1)
a1j = a1j , b1 = b1 ,
(1) (1)
pour i = 2, . . . , n (1.3)
aij = aij − ℓi1 a1j bi = bi − ℓi1 b1
Systèmes d’Equations Linéaires 77
Le système (1.2) contient un sous-système de dimension n − 1 sur lequel on peut répéter la

procédure pour éliminer x2 dans les équations 3 à n. On multiplie la ligne 2 de (1.2) par ℓi2 =
(1) (1)
ai2 /a22 et on la soustrait de la ligne i. Après n − 1 étapes
(A, b) → (A(1) , b(1) ) → (A(2) , b(2) ) → . . . → (A(n−1) , b(n−1) ) =: (R, c)
on obtient un système triangulaire

r11 x1 + r12 x2 + . . . + r1n xn = c1
r22 x2 + . . . + r2n xn = c2
.. .. .. (1.4)
. . .
rnn xn = cn
qui se résoud facilement par “back substitution”
n
X
xn = cn /rnn , xi = (ci − rij xj )/rii pour i = n − 1, . . . , 1. (1.5)
j=i+1
Astuce pour la programmation : après l’élimination, les places de mémoire pour les a21 , a31 , . . .
ne seront plus nécessaires (on sait que ces grandeurs sont nulles) ; on peut donc y stocker les
ℓ21 , ℓ31 , . . . et tout l’algorithme peut se programmer en quelques lignes :
do ir=1,n-1 C ---- BACK SUBSTITUTION ---

do i=ir+1,n x(n)=b(n)/a(n,n)
a(i,ir)=a(i,ir)/a(ir,ir) do i=n-1,1,-1
do j=ir+1,n sum=0.
a(i,j)=a(i,j)-a(i,ir)*a(ir,j) do j=i+1,n
end do sum=sum+a(i,j)*x(j)
b(i)=b(i)-a(i,ir)*b(ir) end do
end do x(i)=(b(i)-sum)/a(i,i)
end do end do
Théorème 1.1 L’élimination de Gauss équivaut à une factorisation
A = LR (1.6)
où    
1 r11 r12 ... r1n
   
 ℓ21 1   r22 ... r2n 
L=
 .. .. .. ,
 R= 
 .. .. 
. (1.7)
 . . .   . . 
ℓn1 ... ℓn,n−1 1 rnn
La formule (1.6) s’appelle “décomposition LR” (left - right) de la matrice A.
Démonstration. En utilisant les matrices

   
1 1
 −ℓ 1  0 1 
 21   
   
 −ℓ
L1 =  31 0 1 ,
 0
L2 =  −ℓ32 1 ,
 ... (1.8)
 .. .. .. ..   .. .. .. .. 
 . . . .  . . . . 
−ℓn1 0 ... 0 1 0 −ℓn2 ... 0 1
le premier pas de l’élimination de Gauss correspond à une multiplication de A avec L1 , le deuxième

avec L2 , etc.,
L1 A = A(1) , L2 A(1) = A(2) , ... , Ln−1 A(n−2) = A(n−1) = R.
Par conséquent,
R = (Ln−1 Ln−2 · . . . · L1 ) · A et A = (Ln−1 Ln−2 · . . . · L1 )−1 · R.
Il reste à montrer que la matrice L de (1.7) est égale à (Ln−1 Ln−2 · . . . · L1 )−1 . Pour ceci, nous
appliquons la même procédure à la matrice L. La multiplication de L avec L1 élimine les éléments
de la première colonne en-dessous de la diagonale, puis la multiplication avec L2 élimine ceux de
la deuxième colonne, etc. Finalement, on obtient (Ln−1 Ln−2 · . . . · L1 ) · L = I =identité, ce qu’il
fallait démontrer.
Calcul du déterminant d’une matrice. La formule (1.6) implique que det A = det L · det R. On
obtient
det A = r11 · . . . · rnn (1.9)
i.e., le déterminant est le produit des pivots.
Résolution de systèmes linéaires. En pratique, on rencontre souvent la situation où il faut résoudre
une suite de systèmes linéaires Ax = b, Ax′ = b′ , Ax′′ = b′′ , etc., possédant tous la même matrice.
Très souvent, on connaı̂t b′ seulement après la résolution du premier système.
C’est la raison pour laquelle on écrit, en général, le programme pour l’élimination de Gauss en
deux sous-programmes :
DEC – calculer la décomposition LR (voir (1.6)) de la matrice;
SOL – résoudre le système Ax = b. D’abord on calcule le vecteur c (voir (1.4)), défini par Lc = b,
puis on résoud le système triangulaire Rx = c.
Pour le problème ci-dessus, on appelle une fois le sous-programme DEC et puis, pour chaque
système linéaire, le sous-programme SOL.
Calcul de l’inverse d’une matrice. Si on choisit pour les b, b′ , b′′ ci-dessus les vecteurs de base
(1, 0, ..., 0)T , (0, 1, ..., 0)T , etc., on obtient pour les x, x′ , x′′ , etc., les colonnes de la matrice inverse
A−1 .
Coût de l’élimination de Gauss. Pour le passage de A à A(1) , on a besoin de
n − 1 divisions (voir (1.1)) et de
(n − 1)2 multiplications et additions (voir (1.3)).
Le calcul de A(2) nécessite n − 2 divisions et (n − 2)2 multiplications et additions, etc. Comme le
travail dû aux divisions est ici négligeable, le coût total de la décomposition LR s’élève à environ
Z
2 2 2 2
n n3
(n − 1) + (n − 2) + . . . + 2 + 1 ≈ x2 dx = opérations
0 3
(opération = multiplication + addition).
En revenant à l’exemple ci-dessus d’une matrice de dimension 20 × 20, l’algorithme de Gauss
nécessite ≈ 2600 opérations, d’un facteur 10−15 fois plus petit que le travail des déterminants.
Exemple numérique. Prenons une matrice 60 × 60 avec coefficients choisies aléatoirement entre
−1 ≤ aij ≤ 1 et calculons A−1 par la méthode de Gauss (en simple précision). Puis, on contrôle
en double précision l’erreur des éléments. Le résultat est présenté en Fig. IV.1 à gauche (noir = 2
décimales justes, blanc = 7 décimales justes). Le résultat semble à désirer!! Mais il donne lieu à
une nouvelle découverte: The Scottish Kilt Phenomenon!!
A−1 = A−1 =
F IG . IV.1: Erreurs des éléments de A−1 d’une matrice aléatoire 60 × 60; sans recherche de pivot
(gauche), avec recherche de pivot (droite) ; noir = 2 décimales justes, blanc = 7 décimales justes
Algorithme avec recherche de pivot.
Exemple 1.2 (Forsythe) Considérons le système
1 · 10−4 · x1 + 1 · x2 = 1.0001
(1.10)
1 · x1 + 1 · x2 = 2 .
La solution exacte est, comme on voit, x1 = 1 et x2 = 1. Appliquons l’élimination de Gauss et

simulons un calcul en virgule flottante avec 3 chiffres significatifs (en base 10).
(1)
Si l’on prend a11 = 1 · 10−4 comme pivot, on obtient ℓ21 = a21 /a11 = 1.00 · 104 , a22 =
1.00−1.00·104 = −1.00·104. L’information contenue dans la valeur de a22 = 1 a tout simplément
(1)
disparue. Il est clair que le reste du calcul est faux. Regardons: b2 = 2.00 − 1.00 · 104 =
(1) (1)
−1.00 · 104 . Par conséquent, x2 = b2 /a22 = 1.00 (exacte!, la première équation n’a pas été
endommagé), mais pour x1 nous obtenons
x1 = (b1 − a12 x2 )/a11 = (1.00 − 1.00 ∗ 1.00)/(1.00 · 10−4 ) = 0.
Le résultat numérique, obtenu pour x1 , est faux.
Nous voyons à cet exemple, qu’il faut éviter qu’un des arr deviendrait trop petit. L’idée est
alors de ramener un des aij 6= 0 à la place du arr par des échanges de lignes (i.e., échange des
équations) et/ou des échanges de colonnes (i.e., échange des xi ), pour le rendre le plus grand
possible. L’algorithme le plus souvent utilisé dans les codes est le suivant :
Recherche partielle de pivot. On ne se contente pas d’un pivot différent de zéro (arr 6= 0), mais
on échange les équations de (0.1) afin que arr soit le plus grand élément (en valeur absolue) des
air , (i = r, r + 1, . . . , n). De cette manière on a toujours |ℓir | ≤ 1. Pour la programmation, il suffit
d’insérer dans le code ci-dessus, après la première ligne, les commandes :

c --- recherche du pivot ----
c --- echange -----
pgval=0.
do j=ir,n
izero=ir
store=a(ir,j)
do i=ir,n
a(ir,j)=a(izero,j)
valabs=abs(a(i,ir))
a(izero,j)=store
if(valabs.gt.pgval)then
end do
pgval=valabs
store=b(ir)
izero=i
b(ir)=a(izero)
end if
b(izero)=store
end do
IV.2 Etude des erreurs ; “Backward Error Analysis”

10−1 errmax 10−1 errmax
sans recherche de pivot avec recherche de pivot
10−2 10−2
10−3 10−3
10−4 10−4
10−5 10−5
10−6 10−6
10−7 10−7
matrice quelconque matrice quelconque
10 −8 n 10−8 n
10 20 30 40 50 10 20 30 40 50
10−1 errmax 10−1 errmax
sans recherche de pivot avec recherche de pivot
10−2 10−2
10−3 10−3
10−4 10−4
10−5 10−5
10−6 10−6
10−7 10−7
matrice orthogonale matrice orthogonale
10 −8 n 10−8 n
10 20 30 40 50 10 20 30 40 50
F IG . IV.2: Erreurs pour 1 million de systèmes linéaires de dimensions 5×5 à 55×55
Les pionniers de l’analyse numérique (Hotelling, von Neumann, Goldstine) dans les années ’40 ont
rencontré d’insurmontables difficultés à analyser les erreurs d’arrondi de la solution des systèmes
linéaires. Ils sont arrivés à la conclusion que les dimensions plus grandes que 10 ou 12 seraient
impossibles. Malgré ces prédictions pessimistes, les résultats numériques n’ont pas été si mauvais.
Faisons une expérience numérique (voir figure IV.2) : pour chaque n = 5, 6, . . . , 55 nous choi-
sissons 2000 matrices aléatoires avec coefficients aij uniformément distribués dans [−1, 1] et des
solutions xi uniformément distribuées dans [−1, 1]. Alors on calcule en double précision les bj
pour cette solution exacte. Ensuite on applique l’algorithme de Gauss, une fois sans recherche
de pivot, et une fois avec recherche de pivot, en simple précision. L’erreur maxi |xnum i − xex
i | de
chaque résultat est représentée par un petit point dans les dessins supérieurs de la figure IV.2. Bien
que nous ne soyons pas surpris par les nombreuses erreurs sans recherche de pivot, quelques cas
demeurent inacceptables à droite ; bon nombre de résultats restent cependant bons !
Faisons une deuxième expérience : une matrice avec aij uniformément distribués dans [−1, 1] pour
j > i est complétée par aji = −aij , pour assurer que Q = (I −A)−1 (I +A) soit orthogonale (Cay-
ley ; voir Γǫoµǫτ ρία II.4). Cette matrice est calculée en double précision, le reste de l’expérience
continue comme auparavant (voir les résultats au bas de la figure IV.2). Cette fois-ci il n’y a pas
d’exception dans la bonne performance de l’algorithme de Gauss avec recherche de pivot.
La “Backward Error Analysis” de Wilkinson.
L’explication théorique de ces phénomènes a été un des grands challenges des années ’50. Il parais-
sait alors difficile d’arriver à un résultat, où même un John von Neumann avait jeté l’éponge !...
L’idée miraculeuse a finalement été publiée par Wilkinson (1961, J. Ass. Comp. Mach. 8) :
Supposons qu’un système de dimension 2 × 2 soit à transformer sur forme triangulaire par un pas
d’élimination

a11 a12 b1 El. par ℓ21 a11 a12 b1
−→
a21 a22 b2 0 a22 − ℓ21 a12 b2 − ℓ21 b1
avec ℓ21 = a21 /a11 .
Première source d’erreur : elle résulte du fait que l’ordinateur calcule avec un faux
ℓb21 = ℓ21 + ǫ où |ǫ| ≤ eps
(car |ℓ21 | ≤ 1 à cause du choix du pivot).

Idée : au lieu de poursuivre les dégâts occasionnés par cette erreur aux calculs ultérieurs et aux
solutions, nous cherchons à modifier les données pour rendre le calcul (théoriquement) correct : si
le a21 du début aurait été égal à a21 + ǫa11 , ce calcul avait été sans erreur !
!
a11 a12 b1 El. par ℓb21 a11 a12 b1
−→ (1) (1) .
a21 + ǫa11 a22 b2 0 a22 b2
(1)
Deuxième source d’erreur. Pour le calcul de a22 = a22 − ℓb21 a12 , il y a une multiplication et une
soustraction à faire ; ensuite, le résultat est placé dans une case de mémoire pour a22 . Les détails du
résultat dépendent de la manière dont le compilateur travaille. Souvent les opérations algébriques
en chaı̂ne se font sur un régistre plus long ; seulement la mise en mémoire provoque une erreur
d’arrondi notable. Sous cette hypothèse1, les deux quantités deviennent
(1) (1) (1) (1) (1) (1)
ab22 = a22 + e1 et bb2 = b2 + e2 , où |e1 | ≤ eps · |a22 | et |e2 | ≤ eps · |b2 |.
Nous transportons à nouveau ces erreurs du côté des données, et le calcul

!
a11 a12 b1 El. par ℓb21 a11 a12 b1
−→ (1) b(1)
a21 + ǫa11 a22 + e1 b2 + e2 0 ab22 b2
1
qui nous est agréable...
est sans erreur. Nous voyons donc que le résultat numérique du système linéaire est le résultat
exact d’un système dont la deuxième ligne a été modifiée par des quantités ≤ a · eps où a =
(k) (k)
max |aij , bi |.
Pour des systèmes de dimensions supérieures, on corrige plusieurs fois les données aij ; d’abord
pour i = 2, . . . , n, j = 1, . . . , n, ensuite pour i = 3, . . . , n, j = 2, . . . , n, etc. Nous arrivons au
célèbre théorème :
Théorème 2.1 (Wilkinson) Soit A une matrice inversible et L, b Rb le résultat numérique de l’élimi-
nation de Gauss (avec recherche de pivot, c.-à-d. |ℓbij | ≤ 1 pour tout i, j). Alors L
bR b =A b avec
 
0 0 0 ... 0 0

1 1 1 ... 1 1 
 
1 2 2 ... 2 2  (k)
|abij − aij | ≤ a · eps · 
1
 . où a = maxi,j,k |aij |. (2.1)
 2 3 ... 3 3 
 
 .. .. .. ... .. .. 
1 2 3 . . . n−1 n−1
Définition 2.2 Un algorithme pour résoudre un problème est numériquement stable (au sens de
“backward analysis”), si le résultat numérique peut être interprété comme un résultat exact pour
des données légèrement perturbées.
Par conséquent, si le résultat est faux, ce n’est pas la faute de la méthode, mais bien celle du
problème. Dans ce cas, on appelle le problème un problème mal conditionné. Nous allons étudier
ces problèmes plus en détail au paragraphe suivant.
Exemple. Calculons la solution (en simple précision) du système Ax = b avec
     
1/2 1/3 1/4 1/5 3511/13860 1/3
 1/3 1/4 1/5 1/6   277/1540   1/11 
     
A=  , b=  ; sol. exacte x =   .
 1/4 1/5 1/6 1/7   40877/291060   1/9 
1/5 1/6 1/7 1/8 3203/27720 1/7
.5000000 .3333333 .2500000 .2000000 .2533189
.3333333 .2500000 .2000000 .1666667 .1798701
.2500000 .2000000 .1666667 .1428571 .1404418
.2000000 .1666667 .1428571 .1250000 .1155483
--------------------------------------
.5000000 .3333333 .2500000 .2000000 .2533189
.6666667 .0277778 .0333333 .0333333 .0109909
.5000000 .0333333 .0416667 .0428571 .0137824
.4000000 .0333333 .0428571 .0450000 .0142208
--------------------------------------
.5000000 .3333333 .2500000 .2000000 .2533189
.6666667 .0333333 .0416667 .0428571 .0137824
.5000000 .8333330 -.0013889 -.0023809 -.0004945
.4000000 1.0000000 .0011905 .0021429 .0004384
--------------------------------------
.5000000 .3333333 .2500000 .2000000 .2533189
.6666667 .0333333 .0416667 .0428571 .0137824
.5000000 .8333330 -.0013889 -.0023809 -.0004945
.4000000 1.0000000 -.8571472 .0001020 .0000146
--------------------------------------
x(1)= .33333951 x(2)= .09086763
x(3)= .11118819 x(4)= .14281443
Les résultats montrent que seulement 3 à 4 décimales sont justes. Mais, pour l’honneur de notre
méthode, nous constatons que les résidus de ces solutions res = Ax − b sont correctes :
res(1)= -.00000003 res(2)= -.00000001
res(3)= .00000000 res(4)= -.00000001
Probl. bien cond. Probl. bien cond. Probl. mal cond.
(A, b) x (A, b) x (A, b) x
b b
xb b b
(A, b) (A, b)
b b
(A, b) xb xb
Alg. stable Alg. instable Alg. stable
F IG . IV.3: Schéma de la “Backward Error Analysis”
Ce problème est mal conditionné, malgré son apparence débonnaire. Nous sommes donc dans la
troisième case du schéma de la Fig. IV.3.
IV.3 La condition d’une matrice

En principe, un problème avec m données et n solutions possède m × n coefficients décrivant la
sensibilité de la n-ème solution par rapport à la m-ème donnée. Devant cette myriade de valeurs,
il est parfois préférable d’exprimer la condition par un seul nombre. On réussira cela à l’aide de
normes de vecteurs et de matrices (recherche initiée par A. Turing 1948).
Rappel sur la norme d’une matrice. Pour une matrice à m lignes et n colonnes, on définit
kAxk
kAk = max kAxk = max , (3.1)
kxk=1 x6=0 kxk
c.-à-d., la norme de A est le plus petit nombre kAk qui possède la propriété
kAxk ≤ kAk · kxk pour tout x ∈ IRn . (3.2)
Evidemment, kAk dépend des normes choisies dans IRn et IRm . Il y a des situations où l’on connaı̂t
des formules explicites pour kAk. Par exemple, si l’on prend la même norme dans les deux espaces
alors,
P
pour kxk1 = ni=1 |xi | , on a
m
X
kAk1 = max |aij | ; (3.3)
j=1,...,n
i=1
Pn
pour kxk2 = ( i=1 |xi |2 )1/2 , on a
q
kAk2 = plus grande valeur propre de AT A; (3.4)
pour kxk∞ = maxi=1,...,n |xi | , on a

n
X
kAk∞ = max |aij | . (3.5)
i=1,...,m
j=1
La norme kAk d’une matrice satisfait toutes les propriétés d’une norme. En plus, elle vérifie
kIk = 1 pour la matrice d’identité et kA · Bk ≤ kAk · kBk.
Après ce rappel sur la norme d’une matrice, essayons d’estimer la condition du problème Ax =
b. Pour ceci, considérons un deuxième système linéaire Abxb = bb avec des données perturbées
abij = aij (1 + ǫij ), |ǫij | ≤ ǫA ,
b
(3.6)
b
i = bi (1 + ǫi ), |ǫi | ≤ ǫb ,
où ǫA et ǫb spécifient la précision des données (par exemple ǫA ≤ eps, ǫb ≤ eps où eps est la
précision de l’ordinateur). Les hypothèses (3.6) impliquent (au moins pour les normes k · k1 et
k · k∞ ) que
kAb − Ak ≤ ǫA · kAk, kbb − bk ≤ ǫb · kbk. (3.7)
Notre premier résultat donne une estimation de kxb − xk, en supposant que (3.7) soit vrai. Un peu
plus loin, on donnera une estimation améliorée valable si (3.6) est satisfait.
Théorème 3.1 Considérons les deux systèmes linéaires Ax = b et Abxb = bb où A est une matrice
inversible. Si (3.7) est vérifié et si ǫA · κ(A) < 1, alors on a
kxb − xk κ(A)
≤ · (ǫA + ǫb ) (3.8)
kxk 1 − ǫA · κ(A)
où κ(A) := kAk · kA−1 k. Le nombre κ(A) s’appelle condition de la matrice A.
Démonstration. De bb − b = Abxb − Ax = (Ab − A)xb + A(xb − x), nous déduisons que

xb − x = A−1 −(Ab − A)xb + (bb − b) . (3.9)
Maintenant, prenons la norme de (3.9), utilisons l’inégalité du triangle, les estimations (3.7), kxbk ≤
kxk + kxb − xk et kbk = kAxk ≤ kAk · kxk. Nous obtenons ainsi

kxb − xk ≤ kA−1 k ǫA · kAk · (kxk + kxb − xk) + ǫb · kAk · kxk .
Ceci donne l’estimation (3.8).
La formule (3.8) montre que pour ǫA · κ(A) ≪ 1, l’amplification maximale de l’erreur des
données sur le résultat est de κ(A).
Propriétés de κ(A). Soit A une matrice inversible. Alors,
a) κ(A) ≥ 1 pour toute A,
b) κ(αA) = κ(A) pour
.
α 6= 0,
c) κ(A) = max kAyk min kAzk.
kyk=1 kzk=1
La propriété (c) permet d’étendre la définition de κ(A) aux matrices de dimension m × n avec
m 6= n.
Démonstration. La propriété (a) est une conséquence de 1 = kIk = kAA−1 k ≤ kAk · kA−1 k. La
propriété (b) est évidente. Pour montrer (c), nous utilisons
−1
−1 kA−1 xk kzk kAzk
kA k = max = max = min .
x6=0 kxk z6=0 kAzk z6=0 kzk
TAB . IV.1: Conditions de matrices de Hilbert et Vandermonde
n 2 4 6 8 10 12
κ(Hn ) 27 2.8 · 104 2.9 · 107 3.4 · 1010 3.5 · 1013 3.8 · 1016
κ(Vn ) 8 5.6 · 102 3.7 · 104 2.4 · 106 1.6 · 108 1.0 · 1010
Exemples de matrices ayant une grande condition. Considérons les matrices Hn (matrice de
Hilbert) et Vn (matrice de Vandermonde) définies par (cj = j/n)
1 n n
Hn = , Vn = ci−1
j .
i + j − 1 i,j=1 i,j=1
Leur condition pour la norme k · k∞ est donnée dans le tableau IV.1. La matrice Vn est précisément
la matrice du problème d’interpolation polynomiale pour noeuds équidistants. La mauvaise con-
dition de cette matrice est liée au mauvais comportement de cette interpolation que nous avons
remcontré au chapitre II.4.
Exemples de matrices ayant une petite condition. Une matrice U est orthogonale si U T U = I.
Pour la norme euclidienne, sa condition vaut 1 car kUk2 = 1 et kU −1 k2 = 1 (l’inverse U −1 = U T
est aussi orthogonale).
Concernant l’interpolation avec des fonctions splines, nous avons rencontré la matrice (voir le
paragraphe II.8, cas équidistant)
  
4 1 

1   
A=  1 4 .1. ..  n (3.10)
h 1 ..
.
..
. 


. .
Le facteur 1/h n’influence pas κ(A). Posons alors h = 1. Avec la formule (3.5), on vérifie
facilement que kAk∞ = 6. Pour estimer kA−1 k∞ , écrivons A sous la forme A = 4(I + N) où I
est l’identité et N contient le reste. On voit que kNk∞ = 1/2. En exprimant A−1 par une série
géométrique, on obtient
1 1
kA−1 k∞ ≤ 1 + kNk∞ + kNk2∞ + kNk3∞ + . . . ≤ .
4 2
Par conséquent, κ∞ (A) ≤ 3 indépendamment de la dimension du système.
IV.4 L’algorithme de Cholesky

Soit B une matrice quelconque et posons
A = BT B , (4.1)
i.e., l’élément aij de A est le produit scalaire des colonnes i et j de B 2 . Alors
A est symétrique (car le prod. scal. est symétrique; ou car AT = B T (B T )T = B T B = A)

(4.2)
2
On appelle cela aussi une Matrice de Gram.
et
A est définie positive (i.e., xT Ax > 0 pour x 6= 0), (4.3)
car xT Ax = xT B T Bx = (Bx)T (Bx) = y T y > 0). Il est nécessaire pour l’inégalité stricte que les
colonnes de B sont linéairement indépendantes 3 . Si Ax = λx, alors xT Ax = λ · xT x > 0, on voit
que chaque valeur propre d’une matrice symétrique et positive définie doit être réelle (voir cours
d’Algèbre) et > 0.
Question: Existe-t-il une “décomposition LR” symétrique
    
a11 a12 a13 a14 ℓ11 ℓ11 ℓ21 ℓ31 ℓ41
a a22 a23  
a24   ℓ21 ℓ22  ℓ22 ℓ32 ℓ42 
 21  
A = L LT ou  =  
 a31 a32 a33 a34   ℓ31 ℓ32 ℓ33  ℓ33 ℓ43 
a41 a42 a43 a44 ℓ41 ℓ42 ℓ43 ℓ44 ℓ44
(4.4)
Il est clair, par (4.1), que A doit être symétrique et définie positive.
Théorème. Pour chaque matrice symétrique et définie-positive existe une décomposition dite “de
Cholesky” 4 (4.4). L’algorithme de Cholesky ci-dessous est toujours numériquement stable. Il
n’est pas nécessaire de faire une recherche de pivot.
Calcul des ℓij .
Pas 1a. Calculons dans (4.4) la valeur de a11 . Elle est
√
a11 = (ℓ11 )2 donc ℓ11 = a11 . (4.5)
Question. Est-on sûr que a11 > 0 ? Oui, il suffit de poser dans la condition (4.3) le vecteur
x = (1, 0, 0, ...)T .
Pas 1b. Calculons dans (4.4) les valeurs de ai1 pour i = 2, 3, 4, .... On obtient
ai1 = ℓi1 · ℓ11 donc (ℓ21 est connu) ℓi1 = ai1 /ℓ11 . (4.6)
La division par ℓ11 ne pose pas de problème, car ℓ11 > 0.
Pas 2a. Calculons dans (4.4) la valeur de a22 . Elle est
q
a22 = (ℓ21 )2 + (ℓ22 )2 donc ℓ22 = a22 − (ℓ21 )2 . (4.7)
Question. Est-on sûr que a22 − (ℓ21 )2 > 0 ? C’est déjà plus difficile. On pose dans la condition
(4.3) le vecteur x = (u, 1, 0, ...)T . Ainsi, xT Ax, que nous savons positif, devient
a11 u2 + 2a21 u + a22 > 0 (4.8)
pour chaque u. Nous obtenons l’information la meilleure, si nous posons pour u la valeur pour
laquelle a11 u2 + 2a21 u it est minimale, i.e., où la dérivée 2a11 u + 2a21 = 0, i.e., u = −a21 /a11 .
Ainsi (4.8) devient, par (4.6) et (4.5),
a221 ℓ221 ℓ211
a22 − = a22 − 2 = a22 − (ℓ21 )2 > 0 .
a11 ℓ11
3
Sinon, la matrice est définie semi-positive.
4
Le “Commandant Cholesky” (1875–1918) entra à l’École Polytechnique à l’âge de vingt ans et en sortit dans
l’arme de l’Artillerie. Affecté à la Section de Géodésie du Service géographique, en juin 1905, il s’y fit remarquer
de suite par une intelligence hors ligne, une grande facilité pour les travaux mathématiques, un esprit chercheur,
des idées originales, parfois même paradoxales, mais toujours empreintes d’une grande élévation de sentiments et
qu’il soutenait avec une extrême chaleur. (...) Cholesky aborda ce problème en apportant dans ses solutions, ... une
originalité marquée. Il imagina pour la résolution des équations de condition par la méthode des moindres carrés un
procédé de calcul très ingénieux ... (copié du Bulletin géodésique No. 1, 1922).
Pas 2b. Calculons dans (4.4) les valeurs de ai2 pour i = 3, 4, .... On obtient
ai2 = ℓi1 · ℓ21 + ℓi2 · ℓ22 donc ℓi2 = (ai2 − ℓi1 · ℓ21 )/ℓ22 . (4.9)
La division par ℓ22 ne pose pas de problème, car ℓ22 > 0.

Pas 3a. Pour la valeur de a33 dans (4.4) on obtient
q
2 2 2
a33 = (ℓ31 ) + (ℓ32 ) + (ℓ33 ) donc ℓ33 = a33 − (ℓ31 )2 − (ℓ32 )2 . (4.10)
QUESTION. Est-on sûr que a33 − (ℓ31 )2 − (ℓ32 )2 > 0 ? Cette fois-ci on va poser dans (4.3) le
vecteur x = (u, v, 1, 0, ...)T , i.e.,
  
a11 a12 a13 u
  
(u v 1 )  a21 a22 a23   v  = a11 u2 + 2a21 uv + . . . + a33 > 0 (4.11)
a31 a32 a33 1
pour tout u et v. On va de nouveau chercher la valeur minimale de cette expression quadratique.
Pour ne pas nous perdre dans les calculs, observons que
    
ℓ11 ℓ11 ℓ21 ℓ31 a11 a12 a13
    
 ℓ21 ℓ22  ℓ22 ℓ32  =  a21 a22 a23  . (4.12)
ℓ31 ℓ32 0 0 a31 a32 ℓ231 + ℓ232
Ainsi, l’expression (4.11) est égale à
  
ℓ11 ℓ21 ℓ31 u
  
y T y + a33 − ℓ231 − ℓ232 > 0 avec y= ℓ22 ℓ32   v  . (4.13)
0 1
Pour ℓ11 u + ℓ21 v = −ℓ31 et ℓ22 v = −ℓ32 nous avons y = 0 et (4.13) devient l’estimation
recherchée.
Algorithme de Cholesky. On continue anisi avec les pas 3b, 4a, 4b, etc., et on obtient l’algorithme
suivant :
for k := 1 q to n do
Pk−1 2
ℓkk := akk − j=1 ℓkj ;
for i := k + 1 to n do
Pk−1
ℓik := (aik − j=1 ℓij ℓkj )/ℓkk .
Coût de cet algorithme. En négligeant les n racines, le nombre d’opérations nécessaires est
d’environ n Z n
X n3
(n − k) · k ≈ (n − x)x dx = .
k=1 0 6
L’algorithme est deux fois plus rapide que la décomposition LR de Gauss.
Solution du système linéraire. Pour résoudre le système Ax = b, on calcule d’abord la décomposition
de Cholesky A = LLT . Puis
T
LL
| {zx} = b ⇒ résoudre successivement les systèmes Lc = b et LT x = c
c
dont les matrices sont triangulaires.
IV.5 Systèmes surdéterminés – méthode des moindres carrés

Considérons un système d’équations linéaires
a11 x1 + a12 x2 + . . . + a1n xn = b1
a21 x1 + a22 x2 + . . . + a2n xn = b2
.. .. .. .. (5.1)
. . . .
am1 x1 + am2 x2 + . . . + amn xn = bm
où m ≥ n (matriciellement: Ax = b avec x ∈ IRn et b ∈ IRm ; A est une matrice m × n).
Evidemment, le système (5.1) ne possède, en général, pas de solution. L’idée est de chercher un
vecteur x tel que
kAx − bk2 → min (5.2)
pour la norme euclidienne. Une justification probabiliste de cette condition sera donnée dans le
paragraphe IV.7. Le nom “méthode des moindres carrés” indique le choix de la norme dans (5.2)
(la somme des carrés des erreurs doit être minimale).
Théorème 5.1 Soit A une matrice m × n (avec m ≥ n) et soit b ∈ IRm . Le vecteur x est solution
de (5.2) si et seulement si
AT Ax = AT b. (5.3)
Les équations du système (5.3) s’appellent “équations normales”.
Démonstration. Les minima de la fonction quadratique
f (x) := kAx − bk2 = (Ax − b)T (Ax − b) = xT AT Ax − 2xT AT b + bT b
sont donnés par 0 = f ′ (x) = 2(xT AT A − bT A).
Interprétation géométrique. L’ensemble E = {Ax | x ∈ IRn } est un sous-espace linéaire de IRm .

Pour un b ∈ IRm arbitraire, x est une solution de (5.2) si et seulement si Ax est la projection
orthogonale de b sur E. Ceci signifie que Ax − b ⊥ Az pour tout z ∈ IRn . On en déduit que
AT (Ax − b) = 0 et on a ainsi établi une deuxième démonstration de (5.3).
Exemple 5.2 Pour étudier le phénomène de la thermo-électricité, on fait l’expérience suivante. On

soude un fil de cuivre avec un fil de constantan de manière à obtenir une boucle fermée. Un point
de soudure est maintenu à température fixe (T0 ≈ 24◦ C), alors que l’on fait varier la température T
de l’autre. Ceci génère une tension U, laquelle est mesurée en fonction de T (voir le tableau IV.2
et la fig. IV.4). Les données du tableau IV.2 sont prises du livre de P.R. Bevington5.
On suppose que cette dépendance obéit à la loi
U = a + bT + cT 2 (5.4)
et on cherche à déterminer les paramètres a, b et c. Les données du tableau IV.2 nous conduisent
au système surdéterminé (n = 3, m = 21)
Ui = a + bTi + cTi2 , i = 1, . . . , 21. (5.5)
En résolvant les équations normales (5.3) pour ce problème, on obtient a = −0.886, b = 0.0352 et
c = 0.598 · 10−4 . Avec ces paramètres, la fonction (5.4) est dessinée dans la fig. IV.4. On observe
une très bonne concordance avec les données.
TAB . IV.2: Tensions mesurées en fonction de la température T
i Ti◦ C Ui i Ti◦ C Ui i Ti◦ C Ui

1 0 −0.89 8 35 0.42 15 70 1.88
2 5 −0.69 9 40 0.61 16 75 2.10
3 10 −0.53 10 45 0.82 17 80 2.31
4 15 −0.34 11 50 1.03 18 85 2.54
5 20 −0.15 12 55 1.22 19 90 2.78
6 25 0.02 13 60 1.45 20 95 3.00
7 30 0.20 14 65 1.68 21 100 3.22
3
U
1
T0 T
0
T
0 50 100
F IG . IV.4: Tension en fonction de la température et schéma de l’expérience
Remarque. Les équations normales (5.3) possèdent toujours au moins une solution (la projection
sur E existe toujours). La matrice AT A est symétrique et non-négative (xT AT Ax = kAxk2 ≥ 0).
Elle est définie positive si les colonnes de A sont linéairement indépendantes (Ax 6= 0 pour x 6= 0).
Dans cette situation, on peut appliquer l’algorithme de Cholesky pour résoudre le système (5.3).
Mais, souvent, il est préférable de calculer la solution directement de (5.2) sans passer par les
équations normales (5.3).
IV.6 Décomposition QR d’une matrice

Dans l’élimination de Gauss, on a multiplié l’équation Ax = b par la matrice triangulaire Ln−1 ·
. . . · L2 · L1 . De cette manière, on a réduit le problème original à Rx = c où R est une matrice
triangulaire supérieure. Malheureusement, la multiplication de Ax − b avec Li ne conserve pas la
norme du vecteur.
Pour résoudre (5.2), nous cherchons une matrice orthogonale Q telle que

T R′ c′
Q (Ax − b) = Rx − c = x− (6.1)
0 c′′
où R′ (une matrice carrée de dimension n) est triangulaire supérieure et (c′ , c′′ )T est la partition
de c = QT b telle que c′ ∈ IRn et c′′ ∈ IRm−n . Comme le produit par une matrice orthogonale ne
5
P.R. Bevington (1969): Data reduction and error analysis for the physical sciences. McGraw-Hill Book Com-
pany).
change pas la norme du vecteur, on a
kAx − bk22 = kQT (Ax − b)k22 = kRx − ck22 = kR′ x − c′ k22 + kc′′ k22 . (6.2)
On obtient alors la solution de (5.2) en résolvant le système
R ′ x = c′ . (6.3)
Le problème consiste à calculer une matrice orthogonale Q (c.-à-d., QT Q = I) et une matrice

triangulaire supérieure R telles que QT A = R ou de façon équivalente
A = QR. (6.4)
Cette factorisation s’appelle la “décomposition QR” de la matrice A. Pour arriver à ce but, on peut
se servir des rotations de Givens (voir exercice ?? du chapitre V) ou des réflexions de Householder.
Réflexions de Householder (1958). Une matrice de la forme
H = I − 2uuT où uT u = 1 (6.5)

Hx
a les propriétés suivantes :
u
• H est une réflexion à l’hyper-plan {x | uT x = 0} x
car Hx = x − u · (2uT x) et Hx + x ⊥ u .
• H est symétrique.
• H est orthogonale, car
H T H = (I − 2uuT )T (I − 2uuT ) = I − 4uuT + 4uuT uuT = I.
En multipliant A avec des matrices de Householder, nous allons essayer de transformer A en une
matrice de forme triangulaire.
L’algorithme de Householder - Businger - Golub. Dans une première étape, on cherche une
matrice H1 = I − 2u1 uT1 (u1 ∈ IRm et uT1 u1 = 1) telle que
 
α1 × ··· ×
 0 × ··· ×
 
H1 A =  .. .. ..  . (6.6)
 . . . 
0 × ··· ×
Si l’on dénote par A1 la première colonne de A, il faut que H1 A1 = α1 e1 = (α1 , 0, . . . , 0)T et on

obtient |α1 | = kH1 A1 k2 = kA1 k2 . La forme particulière de H1 implique que
H1 A1 = A1 − 2u1 · uT1 A1 = α1 e1 .
L’expression uT1 A1 est un scalaire. Par conséquent,
u1 = C · v1 où v1 = A1 − α1 e1 (6.7)
et la constante C est déterminée par ku1 k2 = 1. Comme on a encore la liberté de choisir le signe
de α1 , posons
α1 = −sign(a11 ) · kA1 k2 (6.8)
pour éviter une soustraction mal conditionnée dans le calcul de v1 = A1 − α1 e1 .

Calcul de H1 A. Notons par Aj et (H1 A)j les j èmes colonnes de A et H1 A respectivement. Alors,
on a
2
(H1 A)j = Aj − 2u1 uT1 Aj = Aj − β · v1T Aj · v1 où β= T . (6.9)
v1 v1
Le facteur β peut être calculé à l’aide de
v1T v1 1
β −1 = = AT1 A1 − 2α1 a11 + α12 = −α1 (a11 − α1 ). (6.10)
2 2
Dans une deuxième étape, on applique la procédure précédente à la sous-matrice de dimension

(m − 1) × (n − 1) de (6.6). Ceci donne un vecteur ū2 ∈ IRm−1 et une matrice de Householder
H̄2 = I − 2ū2 ūT2 . En posant u2 = (0, ū2)T , une multiplication de (6.6) par la matrice H2 =
I − 2u2uT2 donne
     
α1 × · · · × α1 × · · · × α1 × × · · · ×
     0 α × ··· ×
 0   0   2 
H2 H1 A = H2  ..  =  .. =
 0 0 × ··· ×
.
 . C   . H̄2 C   : : : :
0 0 0 0 × ··· ×
En continuant cette procédure, on obtient après n étapes (après n − 1 étapes si m = n) une matrice
triangulaire
′
R
H · . . . · H2 H1 A = R = .
| n {z } 0
QT
Ceci donne la décomposition (6.4) avec QT = Hn · . . . · H2 H1 .
Coût de la décomposition QR. La première étape exige le calcul de α1 par la formule (6.8) (≈ m
opérations), le calcul de 2/v1T v1 par la formule (6.10) (travail négligeable) et le calcul de (H1 A)j
pour j = 2, . . . , n par la formule (6.9) (≈ (n − 1) · 2 · m opérations). En tout, cette étape nécessite
environ 2mn opérations. Pour la décomposition QR, on a alors besoin de
2(n2 + (n − 1)2 + . . . + 1) ≈ 2n3 /3 opérations si m = n (matrice carrée);
2m(n + (n − 1) + . . . + 1) ≈ mn2 opérations si m ≫ n.
En comparant encore ce travail avec celui de la résolution des équations normales (≈ mn2 /2
opérations pour le calcul de AT A et ≈ n3 /6 opérations pour la décomposition de Cholesky de
AT A), on voit que la décomposition QR coûte au pire le double.
Remarque. Si les colonnes de la matrice A sont linéairement indépendantes, tous les αi sont
non nuls et l’algorithme de Householder–Businger–Golub est applicable. Une petite modification
(échange des colonnes de A) permet de traiter aussi le cas général.
Concernant la programmation, il est important de ne calculer ni les matrices Hi , ni la matrice
Q. On retient simplement les valeurs αi et les vecteurs vi (pour i = 1, . . . , n) qui contiennent déjà
toutes les informations nécessaires pour la décomposition. Comme pour l’élimination de Gauss,
on écrit deux sous-programmes. DECQR fournit la décomposition QR de la matrice A (c.-à-d.
les αi , vi et la matrice R). Le sous-programme SOLQR calcule QT b et la solution du système
triangulaire R′ x = c′ (voir (6.3)). Le calcul de QT b = Hn · . . . · H2 H1 b se fait avec une formule
analogue à (6.9).
Exemple 6.1 Si les colonnes de A sont “presque” linéairement dépendantes, la résolution du

problème (5.2) à l’aide de la décomposition QR est préférable à celle des équations normales.
Considérons, par exemple,    
1 1 1
   
A =  ǫ 0, b = 0
0 ǫ 0
où ǫ est une petite constante, disons ǫ2 < eps. Avec un calcul exact, on obtient

1 + ǫ2 1 1
AT A = , AT b =
1 1 + ǫ2 1
et la solution est donnée par
1 1
x1 = x2 = 2
= + O(ǫ2 ).
2+ǫ 2
Un calcul en virgule flottante fait disparaı̂tre le ǫ2 dans AT A et cette matrice devient singulière. On
n’obtient pas de solution.
Par contre, l’algorithme de Householder–Businger–Golub donne (en négligeant ǫ2 ) α1 = −1,
v1 = (2, ǫ, 0)T ,. . . et à la fin
   
−1 √−1 −1
√ 
  T 
R= 0 2 · ǫ, Q b =  ǫ/ √2  .
0 0 −ǫ/ 2
La résolution de (6.3) donne une bonne approximation de la solution exacte.
Calcul pour l’exemple 5.2 (dont les données sont multipliées par 100 ; voir équation (7.11) ci-
dessous) :
   
100 0 0 −88 −457 −22912 −1565712 −494
 100 500 2500 −68   0 −3603 −277963 −141 
   
   


100 1000 10000 −52 



0 −3103 −270463 −125  
 100 1500 22500 −33   0 −2603 −257963 −106 
   
  
 100 2000 40000 −14   0 −2103 −240463 −87  
   
 100 2500 62500 2   0 −1603 −217963 −70 
   
 100 3000 90000 20   0 −1103 −190463 −52 
   
   

 100 3500 122500 42 


 0 −603 −157963 −30  
 100 4000 160000 61   0 −103 −120463 −11 
   
   


100 4500 202500 82 



0 396 −77963 9 

 100 5000 250000 103   0 896 −30463 30 
   
  
 100 5500 302500 122   0 1396 22036 49 

   

 100 6000 360000 145 


 0 1896 79536 72 

 100 6500 422500 168   0 2396 142036 95 
   
   


100 7000 490000 188 



0 2896 209536 115 

 100 7500 562500 210   0 3396 282036 137 
   
  
 100 8000 640000 231   0 3896 359536 158 

   
 100 8500 722500 254   0 4396 442036 181 
   
 100 9000 810000 278   0 4896 529536 205 
   
   
 100 9500 902500 300   0 5396 622036 227 
100 10000 1000000 322 0 5896 719536 249
   
−457 −22912 −1565712 −494 −457 −22912 −1565712 −494
 0 13874 1387444 572   0 13874 1387444 572 
   
   


0 0 25324 1 



0 0 −374437 −21  
 0 0 −9816 0   0 0 0 0 
   

 0 0 −39957 −1 

 0 0 0 1 

   
 0 0 −65098 −4   0 0 0 0 
   
 0 0 −85238 −7   0 0 0 −2 
   
   

 0 0 −100379 −5 

 0 0 0 0 

 0 0 −110520 −6   0 0 0 0 
   
   


0 0 −115661 −6 


0 0 0 0 

 0 0 −115802 −5   0 0 0 1 
   

 0 0 −110943 −7 

 0 0 0 0 

   

 0 0 −101083 −4 

 0 0 0 1 

 0 0 −86224 −2   0 0 0 3 
   
   


0 0 −66365 −2 


0 0 0 1 

 0 0 −41506 0   0 0 0 1 
   

 0 0 −11647 0 


 0 0 0 0 

   
 0 0 23212 2   0 0 0 0 
   
 0 0 63071 5   0 0 0 2 
   
   
 0 0 107930 7   0 0 0 1 
0 0 157789 9 0 0 0 0
IV.7 Etude de l’erreur de la méthode des moindres carrés

Supposons d’avoir un système surdéterminé
n
X
aij xj = bi , i = 1, . . . , m . (7.1)
j=1
En pratique, les bi sont des mesures légèrement erronées et il est naturel de les considérer comme
des valeurs plus ou moins aléatoires. L’étude de l’erreur de la solution x, obtenue par la méthode
des moindres carrés, se fait alors dans le cadre de la théorie des probabilités.
Rappel sur la théorie des probabilités. Considérons des variables aléatoires X (dites “con-
tinues”) qui sont spécifiées par une fonction de densité f : IR → IR, c.-à-d., la probabilité de
l’événement que la valeur de X se trouve dans l’intervalle [a, b) est donnée par
Z b
P (a ≤ X < b) = f (x) dx (7.2)
a
∞R
avec f (x) ≥ 0 pour x ∈ IR et −∞ f (x) dx = 1.
On appelle espérance (mathématique) de la variable aléatoire X le nombre réel
Z ∞
µX = E(X) = xf (x) dx, (7.3)
−∞
et variance la valeur
Z ∞ Z ∞
2
σX = Var (X) = (x − µX )2 f (x) dx = x2 f (x) dx − µ2X . (7.4)
−∞ −∞
Exemple 7.1 Si une variable aléatoire satisfait (7.2) avec (voir la fig. IV.5)

1 1 x − µ 2
f (x) = √ · exp − (7.5)
2π · σ 2 σ
alors on dit que la variable aléatoire satisfait la loi normale ou la loi de Gauss – Laplace que l’on
symbolise par N(µ, σ 2 ). On vérifie facilement que µ est l’espérance et σ 2 la variance de cette
variable aléatoire.
La loi normale est parmi les plus importantes en probabilités. Une raison est due au “théorème
de la limite centrale” qui implique que les observations pour la plupart des expériences physiques
obéissent à cette loi.
95 %
µ − 2σ µ−σ µ µ+σ µ + 2σ
F IG . IV.5: Fonction de densité pour la loi normale
Rappelons aussi que n variables aléatoires X1 , . . . , Xn sont indépendantes si, pour tout ai , bi ,
on a n Y
P (ai ≤ Xi < bi , i = 1, . . . , n) = P (ai ≤ Xi < bi ). (7.6)
i=1
Lemme 7.2 Soient X et Y deux variables aléatoires indépendantes avec comme fonctions de
densité f (x) et g(y) respectivement et soient α, β ∈ IR avec α 6= 0. Alors, les variables aléatoires
αX + β et X + Y possèdent les fonctions de densité
Z
1 x − β ∞
f et (f ∗ g)(z) = f (z − y)g(y) dy. (7.7)
|α| α −∞
Leur espérance mathématique est
E(αX + β) = αE(X) + β, E(X + Y ) = E(X) + E(Y ) (7.8)
et leur variance satisfait
Var (αX + β) = α2 Var (X), Var (X + Y ) = Var (X) + Var (Y ). (7.9)
Démonstration. La fonction de densité pour la variable aléatoire αX + β découle de (pour α > 0)

a − β Z Z
b − β (b−β)/α b t − β
P (a ≤ αX + β < b) = P ≤X< = f (x) dx = α−1 f dt.
α α (a−β)/α a α
Les propriétés (7.8) et (7.9) pour αX + β en sont une conséquence directe.
Comme X et Y sont supposées indépendantes, on obtient (en posant z = x + y)
ZZ Z b Z ∞
P (a ≤ X + Y < b) = f (x)g(y) dx dy = f (z − y)g(y) dy dz
a≤x+y<b a −∞
et on trouve la fonction de densité pour X + Y . Un calcul direct donne

Z ∞ Z ∞ Z ∞ Z ∞
E(X + Y ) = z f (z − y)g(y) dy dz = (x + y)f (x)g(y) dy dx = E(X) + E(Y )
−∞ −∞ −∞ −∞
et, de façon similaire, on obtient

Z ∞ Z ∞
Var (X + Y ) = z2 f (z − y)g(y) dy dz − µ2X+Y
−∞ −∞
Z ∞ Z ∞
= (x + y)2f (x)g(y) dy dx − (µX + µY )2 = Var (X) + Var (Y ).
−∞ −∞
Remarque. Si X et Y sont deux variables aléatoires indépendantes qui obéissent à la loi normale,
les variables aléatoires αX + β et X + Y obéissent aussi à cette loi (exercice 13).
Retour au problème (7.1). Pour pouvoir estimer l’erreur du résultat numérique x, faisons les
hypothèses suivantes :
H1: La valeur bi est la réalisation d’une épreuve pour une variable aléatoire Bi . On suppose que
les Bi soient indépendantes et qu’elles obéissent à la loi de Gauss–Laplace avec βi comme
espérance et σi2 comme variance (les βi sont inconnus, mais les σi2 sont supposés connus).
H2: Le système surdéterminé (7.1) possède une solution unique si l’on remplace les bi par les
nombres βi , c.-à-d. qu’il existe un vecteur ξ ∈ IRn tel que Aξ = β où β = (β1 , . . . , βm )T .
Une illustration de cette situation est donnée en Fig. IV.6.
2 2
B2 β3
B2
B3 B3
1 1 β2
B1 β1
0 0 B1
x1 x2 x3 x1 x2 x3
F IG . IV.6: Illustration pour les hypothèses H1 et H2 (les probabilités sont repésentées par un dégradé de
gris).
Motivation de la méthode des moindres carrés par “maximum likelihood”. Par l’hypothèse
H1, la probabilité que Bi soit dans l’intervalle [bi , bi + dbi ) avec dbi (infinimement) petit est
1 1 b − β 2
i i
P (bi ≤ Bi < bi + dbi ) ≈ √ · exp − · dbi .
2π · σi 2 σi
Comme les Bi sont indépendants, la formule (7.6) implique que
m 1 b − β 2
1 Y i i
P (bi ≤ Bi < bi + dbi , i = 1, . . . , m) ≈ √
· exp − · dbi (7.10)
i=1 2π · σi 2 σi
P
1X m
bi − βi 2 1X m
bi − nj=1 aij ξj 2
= C · exp − = C · exp − .
2 i=1 σi 2 i=1 σi
Selon une idée de Gauss (1812), la “meilleure” réponse xi pour les ξi (inconnus) est celle pour
laquelle la probabilité (7.10) est maximale (“maximum likelihood”). Alors, on calcule x1 , . . . , xn
de façon à ce que
m n 2
X bi X aij
− · xj → min . (7.11)
i=1 σi j=1 σi
Si l’on remplace bi /σi par bi et aij /σi par aij , la condition (7.11) est équivalente à (5.2). Par la suite,
nous supposerons que cette normalisation soit déjà effectuée (donc, σi = 1 pour i = 1, . . . , n).
Estimation de l’erreur. La solution de (7.11) est donnée par x = (AT A)−1 AT b. La solution
théorique satisfait ξ = (AT A)−1 AT β. Alors,
m
X
x − ξ = (AT A)−1 AT (b − β) ou xi − ξi = αij (bj − βj )
j=1
où αij est l’élément (i, j) de la matrice (AT A)−1 AT . L’idée est de considérer la valeur xi comme
la réalisation d’une variable aléatoire Xi définie par
m
X m
X
Xi = αij Bj ou Xi − ξ i = αij (Bj − βj ). (7.12)
j=1 j=1
Théorème 7.3 Soient B1 , . . . , Bm des variables aléatoires indépendantes avec βi comme espérance
et σi = 1 comme variance. Alors, la variable aléatoire Xi , définie par (7.12), satisfait
E(Xi ) = ξi et Var (Xi ) = ǫii (7.13)
où ǫii est le ième élément de la diagonale de (AT A)−1 .
Remarque. Les autres éléments de (AT A)−1 sont les covariances de Xi et Xj .

Démonstration. La formule (7.8) donne E(Xi ) = ξi . Pour calculer la variance de Xi , nous
utilisons le fait que Var (Bi ) = 1 et la formule (7.9). Ceci donne avec ei = (0, . . . , 0, 1, 0, . . . , 0)T
que m
X
2 2
σX i
= αij = keTi (AT A)−1 AT k22 = eTi (AT A)−1 AT A(AT A)−1 ei = eTi (AT A)−1 ei = ǫii .
j=1
Exemple 7.4 Pour l’expérience sur la thermo-électricité (voir le paragraphe IV.5), on a supposé
que les mesures bi ont été faites avec une précision correspondant à σi = 0.01. Pour le système
surdéterminé (on écrit x1 , x2 , x3 pour a, b, c et bi pour Ui )
1 Ti T2 bi
· x1 + · x2 + i · x3 = , i = 1, . . . , 21
σi σi σi σi
la matrice (AT A)−1 devient

 
0.356 · 10−4 −0.139 · 10−5 0.113 · 10−7
 
(AT A)−1 =  −0.139 · 10−5 0.765 · 10−7 −0.713 · 10−9  (7.14)
0.113 · 10−7 −0.713 · 10−9 0.713 · 10−11
et on obtient
σX1 = 0.60 · 10−2 , σX2 = 0.28 · 10−3 , σX3 = 0.27 · 10−5 .
Ceci implique qu’avec une probabilité de 95%, la solution exacte (si elle existe) satisfait
a = −0.886 ± 0.012, b = 0.0352 ± 0.0006, c = 0.598 · 10−4 ± 0.054 · 10−4 .

Test de confiance du modèle. Etudions encore si les données sont compatibles avec l’hypothèse
H2.
En utilisant la décomposition QR de la matrice A, le problème surdéterminé Ax = b se trans-
forme en (voir (6.1))
′ ′ ′
R c c
x = ′′ où = QT b. (7.15)
0 c c′′
La grandeur de kc′′ k22 est une mesure de la qualité du résultat numérique. Théoriquement, si l’on a
β à la place de b et ξ à la place de x, cette valeur est nulle.
Notons les éléments de la matrice Q par qij . Alors, les éléments du vecteur c = QT b sont
P Pm
donnés par ci = m ′′
j=1 qji bj et ceux du vecteur c satisfont aussi ci = j=1 qji (bj − βj ). Il est alors
naturel de considérer les variables aléatoires
m
X
Ci = qji(Bj − βj ), i = n + 1, . . . , m. (7.16)
j=1
Pm
Le but est d’étudier la fonction de densité de i=n+1 Ci2 .
Lemme 7.5 Soient B1 , . . . , Bm des variables aléatoires indépendantes satisfaisant la loi normale
N(βi , 1). Alors, les variables aléatoires Cn+1 , . . . , Cm , définies par (7.16), sont indépendantes et
satisfont aussi la loi normale avec
E(Ci ) = 0, Var (Ci ) = 1. (7.17)
Démonstration. Pour voir que les Ci sont indépendants, calculons la probabilité P (ai ≤ Ci <
bi , i = n + 1, . . . , m). Notons par S l’ensemble S = {y ∈ IRm | ai ≤ yi < bi , i = n + 1, . . . , m}
et par C et B les vecteurs (C1 , . . . , Cm )T et (B1 , . . . , Bm )T . Alors, on a
P (ai ≤Ci < bi , i = n + 1, . . . , m) = P (C ∈ S) = P (QT (B − β) ∈ S)
ZZ 1X m
(a) 1 2
= P (B − β ∈ Q(S)) = √ exp − y dy1 . . . dym (7.18)
Q(S) ( 2π)m 2 i=1 i
ZZ 1X m m Z z2
(b) 1 Y bi 1
= √ exp − 2
z dz1 . . . dzm = √ exp − i dzi .
S ( 2π)m 2 i=1 i i=n+1 ai 2π 2
L’identité (a) est une conséquence de l’indépendance des Bi et (b) découle de la transformation
P P
y = Qz, car det Q = 1 et i yi2 = i zi2 (la matrice Q est orthogonale). En utilisant Si = {y ∈
IRm | ai ≤ yi < bi }, on déduit de la même manière que
Z z2
1 bi
P (ai ≤ Ci < bi ) = P (C ∈ Si ) = . . . = √ exp − i dzi . (7.19)
ai 2π 2
Une comparaison de (7.18) avec (7.19) démontre l’indépendance de Cn+1 , . . . , Cm (voir la définition
(7.6)).
Le fait que les Ci satisfont la loi normale N(0, 1) est une conséquence de (7.19).
Théorème 7.6 (Pearson) Soient Y1 , . . . , Yn des variables aléatoires indépendantes qui obéissent
à la loi normale N(0, 1). Alors, la fonction de densité de la variable aléatoire
Y12 + Y22 + . . . + Yn2 (7.20)
est donnée par (voir fig. IV.7)
1
fn (x) = · xn/2−1 · e−x/2 (7.21)
2n/2 · Γ(n/2)
pour x > 0 et par fn (x) = 0 pour x ≤ 0 (“loi de χ2 à n degrés de liberté”). L’espérance de cette
variable aléatoire vaut n et sa variance 2n.
.2 n=3
n=8
.1 n = 18 95 %
0 10 20 30 40
F IG . IV.7: Fonction de densité (7.21)
Démonstration. Considérons d’abord le cas n = 1. Pour 0 ≤ a < b, on a

√ √ √ √
P (a ≤ Y12 < b) = P ( a ≤ Y1 < b) + P (− a ≥ Y1 > − b)
Z √b Z b
1 −x2 /2 1 dt
=2 √ √ ·e dx = √ · e−t/2 · √ ,
a 2π a 2π t
√
ce qui démontre (7.21) pour n = 1 car Γ(1/2) = π.
Pour le cas général, nous procédons par récurrence. Nous utilisons le résultat du Lemme 7.2
qui affirme que la fonction de densité de Y12 +. . .+Yn+1
2
est la convolution de celle de Y12 +. . .+Yn2
2
avec celle de Yn+1 . Le calcul
Z
1 x
(fn ∗ f1 )(x) = √ (x − t)−1/2 e−(x−t)/2 tn/2−1 e−t/2 dt
2 · Γ(1/2) · 2n/2 · Γ(n/2) 0
Z
e−x/2 x
=√ (x − t)−1/2 tn/2−1 dt
2 · Γ(1/2) · 2n/2 · Γ(n/2) 0
Z
x(n+1)/2−1 e−x/2 1
=√ (1 − s)1/2 sn/2−1 ds = fn+1 (x)
2 · Γ(1/2) · 2n/2 · Γ(n/2) 0
nous permet de conclure.
Pour les variables aléatoires Ci de (7.16), ce théorème montre que

m
X
Ci2 (7.22)
i=n+1
est une variable aléatoire ayant comme fonction de densité fm−n (x) (on rappelle qu’après normal-
isation, on a σi = 1 pour les variables aléatoires Bi ).
Appliquons ce résultat à l’exemple du paragraphe IV.5 (voir la formulation (7.11)). Dans ce
cas, on a kc′′ k22 = 25.2 et m − n = 18 degrés de liberté. La fig. IV.7 montre que cette valeur de
kc′′ k22 est suffisamment petite pour être probable.
Si l’on avait travaillé avec le modèle plus simple
U = a + bT (7.23)
(à la place de (5.4)) on aurait trouvé kc′′ k22 = 526.3 et m − n = 19. Cette valeur est trop grande
pour être probable. La conclusion est que, pour les données du tableau IV.2, cette “loi” est à
réfuser !

Numi 4

Transféré par

Droits d'auteur :

Formats disponibles

Numi 4

Transféré par

Informations du document

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

Numi 4

Transféré par

Droits d'auteur :

Formats disponibles

Chapitre IV

Systèmes d’Equations Linéaires

Considérons un système d’équations linéaires (aij , bj donnés)

La méthode de Wanner. Prenons l’exemple

5x2 − 7x3 = 2 ⇒ x2 = 6, x3 = 4 et par (4) x1 = −1.

IV.1 Elimination de Gauss

Le système (1.2) contient un sous-système de dimension n − 1 sur lequel on peut répéter la

(A, b) → (A(1) , b(1) ) → (A(2) , b(2) ) → . . . → (A(n−1) , b(n−1) ) =: (R, c)

on obtient un système triangulaire

do ir=1,n-1 C ---- BACK SUBSTITUTION ---

Théorème 1.1 L’élimination de Gauss équivaut à une factorisation

Démonstration. En utilisant les matrices

le premier pas de l’élimination de Gauss correspond à une multiplication de A avec L1 , le deuxième

Algorithme avec recherche de pivot.

Exemple 1.2 (Forsythe) Considérons le système

La solution exacte est, comme on voit, x1 = 1 et x2 = 1. Appliquons l’élimination de Gauss et

x1 = (b1 − a12 x2 )/a11 = (1.00 − 1.00 ∗ 1.00)/(1.00 · 10−4 ) = 0.

Le résultat numérique, obtenu pour x1 , est faux.

d’insérer dans le code ci-dessus, après la première ligne, les commandes :

IV.2 Etude des erreurs ; “Backward Error Analysis”

ℓb21 = ℓ21 + ǫ où |ǫ| ≤ eps

(car |ℓ21 | ≤ 1 à cause du choix du pivot).

Nous transportons à nouveau ces erreurs du côté des données, et le calcul

Probl. bien cond. Probl. bien cond. Probl. mal cond.

(A, b) x (A, b) x (A, b) x

IV.3 La condition d’une matrice

pour kxk∞ = maxi=1,...,n |xi | , on a

TAB . IV.1: Conditions de matrices de Hilbert et Vandermonde

IV.4 L’algorithme de Cholesky

i.e., l’élément aij de A est le produit scalaire des colonnes i et j de B 2 . Alors

A est symétrique (car le prod. scal. est symétrique; ou car AT = B T (B T )T = B T B = A)

La division par ℓ22 ne pose pas de problème, car ℓ22 > 0.

IV.5 Systèmes surdéterminés – méthode des moindres carrés

Démonstration. Les minima de la fonction quadratique

f (x) := kAx − bk2 = (Ax − b)T (Ax − b) = xT AT Ax − 2xT AT b + bT b

sont donnés par 0 = f ′ (x) = 2(xT AT A − bT A).

Interprétation géométrique. L’ensemble E = {Ax | x ∈ IRn } est un sous-espace linéaire de IRm .

Exemple 5.2 Pour étudier le phénomène de la thermo-électricité, on fait l’expérience suivante. On

Ui = a + bTi + cTi2 , i = 1, . . . , 21. (5.5)

TAB . IV.2: Tensions mesurées en fonction de la température T

i Ti◦ C Ui i Ti◦ C Ui i Ti◦ C Ui

IV.6 Décomposition QR d’une matrice

change pas la norme du vecteur, on a

On obtient alors la solution de (5.2) en résolvant le système

Le problème consiste à calculer une matrice orthogonale Q (c.-à-d., QT Q = I) et une matrice

Réflexions de Householder (1958). Une matrice de la forme

H = I − 2uuT où uT u = 1 (6.5)

H T H = (I − 2uuT )T (I − 2uuT ) = I − 4uuT + 4uuT uuT = I.

Si l’on dénote par A1 la première colonne de A, il faut que H1 A1 = α1 e1 = (α1 , 0, . . . , 0)T et on

L’expression uT1 A1 est un scalaire. Par conséquent,

pour éviter une soustraction mal conditionnée dans le calcul de v1 = A1 − α1 e1 .

Dans une deuxième étape, on applique la procédure précédente à la sous-matrice de dimension

Ceci donne la décomposition (6.4) avec QT = Hn · . . . · H2 H1 .

Exemple 6.1 Si les colonnes de A sont “presque” linéairement dépendantes, la résolution du

IV.7 Etude de l’erreur de la méthode des moindres carrés

Leur espérance mathématique est

E(αX + β) = αE(X) + β, E(X + Y ) = E(X) + E(Y ) (7.8)

et leur variance satisfait

Var (αX + β) = α2 Var (X), Var (X + Y ) = Var (X) + Var (Y ). (7.9)