檢視廣義上細殘量方法的原始碼

佇咧數學上，'''廣義上細殘量方法'''( 一般簡稱'''GMRES''') 是一个求解線性方程組數值解的疊代方法。這个方法利用佇咧 Krylov 空間內底有著上細殘量的向量來逼近解。Arnoldi 疊代方法去予人用來求解這个向量。

GMRES 方法由 Yousef Saad 和 Martin H . Schultz 佇一九八六年提出。

==GMRES 方法==

需要求解的線性方程組記做


: $ Ax=b \ , $。

假設矩陣 _ A _ 是 _ n _ $ \ times $ _ n _ 階的倒反的。進一步，準講 _ b _ 是標準化的，即 | | _ b _ | |=一 ( 佇這篇文章內底，| | ・ | | 是 Euclidean 範數 )。

這个問題的 _ m _ 階 Krylov 子空間是


: $ K _ { m }=\ operatorname { span } \ , \ { b , Ab , A ^ { 二 } b , \ ldots , A ^ { m 影一 } b \ } \ , $。

GMRES 通過予伊殘量 _ Ax _ m − _ b _ 上細的向量 _ x _ m ∈ _ K _ m 來逼近 _ Ax _=_ b _ 的精確解。

猶毋過，向量 _ b _ , _ Ab _ ,…, _ A _ m− 一 _ b _ 差不多是線性相依的。所以，用 Arnoldi 疊代方法求得的這組 _ K _ m 的標準正交基


: $ q _ { 一 } , q _ { 二 } , \ ldots , q _ { m } \ , $

來取代頂懸彼組基。所以乎，向量 _ x _ m ∈ _ K _ m 寫做 _ x _ m=_ Q _ m _ y _ m，其中 _ y _ m ∈'''R'''m 而且 _ Q _ m 是由 _ q _ 一 ,…, _ q _ m 組成的 _ n _ $ \ times $ _ m _ 矩陣。

Arnoldi 過程嘛產生一个 ( _ m _ + 一 ) $ \ times $ _ m _ 階上 Hessenberg 矩陣 $ { \ tilde { H } } _ { m } $ 滿足


: $ AQ _ { m }=Q _ { m + 一 } { \ tilde { H } } _ { m } \ , $。

因為乎 $ Q _ { m } $ 是正交的，阮有


: $ \ | Ax _ { m }-b \ |=\ | { \ tilde { H } } _ { m } y _ { m }-\ beta e _ { 一 } \ | \ , $，

其中


: $ e _ { 一 }=( 一 , 零 , 零 , \ ldots , 零 ) \ , $

是'''R'''m + 一的標準基的第一个向量，並且


: $ \ beta=\ | b-Ax _ { 零 } \ | \ , $，

其中 $ x _ { 零 } $ 是初初向量 ( 通常是零向量 )。所以，求予得殘量


: $ r _ { m }={ \ tilde { H } } _ { m } y _ { m }-\ beta e _ { 一 } $。

的範數上細的 $ x _ { m } $。這是一个 _ m _ 階線性上細平方問題。

這就是 GMRES 方法。佇咧疊代的每一步內底：

一 . 做一步 Arnoldi 疊代方法；
二 . 走揣會使得 | | _ r _ m | | 上細漢的喔 $ y _ { m } $；
三 . 計算 $ x _ { m }=Q _ { m } y _ { m } $；
四 . 若無量無夠細，重複以上的過程。

佇每一步疊代中，著愛計算一改矩陣向量積 _ Aq _ m。對著一般的 _ n _ 階稠密矩陣，這愛計算複雜度差不多二 _ n _ 二改浮點運算。但是對疏疏矩陣，這个計算複雜度減到 O ( _ n _ )。進一步，關於矩陣向量積，佇咧 _ m _ 次疊代中會當進行 O ( _ m _ _ n _ ) 次浮點運算。

==收斂性==

