このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

bicgstab

線形方程式系の求解 — 双共役傾斜安定化法

構文

x = bicgstab(A,b)

x = bicgstab(A,b,tol)

x = bicgstab(A,b,tol,maxit)

x = bicgstab(A,b,tol,maxit,M)

x = bicgstab(A,b,tol,maxit,M1,M2)

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

[x,flag] = bicgstab(___)

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

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

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

説明

x = bicgstab(A,b) は、双共役傾斜安定化法を使用して、線形方程式系 A*x = b を x について解きます。試行が正常に完了した場合、bicgstab は収束を確認するメッセージを表示します。bicgstab が最大反復回数に達しても収束しないか、何らかの理由で実行を停止した場合は、相対残差 norm(b-A*x)/norm(b) と停止時の反復回数を含む診断メッセージが表示されます。

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

x = bicgstab(A,b,tol,maxit) は、使用する最大反復回数を指定します。bicgstab は、maxit 以内の反復で収束しない場合、診断メッセージを表示します。

x = bicgstab(A,b,tol,maxit,M) は前処理行列 M を指定し、y について系 $$A{M}^{-1}y=b$$ を実質的に解くことにより x を計算します。ここで、$$y=Mx$$ です。前処理行列を使用すると、問題の数値的なプロパティと計算の効率を向上させることができます。

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

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

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

[x,flag,relres] = bicgstab(___) は、相対残差 norm(b-A*x)/norm(b) も返します。flag が 0 の場合は、relres <= tol になります。

[x,flag,relres,iter] = bicgstab(___) は、x が計算されたときの反復回数 iter も返します。

[x,flag,relres,iter,resvec] = bicgstab(___) は、2 回の反復ごとに、最初の残差 norm(b-A*x0) を含む残差ノルムのベクトルも返します。

例

線形システムの反復法

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

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

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

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

x = bicgstab(A,b);

bicgstab stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 0.12.

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

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

x = bicgstab(A,b,1e-4,100);

bicgstab 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.051.

許容誤差が緩く、反復回数が多い場合でも、残差誤差はあまり向上しません。反復アルゴリズムがこの方法で停滞する場合、前処理行列が必要であることを示しています。

A の不完全コレスキー分解を計算し、bicgstab への前処理行列入力として L' 因子を使用します。

