Hello, I've got a question about IncompleteCholesky (I know it's "unsupported," but a bunch of updates were made to it recently, which gave me hope). I'm having a weird issue with it. As I understand it, incomplete cholesky as an algorithm should, given a matrix A, produce an approximation to L where LL^T=A. I would expect the "solve" function (or _solve_impl in this case) to give an approximation to A^(-1). However, when I test IncompleteCholesky on a trivial example, _solve_impl seems to be equivalent to multiplying by A rather than by A^(-1). Is this by design?

Regards,
-Avi

Example code:

Code: Select all

    SparseMatrix<T> testmatrix(4,4);
    std::vector<Triplet<T>> triplets;
    for(int i=0;i<4;i++){
        triplets.push_back(Triplet<T>(i,i,4));
    }
    testmatrix.setFromTriplets(triplets.begin(),triplets.end());
    IncompleteCholesky<T> testprecond;
    testprecond.compute(testmatrix);
    Matrix<T,4,1> realsimple,other;
    realsimple=Matrix<T,4,1>::Constant(1);
    testprecond._solve_impl(realsimple,other);
    std::cout<<"Test: "<<other.transpose()<<std::endl;


Okay, I may have some insight on this. There appear to be two issues:

First, in the function _solve_impl, both of the places where the line

Code: Select all

x = m_scale.asDiagonal() * x;

appear, it should instead be something equivalent to 1/m_scale

Second, the paper that's being implemented calls for the inner factorization algorithm to be tried with an initial shift value and then re-run with increasingly large shift values if the factorization fails (the update is shift=max(2*shift, initialShift)). I implemented a somewhat clunky version of this and it seems to make things work; without it I was getting factorization failures.

Should I clean this up and submit it somehow?

Regards,
-Avi

First you'll need to file the bug on the bugzilla page at http://eigen.tuxfamily.org/bz/
List any info that is necessary to recreate the original problem, usually just the version of Eigen, the compiler, and which operating system you are running.
Upload your patch that includes the proposed fixes and the code you used to test the new patch.
Then wait until someone can replicate the bug to get it marked for review.
The time for a fix to get pushed to the main repository will vary.
If you send me a PM with a link to the bugzilla report, I'll look it over and at least try to replicate the bug.

arobinsonmosher wrote:First, in the function _solve_impl, both of the places where the line

Code: Select all

x = m_scale.asDiagonal() * x;

appear, it should instead be something equivalent to 1/m_scale

Right!

Code: Select all

x = m_scale.asDiagonal().inverse() * x;

Second, the paper that's being implemented calls for the inner factorization algorithm to be tried with an initial shift value and then re-run with increasingly large shift values if the factorization fails (the update is shift=max(2*shift, initialShift)). I implemented a somewhat clunky version of this and it seems to make things work; without it I was getting factorization failures.

Should I clean this up and submit it somehow?

Sure, please send a patch on the bugzilla or a pull-request on bitbucket. We should probably also limit the number of tries, for instance, if the entry contains NaN or inf, or if we introduce overflow because of the shift, we must avoid to run into an infinite loop.

I have submitted a bug report with a patch at http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1076

I actually left out my re-try code; it never seemed to generate a factorization that gave good results and I'm not confident enough in my understanding of the algorithm to say whether it's a good idea or not. I did change the response to finding a negative value on the diagonal so that it was a bit better-behaved (avoided printing out to cerr unconditionally and made it possible to use the info() function to check whether it succeeded).