Main Content

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

gmres

線形方程式系の求解 — 一般化最小残差法

説明

x = gmres(A,b) は、一般化最小残差法を使用して、線形方程式系 A*x = bx について解きます。試行が正常に完了した場合、gmres は収束を確認するメッセージを表示します。gmres が最大反復回数に達しても収束しないか、何らかの理由で実行を停止した場合は、相対残差 norm(b-A*x)/norm(b) と停止時の反復回数を含む診断メッセージが表示されます。この構文の場合、関数 gmres は再起動しません。最大反復回数は min(size(A,1),10) です。

x = gmres(A,b,restart)内側の反復restart 回ごとに計算を再起動します。外側の最大反復回数は outer = min(size(A,1)/restart,10) です。gmres は外側の反復ごとに restart 回内側の反復を実行するため、最大合計反復回数は restart*outer です。restartsize(A,1) または [] の場合、関数 gmres は再起動しません。最大反復回数の合計は min(size(A,1),10) です。

x = gmres(A,b,restart,tol) は、このメソッドの許容誤差を指定します。既定の許容誤差は 1e-6 です。

x = gmres(A,b,restart,tol,maxit)外側の最大反復回数の合計を指定します。すなわち、合計反復回数は restart*maxit 以下になります。maxit[] の場合、関数 gmres は既定の設定の min(size(A,1)/restart,10) を使用します。restartsize(A,1) または [] の場合、最大合計反復回数は maxit (restart*maxit ではない) となります。gmres は、最大合計反復回数内で収束しない場合、診断メッセージを表示します。

x = gmres(A,b,restart,tol,maxit,M) は前処理行列 M を指定し、システム M1Ax=M1b を x について効果的に解くことにより x を計算します。前処理行列を使用すると、問題の数値的なプロパティと計算の効率を向上させることができます。

x = gmres(A,b,restart,tol,maxit,M1,M2)M = M1*M2 となるような前処理行列 M の因子を指定します。

x = gmres(A,b,restart,tol,maxit,M1,M2,x0) は解のベクトル x の初期推定を指定します。既定値はゼロのベクトルです。

[x,flag] = gmres(___) は、アルゴリズムが正常に収束したかどうかを示すフラグを返します。flag = 0 の場合、収束は正常に実行されています。この出力構文は、前記のすべての入力引数の組み合わせで使用できます。flag 出力を指定する場合、gmres は診断メッセージを表示しません。

[x,flag,relres] = gmres(___) は、相対残差 norm(M\(b-A*x))/norm(M\b) も返します。これには前処理行列 M が含まれます。flag0 の場合は、relres <= tol になります。

[x,flag,relres,iter] = gmres(___) は、x がベクトル [outer inner] として計算されたときの内側および外側の反復回数も返します。外側の反復回数は範囲 0 <= iter(1) <= maxit 内にあり、内側の反復回数は範囲 0 <= iter(2) <= restart 内にあります。

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

すべて折りたたむ

既定の設定で gmres を使用して正方線形システムを解いてから、解法プロセスで使用される許容誤差と反復回数を調整します。

密度が 50% の乱数スパース行列 A を作成します。また、Ax=b の右辺に乱数のベクトル b を作成します。

rng default
A = sprand(400,400,.5);
A = A'*A;
b = rand(400,1);

gmres を使用して Ax=b を解きます。出力の表示には相対残差誤差 Ax-bb の値が含まれます。

x = gmres(A,b);
gmres stopped at iteration 10 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 10) has relative residual 0.15.

既定では gmres は 10 回の反復と 1e-6 の許容誤差を使用し、アルゴリズムではこの行列についてはその 10 回の反復で収束できません。残差がまだ大きいため、より多くの反復 (または前処理行列) が必要であることを示しています。許容誤差を減らして、アルゴリズムの収束をより簡単にすることもできます。

1e-4 の許容誤差と 100 回の反復を使用して、再度システムを解きます。

tol = 1e-4;
maxit = 100;
x = gmres(A,b,[],tol,maxit);
gmres stopped at iteration 100 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 100) has relative residual 0.052.

