ドキュメンテーション

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

ichol

不完全コレスキー分解

構文

L = ichol(A)
L = ichol(A,opts)

説明

L = ichol(A) は、ゼロ埋め込みを使用して A の不完全コレスキー分解を実行します。

L = ichol(A,opts) は、opts で指定されたオプションを使用して、A の不完全コレスキー分解を実行します。

既定の設定では、ichol は、A の下三角部分を参照し、下三角因子を作成します。

入力引数

A

スパース行列

opts

最大で 5 つのフィールドを含む構造体:

フィールド名まとめ説明
type分解のタイプ。

実行する不完全コレスキーの種類を指定します。このフィールドの有効な値は、'nofill' および 'ict' です。'nofill' バリアントは、ゼロ埋め込みを使用して不完全コレスキー を実行します (IC(0))。'ict' バリアントは、しきい値棄却を使用して不完全コレスキー を実行します (ICT)。既定値は 'nofill' です。

droptol タイプが 'ict' である場合の棄却許容誤差

ICT を実行するときに棄却許容誤差として使用される非負スカラー。ローカル棄却許容誤差よりも絶対値が小さい要素は、まったく棄却されない対角要素を除き、結果の要素から棄却されます。分解のステップ j における局所的な棄却許容誤差は、norm(A(j:end,j),1)*droptol です。'droptol' は、'type''nofill' である場合は無視されます。既定値は 0 です。

michol変更された不完全コレスキーを実行するかどうかの指定

変更された不完全コレスキー (MIC) を実行するかどうかの指定フィールドは 'on' または 'off' です。MIC を実行すると、関係 A*e = L*L'*e が成り立つよう、棄却された要素に応じて対角が補正されます。ここで、e = ones(size(A,2),1) です。既定値は 'off' です。

diagcomp指定した係数で補正された不完全コレスキーの実行

不完全コレスキー係数を構成するときに、グローバル対角シフト alpha として使用される非負スカラー。すなわち、A の不完全コレスキーを実行するのではなく、A + alpha*diag(diag(A)) の分解が行なわれます。既定値は 0 です。

shape参照され、返される三角形の決定

有効な値は 'upper' および 'lower' です。'upper' を指定した場合、 A の上三角だけが参照され、AR'*R で近似されるように、R が構築されます。'lower' を指定した場合、 A の下三角だけが参照され、AL*L' で近似されるように、L が構築されます。既定値は 'lower' です。

すべて折りたたむ

この例では、修正不完全コレスキー分解を生成します。

対称正定行列 A から始めます。

N = 100;
A = delsq(numgrid('S',N));

A は、Dirichlet 境界条件による 100 行 100 列の正方形グリッド上の負の 2 次元 5 点分散ラプラシアンです。A のサイズは 98*98 = 9604 となります (グリッドの境界がディリクレ条件を課すために使用されるので、10,000 ではない)。

埋め込みなしの不完全コレスキー分解は、A が非ゼロを含むため、同じ位置に非ゼロのみを含む分解です。このような分解は、算出が非常に容易です。積 L*L' は一般に、A とは大きく異なりますが、積 L*L' は該当するパターンで A と一致します。

