pcg (MATLAB Function 参考文献)

前提条件付き共役傾斜法

表示

x = pcg(A,b)
pcg(A,b,tol)
pcg(A,b,tol,maxit)
pcg(A,b,tol,maxit,M)
pcg(A,b,tol,maxit,M1,M2)
pcg(A,b,tol,maxit,M1,M2,x0)
pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,...)
[x,flag] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,...) 
[x,flag,relres] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,...)
[x,flag,relres,iter] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,...)
[x,flag,relres,iter,resvec] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,...)

詳細

x = pcg(A,b) は、線形方程式 A*x=b を x について、解きます。 n 行 n 列の係数行列 A は、対称で、正定である必要があり、列ベクトル b は、長さ n である必要があります。A は、afun(x) が、A*x を出力するような関数 afun です。

pcg が収束する場合、その結果に関するメッセージが表示されます。pcg が、設定した最大繰り返し回数に達しても収束しない、または、何らかの理由で、停止した場合、ワーニングメッセージが、相対的な残差 norm(b-A*x)/norm(b) と停止した時点での繰り返し回数とを印刷されます。

pcg(A,b,tol) は、方法のトレランスを設定します。tol が、[] の場合、pcg は、デフォルトの 1e-6 を使います。

pcg(A,b,tol,maxit) は、繰り返しの最大回数を設定します。maxit が、[] の場合、pcg は、デフォルトの min(n,20) を使います。

pcg(A,b,tol,maxit,M) と pcg(A,b,tol,maxit,M1,M2) は、対称な正定前提条件子 M、または、M = M1*M2 を使って、システム inv(M)*A*x = inv(M)*b を x について効率的に解きます。M が、[] の場合、pcg は、前提条件子を適用しません。M は、M\x を戻す関数です。

pcg(A,b,tol,maxit,M1,M2,x0) は、初期推定を指定します。x0 が、[] の場合、pcg は、デフォルトのすべての要素がゼロのベクトルを使います。

pcg(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,...) は、パラメータ p1,p2,... を、関数afun(x,p1,p2,...), m1fun(x,p1,p2,...), m2fun(x,p1,p2,...) に渡します。

[x,flag] = pcg(A,b,tol,maxit,M1,M2,x0) は、収束に関するフラグも戻します。

フラグ
収束

0
pcg は、繰り返し回数が、maxit 以内で、希望するトレランス tol に収束しました。

1
pcg は、maxit 回の繰り返しを行いましたが、収束しませんでした。

2
前提条件子 M は、条件数が良くありません。

3
pcg は、前の処理と同じ結果を出力します。

4

pcg の中で計算されたスカラ量の一つが、計算を続けるには、大きすぎるか、または、小さすぎます。

フラグ	収束
`0`	`pcg` は、繰り返し回数が、`maxit` 以内で、希望するトレランス `tol` に収束しました。
`1`	`pcg` は、`maxit` 回の繰り返しを行いましたが、収束しませんでした。
`2`	前提条件子 `M` は、条件数が良くありません。
`3`	`pcg` は、前の処理と同じ結果を出力します。
`4`	`pcg` の中で計算されたスカラ量の一つが、計算を続けるには、大きすぎるか、または、小さすぎます。

flag が、0 でない場合、戻される解 x は、繰り返し全体に渡って計算される最小ノルム残差です。出力引数 flag が設定された場合は、メッセージは表示されません。

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

[x,flag,relres,iter] = pcg(A,b,tol,maxit,M1,M2,x0) は、x が計算された繰り返し回数もも度します。ここで、0 <= iter <= maxit です。

[x,flag,relres,iter,resvec] = pcg(A,b,tol,maxit,M1,M2,x0) は、各繰り返しでの残差ノルム norm(b-A*x0) を含むベクトルも出力します。

例題

例題 1.

A = gallery('wilk',21);  
b = sum(A,2);
tol = 1e-12;  
maxit = 15;  
M = diag([10:-1:1 1 1:10]);
[x,flag,rr,iter,rv] = pcg(A,b,tol,maxit,M);

また、つぎの行列-ベクトル乗算関数を使います。

function y = afun(x,n)
y = [0; 
     x(1:n-1)] + [((n-1)/2:-1:0)'; 
     (1:(n-1)/2)'].*x + [x(2:n); 
     0];

また、つぎの前提条件子後退解法関数も

function y = mfun(r,n)
y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];

pcg への入力になります。

[x1,flag1,rr1,iter1,rv1] = pcg(@afun,b,tol,maxit,@mfun,...
                               [],[],21);

例題 2.

A = delsq(numgrid('C',25));
b = ones(length(A),1);
[x,flag] = pcg(A,b)

pcg は、デフォルトの繰り返し最大回数 20 回以内で、デフォルトのトレランス 1e-6 に収束しませんでしたので、flag は、1 です。

R = cholinc(A,1e-3);
[x2,flag2,relres2,iter2,resvec2] = pcg(A,b,1e-8,10,R',R)

pcgは、1e-3のドロップ許容誤差をもつ不完全コレスキ分解により前提条件が与えられたときに、6回目の繰り返し(iter2の値)で許容誤差1.2e-9(relres2 の値)に収束するので、flag2は0です。resvec2(1) = norm(b)かつ、resvec2(7) = norm(b-A*x2)です。semilogy(0:iter2,res-vec2/norm(b),'-o')を使って、初期推定(繰り返し回数が0)から始まる各繰り返しでの残差の相対誤差をプロットして、pcgの進行状況を見ることができます。

semilogy(0:iter2,resvec2/norm(b),'-o')
xlabel('iteration number')
ylabel('relative residual')

参考

bicg, bicgstab, cgs, cholinc, gmres, lsqr, minres, qmr, symmlq

@ (関数ハンドル), \ (バックシュラッシュ)

参考文献

Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

pbaspect pchip