Main Content

本页翻译不是最新的。点击此处可查看最新英文版本。

ilu

不完全 LU 分解

说明

示例

[L,U] = ilu(A) 用零填充执行稀疏矩阵 A 的不完全 LU 分解,并返回下三角矩阵 L 和上三角矩阵 U

[L,U,P] = ilu(A) 还返回置换矩阵 P,并满足 LUP*AA*P 的不完全因子。默认情况下,P 是用于不使用主元消去的不完全 LU 分解的单位矩阵。

示例

W = ilu(A) 返回 LU 因子的非零值。输出 W 等于 L + U - speye(size(A))

示例

[___] = ilu(A,options) 使用结构体 options 指定的选项对 A 执行不完全 LU 分解。

例如,通过将 optionstype 字段设置为 "ilutp",您可以使用主元消去执行不完全 LU 分解。然后,通过将 milu 字段设置为 "row""col",可以指定保留修正不完全 LU 分解的行值总和或列值总和。这种选项组合返回置换矩阵 P,使得 LU"row" 选项的 A*P 的不完全因子(其中 U 是列置换的),或 LU"col"P*A 的不完全因子(其中 L 是行置换的)。

示例

全部折叠

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 contains an axes object. The axes object with title Ratio of nonzeros of LU factors with respect to B, xlabel drop tolerance, ylabel nnz(L)+nnz(U)-size(B,1)/nnz(B) contains an object of type line.

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")

Figure contains an axes object. The axes object with title Relative error of the incomplete LU factorization, xlabel drop tolerance, ylabel $||BP-LU|| indexOf F baseline \ ,/ \ ,||B|| indexOf F$ baseline contains an object of type line.

在此示例中,具有阈值调降的不完全 LU 分解的相对误差与调降容差处于相同的数量级(但不能保证一定会发生此情况)。

使用不完全 LU 分解作为 bicgstab 的预条件子来求解线性系统。

加载 west0479,它是一个非对称的 479×479 实稀疏矩阵。

load west0479
A = west0479;

定义 b 以使 Ax=b 的实际解是全为 1 的向量。

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

设置容差和最大迭代次数。

tol = 1e-12;
maxit = 20;

使用 bicgstab 根据请求的容差和迭代次数求解。指定五个输出以返回有关求解过程的信息:

  • x 是计算 A*x = b 所得的解。

  • fl0 是指示算法是否收敛的标志。

  • rr0 是计算的解 x 的相对残差。

  • it0 是计算 x 时所用的迭代序号。

  • rv0b-Ax 的由每个二分之一迭代的残差历史记录组成的向量。

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

bicgstab 未在请求的 20 次迭代内收敛至请求的容差 1e-12,因此 fl01。实际上,bicgstab 的行为太差,因此初始估计值 x0 = zeros(size(A,2),1) 是最佳解,并会返回最佳解(如 it0 = 0 所示)。

为了有助于缓慢收敛,您可以指定预条件子矩阵。由于 A 是非对称的,请使用 ilu 生成预条件子 M=L U。指定调降容差,以忽略值小于 1e-6 的非对角线元。通过指定 LU 作为 bicgstab 的输入,求解预条件方程组 A M-1 M x=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 预条件子产生的相对残差 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")

Figure contains an axes object. The axes object with xlabel Iteration number, ylabel Relative residual contains 3 objects of type line, constantline. These objects represent No preconditioner, ILU preconditioner, Tolerance.

输入参数

全部折叠

输入矩阵,指定为稀疏方阵。

LU 分解选项,指定为包含最多五个字段的结构体:

字段名称摘要描述

type

不完全 LU 分解的类型

type 的有效值为:

  • "nofill"(默认值)- 执行具有 0 填充级别的 ILU 分解(称为 ILU(0))。如果将 type 设置为 "nofill",则仅使用 milu 选项;所有其他字段都将被忽略。

  • "crout" - 执行 ILU 分解的 Crout 版本,称为 ILUC。如果将 type 设置为 "crout",则仅使用 droptolmilu 选项;所有其他字段都将被忽略。

  • "ilutp" - 执行带阈值和主元消去的 ILU 分解,称为 ILUTP。仅可以对此类型执行主元消去。

默认值为 "nofill"

droptol

不完全 LU 分解的调降容差

droptol 的有效值是任何非负标量。U 的非零项满足 abs(U(i,j)) >= droptol*norm(A(:,j)),但对角线元素除外,无论它们是否满足此标准,对角线元素都将保留。在使用主元调整 L 的元之前,将根据局部调降容差检验这些元,以使 L 的非零项满足 abs(L(i,j)) >= droptol*norm(A(:,j))/U(j,j)

默认值为 0,这会生成完全的 LU 分解。

milu

修正不完全 LU 分解的类型

milu 的有效值为:

  • "off" - 不生成修正不完全 LU 分解。

  • "row" - 生成行值总和修正不完全 LU 分解。补偿上三角因子 U 的对角线元素,以便保留原始矩阵的行值总和。即 A*e = L*U*e,其中 e 是由 1 组成的列向量。

  • "col" - 生成列总和修正不完全 LU 分解。补偿上三角因子 U 的对角线元素,以便保留原始矩阵的列总和。即 e'*A = e'*L*U

默认值为 "off"

udiag

是否替换 U 的零对角线元素

udiag 的有效值为 01。如果 udiag1,则上三角因子 U 的任何零对角线元素都将替换为局部调降容差,以尝试避免奇异因子。

默认值为 0

thresh

主元阈值

有效值介于 0(算法选择对角主元)和 1(算法选择列中模最大的条目作为主元)之间。当列中对角线元素的模小于该列中任何次对角线元素模的 thresh 倍时,就会发生主元消去。

默认值为 1

输出参数

全部折叠

下三角因子,以稀疏方阵形式返回。

  • options.type 指定为 "nofill"(默认值)或 "crout" 时,则 L 返回为一个单位下三角矩阵(即,主对角线上元素为 1 的下三角矩阵)。

  • options.type 指定为 "ilutp" 并且 options.milu 指定为 "col" 时,L 返回为一个行置换的单位下三角矩阵。

上三角因子,以稀疏方阵形式返回。

  • options.type 指定为 "nofill"(默认值)或 "crout" 时,U 返回为一个上三角矩阵。

  • options.type 指定为 "ilutp" 并且 options.milu 指定为 "row" 时,U 返回为一个列置换的上三角矩阵。

置换矩阵,以稀疏方阵形式返回。

options.type 指定为 "nofill"(默认值)或 "crout" 时,P 返回为一个与 A 大小相同的单位矩阵,因为这两种类型都不执行主元消去。

LU 因子的非零值,以稀疏方阵形式返回,其中 W = L + U - speye(size(A))

提示

参考

[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.

扩展功能

版本历史记录

在 R2007a 中推出