ドキュメンテーション センター

  • 評価版
  • 製品アップデート

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

gmres

一般化最小残差法 (リスタート付き)

構文

x = gmres(A,b)
gmres(A,b,restart)
gmres(A,b,restart,tol)
gmres(A,b,restart,tol,maxit)
gmres(A,b,restart,tol,maxit,M)
gmres(A,b,restart,tol,maxit,M1,M2)
gmres(A,b,restart,tol,maxit,M1,M2,x0)
[x,flag] = gmres(A,b,...)
[x,flag,relres] = gmres(A,b,...)
[x,flag,relres,iter] = gmres(A,b,...)
[x,flag,relres,iter,resvec] = gmres(A,b,...)

説明

x = gmres(A,b) は線形方程式 A*x = bx について解きます。nn 列の係数行列 A は、正方行列でなければならず、さらに大規模でスパース行列でなければなりません。列ベクトル b の長さは n でなければなりません。A は、afun(x)A*x を返すような関数ハンドル afun になることができます。この構文の場合、関数 gmres はリスタートしません。反復の最大回数は min(n,10) です。

「関数のパラメーター化」では、必要に応じて関数 afun および以下で説明する前提条件関数 mfun に追加のパラメーターを設定する方法が説明されています。

関数 gmres が収束する場合、この影響を示すメッセージが表示されます。関数 gmres が最大反復回数に達したか、何らかの理由で停止したために収束しない場合、相対残差 norm(b-A*x)/norm(b) と、計算が停止または失敗したときの反復回数を示す警告メッセージが返されます。

gmres(A,b,restart) は内側の反復の restart 回ごとに計算をリスタートします。外側の反復最大回数は min(n/restart,10) です。総反復の最大回数は restart*min(n/restart,10) です。restartn または [] の場合、関数 gmres はリスタートしません。総反復回数は min(n,10) です。

gmres(A,b,restart,tol) はこの方法の許容誤差を指定します。tol[] の場合、関数 gmres は既定の設定の 1e-6 を使用します。

gmres(A,b,restart,tol,maxit) は外側の反復最大回数を設定します。すなわち、総反復回数は restart*maxit 以下になります。maxit[] の場合、関数 gmres は既定の設定の min(n/restart,10) を使用します。restartn または [] の場合、総反復回数の最大数は maxit (restart*maxit の代わり) になります。

gmres(A,b,restart,tol,maxit,M) および gmres(A,b,restart,tol,maxit,M1,M2) は前提条件 M または M = M1*M2 を使用し、方程式 inv(M)*A*x = inv(M)*bx について効率的に解きます。M[] の場合、関数 gmres は前提条件子を適用しません。M は、mfun(x)M\x を返すような関数ハンドル mfun です。

gmres(A,b,restart,tol,maxit,M1,M2,x0) は最初の初期推定を指定します。x0[] の場合、関数 gmres は要素がすべてゼロの既定のベクトルを使用します。

[x,flag] = gmres(A,b,...) は以下の収束フラグも返します。

flag = 0

関数 gmres は、maxit 回の外側の反復回数以内で、必要な許容誤差 tol に収束しました。

flag = 1

関数 gmres は、maxit 回反復しましたが収束しませんでした。

flag = 2

前提条件 M は、条件数が良くありません。

flag = 3

関数 gmres は計算を進めていません (連続する 2 回の計算で、結果は厳密に同一でした)。

flag0 でない場合は、返される解 x は、反復全体に渡り計算される最小のノルム残差内にあります。flag 出力が設定されている場合、メッセージは表示されません。

[x,flag,relres] = gmres(A,b,...) は相対残差 norm(b-A*x)/norm(b) も返します。flag0 の場合は、relres <= tol になります。3 番目の出力 relres は、前提条件を使ったシステムの相対残差です。

[x,flag,relres,iter] = gmres(A,b,...)x が計算された時点での外側の反復回数と内側の反復回数の両方も返します。ここで、0 <= iter(1) <= maxit および 0 <= iter(2) <= restart です。

[x,flag,relres,iter,resvec] = gmres(A,b,...) は内側の反復ごとに残差ノルムのベクトルも返します。これらは、前提条件を使ったシステムの残差ノルムです。

行列入力と gmres の使用

A = gallery('wilk',21);  
b = sum(A,2);
tol = 1e-12;  
maxit = 15; 
M1 = diag([10:-1:1 1 1:10]);

x = gmres(A,b,10,tol,maxit,M1);

以下のメッセージが表示されます。

gmres(10) converged at outer iteration 2 (inner iteration 9) to 
a solution with relative residual 3.3e-013

関数ハンドルと gmres の使用

この例では、前述の例 の行列 A に代えて、行列-ベクトル積関数 afun のハンドルを使用し、前提条件 M1 に代えてバックソルブ関数 mfun のハンドルを使用します。この例は、以下を行う関数 run_gmres に含まれています。

  • 関数ハンドル @afun を第 1 引数として使用して、関数 gmres を呼び出す。

  • 関数 afun および関数 mfun を入れ子にした関数として含み、run_gmres のすべての変数を関数 afun および関数 mfun から使用可能にする。

run_gmres のコードを以下に示します。