L = ichol(A);
norm(A-L*L','fro')./norm(A,'fro')
ans = 0.0916
norm(A-(L*L').*spones(A),'fro')./norm(A,'fro')
ans = 4.9606e-17

また、ichol を使用して、しきい値棄却による不完全コレスキー分解を実行することもできます。棄却許容誤差が小さくなると、係数は密になり、積 L*L' はより優れた A の近似となる傾向があります。以下のプロットは、完全コレスキー係数の密度に対する不完全係数の密度の割合と、棄却許容誤差に対してプロットした不完全分解の相対誤差を示します。

n = size(A,1);
ntols = 20;
droptol = logspace(-8,0,ntols);
nrm = zeros(1,ntols);
nz = zeros(1,ntols);
nzComplete = nnz(chol(A,'lower'));
for k = 1:ntols
    L = ichol(A,struct('type','ict','droptol',droptol(k)));
    nz(k) = nnz(L);
    nrm(k) = norm(A-L*L','fro')./norm(A,'fro');
end
figure
loglog(droptol,nrm,'LineWidth',2)
title('Drop tolerance vs norm(A-L*L'',''fro'')./norm(A,''fro'')')

figure
semilogx(droptol,nz./nzComplete,'LineWidth',2)
title('Drop tolerance vs fill ratio ichol/chol')

相対誤差は通常、棄却許容誤差と同じ次元にありますが、必ず同じとは保証できません。

この例では、収束を改善するための前提条件子として、不完全コレスキー分解を使用する方法を説明します。

対称正定行列 A を作成します。

N = 100;
A = delsq(numgrid('S',N));

不完全コレスキー分解を pcg の前提条件子として作成します。定数ベクトルを右側に使用します。ベースラインとして、前提条件子なしで pcg を実行します。

b = ones(size(A,1),1);
tol = 1e-6;
maxit = 100;
[x0,fl0,rr0,it0,rv0] = pcg(A,b,tol,maxit);

fl0 = 1 は、pcg が最大許容反復回数内で相対残差を要求した許容誤差に収めないことを示します。埋め込みなしの不完全コレスキー分解を前提条件子として試行します。

L1 = ichol(A);
[x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L1,L1');

fl1 = 0 は、pcg が反復 59 回で要求した許容誤差に収束することを示します (it1 の値)。この行列はラプラシアンで離散化されていますが、変更済みの不完全コレスキーを使用すると、より良好な前提条件を作成できます。変更済みの不完全コレスキー分解は、演算子の処理を定数ベクトルに保持する近似分解を構築します。すなわち、norm(A*e-L*(L'*e)) は、norm(A-L*L','fro')/norm(A,'fro') がゼロ近傍でない場合でも、e = ones(size(A,2),1) でほぼゼロになります。nofill が既定であるため、この構文ではタイプを指定する必要はありませんが、指定することをお勧めします。

opts.type = 'nofill';
opts.michol = 'on';
L2 = ichol(A,opts);
e = ones(size(A,2),1);
norm(A*e-L2*(L2'*e))
ans = 3.7983e-14
[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L2,L2');

pcg はわずか 38 回の反復で収束します (fl2 = 0)。3 つの収束履歴をすべてプロットすると、以下の収束が得られます。

semilogy(0:maxit,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)');

このプロットは、変更済み不完全コレスキー前提条件がより速く収束を達成することを示します。

また、しきい値棄却による不完全コレスキー分解を試行することもできます。以下のプロットは、各種の棄却許容誤差で構築した前提条件子による pcg の収束を示します。

L3 = ichol(A, struct('type','ict','droptol',1e-1));
[x3,fl3,rr3,it3,rv3] = pcg(A,b,tol,maxit,L3,L3');
L4 = ichol(A, struct('type','ict','droptol',1e-2));
[x4,fl4,rr4,it4,rv4] = pcg(A,b,tol,maxit,L4,L4');
L5 = ichol(A, struct('type','ict','droptol',1e-3));
[x5,fl5,rr5,it5,rv5] = pcg(A,b,tol,maxit,L5,L5'); 

figure
semilogy(0:maxit,rv0./norm(b),'b-','linewidth',2);
hold on
semilogy(0:it3,rv3./norm(b),'b-.','linewidth',2);
semilogy(0:it4,rv4./norm(b),'b--','linewidth',2);
semilogy(0:it5,rv5./norm(b),'b:','linewidth',2);
legend('No Preconditioner','ICT(1e-1)','ICT(1e-2)', ...
   'ICT(1e-3)','Location','SouthEast');

棄却許容誤差 1e-2 で構築した不完全コレスキー前提条件子は、ICT(1e-2) となります。

行列が楕円型偏微分方程式に基づくため、ゼロ埋め込み不完全コレスキーと同じように、しきい値棄却分解も変更により効果を上げることができます (すなわち opts.michol = 'on')。MIC(0) と同様に、変更済みしきい値ベース棄却不完全コレスキーは、前提条件子の処理を定数ベクトルに保持します。すなわち、norm(A*e-L*(L'*e)) はほぼゼロとなります。

この例では、icholdiagcomp オプションの使用方法を説明します。

正定行列の不完全コレスキー分解は、常に存在するわけではありません。以下のコードでは、乱数対称正定行列を構築し、pcg を使用して線形システムを解決しようとします。

S = rng('default');
A = sprandsym(1000,1e-2,1e-4,1);
rng(S);
b = full(sum(A,2));
[x0,fl0,rr0,it0,rv0] = pcg(A,b,1e-6,100);

収束が達成されていないため、不完全コレスキー前提条件子の作成を試行します。

L = ichol(A,struct('type','ict','droptol',1e-3));
Error using ichol
Encountered nonpositive pivot.

ichol が前述のように分解される場合、diagcomp オプションを使用して、シフトされた不完全コレスキー分解を構築できます。すなわち、L を構築するのではなく、L*L'A を近似すること、対角補正の icholL を構築し、L*L' が明示的に M を形成することなく M = A + alpha*diag(diag(A)) を近似することを意味します。対角優位の行列では不完全分解が常に存在するため、alpha を見つけて M を対角優位にすることができます。

alpha = max(sum(abs(A),2)./diag(A))-2
alpha = 62.9341
L1 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha));
[x1,fl1,rr1,it1,rv1] = pcg(A,b,1e-6,100,L1,L1');

ここで、pcg はまだ、指定の反復回数内に必要な許容誤差に収束することができませんが、pcg では、前提条件子なしよりも、この前提条件子のほうが良好な収束が得られています。alpha を小さくすると有効な場合があります。実験をある程度行うと、alpha の適切な値を得ることができます。

alpha = .1;
L2 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha));
[x2,fl2,rr2,it2,rv2] = pcg(A,b,1e-6,100,L2,L2');

これで、pcg が収束し、プロットには各前提条件子の収束履歴が示されます。

semilogy(0:100,rv0./norm(b),'b.');
hold on;
semilogy(0:100,rv1./norm(b),'r.');
semilogy(0:it2,rv2./norm(b),'k.');
legend('No Preconditioner','\alpha \approx 63','\alpha = .1');
xlabel('Iteration Number');
ylabel('Relative Residual');

ヒント

  • このルーチンで得られる係数は、pcgminres など、反復メソッドで解かれる線形等式系の前提条件子として使用できることがあります。

  • ichol が機能するのはスパース正方行列だけです。

参考文献

[1] Saad, Yousef. “Preconditioning Techniques.” Iterative Methods for Sparse Linear Systems. PWS Publishing Company, 1996.

[2] Manteuffel, T.A. “An incomplete factorization technique for positive definite linear systems.” Math. Comput. 34, 473–497, 1980.

参考

| | |