GMRES法

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

 解法 

解くべき方程式を

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

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

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

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

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

アーノルディ過程では、

を満たす(n+1)×n次上ヘッセンベルグ行列が生成される。は直交行列なので、

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

が成り立つ。したがって、は残差

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

以上からGMRES法のアルゴリズムが得られる。すなわち、以下の反復を実行する。
 アーノルディ法を1ステップ計算する
 ||rn||を最小化するを見つける
 を計算する
 残差が十分小さくなるまでこれを繰り返す

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

 収束特性 

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

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

ただ、実際にはGMRES法は良い性能を示すことも多い。これは特定の場合に証明できる。Aが正定値なら、、をそれぞれ最小、最大固有値として、

が成り立つ。

Aが対称正定値なら、をAのユークリッドノルムでの条件数として、

が成り立つ。

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

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

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

 解法の拡張 

 前処理 

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

 リスタート 

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

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

 他の解法との比較 

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

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

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

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

 最小二乗問題の求解 

GMRES法では、

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

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

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

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

が成り立つので、

は三角行列に近く、σが0であれば三角行列である。これを更新するためには、:
としてギブンス回転

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

を得る。

は三角行列である。

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

から容易に解くことができる。実際、gn ∈ Rn、γn ∈ Rとして、ベクトルを

と表すと、これは

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

で与えられる。の更新も容易に行うことができる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. .
 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.

Category:数値線形代数
Category:数値解析
Category:応用数学
Category:計算科学
Category:アルゴリズム