ichol
不完全コレスキー分解
説明
例
不完全コレスキー分解
この例では、修正不完全コレスキー分解を生成します。
対称正定値行列 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')
相対誤差は通常、棄却許容誤差と同じ次元にありますが、必ず同じとは保証できません。
前処理行列としての ichol
の使用
この例では、収束を改善するための前処理行列として、不完全コレスキー分解を使用する方法を説明します。
対称正定値行列 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
が既定であるため、この構文ではタイプを指定する必要はありませんが、指定することをお勧めします。
options.type = 'nofill'; options.michol = 'on'; L2 = ichol(A,options); 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)
となります。
行列が楕円型偏微分方程式に基づくため、ゼロ埋め込み不完全コレスキーと同じように、しきい値棄却分解も変更により効果を上げることができます (すなわち options.michol = 'on'
)。MIC(0)
と同様に、変更済みしきい値ベース棄却不完全コレスキーは、前処理行列の処理を定数ベクトルに保持します。すなわち、norm(A*e-L*(L'*e))
はほぼゼロとなります。
diagcomp
オプションの使用
この例では、ichol
の diagcomp
オプションの使用方法を説明します。
正定値行列の不完全コレスキー分解は、常に存在するわけではありません。以下のコードでは、乱数対称正定値行列を構築し、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
を近似すること、対角補正の ichol
が L
を構築し、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');
入力引数
A
— 入力行列
スパース正方行列
入力行列。スパース正方行列として指定します。
options
— コレスキー分解オプション
構造体
コレスキー分解オプション。最大 5 つのフィールドをもつ構造体として指定します。
フィールド名 | まとめ | 説明 |
---|---|---|
type | 分解のタイプ。 | 実行する不完全コレスキーの種類を指定します。このフィールドの有効な値は、 |
droptol | タイプが 'ict' である場合の棄却許容誤差 | ICT を実行するときに棄却許容誤差として使用される非負スカラー。ローカル棄却許容誤差よりも絶対値が小さい要素は、まったく棄却されない対角要素を除き、結果の要素から棄却されます。分解のステップ |
michol | 変更された不完全コレスキーを実行するかどうかの指定 | 変更された不完全コレスキー (MIC) を実行するかどうかの指定フィールドは |
diagcomp | 指定した係数で補正された不完全コレスキーの実行 | 不完全コレスキー係数を構成するときに、グローバル対角シフト |
shape | 参照され、返される三角形の決定 | 有効な値は |
例: options.type = 'nofill'; options.michol = 'on'; L = ichol(A,options);
ヒント
参照
[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.
拡張機能
スレッドベースの環境
MATLAB® の backgroundPool
を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool
を使用してコードを高速化します。
この関数はスレッドベースの環境を完全にサポートしています。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。
分散配列
Parallel Computing Toolbox™ を使用して、クラスターの結合メモリ上で大きなアレイを分割します。
使用上の注意事項および制限事項:
options
のtype
フィールドを指定する場合、'nofill'
でなければなりません。A
は、実対称または複素エルミートでなければなりません。
詳細については、分散配列を使用した MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。
バージョン履歴
R2011a で導入
MATLAB コマンド
次の MATLAB コマンドに対応するリンクがクリックされました。
コマンドを MATLAB コマンド ウィンドウに入力して実行してください。Web ブラウザーは MATLAB コマンドをサポートしていません。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)