ichol
不完全乔列斯基分解
说明
示例
不完全乔列斯基分解
此示例生成不完全乔列斯基分解。
从一个对称正定矩阵 A
开始:
N = 100;
A = delsq(numgrid('S',N));
A
是带有狄利克雷边界条件的 100×100 方形网格上的二维、五点离散负拉普拉斯矩阵。A
的大小为 98*98 = 9604
(并非 10000,因为网格边框用于施加狄利克雷条件)。
无填充的不完全乔列斯基分解是指在 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
也可用于通过阈值调降生成不完全乔列斯基分解。当调降容差减小时,因子往往会变得更密集,而 L*L'
乘积往往更好地接近 A
。以下绘图显示了不完全分解的相对误差对调降容差的图,以及不完全因子密度与完全乔列斯基因子密度之比。
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
作为预条件子
此示例说明如何使用不完全乔列斯基分解作为预条件子来提高收敛。
创建一个对称正定矩阵 A
。
N = 100;
A = delsq(numgrid('S',N));
创建一个不完全乔列斯基分解,作为 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
未在允许的最大迭代次数中使相对残差趋向于请求的公差。请尝试使用无填充的不完全乔列斯基分解作为预条件子。
L1 = ichol(A); [x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L1,L1');
fl1 = 0
,指示 pcg
收敛至请求的公差,经过 59 次迭代(it1
值)后实现收敛。但是,由于此矩阵是一个离散的拉普拉斯矩阵,因此使用修正不完全乔列斯基分解可创建一个更好的预条件子。修正不完全乔列斯基分解会构造一个近似的分解,该分解会保留算子对常向量的作用。也就是说,对于 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)');
以上绘图显示,修正不完全乔列斯基预条件子创建的收敛更快。
您也可以尝试通过阈值调降使用不完全乔列斯基分解。以下绘图显示了通过各种调降容差构造预条件子时 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
构造的不完全乔列斯基预条件子表示为 ICT(1e-2)
。
与零填充的不完全乔列斯基分解一样,阈值调降分解也会受益于修正(即 options.michol = 'on'
),因为矩阵由一个椭圆偏微分方程生成。与 MIC(0)
一样,修正后的阈值调降不完全乔列斯基分解也将保留预条件子对常向量的作用,也即 norm(A*e-L*(L'*e))
将约为零。
使用 diagcomp
选项
此示例介绍了 ichol
的 diagcomp
选项的用法。
正定矩阵的不完全乔列斯基分解并不总是存在。以下代码构造一个随机的对称正定矩阵,并尝试使用 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);
由于未实现收敛,因此请尝试构造一个不完全乔列斯基预条件子。
L = ichol(A,struct('type','ict','droptol',1e-3));
Error using ichol Encountered nonpositive pivot.
如果 ichol
按以上所示进行分解,您可以使用 diagcomp
选项构造一个偏移的不完全乔列斯基分解。也就是说,带有对角线补偿的 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');
输入参数
A
— 输入矩阵
稀疏方阵
输入矩阵,指定为稀疏方阵。
options
— Cholesky 分解选项
结构体
Cholesky 分解选项,指定为包含最多五个字段的结构体:
字段名称 | 摘要 | 描述 |
---|---|---|
type | 分解的类型 | 指示要执行的不完全乔列斯基分解类型。此字段的有效值为 |
droptol | 类型为 'ict' 时的调降容差 | 执行 ICT 时用作调降容差的非负标量。如果元素的模小于局部调降容差,它将从生成的因子中删除,但对角线元素除外,该元素永不会被删除。分解的第 |
michol | 指示是否执行修正的不完全乔列斯基分解 | 指示是否执行修正的不完全乔列斯基分解 (MIC)。该字段可能为 |
diagcomp | 使用指定的系数执行补偿的不完全乔列斯基分解 | 构造不完全乔列斯基因子时用作全局对角线偏移量 |
shape | 确定引用并返回的三角矩阵 | 有效值为 |
示例: options.type = 'nofill'; options.michol = 'on'; L = ichol(A,options);
参考
[1] Saad, Yousef. “Preconditioning Techniques.” Iterative Methods for Sparse Linear Systems. PWS Publishing Company, 1996.
[2] Manteuffel, T.A. “An incomplete factorization technique for positive definite linear systems.” Math. Comput. 34, 473–497, 1980.
扩展功能
基于线程的环境
使用 MATLAB® backgroundPool
在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool
加快代码运行速度。
此函数完全支持基于线程的环境。有关详细信息,请参阅在基于线程的环境中运行 MATLAB 函数。
分布式数组
使用 Parallel Computing Toolbox™ 在集群的组合内存中对大型数组进行分区。
用法说明和限制:
如果指定
options
的type
字段,则它必须为'nofill'
。A
必须为实对称矩阵或埃尔米特复矩阵。
有关详细信息,请参阅Run MATLAB Functions with Distributed Arrays (Parallel Computing Toolbox)。
版本历史记录
在 R2011a 中推出
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)