Konjugovaná metoda přechodu - Conjugate gradient method

Srovnání konvergence klesání s optimální velikostí kroku (zeleně) a konjugovaným vektorem (červeně) pro minimalizaci kvadratické funkce spojené s daným lineárním systémem. Konjugovaný gradient za předpokladu přesné aritmetiky konverguje maximálně n kroky, kde n je velikost matice systému (zde n = 2).

v matematika, metoda konjugovaného gradientu je algoritmus pro numerické řešení zvláště soustavy lineárních rovnic, jmenovitě ti, jejichž matice je symetrický a pozitivní-definitivní. Metoda sdruženého přechodu je často implementována jako iterační algoritmus, použitelné pro řídký systémy, které jsou příliš velké na to, aby je bylo možné zpracovat přímou implementací nebo jinými přímými metodami, jako je Choleský rozklad. Při numerickém řešení často vznikají velké řídké systémy parciální diferenciální rovnice nebo problémy s optimalizací.

Metodu konjugovaného přechodu lze také použít k řešení neomezeného optimalizace problémy jako minimalizace energie. Byl vyvinut hlavně společností Magnus Hestenes a Eduard Stiefel,^[1][2] kdo to naprogramoval na Z4.^[3]

The biconjugate gradientní metoda poskytuje zobecnění nesymetrických matic. Rozličný metody nelineárního konjugovaného gradientu hledat minima nelineárních rovnic.

Popis problému řešeného přechody konjugátu

Předpokládejme, že chceme vyřešit soustava lineárních rovnic

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

pro vektor X, kde je známo n × n matice A je symetrický (tj., A^T = A), pozitivní-definitivní (tj. X^TSekera > 0 pro všechny nenulové vektory X v Rⁿ), a nemovitý, a b je také známý. Označujeme jedinečné řešení tohoto systému ${ displaystyle mathbf {x} _ {*}}$.

Jako přímá metoda

Říkáme, že dva nenulové vektory u a proti jsou sdružené (s ohledem na A) pokud

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

Od té doby A je symetrický a kladně definitivní, levá strana definuje znak vnitřní produkt

${ displaystyle 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.}$

Dva vektory jsou konjugovány právě tehdy, pokud jsou kolmé vzhledem k tomuto vnitřnímu produktu. Být konjugátem je symetrický vztah: pokud u je konjugován s proti, pak proti je konjugován s u. Předpokládejme to

${ displaystyle P = { mathbf {p} _ {1}, tečky, mathbf {p} _ {n} }}$

je sada n vzájemně konjugované vektory (s ohledem na A). Pak P tvoří a základ pro ${ displaystyle mathbb {R} ^ {n}}$a můžeme vyjádřit řešení X_∗ z ${ displaystyle mathbf {Axe} = mathbf {b}}$ na tomto základě:

${ displaystyle mathbf {x} _ {*} = součet _ {i = 1} ^ {n} alpha _ {i} mathbf {p} _ {i}.}$

Na základě této expanze vypočítáme:

${ displaystyle mathbf {A} mathbf {x} _ {*} = součet _ {i = 1} ^ {n} alpha _ {i} mathbf {A} mathbf {p} _ {i} .}$

Násobení doleva ${ displaystyle mathbf {p} _ {k} ^ { mathsf {T}}}$:

${ displaystyle mathbf {p} _ {k} ^ { mathsf {T}} mathbf {A} mathbf {x} _ {*} = součet _ {i = 1} ^ {n} alpha _ {i} mathbf {p} _ {k} ^ { mathsf {T}} mathbf {A} mathbf {p} _ {i},}$

střídání ${ displaystyle mathbf {Ax _ {*}} = mathbf {b}}$ a ${ displaystyle mathbf {u} ^ { mathsf {T}} mathbf {A} mathbf {v} = langle mathbf {u}, mathbf {v} rangle _ { mathbf {A}} }$:

