Konjugerad gradientmetod

En jämförelse av konvergensen av gradientnedstigning med optimal stegstorlek (i grönt) och konjugerad vektor (i rött) för att minimera en kvadratisk funktion associerad med ett givet linjärt system. Konjugerad gradient, om man antar exakt aritmetik, konvergerar i högst n steg, där n är storleken på systemets matris (här n = 2).

I matematik är den konjugerade gradientmetoden en algoritm för den numeriska lösningen av särskilda system av linjära ekvationer , nämligen de vars matris är positiv-definitiv . Den konjugerade gradientmetoden implementeras ofta som en iterativ algoritm , tillämpbar på glesa system som är för stora för att kunna hanteras av en direkt implementering eller andra direkta metoder som Cholesky-nedbrytningen . Stora glesa system uppstår ofta när man numeriskt löser partiella differentialekvationer eller optimeringsproblem.

Den konjugerade gradientmetoden kan också användas för att lösa obegränsade optimeringsproblem som energiminimering . Den tillskrivs vanligtvis Magnus Hestenes och Eduard Stiefel , som programmerade den på Z4:an och undersökte den mycket.

Den bikonjugerade gradientmetoden ger en generalisering till icke-symmetriska matriser. Olika olinjära konjugerade gradientmetoder söker ett minimum av olinjära optimeringsproblem.

Beskrivning av problemet med konjugerade gradienter

Antag att vi vill lösa systemet med linjära ekvationer

$\mathbf {A} \mathbf {x} =\mathbf {b}$

för vektorn $\mathbf {x}$ , där den kända $n\times n$ matrisen $\mathbf {A}$ är symmetrisk (dvs. A ^T = A ), positiv-definit (dvs x ^T Ax > 0 för alla vektorer som inte är noll $\mathbf {x}$ i R ⁿ ), och real , och $\mathbf {b}$ är också kända. Vi betecknar den unika lösningen för detta system med $\mathbf {x} _{*}$ .

Härledning som en direkt metod

Huvudartikel: Härledning av konjugatgradientmetoden

Den konjugerade gradientmetoden kan härledas från flera olika perspektiv, inklusive specialisering av den konjugerade riktningsmetoden för optimering, och variation av Arnoldi / Lanczos iterationen för egenvärdesproblem . Trots skillnader i deras tillvägagångssätt delar dessa härledningar ett gemensamt ämne - vilket bevisar ortogonaliteten hos residualerna och konjugationen av sökriktningarna. Dessa två egenskaper är avgörande för att utveckla metodens välkända kortfattade formulering.

Vi säger att två vektorer som inte är noll, u och v är konjugerade (med avseende på ${\displaystyle \mathbf {A} } )$ om

$\mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =0.$

Eftersom $\mathbf {A}$ är symmetrisk och positiv-definitiv, definierar den vänstra sidan en inre produkt

$\mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} } :=\langle \mathbf {A} \mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} ^{\mathsf {T}}\mathbf {v} \ rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle .$

Två vektorer är konjugerade om och endast om de är ortogonala med avseende på denna inre produkt. Att vara konjugerat är en symmetrisk relation: om $\mathbf {u}$ är konjugerad till $\mathbf {v}$ , då är $\mathbf {v}$ konjugerad till $\mathbf {u}$ . Anta att

$P=\{\mathbf {p} _{1},\dots ,\mathbf {p} _{n}\}$

är en uppsättning av $n$ ömsesidigt konjugerade vektorer med avseende på $\mathbf {A}$ , dvs $\mathbf {p} _{i}^{\ mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0$ för alla $i\neq j$ . Då $P$ en grund för $\mathbb {R} ^{n}$ , och vi kan uttrycka lösningen $\mathbf {x} _{*}$ av $\mathbf {Ax} =\mathbf {b}$ på denna grund:

$\mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}\Rightarrow \mathbf {A} \mathbf { x} _{*}=\summa _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}.$

Vänstermultiplicera problemet $\mathbf {Ax} =\mathbf {b}$ med vektorn $\mathbf {p} _{k}^{\mathsf {T}}$ ger

