GMRES法のソースを表示

数学において、'''GMRES法'''（GMRESほう、'''generalized minimal residual method'''）は、[[連立一次方程式]]の数値解を求めるための[[反復法 (数値計算)|反復法]]の一種である<ref name="mw">Black, Noel and Moore, Shirley. "Generalized Minimal Residual Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. {{URL|https://mathworld.wolfram.com/GeneralizedMinimalResidualMethod.html}}</ref>。残差を[[クリロフ部分空間]]において最小化することにより、近似解を計算する。ベクトルの計算にはアーノルディ法<ref>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.</ref>が用いられる。ヨセフ・サードとマルティン・H・シュルツにより、1986年に開発された<ref>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.</ref>。

== 解法 ==

解くべき方程式を

{{Indent|<math> Ax = b. \, </math>}}

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

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

{{Indent|<math> K_n = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{n-1}b \}. \, </math>}}

である。GMRES法では、残差''Ax''<sub>''n''</sub> − ''b''を最小化するベクトル''x''<sub>''n''</sub> ∈ ''K''<sub>''n''</sub>によって、''Ax'' = ''b''の厳密解を近似する。

''b'', ''Ab'', …, ''A''<sup>''n''−1</sup>''b''は線形従属に近いため、これを用いる代わりに、[[アーノルディ法]]を用いて''K''<sub>''n''</sub>の基底を構成する正規直交化ベクトル列

{{Indent|<math> q_1, q_2, \ldots, q_n \, </math>}}

を計算する。すなわち、''x''<sub>''n''</sub> ∈ ''K''<sub>''n''</sub>は''y''<sub>''n''</sub> ∈ '''R'''<sup>''n''</sup>を用いて''x''<sub>''n''</sub> = ''Q''<sub>''n''</sub>''y''<sub>''n''</sub>と書くことができる。ただし''Q''<sub>''n''</sub>は''q''<sub>1</sub>, …, ''q''<sub>n</sub>により構成される''m''×''n''行列とする。 

アーノルディ過程では、

{{Indent|<math> AQ_n = Q_{n+1} \tilde{H}_n</math>}}

を満たす(''n''+1)×''n''次上[[ヘッセンベルグ行列]]<math>\tilde{H}_n</math>が生成される。<math>Q_n</math>は直交行列なので、

{{Indent|<math> e_1 = (1,0,0,\ldots,0) \, </math>}}

を標準基底'''R'''<sup>''n''+1</sup>の第1ベクトルとして、

{{Indent|<math> \| Ax_n - b \| = \| \tilde{H}_ny_n - e_1 \|</math>}}

が成り立つ。したがって、<math>x_n</math>は残差

{{Indent|<math> r_n = \tilde{H}_n y_n - e_1. </math>}}

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

以上からGMRES法のアルゴリズムが得られる<ref name="mw"/>。すなわち、以下の反復を実行する。
# アーノルディ法を1ステップ計算する
# ||''r''<sub>''n''</sub>||を最小化する<math> y_n </math>を見つける
# <math> x_n = Q_n y_n </math>を計算する
# 残差が十分小さくなるまでこれを繰り返す

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

== 収束特性 ==

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

一般にはこれは成り立たない。実際、Greenbaum・Pták・Strakošの定理<ref>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.</ref>によれば、任意の単調減少列''a''<sub>1</sub>, …, ''a''<sub>''m''−1</sub>、''a''<sub>''m''</sub> = 0について、上で定義された''r''<sub>''n''</sub>に関して、すべての''n''で||''r''<sub>''n''</sub>|| = ''a''<sub>''n''</sub>となるような行列''A''を見つけることができる。特に、''m''&nbsp;−&nbsp;1回の反復で一定の値を保ちながら、最後の1反復で残差が0になるような行列を見つけることができる。

