Algorithmes et techniques avancees de programmation
4
pages
Français
Documents
Le téléchargement nécessite un accès à la bibliothèque YouScribe Tout savoir sur nos offres
4
pages
Français
Documents
Le téléchargement nécessite un accès à la bibliothèque YouScribe Tout savoir sur nos offres
Niveau: Supérieur, Master
Algorithmes et techniques avancees de programmation Universite Pierre et Marie Curie Master de Sciences et Technologies Mention Mathematiques et Applications Master I : Informatique scientifique en C++ Travaux diriges I. DANAILA extrait du livre : I. Danaila, F. Hecht, O. Pironneau, Simulation numerique en C++, Dunod, 2003 1.1 La methode du gradient conjugue Nous allons montrer dans cette section comment utiliser les classes RNM pour programmer l'algo- rithme du gradient conjugue preconditionne (GCP) pour resoudre le systeme lineaire Ax = b, avec A ? IRn?n une matrice symetrique, definie positive. Rappelons d'abord le principe et les principales caracteristiques de l'algorithme du gradient conjugue (GC) (pour les developpements theoriques de la methode, nous renvoyons a [?, ?] par exemple). Si a, b ? IRn sont deux vecteurs colonnes, on va noter (a, b) = aT b leur produit scalaire euclidien. Considerons la fonctionnelle quadratique E : IRn ? IR E(x) = 1 2 (Ax, x)? (b, x), (1.1) dont le gradient vaut ?E(x) = (Ax ? b). Le minimum de cette fonctionnelle est atteint pour le point x¯ ? IRn qui annule ?E(x), donc verifiant Ax¯ = b. La methode du gradient conjugue fait partie des methodes de descente, qui ont comme principe commun la recherche de x¯ suivant le procede iteratif x0 donne, xi+1 = xi + ?idi, (1.2) avec di ? IRn
Voir
Algorithmes et techniques avancees de programmation Universite Pierre et Marie Curie Master de Sciences et Technologies Mention Mathematiques et Applications Master I : Informatique scientifique en C++ Travaux diriges I. DANAILA extrait du livre : I. Danaila, F. Hecht, O. Pironneau, Simulation numerique en C++, Dunod, 2003 1.1 La methode du gradient conjugue Nous allons montrer dans cette section comment utiliser les classes RNM pour programmer l'algo- rithme du gradient conjugue preconditionne (GCP) pour resoudre le systeme lineaire Ax = b, avec A ? IRn?n une matrice symetrique, definie positive. Rappelons d'abord le principe et les principales caracteristiques de l'algorithme du gradient conjugue (GC) (pour les developpements theoriques de la methode, nous renvoyons a [?, ?] par exemple). Si a, b ? IRn sont deux vecteurs colonnes, on va noter (a, b) = aT b leur produit scalaire euclidien. Considerons la fonctionnelle quadratique E : IRn ? IR E(x) = 1 2 (Ax, x)? (b, x), (1.1) dont le gradient vaut ?E(x) = (Ax ? b). Le minimum de cette fonctionnelle est atteint pour le point x¯ ? IRn qui annule ?E(x), donc verifiant Ax¯ = b. La methode du gradient conjugue fait partie des methodes de descente, qui ont comme principe commun la recherche de x¯ suivant le procede iteratif x0 donne, xi+1 = xi + ?idi, (1.2) avec di ? IRn
- matrice symetrique
- b?
- algorithme sans preconditionnement
- ?hi ?
- ax0 ?
- methode du gradient
- di ?
Universite´ Pierre et Marie Curie Master de Sciences et Technologies Mention Mathe´matiques et Applications
Master I : Informatique scientifique en C++ Travauxdirig´es I. DANAILA
Algorithmesettechniquesavanc´eesde programmation
extrait du livre :I. Danaila, F. Hecht, O. Pironneau, Simulationnum´eriqueenC++,Dunod,2003
1.1Lam´ethodedugradientconjugu´e
Nous allons montrer dans cette section comment utiliser les classesRNMpour programmer l’algo rithmedugradientconjugue´pre´conditionn´e(GCP)pourre´soudrelesyste`meline´aireAx=b, avec n×n A∈IRfiinqiuee´,edt´erteisvyemptorsiicenam.u Rappelonsd’abordleprincipeetlesprincipalescaract´eristiquesdel’algorithmedugradientconjugu´e (GC)(pourlesde´veloppementsthe´oriquesdelam´ethode,nousrenvoyonsa`[?,?] par exemple). n T Sia, b∈IRsont deux vecteurs colonnes, on va noter(a, b) =a bleur produit scalaire euclidien. n Consid´eronslafonctionnellequadratiqueE: IR→IR 1 E(x() =Ax, x)−(b, x),(1.1) 2 dont le gradient vaut∇E(x) = (Ax−b). Le minimum de cette fonctionnelle est atteint pour le n point¯x∈IRqui annule∇E(x),dotna´vcnfiireAx¯ =b. Lam´ethodedugradientconjugue´faitpartiedesscenmt´eethodesdede, qui ont comme principe commun la recherche dex¯erateit´´ed´procfiaviueltns 0i+1ii i xdonne´, x=x+ρ d ,(1.2) i ni avecd∈IRladirection de descenteetρ∈IRlepas de descenteC.emmo´dtninfieesrlpa ram`etresdelam´ethodededescente?Lar´eponseestsimplepourlechoixdupasdedescente,car, i i une foisdetxfixe´s, on peut calculer facilement un pas de descenteoptimal i i (Ax−b, d) i ρ=−(1.3) i i (dAd ,) quir´ealiseleminimumdelafonctionnelleEdans la direction choisie, ou, plus formellement, le i i minimumdelafonctionre´elle,devariabler´eelle:ρ→E(x+ρd). Quant`aladirectiondedescente,plusieurschoixsontpossibles,cequinouspermetdedistinguer: i ii n •lam´ethoaxitno:edederald=e, avec(e)1≤i≤nla base canonique deIR(on peut de´montrer qu’ils’agitenfaitdelame´thodeit´erativedeGaussSeidel);
1
2 i ii •nt:dugradie´mteohedald=g=∇E(x)tivaquon’orlerbs´sabusee(eualavelonafelrdc tionnelle diminue le plus rapidement suivant la direction du gradient) ; i i •dientconjugu´e:malhte´dedoargud=h, avec les directions de descenteconjugue´espar i j rapport`alamatriceAire`adest,c’(Ah ,h) = 0, i=j. L’efficacite´delame´thodedugradientconjugu´er´esidedanssespropri´et´esremarquables,e´nonce´es ici pour l’ite´rationi: i ii •les directions de descentehsont construites de telle manie`re que les gradientsg=Ax−b k j soient tous orthogonaux entre eux,i.e.(gg ,) = 0,∀0≤k < j≤i(ce n’est pas le cas dans la m´ethodedugradient,ou`seulementdeuxgradientssuccessifssontorthogonaux); i10 0i−1 •xellneiseleminr´ealofcnitnomimuedalEsur l’espace vectorielx+V ect{g ,g ,. . . , g}, n de dimensioni; quandi=n, cet espace vectoriel sera exactementIR, ce qui prouve que l’algorithme converge enn;ite´rations au plus i •d’un point de vue pratique, les directions de descentehsedrirtpa`ar,lecualacel`saficostn i gradientsgpusse´lpcevsruetntmeoitresedleeuseotkcgaecsstilethoden´e,etlam´eiaeremtns (voir programme plus loin). Une formulation´eisimptoede l’algorithme (GC) est la suivante :
Algorithme 1.1
Legradientconjugu´epourr´esoudrelesyst`emeline´aireAx=b: 0n soientx∈IRvalela)eeetno´neld(tiairuniε´eprsicila)fi(noee´x 0 0 g=Ax−b 0 0 h=−g –pouri= 0`an i i (hg ,) ρ=−(pas optimal, cf. 1.3) i i (h ,Ah) i+1i ii x=x+ρh(avancement suivanth, cf. 1.2) i+1i i g=g+ρAheitn)c(uclat´liaterduifadgr i+1i+1 (gg ,) i+1i γ=(parame`tre assurant(hh ,) = 0) i i (g ,g) i+1i+1i h=−g+γhc(ucla´tilaterduifdilactredeseoidn)eectn i+1i+1 si(gg ,)< εstop
Sil’algorithme(GC)peutˆetreconsid´er´ecommeuneme´thodeexacte(laconvergence´etantassure´e enne´ticesneetta’adnea`tnemmmco´eisilutntvitare´tiedohte´s),iupluonsaratieleme´ar´gneelts qu’ilconvergeplusrapidement.Pouracc´el´ererlavitessedeconvergence,onpeutco´eproinndntire lesyst`emelin´eaire.L’id´eeiciestdemultiplierlesyste`meinitialparunematriceC, dite de pr´econditionnementonelerduose´rte,e(t`emusysuveaCAx=Cb) par le (GC) classique. Entouterigueur,pourappliquerle(GC),lamatricedusyst`emedoiteˆtresyme´trique,d´efiniepo sitive,conditionquin’estpasassure´emˆemesilamatriceCest syme´trique, de´finie positive (le produitdedeuxmatricessym´etriques,d´efiniepositives,n’estpasn´ecessairementunematrice syme´trique,de´finiepositive).Cetteobservationnousam`ene`afairelesmanipulationsalg´ebriques suivantes pour trouver une formulation simple de l’algorithme (GCP). n×n SoitU∈IRune matrice inversible etz=U xun changement de variable. La fonctionnelleE donn´eepar(1.1)devient:
1 1 −1−1−1−T−1−T E(z) =(U zAU z,)−(b, Uz() =z, zU AU)−(zU b,), 2 2
(1.4)
´ ´ 1.1. LAMETHODE DU GRADIENT CONJUGUE3 −T−1T T−1∗ ∗ avecU= (U() =U)elleionnonctttefiyls´tnsed`neumaseeopaoecrrireC.A z=bde ∗ −T−1∗ −T matriceA=U AUsyme´trique, de´finie positive, et de second membreb=U b. Nous pouvonsappliquerle(GC)pourcesyst`eme,etlesprincipalesquantite´s`acalculerdeviennent(cf. algorithme 1.1) : 0−T−1 0−T T0 0 g=zU AU−U b⇔U g=Ax−b 0 G i+1i−T−1i i+1i−1i ii g=g+ρU AUh⇔G=G+ρA Uh=G+ρAH i H 0 00−1−T0 0 h=−g⇔H=−GU U=−CG C i i−ii iT i (1.5) (hg ,) (U G, U H) (HG ,) ρ=− ⇔ρ=−=− i ii−T−1ii i (Ahh ,) (HAU U, UU H) (H ,AH) i+1ii i+1i i z=z+ρh⇔x=x+ρH i+1i+1−T i+1−T i+1i+1i+1 (gg ,) (U GU G ,) (G ,CG) γ=⇔γ= = i i−T i−i iT i (g ,g) (G, UU G) (CGG ,) i+1i+1i i+1i+1i h=−g+γh⇔H=−CG+γH
−1−T Observons que la matriceC=U Uc’uetlesvetitq,eeinfiisopeuqie´d,tsym´etreselsaue matrice intervenant dans les expressions pre´ce´dentes (la matriceUetspeuocrnotsurriestutilis´eeju n l’algorithme). Si on note(a, b)C= (Ca, b)le produit scalaire induit parCsurIR, l’algorithme (GCP) s’e´crit sous la forme :
Algorithme 1.2
Legradientconjugu´epr´econditionn´eparunematriceCte´mys,iefin´eedquri positive : 0n soientx∈IR,ε,Cn´esdon 0 0 G=Ax−b 0 0 H=−CG – pouri= 0`an i i (HG ,) ρ=− i i (AHH ,) i+1i i x=x+ρH i+1i i G=G+ρAH i+1i+1 (G ,G)C γ= i i (G ,G)C i+1i+1i H=−CG+γH i+1i+1 si(G ,G)C< εstop
−1 Il est clair de ce qui pre´ce`de que le choix«id´lae»pour la matriceC, seraitC=A, l’inverse deAvidemmencecas,´e;adsnedodraugamelth´etasidnoi’l,tliturdiaivneitntcondien´edejugu −1 utile. En pratique, plus la matrice de pre´conditionnement est«proche”»deA, meilleure sera la convergence de l’algorithme. En meˆme temps, la matriceClpsueralˆttediouslaplleetsimp creuse possible, ce qui impose comme choix les plus faciles •C=I, la matrice identite´ (algorithme sans pre´conditionnement); −1 •C=D, l’inverse de la partie diagonale deA.
