ドキュメンテーション

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

cgs

共役傾斜二乗法

構文

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

説明

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

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

関数 cgs が収束する場合、それを伝えるメッセージが表示されます。関数 cgs が設定された最大反復回数に達しても収束しないか、または何らかの理由で停止した場合、警告メッセージが、相対残差 norm(b-A*x)/norm(b) と、計算が停止した反復回数と共に表示されます。

cgs(A,b,tol) は、この手法の許容誤差 tol を指定します。tol[] の場合、関数 cgs は既定の設定の 1e-6 を使用します。

cgs(A,b,tol,maxit) は、最大反復回数 maxit を指定します。maxit[] の場合、関数 cgs は既定の設定の min(n,20) を使用します。

cgs(A,b,tol,maxit,M) および cgs(A,b,tol,maxit,M1,M2) は、前提条件子 M または M = M1*M2 を使用し、方程式系 inv(M)*A*x = inv(M)*bx について効率的に解きます。M[] の場合、関数 cgs は前提条件を適用しません。M は、mfun(x)M\x を返すような関数ハンドル mfun です。

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

[x,flag] = cgs(A,b,...) は、解 x と、cgs の収束に関するフラグを返します。

フラグ

収束

0

関数 cgs は希望する許容誤差 tol に反復回数 maxit 以内で収束しました。

1

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

2

前提条件 M は悪条件です。

3

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

4

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

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

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

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

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

行列入力で cgs を使用する

A = gallery('wilk',21);
b = sum(A,2);
tol = 1e-12;  maxit = 15; 
M1 = diag([10:-1:1 1 1:10]);
x = cgs(A,b,tol,maxit,M1);

以下のメッセージが表示されます。

cgs converged at iteration 13 to a solution with 
relative residual 2.4e-016.

関数ハンドルで cgs を使用する

この例では、前述の例の行列 A に代えて、行列-ベクトル積関数 afun のハンドルを使用し、前提条件 M1 に代えてバックソルブ関数 mfun のハンドルを使用します。この例は、以下を行うファイル run_cgs に含まれています。

  • 第 1 引数として関数ハンドル @afun を指定し、関数 cgs を呼び出す。

  • run_cgs 内のすべての変数を afun および myfun が利用できるようにするため、afun を入れ子関数として含む。

以下に run_cgs のコードを示します。

function x1 = run_cgs
n = 21;
b = afun(ones(n,1));
tol = 1e-12;  maxit = 15; 
x1 = cgs(@afun,b,tol,maxit,@mfun);
 
    function y = afun(x)
        y = [0; x(1:n-1)] + ...
              [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x + ...
              [x(2:n); 0];
    end
 
    function y = mfun(r)
        y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
    end
end

次のように入力した場合、

x1 = run_cgs

MATLAB® ソフトウェアは次の値を返します。

cgs converged at iteration 13 to a solution with 
relative residual 2.4e-016.

前提条件子と cgs の使用

この例は、前提条件の使用方法を示します。

479 行 479 列の非対称実スパース行列 west0479 を読み込みます。

load west0479;
A = west0479;

真の解がすべて 1 のベクトルになるように、b を定義します。

b = full(sum(A,2));

許容誤差と最大反復回数を設定します。

tol = 1e-12;
maxit = 20;

cgs を使用して、要求された許容誤差と反復回数で解を求めます。

[x0,fl0,rr0,it0,rv0] = cgs(A,b,tol,maxit);

fl0 は、cgs が要求した反復回数 20 回以内に要求した許容誤差 1e-12 に収束しなかったため、1 となります。実際、cgs の動作が良好でないため、初期推定 (x0 = zeros(size(A,2),1)) が最適解であり、これは it0 = 0 で示される出力です。MATLAB は残差履歴を rv0 に格納します。

cgs の動作をプロットします。

semilogy(0:maxit,rv0/norm(b),'-o');
xlabel('Iteration number');
ylabel('Relative residual');

プロットでは、解が収束しないことが示されます。前提条件子を使用して結果を改善できます。

A は非対称であるため、ilu で前提条件子を作成します。

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-5));
Error using ilu
There is a pivot equal to zero.  Consider decreasing the 
drop tolerance or consider using the 'udiag' option.

MATLAB は不完全な LU を構成できません。前提条件子では無意味な特異因子が生成されるためです。

エラー メッセージで示されるように、減少させた棄却許容誤差でもう一度試すことができます。

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = cgs(A,b,tol,maxit,L,U);

cgs が相対残差 4.3851e-014 (rr1 の値) を導き出すため、fl1 は 0 です。棄却許容誤差 1e-6 の不完全 LU 分解を前提条件とした場合、3 回目の反復 (it1 の値) で相対残差が指定の許容誤差 1e-12 より小さくなります。出力 rv1(1)norm(b) で、出力 rv1(14)norm(b-A*x2) です。

初期推定 (反復回数が 0) から始まる各反復での相対残差をプロットして、関数 cgs の進行状況を確認できます。

semilogy(0:it1,rv1/norm(b),'-o');
xlabel('Iteration number');
ylabel('Relative residual');

参考文献

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

[2] Sonneveld, Peter, “CGS: A fast Lanczos-type solver for nonsymmetric linear systems,” SIAM J. Sci. Stat. Comput., January 1989, Vol. 10, No. 1, pp. 36–52.

拡張機能

R2006a より前に導入