第 _ m _ 次疊代得著佇 Krylov 子空間 _ K _ m 下的上細殘量。因為每一个子空間包含著後一个子空間內底，所以殘量單調遞減。佇咧第 _ n _ 次疊代後，其中 _ n _ 是矩陣 _ A _ 的階數，Krylov 空間 _ K _ n 是完整的'''R'''n。所以，GMRES 方法達到精確解。毋過，問題是因為：佇咧極少的幾若改疊代了後 ( 相對的是 _ n _ )，向量 _ x _ m 差不多已經是真精確解的一个真好的逼近。

猶毋過，佇一般情形下這是袂發生的。事實上，Greenbaum，Pták 和 Strakoš 的理論說明矣對每一个單調減少的序列 _ a _ 一 ,…, _ a _ n− 一 , _ a _ n=零，會當揣著一个矩陣 _ A _ 對所有 _ m _ 滿足 | | _ r _ m | |=_ a _ m，其中 _ r _ m 是頂頭所定義的殘量。特別的，有可能揣著一个矩陣，會當予𪜶前 _ n _ −  一改疊代的殘量一直保持做常數，毋過干焦佇最後一改疊代時到零。

佇實驗內底，GMRES 方法定定表現了真好。佇特殊的情形下這會使予人證明。若是 _ A _ 是正定的，著


: $ \ | r _ { m } \ | \ leq \ left ( 一-{ \ frac { \ lambda _ { \ mathrm { min } } ( A ^ { T } + A ) } { 二 \ lambda _ { \ mathrm { max } } ( A ^ { T } + A ) } } \ right ) ^ { m / 二 } \ | r _ { 零 } \ | $，

其中 $ \ lambda _ { \ mathrm { min } } ( M ) $ 和 $ \ lambda _ { \ mathrm { max } } ( M ) $ 分別為矩陣 $ M $ 的上小佮上大特徵值。

若是 _ A _ 是對稱的並且是正定的，著


: $ \ | r _ { m } \ | \ leq \ left ( { \ frac { \ kappa _ { 二 } ^ { 二 } ( A ) 影一 } { \ kappa _ { 二 } ^ { 二 } ( A ) } } \ right ) ^ { m / 二 } \ | r _ { 零 } \ | $。

其中 $ \ kappa _ { 二 } ( A ) $ 記為 _ A _ 佇咧 Euclidean 範數下跤的條件數。

一般情形下，其中 _ A _ 是無正定的，著


: $ \ | r _ { m } \ | \ leq \ inf _ { p \ in P _ { m } } \ | p _ { m } ( A ) \ | \ leq \ kappa _ { 二 } ( V ) \ inf _ { p \ in P _ { m } } \ max _ { \ lambda \ in \ sigma ( A ) } | p ( \ lambda ) | \ , $，

其中 _ P _ m 記為次數無超過 _ m _ 而且 _ p _ ( 零 )=一的多項式的集合，_ V _ 是 _ A _ 的譜分解中的矩陣，而且 σ ( _ A _ ) 是 _ A _ 的譜。粗略的講，當 _ A _ 的特徵值聚集佇遠離原點的區域而且 _ A _ 離正規無蓋遠的時，收縮速度較緊。

所有的無等式只界定殘量，毋是真正精差 ( 確定佮彼陣仔疊代 _ x _ m 之間的距離 )。

==GMRES 方法的拓展 ( Restarted GMRES )==

仝其他疊代方法仝款，為著加緊收斂，GMRES 定定結合預處理方法。

疊代的開銷以 O ( _ m _ 二 ) 增加長，其中 _ m _ 是疊代次數。毋過有當時仔，GMRES 方法佇咧 _ k _ 次疊代後重新開始，即 _ x _ k 閣變做轉來初初值。按呢的方法叫做 GMRES ( _ k _ )。

==佮其他解法的比較==