${ displaystyle mathbf {p} _ {k} ^ { mathsf {T}} mathbf {b} = součet _ {i = 1} ^ {n} alpha _ {i} leva langle mathbf {p} _ {k}, mathbf {p} _ {i} right rangle _ { mathbf {A}},}$

pak ${ displaystyle mathbf {u} ^ { mathsf {T}} mathbf {v} = langle mathbf {u}, mathbf {v} rangle}$ a pomocí ${ displaystyle forall i neq k: langle mathbf {p} _ {k}, mathbf {p} _ {i} rangle _ { mathbf {A}} = 0}$ výnosy

${ displaystyle langle mathbf {p} _ {k}, mathbf {b} rangle = alpha _ {k} langle mathbf {p} _ {k}, mathbf {p} _ {k} rangle _ { mathbf {A}},}$

z čehož vyplývá

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

To poskytuje následující metodu řešení rovnice Sekera = b: najít posloupnost n konjugujte směry a poté vypočítejte koeficienty α_k.

Jako iterativní metoda

Pokud zvolíme konjugované vektory p_k pečlivě pak možná nebudeme potřebovat všechny, abychom získali dobrou aproximaci řešení X_∗. Chceme tedy považovat metodu konjugovaného gradientu za iterační metodu. To nám také umožňuje přibližně vyřešit systémy, kde n je tak velká, že přímá metoda by zabrala příliš mnoho času.

Označíme počáteční odhad pro X_∗ podle X₀ (můžeme to předpokládat bez ztráty obecnosti X₀ = 0, jinak zvažte systém Az = b − Sekera₀ namísto). Začínání s X₀ hledáme řešení a v každé iteraci potřebujeme metriku, která nám řekne, zda jsme blíže řešení X_∗ (to je nám neznámé). Tato metrika vychází ze skutečnosti, že řešení X_∗ je také jedinečným minimalizátorem následujícího kvadratická funkce

${ displaystyle 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} ,.}$

Existence jedinečného minimalizátoru je zjevná, protože jeho druhá derivace je dána symetrickou pozitivně-definitivní maticí

${ displaystyle nabla ^ {2} f ( mathbf {x}) = mathbf {A} ,,}$

a že minimalizátor (použijte DF(X) = 0) řeší počáteční problém, který je zřejmý z jeho první derivace

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

To naznačuje použití prvního základního vektoru p₀ být záporem gradientu F v X = X₀. Přechod z F rovná se Sekera − b. Počínaje počátečním odhadem X₀, to znamená, že bereme p₀ = b − Sekera₀. Ostatní vektory v bázi budou konjugovány s přechodem, odtud název metoda konjugovaného gradientu. Všimněte si, že p₀ je také reziduální poskytované tímto počátečním krokem algoritmu.

Nechat r_k být reziduální na kt krok:

${ displaystyle mathbf {r} _ {k} = mathbf {b} - mathbf {Axe} _ {k}.}$

Jak bylo uvedeno výše, r_k je záporný gradient F v X = X_k, takže klesání metoda by vyžadovala pohyb ve směru r_k. Zde však trváme na tom, že pokyny p_k být navzájem konjugovány. Praktickým způsobem, jak to vynutit, je požadavek, aby byl další směr hledání vytvořen z aktuálního zbytkového a všech předchozích směrů vyhledávání.^[4] To dává následující výraz:

${ displaystyle mathbf {p} _ {k} = mathbf {r} _ {k} - součet _ {i$

(vliv omezení konjugace na konvergenci najdete na obrázku v horní části článku). Po tomto směru je další optimální umístění dáno vztahem

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

s

${ displaystyle 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}}}, }$

kde poslední rovnost vyplývá z definice r_k Výraz pro ${ displaystyle alpha _ {k}}$ lze odvodit, pokud jeden nahradí výraz X_k+1 do F a minimalizovat to w.r.t. ${ displaystyle alpha _ {k}}$

${ displaystyle { 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 {zarovnáno}}}$