$\mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} =\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {x} _{*}=\summa _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}=\sum _{i=1}^{n}\alpha _{i}\left\langle \mathbf {p} _{k},\mathbf {p} _{i }\right\rangle _{\mathbf {A} }=\alpha _{k}\left\langle \mathbf {p} _{k},\mathbf {p} _{k}\right\rangle _{\ mathbf {A} }$

och så

$\alpha _{k}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\langle \mathbf {p} _{k},\mathbf { p} _{k}\rangle _{\mathbf {A} }}}.$

Detta ger följande metod för att lösa ekvationen Ax = b : hitta en sekvens av $n$ konjugerade riktningar och beräkna sedan koefficienterna $\alpha _{k}$ .

Som en iterativ metod

Om vi ​​väljer de konjugerade vektorerna $\mathbf {p} _{k}$ noggrant, behöver vi kanske inte alla för att få en bra approximation till lösningen $\mathbf {x} _ {*}$ . Så vi vill betrakta den konjugerade gradientmetoden som en iterativ metod. Detta gör att vi också kan ungefärligen lösa system där n är så stort att den direkta metoden skulle ta för mycket tid.

₀₀ Vi betecknar den initiala gissningen för x _∗ med x ₀ (vi kan utan förlust av generalitet anta att ₀ x = 0 , annars betrakta systemet Az = b − Ax istället). Börjar med x söker vi efter lösningen och i varje iteration behöver vi ett mått för att tala om för oss om vi är närmare lösningen x _∗ (det är okänt för oss). Detta mått kommer från det faktum att lösningen x _∗ också är den unika minimeraren för följande kvadratiska funktion

$f(\mathbf {x} )={\tfrac {1}{2}}\mathbf {x} ^{\mathsf {T}}\mathbf {A} \mathbf {x} -\mathbf {x } ^{\mathsf {T}}\mathbf {b} ,\qquad \mathbf {x} \in \mathbf {R} ^{n}\,.$

Förekomsten av en unik minimerare är uppenbar eftersom dess hessiska matris av andraderivator är symmetrisk positiv-definitiv

$\mathbf {H} (f(\mathbf {x} ))=\mathbf {A} \,,$

och att minimeraren (använd D f ( x )=0) löser det initiala problemet följer av dess första derivata

$\nabla f(\mathbf {x} )=\mathbf {A} \mathbf {x} -\mathbf {b} \,.$

₀₀₀₀₀₀ Detta tyder på att den första basvektorn p är den negativa av gradienten för f vid x = x . Gradienten för f är lika med Ax − b . Om man börjar med en initial gissning x , betyder det att vi tar p = b − Ax . De andra vektorerna i basen kommer att vara konjugerade till gradienten, därav namnet conjugate gradient method . Observera att p också är den rest som tillhandahålls av detta initiala steg i algoritmen.

Låt r _k vara resten i det k: te steget:

$\mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}.$

Som observerats ovan är $\mathbf {r} _{k}$ den negativa gradienten för $f$ vid ${\displaystyle \mathbf {x} _{k}} ,$ så gradienten nedstigningsmetod skulle kräva att röra sig i riktningen r _k . Här insisterar vi dock på att riktningarna $\mathbf {p} _{k}$ måste vara konjugerade till varandra. Ett praktiskt sätt att upprätthålla detta är att kräva att nästa sökriktning byggs av den aktuella resterande och alla tidigare sökriktningar. Konjugationsbegränsningen är en begränsning av ortonormaltyp och därför kan algoritmen ses som ett exempel på Gram-Schmidt-ortonormalisering . Detta ger följande uttryck:

$\mathbf {p} _{k}=\mathbf {r} _{k}-\sum _ {i<k}{\frac {\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {r} _{k}}{\mathbf {p} _{i }^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}}}\mathbf {p} _{i}$

(se bilden överst i artikeln för effekten av konjugationsbegränsningen på konvergens). Efter denna riktning ges nästa optimala läge av

$\mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{ k}$

med

$\alpha _{k}={\frac {\mathbf {p } _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\ mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}},$

där den sista likheten följer av definitionen av $\mathbf {r} _{k}$ . Uttrycket för $\alpha _{k}$ kan härledas om man ersätter uttrycket för x _{k +1} med f och minimerar det med avseende på $\alpha _{k}$

