conjgrad

引数

A

指令した各xについてA*xを計算する 行列または関数またはリスト. Aのそれぞれの型に関する A*x の計算に関して以下に示します.

• 行列.Aが行列の場合, 通常または疎行列が使用可能

• 関数.Aが関数の場合, 以下のヘッダを有する 必要があります :

  function y=A(x)

• リスト.Aがリストの場合, リストの最初の要素には関数を指定し, 添字2から最後までのリストの他の要素には関数の引数を指定します. この関数がコールされた際, xのカレントの値が関数の最初の引数に指定され, その他の引数にはリストで指定されたものが指定されます.

b

右辺ベクトル (大きさ: nx1)

mehtod

スカラー文字列, "pcg", "cgs", "bicg" または "bicgstab" (デフォルト "bicgstab")

tol

相対許容誤差 (デフォルト: 1e-8). 終了基準は残差 r=b-Ax の二乗ノルムを右辺 b の二乗ノルムで割ったものに 基づきます.

maxIter

最大反復回数 (デフォルト: n)

M

プリコンディショナ: 通常または疎行列または M\x を返す関数 (デフォルト: none)

M2

プリコンディショナ: 通常または疎行列または各xの M2\x を返す関数 (デフォルト: none)

x0

初期推定ベクトル (デフォルト: zeros(n,1))

verbose

冗長なログを有効にする場合は1を指定 (デフォルト 0)

x

解ベクトル

flag

conjgrad が反復回数maxi以内に 許容誤差以内に収束した場合は 0, それ以外は 1

err

残差ノルムの最終値 (右辺 b の二乗ノルムを使用)

iter

実行した反復回数

res

相対ノルム残差のベクトル

説明

プリコンディショニング有りまたは無しの 共役勾配法により線形システムAx=b を解きます. プリコンディショナは対称正定行列M,または M=M1*M2となる2つの行列 M1とM2 により定義されます. この場合,この関数はinv(M)*A*x = inv(M)*bを xについて解きます. M, M1 および M2 は, 呼び出し手順y=Milx(x) で対応する左除算y=Mi\xを計算するScilab関数と することができます.

入力引数 method は,使用するソルバを指定します:

• "pcg" プリコンディショナ付き共役勾配法
• "cgs" プリコンディショナ付き二乗共役勾配法
• "bicg" プリコンディショナ付き双共役勾配法
• "bicgstab" プリコンディショナ付き安定化双共役勾配法 (デフォルト)

method="pcg"の場合, A 行列は 対称正定行列(通常または疎行列)または呼び出し手順y=Ax(x)で y=A*xを計算する関数とする必要があります. その他の場合 (method="cgs", "bicg" or "bicgstab"), A は正方行列であることのみ必要です.

良条件および悪条件の問題の例

以下の例では, 2つの線形システムを解きます. 最初の行列は条件数が ~0.02で,アルゴリズムは10回で収束します. これは行列の大きさと同じで,共役勾配法で期待される動作です. 2番目の行列は条件数が1.d-6と小さく,アルゴリズムは収束までに22回と より多くの反復を要します.これがパラメータ maxIter が 30 に設定されている 理由です. "key=value" 構文の他の例については以下を参照ください.

// 良条件問題
A = [ 94  0   0   0    0   28  0   0   32  0
     0   59  13  5    0   0   0   10  0   0
     0   13  72  34   2   0   0   0   0   65
     0   5   34  114  0   0   0   0   0   55
     0   0   2   0    70  0   28  32  12  0
     28  0   0   0    0   87  20  0   33  0
     0   0   0   0    28  20  71  39  0   0
     0   10  0   0    32  0   39  46  8   0
     32  0   0   0    12  33  0   8   82  11
     0   0   65  55   0   0   0   0   11  100];
b = ones(10, 1);
[x, fail, err, iter, res] = conjgrad(A, b, "bicg", 1d-12, 15);
mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)
// 悪条件問題
A = [ 894     0   0     0   0   28  0   0   1000  70000
      0      5   13    5   0   0   0   0   0     0
      0      13  72    34  0   0   0   0   0     6500
      0      5   34    1   0   0   0   0   0     55
      0      0   0     0   70  0   28  32  12    0
      28     0   0     0   0   87  20  0   33    0
      0      0   0     0   28  20  71  39  0     0
      0      0   0     0   32  0   39  46  8     0
      1000   0   0     0   12  33  0   8   82    11
      70000  0   6500  55  0   0   0   0   11    100];