Výsledný algoritmus

Výše uvedený algoritmus poskytuje nejpřímější vysvětlení metody konjugovaného gradientu. Zdánlivě algoritmus, jak je uvedeno, vyžaduje uložení všech předchozích směrů hledání a vektorů zbytků, stejně jako mnoho násobení matice-vektor, a proto může být výpočetně nákladný. Podrobnější analýza algoritmu to však ukazuje r_i je kolmý na r_j , tj. ${ displaystyle mathbf {r} _ {i} ^ { mathsf {T}} mathbf {r} _ {j} = 0}$ , protože i ≠ j. A p_i je A-ortogonální k p_j , tj. ${ displaystyle mathbf {p} _ {i} ^ { mathsf {T}} A mathbf {p} _ {j} = 0}$ , protože i ≠ j. To lze považovat za to, že jak algoritmus postupuje, p_i a r_i rozpětí stejné Krylovský podprostor. Kde r_i - tvoří ortogonální základ s ohledem na standardní vnitřní produkt a - p_i tvoří ortogonální základ s ohledem na vnitřní produkt indukovaný A. Proto X_k lze považovat za projekci X na Krylovském podprostoru.

Algoritmus je pro řešení podrobně popsán níže Sekera = b kde A je skutečná, symetrická, pozitivně určitá matice. Vstupní vektor X₀ může být přibližné počáteční řešení nebo 0. Jedná se o jinou formulaci přesného postupu popsaného výše.

${ displaystyle { begin {aligned} & mathbf {r} _ {0}: = mathbf {b} - mathbf {Axe} _ {0} & { hbox {if}} mathbf {r } _ {0} { text {je dostatečně malý, pak se vrátí}} mathbf {x} _ {0} { text {jako výsledek}} & 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 {je dostatečně malý, pak ukončete smyčku}} & 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 {end repeat}} & { text {return}} mathbf {x} _ {k + 1} { text {jako výsledek}} end {zarovnáno}}}$

Toto je nejčastěji používaný algoritmus. Stejný vzorec pro β_k se také používá u Fletcher – Reevese metoda nelineárního konjugovaného gradientu.

Výpočet alfa a beta

V algoritmu α_k je zvolen tak, že ${ displaystyle mathbf {r} _ {k + 1}}$ je kolmý na r_k. Jmenovatel je zjednodušen od

${ displaystyle 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}}}}$

od té doby ${ displaystyle mathbf {r} _ {k + 1} = mathbf {p} _ {k + 1} - mathbf { beta} _ {k} mathbf {p} _ {k}}$. The β_k je zvolen tak, že ${ displaystyle mathbf {p} _ {k + 1}}$ je konjugován s p_k. Zpočátku, β_k je

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

použitím

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

a rovnocenně

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

čitatel β_k je přepsán jako

${ displaystyle 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}}$

protože ${ displaystyle mathbf {r} _ {k + 1}}$ a r_k jsou kolmé záměrně. Jmenovatel je přepsán jako

${ 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}}$

pomocí toho vyhledejte pokyny p_k jsou konjugovány a opět, že zbytky jsou ortogonální. To dává β v algoritmu po zrušení α_k.

Příklad kódu v MATLAB / GNU oktáva

funkceX =Congrad(A, b, x)r = b - A * X;    p = r;    rsold = r' * r;    pro i = 1: délka (b)        Ap = A * p;        alfa = rsold / (p' * Ap);        X = X + alfa * p;        r = r - alfa * Ap;        rsnew = r' * r;        -li sqrt (rsnew) <1e-10              přestávkakonec        p = r + (rsnew / rsold) * p;        rsold = rsnew;    koneckonec

Numerický příklad

Zvažte lineární systém Sekera = b dána

${ displaystyle 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}},}$

provedeme dva kroky metody konjugovaného přechodu počínaje počátečním odhadem