${\begin{aligned}f(\mathbf {x} _{k+1})&=f(\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{ k})=:g(\alpha _{k})\\g'(\alpha _{k})&{\overset {!}{=}}0\quad \Rightarrow \quad \alpha _{k} ={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{\mathbf {p} _{k} ^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}\,.\end{aligned}}$

Den resulterande algoritmen

Ovanstående algoritm ger den mest enkla förklaringen av den konjugerade gradientmetoden. Till synes kräver algoritmen som sagt lagring av alla tidigare sökriktningar och restvektorer, såväl som många matris-vektormultiplikationer, och kan därför vara beräkningsmässigt dyr. En närmare analys av algoritmen visar dock att $\mathbf {r} _{i}$ är ortogonal mot $\mathbf {r} _{j}$ , dvs ${\displaystyle \mathbf {r} _{i}^{\mathsf {T}}\mathbf {r} _{j}=0} ,$ för i ≠ j. Och $\mathbf {p} _{i}$ är $\mathbf {A}$ -ortogonal mot $\mathbf {p} _{j}$ , dvs ${\displaystyle \mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0} , för i ≠$ j $i \neq j}$ . Detta kan anses vara att allteftersom algoritmen fortskrider, spänner $\mathbf {p} _{i}$ och $\mathbf {r} _{i}$ över samma Krylov-delrum . Där $\mathbf {r} _{i}$ utgör den ortogonala basen med avseende på standardinre produkten, och $\mathbf {p} _{i}$ utgör den ortogonala basen med avseende på till den inre produkten inducerad av $\mathbf {A}$ . Därför $\mathbf {x} _{k}$ betraktas som projektionen av $\mathbf {x}$ på Krylov-underrummet.

0 Algoritmen beskrivs nedan för att lösa Ax = b där $\mathbf {A}$ är en reell, symmetrisk, positiv-definitiv matris. Ingångsvektorn $\mathbf {x} _{0}$ kan vara en ungefärlig initiallösning eller . Det är en annan formulering av den exakta proceduren som beskrivs ovan.

${\begin{aligned}&\mathbf {r} _{0}: =\mathbf {b} -\mathbf {Ax} _{0}\\&{\hbox{if }}\mathbf {r} _{0}{\text{ är tillräckligt liten, returnera sedan }}\mathbf { x} _{0}{\text{ som resultat}}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\text {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{if }}\mathbf {r} _{k+1}{\text{ är tillräckligt liten, avsluta sedan loopen}}\ \&\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\\&{\text{slut upprepa}}\\&{\text{ returnera }}\mathbf {x} _{k+1}{\text{ som resultat}}\end{aligned}}$

Detta är den vanligaste algoritmen. Samma formel för β _k används också i Fletcher–Reeves olinjära konjugerade gradientmetod .

Startar om

Vi noterar att $\mathbf {x} _{1}$ beräknas med metoden för gradientnedstigning som tillämpas på $\mathbf {x} _{0}$ . Inställning av $\beta _{k}=0$ skulle på liknande sätt göra $\mathbf {x} _{k+1}$ beräknad med metoden gradient descent från $\mathbf {x} _{k}$ , dvs kan användas som en enkel implementering av en omstart av de konjugerade gradientiterationerna. Omstarter kan sakta ner konvergensen, men kan förbättra stabiliteten om konjugatgradientmetoden inte beter sig, t.ex. på grund av avrundningsfel .

Explicit restberäkning

Formlerna $\mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf { p} _{k}$ och ${\displaystyle \mathbf {r} _{k}:=\mathbf {b} -\mathbf {Ax} _{k}} , som$ båda håll i exakt aritmetik, gör formlerna $\mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}$ och $\mathbf {r} _{k+1}:=\mathbf {b} - \mathbf {Ax} _{k+1}$ matematiskt ekvivalent. Den förra används i algoritmen för att undvika en extra multiplikation med $\mathbf {A}$ eftersom vektorn $\mathbf {Ap} _{k}$ redan är beräknad för att utvärdera ${ \displaystyle \alpha _{k}}$ . Den senare kan vara mer exakt och ersätter den explicita beräkningen $\mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax } _{k+1}$ för den implicita av rekursionen föremål för ackumulering av avrundningsfel, och rekommenderas därför för en tillfällig utvärdering.

