一般化最小残差法

Method for numerical solution of certain systems of equations

数学において、一般化最小残差法（GMRES）は、非正定値非対称線形方程式の数値解を求める反復法である。この法は、解を最小残差を持つクリロフ部分空間のベクトルで近似する。このベクトルを求めるには、アーノルディ反復法が用いられる。

GMRES法は、 1986年にYousef SaadとMartin H. Schultzによって開発されました^。[1]これは、1975年にPaigeとSaundersによって開発されたMINRES法を一般化・改良したものです。 ^{[2] [3]} MINRES法では行列が対称である必要がありますが、3つのベクトルのみを扱うことができるという利点があります。GMRES法は、 1980年にPeter Pulayによって開発されたDIIS法の特殊なケースです。DIIS法は非線形システムに適用可能です。

方法

任意のベクトルvのユークリッドノルムを と表記する。解くべき（正方）線形方程式系を と表記する。行列 Aはm行m列の逆行列であると仮定する。さらに、 bは正規化されている、すなわち であると仮定する。 ${\displaystyle \|v\|}$ $${\displaystyle Ax=b.}$$ ${\displaystyle \|b\|=1}$

この問題の n 番目のクリロフ部分空間は、初期推定値が与えられた 場合の初期残差です。明らかに、であれば、 $${\displaystyle K_{n}=K_{n}(A,r_{0})=\operatorname {span} \,\{r_{0},Ar_{0},A^{2}r_{0},\ldots ,A^{n-1}r_{0}\}.\,}$$ ${\displaystyle r_{0}=b-Ax_{0}}$ ${\displaystyle x_{0}\neq 0}$ ${\displaystyle r_{0}=b}$ ${\displaystyle x_{0}=0}$

GMRES は、残差のユークリッドノルムを最小化するベクトルによっての正確な解を近似します。 ${\displaystyle Ax=b}$ ${\displaystyle x_{n}\in x_{0}+K_{n}}$ ${\displaystyle r_{n}=b-Ax_{n}}$

ベクトルは線形従属に近い可能性があるので、この基底の代わりに、アーノルディ反復法を使用しての基底を形成する直交ベクトルを見つけます。特に、 です。 ${\displaystyle r_{0},Ar_{0},\ldots A^{n-1}r_{0}}$ ${\displaystyle q_{1},q_{2},\ldots ,q_{n}\,}$ ${\displaystyle K_{n}}$ ${\displaystyle q_{1}=\|r_{0}\|_{2}^{-1}r_{0}}$

したがって、ベクトルはと表すことができます。ここで、は によって形成されるm行n列の行列です。言い換えれば、解のn番目の近似値（すなわち）を求めることは、ベクトル を求めることに帰着します。ベクトル は、以下に示すように、留数を最小化することで決定されます。 ${\displaystyle x_{n}\in x_{0}+K_{n}}$ ${\displaystyle x_{n}=x_{0}+Q_{n}y_{n}}$ ${\displaystyle y_{n}\in \mathbb {R} ^{n}}$ ${\displaystyle Q_{n}}$ ${\displaystyle q_{1},\ldots ,q_{n}}$ ${\displaystyle x_{n}}$ ${\displaystyle y_{n}}$

アーノルディ過程は、( )行上ヘッセンベルク行列も構築します。この行列 は、 の計算を簡略化するために使用される等式を満たします（§ 最小二乗問題の解法を参照）。対称行列の場合、実際には対称三重対角行列が得られ、MINRES法が用いられることに注意してください。 ${\displaystyle {\tilde {H}}_{n}}$ ${\displaystyle n+1}$ ${\displaystyle n}$ $${\displaystyle AQ_{n}=Q_{n+1}{\tilde {H}}_{n}\,}$$ ${\displaystyle y_{n}}$

