ドキュメンテーション

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

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 に含まれています。

  • 第 1 引数として関数ハンドル @afun を指定し、関数 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 を再起動しない前提条件子の使用方法を示します。

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

load west0479;
A = west0479;

許容誤差と最大反復回数を設定します。

tol = 1e-12;
maxit = 20;

真の解がすべて 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 に格納します。

gmres の動作をプロットします。

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

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

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 を構築できません。これは前提条件子では無意味です。

エラー メッセージに示されるように、棄却許容誤差を小さくしてもう一度試してください。

[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))) です。

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

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

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

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

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

load west0479;
A = west0479;

真の解がすべて 1 のベクトルになるように、b を定義します。

b = full(sum(A,2));

前の例と同じように、不完全な LU 前提条件子を構築します。

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

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

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 となります。

figure
semilogy(1:1/3:6,rv3/norm(b),'-o');
h1 = gca;
h1.XTick = [1:1/3:6];
h1.XTickLabel = ['1';' ';' ';'2';' ';' ';'3';' ';' ';'4';' ';' ';'5';' ';' ';'6';];
title('gmres(3)')
xlabel('Iteration number');
ylabel('Relative residual');

figure
semilogy(1:1/4:3,rv4/norm(b),'-o');
h2 = gca;
h2.XTick = [1:1/4:3];
h2.XTickLabel = ['1';' ';' ';' ';'2';' ';' ';' ';'3'];
title('gmres(4)')
xlabel('Iteration number');
ylabel('Relative residual');

figure
semilogy(1:1/5:2.8,rv5/norm(b),'-o');
h3 = gca;
h3.XTick = [1:1/5:2.8];
h3.XTickLabel = ['1';' ';' ';' ';' ';'2';' ';' ';' ';' '];
title('gmres(5)')
xlabel('Iteration number');
ylabel('Relative residual');

参考文献

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.

参考

| | | | | | | | |

R2006a より前に導入

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