MATLAB中ichol函数用法

目录

语法

说明

示例

不完全 Cholesky 分解

使用 ichol 作为预条件子

使用 diagcomp 选项

---

        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 仅适用于稀疏方阵。