の列は直交しているので、 はの標準基底の最初のベクトルであり、は最初の試行残差ベクトル（通常は）である。したがって、 は残差のユークリッドノルムを最小化することで求められる。これはサイズnの線形最小二乗問題である。 ${\displaystyle Q_{n}}$ $${\displaystyle {\begin{aligned}\left\|r_{n}\right\|&=\left\|b-Ax_{n}\right\|\\&=\left\|b-A(x_{0}+Q_{n}y_{n})\right\|\\&=\left\|r_{0}-AQ_{n}y_{n}\right\|\\&=\left\|\beta q_{1}-AQ_{n}y_{n}\right\|\\&=\left\|\beta q_{1}-Q_{n+1}{\tilde {H}}_{n}y_{n}\right\|\\&=\left\|Q_{n+1}(\beta e_{1}-{\tilde {H}}_{n}y_{n})\right\|\\&=\left\|\beta e_{1}-{\tilde {H}}_{n}y_{n}\right\|\end{aligned}}}$$ $${\displaystyle e_{1}=(1,0,0,\ldots ,0)^{T}\,}$$ ${\displaystyle \mathbb {R} ^{n+1}}$ $${\displaystyle \beta =\|r_{0}\|\,,}$$ ${\displaystyle r_{0}}$ ${\displaystyle b}$ ${\displaystyle x_{n}}$ $${\displaystyle r_{n}={\tilde {H}}_{n}y_{n}-\beta e_{1}.}$$

これによりGMRES法が生成される。第回の反復では次のようになる。 ${\displaystyle n}$

1. アーノルディ法で計算します。 ${\displaystyle q_{n}}$
2. を最小化する を見つけます。 ${\displaystyle y_{n}}$ ${\displaystyle \|r_{n}\|}$
3. 計算する; ${\displaystyle x_{n}=x_{0}+Q_{n}y_{n}}$
4. 残差がまだ十分に小さくない場合は繰り返します。

各反復処理において、行列ベクトル積を計算する必要があります。これは、サイズ の一般的な密行列では約 回の浮動小数点演算を必要としますが、疎行列ではコストがまで減少する可能性があります。行列ベクトル積に加えて、n回目の反復処理では浮動小数点演算を計算する必要があります。 ${\displaystyle Aq_{n}}$ ${\displaystyle 2m^{2}}$ ${\displaystyle m}$ ${\displaystyle O(m)}$ ${\displaystyle O(nm)}$

収束

n回目の反復は、クリロフ部分空間 における残差を最小化します。すべての部分空間は次の部分空間に含まれるため、残差は増加しません。m回の反復（mは行列Aのサイズ）の後、クリロフ空間K _mはR ^m全体となり、したがって GMRES 法は厳密解に到達します。ただし、この手法の考え方は、（ mに対して）少数の反復後には、ベクトルx _{n が}既に厳密解の良い近似値になっているというものです。 ${\displaystyle K_{n}}$

これは一般には起こりません。実際、Greenbaum、Pták、Strakošの定理によれば、任意の非増加列a ₁ , ..., a _{m −1} , a _m = 0 に対して、すべてのnに対して ‖ r _n ‖ = a _nとなるような行列Aが見つかります。ここで、 r _nは上で定義した残差です。特に、残差がm  − 1 回の反復処理で一定のままで、最後の反復処理でのみゼロになるような行列が見つかる可能性があります。

しかし実際には、GMRESは多くの場合良好な性能を示します。これは特定の状況で証明できます。A の対称部分、つまり が正定値 である場合、となります。ここで、 と はそれぞれ行列 の最小固有値と最大固有値を表します。^[4] ${\displaystyle (A^{T}+A)/2}$ $${\displaystyle \|r_{n}\|\leq \left(1-{\frac {\lambda _{\min }^{2}(1/2(A^{T}+A))}{\lambda _{\max }(A^{T}A)}}\right)^{n/2}\|r_{0}\|,}$$ ${\displaystyle \lambda _{\mathrm {min} }(M)}$ ${\displaystyle \lambda _{\mathrm {max} }(M)}$ ${\displaystyle M}$

Aが対称かつ正定値である場合、 が成り立ちます。ここで はユークリッドノルムにおけるAの条件数を表します。 $${\displaystyle \|r_{n}\|\leq \left({\frac {\kappa _{2}(A)^{2}-1}{\kappa _{2}(A)^{2}}}\right)^{n/2}\|r_{0}\|.}$$ ${\displaystyle \kappa _{2}(A)}$

Aが正定値でない一般的なケースでは、次式が成り立ちます。ここで、 P _{n は}p (0) = 1 を満たす n 次以下の多項式の集合、VはAのスペクトル分解に現れる行列、σ ( A ) はAのスペクトルです。大まかに言えば、これはAの固有値が原点から離れて密集し、Aが正規分布からそれほど離れていないときに、高速収束が起こることを示しています。^[5] $${\displaystyle {\frac {\|r_{n}\|}{\|b\|}}\leq \inf _{p\in P_{n}}\|p(A)\|\leq \kappa _{2}(V)\inf _{p\in P_{n}}\max _{\lambda \in \sigma (A)}|p(\lambda )|,\,}$$

