gmres

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

ページ内をすべて折りたたむ

構文

x = gmres(A,b)

x = gmres(A,b,restart)

x = gmres(A,b,restart,tol)

x = gmres(A,b,restart,tol,maxit)

x = gmres(A,b,restart,tol,maxit,M)

x = gmres(A,b,restart,tol,maxit,M1,M2)

x = gmres(A,b,restart,tol,maxit,M1,M2,x0)

[x,flag] = gmres(___)

[x,flag,relres] = gmres(___)

[x,flag,relres,iter] = gmres(___)

[x,flag,relres,iter,resvec] = gmres(___)

説明

x = gmres(A,b) は、一般化最小残差法を使用して、連立線形方程式 A*x = b を x について解きます。試行が正常に完了した場合、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 です。restart が size(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) を使用します。restart が size(A,1) または [] の場合、最大合計反復回数は maxit (restart*maxit ではない) となります。gmres は、最大合計反復回数内で収束しない場合、診断メッセージを表示します。

例

x = gmres(A,b,restart,tol,maxit,M) は前処理行列 M を指定し、方程式 $M^{- 1} A x = M^{- 1} b$ を実質的に解くことにより 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 が含まれます。flag が 0 の場合は、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 = sprandn(400,400,0.5) + 12*speye(400);
b = rand(400,1);

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

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.35.

既定では 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.0045.

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

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

restart = 100;
maxit = 20;
x = gmres(A,b,restart,tol,maxit);

gmres(100) converged at outer iteration 2 (inner iteration 75) to a solution with relative residual 9.3e-05.

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

再起動なしの `gmres` での前処理行列の使用

ライブスクリプトを開く

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

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

load west0479
A = west0479;

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

x は A*x = b の計算された解です。
fl0 はアルゴリズムが収束したかどうかを示すフラグです。
rr0 は計算解 x の相対残差です。
it0 は、x が計算されたときの内側と外側の反復回数を示す 2 要素ベクトル [outer inner] です。
rv0 は $‖ b - Ax ‖$ の残差履歴のベクトルです。

[x,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 = L U$ を生成します。棄却許容誤差を指定して、1e-6 よりも小さい値をもつ非対角エントリを無視します。gmres への入力として L および U を指定して、前処理された方程式 $M^{- 1} A x = M^{- 1} b$ を解きます。前処理行列を指定すると、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 = 
1.3557e-13

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')

Figure contains an axes object. The axes object with xlabel Iteration number, ylabel Relative residual contains 3 objects of type line, constantline. These objects represent No preconditioner, ILU preconditioner, Tolerance.

再起動ありの `gmres` での前処理行列の使用

ライブスクリプトを開く

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

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

load west0479
A = west0479;

$Ax = b$ に対する真の解がすべて 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

fl3、fl4 および 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

Figure contains an axes object. The axes object with title gmres(N) for N = 3, 4, 5, xlabel Outer iteration number, ylabel Relative residual contains 4 objects of type line, constantline. These objects represent gmres(3), gmres(4), gmres(5), Tolerance.

一般的に、内側の反復回数が増えると、外側の反復ごとに 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 と等価のベクトルとして初期推定を 2 番目の解に指定します。

maxit = 200;
x1 = gmres(A,b,[],[],maxit);

gmres converged at iteration 27 to a solution with relative residual 9.5e-07.

x0 = 0.99*e;
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
     0     0     0     0     0     0     0     0     0     1     0     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     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     3     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
      ⋮

ウィルキンソン行列の構造は特殊なため、演算 A*x を関数ハンドルで表すことができます。A がベクトルを乗算する場合、結果のベクトルのほとんどの要素はゼロとなります。結果の非ゼロ要素は、A の非ゼロの三重対角要素に対応します。さらに、主対角のみに 1 と等しくない非ゼロ要素があります。

式 $Ax$ は次のようになります。

$Ax = [\begin{array}{c} 10 & 1 & 0 & \dots & \dots & \dots & 0 & 0 \\ 1 & 9 & 1 & 0 & 0 \\ 0 & 1 & 8 & 1 & 0 & ⋮ \\ ⋮ & 0 & 1 & 7 & 1 & 0 \\ 0 & 1 & 6 & 1 & 0 & ⋮ \\ ⋮ & 0 & 1 & 5 & 1 & 0 \\ 0 & 1 & 4 & 1 & 0 & ⋮ \\ ⋮ & 0 & 1 & 3 & ⋱ & 0 \\ 0 & 0 & ⋱ & ⋱ & 1 \\ 0 & 0 & \dots & \dots & \dots & 0 & 1 & 10 \end{array}] [\begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5} \\ ⋮ \\ ⋮ \\ x_{21} \end{array}] = [\begin{array}{c} 10 x_{1} + x_{2} \\ x_{1} + 9 x_{2} + x_{3} \\ x_{2} + 8 x_{3} + x_{4} \\ ⋮ \\ x_{19} + 9 x_{20} + x_{21} \\ x_{20} + 10 x_{21} \end{array}]$ .

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

