luinc (MATLAB Function 参考文献)

不完全な LU 行列分解

表示

luinc(X,'0')
[L,U] = luinc(X,'0')
[L,U,P] = luinc(X,'0') 
luinc(X,droptol)
luinc(X,options)
[L,U] = luinc(X,options)
[L,U] = luinc(X,droptol)
[L,U,P] = luinc(X,options)
[L,U,P] = luinc(X,droptol)

詳細

luinc は、単位下三角行列、上三角行列、置換行列を作ります。

luinc(X,'0') は、正方スパース行列のレベル0の不完全な LU 分解を計算します。三角因子は、オリジナルのスパース行列 X の置換と同じスパースパターンをもち、それらの積はスパースパターン上で置換された X と一致します。luinc(X,'0') は、因子の下三角部分と同じ行列にある上三角因子を出力します。置換の情報は失われていますが、nnz(luinc(X,'0')) = nnz(X) で、キャンセルによるゼロの例外があります。

[L,U] = luinc(X,'0') は、L に置換行列と下三角行列の積を、U に上三角行列を出力します。L, U,X のスパースパターンを比較することはできませんが、非ゼロの数はキャンセルによるL と U のゼロを除いて、保持されます。

  nnz(L)+nnz(U) = nnz(X)+n, ここで、X は、 n 行 n 列です。

積 L*U は、スパースパターン上で、X と一致します。(L*U).*spones(X)-X は、eps程度の要素をもちます。

[L,U,P] = luinc(X,'0') は、L に単位下三角行列、Uに上三角行列、P に置換行列を出力します。Lは、置換された X の下三角行列と同じスパースパターンをもちます。

  spones(L) = spones(tril(P*X))

P*X がゼロである Lの対角の 1 の可能な例外と、P*X が 0 でないキャンセルによる L のゼロの例外があります。Uは、P*X の上三角と同じスパースパターンをもちます。

  spones(U) = spones(triu(P*X))

P*X が 0でない部分でのキャンセルによる U の 0の例外があります。積 L*U は、スパースパターン上で置換行列P*X と丸め誤差内で一致します。(L*U).*spones(P*X)-P*Xは、eps 程度の要素をもちます。

luinc(X,droptol) は、ドロップトレランスを使ってスパース行列の不完全な LU 分解を計算します。droptol は、非負のスカラでなくてはいけません。luinc(X,droptol) は、lu(X) により出力される完全なLU因子の近似を行います。この近似は、ドロップトレランスが、0になるまで、ドロップトレランスを小さい値に変更する毎に改善されます。このとき、lu(X)で行われるような完全な LU分解が行われます。

三角行列の不完全因子の各列 j が計算されるに連れて、ローカルなドロップトレランス（ドロップトレランスとX の対応する列のノルムの積）

  droptol*norm(X(:,j))

よりも小さな要素は、適切な因子からドロップされます。

このドロップ則の唯一の例外は、上三角行列因子の対角要素で、これは特異因子を除くように保持されています。

luinc(X,options) は、droptol, milu, udiag, thresh の任意の組み合わせで使われる4つまでのフィールドをもつ構造体を指定します。optionsの他のフィールドは無視されます。

droptol は、不完全な因子分解のドロップトレランスです。

milu が 1 の場合、luinc は、列のドロップされた要素を上三角行列の対角要素から引いた、修正不完全LU 分解を作成します。デフォルト値は 0 です。

udiag が 1 の場合、上三角行列の対角上のゼロは、ローカルなドロップトレランスで置き換えられます。デフォルトは、0 です。

thresh は、0（対角ピボットを行う）と 1 の間のピボットスレッシュホールドで、デフォルトで、ピボットを行うために列中の最大要素を常に選択します。threshについては、lu で詳しく記述しています。

luinc(X,options) は、optionsが唯一のフィールドとして、droptol をもつ場合、luinc(X,droptol)と同じです。

[L,U] = luinc(X,options) は、Lに単位下三角単位行列の置換を、U に上三角行列を出力します。積L*U は、 X の近似です。luinc(X,options)は、同じ行列の中に因子の下三角部と上三角部を出力します。置換の情報は失われています。

[L,U] = luinc(X,options) は、optionsが、唯一のフィールドとして、droptol をもつ場合、luinc(X,droptol)と同じです。

[L,U,P] = luinc(X,options) は、L に単位下三角行列、Uに上三角行列、P に置換行列を出力します。Uの非ゼロ要素は、つぎの式を満足します。

  abs(U(i,j)) >= droptol*norm((X:,j)),

例外として、この規範を満たさないにも拘わらず、保持されている対角要素があります。L の要素は、ピボットによってスケーリングされる前にローカルなドロップトレランスに対してテストされるので、L の非ゼロ要素に対して、以下のようになります。

abs(L(i,j)) >= droptol*norm(X(:,j))/U(j,j).

積 L*U は、置換されたP*X への近似です。

[L,U,P] = luinc(X,options) は、options が唯一のフィールドとして、droptol をもつ場合、[L,U,P] = luinc(X,droptol) と同じです。

注意

不完全分解は、大きなスパース線形方程式システムを解くための前提条件として有効です。下三角行列は主対角要素すべてに1をもちますが、上三角行列の対角上の1つの0は行列を特異にします。ドロップトレランスをもつ不完全分解は、上三角行列の対角上に0がある場合、ワーニングメッセージを表示します。同様に、ゼロ対角を置き換えるためにオプションudiagを使うと、問題の徴候がなくなるだけで、それを解くことはしません。前提条件行列は特異ではありませんが、使用するのに有効ではなく、ワーニングメッセージが出力されます。

制限

luinc(X,'0') は、正方行列でのみ有効です。

例題

スパース行列について考え、LU分解を計算します。

load west0479;
S = west0479;
LU = lu(S);

レベル 0 の不完全 LU分解を計算します。

[L,U,P] = luinc(S,'0');
D = (L*U).*spones(P*S)-P*S;

spones(U) と spones(triu(P*S)) は、同じです。

code>spones(L) と spones(tril(P*S)) は、対角上の 73 箇所で一致しません。ここでは、L は 1で、P*S は 0 です。(206,113)でも一致せず、キャンセルのため L は 0で、P*S は -1 です。Dは、 eps 程度の要素をもちます。

[IL0,IU0,IP0] = luinc(S,0);
[IL1,IU1,IP1] = luinc(S,1e-10);
                .
                .
                .

0 のドロップトレランスは、完全な LU分解を行います。ドロップトレランスを大きくすると、因子のスパース性を増加させます(非ゼロを減少させます）が、下の2番目の図のドロップトレランスとnorm(L*U-P*S,1)/norm(S,1) のプロットのように、因子の誤差を増加させます。

アルゴリズム

luinc(X,'0') は、部分ピボットによる LU分解の変形である "KJI" を基本にしています。X の非ゼロである位置に対してのみアップデートが行われています。

luinc(X,droptol) と luinc(X,options)は、スパース行列に対する列指向の lu に基づいています。

参考

lu, cholinc, bicg

参考文献

Saad, Yousef, Iterative Methods for Sparse Linear Systems, PWS Publishing Company, 1996, Chapter 10 - Preconditioning Techniques.

lu magic