[x, fail, err, iter, res] = conjgrad(A, b, method="pcg", maxIter=30, tol=1d-12);
mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)

Aに疎行列または関数またはリストを指定する例

以下の例は,この手法が疎行列を同様に処理できることを示すものです. また,この例は右辺を計算する関数を"conjgrad"のプリミティブに指定する ケースも示します. この例の最後のケースは, リストをプリミティブに指定すた場合です.

// 良条件問題
A = [ 94  0   0   0    0   28  0   0   32  0
     0   59  13  5    0   0   0   10  0   0
     0   13  72  34   2   0   0   0   0   65
     0   5   34  114  0   0   0   0   0   55
     0   0   2   0    70  0   28  32  12  0
     28  0   0   0    0   87  20  0   33  0
     0   0   0   0    28  20  71  39  0   0
     0   10  0   0    32  0   39  46  8   0
     32  0   0   0    12  33  0   8   82  11
     0   0   65  55   0   0   0   0   11  100];
b = ones(10, 1);
// Aを疎行列に変換
Asparse=sparse(A);
[x, fail, err, iter, res] = conjgrad(Asparse, b, "cgs", maxIter=30, tol=1d-12);
mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)
// 右辺を計算する関数を定義
function y=Atimesx(x)
  A = [ 94  0   0   0    0   28  0   0   32  0
       0   59  13  5    0   0   0   10  0   0
       0   13  72  34   2   0   0   0   0   65
       0   5   34  114  0   0   0   0   0   55
       0   0   2   0    70  0   28  32  12  0
       28  0   0   0    0   87  20  0   33  0
       0   0   0   0    28  20  71  39  0   0
       0   10  0   0    32  0   39  46  8   0
       32  0   0   0    12  33  0   8   82  11
       0   0   65  55   0   0   0   0   11  100];
  y = A*x
endfunction
// スクリプトAtimesx をプリミティブに指定
[x, fail, err, iter, res] = conjgrad(Atimesx, b, maxIter=30, tol=1d-12);
mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)
// 右辺を計算する関数を定義
function y=Atimesxbis(x, A)
  y = A*x
endfunction
// プリミティブへのリストを指定
Alist = list(Atimesxbis, Asparse);
[x, fail, err, iter, res] = conjgrad(Alist, b, maxIter=30, tol=1d-12);
mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)

key=value 構文の例

以下の例は"key=value"構文で引数を指定する方法を示すものです. これにより,位置を固定せずに引数を指定でき, つまり, 引数リストにおける順序に依存せずに, 引数を設定できます. 有効なキーはオプション引数の名前,つまり, tol, maxIter, %M, %M2, x0, verbose です. 以下の例では, verbose オプションが maxIterオプションの 前に指定されていることに注意してください. "key=value"構文ではない場合は引数の位置が固定されるため, maxIter は先, verbose は後で使用する必要があります.

// key=value 構文で指定される引数の例
A = [100 1; 1 10];
b = [101; 11];
[xcomputed, flag, err, iter, res] = conjgrad(A, b, verbose=1);
// key=value 構文の場合, 順序は関係ありません
[xcomputed, flag, err, iter, res] = conjgrad(A, b, verbose=1, maxIter=0);

参照

• backslash — (\) 左行列除算.
• qmr — プリコンディショナ付きのQuasi Minimal Residual法
• gmres — Generalized Minimum RESidual 法

参考文献

PCG

"Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra, Eijkhout, Pozo, Romine, and Van der Vorst, SIAM Publications, 1993, ftp netlib2.cs.utk.edu/linalg/templates.ps

"Iterative Methods for Sparse Linear Systems, Second Edition", Saad, SIAM Publications, 2003, ftp ftp.cs.umn.edu/dept/users/saad/PS/all_ps.zip

CGS

"CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear systems" by Peter Sonneveld.

Original article

Article on ACM

Some theory around CGS

BICG

"Numerical Recipes: The Art of Scientific Computing." (third ed.) by William Press, Saul Teukolsky, William Vetterling, Brian Flannery.

http://apps.nrbook.com/empanel/index.html?pg=87

Article on ACM

Some theory around BICG

BICGSTAB

"Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems" by Henk van der Vorst. 339

Original article

Article on ACM

Some theory around BICG

履歴