Metoda gradientului conjugat
În analiza numerică , metoda gradientului conjugat este un algoritm pentru rezolvarea sistemelor de ecuații liniare a căror matrice este simetric pozitivă definită . Această metodă, imaginată în 1950 simultan de Cornelius Lanczos , Eduard Stiefel și Magnus Hestenes , este o metodă iterativă care converge într-un număr finit de iterații (cel mult egal cu dimensiunea sistemului liniar). Cu toate acestea, marele său interes practic din punctul de vedere al timpului de calcul provine din faptul că o inițializare inteligentă (numită „precondiționare”) face posibilă conducerea în doar câteva pasaje la o estimare foarte apropiată de soluția exactă, aceasta de aceea, în practică, una se limitează la un număr de iterații mult mai mic decât numărul necunoscutelor.
Metoda gradientului biconjugat prevede o generalizare pentru matrici nesimetrice.
Principiu
Obiectivul este de a minimiza funcția în care A este o matrice pătrată simetrică definită pozitivă de mărimea n.
f:X↦12(LAX,X)-(b,X){\ displaystyle f: x \ mapsto {\ frac {1} {2}} (\ mathbf {A} x, x) - (b, x)}
Calculul arată că o soluție a problemei este soluția sistemului : într-adevăr, avem .
LAX=b{\ displaystyle \ mathbf {A} x = b}∇f(X)=LAX-b{\ displaystyle \ nabla f \ left (x \ right) = \ mathbf {A} xb}
Intuitiv, funcția f poate fi deci văzută ca un antiderivant al reziduului . Prin anularea gradientului lui f , obținem vectorul x care minimizează eroarea.
LAX-b{\ displaystyle \ mathbf {A} xb}
Metoda gradientului conjugat văzută ca o metodă directă
Reamintim că doi vectori diferiți de zero u și v sunt conjugați față de A dacă
tuTLAv=0.{\ displaystyle u ^ {\ mathrm {T}} \ mathbf {A} v = 0.}Știind că A este simetric pozitiv definit, deducem un produs scalar
⟨tu,v⟩LA: =⟨LAtu,v⟩=⟨tu,LATv⟩=⟨tu,LAv⟩=tuTLAv.{\ displaystyle \ langle u, v \ rangle _ {\ mathbf {A}}: = \ langle \ mathbf {A} {u}, {v} \ rangle = \ langle {u}, \ mathbf {A} ^ { \ mathrm {T}} {v} \ rangle = \ langle {u}, \ mathbf {A} {v} \ rangle = {u} ^ {\ mathrm {T}} \ mathbf {A} {v}.}Doi vectori sunt conjugați dacă sunt, prin urmare, ortogonali pentru acest produs scalar.
Conjugarea este o relație simetrică: dacă u este conjugat cu v pentru A , atunci v este conjugat cu u .
Să presupunem că { p k } este o succesiune de n direcții conjugate în perechi. Atunci { p k } formează o bază a lui R n , deci soluția x * a lui A x = b în această bază:
X∗=∑eu=1nuαeupeu{\ displaystyle x _ {*} = \ sum _ {i = 1} ^ {n} \ alpha _ {i} p_ {i}}Coeficienții sunt dați de
b=LAX∗=∑eu=1nuαeuLApeu.{\ displaystyle {b} = \ mathbf {A} {x} _ {*} = \ sum _ {i = 1} ^ {n} \ alpha _ {i} \ mathbf {A} {p} _ {i} .}
pkTb=pkTLAX∗=∑eu=1nuαeupkTLApeu=αkpkTLApk.{\ displaystyle {p} _ {k} ^ {\ mathrm {T}} {b} = {p} _ {k} ^ {\ mathrm {T}} \ mathbf {A} {x} _ {*} = \ sum _ {i = 1} ^ {n} \ alpha _ {i} {p} _ {k} ^ {\ mathrm {T}} \ mathbf {A} {p} _ {i} = \ alpha _ { k} {p} _ {k} ^ {\ mathrm {T}} \ mathbf {A} {p} _ {k}.}(pentru că sunt conjugate două câte două)
∀eu≠k,peu,pk{\ displaystyle \ forall i \ neq k, p_ {i}, p_ {k}}
αk=pkTbpkTLApk=⟨pk,b⟩⟨pk,pk⟩LA=⟨pk,b⟩‖pk‖LA2.{\ displaystyle \ alpha _ {k} = {\ frac {{p} _ {k} ^ {\ mathrm {T}} {b}} {{p} _ {k} ^ {\ mathrm {T}} \ mathbf {A} {p} _ {k}}} = {\ frac {\ langle {p} _ {k}, {b} \ rangle} {\, \, \, \ langle {p} _ {k} , {p} _ {k} \ rangle _ {\ mathbf {A}}}} = {\ frac {\ langle {p} _ {k}, {b} \ rangle} {\, \, \, \ | {p} _ {k} \ | _ {\ mathbf {A}} ^ {2}}}.}
Avem astfel ideea de ghidare a metodei de rezolvare a sistemului A x = b : găsiți o succesiune de n direcții conjugate și calculați coeficienții α k .
Metoda gradientului conjugat văzută ca o metodă iterativă
Alegând corect direcțiile conjugate p k , nu este necesar să le determinăm pe toate pentru a obține o bună aproximare a soluției x * . Este astfel posibil să considerăm metoda gradientului conjugat ca metodă iterativă. Această alegere face posibilă luarea în considerare a rezoluției unor sisteme foarte mari, unde calculul tuturor direcțiilor ar fi fost foarte lung.
Considerăm astfel un prim vector x 0 , pe care îl putem presupune a fi zero (în caz contrar, trebuie să considerăm sistemul A z = b - A x 0 ). Algoritmul va consta, începând de la x 0 , să „se apropie” de soluția necunoscută x * , ceea ce presupune definirea unei metrici. Această valoare provine din faptul că soluția x * este minimizatorul unic al formei pătratice :
f(X)=12XTLAX-XTb,X∈Rnu.{\ displaystyle f (\ mathbf {x}) = {\ frac {1} {2}} x ^ {\ mathrm {T}} \ mathbf {A} xx ^ {\ mathrm {T}} b, \ quad x \ in \ mathbb {R} ^ {n}.}Astfel, dacă f ( x ) scade după o iterație, atunci ne apropiem de x * .
Prin urmare, acest lucru sugerează luarea primei direcții p 1 ca opusă gradientului de la f la x = x 0 . Gradientul valorează A x 0 - b = - b , conform primei noastre ipoteze. Următorii vectori ai bazei vor fi astfel conjugați cu gradientul, de unde și denumirea de „metodă de gradient conjugat”.
Să r k să fie reziduul de la k - lea iterație:
rk=b-LAXk.{\ displaystyle r_ {k} = b- \ mathbf {A} x_ {k}. \,}Rețineți că r k este opusul gradientului lui f la x = x k , deci algoritmul de gradient indică deplasarea în direcția r k . Se reamintește că direcțiile p k sunt conjugate două câte două. De asemenea, dorim ca următoarea direcție să fie construită din reziduul curent și direcțiile construite anterior, ceea ce este o presupunere rezonabilă în practică.
Constrângerea de conjugare este o constrângere de ortonormalitate, deci problema împărtășește similarități cu metoda Gram-Schmidt .
Avem astfel
pk+1=rk-∑eu≤kpeuTLArkpeuTLApeupeu{\ displaystyle {p} _ {k + 1} = {r} _ {k} - \ sum _ {i \ leq k} {\ frac {{p} _ {i} ^ {\ mathrm {T}} \ mathbf {A} {r} _ {k}} {{p} _ {i} ^ {\ mathrm {T}} \ mathbf {A} {p} _ {i}}} {p} _ {i}}Urmând această direcție, următorul punct este dat de
Xk+1=Xk+αk+1pk+1{\ displaystyle {x} _ {k + 1} = {x} _ {k} + \ alpha _ {k + 1} {p} _ {k + 1}}
etapa
α k +1 este determinată astfel încât să se minimizeze :
g(α)=f(Xk+αpk+1){\ displaystyle g (\ alpha) = f ({x} _ {k} + \ alpha {p} _ {k + 1})}
g(α)=12α2pk+1TLApk+1+αpk+1T(LAXk-b)+vs.onustlanute{\ displaystyle g (\ alpha) = {\ frac {1} {2}} \ alpha ^ {2} {p} _ {k + 1} ^ {\ mathrm {T}} \ mathbf {A} {p} _ {k + 1} + \ alpha {p} _ {k + 1} ^ {\ mathrm {T}} (\ mathbf {A} {x} _ {k} -b) + constant}
minimul de
g este atins și ca
A este pozitiv definit ,
dgdα(αk+1)=0{\ displaystyle {\ frac {dg} {d \ alpha}} (\ alpha _ {k + 1}) = 0}pk+1TLApk+1>0{\ displaystyle {p} _ {k + 1} ^ {\ mathrm {T}} \ mathbf {A} {p} _ {k + 1}> 0}
prin urmare:
αk+1=pk+1T(b-LAXk)pk+1TLApk+1=pk+1Trkpk+1TLApk+1{\ displaystyle \ alpha _ {k + 1} = {\ frac {{p} _ {k + 1} ^ {\ mathrm {T}} (b- \ mathbf {A} {x} _ {k})} {{p} _ {k + 1} ^ {\ mathrm {T}} \ mathbf {A} {p} _ {k + 1}}} = {\ frac {{p} _ {k + 1} ^ { \ mathrm {T}} {r} _ {k}} {{p} _ {k + 1} ^ {\ mathrm {T}} \ mathbf {A} {p} _ {k + 1}}}}Algoritm
Pentru a începe recurența, este necesar să porniți de la o estimare inițială x 0 a vectorului căutat x ; și numărul de iterații N necesare astfel încât (unde ε este un număr pozitiv arbitrar aproape de zero) să depindă de x 0 ales. Din păcate, atât sigur , cât și general (adică eficiente pentru toate tipurile de matrice simetrice pozitive) metodele de „ precondiționare ” pentru formarea unui x 0 corect sunt ele însele costisitoare din punct de vedere al calculului. În practică, intuiția fizică, ghidată de natura fizică a problemei de rezolvat, sugerează uneori o inițializare eficientă: aceste idei au dat naștere de mai bine de treizeci de ani la o abundentă literatură de specialitate.
‖XNU-X‖<ε{\ displaystyle \ | x_ {N} -x \ | <\ varepsilon}
Algoritm iterativ pseudo-cod
Algoritmul de mai jos rezolvă A x = b , unde A este o matrice definită reală, simetrică și pozitivă. Vectorul de intrare x 0 poate fi o aproximare a soluției inițiale sau 0.
r0: =b-LAX0p0: =r0k: =0repetaαk: =rkTrkpkTLApkXk+1: =Xk+αkpkrk+1: =rk-αkLApkdacă rk+1 este suficient de mic, apoi ieșim din buclăβk: =rk+1Trk+1rkTrkpk+1: =rk+1+βkpkk: =k+1sfârșitul repetăriiRezultatul este Xk+1{\ displaystyle {\ begin {align} & \ mathbf {r} _ {0}: = \ mathbf {b} - \ mathbf {Ax} _ {0} \\ & \ mathbf {p} _ {0}: = \ mathbf {r} _ {0} \\ & k: = 0 \\ & {\ hbox {repeat}} \\ & \ qquad \ alpha _ {k}: = {\ frac {\ mathbf {r} _ { k} ^ {\ mathsf {T}} \ mathbf {r} _ {k}} {\ mathbf {p} _ {k} ^ {\ mathsf {T}} \ mathbf {Ap} _ {k}}} \ \ & \ qquad \ mathbf {x} _ {k + 1}: = \ mathbf {x} _ {k} + \ alpha _ {k} \ mathbf {p} _ {k} \\ & \ qquad \ mathbf { r} _ {k + 1}: = \ mathbf {r} _ {k} - \ alpha _ {k} \ mathbf {Ap} _ {k} \\ & \ qquad {\ hbox {si}} r_ {k + 1} {\ hbox {este suficient de mic, așa că ieșim din buclă}} \\ & \ qquad \ beta _ {k}: = {\ frac {\ mathbf {r} _ {k + 1} ^ {\ mathsf {T}} \ mathbf {r} _ {k + 1}} {\ mathbf {r} _ {k} ^ {\ mathsf {T}} \ mathbf {r} _ {k}}} \\ & \ qquad \ mathbf {p} _ {k + 1}: = \ mathbf {r} _ {k + 1} + \ beta _ {k} \ mathbf {p} _ {k} \\ & \ qquad k: = k + 1 \\ & {\ hbox {sfârșitul repetării}} \\ & {\ hbox {Rezultatul este}} \ mathbf {x} _ {k + 1} \ end {align}}}Convergenţă
Putem arăta următorul rezultat la convergența algoritmului:
‖X∗-Xk‖LA⩽2(κ(LA)-1κ(LA)+1)k‖X∗-X0‖LA,{\ displaystyle \ | x ^ {*} - x_ {k} \ | _ {\ mathbf {A}} \ leqslant 2 \ left ({\ frac {{\ sqrt {\ kappa (\ mathbf {A})}} -1} {{\ sqrt {\ kappa (\ mathbf {A})}} + 1}} \ right) ^ {k} \ | x ^ {*} - x_ {0} \ | _ {\ mathbf {A }},}unde desemnează condiționarea matricei șiκ(LA){\ displaystyle {\ sqrt {\ kappa (\ mathbf {A})}}}‖z‖LA=zTLAz.{\ displaystyle \ | z \ | _ {\ mathbf {A}} = {\ sqrt {z ^ {T} \ mathbf {A} z}}.}
Prin urmare, metoda gradientului conjugat are o convergență superliniară, care poate fi subminată prin condiționarea slabă a matricei. Cu toate acestea, rămâne mai bun decât algoritmii cu direcția unei pante mai abrupte .
Solver
-
(ro) M1CG1 - Un rezolvator de sisteme liniare simetrice prin iterații de gradient conjugat, folosind / construind un precondiționator BFGS / ℓ-BFGS . Scris în Fortran -77. Solverul are interesul de a oferi posibilitatea de a construi un precondiționator BFGS sau ℓ-BFGS (in) , care va putea fi util pentru rezoluția unui sistem liniar cu o matrice apropiată și un al doilea membru diferit.
Note și referințe
Note
-
Magnus Hestenes și Eduard Stiefel , „ Methods of Conjugate Gradients for Solving Linear Systems ”, Journal of Research of the National Bureau of Standards , vol. 49, nr . 6,1952( citiți online [PDF] )
-
Potrivit lui Dianne O'Leary (cf. bibliografie), articolul lui J. Meijerink și H. van der Vorst , „ O metodă de soluție iterativă pentru sistemele liniare ale căror matrice de coeficienți este simetrică o matrice M ”, Matematică de calcul , nr . 31,1977, p. 148-162marchează un pas decisiv în ideea precondiționării unui sistem liniar înainte de a aplica algoritmul acestuia. Acest articol pionier propunea precondiționarea prin descompunerea incompletă a LU . Urmat, printre alte condiții preliminare, de SOR întrerupt ( M. DeLong și J. Ortega , „ SRG ca o precondiție ” , Matematica numerică aplicată , nr . 18,1995, p. 431-440), și prin metoda Gauss-Seidel întreruptă ( Y. Saad și Gene Golub ( ed. ), precondiționari paraleli pentru matrice generale rare , Progrese recente în metodele iterative , Springer Verlag,1994, 165-199 p.). O privire de ansamblu asupra diferitelor tehnici poate fi găsită, printre altele: A. Bruaset , O anchetă a metodelor iterative precondiționate , Longman Scientific & Technical, col. "Note de cercetare Pitman în matematică",1995 ; J. Erhel și K. Burrage , Despre performanța diverselor strategii adaptive precondiționate GMRES , INRIA / IRISA,1997, etc.
Bibliografie
- Philippe Ciarlet , Introducere în analiza și optimizarea matricii numerice , Masson, col. „Matematică. Aplic. pentru Maeștri,1985( reeditare 2001) ( ISBN 2-225-68893-1 )
- Dianne P. O'Leary ( dir. ), Metode lineare și neliniare ale gradientului conjugat , AMS-SIAM,1996, „Conjugați gradientul și algoritmii KMP conexi: începuturile”
Vezi și tu
Articole similare
linkuri externe
<img src="https://fr.wikipedia.org/wiki/Special:CentralAutoLogin/start?type=1x1" alt="" title="" width="1" height="1" style="border: none; position: absolute;">