${ displaystyle mathbf {x} _ {0} = { začátek {bmatrix} 2 1 konec {bmatrix}}}$

za účelem nalezení přibližného řešení systému.

Řešení

Přesné řešení je

${ displaystyle mathbf {x} = { begin {bmatrix} { frac {1} {11}} { frac {7} {11}} end {bmatrix}} cca { začátek {bmatrix} 0,0909 0,6364 end {bmatrix}}}$

Naším prvním krokem je výpočet zbytkového vektoru r₀ spojený s X₀. Tento zbytek se počítá ze vzorce r₀ = b - Sekera₀, a v našem případě se rovná

${ displaystyle 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}.}$

Jelikož se jedná o první iteraci, použijeme zbytkový vektor r₀ jako náš počáteční směr hledání p₀; způsob výběru p_k se v dalších iteracích změní.

Nyní vypočítáme skalár α₀ pomocí vztahu

${ displaystyle 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}}.}$

Nyní můžeme počítat X₁ pomocí vzorce

${ displaystyle 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}} = { begin {bmatrix} 0,2356 0,3384 end {bmatrix} }.}$

Tento výsledek dokončuje první iteraci, výsledkem je „vylepšené“ přibližné řešení systému, X₁. Nyní můžeme pokračovat a vypočítat další zbytkový vektor r₁ pomocí vzorce

${ displaystyle 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}} = { begin {bmatrix} -0,2810 0,7492 end {bmatrix}}.}$

Naším dalším krokem v tomto procesu je výpočet skaláru β₀ který bude nakonec použit k určení dalšího směru hledání p₁.

${ displaystyle beta _ {0} = { frac { mathbf {r} _ {1} ^ { mathsf {T}} mathbf {r} _ {1}} { mathbf {r} _ {0 } ^ { mathsf {T}} mathbf {r} _ {0}}} = { frac {{ begin {bmatrix} -0,2810 & 0,7492 end {bmatrix}} { begin {bmatrix} -0,2810 0,7492 end {bmatrix}}}} {{ start {bmatrix} -8 & -3 end {bmatrix}} { begin {bmatrix} -8 - 3 end {bmatrix}}}} = 0,0088. }$

Nyní pomocí tohoto skaláru β₀, můžeme vypočítat další směr hledání p₁ pomocí vztahu

${ displaystyle mathbf {p} _ {1} = mathbf {r} _ {1} + beta _ {0} mathbf {p} _ {0} = { begin {bmatrix} -0,2810 0,7492 end {bmatrix}} + 0,0088 { begin {bmatrix} -8 - 3 end {bmatrix}} = { begin {bmatrix} -0,3511 0,7229 end {bmatrix}}.}$

Nyní vypočítáme skalár α₁ pomocí našich nově získaných p₁ pomocí stejné metody, jaká byla použita pro α₀.

${ displaystyle alpha _ {1} = { frac { mathbf {r} _ {1} ^ { mathsf {T}} mathbf {r} _ {1}} { mathbf {p} _ {1 } ^ { mathsf {T}} mathbf {Ap} _ {1}}} = { frac {{ begin {bmatrix} -0,2810 & 0,7492 end {bmatrix}} { begin {bmatrix} -0,2810 0,7492 end {bmatrix}}} {{ start {bmatrix} -0,3511 a 0,7229 end {bmatrix}} { begin {bmatrix} 4 a 1 1 & 3 end {bmatrix}} { begin {bmatrix} -0,3511 0,7229 end {bmatrix}}}} = 0,4122.}$

Nakonec zjistíme X₂ pomocí stejné metody jako ta, která byla použita k vyhledání X₁.

${ displaystyle mathbf {x} _ {2} = mathbf {x} _ {1} + alpha _ {1} mathbf {p} _ {1} = { 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}}.}$

Výsledek, X₂, je „lepším“ přiblížením k řešení systému než X₁ a X₀. Pokud by se v tomto příkladu měla použít přesná aritmetika namísto omezené přesnosti, pak by bylo přesného řešení teoreticky dosaženo po n = 2 iterace (n pořadí systému).