En norm för restvärdet används vanligtvis för stoppkriterier. Normen för den explicita residualen $\mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+ 1}$ ger en garanterad noggrannhetsnivå både i exakt aritmetik och i närvaro av avrundningsfel, där konvergensen stagnerar naturligt. Däremot implicita residual $\mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{ k}\mathbf {Ap} _{k}$ är känd för att bli mindre i amplitud långt under nivån för avrundningsfel och kan därför inte användas för att bestämma stagnationen av konvergens.

Beräkning av alfa och beta

I algoritmen är α _k vald så att $\mathbf {r} _{k+1}$ är ortogonal mot $\mathbf {r} _{k}$ . Nämnaren förenklas från

$\alpha _{k}={\frac {\mathbf {r} _{k}^ {\mathsf {T}}\mathbf {r} _{k}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}} }={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T} }\mathbf {Ap} _{k}}}$

eftersom $\mathbf {r} _{k+1}=\mathbf {p} _{k+1}-\mathbf {\beta } _{ k}\mathbf {p} _{k}$ . β { _k är vald så att $\mathbf {p} _{k+1}$ är konjugerat till $\displaystyle \mathbf {p} _{k}}$ . Inledningsvis är βk _{_}

$\beta _{k}=-{\frac {\mathbf {r} _{k+1}^{\mathsf {T }}\mathbf {A} \mathbf {p} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}} }$

använder sig av

$\mathbf {r} _{k+1}=\mathbf {r} _{k}-\alpha _{k}\mathbf {A} \ mathbf {p} _{k}$

och motsvarande

$\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}} (\mathbf {r} _{k}-\mathbf {r} _{k+1}),$

täljaren för β _k skrivs om som

$\mathbf {r} _{ k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k +1}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})=-{\frac {1}{\alpha _{k}} }\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}$

eftersom $\mathbf {r} _{k+1}$ och $\mathbf {r} _{k}$ är ortogonala till sin design. Nämnaren skrivs om som

${ \displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}=(\mathbf {r} _{k}+\beta _{k -1}\mathbf {p} _{k-1})^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k }}}\mathbf {r} _{k}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})={\frac {1} {\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}$

med att sökriktningarna pk är konjugerade och återigen att residualerna är ortogonala _. Detta ger β i algoritmen efter att ha avbrutit α _k .

Exempelkod i MATLAB / GNU Octave

   
          
      
        

       
            
              
              
              
            
           
            
        
                
          
    
 funktion  x  =  konjgrad  (  A, b, x  )  r  =  b  -  A  *  x  ;  p  =  r  ;  rsold  =  r  '  *  r  ;  för  i  =  1  :  längd  (  b  )  Ap  =  A  *  p  ;  alpha  =  rsold  /  (  p  '  *  Ap  );  x  =  x  +  alfa  *  p  ;  r  =  r  -  alfa  *  Ap  ;  rsnew  =  r  '  *  r  ;  if  sqrt  (  rsnew  )  <  1e-10  break  end  p  =  r  +  (  rsnew  /  rsold  )  *  p  ;  rsold  =  rsnew  ;  slutändan  _ 

Numeriskt exempel

Betrakta det linjära systemet Ax = b givet av

$\mathbf {A} \mathbf {x} ={\begin{bmatrix}4&1\\1&3\end{bmatrix} }{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}},$

vi kommer att utföra två steg av den konjugerade gradientmetoden som börjar med den första gissningen

$\mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}$

för att hitta en ungefärlig lösning på systemet.

Lösning

För referens är den exakta lösningen

$\mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\\\{\frac {7}{11 }}\end{bmatrix}}\approx {\begin{bmatrix}0.0909\\\\0.6364\end{bmatrix}}$

₀₀₀₀ Vårt första steg är att beräkna restvektorn r associerad med x . Denna residual beräknas från formeln r = b - Ax , och är i vårt fall lika med

$\mathbf {r} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}-{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{ bmatrix}2\\1\end{bmatrix}}={\begin{bmatrix}-8\\-3\end{bmatrix}}=\mathbf {p} _{0}.$