許容誤差が緩く、反復回数が多い場合でも、残差誤差はあまり向上しません。反復アルゴリズムがこの方法で停滞する場合、前処理行列が必要であることを示しています。ただし、gmres には、内側の反復回数を制御する入力もあります。内側の反復回数の値を指定することで、gmres で外側の反復ごとに実行する作業量が増えします。

restart 値に 100 および maxit 値に 250 を使用してシステムをもう一度解きます。100 回の反復を一度に実行するのではなく、gmres は再起動の間に 100 回の反復を実行し、これを 250 回反復します。

restart = 100;
maxit = 250;
x = gmres(A,b,restart,tol,maxit);
gmres(100) converged at outer iteration 134 (inner iteration 39) to a solution with relative residual 0.0001.

この場合に、gmres に対して大きい restart 値を指定すると、許容反復回数内の解に収束できます。ただし、restart 値を大きくすると、A も大きい場合は大量のメモリが使用される可能性があります。

再起動なしの gmres で前処理行列を使用して線形システムを解いた結果を調べます。

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

load west0479
A = west0479;

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

gmres を使用して、要求された許容誤差と反復回数で解を求めます。解法プロセスに関する情報を返す出力を 5 つ指定します。

  • x0A*x0 = b の計算された解です。

  • fl0 はアルゴリズムが収束したかどうかを示すフラグです。

  • rr0 は計算解 x0 の残差です。

  • it0 は、x0 が計算されたときの内側と外側の反復回数を示す 2 要素ベクトル [inner outer] です。

  • rv0Ax-b の残差履歴のベクトルです。