$Ax = [\begin{array}{c} 0 + 10 x_{1} + x_{2} \\ x_{1} + 9 x_{2} + x_{3} \\ x_{2} + 8 x_{3} + x_{4} \\ ⋮ \\ x_{19} + 9 x_{20} + x_{21} \\ x_{20} + 10 x_{21} + 0 \end{array}]$ = $[\begin{array}{c} 0 \\ x_{1} \\ ⋮ \\ x_{20} \end{array}] + [\begin{array}{c} 10 x_{1} \\ 9 x_{2} \\ ⋮ \\ 10 x_{21} \end{array}] + [\begin{array}{c} x_{2} \\ ⋮ \\ x_{21} \\ 0 \end{array}]$

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

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

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

ここで、gmres に A*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.

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

afun(x1)

ローカル関数

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

入力引数

すべて折りたたむ

`A` — 係数行列
行列 | 関数ハンドル

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

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

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

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

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

`b` — 線形方程式の右辺
列ベクトル

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

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

`restart` — 再起動前の内側の反復の回数
`[]` または `size(A,1)` (既定値) | 正のスカラー整数

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

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

データ型: double

`tol` — メソッドの許容誤差
`[]` または `1e-6` (既定値) | 正のスカラー

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

データ型: single | double

`maxit` — 外側の最大反復回数
正のスカラー整数

外側の最大反復回数。正のスカラー整数として指定します。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) になります。

データ型: single | double

`M`, `M1`, `M2` — 前処理行列 (個別の引数として)
`eye(size(A))` (既定値) | 行列 | 関数ハンドル

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

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

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

オプションで、M、M1 または M2 のいずれかを、行列ではなく関数ハンドルとして指定できます。関数ハンドルは、前処理行列全体を作成するのではなく行列ベクトル演算を行い、計算をより効率化します。

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

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

`x0` — 初期推定
`[]` またはゼロの列ベクトル (既定値) | 列ベクトル

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

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

出力引数

すべて折りたたむ

`x` — 線形システムの解
列ベクトル

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

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

データ型: single | double

`flag` — 収束フラグ
スカラー

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

フラグの値	収束
`0`	成功 — `gmres` は希望する許容誤差 `tol` に反復回数 `maxit` 以内で収束しました。
`1`	失敗 — `gmres` は `maxit` 回の反復を行いましたが収束しませんでした。
`2`	失敗 — 前処理行列 `M` または `M = M1*M2` が悪条件になっています。
`3`	失敗 — 連続する 2 回の反復が同じであり、`gmres` は計算を進めていません。
`4`	失敗 — `gmres` アルゴリズムで計算されたスカラー量の 1 つが、計算を続けるには大きすぎるかまたは小さすぎます。

`relres` — 相対残差誤差
スカラー

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

`iter` — 外側と内側の反復回数
ベクトル

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

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

データ型: double

`resvec` — 残差誤差
ベクトル

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

詳細

すべて折りたたむ

一般化最小残差法

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

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

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

内側の反復と外側の反復

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

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

たとえば、restart を指定しないと、最大反復回数は maxit の値によって決定され、gmres は再起動しません。

If the restart argument is not specified and the value of maxit is 4, then gmres performs a total of 4 iterations.

ただし、restart を指定すると、関数 gmres は、外側の反復 (maxit で指定) ごとに内側の反復 (restart で指定) を数回実行します。この場合、最大合計反復回数は restart*maxit です。

If the restart argument is specified as 4, and the maxit argument is specified as 2, then gmres performs 4 inner iterations for each outer iteration for a total of 8 iterations.

ヒント

ほとんどの反復メソッドにおける収束は、係数行列の条件数 cond(A) に依存します。A が正方行列の場合、equilibrate を使用して条件数を改善することができ、それ自体でほとんどの反復ソルバーが収束しやすくなります。ただし、equilibrate を使用することでも、平衡化した行列 B = R*P*A*C を次に因子分解する際に前処理行列の品質が向上します。
dissect や symrcm などの行列の並べ替え関数を使用して係数行列の行と列を並べ替え、係数行列が因子分解されて前処理行列が生成される際に非ゼロの数を最小化できます。これによって、後で前処理を使用した線形システムを解くのに必要なメモリと時間を削減できます。

