Télécharge Notes sur le thème de l'algorithmique numérique - 1° partie et plus Notes au format PDF de Application informatique sur Docsity uniquement!
Généralités
Normes et rayon spectral:
Définition norme vectorielle :
x x R R
n , vérifiant :
- x 0 et x 0 x 0
- . x .x
- x y x y
Exemple de norme vectorielle :
X 1 xi
i ,
2 2 ^
X (^) Supx
Proposition : En dimension finie, toutes ces normes sont équivalentes ; en pratique, on choisit celle qui nous arrange.
Définition norme matricielle :
C'est une norme dans l'espace vectorielle Mn(R), avec en plus A B A.B.
Norme matricielle induite : A partir d'une norme vectorielle, on construit une norme matricielle induite :
Sup AX X
AX
A Sup X 0 X 1 ou X 1
Notons de plus que I 1.
En particulier pour la norme 2, on a la relation : I AA
t
2
avec le rayon spectral.
Propriétés :
- AX A.X et AB A.B
- Cas de A symétrique : (^) i i
A A max
2
- Dans le cas général, (^) max 2
A , maximum des valeurs singulières i iavec i les valeurs
propres (positives) de la matrice symétrique A A
t
. (?)
- Soit A une matrice carré quelconque n x n. Pour toutes normes de matrice induite, on a A A.
Conditionnement d’une matrice inversible :
Considérons un système linéaire AX=b. L'expérience de Wilson met un évidence une instabilité des calculs dans certains cas. Si l'on modifie sensiblement les paramètres de la matrice A, ou du second membre b, on peut trouver des résultats complètement différends!
1 .
Cond A A A
- Cond A Cond A
- Cond A 1 avec Cond Id 1
- Soit une matrice Q orthogonale : 2 2
QX X et Cond Q 1
- Les bons conditionnements sont les petits conditionnements (au mieux 1 pour les matrices orthogonales).
- Pour A symétrique :
min
max 2
Cond A ; dans la cas général :
min
max 2
Cond A
Théorèmes sur le conditionnement :
Théorème : AX Bet A XX bb. On a
b
b Cond A X
X
.
Théorème : AX Bet AA XX b. On a
A
A
Cond A X X
X
.
Remarque : Le conditionnement traduit l'amplification des grandeurs relatives.
Inverse numérique d’une matrice :
Soit A une matrice inversible et A
M est un inverse de A du point de vue numérique si :
1
A
M A
est petit, ce qui correspond à M=A-
M A
1
est petit, ce qui correspond à A=M-
- AM^ ^ Id et^ MAId sont petits, ce qui correspond à AM Id 0 etMA-Id 0
On s'adapte au calcul que l'on veut faire. Notons que l'on résout rarement un système linéaire en calculant
l'inverse X A b
1
On veut traiter tous les seconds membres en une seule passe ; pour cela on considère la matrice n x 2n suivante :
qui après calculs donne
.
On donne ici un exemple d’écriture en langage de description algorithmique :
Procédure pivot ( k , l : entiers ; OK : booléens ) Début
Max := A k , k;
l := k ; Pour i := k+1 à n faire
Si max < A i ,k
Alors Début
Max := A i , k;
l := i ; Fin ; OK := Max > ; Fin ;
Procédure permuter ( k , l : entiers ) Début Si k l Alors Pour j := k à 2n faire Début
c := A^ k^ ,^ j;
A^ k^ , j: = A l, j;
A^ l^ , j := c ;
Fin ; Fin ;
On donne maintenant le programme principal, réalisant le calcul de l’inverse.
Début Initialiser A et lire ;
/* triangularisation */ Pour k := 1 à n faire Début Pivot (k , l , OK ) ; Si non(OK) Alors Début Ecrire « matrice non inversible » ; Exit echec ; Fin ; Permuter ( k , l ) ; Pour j := k+1 à 2n faire
A k, j : A k,j /A k,k;
Pour i := k+1 à n faire Pour j := k+1 à 2n faire
A i , j : Ai ,k A i,k *A k,j;
Fin ;
/* n remontées */ Pour k := 1 à n faire Pour i := n-1 à 1 par pas de -1 faire Début s := 0 ; Pour j := i+1 à n faire
s : sA i ,j *A j,nk;
A i , n k :Ai ,nk s;
Fin ; Fin ;
Remarques :
- Choix du pivot : pour des raisons de stabilité numérique, on divise toujours par le terme de la colonne k qui
a la plus grande valeur absolue (Cf. Pivot partiel maximum).
- Recherche du pivot : si la plus grande valeur absolue est inférieure à la précision machine , alors on arrête
en disant que la matrice n’est pas inversible à près (test d’inversibilité informatique). Cependant, elle peut très bien être inversible du point de vue mathématique!
Factorisation A=LU
Soit A une matrice régulière, de dimension n.
Lien avec la méthode de Gauss
On applique la méthode de Gauss sans pivoter. U est une matrice triangulaire supérieure réelle, telle que
J A U k n
k 1 ,
.
On donne la définition des matrices
,
1 ,
n k
Jk k k
avec k k
lk l k a
a
,
,
, pour l k.
On démontre que l’inverse des J (^) k existe et que (^)
k n
L Jk 1 ,
1 est une matrice triangulaire inférieure, avec des 1
sur la diagonale. En fait, on a
,
1 ,
1
n k
Jk k k
.
Calcul algorithmique des coefficients de L et de U
On procède par identification en utilisant les propriétés de L (triangulaire inférieure avec une diagonale de 1) et
de U (triangulaire supérieure). On cherche un algorithme par ligne donnant les
li (^) , jpour 1 ji 1 etlesui,jpour 1 j n.
Pour i := 1 à n faire Début Pour j := 1 à i-1 faire
Matrice de permutation élémentaire i j :
lignej
lignei
1
1
1 0
1
1
0 1
1
1
,
Pi j
Le produit PA permute les lignes i et j de la matrice A. Le produit AP permute les colonnes i et j de la matrice A.
(à compléter avec le cours : obtention de la factorisation et utilisation de la factorisation )
Factorisation A=BB
t : Cholesky
Cette factorisation s'applique pour des matrices A symétriques et positives.
Théorème : Soit A une matrice SDP de dimension n. Il existe une et une seule matrice B triangulaire inférieure,
telle que A=BB
t .
Factorisation A=LDL
t : Crout
Théorème : Soit A une matrice SDP de dimension n. Il existe une et une seule matrice L triangulaire inférieure et
D diagonale positive, telles que A=LDLt.
Remarque : Crout a le même comportement numérique que Cholesky. Cependant, il est plus utilisé car
généralisable au cas complexe.
Soit A une matrice symétrique définie positive (quelconque ?).
L est une matrice triangulaire inférieure de dimension n dont la diagonale se compose de 1.
D est une matrice diagonale de dimension n.
L est en fait la matrice B pour laquelle les termes d'une colonne sont divisés par les coefficients b (^) ii. D est la
matrice dont les coefficients diagonaux sont les b (^) ii
² .
Algorithme par ligne pour l’obtention de la factorisation
On écrit
li, j
L et
d n
d
D
1
.
Il s'agit de calculer les li,j et les di. On a (^) ij j
j
k
a (^) i j (^) likdkljkl d
1
0
,.
On établit un algorithme par ligne : Pour la 1ère^ ligne, on écrit d 1 a 1 , 1. Supposons maintenant le calcul
effectué jusqu'à la ligne i-1 ; le calcul de la ligne i est donné par
j
j
k
ik k jk
i j ij d
l dl
l a
1
1 , , pour^ j^ ^ i^1 et
par (^)
1
1
2 ,
j
k
d (^) i aii lik dk. Ainsi une CNS pour cette factorisation est que les d (^) i ne soient pas nuls.
On donne l’algorithme correspondant :
Pour i de 1 à n faire Début Pour j de 1 à i-1 faire
(^) j
j
k
li (^) j: aij likdkljk /d
1
1
, ,
/* li (^) , i: 1 */
^ ^
1
1
2 : ,
j
k
d (^) i aii lik dk ;
Fin ;
Ce qui correspond en effet à une factorisation de Crout, car la matrice D diest bien diagonale et la matrice
L est bien triangulaire inférieure, par construction.
On donne l’algorithme de calcul de L et de D par colonne, qui procède en écrasant la partie triangulaire
inférieure de A, (les 1 de la diagonale de L ne sont pas stockés, car on préfère y stocker D) :
Début Pour k de 1 à n-1 faire Pour j de k+1 à n faire Début
d : aj ,k ak,k ;
Pour i de j à n faire a (^) i , j: ai,jai,k.d ;
a (^) j , k:d ;
Fin ; Fin ;
La structure donnée d’implémentation serait un vecteur à une dimension stockant consécutivement les colonnes
de la partie triangulaire inférieure de A.
Remarque fondamentale : ai (^) , jest modifié par a (^) i,k et a (^) j,k.
Etude du remplissage lors de la factorisation de matrices SDP creuses
Considérons une matrice A symétrique creuse. On utilise une structure de donnée type profil, il est important de
savoir où se trouvent les termes à priori non nuls de L.
Une matrice est dite creuse si elle contient « un grand nombre » de termes nuls. Plus précisément, on considère
que la matrice est creuse à partir du moment, où son pourcentage de termes nuls permet un gain informatique par
rapport à un traitement classique qui la considérerait pleine.
Définitions :
- Un terme de remplissage apparaît en (i,j) au cours de la factorisation si et seulement si ai (^) , j 0 et
li (^) , j 0. On reformule le problème comme la conservation du creux initial de la matrice.
- li,j est logiquement non nul si ai (^) , j 0 ouk, 1 kj 1 telqueli,klj,k 0. C’est à dire, à la k éme
itération 1 k n 1 un terme ai , jde la matrice A qui est nul, devient non nul ssi
a (^) j , k 0 etai,k 0 avec k j 1. Preuve (algorithme par colonne): kk
ik jk i j ij a
a a a a ,
0
, ,
0
, 0
,
- Les li,j logiquement non nuls se trouvent à l'intérieur du profil de A.
Exemple de remplissage total :
x x x x x
x x x x
x x x
x x
x
E
x x
x x
x x
x x
x
E 1 2
Le coefficient (3,2) dépend uniquement des coefficients (3,1) et (2,1) qui sont non nuls, par conséquent le
coefficient (3,2) devient logiquement non nul : effet de remplissage. Ensuite, le coefficient (4,2) qui dépend des
coefficients (4,1) et (2,1) devient logiquement non nul. (…) Le coefficient (5,3) dépend des coefficients (5,k) et
(3,k) pour k variant de 1 à 3-1=2, etc.
Remarque : Ainsi, on voit que le phénomène dépend de l’ordre selon lequel on exécute les itérations, c’est à dire
de la numérotation des inconnues du système! Si l’on permute les indices 5 et 1 dans l’exemple précédent, on
obtient la matrice suivante :
x x x x x
x
x
x
x
E
3 qui ne produit pas de remplissage!
La question qui se pose maintenant est de savoir si l’on peut trouver une matrice de permutation P, telle que la
factorisation de la matrice PAPt^ crée le moins de remplissage possible? Il n’existe pas de méthode exacte en
temps raisonnable et on va utiliser des méthodes approchées (heuristiques).
Caractérisation du remplissage par les graphes
On utilise un graphe afin de déterminer quels sont les termes logiquement non nuls qui vont apparaître au cours
du calcul.
Le graphe associé à une matrice A symétrique d'ordre n possède n sommets que sont les indices de 1 à n des
coefficients de la matrice. Un coefficient non nul est repéré dans le graphe par un arc entre les indices de ligne et
de colonne de ce coefficient.
Exemple : Graphe associé à la matrice E 1 et à la permutation identité :
Considérons une permutation donnée telle que x i iet plaçons nous à la k
éme itération :
- trouver le monotone adjacent de^ x^ kqui est l’ensemble de tous les sommets liés à^ x^ kavec un numéro plus
grand, Madj xk x (^) j, jket telque xj,xksoitunearêtepossible
- rajouter une arête de remplissage pour tout couple de tels sommets non liés entre eux.
- On fait cela dans l’ordre des k croissants de 1 à n-1. A la fin, on obtient le graphe dit d’élimination.
Exemple :
- Madj x 1 x 2 , x 3 ,x 4 ,x 5 remplissage total
- Madj x 2 pas de remplissage
- Madj x 3 pas de remplissage
- Madj x 4 pas de remplissage
Profil ligne
Considérons une matrice SDP A de dimension n.
Définition : Pour chaque ligne, on considère la colonne c (^) icorrespondant au premier terme non nul sur la ligne.
La contribution de la ligne i au profil est Pi a (^) i, j,cijiqui sera effectivement stockée par ordre
croissant de colonne (soit i ci 1 coefficients). Le profil ligne est l’union des Pi pour i variant de 1 à n. On a
besoin d’un vecteur de dimension n donnant pour tout i la valeur de ci.
La taille du profil est (^)
n
i
m i ci 1
Propriété fondamentale : Le profil ligne se conserve par factorisation (statique). Le remplissage est inclus dans
le profil ligne. Par conséquent, on va chercher à trouver des numérotations qui tendent à minimiser la taille du
profil ligne.