₀₀ Eftersom detta är den första iterationen kommer vi att använda restvektorn r som vår initiala sökriktning p ; metoden för att välja p _k kommer att ändras i ytterligare iterationer.

Vi beräknar nu skalären α ₀ med hjälp av sambandet

$\alpha _{0}={\frac {\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}{\mathbf {p} _{0}^{\mathsf {T}}\mathbf {Ap} _{0}}}={\frac {{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end {bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3 \end{bmatrix}}}}={\frac {73}{331}}\approx 0,2205$

Vi kan nu beräkna x ₁ med formeln

$\mathbf {x} _{1}=\mathbf {x} _{0}+\alpha _{0}\mathbf {p} _{0}={\begin{bmatrix}2\\1\ end{bmatrix}}+{\frac {73}{331}}{\begin{bmatrix}-8\\-3\end{bmatrix}}\approx {\begin{bmatrix}0.2356\\0.3384\end{bmatrix}} }}.$

Detta resultat fullbordar den första iterationen, resultatet är en "förbättrad" ungefärlig lösning på systemet, x ₁ . Vi kan nu gå vidare och beräkna nästa restvektor r ₁ med hjälp av formeln

$\mathbf {r} _{1}=\mathbf {r} _{0}-\alpha _{0}\mathbf {A} \mathbf {p} _{0}={\begin{bmatrix} -8\\-3\end{bmatrix}}-{\frac {73}{331}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\- 3\end{bmatrix}}\approx {\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.$

Vårt nästa steg i processen är att beräkna den skalära β ₀ som så småningom kommer att användas för att bestämma nästa sökriktning p ₁ .

$\beta _{0}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {r} _{0 }^{\mathsf {T}}\mathbf {r} _{0}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}- 0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}=0.0088 .$

Nu, med hjälp av denna skalära β ₀ , kan vi beräkna nästa sökriktning p ₁ med hjälp av relationen

$\mathbf {p} _{1}=\mathbf {r} _{1}+\beta _{0}\mathbf {p} _{0}\approx {\begin{bmatrix}-0.2810\\ 0.7492\end{bmatrix}}+0.0088{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}.$

Vi beräknar nu skalären α ₁ med vår nyförvärvade p ₁ med samma metod som den som används för α ₀ .

$\alpha _{1}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {p} _{1 }^{\mathsf {T}}\mathbf {Ap} _{1}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}- 0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-0.3511&0.7229\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix} }-0.3511\\0.7229\end{bmatrix}}}}=0.4122.$

Slutligen hittar vi x ₂ med samma metod som den som användes för att hitta x ₁ .

