Algorithmes et techniques avancees de programmation

Marie Curie

Description

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

• matrice symetrique


• b?


• algorithme sans preconditionnement


• ?hi ?


• ax0 ?


• methode du gradient


• di ?

Sujets

• Matrice symétrique

Informations

• Publié par: cidur
• Nombre de lectures: 131
• Langue: Français

Extrait

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 n i avecd∈IRladirection de descenteetρ∈IRlepas de descenteC.emmo´dtninfieesrlpa 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 i i n •lam´ethoaxitno:edederald=e, avec(e) 1≤i≤n la base canonique deIR(on peut de´montrer qu’ils’agitenfaitdelame´thodeit´erativedeGaussSeidel);

1

2 i ii •nt:dugradie´mteohedald=g=∇E(x)tivaquon’orlerbs´sabusee(eualavelonafelrdc 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`adest,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: 0 n 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−T 0 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+1 i 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 : 0 n 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,tliturdiaivneitntcondien´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.