pcg

引数

A

行列, または指定したxについて A*xを計算する関数またはリスト. 以下にAの型に応じたA*xの計算に関する説明を示します.

• 行列.Aが行列の場合, 通常の行列 または疎行列とすることができます

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

  function y=A(x)

• list. Aがリストの場合, リストの最初の要素は関数で,リストのその他の 要素(インデックス2から末尾まで)は関数の引数となります. 関数がコールされる時, xのカレントの値は関数に最初の引数として渡されます. 他の引数はリストに指定された値です.

b

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

tol

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

maxIter

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

M

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

M2

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

x0

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

verbose

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

x

解ベクトル

flag

maxi回の反復回数内に pcgが指定した 許容誤差に収束する場合に 0, そうでない場合に 1

err

最終残差相対ノルム (右辺bの2次ノルムが使用されます)

iter

実行された反復回数

res

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

説明

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

A行列は対称正定行列 (通常の行列または疎行列)または y=A*xを計算する 呼び出し手順y=Ax(x)を有する関数と する必要があります.

良好な条件または悪い条件の問題の例

以下の例では,2つの線形システムを解きます. 最初の行列の条件数は ~0.02 に等しくなり, アルゴリズムはちょうど10回の反復で収束します. これが行列の大きさの場合,共役勾配法で指定される動作となります. 後者は,条件数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]=pcg(A,b,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]=pcg(A,b,maxIter=30,tol=1d-12);
mprintf("      fail=%d, iter=%d, errrel=%e\n",fail,iter,err)

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

以下の例では,疎行列を同時に処理する方法を示します. 右辺を計算する関数が "pcg" プリミティブに指定される 場合も示します. この例における最後のケースでは, リストがプリミティブに指定される場合です.

//良い条件の例
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]=pcg(Asparse,b,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]=pcg(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]=pcg(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. 以下の例では, maxIterオプションの前にverboseオプションを指定します. "key=value" 構文ではなく, 位置を指定する引数の場合は, maxIterを最初にverboseを後に指定する必要があります.

// 引数をkey=value構文で指定する例
A=[100,1;1,10];
b=[101;11];
[xcomputed, flag, err, iter, res]=pcg(A,b,verbose=1);
// key=value構文では, 順番は関係ありません
[xcomputed, flag, err, iter, res]=pcg(A,b,verbose=1,maxIter=0);

参照

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

参考文献

"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

履歴