Livre_silo 30 août 2013 16:32 Page 182
¨
©
¨
©
¨
©
¨
©
C o p y r i g h t E y r o l l e s
182
Informatique pour tous
POUR ALLER PLUS LOIN Presque toutes les matrices ont une décomposition LU .
Tout d’abord, si on prend des coefficients aléatoires (en différents sens raisonnables), on
obtient avec probabilité 1 un système de Cramer (une matrice inversible).
Encore mieux : on n’aura en général pas de problème de pivot ! Plus précisément, si les
mineurs principaux (les matrices (k, k) « en haut à gauche » extraites de la matrice initiale)
sont tous inversibles, alors à chaque étape k du pivot ⁵, le coefficient présent en position
(k, k) est non nul, et on peut l’utiliser pour pivoter.
Ainsi, partant d’une matrice A, on trouve une matrice L 1 triangulaire inférieure, telle que
U = L 1 A soit triangulaire supérieure. Si on note L = L
−1
1 , on a la décomposition A = LU
(L pour Lower et U pour Upper).
Si par malheur on trouve à l’étape k le coefficient a k,k qui est nul, alors on prend le plus
petit j > k tel que a j,k soit non nul et on l’utilise comme pivot pour placer des zéros dessous et à droite par des transvections sur les colonnes. On parvient ainsi, modulo quelques
dernières dilatations, à construire deux matrices L et U respectivement triangulaires inférieure et supérieure, telles que A = LTσU , avec Tσ une matrice de permutation : c’est la
décomposition de Bruhat ⁶.
Les matrices L et U ne sont pas uniques, mais la permutation σ l’est. Dans la partition de
GLn(R) en n! composantes, celles associées à l’identité — c’est-à-dire celles possédant une
décomposition LU — constituent la « grosse cellule » : c’est un ouvert dense.
7.2 Mise en œuvre
7.2.1 Découper le travail
On commence par réfléchir aux bons outils, à savoir les fonctions et programmes auxiliaires à l’aide desquels l’écriture du programme principal sera quasiment une traduction
en anglais de l’algorithme ! Il faudra déléguer les opérations suivantes :
• la recherche d’un pivot ;
• les échanges de lignes ;
• les transvections.
Le premier programme prend en entrée une matrice A et un indice i. Il doit renvoyer un
indice j ⩾ i tel que |a j,i | soit maximale :
def chercher_pivot(A, i):
n = len(A) # le nombre de lignes
j = i # la ligne du maximum provisoire
for k in range(i+1, n):
if abs(A[k][i]) > abs(A[j][i]):
j = k # un nouveau maximum provisoire
return j # en faisant bien attention à l'indentation :-)
5. Le pivot normal et non partiel.
6. Les analystes numériciens ont tendance à privilégier l’ordre LU P (la matrice de permutation en dernier).
¨
©
¨
©
¨
©
¨
©
C o p y r i g h t E y r o l l e s
182
Informatique pour tous
POUR ALLER PLUS LOIN Presque toutes les matrices ont une décomposition LU .
Tout d’abord, si on prend des coefficients aléatoires (en différents sens raisonnables), on
obtient avec probabilité 1 un système de Cramer (une matrice inversible).
Encore mieux : on n’aura en général pas de problème de pivot ! Plus précisément, si les
mineurs principaux (les matrices (k, k) « en haut à gauche » extraites de la matrice initiale)
sont tous inversibles, alors à chaque étape k du pivot ⁵, le coefficient présent en position
(k, k) est non nul, et on peut l’utiliser pour pivoter.
Ainsi, partant d’une matrice A, on trouve une matrice L 1 triangulaire inférieure, telle que
U = L 1 A soit triangulaire supérieure. Si on note L = L
−1
1 , on a la décomposition A = LU
(L pour Lower et U pour Upper).
Si par malheur on trouve à l’étape k le coefficient a k,k qui est nul, alors on prend le plus
petit j > k tel que a j,k soit non nul et on l’utilise comme pivot pour placer des zéros dessous et à droite par des transvections sur les colonnes. On parvient ainsi, modulo quelques
dernières dilatations, à construire deux matrices L et U respectivement triangulaires inférieure et supérieure, telles que A = LTσU , avec Tσ une matrice de permutation : c’est la
décomposition de Bruhat ⁶.
Les matrices L et U ne sont pas uniques, mais la permutation σ l’est. Dans la partition de
GLn(R) en n! composantes, celles associées à l’identité — c’est-à-dire celles possédant une
décomposition LU — constituent la « grosse cellule » : c’est un ouvert dense.
7.2 Mise en œuvre
7.2.1 Découper le travail
On commence par réfléchir aux bons outils, à savoir les fonctions et programmes auxiliaires à l’aide desquels l’écriture du programme principal sera quasiment une traduction
en anglais de l’algorithme ! Il faudra déléguer les opérations suivantes :
• la recherche d’un pivot ;
• les échanges de lignes ;
• les transvections.
Le premier programme prend en entrée une matrice A et un indice i. Il doit renvoyer un
indice j ⩾ i tel que |a j,i | soit maximale :
def chercher_pivot(A, i):
n = len(A) # le nombre de lignes
j = i # la ligne du maximum provisoire
for k in range(i+1, n):
if abs(A[k][i]) > abs(A[j][i]):
j = k # un nouveau maximum provisoire
return j # en faisant bien attention à l'indentation :-)
5. Le pivot normal et non partiel.
6. Les analystes numériciens ont tendance à privilégier l’ordre LU P (la matrice de permutation en dernier).