これらの不等式はすべて、実際の誤差ではなく残差、つまり現在の反復x _nと正確な解の間の距離のみを制限します。

方法の拡張

他の反復​​法と同様に、GMRES は通常、収束を高速化するために前処理法と組み合わせられます。

反復のコストはO( n^ ² )で増加する（nは反復回数）。そのため、この手法は、例えばk回の反復後に、初期推定値xk_を用いて再計算されることがある。この手法はGMRES( k )またはRestarted GMRESと呼ばれる。非正定値行列の場合、再計算後の部分空間が以前の部分空間に近い場合が多いため、この手法は収束が停滞する可能性がある。

GMRES法とリスタートGMRES法の欠点は、GCROT法やGCRODR法などのGCRO型法におけるクリロフ部分空間の再利用によって解決される。^[6] GMRES法におけるクリロフ部分空間の再利用は、線形システムのシーケンスを解く必要がある場合に収束を高速化することもできる。^[7]

他のソルバーとの比較

対称行列の場合、アーノルディ反復法はランチョス反復法に帰着します。対応するクリロフ部分空間法は、ペイジとサンダースによる最小残差法（MinRes）です。非対称行列の場合とは異なり、MinRes法は3項漸化式で与えられます。一般行列の場合、GMRES法のように短い漸化式で与えられ、かつ残差ノルムを最小化するクリロフ部分空間法は存在しないことが示されます。

非対称ランチョス反復法を基盤とする別の手法、特にBiCG法が注目される。これらの手法では3項漸化式を用いるが、残差の最小化は達成されず、したがって残差は単調に減少しない。収束性も保証されていない。

3番目のクラスは、CGSやBiCGSTABなどの手法によって形成されます。これらも3項再帰関係（したがって最適性なし）を扱い、収束せずに途中で終了することもあります。これらの手法の背後にある考え方は、反復列の生成多項式を適切に選択することです。

これら3つのクラスのいずれも、すべての行列に対して最適というわけではありません。あるクラスが他のクラスよりも優れた性能を示す例が常に存在します。そのため、実際には複数のソルバーを試用し、特定の問題に対してどれが最適かを確認します。

最小二乗問題を解く

GMRES 法の 1 つの部分は、を最小化するベクトルを見つけることです。は ( n  + 1) 行n列の行列である ことに注意してください。したがって、これはn個の未知数に対してn +1 個の方程式の制約付き線形システムを生成します。 ${\displaystyle y_{n}}$ $${\displaystyle \left\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\right\|.}$$ ${\displaystyle {\tilde {H}}_{n}}$

最小値はQR分解を用いて計算できます。( n  + 1)行( n  + 1)列の直交行列Ω _nと、( n  + 1)行n列の上三角行列を 求めます。この三角行列は列数より行数が1つ多いため、最下行はゼロになります。したがって、この行列は次のように分解できます。ここで、はn行n列（つまり正方行列）の三角行列です。 ${\displaystyle {\tilde {R}}_{n}}$ $${\displaystyle \Omega _{n}{\tilde {H}}_{n}={\tilde {R}}_{n}.}$$ $${\displaystyle {\tilde {R}}_{n}={\begin{bmatrix}R_{n}\\0\end{bmatrix}},}$$ ${\displaystyle R_{n}}$

