conjgrad
共役勾配ソルバ
呼び出し手順
[x, flag, err, iter, res] = conjgrad(A, b [, method [, tol [, maxIter [, M [, M2 [, x0 [, verbose]]]]]]]) [x, flag, err, iter, res] = conjgrad(A, b [, method [, key=value,...]])
引数
- A
指令した各
xについてA*xを計算する 行列または関数またはリスト. Aのそれぞれの型に関する A*x の計算に関して以下に示します.行列.Aが行列の場合, 通常または疎行列が使用可能関数.Aが関数の場合, 以下のヘッダを有する 必要があります :function y=A(x)
リスト.Aがリストの場合, リストの最初の要素には関数を指定し, 添字2から最後までのリストの他の要素には関数の引数を指定します. この関数がコールされた際, xのカレントの値が関数の最初の引数に指定され, その他の引数にはリストで指定されたものが指定されます.
- b
右辺ベクトル (大きさ: nx1)
- method
スカラー文字列, "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);
参照
参考文献
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.
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
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
履歴
| バージョン | 記述 |
| 5.5.0 | 導入 |
| Report an issue | ||
| << Linear Equations (Iterative Solvers) | Linear Equations (Iterative Solvers) | gmres >> |