L = ichol(A);
x = bicgstab(A,b,1e-4,100,L');

bicgstab converged at iteration 30.5 to a solution with relative residual 4.9e-05.

前処理行列を使用すると、bicgstab が収束できるまで問題の数値特性が向上します。

前処理行列を使用した bicgstab の使用

線形システムを解くために前処理行列を bicgstab と使用する効果を調べます。

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

load west0479
A = west0479;

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

• x は A*x = b の計算された解です。

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

• rr0 は計算解 x の相対残差です。

• it0 は x が計算されたときの反復回数です。

• rv0 は $\Vert \mathit{b}-\mathrm{Ax}\Vert$ の残差履歴のベクトルです。

[x,fl0,rr0,it0,rv0] = bicgstab(A,b,tol,maxit); 
fl0

fl0 = 1

rr0

rr0 = 1

it0

it0 = 0

fl0 は、bicgstab が要求した反復回数 20 回以内に要求した許容誤差 1e-12 に収束しなかったため 1 となります。実際、bicgstab の動作が良好でないために初期推定 x0 = zeros(size(A,2),1) が最適解であり、it0 = 0 で示されるとおり、これが返されます。

遅い収束への対応として、前処理行列を指定できます。A は非対称であるため、ilu を使用して前処理行列 $\mathit{M}=\mathit{L}\mathit{U}$ を生成します。棄却許容誤差を指定して、1e-6 よりも小さい値をもつ非対角エントリを無視します。bicgstab への入力として L および U を指定して、前処理された系 ${\mathit{A}\mathit{M}}^{-1}\mathit{M}\mathit{x}=\mathit{b}$ を解きます。

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

fl1 = 0

rr1

rr1 = 5.9892e-14

it1

it1 = 3

ilu 前処理行列を使用すると、3 回目の反復で 1e-12 の指定の許容誤差より少ない相対残差が生成されます。出力 rv1(1) は norm(b)、出力 rv1(end) は norm(b-A*x1) になります。

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

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

初期推定の指定

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

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

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

bicgstab を使用して $\mathrm{Ax}=\mathit{b}$ を 2 回解きます。1 回は既定の初期推定、もう 1 回は解の適切な初期推定を使用します。両方の解に対して 50 回の反復と既定の許容誤差を使用します。すべての要素が 0.99 と等価のベクトルとして初期推定を 2 番目の解に指定します。

maxit = 50;
x1 = bicgstab(A,b,[],maxit);

bicgstab converged at iteration 20.5 to a solution with relative residual 9.3e-07.

x0 = 0.99*e;
x2 = bicgstab(A,b,[],maxit,[],[],x0);

bicgstab converged at iteration 4 to a solution with relative residual 8.7e-07.

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

中間結果を返す

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

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

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

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

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

bicgstab に、係数行列 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 がベクトルを乗算する場合、結果のベクトルのほとんどの要素はゼロとなります。結果の非ゼロ要素は、A の非ゼロの三重対角要素に対応します。さらに、主対角のみに 1 と等しくない非ゼロ要素があります。

式 $\mathrm{Ax}$ は次のようになります。

$\mathrm{Ax}=\left[\begin{array}{cccccccccc}10& 1& 0& \cdots & & \cdots & & \cdots & 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& \ddots & 0\\ 0& & & & & & 0& \ddots & \ddots & 1\\ 0& 0& \cdots & & \cdots & & \cdots & 0& 1& 10\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\\ {\mathit{x}}_3\\ {\mathit{x}}_4\\ {\mathit{x}}_5\\ ⋮\\ \\ ⋮\\ \\ {\mathit{x}}_{21}\end{array}\right]=\left[\begin{array}{c}10{\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+8{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$.

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

$\mathrm{Ax}=\left[\begin{array}{c}0+10{\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+8{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}+0\end{array}\right]$=$\left[\begin{array}{c}0\\ {\mathit{x}}_1\\ ⋮\\ {\mathit{x}}_{20}\end{array}\right]+\left[\begin{array}{c}10{\mathit{x}}_1\\ 9{\mathit{x}}_2\\ ⋮\\ 10{\mathit{x}}_{21}\end{array}\right]+\left[\begin{array}{c}{\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{21}\\ 0\end{array}\right]$

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

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

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

ここで、bicgstab に A*x を計算する関数ハンドルを与えて、線形システム $\mathrm{Ax}=\mathit{b}$ を解きます。1e-12 の許容誤差と 50 回の反復を使用します。

b = ones(21,1);
tol = 1e-12;  
maxit = 50;
x1 = bicgstab(@afun,b,tol,maxit)

bicgstab converged at iteration 11.5 to a solution with relative residual 5.2e-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 — 係数行列

行列 | 関数ハンドル

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

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

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

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

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

b — 線形方程式の右辺

列ベクトル

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

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

tol — メソッドの許容誤差

[] または 1e-6 (既定値) | 正のスカラー

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

データ型: double

maxit — 最大反復回数

[] または min(size(A,1),20) (既定値) | 正のスカラー整数

最大反復回数。正のスカラー整数として指定します。maxit の値を増やして、bicgstab が許容誤差 tol を満たすためにより多くの反復を行えるようにします。一般に、tol の値が小さいほど、計算を正常に完了するための反復がより多く必要であるということを意味します。

M, M1, M2 — 前処理行列 (個別の引数として)

eye(size(A)) (既定値) | 行列 | 関数ハンドル

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

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

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

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

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

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

x0 — 初期推定

[] またはゼロの列ベクトル (既定値) | 列ベクトル

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

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

出力引数

x — 線形システムの解

列ベクトル

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

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

flag — 収束フラグ

スカラー

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

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

relres — 相対残差誤差

スカラー

相対残差誤差。スカラーとして返されます。相対残差誤差 relres = norm(b-A*x)/norm(b) は解の正確さを示しています。計算が許容誤差 tol に反復回数 maxit 以内で収束する場合、relres <= tol になります。

データ型: double

iter — 反復回数

スカラー

反復回数。スカラーとして返されます。この出力は、x に対する計算解が計算された時点での反復回数を示します。bicgstab の外側の反復ごとに内側の反復が 2 つ含まれるため、iter を反復回数の小数として返すことができます。

データ型: double

resvec — 残差誤差

ベクトル

残差誤差。ベクトルとして返されます。残差誤差 norm(b-A*x) によって、アルゴリズムが x の特定の値に対する収束にどの程度近いかがわかります。resvec の要素数は反復回数の半分と同じです。resvec の内容を調べて、tol または maxit の値を変更するかどうかを決定することができます。

データ型: double

詳細

双共役傾斜安定化法

双共役傾斜安定化 (BiCGSTAB) アルゴリズムは、再起動ありの GMRES ステップを使用して BiCG ステップにおける不規則な収束動作を緩和することにより、BiCG アルゴリズムを改善するために開発されました。

BiCGSTAB の GMRES(1) 修復ステップは次数 1 の最小残渣多項式のみを使用するため、これらの低次元ステップが不十分な場合、このメソッドは機能しない場合があります。高次の多項式の使用により、BiCGSTABL などのその他のメソッドが開発されました[1]。

ヒント

• ほとんどの反復メソッドにおける収束は、係数行列の条件数 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] van der Vorst, H.A., "BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems," SIAM J. Sci. Stat. Comput., March 1992, Vol. 13, No. 2, pp. 631–644.

拡張機能

スレッドベースの環境

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

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

GPU 配列

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

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

• 入力 A がスパース行列の場合:

  • 1 つのスパース前処理行列 M のみがサポートされています。

  • 2 つの前処理行列 M1 と M2 を使用する場合、どちらの前処理行列も関数でなければなりません。

  • 出力パラメーター iter と resvec は反復全体のみを記録します。

  • bigcstab は停滞 (フラグ 3) を検出しません。代わりに、収束の失敗 (フラグ 1) をレポートします。

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

分散配列

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

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

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

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

バージョン履歴

R2006a より前に導入

参考

bicg | bicgstabl | cgs | gmres | lsqr | ilu | minres | pcg | qmr | symmlq | mldivide

トピック

• 線形方程式の反復法