ドキュメンテーション センター

  • 評価版
  • 製品アップデート

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

pcg

前処理を使った共役勾配法

構文

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)
[x,flag] = pcg(A,b,...)
[x,flag,relres] = pcg(A,b,...)
[x,flag,relres,iter] = pcg(A,b,...)
[x,flag,relres,iter,resvec] = pcg(A,b,...)

説明

x = pcg(A,b) は、線形方程式 A*x=bx について解きます。nn 列の係数行列 A は、対称かつ正定でなければならず、さらに大規模でスパース行列でなければなりません。列ベクトル b の長さは n でなければなりません。A を、afun(x)A*x を返すような関数ハンドル afun になるように指定することもできます。

「関数のパラメーター化」では、必要に応じて関数 afun および以下で説明する前提条件関数 mfun に追加のパラメーターを設定する方法が説明されています。

関数 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)*bx について効率的に解きます。M[] の場合、関数 pcg は前提条件を適用しません。M は、mfun(x)M\x を返すような関数ハンドル mfun です。

pcg(A,b,tol,maxit,M1,M2,x0) は、初期推定を指定します。x0[] の場合、関数 pcg は既定の設定のゼロベクトルを使用します。

[x,flag] = pcg(A,b,...) は、収束に関するフラグも返します。

フラグ

収束

0

関数 pcg maxit の反復の中で、希望する許容誤差 tol に収束しました。

1

関数 pcg は、maxit 回の反復に対して、収束しませんでした。

2

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

3

関数 pcg の出力結果に変化が見られません (連続する 2 回の計算で、結果がまったく同じでした)。

4

関数 pcg の実行中に計算されたスカラー量の 1 つが、計算を続けるには大き過ぎるか、または小さ過ぎます。

flag0 でない場合は、返される解 x は、反復全体に渡り計算される最小のノルム残差内にあります。flag 出力が設定されている場合、メッセージは表示されません。

[x,flag,relres] = pcg(A,b,...) は相対残差の推定 norm(b-A*x)/norm(b) も返します。flag0 の場合は、relres <= tol になります。

[x,flag,relres,iter] = pcg(A,b,...) は、x が計算されたときの反復回数も返します。ここで、0 <= iter <= maxit です。

[x,flag,relres,iter,resvec] = pcg(A,b,...) は、各反復で、norm(b-A*x0) を含めて、残差ノルムのベクトルも返します。

大規模行列と pcg の使用

この例では、行列入力および関数ハンドルで pcg を使用する方法を示します。

n1 = 21; 
A = gallery('moler',n1);  
b1 = sum(A,2);
tol = 1e-6;  
maxit = 15;  
M1 = spdiags((1:n1)',0,n1,n1);
[x1,flag1,rr1,iter1,rv1] = pcg(A,b1,tol,maxit,M1);

行列 A の代わりに、以下の関数を使用することもできます。

function y = applyMoler(x)
y = x;
y(end-1:-1:1) = y(end-1:-1:1) - cumsum(y(end:-1:2));
y(2:end) = y(2:end) - cumsum(y(1:end-1));

この関数を使用すると、全体の行列 A を保存する必要がないため、以下のように大規模システムをより効率的に解くことができます。

n2 = 21;
b2 = applyMoler(ones(n2,1));
tol = 1e-6;
maxit = 15;
M2 = spdiags((1:n2)',0,n2,n2);
[x2,flag2,rr2,iter2,rv2] = pcg(@applyMoler,b2,tol,maxit,M2);

前提条件子と pcg の使用

この例では、pcg で前処理行列を使用する方法を示します。

  1. 以下のように入力行列を作成し、pcg を使用してシステムを解きます。

    A = delsq(numgrid('S',100));
    b = ones(size(A,1),1);
    [x0,fl0,rr0,it0,rv0] = pcg(A,b,1e-8,100);

    関数 pcg は、要求された最大反復回数 100 回以内で、要求された許容誤差 1e-8 に収束しなかったため、fl0 は、1 です。前提条件子を使用することにより、システムをより迅速に収束させることができます。

  2. ゼロ充てんの不完全コレスキー分解を作成するには、以下のように入力引数を 1 つだけもつ ichol を使用します。

    L = ichol(A);
    [x1,fl1,rr1,it1,rv1] = pcg(A,b,1e-8,100,L,L’);

    ゼロ充てんの不完全コレスキー分解による前提条件前の場合、pcg が、77 回目の反復 (it1 の値) で要求許容誤差 1e-8 未満となる相対残差 9.8e-09 (rr1 の値) を導き出すので、fl10 です。rv1(1) = norm(b)rv1(78) = norm(b-A*x1) があります。

    前述の行列は、ディリクレ境界条件付きの 100 x 100 グリッドでのラプラシアンの離散化を表します。これは、修正不完全コレスキー前提条件子を使用することで、より効果的に実行できる可能性があることを意味します。

  3. 修正不完全コレスキー前提条件子を作成するには、micholオプションを以下のように使用します。

    L = ichol(A,struct(‘michol’,’on’);
    [x2,fl2,rr2,it2,rv2] = pcg(A,b,1e-8,100,L,L’);

    この場合、47 回の反復で収束できます。

  4. 以下のように、初期推定 (反復数 0) から始めて各残差履歴をプロットすることで、前提条件が pcg の収束レートにどのように影響するかを確認できます。

    figure;
    semilogy(0:it0,rv0/norm(b),'b.');
    hold on;
    semilogy(0:it1,rv1/norm(b),'r.');
    semilogy(0:it2,rv2/norm(b),'k.');
    legend('No Preconditioner','IC(0)','MIC(0)');
    xlabel('iteration number');
    ylabel('relative residual');
    hold off;
    

参考文献

[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.

参考

| | | | | | | | | |

この情報は役に立ちましたか?