Vlastnosti konvergence

Metodu konjugovaného gradientu lze teoreticky považovat za přímou metodu, protože vytváří přesné řešení po konečném počtu iterací, který není větší než velikost matice, při absenci chyba zaokrouhlení. Metoda konjugovaného gradientu je však nestabilní s ohledem na i malé odchylky, např. Většina směrů není v praxi konjugovaná a přesné řešení se nikdy nezíská. Naštěstí lze metodu konjugovaného gradientu použít jako iterační metoda protože poskytuje monotónně zlepšující aproximace ${ displaystyle mathbf {x} _ {k}}$ k přesnému řešení, které může dosáhnout požadované tolerance po relativně malém počtu iterací (ve srovnání s velikostí problému). Zlepšení je obvykle lineární a jeho rychlost je určena číslo podmínky ${ displaystyle kappa (A)}$ matice systému ${ displaystyle A}$: větší ${ displaystyle kappa (A)}$ je, tím pomalejší je zlepšení.^[5]

Li ${ displaystyle kappa (A)}$ je velký, předběžná úprava se používá k nahrazení původního systému ${ displaystyle mathbf {Ax} - mathbf {b} = 0}$ s ${ displaystyle mathbf {M} ^ {- 1} ( mathbf {Axe} - mathbf {b}) = 0}$ takhle ${ displaystyle kappa ( mathbf {M} ^ {- 1} mathbf {A})}$ je menší než ${ displaystyle kappa ( mathbf {A})}$, viz. níže.

Věta o konvergenci

Definujte podmnožinu polynomů jako

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

kde ${ displaystyle Pi _ {k}}$ je sada polynomy maximálního stupně ${ displaystyle k}$.

Nechat ${ displaystyle left ( mathbf {x} _ {k} right) _ {k}}$ být iterativní aproximace přesného řešení ${ displaystyle mathbf {x} _ {*}}$a definovat chyby jako ${ displaystyle mathbf {e} _ {k}: = mathbf {x} _ {k} - mathbf {x} _ {*}}$Nyní lze míru konvergence přiblížit jako ^[6]

${ displaystyle { 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}} vpravo) ^ {k} vlevo | mathbf {e} _ {0} vpravo | _ { mathbf {A}} ,, end {zarovnáno}}}$

kde ${ displaystyle sigma ( mathbf {A})}$ označuje spektrum, a ${ displaystyle kappa ( mathbf {A})}$ označuje číslo podmínky.

Důležitý limit, když ${ displaystyle kappa ( mathbf {A})}$ má sklony k ${ displaystyle infty}$

${ displaystyle { frac {{ sqrt { kappa ( mathbf {A})}} - 1} {{ sqrt { kappa ( mathbf {A})}} + 1}} přibližně 1- { frac {2} { sqrt { kappa ( mathbf {A})}}} quad { text {pro}} quad kappa ( mathbf {A}) gg 1 ,.}$

Tento limit ukazuje rychlejší míru konvergence ve srovnání s iterativními metodami Jacobi nebo Gauss – Seidel které měřítko jako ${ displaystyle cca 1 - { frac {2} { kappa ( mathbf {A})}}}$.

Předpřipravená metoda konjugovaného gradientu

Většinou, předběžná úprava je nutné k zajištění rychlé konvergence metody konjugovaného gradientu. Předpřipravená metoda konjugovaného přechodu má následující podobu:^[7]

${ displaystyle mathbf {r} _ {0}: = mathbf {b} - mathbf {Axe} _ {0}}$

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

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

${ displaystyle k: = 0 ,}$

opakovat

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

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

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

-li r_k+1 je dostatečně malý pak výstupní smyčka skončit, pokud

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

${ displaystyle 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}}$

${ displaystyle k: = k + 1 ,}$

konec opakování

