Décompositions LU et Choleski
(Programmation avec Maple)
Préparation à la nouvelle épreuve d’informatique
de l’École Polytechnique.
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Énoncé
Énoncé
Ce document vise a la préparation a la nouvelle épreuve d’informatique de l’École Polytechnique.
Les caractéristiques de cette épreuve peuvent être consultées à l’adresse :
http ://www.enseignement.polytechnique.fr/informatique/concours/
Parmi les langages de programmation possibles, on a choisi Maple, dont on n’a utilisé que les fonc-
tionnalités de base, conformément aux demandes des concepteurs de l’épreuve.
I. La décomposition LU
Soit A une matrice inversible de Mn(K). Une décomposition LU de A est une égalité A = LU ,
ou L est une matrice triangulaire inférieure (L pour “Low”) a diagonale unité (tous les coefficients
diagonaux valent 1), et où U est une matrice triangulaire supérieure (U pour “Up”). La matrice U
est nécessairement inversible (donc ses coefficients diagonaux sont non nuls.)
Par exemple, A =
2 −3
2 −3
−2
4 −9 −2
5
−2
1 −1
2
3
5 −4
=
0 0
0
1
0 0
1
−1
1 0
3
2
−1 −2 1 1
1 −1
2 −3
1
0 −1 −2
2
0
0
2
0 −5
0
0
Pour tout indice k compris entre 1 et n, on notera Ak la matrice extraite de A formée à l’intersection
des k premieres lignes et des k premieres colonnes de A. On notera ∆k le déterminant de Ak. Les ∆k
sont appelés les mineurs principaux de A.
Proposition (Existence et unicité de la décomposition LU )
Une matrice A de GLn(K) possède une décomposition LU si et seulement si ses mineurs princi-
paux sont non nuls. Cette décomposition est alors unique. [Démonstration]
Dans la suite de cette partie, on suppose que la matrice A possède une décomposition LU et on voit
comment mettre en œuvre le calcul des matrices L et U .
On note aij, ‘ij, et uij les termes généraux de A, L, U .
1. Ecrire les égalités donnant aik en fonction des ‘ij (avec j ≤ i) et ujk (avec j ≤ k).
En déduire les expressions :
– De uij, pour i ≤ j, en fonction de aij, de ‘ik (k < i), de ukj (k < i).
– De ‘ij, pour i > j, en fonction de aij, de ‘ik (k < j), de ukj (k ≤ j).
Montrer comment les égalités obtenues permettent de calculer de proche en proche (et on
précisera dans quel ordre) tous les coefficients de L et de U .
[ S ]
2. En déduire une procédure LU1 calculant la décomposition d’une matrice A (supposée carrée
inversible d’ordre n, et possédant une telle décomposition).
La syntaxe sera LU1(A, ‘L‘, ‘U ‘, n) et les deux matrices L et U seront placées dans les noms de
variable homonymes passés en argument.
[ S ]
c(cid:13)EduKlub S.A.
Page 1
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Énoncé
3. Dans un souci d’économie, on décide de former une seule matrice B, carrée d’ordre n, dont le
terme général bij est donné par : bij = uij si i ≤ j, et bij = ‘ij si i > j.
(a) Montrer que B peut être formée de la manière suivante, dans l’ordre des j croissants :
∀i ∈ {1, . . . , n},
bij = aij −
min(i,j)−1
P
k=1
bik bkj suivi de bij ←
bij
bjj
si i > j.
Montrer qu’il est possible, dans le programme, de n’utiliser qu’une variable pour désigner
la matrice initiale A, puis pour former petit à petit la matrice B.
[ S ]
(b) En déduire une procédure LU2 décomposant la matrice A, et plaçant le résultat dans la
matrice A elle-même (les coefficients sur et au-dessus de la diagonale étant ceux de U et
les coefficients sous la diagonale étant ceux de L.) [ S ]
4. Décomposition P A = LU .
On a vu que l’existence de la décomposition A = LU est soumise à des conditions sur la matrice
inversible A (la première d’entre elles étant que a11 doit être non nul.)
Chacune de ces conditions équivaut à la non-nullité d’un coefficient utilisé comme pivot. Si un
tel pivot est nul, la décomposition est impossible, à moins qu’on ne puisse continuer moyennant
un échange de lignes.
Si ce pivot est “presque nul”, on utilisera encore un échange de lignes (sinon la division par le
pivot conduit à une imprécision importante.)
Ces considérations conduisent à envisager la décomposition LU non pas de la matrice A initiale,
mais d’une matrice B = P A, où P est une matrice de permutation.
Soit P un élément de Mn(K), dont le coefficient d’indice (i, j) est noté pij.
On dit que P est une matrice de permutation s’il existe une permutation σ de {1, 2, . . . , n}
telle que pij = δσ(i)j (notations de Kronecker), pour tous indices i et j.
Pour toute matrice A de Mn(K), la matrice B = P A se déduit de A par la permutation σ
appliquée aux lignes de A : pour tout i de {1, . . . , n}, la ligne d’indice i de B = P A est en effet
la ligne d’indice σ(i) de A.
Proposition (Existence de la décomposition P A = LU )
Pour tout matrice A de GLn(K), il existe au moins une matrice de permutation P telle que
la matrice P A admette une décomposition LU .
[Démonstration]
On a vu dans la question précédente que la décomposition LU de A pouvait s’effectuer di-
rectement sur cette matrice. Cette décomposition s’effectue en n étapes successives (une par
colonne.)
Pour tout j de {1, . . . , n}, l’étape j se déroule en deux phases :
– Pour tout i de {1, . . . , n} on réalise l’affectation aij ← aij −
aikakj.
On modifie ainsi tous les coefficients de la j-ème colonne.
min(i,j)−1
P
k=1
Le pivot est maintenant le nouveau coefficient diagonal p = ajj.
c(cid:13)EduKlub S.A.
Page 2
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Énoncé
– Pour tout i de {j + 1, . . . , n} (donc sous la diagonale) on divise aij par le pivot p.
Plaçons nous apres la premiere phase de l’étape j. Notons A0 l’état actuel de la matrice dans
laquelle s’effectue la décomposition, et A l’état initial de cette matrice.
Donnons-nous un indice de ligne i ≥ j.
Notons L1, . . . , Ln les lignes de A et L0
Il n’est pas difficile de se convaincre (en remontant les calculs effectués pour en arriver à cette
phase de la décomposition) que la ligne L0
i ne dépend que des lignes L1, . . . , Lj−1, Li. Notons
par exemple L0
n celles de A0.
1, . . . , L0
i = ϕ(L1, . . . , Lj−1, Li).
On se rend compte également que cette fonction ϕ est indépendante de l’indice i ≥ j.
Cela est valable en particulier pour l’indice j lui-même.
On en déduit que si on veut que soient échangées les lignes L0
il suffit d’échanger les lignes Lj et Li dans la matrice initiale. Il est donc possible de procéder,
a l’issue de la premiere phase de l’étape j, à l’échange de L0
i (où i > j) de
façon a rendre maximum (en valeur absolue) le pivot p qui sera utilisé dans la deuxieme phase
de cette étape.
i (à ce stade de la décomposition)
j avec une ligne L0
j et L0
Il suffit alors de mémoriser les échanges de lignes effectués tout au long de la décomposition. On
aura ainsi obtenu la décomposition LU de la matrice P A, où P est la matrice de permutation
résumant les échanges successifs de lignes.
Si cette méthode échoue, c’est-a-dire si le pivot maximum obtenu a un moment donné est nul,
c’est que A n’a pas de décomposition P A = LU : c’est qu’elle n’est pas inversible.
(a) Écrire une procédure LU décomposant la matrice A, et plaçant le résultat dans la matrice
A elle-même (les coefficients sur et au-dessus de la diagonale étant ceux de U et les
coefficients sous la diagonale étant ceux de L.)
La procédure LU se distinguera de LU2 par la recherche d’un pivot de valeur absolue
maximum sur chaque colonne (suivie éventuellement d’un échange de lignes.)
On utilisera la syntaxe LU(A, ‘s‘, n) : dans cette écriture la variable s recevra la per-
mutation représentant la succession des échanges de lignes, sous la forme d’un tableau
unidimensionnel [s(1), . . . , s(n)].
[ S ]
(b) Montrer comment ce qui précede permet de ramener la résolution d’un systeme linéaire
AX = B (ou A est carrée d’ordre n et inversible, et ou B est un vecteur colonne de
hauteur n) a la résolution de deux systemes triangulaires successifs.
Écrire une fonction syst renvoyant la solution X du système AX = B, sous la forme d’un
vecteur de taille n. La syntaxe sera syst(M, B, s, n), où :
– M est la matrice carrée contenant la décomposition LU .
– B est un vecteur représentant les seconds membres.
– s est un vecteur représentant la permutation P telle que P A = LU .
[ S ]
c(cid:13)EduKlub S.A.
Page 3
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Énoncé
(c) Vérifier sur un exemple numérique de résolution de système (en imposant par exemple une
valeur très faible pour a11 dans A) que la décomposition obtenue par LU peut se révéler
beaucoup plus précise que celle qui est obtenue par LU2.
[ S ]
(d) Écrire une fonction inv inversant une matrice carrée A d’ordre n.
[ S ]
II. La décomposition de Choleski
Les matrices considérées ici sont à coefficients réels. On identifiera un vecteur u de Rn avec la matrice
uni-colonne U de ses coordonnées dans la base canonique (e).
Le produit scalaire canonique de R peut donc s’écrire : < u, v >= tU V =
n
P
k=1
uk vk.
Soit S une matrice carrée d’ordre n, symétrique et à coefficients réels. On sait que S est diagonalisable
sur R, et que ses différents sous-espaces propres sont orthogonaux deux à deux.
En particulier, il existe une base orthonormée (ε) = ε1, . . . , εn formée de vecteurs propres de S.
Proposition (matrices symétriques définies positives)
Soit S une matrice carrée d’ordre n, symétrique et à coefficients réels.
Les conditions suivantes sont équivalentes :
– Toutes les valeurs propres de S sont strictement positives.
– Pour tout vecteur non nul X de Rn, tXSX > 0.
– Il existe une matrice inversible M telle que S = M tM .
Si ces conditions sont réalisées, on dit S est définie positive.
[Démonstration]
La matrice inversible M telle que S = M tM n’est pas unique, à moins de lui en demander “un peu
plus”. C’est l’objet de la proposition suivante, qui introduit la décomposition de Choleski.
Proposition (décomposition de Choleski)
Soit S une matrice symétrique définie positive. Il existe une unique matrice L triangulaire
inférieure à coefficients diagonaux strictement positifs, telle que S = L tL.
Cette écriture est appelée décomposition de Choleski de S. [Démonstration]
1. Soit S = L tL la décomposition de Choleski d’une matrice S définie positive.
On note sij et ‘ij les coefficients d’indice i, j de S et L.
Établir un système d’égalités permettant de calculer les coefficients ‘ij.
[ S ]
2. Écrire une fonction choleski prenant en argument une matrice carrée M à coefficients réels,
supposée définie positive, et renvoyant la matrice L de la décomposition de M .
[ S ]
c(cid:13)EduKlub S.A.
Page 4
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Corrigé
Corrigé
I. La décomposition LU
1. Le terme général aij de la matrice A = LU s’écrit aij =
Or ‘ik = 0 quand k > i et ukj = 0 quand k > j.
La somme est donc limitée à k = min(i, j) : aij =
n
P
k=1
‘ik ukj.
‘ik ukj. On en déduit :
min(i,j)
P
k=1
i−1
P
k=1
i
P
k=1
i−1
P
k=1
Si i ≤ j, alors aij =
‘ik ukj =
‘ik ukj + uij ⇒ uij = aij −
‘ik ukj
(1)
Si i > j : aij =
j
P
k=1
‘ik ukj =
j−1
P
k=1
‘ik ukj + ‘ij ujj ⇒ ‘ij = 1
ujj
(cid:16)
aij −
j−1
P
k=1
(cid:17)
‘ik ukj
(2)
Ces égalités permettent de former L et U , ligne par ligne, ou colonne par colonne.
– Pour j = 1, l’égalité (1) donne u11 = a11 (non nul car a11 = ∆1.)
Puis (2) donne la sous-diagonale de la colonne 1 de L : ∀ i ∈ {2, . . . , n}, ‘i1 =
ai1
u11
.
– Soit j ≥ 2. Supposons connues les colonnes 1 à j −1 de L et U .
Quand i varie de 1 a j, l’égalité (1) donne la j-eme colonne de U .
Quand i varie de j +1 a n, l’égalité (2) donne la j-eme colonne de L.
La division par ujj est possible car par hypothèse la décomposition A = LU existe.
En tout cas, tout nouveau coefficient est obtenu en fonction de coefficients déjà calculés !
On peut donc calculer de proche en proche tous les coefficients de U et L, en procédant de la
colonne 1 a la colonne n, et sur chaque colonne de la ligne 1 a la ligne n.
[ Q ]
local i,j,k,c:
L:=matrix(n,n,0); U:=matrix(n,n,0); # matrices nulles d’ordre n
for j to n do
boucle de formation des matrices L et U
for i to j do
c:=A[i,j];
for k to i-1 do
2. La procédure LU1 est la traduction des calculs précédents :
> LU1:=proc(A::array,L::name,U::name,n::integer)
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
od;
L[j,j]:=1;
for i from j+1 to n do
c:=A[i,j];
for k to j-1 do
od;
U[i,j]:=c;
c:=c-L[i,k]*U[k,j];
c:=c-L[i,k]*U[k,j]
calcul de la colonne j de U
initialise le coefficient d’indice i, j
boucle de calcul de U [i, j]
place le coefficient dans U
coefficient diagonal d’indice j de L
calcul de la colonne j de L
initialise le coefficient d’indice i, j
boucle de calcul de L[i, j]
Publicité
c(cid:13)EduKlub S.A.
Page 5
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Corrigé
od;
od;
L[i,j]:=c/U[j,j];
>
>
>
>
> end:
On teste LU1 sur la matrice A citée comme exemple dans l’énoncé. On affiche ensuite L et U ,
et on vérifie l’égalité A = LU :
place le coefficient dans L
od;
> A:=matrix([[2,-3,1,-1],[-2,2,-3,2],[4,-9,-2,3],[-2,5,5,-4]]):
> LU1(A,’L’,’U’,4): eval(L), eval(U), evalm(L&*U);
0 0
0
1
0 0
1
−1
2
1 0
3
−1 −2 1 1
,
1 −1
2 −3
1
0 −1 −2
2
0
0
2
0 −5
0
0
,
2 −3
2 −3
−2
4 −9 −2
5
−2
1 −1
2
3
5 −4
[ Q ]
3. Décomposition “in situ”.
(a) Dans les égalités (1) et (2) de la première question, tous les coefficients de L et U peuvent
être renommés à l’aide de la convention suivante :
Pour tous indices i, j notons bij = uij si i ≤ j et bij = ‘ij si i > j.
On en déduit que les égalités (1) et (2) peuvent être réécrites sous la forme :
j−1
P
k=1
(2) Si i > j, alors bij = 1
bjj
(1) Si i ≤ j, alors bij = aij −
i−1
P
k=1
bik bkj
aij −
(cid:16)
Il revient au même d’écrire, dans l’ordre croissant des indices de colonne j :
(cid:17)
bik bkj
Pour tout i de {1, . . . , n}, bij = aij −
bik bkj, suivi de bij ←
min(i,j)−1
P
k=1
bij
bii
si i > j.
Chaque aij n’est utilisé que lors du calcul de bij. On peut donc remplacer au fur et à
mesure les aij par les bij et ainsi n’utiliser que la variable A.
[ Q ]
(b) Voici la procédure LU2, qui se déduit très simplement des considérations précédentes :
A[i,j]:=A[i,j]-A[i,k]*A[k,j]
for k to min(i,j)-1 do
for i to n do
> LU2:=proc(A::array,n)
local i,j,k,p;
>
for j to n do
>
>
>
>
>
>
>
>
>
>
>
> end:
od;
fi;
od;
A[i,j]:=A[i,j]/p
od;
p:= A[j,j];
for i from j+1 to n do
calcul de la colonne j
calcul du terme d’indice i, j
c’est le pivot de la colonne j
modification des termes sous-diagonaux
c(cid:13)EduKlub S.A.
Page 6
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Corrigé
On décompose la même matrice. Le résultat, dans A, est une superposition de L, U comme
on le voit en rappelant ces matrices obtenues avec LU1 :
> LU2(A,4): eval(A), eval(L), eval(U);
1 −1
2 −3
1
−1 −1 −2
2
2
2
3
1 −5
−1 −2
,
0 0
0
1
0 0
1
−1
2
1 0
3
−1 −2 1 1
,
1 −1
2 −3
1
0 −1 −2
2
0
0
2
0 −5
0
0
[ Q ]
4. Décomposition P A = LU .
(a) Voici la procédure LU :
od;
for i to n do
valeur et ligne initiales du pivot
boucle de recherche du meilleur pivot
for k to min(i,j)-1 do
A[i,j]:=A[i,j]-A[i,k]*A[k,j]
if abs(A[k,j])>abs(p) then p:=A[k,j]; i:=k; fi;
od;
p:=A[j,j]; i:=j;
for k from j+1 to n do
par défaut, s = permutation identité
de la colonne j = 1 à la colonne j = n
de haut en bas sur la colonne j
boucle de calcul de A[i, j]
local i,j,k,p,q,t;
s:=array(1..n);
for i to n do s[i]:=i od;
for j to n do
> LU:=proc(A::array,s::name,n::integer)
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> end:
On reprend la matrice utilisée pour illustrer le fonctionnement de LU1 et LU2.
od;
t:=s[j]; s[j]:=s[i]; s[i]:=t; # mémorise l’échange des lignes
fi;
for k from j+1 to n do A[k,j]:=A[k,j]/p od; # divise par le pivot
si le pivot n’est pas diagonal
boucle d’échange des lignes j et i
t:=A[j,k]: A[j,k]:=A[i,k]: A[i,k]:=t
si pivot=0, matrice non inversible
ERROR("matrice non inversible")
fi;
if i>j then
od;
if p=0 then
for k to n do
od;
On la décompose avec LU, et on affiche la permutation des lignes, et le résultat de la
décomposition (L et U sont superposées dans la variable initiale A.)
> A:=matrix([[2,-3,1,-1],[-2,2,-3,2],[4,-9,-2,3],[-2,5,5,-4]]):
> LU(A,’s’,4); eval(s),eval(A);
c(cid:13)EduKlub S.A.
Page 7
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Corrigé
[3, 2, 4, 1],
4
3
−2
−9
7/2
−1/2 −5/2 −4
−1/2 −1/5
16/5 −9/5
1/2 −3/5 −1/8 −5/8
Ce résultat montre que ce n’est pas A qui a été décomposée, mais B = P A, avec :
B = P A =
0 0 1 0
0 1 0 0
0 0 0 1
1 0 0 0
2 −3
2 −3
−2
4 −9 −2
5
−2
1 −1
2
3
5 −4
=
4 −9 −2
2 −3
−2
−2
Publicité
5
2 −3
3
2
5 −4
1 −1
Cela est confirmé en appliquant la procédure LU1 à la matrice B :
> B:=matrix([[4,-9,-2,3],[-2,2,-3,2],[-2,5,5,-4],[2,-3,1,-1]]):
> LU1(B,’L’,’U’,4); eval(L), eval(U);
[ Q ]
0
1
0
1
0
−1/2
−1/2 −1/5
0
1/2 −3/5 −1/8 1
0
0
1
,
−2
4 −9
0 −5/2 −4
0
0
3
7/2
16/5 −9/5
−5/8
0
0
0
(b) Il existe une matrice de permutation P , et deux matrices L, U telles que P A = LU .
Le systeme AX = B est alors équivalent a P AX = P B, c’est-a-dire a LU X = P B.
Si P est associée à la permutation σ, alors B0 = P B = (bσ(1), bσ(2), . . . , bσ(n)).
Dans ces conditions, trouver X c’est résoudre LY = B0 puis U X = Y .
Ainsi le systeme initial se ramene a deux systemes triangulaires successifs.
Posons X = (x1, x2 . . . , xn) et Y = (y1, y2, . . . , yn). Le système LY = B0 s’écrit :
1
‘21
0
1
0
0
‘31
...
‘n1
‘32
...
‘n2
1
. . .
. . .
. . .
. . .
. . .
. . .
‘n,n−1
0
0
0
0
1
=
y1
y2
y3
...
yn
bσ(1)
bσ(2)
bσ(3)
...
bσ(n)
⇔
y1 = bσ(1)
y2 = bσ(2) − ‘21y1
y3 = bσ(3) − ‘31y1 − ‘32y2
. . . = . . .
yn = bσ(n) −
‘njyj
n−1
P
j=1
Ensuite le système U X = Y s’écrit :
u11 u12 u13
u22 u23
0
0
0
u33
...
...
. . .
. . .
0
0
xn = yn/unn
. . . u1n
. . . u2n
. . . u3n
...
. . .
unn
0
x1
x2
x3
...
xn
=
y1
y2
y3
...
yn
⇔
xn−1 =
. . . = . . .
1
u11
x1 =
yn−1 − un−1,nyn
un−1,n−1
(cid:16)
y1 −
n
P
j=2
(cid:17)
u1j xj
Voici le listing de syst. La variable X désigne successivement les seconds membres (après
permutation), la solution Y du premier système, puis la solution finale.
c(cid:13)EduKlub S.A.
Page 8
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Corrigé
t:=0; for j to i-1 do t:=t+A[i,j]*X[j] od; X[i]:=X[i]-t;
initialisation
permutation sur le second membre
résolution du système LY = B
local X,i,j,t;
X:=array(1..n);
for i to n do X[i]:=B[s[i]] od;
for i from 2 to n do
> syst:=proc(A::array,B::array,s::array,n::integer)
>
>
>
>
>
>
>
>
>
>
> end:
On reprend ici une matrice plusieurs fois utilisée avec les seconds membres (1, 6, 5, 2).
t:=0; for j from i+1 to n do t:=t+A[i,j]*X[j] od;
X[i]:=(X[i]-t)/A[i,i];
od;
for i from n to 1 by -1 do
résolution du système U X = Y
renvoie la solution du système
od; eval(X);
On constate que la solution du système est X = (−17, −10, −2, −7).
> A:=matrix([[2,-3,1,-1],[-2,2,-3,2],[4,-9,-2,3],[-2,5,5,-4]]):
> AA:=copy(A): B:=array([1,6,5,2]):
> LU(A,’s’,4); X:=syst(A,B,s,4);
X := [−17, −10, −2, −7]
On vérifie maintenant que la solution est correcte. Le produit AX redonne bien B.
> evalm(AA&*X);
[1, 6, 5, 2]
[ Q ]
(c) On reprend l’exemple précédent, mais on remplace A[1, 1] = 2 par
1
1000000
> A:=matrix([[2,-3,1,-1],[-2,2,-3,2],[4,-9,-2,3],[-2,5,5,-4]]):
> A[1,1]:=1/1000000: A1:=copy(A): A2:=copy(A):
> B:=array([1,6,5,2]):
On résout le systeme AX = B apres LU2 et LU de façon exacte.
.
On constate (c’est normal) que les résultats sont les mêmes.
On calcule aussi une valeur approchée de la solution exacte (12 chiffres significatifs.)
> LU(A1,’s’,4); X:=syst(A1,B,s,4);
> LU2(A2,4); X:=syst(A2,B,array([1,2,3,4]),4);
> map(evalf,X);
X :=
h
−
340000000
61999979
, −
−76000062
61999979
, −
487999736
61999979
, −
653999603
61999979
X :=
h
−
340000000
61999979
, −
−76000062
61999979
, −
487999736
61999979
, −
653999603
61999979
i
i
[−5.483872825, −1.225807867, 7.870966150, 10.54838427]
c(cid:13)EduKlub S.A.
Page 9
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Corrigé
On met maintenant tous les coefficients de A au format “virgule flottante”.
On calcule la solution du systeme apres LU puis après LU2.
On constate des erreurs d’arrondi assez sensibles dans le premier cas, alors que la fonction
LU donne un résultat extrêmement précis. On remarque que dans ce dernier cas, on a dû
Publicité
préciser que la permutation s est égale à l’identité.
> A:=map(evalf,A): A1:=copy(A): A2:=copy(A):
> LU(A1,’s’,4); X:=syst(A1,B,s,4);
> LU2(A2,4); X:=syst(A2,B,array([1,2,3,4]),4);
X := [−5.483872825, −1.225807867, 7.870966153, 10.54838427]
X := [−5.490000000, −1.225966469, 7.872628719, 10.55052264]
[ Q ]
(d) Voici la fonction inv. L’idée est évidemment de résoudre le systeme AX = B, ou les
second membres B sont successivement l’une des n colonnes de la matrice identité.
local B,M,AA,s,i,j,t;
M:=array(1..n,1..n); AA:=copy(A); # initialisations
LU(AA,’s’,n);
for j to n do
> inv:=proc(A::array,n::integer)
>
>
>
>
>
>
>
>
>
> end:
B:=vector(n,0); B[j]:=1;
B:=syst(AA,B,s,n);
for i to n do M[i,j]:=B[i] od;
od;
RETURN(eval(M));
décomposition P A = LU
de la colonne 1 à la colonne n
j-ième colonne de l’identité
résolution du système
j-ième colonne de l’inverse de A
renvoie la matrice inverse de M
On reprend la matrice A ayant servi d’exemple jusqu’ici, et on calcule son inverse avec
notre fonction inv, puis avec la fonction intégrée de Maple. On constate que les résultats
sont identiques.
> A:=matrix([[2,-3,1,-1],[-2,2,-3,2],[4,-9,-2,3],[-2,5,5,-4]]):
> inv(A,4), evalm(A&^(-1));
[ Q ]
−
−
−
−
21
20
4
5
9
10
8
5
−
7
4
−
1
2
−1
−1 −
−1 −
−
13
20
2
5
3
10
1
5
−
−
11
10
3
5
1
5
1
5
,
−
−
−
−
21
20
4
5
9
10
8
5
−
−
7
4
−
1
2
−1
−
13
20
2
5
3
10
1
5
−
−
11
10
3
5
1
5
1
5
−
c(cid:13)EduKlub S.A.
Page 10
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Corrigé
II. La décomposition de Choleski
1. L’égalité L tL = S permet par identification de calculer les ‘i,j :
n
P
k=1
S = L tL ⇔ ∀ i, j ∈ {1, . . . , n}, si,j =
[L]i,k[ tL]k,j =
n
P
k=1
‘i,k‘j,k =
min(i,j)
P
k=1
‘i,k‘j,k
La somme est limitée supérieurement à min(i, j) car L est triangulaire inférieure.
D’autre part, on sait que la matrice S est symétrique. L’identification dans S = L tL peut
donc se limiter, sans perdre aucune généralitén, au cas i ≥ j (coefficients diagonaux et sub-
diagonaux). Pour tout entier j fixé (compris entre 1 et n) on trouve :
– (1) : Pour i = j :
j
P
k=1
j,k = sj,j ⇒ ‘ 2
‘ 2
j,j = sj,j −
j−1
P
k=1
‘ 2
j,k
– (2) : Pour j + 1 ≤ i ≤ n :
j
P
k=1
‘i,k‘j,k = si,j ⇒ ‘i,j =
(cid:16)
1
‘j,j
si,j −
‘i,k‘j,k
(cid:17)
.
j−1
P
k=1
Les coefficients de L peuvent donc être calculés colonne par colonne.
Plus précisément, lors du calcul de la j-ème colonne de L :
– La phase (1) donne le j-ème coefficient diagonal ‘j,j > 0 (si le second membre n’est pas
strictement positif, c’est que S n’est pas définie positive). On constate que ‘j,j est obtenu
en fonction de sj,j et de coefficients déja calculés dans L (sur la même ligne, mais sur les
colonnes précédentes).
– La phase (2) donne les termes subdiagonaux de la j-ème colonne de L. Chaque ‘i,j est obtenu
en fonction de si,j (en même position dans S), du coefficient diagonal ‘j,j (qui vient d’être
calculé) et de termes déjà connus dans L (sur les lignes i, j, mais sur des colonnes d’indice
inférieur à j).
[ Q ]
2. Voici la fonction choleski
t:=t-L[j,k]^2;
local L,i,j,k,t:
L:=matrix(n,n,0);
for j to n do
t:=S[j,j];
for k to j-1 do
> choleski:=proc(S::array,n::integer)
>
>
>
>
>
>
>
>
>
>
>
>
>
>
fi;
t:=sqrt(t); L[j,j]:=t;
for i from j+1 to n do
t:=S[i,j];
for k to j-1 do
od;
if t<=0 then
ERROR("matrice non def+")
matrice nulle d’ordre n
calcul de L, colonne par colonne
j-ème coefficient diagonal de S
calcule le j-ème coefficient diagonal de L
erreur si S n’est pas définie positive
place le j-ème coefficient diagonal de L
coefficients sous-diagonaux, colonne j
c(cid:13)EduKlub S.A.
Page 11
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Corrigé
t:=t-L[i,k]*L[j,k];
od;
L[i,j]:=t/L[j,j];
od;
place le coefficient d’indice (i, j) dans L
od;
RETURN(eval(L));
renvoie la matrice L
>
>
>
>
>
>
> end:
On crée une matrice symétrique S :
S:=matrix([[1,-2,4],[-2,13,-11],[4,-11,21]]);
S :=
4
1
−2
4 −11
−2
13 −11
21
On décompose la matrice S puis on vérifie que le résultat obtenu est correct :
> L:=choleski(S,3): eval(L),evalm(L&*linalg[transpose](L));
0
0
1
0
3
−2
4 −1 2
,
[ Q ]
4
1
−2
4 −11
−2
13 −11
21
c(cid:13)EduKlub S.A.
Page 12
Tous droits de l’auteur des œuvres réservés. Sauf autorisation, la reproduction ainsi que toute utilisation des œuvres autre que la consultation
individuelle et privée sont interdites.
www.klubprepa.net
Jean-Michel Ferrard
Algorithmique avec Maple
Décompositions LU et Choleski
Quelques démonstrations
Quelques démonstrations
1. Existence et unicité de la décomposition LU
Proposition [Retour Énoncé]
Une matrice A de GLn(K) possède une décomposition LU si et seulement si ses mineurs princi-
paux sont non nuls. Cette décomposition est alors unique.
Démonstration
– Unicité de la décomposition LU
Supposons que A ait les décompositions A = LU et A = L0U 0. Alors L0−1L = U U 0−1.
Mais les matrices triangulaires supérieures inversibles forment un sous-groupe de GLn(K), de même
que les matrices triangulaires inférieures à diagonale unité. On en déduit que la matrice L0−1L =
U U 0−1 est a diagonale unité et qu’elle est a la fois triangulaire inférieure et triangulaire supérieure
donc diagonale : ce ne peut être que la matrice In.
Il en découle L = L0 et U = U 0. On ainsi prouvé l’unicité (si existence.)
– Existence de la décomposition LU (sens direct)
On suppose que la matrice inversible A possède une décomposition LU .
On veut montrer que toutes ses sous-matrices principales sont inversibles.
Pour tout k de {1, . . . , n}, on utilise une décomposition par blocs de A, L, U :
A =
(cid:19)
(cid:18) Ak A0
k
A00
k A000
k
, L =
(cid:19)
(c...