$\mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}\approx {\begin{bmatrix}0.2356\\0.3384 \end{bmatrix}}+0.4122{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}={\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.$

₀ Resultatet, x ₂ , är en "bättre" approximation av systemets lösning än x ₁ och x . Om exakt aritmetik skulle användas i detta exempel istället för begränsad precision, så skulle den exakta lösningen teoretiskt ha uppnåtts efter n = 2 iterationer ( n är systemets ordning).

Konvergensegenskaper

Den konjugerade gradientmetoden kan teoretiskt ses som en direkt metod, eftersom den i frånvaro av avrundningsfel producerar den exakta lösningen efter ett ändligt antal iterationer, som inte är större än matrisens storlek. I praktiken erhålls aldrig den exakta lösningen eftersom konjugatgradientmetoden är instabil med avseende på även små störningar, t.ex. är de flesta riktningar i praktiken inte konjugerade, på grund av en degenerativ karaktär av att generera Krylov-underrymden.

Som en iterativ metod förbättrar den konjugerade gradientmetoden monotont (i energinormen) approximationerna $\mathbf {x} _{k}$ till den exakta lösningen och kan nå den erforderliga toleransen efter en relativt liten (jämfört med problemets storlek) antal iterationer. Förbättringen är vanligtvis linjär och dess hastighet bestäms av villkorsnumret $\kappa (A)$ för systemmatrisen $A$ : desto större $\kappa ( A)$ är, desto långsammare är förbättringen.

Om $\kappa (A)$ är stor, används förkonditionering vanligtvis för att ersätta det ursprungliga systemet $\mathbf {Ax} -\mathbf {b} =0$ med $\mathbf {M} ^{-1}(\mathbf {Ax} -\mathbf {b} )=0$ så att $\kappa (\mathbf {M} ^{-1}\mathbf {A} )$ är mindre än $\kappa (\mathbf {A} )$ , se nedan.

Konvergenssats

Definiera en delmängd av polynom som

$\Pi _{k}^{*}:=\left\lbrace \ p\in \Pi _{k}\ :\ p(0)=1\ \right\rbrace \,,$

där $\Pi _{k}$ är mängden polynom med maximal grad $k$ .

Låt $\left(\mathbf {x} _{k}\right)_{k}$ vara de iterativa approximationerna av den exakta lösningen $\mathbf {x} _{* }$ , och definiera felen som $\mathbf {e} _{k}:=\mathbf {x} _{k}-\mathbf {x} _{*}$ . Nu kan konvergenshastigheten approximeras som

${\begin{aligned}\left\|\mathbf {e} _{k}\right \|_{\mathbf {A} }&=\min _{p\in \Pi _{k}^{*}}\left\|p(\mathbf {A} )\mathbf {e} _{0 }\right\|_{\mathbf {A} }\\&\leq \min _{p\in \Pi _{k}^{*}}\,\max _{\lambda \in \sigma (\ mathbf {A} )}|p(\lambda )|\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\left({\ frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\ \left \|\mathbf {e} _{0}\right\|_{\mathbf {A} }\,,\end{aligned}}$

där $\sigma (\mathbf {A} )$ anger spektrumet och $\kappa (\mathbf {A} )$ anger villkorsnumret .

Observera att den viktiga gränsen när $\kappa (\mathbf {A} )$ tenderar att $\infty$

${\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\approx 1-{ \frac {2}{\sqrt {\kappa (\mathbf {A} )}}}\quad {\text{for}}\quad \kappa (\mathbf {A} )\gg 1\,.$

Denna gräns visar en snabbare konvergenshastighet jämfört med de iterativa metoderna för Jacobi eller Gauss–Seidel som skalas som $\approx 1-{\frac {2}{\kappa (\mathbf {A) } )}}$ .

Inget avrundningsfel antas i konvergenssatsen, men konvergensgränsen är allmänt giltig i praktiken som teoretiskt förklaras av Anne Greenbaum .

Praktisk konvergens

Om det initieras slumpmässigt är det första steget av iterationer ofta det snabbaste, eftersom felet elimineras inom Krylov-underrummet som initialt återspeglar ett mindre effektivt villkorsnummer. Det andra steget av konvergens är vanligtvis väl definierat av den teoretiska konvergensen bunden med ${\textstyle {\sqrt {\kappa (\mathbf {A} )}}} ,$ men kan vara superlinjär, beroende på en fördelning av spektrumet för matrisen $A$ och spektralfördelningen av felet. I det sista steget uppnås den minsta möjliga noggrannheten och konvergensen stannar eller metoden kan till och med börja divergera. I typiska vetenskapliga datortillämpningar i flyttalformat med dubbel precision för matriser av stora storlekar, använder den konjugerade gradientmetoden ett stoppkriterie med en tolerans som avslutar iterationerna under det första eller andra steget.

Den förkonditionerade konjugerade gradientmetoden

Se även: Preconditioner

I de flesta fall är förkonditionering nödvändig för att säkerställa snabb konvergens av konjugatgradientmetoden. Om $\mathbf {M} ^{-1}$ är symmetrisk positiv-definitiv och $\mathbf {M} ^{-1}\mathbf {A}$ har en bättre villkorsnummer än $\mathbf {A}$ , kan en förkonditionerad konjugerad gradientmetod användas. Den har följande form:

$\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}$

$\mathbf { z} _{0}:=\mathbf {M} ^{-1}\mathbf {r} _{0}$

$\mathbf {p} _{0}:=\mathbf {z } _{0}$

$k:=0\,$

upprepa

$\alpha _{k}:={\frac { \mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}$

$\mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{ k}\mathbf {p} _{k}$

$\mathbf {r} _{k+1}:=\mathbf {r} _{k }-\alpha _{k}\mathbf {Ap} _{k}$

om r _{k +1} är tillräckligt liten, gå ur loopslut om

${\displaystyle \mathbf { z} _{k+1}:=\mathbf {M} ^{-1}\mathbf {r} _{k+1}} β$

$\beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k+1}}{\mathbf { r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}}$

${\displaystyle \mathbf {p} _ {k+1}:=\mathbf {z} _{k+1}+\beta _{k}\mathbf {p} _{k}} k := k + 1 {\displaystyle$

$1 \,}$

avsluta upprepning

Resultatet är x _{k +1}

Ovanstående formulering är ekvivalent med att tillämpa den vanliga konjugatgradientmetoden på det förkonditionerade systemet

$\mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\ mathsf {T}}\mathbf {\hat {x}} =\mathbf {E} ^{-1}\mathbf {b}$

