Télécharge Notes sur le thème de l'algorithmique numérique - 2° partie et plus Notes au format PDF de Application informatique sur Docsity uniquement!
On adapte l’algorithme de factorisation :
j
j
k c c
li j aij likdkljk d i j
1
max ,
, ainsi on diminue le volume de
calcul.
Optimisation du profil ligne
Critère approché : On va chercher une numérotation des sommets du grapheG minimisant le profil A. Pour
cela, on cherche à diminuer la taille m. On utilise un critère local : à i fixé, on va essayer d’obtenir c (^) iproche
autant que possible de i, ce qui va induire un tassement vers la diagonale.
Le problème en terme de numérotation surG est le suivant : « numéroter ses voisins avec des numéros les plus
proches possibles de i ».
Numérotation en couches :
On partitionne l’ensemble des sommets du graphe par rapport à leur distance à un sommet de départ x. On numérote continûment selon les couches et le plus grand écart varie comme un O(longueur d’une couche). Le
nombre de sommets étant fixé, on cherche le sommet x de départ de telle manière que le nombre de couche
h x générés à partir de x soit maximum ; ceci aura tendance à minimiser la longueur des couches.
Exemple : prenons x = 1, on a h(x) = 1.
Algorithme :
- Choisir x quelconque
- Construire la structure en couche à partir de x
- Choisir un sommet y de la couche terminale ayant le moins de voisins et construire les couches à partir de ce
nouveau sommet, nécessairement h y h x
- Si h(x) > h(y) Alors réitérer, Sinon, on est dans le cas où h(x) = h(y), et y est le point de départ.
A partir de l’exemple précédent, on choisit y = 2 et on obtient h(y) = 2. En réitérant de nouveau, on obtiendra
pas mieux. Donc on déduit x 2.
On applique le changement de numérotation :
Ancien numéro i 1 2 3 4 5
Nouveau numéro i^2 1 3 4
Ce qui donne la transformation suivante :
x x
x x
x x
x x
x
E
x x
x x
x x
x x
x
E
Limitation du remplissage
Théorème : CNS de création d’une arête de remplissage
SoitG le graphe numéroté associé à la matrice A.
Une arête de remplissage xi , xjest créée si et seulement si :
- l’arête est déjà présente dansG ;
- il existe un chemin allant de xi à x (^) jtel que les numéros intermédiaires soient plus petits que ceux des
extrémités.
Idée de numérotation : Pour mettre en défaut ce théorème, on utilise des séparateurs topologiques. Cette méthode consiste à séparer deux sous-ensembles A et B par un ensemble C, appelé séparateur ; de sorte que tout
chemin allant d’un élément de A vers un élément de B passe au moins par un sommet de C. Si on numérote les
sommets de C avec des numéros plus grands que ceux affectés à A et B, alors il n’y aura pas d’arête de remplissage créées entre les éléments de A et ceux de B.
On applique ce principe de numérotation récursivement par bloc, jusqu’à tomber sur une grille qui n’est plus séparable. Ce faisant on limite excellemment le remplissage, mais le profil risque d’être mauvais.
Remarque : Il faut une autre structure de donnée que celle du profil non adapté. Par ailleurs, il y a aussi une indépendance des calculs exploitable sur une machine à architecture parallèle.
Exemple : On reprend l’exemple précédent…
On effectue le changement de numérotation :
x x x x x
x
x
x
x
E
x x
x x
x x
x x
x
E
1 0 5 et
l’on peut vérifier que cette dernière matrice ne se remplie pas du tout!
Utilisation de la factorisation pour résoudre Ax=b
Il s’agit de résoudre Ax ben utilisant la factorisation de Cholesky - CROUT : ALDLX b t . On procède
en trois étapes.
t et résolvons Ly b, dont la complexité est en 2
2 n (matrice triangulaire).
On écrit un algorithme utilisant un accès par ligne dans la matrice L.
y 1 : b (^1) ;
Pour i de 2 à n faire
1
1
,
i
j
yi bi lijyj ;
- Résolution du système diagonal : Posons z Lx t et résolvons Dz y(par n divisions).
A B
C
B
A
C
Systèmes linéaires Ax=b : méthodes itératives
Principe
Méthode itératives classiques :
On cherche à obtenir une suite de vecteur X nnconvergeante vers X tellequeAX b.
En posant A MN où M est une matrice inversible, on propose la formule de récurrence
MX (^) n 1 NXn bou encore X (^) n M NXn M b 1 1 1
^ . Si^ X^ est solution du problème, il vient
que ^ ^ ^ ^ ^0
1 X X M N X X
n (^) n
. B ^ M N
1 est appelée la matrice d'itération.
Ainsi pour que X n X, il est nécessaire que ^ ^0
1
n M N.
Remarque : Le cadre normale de convergence des méthodes itératives classiques est : « A matrice symétrique définie positive ».
Vitesse de convergence :
Théorème : B n 0 ssi B 1
n
Vitesse de convergence : Considérons la matrice d'itération B symétrique, alors B B
2
et l'on obtient
le résultat suivant :
2 0 2 x x B. x x
k
k ^ .
En conclusion, le rayon spectral de la matrice d'itération mesure la vitesse de convergence. La convergence est d'autant plus rapide que le rayon spectral est plus petit. On admettra que ce résultat ce généralise aux matrices quelconques.
Méthode Jacobi
Notation :
E
D
F
A
On pose A MNavec MDetNEFselon le schéma ci-dessus. Ce qui donne la formule de
récurrence suivante : X n D E F Xn D b
1 1 1
^ .
Soit i n
i x (^) k xk 1 , le vecteur itéré, on donne l’algorithme donnant x (^) k 1 en fonction de x (^) k.
Pour i de 1 à n faire
i
j k
n
ji
ij
j k
i
j
ij ii
i k a x a x b a
x 1
,
1
1
, ,
1
Théorème : Soit A une matrice symétrique, à diagonale strictement dominante, c'est à dire
i j
ai (^) ,i ai,j
alors la méthode est convergente car 1
1
D E F.
A priori, cette méthode converge pour une matrice A SDP.
Méthode Gauss-Seidel
Boucle directe
On reprend le schéma de matrice précédent et on pose M DE et N F. La formule de récurrence
s'écrit : D E Xn 1 FXnb.
On calcule i xk (^) 1 , l'i ème (^) composante du vecteur X (^) k 1 : (^) i
n
ji
j ij k
j k
i
j
ij
i
a i ixk a x a x b
1
1 ,
1
1
, 1 ,.
On propose l’algorithme suivant :
Pour i de 1 à n faire
i
n
ji
j ij k
j k
i
j
ij ii
i k a x a x b a
x 1
1 ,
1
1
, ,
1
On procède par écrasement du vecteur X (^) k:
t
j k
i k
j xk j i x x i j n
,^1 ^ ^1 , ^1
calculencours
1 déjàcalculés
Boucle rétrograde
Cette fois, on pose M^ ^ DFet N E.
On propose l’algorithme suivant :
Pour i de n à 1 faire
i
n
ji
j ij k
j k
i
j
ij ii
i k a x a x b a
x 1
, 1
1
1
, ,
1
Conclusion
Théorème : Soit A une matrice symétrique définie positive, alors la méthode est convergente
car ^ ^1
1
D E F.
Comparaison : Si A est SDP les méthodes de Jacobi et de Gauss-Seidel convergent, et Gauss-Seidel converge
plus vite que Jacobi car J L J D E F L D E F D F E
1 1 1
(^2 ) avec et ou
Méthode de relaxation
Pour cette méthode, on pose E w
D
M et D F w
w N M A
, avec w 0. On donne la
matrice d’itération : (^)
D F w
w E w
D
L (^) w
1
. Ainsi pour w 1 , on se ramène à la méthode de
Gauss-Seidel.
Théorème : Soit A une matrice SDP. Il y a convergence si et seulement si Lw 1. On démontre que ceci est
vrai si et seulement si w 0 , 2 .
Remarque :
- w 1 : Gauss-Seidel
- w 1 : Sous-relaxation
- w 1 : Sur-relaxation
Méthodes du gradient
1
Principe : méthode de la plus grande descente
Soit A une matrice symétrique définie positive.
Théorème : Soit J X AX, X 2 b,X. Alors X est solution de AX b, si et seulement si X réalise
le minimum de J.
En dimension 1, A a J x ax bx
2 et (équation de parabole) ; d'où le minimum de J est a
b , ce qui
correspond bien à la solution de l'équation linéaire ax b.
Formule de récurrence : Prenons un X 0 quelconque, et supposons que nous ayons obtenu par itération X (^) k.
Considérons la courbe de niveau k X tel queJ X J Xk, pour calculer X k+1, on se déplace du point X k,
appartenant là a courbe de niveau k, dans la direction de son gradient : g (^) k gradAXkb 2
. On donne
ainsi la formule de récurrence : Xk 1 Xk kgk. Il reste maintenant à choisir la valeur de k, de sorte que
J soit le plus petit possible. On montre que J Xk k Agk,gk 2 k gk,gk J Xk
2
1 ^ , donc il faut
prendre
k k
k k k Ag g
g g
,
On résume le passage de Xk à X (^) k+1 :
k k k Ag g
g g
,
Algorithme du gradient conjugué
Dans cette méthode, on va utiliser un ensemble de directions conjugués. g (^) k est la meilleure direction locale,
mais pas à long terme, sinon l'algorithme convergerait d'un coup! D'où l'idée, d'apporter une correction pour le
choix de la direction.
On pose Xk 1 Xk kdkavec dk gk k 1 dk 1. En résolvant le système obtenu, par annulation des
dérivées partielles pour la fonction J (qui doit être minimum), on calcule les valeurs des coefficients ket k 1.
On va maintenant résumer la méthode :
- Initialisation : X 0 quelconque, g 0 AX 0 b, et d 0 g 0.
- Passage de l'indice k à l'indice k+1 : On suppose que l'on connaît X (^) k , dk,etgk. On calcule
k k
k k k Ad d
g g
,
, d'où Xk 1 Xk kdk.
1 Cette partie a été rédigée en privilégiant les notes de cours plutôt que celle de TD, où les notations et les méthodes diffèrent sur quelques points.
- On prépare la suite en calculant : g (^) k 1 gk (^) kAdkou AXk 1 b; k k
k k k g g
g g
,
1 ,^ 1
et
dk 1 gk 1 kd k.
Remarque : On démontre que la correction apportée dans cette méthode est optimale.
Formule par récurrence : On démontre après coup que cette méthode donne la solution de l'équation linéaire au bout de N itérations! Par conséquent, on pourrait la considérer comme une méthode directe, mais on préfère la
prendre comme une méthode itérative, et s'arrêter pour k Nquand la précision est satisfaisante.
Préconditionnement :
Convergence de la formule :
On établit que :
A
k
k (^) A
X X
Cond A
Cond A X X
. Notons que plus Cond(A) est petit, plus la
convergence est rapide. En particulier, si Cond A 1 , alors on converge d'un coup, c'est à dire X 0 est déjà
solution! Rappelons : Cond A 1. Ici, on voit que la vitesse de convergence va dépendre de la matrice A, et
plus particulièrement de son conditionnement. On va alors introduire la notion de préconditionnement.
Transformation :
On cherche à résoudre un autre système AY b
ayant même solution, mais pour lequel la matrice A
possède
un meilleur conditionnement.
On considère t C E.E. On pose t A E AE ..
et b E b
. Ainsi AX b AY b
pour
Y E X t .
Remaque : Les algorithmes passant par l'intermédiaire de Y sont dits non transformés. Lorsque Y n'apparaît plus dans le calcul l'algorithme est dit transformé.
Algorithme du gradient préconditionné transformé :
Prenons un X 0 quelconque, et supposons que nous ayons obtenu par itération Xk. On pose g (^) k AXkbet
hk C g k 1 , ce qui revient à résoudre Chk gk. Résolution d'un système linéaire! Certes. Mais simple,
puisqu'on a supposé auparavant, que l'on avait une factorisation de Cholesky de la matrice C, c'est à dire
t C E.E (algorithme de complexité n²). Ensuite, on calcule k k
k k k Ah h
g h
,
. Et on passe de Xk à Xk+1 en
faisant Xk Xk khk
Ainsi, on obtient le résultat fondamental suivant
A
k
k (^) A
X X
CondA
CondA X X
. Il faut à présent
travailler pour que le conditionnement de à soit meilleur que celui de A. Or, nous avons posé
t A E AE
par où nous pouvons déduire que Cond A CondC A
, car ces deux matrices ont les mêmes valeurs
propres.
Choix de C (méthode heuristique) :
- Factorisation incomplète de Cholesky : Lorsque l'on pose t A BB, il peut sur venir un phénomène de
remplissage! D'où l'idée de réaliser une factorisation incomplète A EE R
t . On ne calcule dans la matrice E que les l (^) i , jtels queai,j 0 , pour lesquels donc on évite le phénomène de remplissage.
Si l'on pose C BB A t , on obtient un conditionnement de 1 (le meilleur) car C AId 1 ; cependant il peut se produire un phénomène de remplissage peut souhaitable. Et si l'on pose C Id, le conditionnement n'est pas amélioré par rapport à celui de A.
Problèmes au moindres carrés (linéaire)
Régression linéaire
Etant donné une série de M points xi , yi, on cherche une relation de la forme y i axib. On cherche a et
b qui minimise l'erreur
M
erreur ei
(^22)
avec e i yiax ib.
Cas général
On étudie un cas plus général
N
f ti xj. jti , avec N inconnues x j. On a M équations avec M N.
On pose
j i M N A t
et b ft iM. L'erreur est Ax b. On cherche les x jqui minimise cette
erreur
i M j N
Ax b aij xj bi 1
2
1
2 2 .
Méthode des équations normales
Théorème
Le vecteur x^ réalise le minimum de Ax b 2 si et seulement si x^ est solution de ^ A Ax^ Ab
t t (équation
normale). Il existe toujours au moins une solution. La solution est unique si et seulement si le rang de A est N.
Défauts numériques On appliquera la méthode des équations normales pour un nombre faible d'inconnues 2, 3, 4 ou 5. Un premier défaut est l'imprécision du calcul engendré par le produit matriciel A A t
. Par ailleurs, si la matrice A est creuse, A A t peut très bien ne pas l'être…
Méthode de la factorisation QR (rectangulaire)
Définition ( M > N )
AM (^) N QMMRN N, avec Q une matrice orthogonale et (^)
O
R
R
NN , avec R une matrice triangulaire
supérieure.
Théorème : minimiser 2 Ax b pour N supérieur à 3 ou 4
Soit M
t c
c Q b
2
1 avec c 1 de dimension N. Le vecteur x minimise 2 Ax b si et seulement si x est
solution de Rx c 1.
Obtention de la factorisation A = QR par Householder
Définition des matrices de Householder
Soit u un vecteur unitaire 1
2 u de dimension m. On définit la matrice de Householder de dimensio m m
t
H u Idm 2 u.u. H uest symétrique et orthogonale H u 2 1 . On note également
t
u m
vv H Id
en posant v 2 .u.
Propriété fondamentale des matrices de Householder
Soit a un vecteur de dimension m, non nul. Il existe H, une matrice de Householder, telle que O m
H a
avec un réel.
Preuve :
2 2
2
a
- On choisit le signe de : le signe opposé à a 1.
- On pose. 1 0 2
a .
a v
.
Principe de l'algorithme de factorisation
On veut obtenir la matrice RM (^) N à partir de la matrice A A 0. On procède colonne par colonne.
On rappelle la géométrie de
O
x
R
1
.
- Occupons nous de la première colonne de A 0. Il existe une matrice de Householder H 1 de dimension m
telle que H 1 A 0 A 1 , où la première colonne de A 1 est
1
.
