ilu

関数 ilu では、ゼロで埋める分解 (ILU(0))、Crout バージョン (ILUC)、棄却とピボットのしきい値をもつ分解 (ILUTP) の 3 種類の不完全 LU 分解を行います。

既定では、ilu はスパース行列入力のゼロで埋める不完全 LU 分解を実行します。たとえば、7840 個の非ゼロ要素をもつスパース行列の完全分解と不完全分解を求めるとします。その完全 LU 因子は 126,478 個の非ゼロ要素をもち、ゼロ埋め込みを使用するその不完全 LU 因子は、A と同じ数の 7840 個の非ゼロ要素をもちます。

A = gallery("neumann",1600) + speye(1600);
n = nnz(A)

n = 7840

n = nnz(lu(A))

n = 126478

n = nnz(ilu(A))

n = 7840

ゼロで埋める分解ではその LU 因子で入力行列のスパース パターンが保持されるため、分解の相対誤差は、A の非ゼロ要素のパターンでは基本的にゼロです。

[L,U] = ilu(A);
e = norm(A-(L*U).*spones(A),"fro")/norm(A,"fro")

e = 4.8874e-17

ただし、このようなゼロで埋められた因子の積は、元の行列を十分に近似していません。

e = norm(A-L*U,"fro")/norm(A,"fro")

e = 0.0601

精度を高めるには、非スパースを使う他のタイプの不完全 LU 分解を使用します。たとえば、棄却許容誤差が 1e-6 の Crout バージョンを使用します。

options.droptol = 1e-6;
options.type = "crout";
[L,U] = ilu(A,options);

Crout バージョンの不完全分解の LU 因子は 51,482 個の非ゼロ要素をもちます (完全分解より少ない)。非スパースを使用した不完全 LU 因子の積のほうが、元の行列を良く近似しています。

n = nnz(ilu(A,options))

n = 51482

e = norm(A-L*U,"fro")./norm(A,"fro")

e = 9.3040e-07

比較のために、同じ入力行列 A について棄却とピボットのしきい値をもつ不完全分解を行うと、Crout バージョンと同様の結果になります。

options.type = "ilutp";
[L,U,P] = ilu(A,options);
n = nnz(ilu(A,options))

n = 51541

norm(P*A-L*U,"fro")./norm(A,"fro")

ans = 9.4960e-07

不完全 LU 分解の棄却許容誤差を変更してスパース行列を因数分解します。

行列 west0479 を読み込みます。これは 479 行 479 列の実数値の非対称スパース行列です。condest を使用して行列の条件数を推定します。

load west0479
A = west0479;
c1 = condest(A)

c1 = 1.4244e+12

equilibrateを使用して行列の条件数を改善します。元の行列 A を並べ替えて再スケーリングし、より良い条件数をもち、かつ対角要素が 1 と –1 のみの新しい行列 B = R*P*A*C を作成します。

[P,R,C] = equilibrate(A);
B = R*P*A*C;
c2 = condest(B)

c2 = 5.1036e+04

棄却とピボットのしきい値をもち、行の和を保持する B の不完全 LU 分解を実行するようにオプションを指定します。比較のために、まず、非スパースの棄却許容誤差としてゼロを指定して、完全な LU 分解を行います。

options.type = "ilutp";
options.milu = "row";
options.droptol = 0;
[L,U,P] = ilu(B,options);

この因数分解では、入力行列 B の近似の精度が非常に高くなりますが、因子は B より大幅に密になります。

e = norm(B*P-L*U,"fro")/norm(B,"fro")

e = 1.0345e-16

nLU = nnz(L)+nnz(U)-size(B,1)

nLU = 19118

nB = nnz(B)

nB = 1887

棄却許容誤差を変更して、スパース性を上げて B の近似の精度を下げた不完全 LU 因子を取得できます。たとえば、以下のプロットは、棄却許容誤差に対してプロットした入力行列の密度に対する不完全因子の密度の割合と、不完全分解の相対誤差を示します。

ntols = 20;
tau = logspace(-6,-2,ntols);
e = zeros(1,ntols);
nLU = zeros(1,ntols);
for k = 1:ntols
    options.droptol = tau(k);
    [L,U,P] = ilu(B,options);
    nLU(k) = nnz(L)+nnz(U)-size(B,1);
    e(k) = norm(B*P-L*U,"fro")/norm(B,"fro");
end
figure
semilogx(tau,nLU./nB,LineWidth=2)
title("Ratio of nonzeros of LU factors with respect to B")
xlabel("drop tolerance")
ylabel("nnz(L)+nnz(U)-size(B,1)/nnz(B)",FontName="FixedWidth")

figure
loglog(tau,e,LineWidth=2)
title("Relative error of the incomplete LU factorization")
xlabel("drop tolerance")
ylabel("$||BP-LU||_F\,/\,||B||_F$",Interpreter="latex")

この例では、しきい値棄却による不完全 LU 分解の相対誤差は、棄却許容誤差と同じ次元にあります (常に同じとは限りません)。

不完全 LU 分解をbicgstabの前処理行列として使用して線形システムを解きます。

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

load west0479
A = west0479;

$\mathit{Ax}=\mathit{b}$ に対する真の解がすべて 1 のベクトルになるように、b を定義します。

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

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

tol = 1e-12;
maxit = 20;

bicgstab を使用して、要求された許容誤差と反復回数で解を求めます。解法プロセスに関する情報を返す出力を 5 つ指定します。

[x,fl0,rr0,it0,rv0] = bicgstab(A,b,tol,maxit); 
fl0

fl0 = 1

rr0

rr0 = 1

it0

it0 = 0

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

遅い収束への対応として、前処理行列を指定できます。A は非対称であるため、ilu を使用して前処理行列 $\mathit{M}=\mathit{L}\mathit{U}$ を生成します。棄却許容誤差を指定して、1e-6 よりも小さい値をもつ非対角エントリを無視します。bicgstab への入力として L および U を指定して、前処理された系 ${\mathit{A}\mathit{M}}^{-1}\mathit{M}\mathit{x}=\mathit{b}$ を解きます。

options = struct("type","ilutp","droptol",1e-6);
[L,U] = ilu(A,options);
[x_precond,fl1,rr1,it1,rv1] = bicgstab(A,b,tol,maxit,L,U);
fl1

fl1 = 0

rr1

rr1 = 5.9892e-14

it1

it1 = 3

ilu 前処理行列を使用すると、3 回目の反復で 1e-12 の要求された許容誤差より少ない相対残差 rr1 が生成されます。出力 rv1(1) は norm(b)、出力 rv1(end) は norm(b-A*x1) になります。

半分の反復ごとでの相対残差をプロットして、bicgstab の進行状況を確認できます。指定された許容誤差のラインと共に、それぞれの解の残差履歴をプロットします。

semilogy((0:numel(rv0)-1)/2,rv0/norm(b),"-o")
hold on
semilogy((0:numel(rv1)-1)/2,rv1/norm(b),"-o")
yline(tol,"r--");
legend("No preconditioner","ILU preconditioner","Tolerance",Location="East")
xlabel("Iteration number")
ylabel("Relative residual")