ilu
不完全 LU 分解
说明
示例
ilu 函数提供三种类型的不完全 LU 分解：零填充分解 (ILU(0))、Crout 版本 (ILUC)，以及具有阈值调降和主元消去的分解 (ILUTP)。
默认情况下，ilu 执行稀疏矩阵输入的零填充不完全 LU 分解。例如，求具有 7840 个非零值的稀疏矩阵的完全和不完全分解。其完全 LU 因子有 126,478 个非零值，其具有零填充的不完全 LU 因子有 7840 个非零值，与 A 的数量相同。
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 进行置换和重新缩放，以创建一个新矩阵 B = R*P*A*C，它具有更好的条件数且对角线元只有 1 和 -1。
[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
您可以通过更改调降容差以获得不完全 LU 因子，这些因子更稀疏，但在逼近 B 时不太准确。例如，以下绘图显示了不完全因子的密度与输入矩阵的密度之比，以及不完全分解的相对误差，它们分别相对于调降容差的图。
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 的预条件子来求解线性系统。
加载 west0479，它是一个非对称的 479×479 实稀疏矩阵。
load west0479
A = west0479;定义 b 以使 的实际解是全为 1 的向量。
b = full(sum(A,2));
设置容差和最大迭代次数。
tol = 1e-12; maxit = 20;
使用 bicgstab 根据请求的容差和迭代次数求解。指定五个输出以返回有关求解过程的信息：
x是计算A*x = b所得的解。fl0是指示算法是否收敛的标志。rr0是计算的解x的相对残差。it0是计算x时所用的迭代序号。rv0是 的由每个二分之一迭代的残差历史记录组成的向量。
[x,fl0,rr0,it0,rv0] = bicgstab(A,b,tol,maxit); fl0
fl0 = 1
rr0
rr0 = 1
it0
it0 = 0
bicgstab 未在请求的 20 次迭代内收敛至请求的容差 1e-12，因此 fl0 为 1。实际上，bicgstab 的行为太差，因此初始估计值 x0 = zeros(size(A,2),1) 是最佳解，并会返回最佳解（如 it0 = 0 所示）。
为了有助于缓慢收敛，您可以指定预条件子矩阵。由于 A 是非对称的，请使用 ilu 生成预条件子 。指定调降容差，以忽略值小于 1e-6 的非对角线元。通过指定 L 和 U 作为 bicgstab 的输入，求解预条件方程组 。
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 预条件子产生的相对残差 rr1 小于请求的容差 1e-12。输出 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")
输入参数
输出参量
参考
[1] Saad, Y. "Preconditioning Techniques", chap. 10 in Iterative Methods for Sparse Linear Systems. The PWS Series in Computer Science. Boston: PWS Pub. Co, 1996.
