目录
ichol函数的功能是对矩阵进行不完全 Cholesky 分解。
语法
L = ichol(A)
L = ichol(A,options)
说明
L = ichol(A) 通过零填充执行 A 的不完全 Cholesky 分解。A 必须为稀疏方阵。
L = ichol(A,options) 使用结构图 options 指定的选项对 A 执行不完全 Cholesky 分解。
默认情况下,ichol 引用 A 的下三角并生成下三角因子。
示例
不完全 Cholesky 分解
此示例生成不完全 Cholesky 分解。从一个对称正定矩阵A开始:
N = 100;
A = delsq(numgrid('S',N));
A 是带有 Dirichlet 边界条件的 100×100 方形网格上的二维、五点离散负 Laplacian 矩阵。A 的大小为 98*98 = 9604(并非 10000,因为网格边框用于施加 Dirichlet 条件)。
无填充的不完全 Cholesky 分解是指在 A 包含非零值的相同位置仅包含非零值的分解。此分解的计算非常容易。尽管 L*L' 乘积通常与 A 完全不同,但 L*L' 乘积将与 A 在其向上舍入模式上匹配。
L = ichol(A);
norm(A-L*L','fro')./norm(A,'fro')
ans = 0.0916
norm(A-(L*L').*spones(A),'fro')./norm(A,'fro')
ans = 4.9606e-17
ichol 也可用于通过阈值调降生成不完全 Cholesky 分解。当调降容差减小时,因子往往会变得更密集,而 L*L' 乘积往往更好地接近 A。以下绘图显示了不完全分解的相对误差对调降容差的图,以及不完全因子密度与完全 Cholesky 因子密度之比。
n = size(A,1);
ntols = 20;
droptol = logspace(-8,0,ntols);
nrm = zeros(1,ntols);
nz = zeros(1,ntols);
nzComplete = nnz(chol(A,'lower'));
for k = 1:ntols
L = ichol(A,struct('type','ict','droptol',droptol(k)));
nz(k) = nnz(L);
nrm(k) = norm(A-L*L','fro')./norm(A,'fro');
end
figure
loglog(droptol,nrm,'LineWidth',2)
title('Drop tolerance vs norm(A-L*L'',''fro'')./norm(A,''fro'')')
如图所示:
figure
semilogx(droptol,nz./nzComplete,'LineWidth',2)
title('Drop tolerance vs fill ratio ichol/chol')
如图所示:
该相对误差通常与调降容差的量级相当,但不能保证一定如此。
使用 ichol
作为预条件子
此示例说明如何使用不完全 Cholesky 分解作为预条件子来提高收敛。创建一个对称正定矩阵 A
。
N = 100;
A = delsq(numgrid('S',N));
创建一个不完全 Cholesky 分解,作为 pcg 的预条件子。在右侧使用一个常向量。先在没有预条件子的情况下执行 pcg,以作为基准。
b = ones(size(A,1),1);
tol = 1e-6;
maxit = 100;
[x0,fl0,rr0,it0,rv0] = pcg(A,b,tol,maxit);
请注意,fl0 = 1 指示 pcg 未在允许的最大迭代次数中使相对残差趋向于请求的公差。请尝试使用无填充的不完全 Cholesky 分解作为预条件子。
L1 = ichol(A);
[x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L1,L1');
fl1 = 0,指示 pcg 收敛至请求的公差,经过 59 次迭代(it1 值)后实现收敛。但是,由于此矩阵是一个离散的 Laplacian 矩阵,因此使用修正不完全 Cholesky 分解可创建一个更好的预条件子。修正不完全 Cholesky 分解会构造一个近似的分解,该分解会保留算子对常向量的作用。也就是说,对于 e = ones(size(A,2),1),norm(A*e-L*(L'*e)) 将约为零,即使 norm(A-L*L','fro')/norm(A,'fro') 不接近零。不必为此语法指定类型,因为 nofill 是默认值,但这是一种好的做法。
options.type = 'nofill';
options.michol = 'on';
L2 = ichol(A,options);
e = ones(size(A,2),1);
norm(A*e-L2*(L2'*e))
ans = 3.7983e-14
[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L2,L2');
pcg 收敛 (fl2 = 0),但仅用了 38 次迭代。绘制所有三种收敛历史记录将显示以下收敛情况。
semilogy(0:maxit,rv0./norm(b),'b.');
hold on
semilogy(0:it1,rv1./norm(b),'r.');
semilogy(0:it2,rv2./norm(b),'k.');
legend('No Preconditioner','IC(0)','MIC(0)');
如图所示:
以上绘图显示,修正不完全 Cholesky 预条件子创建的收敛更快。
也可以尝试通过阈值调降使用不完全 Cholesky 分解。以下绘图显示了通过各种调降容差构造预条件子时 pcg 的收敛。
L3 = ichol(A, struct('type','ict','droptol',1e-1));
[x3,fl3,rr3,it3,rv3] = pcg(A,b,tol,maxit,L3,L3');
L4 = ichol(A, struct('type','ict','droptol',1e-2));
[x4,fl4,rr4,it4,rv4] = pcg(A,b,tol,maxit,L4,L4');
L5 = ichol(A, struct('type','ict','droptol',1e-3));
[x5,fl5,rr5,it5,rv5] = pcg(A,b,tol,maxit,L5,L5');
figure
semilogy(0:maxit,rv0./norm(b),'b-','linewidth',2);
hold on
semilogy(0:it3,rv3./norm(b),'b-.','linewidth',2);
semilogy(0:it4,rv4./norm(b),'b--','linewidth',2);
semilogy(0:it5,rv5./norm(b),'b:','linewidth',2);
legend('No Preconditioner','ICT(1e-1)','ICT(1e-2)', ...
'ICT(1e-3)','Location','SouthEast');
如图所示:
请注意,通过调降容差 1e-2 构造的不完全 Cholesky 预条件子表示为 ICT(1e-2)。
与零填充的不完全 Cholesky 分解一样,阈值调降分解也会受益于修正(即 options.michol = 'on'),因为矩阵由一个椭圆偏微分方程生成。与 MIC(0) 一样,修正后的阈值调降不完全 Cholesky 分解也将保留预条件子对常向量的作用,也即 norm(A*e-L*(L'*e)) 将约为零。
使用 diagcomp
选项
此示例介绍了 ichol 的 diagcomp 选项的用法。
正定矩阵的不完全 Cholesky 分解并不总是存在。以下代码构造一个随机的对称正定矩阵,并尝试使用 pcg 对线性系统求解。
S = rng('default');
A = sprandsym(1000,1e-2,1e-4,1);
rng(S);
b = full(sum(A,2));
[x0,fl0,rr0,it0,rv0] = pcg(A,b,1e-6,100);
由于未实现收敛,因此请尝试构造一个不完全 Cholesky 预条件子。
L = ichol(A,struct('type','ict','droptol',1e-3));
Error using ichol
Encountered nonpositive pivot.
如果 ichol 按以上所示进行分解,您可以使用 diagcomp 选项构造一个偏移的不完全 Cholesky 分解。也就是说,带有对角线补偿的 ichol 会构造 L 以使 L*L' 在未显式构造 M 的情况下约为 M = A + alpha*diag(diag(A)),而不是构造 L 以使 L*L' 约为 A。由于对于对角线占优的矩阵始终可以进行不完全分解,因此可以求得使 M 在对角线上占优的 alpha。
alpha = max(sum(abs(A),2)./diag(A))-2
alpha = 62.9341
L1 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha));
[x1,fl1,rr1,it1,rv1] = pcg(A,b,1e-6,100,L1,L1');
此处的 pcg 仍然无法在所需的迭代次数内收敛至所需公差,但正如以下绘图所示,与没有预条件子相比,收敛更适用于带有此预条件子的 pcg。选择更小的 alpha 可能会有帮助。通过试验,我们可以为 alpha 设置合适的值。
alpha = .1;
L2 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha));
[x2,fl2,rr2,it2,rv2] = pcg(A,b,1e-6,100,L2,L2');
现在,pcg 将会收敛,并且绘图可以显示每个预条件子的收敛历史记录。
semilogy(0:100,rv0./norm(b),'b.');
hold on;
semilogy(0:100,rv1./norm(b),'r.');
semilogy(0:it2,rv2./norm(b),'k.');
legend('No Preconditioner','\alpha \approx 63','\alpha = .1');
xlabel('Iteration Number');
ylabel('Relative Residual');
如图所示:
提示
-
此例程提供的因子可能很有用,可用作通过 pcg 或 minres 等迭代方法求解的线性系统的预条件子。
-
ichol 仅适用于稀疏方阵。