ilu
不完全 LU 分解
说明
示例
不同类型的不完全 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 分解的调降容差
更改不完全 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 分解的相对误差与调降容差处于相同的数量级(但不能保证一定会发生此情况)。
使用 ilu
作为预条件子来求解线性系统
使用不完全 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")
输入参数
A
— 输入矩阵
稀疏方阵
输入矩阵,指定为稀疏方阵。
options
— LU 分解选项
结构体
LU 分解选项,指定为包含最多五个字段的结构体:
字段名称 | 摘要 | 描述 |
---|---|---|
| 不完全 LU 分解的类型 |
默认值为 |
| 不完全 LU 分解的调降容差 |
默认值为 |
| 修正不完全 LU 分解的类型 |
默认值为 |
| 是否替换 U 的零对角线元素 |
默认值为 |
| 主元阈值 | 有效值介于 默认值为 |
输出参量
L
— 下三角因子
稀疏方阵
下三角因子,以稀疏方阵形式返回。
当
options.type
指定为"nofill"
(默认值)或"crout"
时,则L
返回为一个单位下三角矩阵(即,主对角线上元素为 1 的下三角矩阵)。当
options.type
指定为"ilutp"
并且options.milu
指定为"col"
时,L
返回为一个行置换的单位下三角矩阵。
U
— 上三角因子
稀疏方阵
上三角因子,以稀疏方阵形式返回。
当
options.type
指定为"nofill"
(默认值)或"crout"
时,U
返回为一个上三角矩阵。当
options.type
指定为"ilutp"
并且options.milu
指定为"row"
时,U
返回为一个列置换的上三角矩阵。
P
— 置换矩阵
稀疏方阵
置换矩阵,以稀疏方阵形式返回。
当 options.type
指定为 "nofill"
(默认值)或 "crout"
时,P
返回为一个与 A
大小相同的单位矩阵,因为这两种类型都不执行主元消去。
W
— LU 因子的非零值
稀疏方阵
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.
扩展功能
基于线程的环境
使用 MATLAB® backgroundPool
在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool
加快代码运行速度。
此函数完全支持基于线程的环境。有关详细信息,请参阅在基于线程的环境中运行 MATLAB 函数。
分布式数组
使用 Parallel Computing Toolbox™ 在集群的组合内存中对大型数组进行分区。
用法说明和限制:
如果在结构体数组
options
中包含字段type
,则必须将其设置为'nofill'
。
有关详细信息,请参阅Run MATLAB Functions with Distributed Arrays (Parallel Computing Toolbox)。
版本历史记录
在 R2007a 中推出
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)