GMRES法

数学において、GMRES法（GMRESほう、generalized minimal residual method）は、連立一次方程式の数値解を求めるための反復法の一種である[1]。残差をクリロフ部分空間において最小化することにより、近似解を計算する。ベクトルの計算にはアーノルディ法[2]が用いられる。ヨセフ・サードとマルティン・H・シュルツにより、1986年に開発された[3]。

解法

解くべき方程式を

A x = b . {\displaystyle Ax=b.\,}

とする。ただし行列Aは正則なm次正方行列、bは正規化された（ユークリッドノルム||·||に関して||b|| = 1となる）ベクトルとする。

この問題のn次クリロフ部分空間は

K n = span { b , A b , A 2 b , … , A n − 1 b } . {\displaystyle K_{n}=\operatorname {span} \,\{b,Ab,A^{2}b,\ldots ,A^{n-1}b\}.\,}

である。GMRES法では、残差Axn − bを最小化するベクトルxn ∈ Knによって、Ax = bの厳密解を近似する。

b, Ab, …, An−1bは線形従属に近いため、これを用いる代わりに、アーノルディ法を用いてKnの基底を構成する正規直交化ベクトル列

q 1 , q 2 , … , q n {\displaystyle q_{1},q_{2},\ldots ,q_{n}\,}

を計算する。すなわち、xn ∈ Knはyn ∈ Rnを用いてxn = Qnynと書くことができる。ただしQnはq1, …, qnにより構成されるm×n行列とする。

アーノルディ過程では、

A Q n = Q n + 1 H ~ n {\displaystyle AQ_{n}=Q_{n+1}{\tilde {H}}_{n}}

を満たす(n+1)×n次上ヘッセンベルグ行列 H ~ n {\displaystyle {\tilde {H}}_{n}} が生成される。 Q n {\displaystyle Q_{n}} は直交行列なので、

e 1 = ( 1 , 0 , 0 , … , 0 ) {\displaystyle e_{1}=(1,0,0,\ldots ,0)\,}

を標準基底Rn+1の第1ベクトルとして、

‖ A x n − b ‖ = ‖ H ~ n y n − e 1 ‖ {\displaystyle \|Ax_{n}-b\|=\|{\tilde {H}}_{n}y_{n}-e_{1}\|}

が成り立つ。したがって、 x n {\displaystyle x_{n}} は残差

r n = H ~ n y n − e 1 . {\displaystyle r_{n}={\tilde {H}}_{n}y_{n}-e_{1}.}

のノルムを最小化することにより計算できる。これはn次の線形最小二乗問題である。

以上からGMRES法のアルゴリズムが得られる[1]。すなわち、以下の反復を実行する。

  1. アーノルディ法を1ステップ計算する
  2. ||rn||を最小化する y n {\displaystyle y_{n}} を見つける
  3. x n = Q n y n {\displaystyle x_{n}=Q_{n}y_{n}} を計算する
  4. 残差が十分小さくなるまでこれを繰り返す

各反復で、行列ベクトル積Aqnを計算する必要がある。これはm次の一般密行列では約2m2回の浮動小数点数演算を要するが、疎行列ではO(m)とすることができる。行列ベクトル積の他に、第nステップの計算にはO(n m)の浮動小数演算が必要である。

収束特性

第nステップでは、クリロフ部分空間Kn内での残差が最小化される。部分空間は常に次の部分空間に含まれるので、残差は単調に減少する。Aの次数をmとすると、m回目の反復の後ではクリロフ部分空間KmはRmと等しいので、GMRES法は厳密解に到達する。しかし、アイデアの骨子は、（mと比較して）少ない回数の反復でベクトルxnが十分な解の近似になる点にある。

一般にはこれは成り立たない。実際、Greenbaum・Pták・Strakošの定理[4]によれば、任意の単調減少列a1, …, am−1、am = 0について、上で定義されたrnに関して、すべてのnで||rn|| = anとなるような行列Aを見つけることができる。特に、m − 1回の反復で一定の値を保ちながら、最後の1反復で残差が0になるような行列を見つけることができる。