QR分解は、ヘッセンベルク行列の差がゼロの行と列のみであるため、反復ごとに簡単に更新できます。 ここで、 h _n+1 = ( h _{1, n +1} , ..., h _{n +1, n +1} ) ^{T です}。これは、ゼロと乗法単位行列の行を追加したΩ _nをヘッセンベルク行列に前置すると、ほぼ三角行列になることを意味します。σ がゼロの場合、これは三角行列になります。これを解決するには、ギブンズ回転が必要です。このギブンズ回転 により、次の式が形成されます 。実際、 はを持つ三角行列です。 $${\displaystyle {\tilde {H}}_{n+1}={\begin{bmatrix}{\tilde {H}}_{n}&h_{n+1}\\0&h_{n+2,n+1}\end{bmatrix}},}$$ $${\displaystyle {\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{n+1}\\0&\rho \\0&\sigma \end{bmatrix}}}$$ $${\displaystyle G_{n}={\begin{bmatrix}I_{n}&0&0\\0&c_{n}&s_{n}\\0&-s_{n}&c_{n}\end{bmatrix}}}$$ $${\displaystyle c_{n}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\text{and}}\quad s_{n}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}.}$$ $${\displaystyle \Omega _{n+1}=G_{n}{\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}.}$$ $${\displaystyle \Omega _{n+1}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{n+1}\\0&r_{n+1,n+1}\\0&0\end{bmatrix}}}$$ ${\textstyle r_{n+1,n+1}={\sqrt {\rho ^{2}+\sigma ^{2}}}}$

QR分解が与えられれば、最小化問題は次のように簡単に解くことができます。ベクトルをg _n ∈ R ⁿおよびγ _n ∈ Rで 表すと、次のようになります。 この式を最小化するベクトルyは次のように与えられます 。ここでも、ベクトルは簡単に更新できます。^[8] $${\displaystyle {\begin{aligned}\left\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\right\|&=\left\|\Omega _{n}({\tilde {H}}_{n}y_{n}-\beta e_{1})\right\|\\&=\left\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\right\|.\end{aligned}}}$$ ${\displaystyle \beta \Omega _{n}e_{1}}$ $${\displaystyle {\tilde {g}}_{n}={\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}}$$ $${\displaystyle {\begin{aligned}\left\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\right\|&=\left\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\right\|\\&=\left\|{\begin{bmatrix}R_{n}\\0\end{bmatrix}}y_{n}-{\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}\right\|.\end{aligned}}}$$ $${\displaystyle y_{n}=R_{n}^{-1}g_{n}.}$$ ${\displaystyle g_{n}}$

サンプルコード

通常のGMRES（MATLAB / GNU Octave）

関数[x, e] = gmres ( A, b, x, max_iterations, しきい値) n = length ( A ); m = max_iterations ;        % x を初期ベクトルとして使用しますr = b - A * x ;        beta = norm ( r ); % 相対残差が追跡され、収束チェックのために入力しきい値と比較されますb_norm = norm ( b );       % 1次元ベクトルを初期化するsn = zeros ( m , 1 ); cs = zeros ( m , 1 ); g = zeros ( m + 1 , 1 ); g ( 1 ) = beta ; Q (:, 1 ) = r / beta ; k = 1 : mの場合                         % arnoldi を実行します[ H ( 1 : k + 1 , k ), Q (:, k + 1 )] = arnoldi ( A , Q , k ); % H 行目の最後の要素を削除し (したがって R をその場で計算)、回転行列を更新します[ H ( 1 : k + 1 , k ), cs ( k ), sn ( k )] = apply_givens_rotation ( H ( 1 : k + 1 , k ), cs , sn , k ); % ギブンズ回転を適用して、新しく拡張された g を計算しますg ( k + 1 ) = - sn ( k ) * g ( k ); g ( k ) = cs ( k ) * g ( k );                                  % 最小値では、||r_k|| = |g(k+1)| current_relative_tol = abs ( g ( k + 1 )) / b_norm ; if ( current_relative_tol <= threshold ) break ; end end % しきい値に達していない場合は、この時点で k = m (m+1 ではない) % 結果を計算しますy = H ( 1 : k , 1 : k ) \ g ( 1 : k ); x = x + Q (:, 1 : k ) * y ; end                               %----------------------------------------------------% % アーノルディ関数 % %----------------------------------------------------% function [h, q] = arnoldi ( A, Q, k ) q = A * Q (:, k ); % クリロフベクトルfor i = 1 : k % ヘッセンベルグ行列を維持した修正グラム・シュミットh ( i ) = q ' * Q (:, i ); q = q - h ( i ) * Q (:, i ); end h ( k + 1 ) = norm ( q ); q = q / h ( k + 1 ); end                                     %----------------------------------------------------------------------% % H 列にギブンズ回転を適用 % %---------------------------------------------------------------------% function [h, cs_k, sn_k] = apply_givens_rotation ( h, cs, sn, k ) % i 番目の列に適用for i = 1 : k - 1 temp = cs ( i ) * h ( i ) + sn ( i ) * h ( i + 1 ) ; h ( i + 1 ) = - sn ( i ) * h ( i ) + cs ( i ) * h ( i + 1 ) ; h ( i ) = temp ; end                                 % 回転の次の sin cos 値を更新します[ cs_k , sn_k ] = givens_rotation ( h ( k ), h ( k + 1 ) );        % H(i + 1, i) を消去h ( k ) = cs_k * h ( k ) + sn_k * h ( k + 1 ); h ( k + 1 ) = 0.0 ;終了                %%----ギブンズ回転行列を計算します----%% function [cs, sn] = givens_rotation ( v1, v2 ) % if (v1 == 0) % cs = 0; % sn = 1; % else t = sqrt ( v1 ^ 2 + v2 ^ 2 ); % cs = abs(v1) / t; % sn = cs * v2 / v1; cs = v1 / t ; % see http://www.netlib.org/eispack/comqr.f sn = v2 / t ; % end end                 