對稱矩陣來講，Arnoldi 疊代方法變成講 Lanczos 疊代方法。對應的 Krylov 子空間方法叫做 Paige 和 Saunders 的上細殘量方法 ( MinRes )。無像非對稱的情形，MinRes 方法由三項循環關係 ( three-term recurrence relation ) 給出，並且仝款 GMRES 仝款，殘量的範數上細。毋過對一般矩陣，Krylov 子空間方法袂當由短的循環關係 ( short recurrence relation ) 給出。

另外一類方法的由非對稱 Lanczos 疊代方法予出，特別的所在是 BiCG 方法。這利用矣 three-term recurrence relation，但是𪜶無達到上細的殘量，所以對遮的方法殘量袂單調遞減。斂性是袂當保證的。

第三類辦法由 CGS 和 BiCGSTAB 給出。這寡嘛由 three-term recurrence relation 給出 ( 所以，非上優 )。而且可能過早的終止疊代矣而且無達到收斂的目的。遮的方法的想法是合適的選擇共疊代序列所產生的多項式。

所有的對所有矩陣，這三類方法攏毋是上好的；總是有一類的方法予另外一類。因而，各種解法應該愛做實際的試驗，來決定對於予定的問題佗一種是上優的。

==求解上細漢平常時問題==

GMRES 方法其中一部份是求解向量 $ y _ { m } $ 予得


: $ \ | { \ tilde { H } } _ { m } y _ { m }-e _ { 一 } \ | \ , $

上細漢。這會當通過計算 QR 分解來實現：揣著一个 ( _ m _ + 一 ) $ \ times $ ( _ m _ + 一 ) 階正交矩陣 Ωm 佮一个 ( _ m _ + 一 ) $ \ times $ _ m _ 上三角矩陣 $ { \ tilde { R } } _ { m } $ 滿足


: $ \ Omega _ { m } { \ tilde { H } } _ { m }={ \ tilde { R } } _ { m } $。

三角矩陣的行數比列加一，所以伊的上尾一逝由零組成。所以，伊會當分解掉去


: $ { \ tilde { R } } _ { m }={ \ begin { bmatrix } R _ { m } \ \ 零 \ end { bmatrix } } $，

其中 $ R _ { m } $ 是一个 _ m _ $ \ times $ _ m _ 階三角 ( 方 ) 矩陣。

QR 分解會當簡單的進行落去 ( update )，對一步來疊代到後一步疊代。因為逐改 Hessenberg 矩陣只要一行零元素佮一列元素頂懸有所無仝：


: $ { \ tilde { H } } _ { m + 一 }={ \ begin { bmatrix } { \ tilde { H } } _ { m } & h _ { m } \ \ 零 & h _ { m + 一 , m } \ end { bmatrix } } $，

其中 _ h _ m=( _ h _ 一 _ m _ ,…_ h _ mm ) T。這意味對，Hessenberg 矩陣倒乘 Ωm 的擴大矩陣 ( 通過閣上零元素佮單位元素 )，所得的是類似三角矩陣的矩陣：


: $ { \ begin { bmatrix } \ Omega _ { m } & 零 \ \ 零 & 一 \ end { bmatrix } } { \ tilde { H } } _ { m + 一 }={ \ begin { bmatrix } R _ { m } & r _ { k } \ \ 零 & \ rho \ \ 零 & \ sigma \ end { bmatrix } } $

這个矩陣會當三角化，若是 σ 為零。為著欲修正這个矩陣，愛需要做 Givens 旋轉


: $ G _ { m }={ \ begin { bmatrix } I _ { m 影一 } & 零 & 零 \ \ 零 & c _ { m } & s _ { m } \ \ 零 &-s _ { m } & c _ { m } \ end { bmatrix } } $

其中


: $ c _ { m }={ \ frac { \ rho } { \ sqrt { \ rho ^ { 二 } + \ sigma ^ { 二 } } } } \ quad { \ mbox { and } } \ quad s _ { m }={ \ frac { \ sigma } { \ sqrt { \ rho ^ { 二 } + \ sigma ^ { 二 } } } } $。