ただ、実際にはGMRES法は良い性能を示すことも多い。これは特定の場合に証明できる。Aが正定値なら、 λ m i n {\displaystyle \lambda _{\mathrm {min} }} 、 λ m a x {\displaystyle \lambda _{\mathrm {max} }} をそれぞれ最小、最大固有値として、

‖ r n ‖ ≤ ( 1 − λ m i n ( A ⊤ + A ) 2 λ m a x ( A T + A ) ) n / 2 ‖ r 0 ‖ {\displaystyle \|r_{n}\|\leq \left(1-{\frac {\lambda _{\mathrm {min} }(A^{\top }+A)}{2\lambda _{\mathrm {max} }(A^{T}+A)}}\right)^{n/2}\|r_{0}\|}

が成り立つ。

Aが対称正定値なら、 κ 2 ( A ) {\displaystyle \kappa _{2}(A)} をAのユークリッドノルムでの条件数として、

‖ r n ‖ ≤ ( κ 2 2 ( A ) − 1 κ 2 2 ( A ) ) n / 2 ‖ r 0 ‖ {\displaystyle \|r_{n}\|\leq \left({\frac {\kappa _{2}^{2}(A)-1}{\kappa _{2}^{2}(A)}}\right)^{n/2}\|r_{0}\|}

が成り立つ。

Aが正定値でない場合には、Pnをp(0) = 1を満たす高々n次の多項式集合、VをAのスペクトル分解に現れる行列、σ(A)をAのスペクトルとして、

‖ r n ‖ ≤ inf p ∈ P n ‖ p n ( A ) ‖ ≤ κ 2 ( V ) inf p ∈ P n max λ ∈ σ ( A ) | p ( λ ) | {\displaystyle \|r_{n}\|\leq \inf _{p\in P_{n}}\|p_{n}(A)\|\leq \kappa _{2}(V)\inf _{p\in P_{n}}\max _{\lambda \in \sigma (A)}|p(\lambda )|}

を得る。おおざっぱに言えば、これはAの固有値が0から遠くかつ密集しており、Aが正規行列からそれほど離れていない場合に、速く収束することを意味している[5]。

これらの不等式は誤差、つまり現在の反復ベクトルxnと真の解との距離ではなく、残差に関するものである。

解法の拡張

前処理

他の反復法と同様に、GMRES法も通常は収束を速める目的で前処理と組み合わせて使用される[6][7]。

リスタート

nを反復回数として、反復のコストはO(n2)である。すなわち、nとして、連立方程式の本数Nを用いると、直接解法と同程度の求解コストがかかることになる。したがって、GMRES法ではk回の反復の後に、xkを初期ベクトルとして、リスタートを行うことがある。これをGMRES(k)と呼ぶ。しかしながら、リスタートを行うと、収束は非常に遅くなることが知られている。

Flexible GMRES法

前処理をより効率的に行えるよう、アルゴリズムに変更を加えた手法である[8]。前処理に反復解法を用いることができる(前処理にGMRES自身を用いるなど)。

他の解法との比較

アーノルディ法は対称行列の場合にはランチョス法[9]に帰着する。対応するクリロフ部分空間法はPaige・SaundersによるMINRES法（minimal residual method[10][11][12]）である。非対称な場合と異なり、MINRES法は3項漸化式で与えられる。一般行列の場合には、GMRES法のように短い漸化式で残差ノルムを最小化するようなクリロフ部分空間法は存在しないことが知られている。

別な系統として、非対称ランチョス法、特に双共役勾配法[13]に基づく解法がある。これらは3項漸化式を用いるが、最小残差は計算しない。したがって残差は単調には減少せず、収束も保証されない。

さらに、CGS法（conjugate gradient squared method[14]）やBiCGSTAB法（stabilized biconjugate gradient method[15]）などによる解法がある。これらも3項漸化式に基づいているため、最適性は保証されず、解に収束する前に終了することがある。これらの解法のアイデアは、反復毎の生成多項式を適切に選択する点にある。

