ichol作为cholinc替代品:nonpositive pivot

Question

在Matlab 2012中,该cholinc命令被标记为已过时.警告消息说它将被替换为ichol.直到现在我才使用cholinc(A,droptol),通常是droptol=1E-15.在我尝试使用的新版本中ichol(A,struct('droptol',droptol,'type','ict')),大部分时间都可以使用,但有时我会收到警告信息

Error using ichol
Encountered nonpositive pivot.

这是一个根本问题(即问题cholinc是否已经但没有报告)或者是否有办法以ichol前所未有的方式行事cholinc？

Answer 1

该误差表明不完整的Cholesky方法已经破坏,这是对称正定但不是对角占优矩阵的众所周知的可能性.也就是说,即使矩阵具有(完全)Cholesky分解,它也可能没有不完整的Cholesky分解.

cholinc不会遭受破坏,因为它不是真正不完整的Cholesky.相反,它调用luinc没有旋转,抛出L然后缩放得到的U以获得一种不完整的Cholesky因子(参见文档cholinc,算法部分的第一段).您可以获得与cholinc使用非常类似的因子ilu(注意luinc也是过时的).

[L,U] = ilu(A,struct('type','ilutp','droptol',droptol,'thresh',0));
R = diag(sqrt(abs(diag(U))))\U;
% Essentially the same as cholinc.

ichol如果可能,强烈建议使用.请注意,您可以使用该'diagcomp'选项来尝试防止分解中断,但找到有效参数alpha可能需要进行实验.有关ichol示例,请参阅doc .什么ichol时候没有崩溃,它利用对称性往往会快得多.此外,它倾向于产生比较稀疏的因素cholinc导致更快地应用因子作为预处理器,这可能意味着更快的解决方案时间pcg.例如,

M = delsq(numgrid('L',200));
tic; R1 = ichol(M,struct('type','ict','droptol',1e-2,'shape','upper')); toc
% Elapsed time is 0.013809 seconds.
nnz(R1)
% ans = 145632
tic; R = cholinc(M,1e-2); toc
% Elapsed time is 0.167155 seconds.
nnz(R)
% ans = 173851

公平地说,cholinc上面的时间包括发送警告的时间,但是只有一个警告的tic/toc表示时间处于这个特定计算的噪声中.

最后,请注意,默认情况下,ichol引用输入矩阵的下三角形并返回下三角因子.优先选择较低的三角形因子可以显着提高性能:

tic; L = ichol(M,struct('type','ict','droptol',1e-2)); toc
% Elapsed time is 0.008895 seconds.

最后一点,上面提到的1e-15的跌落容差非常小.如果是那种你正在使用公差的,你可能会更好使用一个完整的因式分解,例如chol,ldl或lu.