[x0,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 0.7603
it0
it0 = 1×2

     1    20

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

遅い収束への対応として、前処理行列を指定できます。A は非対称であるため、ilu を使用して前処理行列 M=LU を生成します。棄却許容誤差を指定して、1e-6 よりも小さい値をもつ非対角エントリを無視します。gmres への入力として L および U を指定して、前処理されたシステム M-1Ax=M-1b を解きます。前処理行列を指定すると、gmres は出力 rr1 および rv1 の前処理されたシステムの残余ノルムを計算します。

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

     1     6

ilu 前処理行列を使用すると、6 回目の外側の反復で 1e-12 の指定の許容誤差より少ない相対残差が生成されます。最初の残差 rv1(1)norm(U\(L\b)) です。ここでは M = L*U です。最後の残差 rv1(end)norm(U\(L\(b-A*x1))) です。

各反復での相対残差をプロットして、gmres の進行状況を確認できます。指定された許容誤差のラインと共に、それぞれの解の残差履歴をプロットします。

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(U\(L\b)),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

再起動ありの gmres で前処理行列を使用します。

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

load west0479
A = west0479;

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

b = sum(A,2);

1e-6 の棄却許容誤差で不完全な LU 前処理行列を作成します。

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

再起動ありの gmres を使用するメリットは、このメソッドを実行するために必要なメモリの量を制限できることです。再起動なしの gmres では、Krylov サブスペースの基底を格納するために maxit 個のベクトルのストレージが必要です。また、gmres は、各ステップで以前のベクトルすべてに対して直交させなければなりません。再起動により、使用するワークスペースの量と外側反復ごとに実行する作業量を制限することができます。

不完全な LU 因子を前処理行列として使用して、gmres(3)gmres(4)gmres(5) を実行します。1e-12 の許容誤差と、20 回の外側の最大反復を使用します。

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

fl3fl4 および fl5 はすべて 0 です。その理由は gmres が再起動されるたびに、指定の許容誤差 1e-12 未満になる相対残差を導き出すからです。

次のプロットは、各再起動ありの gmres メソッドの収束履歴を示します。gmres(3) は、外側反復 5、内側反復 3 で収束します (it3 = [5, 3])。これは、外側反復 6、内側反復 0 と同じであるため、最後の目盛り上のマーク 6 になります。

semilogy(1:1/3:6,rv3/norm(U\(L\b)),'-o');
h1 = gca;
h1.XTick = (1:6);
title('gmres(N) for N = 3, 4, 5')
xlabel('Outer iteration number');
ylabel('Relative residual');
hold on
semilogy(1:1/4:3,rv4/norm(U\(L\b)),'-o');
semilogy(1:1/5:2.8,rv5/norm(U\(L\b)),'-o');
yline(tol,'r--');
hold off
legend('gmres(3)','gmres(4)','gmres(5)','Tolerance')
grid on

一般的に、内側の反復回数が増えると、外側の反復ごとに gmres が行う作業量が増え、収束も迅速になります。

gmres に解の初期推定を指定する効果を調べます。

三重対角スパース行列を作成します。x の想定される解が 1 のベクトルとなるよう、Ax=b の右辺のベクトルとして各行の合計を使用します。

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

gmres を使用して Ax=b を 2 回解きます。1 回は既定の初期推定、もう 1 回は解の適切な初期推定を使用します。両方の解に対して 200 回の反復を使用し、すべての要素が 0.99 と等価のベクトルとして初期推定を指定します。

maxit = 200;
x1 = gmres(A,b,[],[],maxit);
gmres converged at iteration 27 to a solution with relative residual 9.5e-07.
x0 = 0.99*ones(size(A,2),1);
x2 = gmres(A,b,[],[],maxit,[],[],x0);
gmres converged at iteration 7 to a solution with relative residual 6.7e-07.

この場合、初期推定を指定すると gmres をより迅速に収束させることができます。

中間結果を返す

for ループで gmres を呼び出して、初期推定を使用して中間結果を取得することもできます。ソルバーを呼び出すたびに、数回の反復が行われ、計算された解が格納されます。その後、その解を次の反復のバッチに対する初期ベクトルとして使用します。

たとえば、次のコードは 100 回の反復を 4 回実行し、for ループを通過するたびに、解のベクトルを格納します。

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = gmres(A,b,[],tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

X(:,k) は、for ループの反復 k で計算された解のベクトルであり、R(k) はその解の相対残差です。

gmres に、係数行列 A の代わりに A*x を計算する関数ハンドルを与えて線形システムを解きます。

gallery で生成されたウィルキンソンのテスト行列の 1 つは 21 行 21 列の三重対角行列です。行列をプレビューします。

A = gallery('wilk',21)
A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

ウィルキンソン行列の構造は特殊なため、演算 A*x を関数ハンドルで表すことができます。A の各行は x の要素を乗算するため、いくつかの結果のみが非ゼロです (三重対角行列の非ゼロに対応)。さらに、主対角のみに 1 と等しくない非ゼロ要素があります。式 Ax は次のようになります。

[1010001910001810017100161001510014100130001000110][x1x2x3x4x5x21]=[10x1+x2x1+9x2+x3x2+8x3+x4x19+9x20+x21x20+10x21]

結果のベクトルは、3 つのベクトルの合計として記述できます。

[0+10x1+x2x1+9x2+x3x2+8x3+x4x19+9x20+x21x20+10x21+0]=[0x1x20]+[10x19x210x21]+[x2x210]

MATLAB® で、これらのベクトルを作成および合算することにより A*x の値を与える関数を記述します。

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

(この関数は、ローカル関数として例の最後に保存されています)

ここで、gmresA*x を計算する関数ハンドルを与えて、線形システム Ax=b を解きます。許容誤差 1e-12、15 回の外側の反復、および再起動前の 10 回の内側の反復を使用します。

b = ones(21,1);
tol = 1e-12;  
maxit = 15;
restart = 10;
x1 = gmres(@afun,b,restart,tol,maxit)
gmres(10) converged at outer iteration 5 (inner iteration 10) to a solution with relative residual 5.3e-13.
x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

afun(x1) が 1 のベクトルを生成していることを確認します。

afun(x1)
ans = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

ローカル関数

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

入力引数

すべて折りたたむ

係数行列。正方行列または関数ハンドルとして指定します。この行列は線形システム A*x = b の係数行列です。一般に、A は大規模なスパース行列、または大規模なスパース行列と列ベクトルの積を返す関数ハンドルです。

関数ハンドルとしての A の指定

係数行列を、行列ではなく関数ハンドルとして指定して、計算時のメモリを節約できます。この関数ハンドルは、係数行列全体を作成するのではなく、行列とベクトルの積を返し、計算をより効率的にします。

関数ハンドルを使用するには、関数シグネチャ function y = afun(x) を使用します。関数のパラメーター化では、必要な場合に関数 afun に追加のパラメーターを指定する方法を説明しています。関数呼び出し afun(x)A*x の値を返さなければなりません。

データ型: double | function_handle
複素数のサポート: あり

線形方程式の右辺。列ベクトルとして指定します。b は、size(A,1) と同じ長さの列ベクトルでなければなりません。

データ型: double
複素数のサポート: あり

再起動前の内側の反復回数。スカラー整数として指定します。この入力を maxit 入力と共に使用して、最大反復回数 restart*maxit を制御します。restart[] または size(A,1) の場合、関数 gmres は再起動しません。最大合計反復回数は maxit です。

restart 値が大きいと、通常は収束動作が向上しますが、時間とメモリの要件も高くなります。

データ型: double

メソッドの許容誤差。正のスカラーとして指定します。この入力を使用して計算の精度と実行時間とをトレードオフします。成功するには、gmres は、許容される反復回数内に許容誤差を満たさなければなりません。tol を小さい値にすると、解がより正確でなければ計算が成功しないことを意味します。

データ型: double

外側の最大反復回数。正のスカラー整数として指定します。maxit の値を増やして、gmres が許容誤差 tol を満たすためにより多くの反復を行えるようにします。一般的に、tol の値を小さくすると、計算を正常に実行するために反復回数を増やす必要があります。

restart 入力も指定した場合、合計反復回数は restart*maxit となります。指定しない場合、合計反復回数は maxit です。

maxit の既定値は、restart を指定したかどうかによって異なります。

  • restart を指定しない場合、あるいは [] または size(A,1) として指定した場合、maxit の既定値は min(size(A,1),10) となります。

  • restart を範囲 1 <= restart < size(A,1) 内の値として指定した場合、maxit の既定値は min(ceil(size(A,1)/restart),10) になります。

データ型: double

前処理行列。行列または関数ハンドルの個別の引数として指定します。前処理行列 M またはその行列の因子 M = M1*M2 を指定して線形システムの数値的側面を改善し、gmres での迅速な収束を容易にします。不完全行列因数分解関数 ilu および ichol を使用して前処理行列を生成できます。また、因数分解の前に equilibrate を使用して、係数行列の条件数を改善できます。前処理行列の詳細については、線形方程式の反復法を参照してください。

gmres は指定なし前処理行列を単位行列として扱います。

関数ハンドルとしての M の指定

MM1、または M2 のいずれも、行列ではなく関数ハンドルとして指定して、計算でメモリを節約できます。関数ハンドルは、前処理行列全体を作成するのではなく行列ベクトル演算を行い、計算をより効率化します。

関数ハンドルを使用するには、関数シグネチャ function y = mfun(x) を使用します。関数のパラメーター化では、必要な場合に関数 mfun に追加のパラメーターを指定する方法を説明しています。関数呼び出し mfun(x)M\x または M2\(M1\x) の値を返さなければなりません。

データ型: double | function_handle
複素数のサポート: あり

初期推定。size(A,2) に等しい長さの列ベクトルとして指定します。gmres に、既定値のゼロのベクトルではなく、より妥当な初期推定 x0 を与えることができれば、計算時間が節約でき、アルゴリズムがより速く収束するようにできます。

データ型: double
複素数のサポート: あり

出力引数

すべて折りたたむ

線形システムの解。ベクトルとして返されます。この出力は線形システム A*x = b の近似解を与えます。計算が正常に完了した場合 (flag = 0)、relrestol 以下になります。

計算が正常に完了しない場合 (flag ~= 0) は常に、gmres によって返される解 x は、すべての反復にわたり計算された最小残差ノルムをもちます。

収束フラグ。次の表のいずれかのスカラー値として返されます。収束フラグは、計算が正常に実行されたかどうかを示し、さまざまな失敗の形式を区別します。

フラグの値

収束

0

成功 — gmres は希望する許容誤差 tol に反復回数 maxit 以内で収束しました。

1

失敗 — gmresmaxit 回の反復を行いましたが収束しませんでした。

2

失敗 — 前処理行列 M または M = M1*M2 が悪条件になっています。

3

失敗 — 連続する 2 回の反復が同じであり、gmres は計算を進めていません。

4

失敗 — gmres アルゴリズムで計算されたスカラー量の 1 つが、計算を続けるには大きすぎるかまたは小さすぎます。

相対残差誤差。スカラーとして返されます。相対残差誤差 relres = norm(M\(b-A*x))/norm(M\b) は解の正確さを示しています。gmres では相対残差計算に前処理行列 M を含みますが、他のほとんどの反復ソルバーには含まれません。計算が許容誤差 tol に反復回数 maxit 以内で収束する場合、relres <= tol になります。

データ型: double

外側と内側の反復回数。2 要素ベクトル [outer inner] として返されます。この出力は、x の解が計算されたときの内側と外側の反復回数を示しています。

  • restart が指定されていない場合、あるいは [] または size(A,1) として指定されている場合、outer = 1 とすべての反復は内側の反復と見なされます。

  • restart に範囲 1 <= restart < size(A,1) 内の値が指定されている場合、外側の反復回数は 0 <= outer <= maxit 範囲内となり、内側の反復回数は 0 <= inner <= restart 範囲内となります。

データ型: double

残差誤差。ベクトルとして返されます。残差誤差 norm(M\(b-A*x)) によって、アルゴリズムが x の特定の値に対する収束にどの程度近いかがわかります。gmres では相対残差計算に前処理行列 M を含みますが、他のほとんどの反復ソルバーには含まれません。resvec の要素数は、反復回数の合計と等価です (restart が使用されている場合、これは最大で restart*maxit です)。resvec の内容を調べて、restarttol または maxit の値を変更するかどうかを決定することができます。

データ型: double

詳細

すべて折りたたむ

一般化最小残差法

一般化最小残差 (GMRES) アルゴリズムは、最小残差 (MINRES) アルゴリズムを非対称行列に拡張するために開発されました。

共役勾配 (CG) 手法と同様に、GMRES アルゴリズムは直交シーケンスを計算しますが、GMRES はシーケンス内に以前のベクトルをすべて保存する必要があります。このように以前のベクトルを保存した場合、確認しないでおくと、大量のメモリが使用される可能性があります。アルゴリズムの "再起動あり" のバージョンは、中間シーケンスを定期的にクリアし、結果を別の反復の初期値として使用して、これらのシーケンスの保存を制御できます。

適切なパフォーマンスにするには、適切な "restart" 値を選択することが重要ですが、適切な値の選択はほとんどが経験によります。再起動前の反復回数が少なすぎると、アルゴリズムの収束が非常に遅くなるか、まったく収束できない可能性があります。一方、restart 値が大きいと、アルゴリズムのストレージ要件が高くなり、不必要な作業を行う可能性があります[1]

内側の反復と外側の反復

"内側の反復" とは、再起動前に gmres が実行する反復です。引数 restart を使用して内側の反復回数を指定できます。

gmres が再起動するたびに、"外側の反復" 回数が増えていきます。外側の反復回数は引数 maxit を使用して指定できます。外側の反復回数の既定値は min(size(A,1)/restart,10) です。

最大合計反復回数は restart*maxit です。

ヒント

  • 最大反復のメソッドの収束は、係数行列の条件数 cond(A) に依存します。equilibrate を使用して A の条件数を改善することができ、それ自体で最大反復のソルバーが収束しやすくなります。ただし、equilibrate を使用することでも、平衡化した行列 B = R*P*A*C を次に因子分解する際に前処理行列の品質が向上します。

  • dissectsymrcm などの行列の並べ替え関数を使用して係数行列の行と列を並べ替え、係数行列が因子分解されて前処理行列が生成される際に非ゼロの数を最小化できます。これによって、後で前処理を使用した線形システムを解くのに必要なメモリと時間を削減できます。

参照

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

[2] Saad, Yousef 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 より前に導入