これらのいずれも、すべての行列に対して万能というわけではない。したがって、実際には与えられた問題に対して複数の解法を試みる必要がある。

最小二乗問題の求解

GMRES法では、

‖ H ~ n y n − e 1 ‖ {\displaystyle \|{\tilde {H}}_{n}y_{n}-e_{1}\|}

を最小化するベクトル y n {\displaystyle y_{n}} を見つける必要がある。 これはQR分解を計算することで実現できる。すなわち、

Ω n H ~ n = R ~ n . {\displaystyle \Omega _{n}{\tilde {H}}_{n}={\tilde {R}}_{n}.}

を満たすn+1次正方直交行列Ωnとn+1×n次上三角行列 R ~ n {\displaystyle {\tilde {R}}_{n}} を計算する。三角行列の行数は列数より1多いので、最下行はすべて零である。したがって、 R n {\displaystyle R_{n}} をn×n次三角行列として、これを

R ~ n = [ R n 0 ] , {\displaystyle {\tilde {R}}_{n}={\begin{bmatrix}R_{n}\\0\end{bmatrix}},}

と分解することができる。

QR分解は、零の行と1列分の値が異なるだけなので、簡単に更新することができる。つまりhn = (h1n, … hnn)Tとして

H ~ n + 1 = [ H ~ n h n 0 h n + 1 , n ] {\displaystyle {\tilde {H}}_{n+1}={\begin{bmatrix}{\tilde {H}}_{n}&h_{n}\\0&h_{n+1,n}\end{bmatrix}}}

が成り立つので、

[ Ω n 0 0 1 ] H ~ n + 1 = [ R n r k 0 ρ 0 σ ] {\displaystyle {\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{k}\\0&\rho \\0&\sigma \end{bmatrix}}}

は三角行列に近く、σが0であれば三角行列である。これを更新するためには、: c n = ρ ρ 2 + σ 2 and s n = σ ρ 2 + σ 2 {\displaystyle c_{n}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\mbox{and}}\quad s_{n}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}} としてギブンス回転

G n = [ I n − 1 0 0 0 c n s n 0 − s n c n ] {\displaystyle G_{n}={\begin{bmatrix}I_{n-1}&0&0\\0&c_{n}&s_{n}\\0&-s_{n}&c_{n}\end{bmatrix}}}

を行えばよい。これにより、

Ω n + 1 = G n [ Ω n 0 0 1 ] {\displaystyle \Omega _{n+1}=G_{n}{\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}}

を得る。

Ω n + 1 H ~ n + 1 = [ R n r n 0 r n n 0 0 ] with r n n = ρ 2 + σ 2 {\displaystyle \Omega _{n+1}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{n}\\0&r_{nn}\\0&0\end{bmatrix}}\quad {\text{with}}\quad r_{nn}={\sqrt {\rho ^{2}+\sigma ^{2}}}}

は三角行列である。

QR分解が与えられると、最小化問題は

‖ H ~ n y n − e 1 ‖ = ‖ Ω n ( H ~ n y n − e 1 ) ‖ = ‖ R ~ n y n − Ω n e 1 ‖ {\displaystyle \|{\tilde {H}}_{n}y_{n}-e_{1}\|=\|\Omega _{n}({\tilde {H}}_{n}y_{n}-e_{1})\|=\|{\tilde {R}}_{n}y_{n}-\Omega _{n}e_{1}\|}

から容易に解くことができる。実際、gn ∈ Rn、γn ∈ Rとして、ベクトル Ω n e 1 {\displaystyle \Omega _{n}e_{1}} を

g ~ n = [ g n γ n ] {\displaystyle {\tilde {g}}_{n}={\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}}

と表すと、これは

‖ H ~ n y n − e 1 ‖ = ‖ R ~ n y n − Ω n e 1 ‖ = ‖ [ R n 0 ] y − [ g n γ n ] ‖ {\displaystyle \|{\tilde {H}}_{n}y_{n}-e_{1}\|=\|{\tilde {R}}_{n}y_{n}-\Omega _{n}e_{1}\|=\left\|{\begin{bmatrix}R_{n}\\0\end{bmatrix}}y-{\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}\right\|}

