GMRES法

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

解くべき方程式を

$Ax=b.\,$

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

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

$K_{n}=\operatorname {span} \,\{b,Ab,A^{2}b,\ldots ,A^{n-1}b\}.\,$

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

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

$q_{1},q_{2},\ldots ,q_{n}\,$

を計算する。すなわち、x_n ∈ K_nはy_n ∈ Rⁿを用いてx_n = Q_ny_nと書くことができる。ただしQ_nはq₁, …, q_nにより構成されるm×n行列とする。

アーノルディ過程では、

$AQ_{n}=Q_{n+1}{\tilde {H}}_{n}$

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

$e_{1}=(1,0,0,\ldots ,0)\,$

を標準基底Rⁿ⁺¹の第1ベクトルとして、

$\|Ax_{n}-b\|=\|{\tilde {H}}_{n}y_{n}-e_{1}\|$

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

$r_{n}={\tilde {H}}_{n}y_{n}-e_{1}.$

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

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

アーノルディ法を1ステップ計算する
||r_n||を最小化する $y_{n}$ を見つける
$x_{n}=Q_{n}y_{n}$ を計算する
残差が十分小さくなるまでこれを繰り返す

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

収束特性

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

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

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

$\|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が対称正定値なら、 $\kappa _{2}(A)$ をAのユークリッドノルムでの条件数として、

$\|r_{n}\|\leq \left({\frac {\kappa _{2}^{2}(A)-1}{\kappa _{2}^{2}(A)}}\right)^{n/2}\|r_{0}\|$

が成り立つ。

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

$\|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]。

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

解法の拡張

前処理

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

リスタート

nを反復回数として、反復のコストはO(n²)である。すなわち、nとして、連立方程式の本数Nを用いると、直接解法と同程度の求解コストがかかることになる。したがって、GMRES法ではk回の反復の後に、x_kを初期ベクトルとして、リスタートを行うことがある。これを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法では、

$\|{\tilde {H}}_{n}y_{n}-e_{1}\|$

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

$\Omega _{n}{\tilde {H}}_{n}={\tilde {R}}_{n}.$

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

${\tilde {R}}_{n}={\begin{bmatrix}R_{n}\\0\end{bmatrix}},$

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

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

${\tilde {H}}_{n+1}={\begin{bmatrix}{\tilde {H}}_{n}&h_{n}\\0&h_{n+1,n}\end{bmatrix}}$

が成り立つので、

${\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}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\mbox{and}}\quad s_{n}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}$ としてギブンス回転

$G_{n}={\begin{bmatrix}I_{n-1}&0&0\\0&c_{n}&s_{n}\\0&-s_{n}&c_{n}\end{bmatrix}}$

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

$\Omega _{n+1}=G_{n}{\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}$

を得る。

$\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分解が与えられると、最小化問題は

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

から容易に解くことができる。実際、g_n ∈ R_n、γ_n ∈ Rとして、ベクトル $\Omega _{n}e_{1}$ を

${\tilde {g}}_{n}={\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}$

と表すと、これは

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

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

補足

1 2 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
↑ 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.
↑ 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.
↑ 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.
↑ Trefethen & Bau, Thm 35.2
↑ 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.
↑ 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.
↑ 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.
↑ Black, Noel and Moore, Shirley. "Minimal Residual Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. mathworld.wolfram.com/MinimalResidualMethod.html
↑ 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.
↑ Black, Noel and Moore, Shirley. "Biconjugate Gradient Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. mathworld.wolfram.com/BiconjugateGradientMethod.html
↑ 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
↑ 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.
↑ Stoer and Bulirsch, §8.7.2