var

$\mathbf {EE} ^{\mathsf {T}}=\mathbf {M} ,\qquad \mathbf {\hat {x}} =\mathbf {E} ^{\mathsf {T}}\mathbf {x} .$

Den Cholesky sönderdelningen av förkonditioneringsmedlet måste användas för att bibehålla systemets symmetri (och positiva bestämdhet). Denna nedbrytning behöver dock inte beräknas, och det räcker med att känna till $\mathbf {M} ^{-1}$ . Det kan visas att $\mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf { T}}$ har samma spektrum som $\mathbf {M} ^{-1}\mathbf {A}$ .

Förkonditioneringsmatrisen M måste vara symmetrisk positiv-definitiv och fixerad, dvs den kan inte ändras från iteration till iteration. Om något av dessa antaganden på förkonditioneringsmedlet kränks, kan beteendet hos den förkonditionerade konjugerade gradientmetoden bli oförutsägbart.

Ett exempel på en vanlig förkonditioneringsmedel är den ofullständiga Cholesky-faktoriseringen .

Den flexibla förkonditionerade konjugatgradientmetoden

I numeriskt utmanande applikationer används sofistikerade förkonditioneringsmedel, vilket kan leda till variabel förkonditionering, växlande mellan iterationerna. Även om förkonditioneraren är symmetrisk positiv-definitiv vid varje iteration, gör det faktum att den kan ändras argumenten ovan ogiltiga, och i praktiska tester leder det till en betydande nedgång av konvergensen av algoritmen som presenteras ovan. Använder Polak-Ribière -formeln

$\beta _{k}:={\frac {\mathbf {r} _{k+1}^ {\mathsf {T}}\left(\mathbf {z} _{k+1}-\mathbf {z} _{k}\right)}{\mathbf {r} _{k}^{\mathsf { T}}\mathbf {z} _{k}}}$

istället för Fletcher-Reeves formel

$\beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T }}\mathbf {z} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}}$

kan dramatiskt förbättra konvergensen i detta fall. Denna version av den förkonditionerade konjugatgradientmetoden kan kallas flexibel, eftersom den tillåter variabel förkonditionering. Den flexibla versionen har också visat sig vara robust även om förkonditioneraren inte är symmetrisk positiv definit (SPD).

Implementeringen av den flexibla versionen kräver lagring av en extra vektor. För en fast SPD-förkonditionering, $\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k}=0,$ så båda formlerna för β _k är ekvivalenta i exakt aritmetik, dvs utan avrundningsfelet .

Den matematiska förklaringen av metodens bättre konvergensbeteende med Polak–Ribière- formeln är att metoden är lokalt optimal i detta fall, i synnerhet konvergerar den inte långsammare än den lokalt optimala brantaste nedstigningsmetoden.

Mot. den lokalt optimala brantaste nedstigningsmetoden

I både den ursprungliga och den förkonditionerade konjugerade gradientmetoden behöver man bara ställa in $\beta _{k}:=0$ för att göra dem lokalt optimala, med hjälp av linjesökning , metoderna för brantast nedstigning . Med denna substitution är vektorer p alltid desamma som vektorer z , så det finns inget behov av att lagra vektorer p . Således är varje iteration av dessa brantaste nedstigningsmetoder lite billigare jämfört med den för konjugerade gradientmetoder. De senare konvergerar dock snabbare, såvida inte en (mycket) variabel och/eller icke-SPD förkonditionerare används, se ovan.