function x1 = run_gmres
n = 21;
b = afun(ones(n,1));
tol = 1e-12;  maxit = 15; 
x1 = gmres(@afun,b,10,tol,maxit,@mfun);
 
    function y = afun(x)
        y = [0; x(1:n-1)] + ...
              [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x + ...
              [x(2:n); 0];
    end
 
    function y = mfun(r)
        y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
    end
end

次のように入力した場合、

x1 = run_gmres;

MATLAB® ソフトウェアが以下のメッセージを表示します。

gmres(10) converged at outer iteration 2 (inner iteration 10) 
to a solution with relative residual 1.1e-013.

再起動なしの前提条件の使用

この例は、gmresを再起動しない前提条件の使用方法を示します。

  1. 479 行 479 列の非対称実スパース行列 west0479 を読み込みます。

    load west0479;
    A = west0479;
  2. 許容誤差と反復の最大数を設定します。

    tol = 1e-12; maxit = 20;
    
  3. 真の解がすべて 1 のベクトルになるように、b を定義します。

    b = full(sum(A,2));
    [x0,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);

    fl0 は、gmres が要求した反復回数 20 回以内に要求した許容誤差 1e-12 に収束しなかったため、1 となります。gmres が返す最適な近似解は、最後の 1 つです (it0(2) = 20 を参照)。MATLAB は残差履歴を rv0 に格納します。

  4. gmres の挙動をプロットします。

    semilogy(0:maxit,rv0/norm(b),'-o');
    xlabel('Iteration number');
    ylabel('Relative residual');

    このプロットは、解が徐々に収束することを示します。前提条件が結果を改善することがあります。

  5. ilu を使用して前提条件を作成します。これは、A が非対称であるためです。

    [L,U] = ilu(A,struct('type','ilutp','droptol',1e-5));
    Error using ilu
    There is a pivot equal to zero. Consider decreasing 
    the drop tolerance or consider using the 'udiag' option.
    

    メモ: MATLAB は不完全な LU を構成できません。前提条件子では無意味な特異因子が生成されるためです。

  6. エラー メッセージに示されるように、減少させた棄却許容誤差でもう一度試してください。

    [L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
    [x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);

    gmres9.5436e-14 (rr1 の値) に相対残差を導き出すため、fl1 は 0 です。相対残差は、棄却許容誤差 1e-6 で不完全 LU 分解を前提条件とした場合、6 回目の反復で指定の許容誤差 1e-12 より小さくなります (it1(2) の値)。出力 rv1(1)norm(M\b) です (M = L*U)。出力 rv1(7)norm(U\(L\(b-A*x1))) です。

  7. 初期推定 (反復回数が 0) から始まる各反復での残差の相対残差をプロットして、関数 gmres の進行状況を確認します。

    semilogy(0:it1(2),rv1/norm(b),'-o');
    xlabel('Iteration number');
    ylabel('Relative residual');

再起動ありの前提条件の使用

この例は、gmresを再起動する前提条件の使用方法を示します。

  1. 479 行 479 列の非対称実スパース行列 west0479 を読み込みます。

    load west0479;
    A = west0479;
  2. b を定義して真の解がすべてのベクトルになるようにします。

    b = full(sum(A,2));
  3. 前の例と同じように、不完全な LU 前提条件を構築します。

    [L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));

    再起動ありの gmres を使用するメリットは、このメソッドを実行するのに必要なメモリの量を制限できることです。再起動なしの場合、gmres では、Krylov サブスペースの基底を格納するのに maxit ベクトルのストレージが必要となります。また、gmres は、各ステップで以前のベクトルすべてに対して直交させなければなりません。再起動により、使用するワークスペースの量と外側反復ごとに実行する作業量を制限することができます。前提条件 gmres が前述の 6 反復内で収束した場合でも、アルゴリズムでは最大 12 の基底ベクトルが許可されているため、そのスペースがすべて割り当てられます。

  4. gmres(3)gmres(4)、および gmres(5) を実行します。

    tol = 1e-12; maxit = 20;
    re3 = 3;
    [x3,fl3,rr3,it3,rv3] = gmres(A,b,re3,tol,maxit,L,U);
    re4 = 4;
    [x4,fl4,rr4,it4,rv4] = gmres(A,b,re4,tol,maxit,L,U);
    re5 = 5;
    [x5,fl5,rr5,it5,rv5] = gmres(A,b,re5,tol,maxit,L,U);

    fl3fl4、および fl5 はすべて 0 です。その理由はgmres が再起動されるたびに、相対残差が指定の許容誤差 1e-12 未満になるように導き出すからです。以下のプロットは、再起動されたそれぞれの gmres メソッドの収束履歴を示します。gmres(3) は、外側反復 5、内側反復 3 で収束します (it3 = [5, 3])。これは、外側反復 6、内側反復 0 と同じであるため、最終目盛りでマーキング 6 となります。

参照

Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems:Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

Saad, Youcef and Martin H. Schultz, "GMRES:A generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM J. Sci.Stat.Comput., July 1986, Vol. 7, No. 3, pp. 856-869.

参考

| | | | | | | | | | |

この情報は役に立ちましたか?