Incomplete Cholesky factorization algorithm

Question

I want to implement incomplete Cholesky factorization to precondition, the algorithm I refer from incomplete Cholesky factorization,

function a = ichol(a)
    n = size(a,1);

    for k=1:n
        a(k,k) = sqrt(a(k,k));
        for i=(k+1):n
            if (a(i,k)!=0)
                a(i,k) = a(i,k)/a(k,k);            
            endif
        endfor
        for j=(k+1):n
            for i=j:n
                if (a(i,j)!=0)
                    a(i,j) = a(i,j)-a(i,k)*a(j,k);  
                endif
            endfor
        endfor
    endfor

    for i=1:n
        for j=i+1:n
            a(i,j) = 0;
        endfor
    endfor            
endfunction

but it seems not right.

The problem is the algorithm cannot guarantee that a[k,k] is a positive integer. When a[k,k] is not a positive integer, the algorithm will incur an error.

I want to ask for the correct version of Cholesky decomposition and some relevant references.

Answer 1

Cholesky decomposition is applicable to positive-definite matrices (for positive-semidefinite the decomposition exists, but is not unique). The positive-definiteness, is what ensures that a[k,k] is a positive number and sqrt is ok (see, for example, a Wiki explanation on that).

A similar story is happening with an incomplete Cholesky factorization, as its applicability is also limited to positive-definite matrices.

The fact that you encounter a problem is due to the fact that the matrix you are trying to factorize does not possess the required properties for Cholesky factorizations family (most likely) or the matrix loses its positive-definiteness during the numerical procedure (very unlikely).

Thus, the code above for the incomplete Cholesky factorization is correct.

In your situation, I would suggest looking into, say, incomplete LU factorization, which is a bit heavier but has fewer expectations of the matrices it is applied to. It is a common factorization to be used in many types of preconditioners, which involve matrices that are not positive-definite.

Alternatively, if you can do changes to the matrix you are factorizing (say, you are working on a precondition) that will enforce its positive-definiteness, it may be an even better choice.

Answer 2

When fed a symmetric positive definite matrix, the full Cholesky factorization is guaranteed to produce a valid real-valued result, issues of finite precision aside. Such issues are rare due to its stability, which is part of why it's so popular.

One likely reason you're not seeing this behavior in your tests is that you're NOT running the full Cholesky factorization. Incomplete Cholesky just doesn't have this guarantee.

This is perhaps the biggest weakness of the incomplete factorization: in order to reliably use it (breakdown pretty much guaranteed on problems that aren't tiny), some heuristic strategy to deal with this is necessary.

So there are a few well-researched ways to deal with this, including a positive global diagonal shift and something called the Jennings-Malik strategy. For lots of great info, see On positive semidefinite modification schemes for incomplete Cholesky factorization (SIAM 2014).