Konjugerad gradientmetod som optimal återkopplingskontroll för dubbelintegrator

Konjugatgradientmetoden kan också härledas med hjälp av optimal kontrollteori . I detta tillvägagångssätt faller den konjugerade gradientmetoden ut som en optimal återkopplingskontroller ,

$$u=k(x,v):=-\gamma _{a}\nabla f(x)-\ gamma _{b}v$$

för dubbelintegratorsystemet ,

$${\dot {x}}=v,\quad {\dot {v}}=u$$

Storheterna $\gamma _{a}$ och $\gamma _{b}$ är variabla återkopplingsförstärkningar.

Konjugera gradient på normalekvationerna

Den konjugerade gradientmetoden kan appliceras på en godtycklig n - by- m matris genom att applicera den på normala ekvationer A ^T A och höger sida vektor A ^T b , eftersom A ^T A är en symmetrisk positiv-semidefinit matris för vilken A som helst . Resultatet är konjugerad gradient på normalekvationerna (CGNR).

A ^T Axe = A ^T b

Som en iterativ metod är det inte nödvändigt att bilda A ^T A explicit i minnet utan endast att utföra matris-vektor och transponera matris-vektor multiplikationer. Därför är CGNR särskilt användbart när A är en gles matris eftersom dessa operationer vanligtvis är extremt effektiva. Men nackdelen med att bilda de normala ekvationerna är att villkorstalet κ( AT A ) är lika med κ ² ( A ) och därför kan konvergenshastigheten för CGNR vara långsam och kvaliteten på den ungefärliga lösningen kan vara känslig för ^avrundning fel. Att hitta en bra förkonditionering är ofta en viktig del av att använda CGNR-metoden.

Flera algoritmer har föreslagits (t.ex. CGLS, LSQR). LSQR- algoritmen påstås ha den bästa numeriska stabiliteten när A är dåligt konditionerad, dvs. A har ett stort villkorsnummer .

Konjugerad gradientmetod för komplexa hermitiska matriser

Den konjugerade gradientmetoden med en trivial modifiering kan utökas till att lösa, givet komplexvärderad matris A och vektor b, systemet av linjära ekvationer $\mathbf {A} \mathbf {x} =\mathbf {b }$ för vektorn x med komplext värde, där A är hermitisk (dvs. A' = A) och positiv-definitiv matris , och symbolen ' anger den konjugata transponeringen med stilen MATLAB / GNU Octave . Den triviala modifieringen är helt enkelt att ersätta den konjugerade transponeringen med den verkliga transponeringen överallt. Denna substitution är bakåtkompatibel, eftersom konjugerad transponering förvandlas till verklig transponering på vektorer och matriser med realt värde. Ovanstående exempelkod i MATLAB/GNU Octave fungerar alltså redan för komplexa hermitiska matriser som inte behövde modifieras.

Se även

• Biconjugate gradient method (BiCG)
• Konjugerad restmetod
• Gaussisk trosutbredning
• Iterativ metod: Linjära system
• Krylov underrum
• Icke-linjär konjugerad gradientmetod
• Förkonditionering
• Gles matris-vektor multiplikation

Vidare läsning

•   Atkinson, Kendell A. (1988). "Avsnitt 8.9". En introduktion till numerisk analys (2:a uppl.). John Wiley och söner. ISBN 978-0-471-50023-0 .
•   Avriel, Mordecai (2003). Icke-linjär programmering: analys och metoder . Dover Publishing. ISBN 978-0-486-43227-4 .
•   Golub, Gene H.; Van Loan, Charles F. (2013). "Kapitel 11". Matrix Computations (4:e upplagan). Johns Hopkins University Press. ISBN 978-1-4214-0794-4 .
•   Saad, Yousef (2003-04-01). "Kapitel 6" . Iterativa metoder för glesa linjära system (2:a uppl.). SIAM. ISBN 978-0-89871-534-7 .
• Gérard Meurant: "Detektion och korrigering av tysta fel i konjugatgradientalgoritmen", Numerical Algorithms, vol.92 (2023), s.869-891. url= https://doi.org/10.1007/s11075-022-01380-1

externa länkar

• "Conjugate gradients, method of" , Encyclopedia of Mathematics , EMS Press , 2001 [1994]