Výsledek je X_k+1

Výše uvedená formulace je ekvivalentní aplikaci metody konjugovaného gradientu bez předběžné úpravy systému^[1]

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

kde

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

Matice preconditioneru M musí být symetrický kladně definitivní a pevný, tj. nemůže se měnit z iterace na iteraci. Pokud dojde k porušení některého z těchto předpokladů na předkondicionéru, může se chování metody předpřipraveného konjugovaného přechodu stát nepředvídatelným.

Příklad běžně používaného kondicionér je neúplná Choleského faktorizace.^[8]

Metoda flexibilního předpřipraveného konjugovaného gradientu

V numericky náročných aplikacích se používají sofistikované předkondicionéry, které mohou vést k variabilním předběžným úpravám, které se mezi iteracemi mění. I když je preconditioner symetrický kladně definitivní při každé iteraci, skutečnost, že se může změnit, činí výše uvedené argumenty neplatnými a v praktických testech vede k výraznému zpomalení konvergence výše uvedeného algoritmu. Za použití Polak – Ribière vzorec

${ displaystyle 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}}}}$

místo Fletcher – Reeves vzorec

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

může v tomto případě výrazně zlepšit konvergenci.^[9] Tuto verzi předem připravené metody konjugovaného přechodu lze volat^[10] flexibilní, protože umožňuje variabilní předběžné úpravy. Je zobrazena také flexibilní verze^[11] být robustní, i když předkondicionér není symetrický pozitivní určitý (SPD).

Implementace flexibilní verze vyžaduje uložení dalšího vektoru. Pro pevný předkondicionér SPD ${ displaystyle mathbf {r} _ {k + 1} ^ { mathsf {T}} mathbf {z} _ {k} = 0,}$ takže oba vzorce pro β_k jsou ekvivalentní v přesné aritmetice, tj. bez chyba zaokrouhlení.

Matematické vysvětlení lepšího konvergenčního chování metody s Polak – Ribière vzorec je, že metoda je místně optimální v tomto případě zejména nekonverguje pomaleji než lokálně optimální metoda nejstrmějšího sestupu.^[12]

Příklad kódu v MATLAB / GNU Octave

