ichol 作为 cholinc 替代品：非正向枢轴

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的文档，算法部分的第一段）。您可以使用ilu 获得与cholinc 非常相似的因子（注意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 的文档。当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。