通過這 Givens 旋轉，所以阮構造


: $ \ Omega _ { m + 一 }=G _ { m } { \ begin { bmatrix } \ Omega _ { m } & 零 \ \ 零 & 一 \ end { bmatrix } } $。

事實上，


: $ \ Omega _ { m + 一 } { \ tilde { H } } _ { m + 一 }={ \ begin { bmatrix } R _ { m } & r _ { m } \ \ 零 & r _ { mm } \ \ 零 & 零 \ end { bmatrix } } \ quad { \ text { 其中 } } \ quad r _ { mm }={ \ sqrt { \ rho ^ { 二 } + \ sigma ^ { 二 } } } $

是一个三角矩陣。

給出了 QR 分解，上細漢的問題就較𠢕解決矣。注意著


: $ \ | { \ tilde { H } } _ { m } y _ { m }-e _ { 一 } \ |=\ | \ Omega _ { m } ( { \ tilde { H } } _ { m } y _ { m }-e _ { 一 } ) \ |=\ | { \ tilde { R } } _ { m } y _ { m }-\ Omega _ { m } e _ { 一 } \ | $。

記 $ \ Omega _ { m } e _ { 一 } $ 為


: $ { \ tilde { g } } _ { m }={ \ begin { bmatrix } g _ { m } \ \ \ gamma _ { m } \ end { bmatrix } } $

其中 _ g _ m ∈'''R'''m 和 γm ∈'''R'''，著


: $ \ | { \ tilde { H } } _ { m } y _ { m }-e _ { 一 } \ |=\ | { \ tilde { R } } _ { m } y _ { m }-\ Omega _ { m } e _ { 一 } \ |=\ left \ | { \ begin { bmatrix } R _ { m } \ \ 零 \ end { bmatrix } } y-{ \ begin { bmatrix } g _ { m } \ \ \ gamma _ { m } \ end { bmatrix } } \ right \ | $。

會當予這个表達式上細的向量 _ y _ 為


: $ y _ { m }=R _ { m } ^ { 影一 } g _ { m } $。

閣一改，向量 $ g _ { m } $ 會當簡單的來進行落去喔 ( update )。

==Example code==

===Regular GMRES ( python 三 )===

==註記==

==參考==

* A . Meister , _ Numerik linearer Gleichungssysteme _ , 二 nd edition , Vieweg 兩千空五 , ISBN 九百七十八追三分五百二十八分一鋪三千一百三十五刣七 .
* Y . Saad , _ Iterative Methods for Sparse Linear Systems _ , 二 nd edition , Society for Industrial and Applied Mathematics , 兩千空三 . ISBN 九百七十八追空抹八九千八百七十一孵五百三十四刣七 .
* Y . Saad and M . H . Schultz , " GMRES : A generalized minimal residual algorithm for solving nonsymmetric linear systems " , _ SIAM J . Sci . Stat . Comput . _ ,'''七''': 八百五十六交八百六十九 , 一千九百八十六 . doi : 十 . 九十五空七千空五十八分之一千一百三十七 .
* J . Stoer and R . Bulirsch , _ Introduction to numerical analysis _ , 三 rd edition , Springer , New York , 兩千空二 . ISBN 九百七十八追空九三百八十七石得九九五千四百五十二孵三 .
* Lloyd N . Trefethen and David Bau , III , _ Numerical Linear Algebra _ , Society for Industrial and Applied Mathematics , 一千九百九十七 . ISBN 九百七十八追空抹八七九千八百七十一孵三百六十一孵九 .
* Dongarra et al . , _ Templates for the Solution of Linear Systems : Building Blocks for Iterative Methods _ , 二 nd Edition , SIAM , Philadelphia , 一千九百九十四
* https : / / github . com / J-N-ch / GMRES \ _ py \ _ restart

[[分類: 待校正]]