3.4 Autres types de factorisation
85
avec α ∈ R
+ , v
T ∈ R
i−1 et cherchons une factorisation de Ai de la forme
Ai = H
T
i Hi =
H
T
i−1
0
h
T
β
Hi−1 h
0
T
β
.
Par identification avec les ´ el´ ements de Ai, on obtient les ´ equations H
T
i−1 h = v et
h
T h + β
2 = α. Le vecteur h est ainsi d´ etermin´ e de fa¸ con unique, puisque H
T
i−1 est
inversible. De plus, en utilisant les propri´ et´ es des d´ eterminants, on a
0 < d´ et(Ai) = d´ et(H
T
i ) d´ et(Hi) = β
2 (d´ et(Hi−1))
2 ,
ce qui implique que β est un nombre r´ eel. Par cons´ equent, β =
√
α − h T h est
l’´ el´ ement diagonal cherch´ e, ce qui conclut la preuve par r´ ecurrence.
Montrons maintenant les formules (3.42). Le fait que h11 =
√
a11 est une cons´ equence imm´ ediate de la r´ ecurrence au rang i = 1. Pour un i quelconque, les relations (3.42)1 sont les formules de “remont´ ee” pour la r´ esolution du syst` eme lin´ eaire
H
T
i−1 h = v = [a1i, a2i, . . . , ai−1,i]
T , et les formules (3.42)2 donnent β =
√
α − h T h,
o` u α = aii.
3
L’algorithme correspondant `
a (3.42) n´ ecessite environ (n
3 /3) flops, et il est
stable par rapport `
a la propagation des erreurs d’arrondi. On peut en effet
montrer que la matrice triangulaire sup´ erieure ˜
H est telle que ˜
H
T ˜
H = A + δA,
o` u δA est une matrice de perturbation telle que δA 2 ≤ 8n(n + 1)uA 2 ,
quand on consid` ere les erreurs d’arrondi et qu’on suppose 2n(n + 1)u ≤ 1 −
(n + 1)u (voir [Wil68]).
Dans la d´ ecomposition de Cholesky, il est aussi possible de stocker la matrice H
T dans la partie triangulaire inf´ erieure de A, sans allocation suppl´ ementaire de m´ emoire. En proc´ edant ainsi, on conserve `
a la fois A et la partie
factoris´ ee. On peut en effet stocker la matrice A dans le bloc triangulaire sup´ erieur puisque A est sym´ etrique et que ses termes diagonaux sont donn´ es par
a 11 = h
2
11 , a ii = h
2
ii +
i−1
k=1 h
2
ik , i = 2, . . . , n.
Un exemple d’impl´ ementation de la d´ ecomposition de Cholesky est propos´ e
dans le Programme 7.
Programme 7 - chol2 : Factorisation de Cholesky
function [A]=chol2(A)
% CHOL2 Factorisation de Cholesky d’une matrice A sym. def. pos.
% R=CHOL2(A) renvoie une matrice triangulaire sup´ erieure R telle
% que R’*R=A.
[n,m]=size(A);
if n ˜= m, error(’Seulement les syst` emes carr´ es’); end
for k=1:n-1
if A(k,k) <= 0, error(’Pivot nul ou n´ egatif’); end
A(k,k)=sqrt(A(k,k)); A(k+1:n,k)=A(k+1:n,k)/A(k,k);
for j=k+1:n, A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k); end
end
85
avec α ∈ R
+ , v
T ∈ R
i−1 et cherchons une factorisation de Ai de la forme
Ai = H
T
i Hi =
H
T
i−1
0
h
T
β
Hi−1 h
0
T
β
.
Par identification avec les ´ el´ ements de Ai, on obtient les ´ equations H
T
i−1 h = v et
h
T h + β
2 = α. Le vecteur h est ainsi d´ etermin´ e de fa¸ con unique, puisque H
T
i−1 est
inversible. De plus, en utilisant les propri´ et´ es des d´ eterminants, on a
0 < d´ et(Ai) = d´ et(H
T
i ) d´ et(Hi) = β
2 (d´ et(Hi−1))
2 ,
ce qui implique que β est un nombre r´ eel. Par cons´ equent, β =
√
α − h T h est
l’´ el´ ement diagonal cherch´ e, ce qui conclut la preuve par r´ ecurrence.
Montrons maintenant les formules (3.42). Le fait que h11 =
√
a11 est une cons´ equence imm´ ediate de la r´ ecurrence au rang i = 1. Pour un i quelconque, les relations (3.42)1 sont les formules de “remont´ ee” pour la r´ esolution du syst` eme lin´ eaire
H
T
i−1 h = v = [a1i, a2i, . . . , ai−1,i]
T , et les formules (3.42)2 donnent β =
√
α − h T h,
o` u α = aii.
3
L’algorithme correspondant `
a (3.42) n´ ecessite environ (n
3 /3) flops, et il est
stable par rapport `
a la propagation des erreurs d’arrondi. On peut en effet
montrer que la matrice triangulaire sup´ erieure ˜
H est telle que ˜
H
T ˜
H = A + δA,
o` u δA est une matrice de perturbation telle que δA 2 ≤ 8n(n + 1)uA 2 ,
quand on consid` ere les erreurs d’arrondi et qu’on suppose 2n(n + 1)u ≤ 1 −
(n + 1)u (voir [Wil68]).
Dans la d´ ecomposition de Cholesky, il est aussi possible de stocker la matrice H
T dans la partie triangulaire inf´ erieure de A, sans allocation suppl´ ementaire de m´ emoire. En proc´ edant ainsi, on conserve `
a la fois A et la partie
factoris´ ee. On peut en effet stocker la matrice A dans le bloc triangulaire sup´ erieur puisque A est sym´ etrique et que ses termes diagonaux sont donn´ es par
a 11 = h
2
11 , a ii = h
2
ii +
i−1
k=1 h
2
ik , i = 2, . . . , n.
Un exemple d’impl´ ementation de la d´ ecomposition de Cholesky est propos´ e
dans le Programme 7.
Programme 7 - chol2 : Factorisation de Cholesky
function [A]=chol2(A)
% CHOL2 Factorisation de Cholesky d’une matrice A sym. def. pos.
% R=CHOL2(A) renvoie une matrice triangulaire sup´ erieure R telle
% que R’*R=A.
[n,m]=size(A);
if n ˜= m, error(’Seulement les syst` emes carr´ es’); end
for k=1:n-1
if A(k,k) <= 0, error(’Pivot nul ou n´ egatif’); end
A(k,k)=sqrt(A(k,k)); A(k+1:n,k)=A(k+1:n,k)/A(k,k);
for j=k+1:n, A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k); end
end