ただ、実際にはGMRES法は良い性能を示すことも多い。これは特定の場合に証明できる。''A''が正定値なら、<math>\lambda_{\mathrm{min}}</math>、<math>\lambda_{\mathrm{max}}</math>をそれぞれ最小、最大固有値として、
{{Indent|<math> \|r_n\| \leq \left( 1-\frac{\lambda_{\mathrm{min}}(A^\top + A)}{2 \lambda_{\mathrm{max}}(A^T + A)} \right)^{n/2} \|r_0\|</math>}}
が成り立つ。

''A''が対称正定値なら、<math>\kappa_2(A)</math>を''A''のユークリッドノルムでの条件数として、
{{Indent|<math> \|r_n\| \leq \left( \frac{\kappa_2^2(A)-1}{\kappa_2^2(A)} \right)^{n/2} \|r_0\|</math>}}
が成り立つ。

''A''が正定値でない場合には、''P''<sub>''n''</sub>を''p''(0) = 1を満たす高々''n''次の多項式集合、''V''を''A''のスペクトル分解に現れる行列、σ(''A'')を''A''のスペクトルとして、
{{Indent|<math> \|r_n\| \le \inf_{p \in P_n} \|p_n(A)\| \le \kappa_2(V) \inf_{p \in P_n} \max_{\lambda \in \sigma(A)} |p(\lambda)|</math>}}
を得る。おおざっぱに言えば、これは''A''の固有値が0から遠くかつ密集しており、''A''が正規行列からそれほど離れていない場合に、速く収束することを意味している<ref>Trefethen & Bau, Thm 35.2</ref>。

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

== 解法の拡張 ==

=== 前処理 ===

他の反復法と同様に、GMRES法も通常は収束を速める目的で[[前処理行列|前処理]]と組み合わせて使用される<ref>Baglama, J., Calvetti, D., Golub, G. H., & Reichel, L. (1998). Adaptively preconditioned GMRES algorithms. SIAM Journal on Scientific Computing, 20(1), 243-269.</ref><ref>Burrage, K., & Erhel, J. (1998). On the performance of various adaptive preconditioned GMRES strategies. Numerical linear algebra with applications, 5(2), 101-121.</ref>。

=== リスタート ===

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

=== Flexible GMRES法 ===
前処理をより効率的に行えるよう、アルゴリズムに変更を加えた手法である<ref>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.
</ref>。前処理に反復解法を用いることができる(前処理にGMRES自身を用いるなど)。

== 他の解法との比較 ==

アーノルディ法は対称行列の場合には[[ランチョス法]]<ref>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.</ref>に帰着する。対応するクリロフ部分空間法はPaige・SaundersによるMINRES法（minimal residual method<ref>Black, Noel and Moore, Shirley. "Minimal Residual Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. {{URL|https://mathworld.wolfram.com/MinimalResidualMethod.html}}</ref><ref>Paige, C. and Saunders, M. "Solution of Sparse Indefinite Systems of Linear Equations." SIAM J. Numer. Anal. 12, 617-629, 1975.</ref><ref>Fong, D. C. L., & Saunders, M. (2012). CG versus MINRES: An empirical comparison. Sultan Qaboos University Journal for Science [SQUJS], 17(1), 44-62.</ref>）である。非対称な場合と異なり、MINRES法は3項漸化式で与えられる。一般行列の場合には、GMRES法のように短い漸化式で残差ノルムを最小化するようなクリロフ部分空間法は存在しないことが知られている。

別な系統として、非対称ランチョス法、特に[[双共役勾配法]]<ref>Black, Noel and Moore, Shirley. "Biconjugate Gradient Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. {{URL|https://mathworld.wolfram.com/BiconjugateGradientMethod.html}}</ref>に基づく解法がある。これらは3項漸化式を用いるが、最小残差は計算しない。したがって残差は単調には減少せず、収束も保証されない。

