Algorithmes et techniques avancees de programmation

De
Publié par

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 ?


Publié le : vendredi 8 juin 2012
Lecture(s) : 110
Source : math.univ-lyon1.fr
Nombre de pages : 4
Voir plus Voir moins
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 AIRinqiuee´,edt´erteisvyemptorsiicenam.u Rappelonsdabordleprincipeetlesprincipalescaract´eristiquesdelalgorithmedugradientconjugu´e (GC)(pourlesde´veloppementsthe´oriquesdelam´ethode,nousrenvoyonsa`[?,?] par exemple). n T Sia, bIRsont deux vecteurs colonnes, on va noter(a, b) =a bleur produit scalaire euclidien. n Consid´eronslafonctionnellequadratiqueE: IRIR 1 E(x() =Ax, x)(b, x),(1.1) 2 dont le gradient vautE(x) = (Axb). Le minimum de cette fonctionnelle est atteint pour le n point¯xIRqui annuleE(x),dotna´vcnireAx¯ =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 avecdIRladirection de descenteetρIRlepas de descenteC.emmo´dtnineesrlpa ram`etresdelam´ethodededescente?Lar´eponseestsimplepourlechoixdupasdedescente,car, i i une foisdetxfixe´s, on peut calculer facilement un pas de descenteoptimal i i (Axb, 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)1inla base canonique deIR(on peut de´montrer quilsagitenfaitdelame´thodeit´erativedeGaussSeidel);
1
2 i ii nt:dugradie´mteohedald=g=E(x)tivaquonorlerbs´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. Lefcacite´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=Axb k j soient tous orthogonaux entre eux,i.e.(gg ,) = 0,0k < ji(ce n’est pas le cas dans la m´ethodedugradient,ou`seulementdeuxgradientssuccessifssontorthogonaux); i10 0i1 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 soientxIRvalela)eeetno´neld(tiairuniε´eprsicila)(noee´x 0 0 g=Axb 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
Silalgorithme(GC)peutˆetreconsid´er´ecommeuneme´thodeexacte(laconvergence´etantassure´e enne´ticesneettaadnea`tnemmmco´eisilutntvitare´tiedohte´s),iupluonsaratieleme´ar´gneelts quilconvergeplusrapidement.Pouracc´el´ererlavitessedeconvergence,onpeutco´eproinndntire lesyst`emelin´eaire.Lid´eeiciestdemultiplierlesyste`meinitialparunematriceC, dite de pr´econditionnementonelerduose´rte,e(t`emusysuveaCAx=Cb) par le (GC) classique. Entouterigueur,pourappliquerle(GC),lamatricedusyst`emedoiteˆtresyme´trique,d´eniepositive,conditionquinestpasassure´emˆemesilamatriceCest syme´trique, de´finie positive (le produitdedeuxmatricessym´etriques,d´eniepositives,nestpasn´ecessairementunematrice syme´trique,de´niepositive).Cetteobservationnousam`ene`afairelesmanipulationsalg´ebriques suivantes pour trouver une formulation simple de l’algorithme (GCP). n×n SoitUIRune matrice inversible etz=U xun changement de variable. La fonctionnelleE donn´eepar(1.1)devient:
1 1 111T1T E(z) =(U zAU z,)(b, Uz() =z, zU AU)(zU b,), 2 2
(1.4)
´ ´ 1.1. LAMETHODE DU GRADIENT CONJUGUE3 T1T T1∗ ∗ avecU= (U() =U)elleionnonctttefiyls´tnsed`neumaseeopaoecrrireC.A z=bde ∗ −T1∗ −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) : 0T1 0T T0 0 g=zU AUU bU g=Axb   0 G i+1iT1i i+1i1i ii g=g+ρU AUhG=G+ρA Uh=G+ρAH   i H 0 001T0 0 h=gH=GU U=CG   C i iii iT i (1.5) (hg ,) (U G, U H) (HG ,) ρ=− ⇔ρ==i iiT1ii i (Ahh ,) (HAU U, UU H) (H ,AH) i+1ii i+1i i z=z+ρhx=x+ρH i+1i+1T i+1T i+1i+1i+1 (gg ,) (U GU G ,) (G ,CG) γ=γ= = i iT ii iT i (g ,g) (G, UU G) (CGG ,) i+1i+1i i+1i+1i h=g+γhH=CG+γH
1T Observons que la matriceC=U Ucuetlesvetitq,eeinisopeuqie´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,ien´eedquri positive : 0n soientxIR,ε,Cn´esdon 0 0 G=Axb 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´etasidnoil,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.
Soyez le premier à déposer un commentaire !

17/1000 caractères maximum.