参照

• 共役勾配法

参考文献

1. Saad, Youcef; Schultz, Martin H. (1986). 「GMRES: 非対称線形システムを解くための一般化最小残差アルゴリズム」 . SIAM Journal on Scientific and Statistical Computing . 7 (3): 856– 869. doi :10.1137/0907058. ISSN  0196-5204.
2. ペイジとサンダース、「スパース不定値線形方程式の解」、SIAM J. Numer. Anal.、第12巻、617ページ（1975年）https://doi.org/10.1137/0712047
3. ニファ、ナウファル (2017). Solveurs Performants pour l'optimisation sous contraintes enification de paramètres [パラメータ同定問題における制約付き最適化のための効率的なソルバー] (論文) (フランス語)。
4. Eisenstat, Elman & Schultz 1983, Thm 3.3. 注：GCRのすべての結果はGMRESにも当てはまります。Saad & Schultz 1986を参照。
5. Trefethen, Lloyd N.; Bau, David, III. (1997).数値線形代数. フィラデルフィア: 産業応用数学協会. 定理35.2. ISBN 978-0-89871-361-9。{{cite book}}: CS1 maint: multiple names: authors list (link)
6. Amritkar, Amit; de Sturler, Eric; Świrydowicz, Katarzyna; Tafti, Danesh; Ahuja, Kapil (2015). 「CFDアプリケーションのためのKrylov部分空間のリサイクルと新しいハイブリッドリサイクルソルバー」. Journal of Computational Physics . 303 : 222. arXiv : 1501.03358 . Bibcode :2015JCoPh.303..222A. doi :10.1016/j.jcp.2015.09.040. S2CID  2933274.
7. Gaul, André (2014).線形システムのシーケンスに対するKrylov部分空間法のリサイクル（Ph.D.）. ベルリン工科大学. doi :10.14279/depositonce-4147.
8. Stoer, Josef; Bulirsch, Roland (2002).数値解析入門. 応用数学テキスト（第3版）. ニューヨーク: Springer. §8.7.2. ISBN 978-0-387-95452-3。

• マイスター、アンドレアス。ヴェメル、クリストフ (2005)。Numeric リニア Gleichungssysteme。ヴィースバーデン: Vieweg。ISBN 978-3-528-13135-7。
• Saad, Y. (2003).疎線形システムのための反復法（第2版）. フィラデルフィア: SIAM. ISBN 978-0-89871-534-7。
• アイゼンスタット, スタンレー C.; エルマン, ハワード C.; シュルツ, マーティン H. (1983). 「非対称線形方程式系に対する変分反復法」. SIAM Journal on Numerical Analysis . 20 (2): 345– 357. doi :10.1137/0720023. ISSN  0036-1429.
• ドンガラ他著『線形システムの解法テンプレート：反復法の構成要素』第2版、SIAM、フィラデルフィア、1994年
• Imankulov, Timur; Lebedev, Danil; Matkerim, Bazargul; Daribayev, Beimbet; Kassymbek, Nurislam (2021-10-08). 「多孔質媒体における多相多成分流の数値シミュレーション：ニュートン法に基づく効率解析」. Fluids . 6 (10): 355. doi : 10.3390/fluids6100355 . ISSN  2311-5521.

Retrieved from "https://en.wikipedia.org/w/index.php?title=Generalized_minimal_residual_method&oldid=1318400507"