一般化最小残差法

数学において、一般化最小残差法 (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\|bA(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 {分} }(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は反復回数です。したがって、この方法は、 x _{k を}初期推定値として、たとえばk回の反復の後に再開されることがあります。結果として得られる方法は、GMRES( k ) または Restarted GMRES と呼ばれます。非正定値行列の場合、再開されたサブスペースは以前のサブスペースに近いことが多いため、この方法では収束が停滞する可能性があります。

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

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

別のクラスの方法は、非対称 Lanczos 反復法、特にBiCG 法に基づいています。これらは 3 項の再帰関係を使用しますが、最小残差は達成されないため、これらの方法では残差が単調に減少しません。収束も保証されません。

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

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

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}{\チルダ {H}}_{n}={\チルダ {R}}_{n}.}$$$${\displaystyle {\tilde {R}}_{n}={\begin{bmatrix}R_{n}\\0\end{bmatrix}},}$$${\displaystyle R_{n}}$

QR 分解は、ヘッセンベルク行列の差が 0 の行と 1 列のみであるため、反復ごとに簡単に更新できます。 ここで、h _n+1 = ( h _{1, n +1} , ..., h _{n +1, n +1} ) ^Tです。これは、ヘッセンベルク行列に Ω _{n を}事前に乗算し、0 と乗法単位行列の行を追加すると、ほぼ三角行列になることを意味します。σ がゼロの場合、これは三角行列になります。これを解決するには、ギブンズ回転が必要です 。 このギブンズ回転により、次の式が形成されます。 実際、 は の三角行列です。 $${\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}}$

関数[x, e] = gmres ( A, b, x, max_iterations, しきい値) n = length ( A ); m = max_iterations ; 
    
    

  % x を初期ベクトルとして使用します
r = b - A * x ;        

  b_norm = norm ( b ); error = norm ( r ) / b_norm ;  
      

  % 1D ベクトルを初期化します
sn = zeros ( m , 1 ); cs = zeros ( m , 1 ); %e1 = zeros(n, 1); e1 = zeros ( m + 1 , 1 ); e1 ( 1 ) = 1 ; e = [ error ]; r_norm = norm ( r ); Q (:, 1 ) = r / r_norm ; % 注: これは上記のセクション「方法」のベータ スカラーではなく、% ベータ スカラーに e1 を掛けたものですbeta = r_norm * e1 ; k = 1の場合: m     
     
  
     
    
    
    
      
  
  
      
     

    % arnoldi を実行します
[ H ( 1 : k + 1 , k ), Q (:, k + 1 )] = arnoldi ( A , Q , k ); % H 行目の最後の要素を削除し、回転行列を更新します[ H ( 1 : k + 1 , k ), cs ( k ), sn ( k )] = apply_givens_rotation ( H ( 1 : k + 1 , k ), cs , sn , k ); % 残差ベクトルを更新しますbeta ( k + 1 ) = - sn ( k ) * beta ( k ); beta ( k ) = cs ( k ) * beta ( k ); error = abs ( beta ( k + 1 )) / b_norm ;           
    
    
            
    
    
          
            
                

    % エラーを保存します
e = [ e ; error ];       

    if ( error <= threshold ) break ; end end % しきい値に達していない場合、この時点で k = m になります (m+1 ではありません) % 結果を計算しますy = H ( 1 : k , 1 : k ) \ beta ( 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 ; % http://www.netlib.org/eispack/comqr.f を参照してくださいsn = v2 / t ; % end end 




        


          
        

• 共役勾配法