ichol as cholinc replacement: nonpositive pivot

Question

In Matlab 2012 the cholinc command is marked as obsolete. The warning message says it is to be replaced by ichol. Until now I used cholinc(A,droptol), typically with droptol=1E-15. In the new version I tried to use ichol(A,struct('droptol',droptol,'type','ict')), which works most of the time, but sometimes I get a warning message

Error using ichol
Encountered nonpositive pivot.

Is this a fundamental problem (i.e. a problem cholinc also had but didn't report) or is there a way to make ichol behave in the way cholinc formerly did?

Answer 1

The error indicates that the incomplete Cholesky method has broken down, which is a well-known possibility for symmetric positive-definite, but not diagonally dominant matrices. That is, even if a matrix has a (complete) Cholesky factorization, it may not have an incomplete Cholesky factorization.

cholinc doesn't suffer breakdown because it is not a true incomplete Cholesky. Rather, it calls luinc with no pivoting, throws out L and then scales the resulting U to obtain a kind of incomplete Cholesky factor (see the doc for cholinc, first paragraph of the Algorithms section). You can obtain a factor very similar to that from cholinc by using ilu (note that luinc is also obsolete).

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

Use of ichol is strongly encouraged, if possible. Note that you can use the 'diagcomp' option to try to prevent breakdown in your factorization, but finding an effective parameter alpha may require experimentation. See the doc for ichol for an example. When ichol does not breakdown, it tends to be much faster as it exploits symmetry. Also, it tends to produce sparser factors than cholinc which results in faster application of the factor as a preconditioner, which can mean faster solution times with pcg. For example,

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

To be fair, the timing of cholinc above includes time to dispatch the warning, but a tic/toc of just a single warning shows that time to be in the noise of this particular computation.

Finally, note that by default, ichol references the lower triangle of the input matrix and returns a lower triangular factor. Preferring lower triangular factors can give a notable performance increase:

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

As a final remark, the drop tolerance of 1e-15 you reference above is very small. If that is the kind of tolerance you are using, you might be better off using a complete factorization such as chol, ldl, or lu.