Konjugert gradientmetode
I numerisk analyse er konjugatgradientmetoden en algoritme for å løse systemer av lineære ligninger hvis matrise er symmetrisk positiv bestemt . Denne metoden, forestilt samtidig i 1950 av Cornelius Lanczos , Eduard Stiefel og Magnus Hestenes , er en iterativ metode som konvergerer i et endelig antall iterasjoner (høyst lik dimensjonen til det lineære systemet). Imidlertid kommer den store praktiske interessen fra beregningstidspunktet fra det faktum at en smart initialisering (kalt "forhåndskonditionering") gjør det mulig å føre på bare noen få passasjer til et estimat veldig nær den nøyaktige løsningen, dette er det i praksis en er begrenset til et antall iterasjoner som er mye lavere enn antall ukjente.
Den biconjugate gradientmetode gir en generalisering for usymmetriske matriser.
Prinsipp
Målet er å minimere funksjonen der A er en positiv bestemt symmetrisk firkantmatrise av størrelse n.
f:x↦12(PÅx,x)-(b,x){\ displaystyle f: x \ mapsto {\ frac {1} {2}} (\ mathbf {A} x, x) - (b, x)}
Beregning viser at en løsning på problemet er løsningen på systemet : ja, det har vi .
PÅx=b{\ displaystyle \ mathbf {A} x = b}∇f(x)=PÅx-b{\ displaystyle \ nabla f \ left (x \ right) = \ mathbf {A} xb}
Intuitivt kan funksjonen f derfor sees på som et antiderivativ for resten . Ved å kansellere gradienten til f , får vi vektoren x som minimerer feilen.
PÅx-b{\ displaystyle \ mathbf {A} xb}
Konjugatgradientmetoden sett på som en direkte metode
Vi husker at to ikke-nullvektorer u og v er konjugerte med hensyn til A if
uTPÅv=0.{\ displaystyle u ^ {\ mathrm {T}} \ mathbf {A} v = 0.}Når vi vet at A er symmetrisk positivt bestemt, utleder vi et skalarprodukt
⟨u,v⟩PÅ: =⟨PÅu,v⟩=⟨u,PÅTv⟩=⟨u,PÅv⟩=uTPÅv.{\ 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}.}To vektorer er konjugerte hvis de derfor er ortogonale for dette skalære produktet.
Bøyningen er et symmetrisk forhold: hvis u er konjugert til v for A , så blir v konjugert til u .
Anta at { p k } er en sekvens av n retninger konjugert i par. Deretter { p k } danner en basis for R n , slik at løsningen x * av A x = b i denne basis:
x∗=∑Jeg=1ikkeαJegsJeg{\ displaystyle x _ {*} = \ sum _ {i = 1} ^ {n} \ alpha _ {i} p_ {i}}Koeffisientene er gitt av
b=PÅx∗=∑Jeg=1ikkeαJegPÅsJeg.{\ displaystyle {b} = \ mathbf {A} {x} _ {*} = \ sum _ {i = 1} ^ {n} \ alpha _ {i} \ mathbf {A} {p} _ {i} .}
skTb=skTPÅx∗=∑Jeg=1ikkeαJegskTPÅsJeg=αkskTPÅsk.{\ 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}.}(fordi er konjugert to og to)
∀Jeg≠k,sJeg,sk{\ displaystyle \ forall i \ neq k, p_ {i}, p_ {k}}
αk=skTbskTPÅsk=⟨sk,b⟩⟨sk,sk⟩PÅ=⟨sk,b⟩‖sk‖PÅ2.{\ 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}}}.}
Vi har altså den ledende ideen om metoden for å løse systemet A x = b : finn en sekvens av n konjugerte retninger, og beregne koeffisientene α k .
Konjugatgradientmetoden sett på som en iterativ metode
Ved å velge konjugatretningene p k riktig , er det ikke nødvendig å bestemme dem alle for å oppnå en god tilnærming av løsningen x * . Det er således mulig å betrakte konjugatgradientmetoden som en iterativ metode. Dette valget gjør det dermed mulig å vurdere oppløsningen til veldig store systemer, der beregningen av alle retninger ville ha vært veldig lang.
Vi anser altså en første vektor x 0 , som vi kan anta å være null (ellers må vi vurdere systemet A z = b - A x 0 ). Algoritmen vil, fra x 0 , bestå i å "komme nærmere" den ukjente løsningen x * , som antar definisjonen av en beregning. Denne beregningen kommer fra det faktum at løsningen x * er den unike minimiseringen av kvadratisk form :
f(x)=12xTPÅx-xTb,x∈Rikke.{\ displaystyle f (\ mathbf {x}) = {\ frac {1} {2}} x ^ {\ mathrm {T}} \ mathbf {A} xx ^ {\ mathrm {T}} b, \ quad x \ in \ mathbb {R} ^ {n}.}Dermed, hvis f ( x ) avtar etter en iterasjon, så nærmer vi oss x * .
Dette antyder derfor å ta den første retningen p 1 som det motsatte av gradienten fra f til x = x 0 . Gradienten er verdt A x 0 - b = - b , ifølge vår første hypotese. Følgende vektorer av basen vil således bli konjugert til gradienten, derav navnet "konjugert gradientmetode".
La r k være residiet i den k- te iterasjon:
rk=b-PÅxk.{\ displaystyle r_ {k} = b- \ mathbf {A} x_ {k}. \,}Merk at r k er det motsatte av gradienten til f ved x = x k , så gradientalgoritmen indikerer å bevege seg i retningen r k . Det minnes at retningen p k er konjugert to og to. Vi ønsker også at følgende retning skal konstrueres fra strømresten og retningene som tidligere er konstruert, noe som er en rimelig antagelse i praksis.
Konjugasjonsbegrensningen er en orthonormalitetsbegrensning, så problemet deler likheter med Gram-Schmidt-metoden .
Det har vi altså
sk+1=rk-∑Jeg≤ksJegTPÅrksJegTPÅsJegsJeg{\ 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}}Etter denne retningen er det neste punktet gitt av
xk+1=xk+αk+1sk+1{\ displaystyle {x} _ {k + 1} = {x} _ {k} + \ alpha _ {k + 1} {p} _ {k + 1}}
trinnet
α k +1 blir bestemt for å minimere :
g(α)=f(xk+αsk+1){\ displaystyle g (\ alpha) = f ({x} _ {k} + \ alpha {p} _ {k + 1})}
g(α)=12α2sk+1TPÅsk+1+αsk+1T(PÅxk-b)+vs.oikkestpåikkete{\ 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) + konstant}
minimum
g er nådd, og som
A er positiv ,
dgdα(αk+1)=0{\ displaystyle {\ frac {dg} {d \ alpha}} (\ alpha _ {k + 1}) = 0}sk+1TPÅsk+1>0{\ displaystyle {p} _ {k + 1} ^ {\ mathrm {T}} \ mathbf {A} {p} _ {k + 1}> 0}
derfor :
αk+1=sk+1T(b-PÅxk)sk+1TPÅsk+1=sk+1Trksk+1TPÅsk+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}}}}Algoritme
For å starte gjentakelsen er det nødvendig å starte fra et innledende estimat x 0 av den søkte vektoren x ; og antall iterasjoner N nødvendig slik at (der ε er et positivt tall vilkårlig nær null) avhenger av valgt x 0 . Dessverre er både sikre og generelle (dvs. effektive for alle slags positive symmetriske matriser) " forhåndskondisjonerings " -metoder for å danne en korrekt x 0 i seg selv beregningsdyr. I praksis antyder fysisk intuisjon, styrt av den fysiske arten av problemet som skal løses, noen ganger en effektiv initialisering: disse ideene har gitt opphav til en rikelig spesialisert litteratur i mer enn tretti år.
‖xIKKE-x‖<ε{\ displaystyle \ | x_ {N} -x \ | <\ varepsilon}
Pseudokode iterativ algoritme
Algoritmen nedenfor løser A x = b , hvor A er en reell, symmetrisk og positiv bestemt matrise. Inngangsvektoren x 0 kan være en tilnærming til den opprinnelige løsningen eller 0.
r0: =b-PÅx0s0: =r0k: =0gjentaαk: =rkTrkskTPÅskxk+1: =xk+αkskrk+1: =rk-αkPÅskhvis rk+1 er liten nok, så går vi ut av løkkenβk: =rk+1Trk+1rkTrksk+1: =rk+1+βkskk: =k+1slutten av gjentakelsenResultatet er 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 {er liten nok, så vi går ut av løkken}} \\ & \ 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 {slutten av gjentakelse}} \\ & {\ hbox {Resultatet er}} \ mathbf {x} _ {k + 1} \ end {justert}}}Konvergens
Vi kan vise følgende resultat på konvergensen av algoritmen:
‖x∗-xk‖PÅ⩽2(κ(PÅ)-1κ(PÅ)+1)k‖x∗-x0‖PÅ,{\ 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 }},}hvor betegner kondisjonering av matrisen ogκ(PÅ){\ displaystyle {\ sqrt {\ kappa (\ mathbf {A})}}}‖z‖PÅ=zTPÅz.{\ displaystyle \ | z \ | _ {\ mathbf {A}} = {\ sqrt {z ^ {T} \ mathbf {A} z}}.}
Konjugatgradientmetoden har derfor en superlinjær konvergens, som kan undergraves av dårlig kondisjonering av matrisen. Imidlertid er det fortsatt bedre enn algoritmene med retning av brattere skråning .
Løser
-
(no) M1CG1 - En løser av symmetriske lineære systemer ved konjugerte gradient-iterasjoner, ved bruk / bygging av en BFGS / ℓ-BFGS-for-conditioner . Skrevet i Fortran -77. Løseren har interessen i å tilby muligheten for å bygge en for- condition BFGS eller ℓ-BFGS (in) , som vil være nyttig for oppløsningen av et lineært system med en nær matrise og et annet annet medlem.
Merknader og referanser
Merknader
-
Magnus Hestenes og Eduard Stiefel , “ Methods of Conjugate Gradients for Solving Linear Systems ”, Journal of Research of the National Bureau of Standards , vol. 49, n o 6,1952( les online [PDF] )
-
I følge Dianne O'Leary (jf. Bibliografi), artikkelen av J. Meijerink og H. van der Vorst , “ En iterativ løsningsmetode for lineære systemer der koeffisientmatrisen er symmetrisk en M-matrise ”, Matematikk av Computation , n o 311977, s. 148-162markerer et avgjørende trinn i ideen om å forkondisjonere et lineært system før algoritmen brukes på det. Denne banebrytende artikkelen foreslo forbehandling av ufullstendig LU-dekomponering . Etterfulgt blant annet forbehandling av SOR avbrutt ( M. DeLong og J. Ortega , " SRG som en forutsetning ," Applied Mathematics Numerisk , n o 181995, s. 431-440), og av Interrupted Gauss-Seidel Method ( Y. Saad og Gene Golub ( red. ), Parallelle forutsetninger for generelle sparsomme matriser , Recent Advances in Iterative Methods , Springer Verlag,1994, 165-199 s.). En oversikt over de forskjellige teknikkene finnes blant annet i: A. Bruaset , A survey of preconditioned iterative methods , Longman Scientific & Technical, coll. "Pitman Research Notes in Mathematics",1995 ; J. Erhel og K. Burrage , om utførelsen av ulike adaptive forhåndskonditionerte GMRES-strategier , INRIA / IRISA,1997, etc.
Bibliografi
- Philippe Ciarlet , Introduksjon til numerisk matriseanalyse og optimalisering , Masson, koll. "Matte. Appl. for mestrene,1985( opptrykk 2001) ( ISBN 2-225-68893-1 )
- Dianne P. O'Leary ( dir. ), Lineære og ikke-lineære konjugatgradientrelaterte metoder , AMS-SIAM,1996, "Konjugatgradient og relaterte KMP-algoritmer: begynnelsen"
Se også
Relaterte artikler
Eksterne linker
<img src="https://fr.wikipedia.org/wiki/Special:CentralAutoLogin/start?type=1x1" alt="" title="" width="1" height="1" style="border: none; position: absolute;">