と書ける。これを最小化するベクトルyは

y n = R n − 1 g n {\displaystyle y_{n}=R_{n}^{-1}g_{n}}

で与えられる。 g n {\displaystyle g_{n}} の更新も容易に行うことができる[16]。

補足

  1. ^ a b Black, Noel and Moore, Shirley. "Generalized Minimal Residual Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. mathworld.wolfram.com/GeneralizedMinimalResidualMethod.html
  2. ^ W. E. Arnoldi (1951), "The principle of minimized iterations in the solution of the matrix eigenvalue problem," Quarterly of Applied Mathematics, volume 9, pages 17–29.
  3. ^ Saad, Y., & Schultz, M. H. (1986). GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing, 7(3), 856-869.
  4. ^ Greenbaum, A., Pták, V., & Strakoš, Z. E. K. (1996). Any nonincreasing convergence curve is possible for GMRES. SIAM journal on matrix analysis and applications, 17(3), 465-469.
  5. ^ Trefethen & Bau, Thm 35.2
  6. ^ Baglama, J., Calvetti, D., Golub, G. H., & Reichel, L. (1998). Adaptively preconditioned GMRES algorithms. SIAM Journal on Scientific Computing, 20(1), 243-269.
  7. ^ Burrage, K., & Erhel, J. (1998). On the performance of various adaptive preconditioned GMRES strategies. Numerical linear algebra with applications, 5(2), 101-121.
  8. ^ Frayssé, V., Giraud, L., & Gratton, S. (1998). A set of Flexible-GMRES routines for real and complex arithmetics. Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique TR/PA/98/20.
  9. ^ Lanczos, C. (1950). "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators". Journal of Research of the National Bureau of Standards. 45 (4): 255–282.
  10. ^ Black, Noel and Moore, Shirley. "Minimal Residual Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. mathworld.wolfram.com/MinimalResidualMethod.html
  11. ^ Paige, C. and Saunders, M. "Solution of Sparse Indefinite Systems of Linear Equations." SIAM J. Numer. Anal. 12, 617-629, 1975.
  12. ^ Fong, D. C. L., & Saunders, M. (2012). CG versus MINRES: An empirical comparison. Sultan Qaboos University Journal for Science [SQUJS], 17(1), 44-62.
  13. ^ Black, Noel and Moore, Shirley. "Biconjugate Gradient Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. mathworld.wolfram.com/BiconjugateGradientMethod.html
  14. ^ Black, Noel and Moore, Shirley. "Conjugate Gradient Squared Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. mathworld.wolfram.com/ConjugateGradientSquaredMethod.html
  15. ^ Van der Vorst, H. A. (1992). Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM Journal on scientific and Statistical Computing, 13(2), 631-644.
  16. ^ Stoer and Bulirsch, §8.7.2

参考文献

  • A. Meister, Numerik linearer Gleichungssysteme, 2nd edition, Vieweg 2005, ISBN 978-3-528-13135-7.
  • Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd edition, en:Society for Industrial and Applied Mathematics, 2003. ISBN 978-0-89871-534-7.
  • Y. Saad and M.H. Schultz, "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems", SIAM J. Sci. Stat. Comput., 7:856-869, 1986. doi:10.1137/0907058.
  • J. Stoer and R. Bulirsch, Introduction to numerical analysis, 3rd edition, Springer, New York, 2002. ISBN 978-0-387-95452-3.
  • Lloyd N. Trefethen and David Bau, III, Numerical Linear Algebra, Society for Industrial and Applied Mathematics, 1997. ISBN 978-0-89871-361-9.
  • Dongarra et al. , Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition, SIAM, Philadelphia, 1994.
  • 表
  • 話
  • 編
  • 歴
連立一次方程式
ベクトル
ベクトル空間
計量ベクトル空間
行列・線型写像
演算・操作
不変量
クラス
行列式
多重線型代数
数値線形代数
基本的な概念
ソフトウェア
ライブラリ
反復法・技法
人物
行列値関数
その他
カテゴリ カテゴリ