さらに、CGS法（conjugate gradient squared method<ref>Black, Noel and Moore, Shirley. "Conjugate Gradient Squared Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. {{URL|https://mathworld.wolfram.com/ConjugateGradientSquaredMethod.html}}</ref>）やBiCGSTAB法（stabilized biconjugate gradient method<ref>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.</ref>）などによる解法がある。これらも3項漸化式に基づいているため、最適性は保証されず、解に収束する前に終了することがある。これらの解法のアイデアは、反復毎の生成多項式を適切に選択する点にある。

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

== 最小二乗問題の求解 ==

GMRES法では、
{{Indent|<math> \| \tilde{H}_n y_n - e_1 \|</math>}}
を最小化するベクトル<math>y_n</math>を見つける必要がある。
これは[[QR分解]]を計算することで実現できる。すなわち、
{{Indent|<math> \Omega_n \tilde{H}_n = \tilde{R}_n. </math>}}
を満たす''n''+1次正方[[直交行列]]Ω<sub>''n''</sub>と''n''+1×''n''次上[[三角行列]]<math>\tilde{R}_n</math>を計算する。三角行列の行数は列数より1多いので、最下行はすべて零である。したがって、<math>R_n</math>を''n''×''n''次三角行列として、これを
{{Indent|<math> \tilde{R}_n = \begin{bmatrix} R_n \\ 0 \end{bmatrix}, </math>}}
と分解することができる。

QR分解は、零の行と1列分の値が異なるだけなので、簡単に更新することができる。つまり''h''<sub>''n''</sub> = (''h''<sub>1''n''</sub>, … ''h''<sub>''nn''</sub>)<sup>T</sup>として
{{Indent|<math>\tilde{H}_{n+1} = \begin{bmatrix} \tilde{H}_n & h_n \\ 0 & h_{n+1,n} \end{bmatrix}</math>}}
が成り立つので、
{{Indent|<math> \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} </math>}}
は三角行列に近く、σが0であれば三角行列である。これを更新するためには、:<math> c_n = \frac{\rho}{\sqrt{\rho^2+\sigma^2}} \quad\mbox{and}\quad s_n = \frac{\sigma}{\sqrt{\rho^2+\sigma^2}}</math>
として[[ギブンス回転]]
{{Indent|<math> G_n = \begin{bmatrix} I_{n-1} & 0 & 0 \\ 0 & c_n & s_n \\ 0 & -s_n & c_n \end{bmatrix} </math>}}
を行えばよい。これにより、
{{Indent|<math> \Omega_{n+1} = G_n \begin{bmatrix} \Omega_n & 0 \\ 0 & 1 \end{bmatrix}</math>}}
を得る。
{{Indent|<math> \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} </math>}}
は三角行列である。

QR分解が与えられると、最小化問題は
{{Indent|<math> \| \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 \|</math>}}
から容易に解くことができる。実際、''g''<sub>''n''</sub> ∈ '''R'''<sub>''n''</sub>、γ<sub>''n''</sub> ∈ '''R'''として、ベクトル<math>\Omega_ne_1</math>を
{{Indent|<math> \tilde{g}_n = \begin{bmatrix} g_n \\ \gamma_n \end{bmatrix} </math>}}
と表すと、これは
{{Indent|<math> \| \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\|</math>}}
と書ける。これを最小化するベクトル''y''は
{{Indent|<math> y_n = R_n^{-1} g_n</math>}}
で与えられる。<math>g_n</math>の更新も容易に行うことができる<ref>Stoer and Bulirsch, §8.7.2</ref>。

== 補足 ==
{{reflist|2}}

== 参考文献 ==
{{div col|rules=yes}}
* 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.
* [http://www.netlib.org/linalg/html_templates/node29.html#SECTION00734000000000000000|J. Dongarra et al. , ''Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods''], 2nd Edition, SIAM, Philadelphia, 1994.
{{div col end}}
{{linear algebra}}
[[Category:数値線形代数]]
[[Category:数値解析]]
[[Category:応用数学]]
[[Category:計算科学]]
[[Category:アルゴリズム]]