参照

[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.

拡張機能

すべて展開する

スレッドベースの環境
MATLAB® の `backgroundPool` を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の `ThreadPool` を使用してコードを高速化します。

この関数はスレッドベースの環境を完全にサポートしています。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

GPU 配列
Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

gmres 関数は GPU 配列入力をサポートしますが、次の使用上の注意および制限があります。

入力 A がスパース行列の場合:
- スパース前処理行列では、resvec は前処理行列を適用せずに残差ノルム norm(B-A*X) を返します。
- M1 と M2 の 2 つの前処理行列を使用する場合、それらは下三角行列と上三角行列であるか、または両方とも関数ハンドルである必要があります。関数ハンドルの代わりに下三角および上三角の前処理行列を使用すると、計算速度が大幅に向上します。
- GPU 配列では、gmres は停滞 (フラグ 3) を検出しません。代わりに、収束の失敗 (フラグ 1) をレポートします。

詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

分散配列
Parallel Computing Toolbox™ を使用して、クラスターの結合メモリ上で大きなアレイを分割します。

使用上の注意および制限:

M1 が関数である場合は、各行に個別に適用されます。

詳細については、分散配列を使用した MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

すべて展開する

R2025a: 単精度の入力引数の指定

以下の引数は単精度として指定できます。

A — 係数行列
b — 線形方程式の右辺
M、M1、M2 — 前処理行列
x0 — 初期推定

これらの引数のいずれかを単精度として指定すると、関数は単精度で計算を行い、線形システムの解、相対残差誤差、および残差誤差の出力を single 型として返します。計算速度を向上させるには、関数ハンドル出力を含むすべての引数を同じ精度として指定します。

参考

lsqr | minres | pcg | qmr | symmlq | mldivide | ilu

トピック

線形方程式の反復法

gmres

構文

説明

例

線形システムの反復法

再起動なしの gmres での前処理行列の使用

再起動ありの gmres での前処理行列の使用

初期推定の指定

数値行列に代わる関数ハンドルの使用

入力引数

A — 係数行列 行列 | 関数ハンドル

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

b — 線形方程式の右辺 列ベクトル

restart — 再起動前の内側の反復の回数 [] または size(A,1) (既定値) | 正のスカラー整数

tol — メソッドの許容誤差 [] または 1e-6 (既定値) | 正のスカラー

maxit — 外側の最大反復回数 正のスカラー整数

M, M1, M2 — 前処理行列 (個別の引数として) eye(size(A)) (既定値) | 行列 | 関数ハンドル

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

x0 — 初期推定 [] またはゼロの列ベクトル (既定値) | 列ベクトル

出力引数

x — 線形システムの解 列ベクトル

flag — 収束フラグ スカラー

relres — 相対残差誤差 スカラー

iter — 外側と内側の反復回数 ベクトル

resvec — 残差誤差 ベクトル

詳細

一般化最小残差法

内側の反復と外側の反復

ヒント

参照

拡張機能

スレッドベースの環境 MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

GPU 配列 Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

分散配列 Parallel Computing Toolbox™ を使用して、クラスターの結合メモリ上で大きなアレイを分割します。

バージョン履歴

R2025a: 単精度の入力引数の指定

参考

トピック

再起動なしの `gmres` での前処理行列の使用

再起動ありの `gmres` での前処理行列の使用

`A` — 係数行列
行列 | 関数ハンドル

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

`b` — 線形方程式の右辺
列ベクトル

`restart` — 再起動前の内側の反復の回数
`[]` または `size(A,1)` (既定値) | 正のスカラー整数

`tol` — メソッドの許容誤差
`[]` または `1e-6` (既定値) | 正のスカラー

`maxit` — 外側の最大反復回数
正のスカラー整数

`M`, `M1`, `M2` — 前処理行列 (個別の引数として)
`eye(size(A))` (既定値) | 行列 | 関数ハンドル

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

`x0` — 初期推定
`[]` またはゼロの列ベクトル (既定値) | 列ベクトル

`x` — 線形システムの解
列ベクトル

`flag` — 収束フラグ
スカラー

`relres` — 相対残差誤差
スカラー

`iter` — 外側と内側の反復回数
ベクトル

`resvec` — 残差誤差
ベクトル

スレッドベースの環境
MATLAB® の `backgroundPool` を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の `ThreadPool` を使用してコードを高速化します。

GPU 配列
Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

分散配列
Parallel Computing Toolbox™ を使用して、クラスターの結合メモリ上で大きなアレイを分割します。