funkce[x, k] =cgp(x0, A, C, b, mit, stol, bbA, bbC)% Synopse:% x0: počáteční bod% A: Matice A systému Ax = b% C: Předběžná matice může být vlevo nebo vpravo% mit: Maximální počet iterací% stol: tolerance normy zbytku% bbA: Black Box, který počítá produkt matice-vektor pro A * u% bbC: Black Box, který počítá:% pro levostranný kondicionér: ha = C  ra% pro pravý kondicionér: ha = C * ra% x: Odhadovaný bod řešení% k: Počet provedených iterací %% Příklad:% tic; [x, t] = cgp (x0, S, speye (1), b, 3000, 10 ^ -8, @ (Z, o) Z * o, @ (Z, o) o); toc% Uplynulý čas je 0,550190 sekund.%% Reference:% Métodos iterativos tipo Krylov para sistema lineales% B. Molina y M. Raydan - {{ISBN | 908-261-078-X}}        -li nargin <8, chyba ('Nedostatek vstupních argumentů. Zkuste pomoc.'); konec;        -li isempty (A), error ('Input matrix A must not be empty.'); konec;        -li isempty (C), error ('Input preconditioner matrix C must not be empty.'); konec;        X = x0;        ha = 0;        hp = 0;        hpp = 0;        ra = 0;        rp = 0;        rpp = 0;        u = 0;        k = 0;        ra = b - bbA(A, x0); % <--- ra = b - A * x0;        zatímco norma (ra, inf)> stol                ha = bbC(C, ra); % <--- ha = C  ra;                k = k + 1;                -li (k == mit), Varování('GCP: MAXIT', „bylo dosaženo, žádná konverze.“); vrátit se; konec;                hpp = hp;                rpp = rp;                hp = ha;                rp = ra;                t = rp' * hp;                -li k == 1                        u = hp;                jinýu = hp + (t / (rpp '* hpp)) * u;                konec;                Au = bbA (A, u); % <--- Au = A * u;                a = t / (u '* Au);                X = X + A * u;                ra = rp - A * Au;        konec;

Vs. lokálně optimální metoda nejstrmějšího klesání

V původní i předpřipravené konjugované gradientní metodě je potřeba pouze nastavit ${ displaystyle beta _ {k}: = 0}$ aby byly lokálně optimální pomocí vyhledávání linek, nejstrmější sestup metody. S touto substitucí vektory p jsou vždy stejné jako vektory z, takže není třeba ukládat vektory p. Tedy každá jejich iterace nejstrmější sestup metody je o něco levnější ve srovnání s metodami s konjugovaným gradientem. Ty však konvergují rychleji, pokud nejde o (vysoce) proměnnou a / nebo non-SPD kondicionér se používá, viz výše.

Odvození metody

Metodu konjugovaného gradientu lze odvodit z několika různých perspektiv, včetně specializace metody směru konjugátu pro optimalizaci a variace Arnoldi /Lanczos iterace pro vlastní číslo problémy. Navzdory rozdílům v jejich přístupech sdílejí tyto derivace společné téma - což dokazuje ortogonalitu zbytků a konjugaci směrů hledání. Tyto dvě vlastnosti jsou zásadní pro vývoj dobře známé stručné formulace metody.

Metodu konjugovaného gradientu lze také odvodit pomocí teorie optimální regulace.^[13] V tomto přístupu metoda konjugovaného přechodu vypadne jako optimální zpětnovazební regulátor,

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

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

${ displaystyle gamma _ {a}}$

${ displaystyle gamma _ {b}}$

Konjugujte gradient na normálních rovnicích

Metodu sdruženého přechodu lze použít na libovolné n-podle-m matice jeho aplikací na normální rovnice A^TA a pravá strana vektor A^Tb, od té doby A^TA je symetrický kladně semidefinitní matice pro všechny A. Výsledkem je konjugovaný gradient na normálních rovnicích (CGNR).

A^TSekera = A^Tb

Jako iterační metodu není nutné vytvářet A^TA explicitně v paměti, ale pouze k provádění násobení matice-vektor a transpozice matice-vektor. Proto je CGNR zvláště užitečný, když A je řídká matice protože tyto operace jsou obvykle extrémně efektivní. Nevýhodou formování normálních rovnic je však to, že číslo podmínky κ (A^TA) se rovná κ²(A), a tak rychlost konvergence CGNR může být pomalá a kvalita přibližného řešení může být citlivá na chyby zaokrouhlování. Nalezení dobrého kondicionér je často důležitou součástí používání metody CGNR.

Bylo navrženo několik algoritmů (např. CGLS, LSQR). Algoritmus LSQR má údajně nejlepší numerickou stabilitu, když A je špatně podmíněný, tj. A má velký číslo podmínky.

Viz také

Reference

Další čtení

• Atkinson, Kendell A. (1988). „Oddíl 8.9“. Úvod do numerické analýzy (2. vyd.). John Wiley and Sons. ISBN  978-0-471-50023-0.
• Avriel, Mordechaj (2003). Nelineární programování: Analýza a metody. Dover Publishing. ISBN  978-0-486-43227-4.
• Golub, Gene H .; Van Loan, Charles F. (1996-10-15). „Kapitola 10“. Maticové výpočty (3. vyd.). Johns Hopkins University Press. ISBN  978-0-8018-5414-9.
• Saad, Yousef (01.04.2003). "Kapitola 6". Iterační metody pro řídké lineární systémy (2. vyd.). SIAM. ISBN  978-0-89871-534-7.

externí odkazy

• "Konjugované přechody, metoda", Encyclopedia of Mathematics, Stiskněte EMS, 2001 [1994]