Bug 567 - ConjugateGradient returns nans for zero right-hand side
ConjugateGradient returns nans for zero right-hand side
Status: RESOLVED FIXED
Product: Eigen
Classification: Unclassified
Component: Sparse
unspecified
All All
: Normal major
Assigned To: Nobody
:
Depends on:
Blocks:
  Show dependency treegraph
 
Reported: 2013-03-18 18:01 UTC by Clemens Hofreither
Modified: 2013-03-20 16:29 UTC (History)
2 users (show)



Attachments

Description Clemens Hofreither 2013-03-18 18:01:43 UTC
The following code:

////////////////////////////////

#include <iostream>
#include "Eigen/Dense"
#include "Eigen/Sparse"

using namespace Eigen;
using namespace std;

int main()
{
    SparseMatrix<double> A(3, 3);
    for (int i = 0; i < 3; ++i)
        A.insert(i,i) = 1.0;

    VectorXd b = VectorXd::Zero(3);

    VectorXd x = ConjugateGradient< SparseMatrix<double> >().compute( A ).solve( b );

    cout << x << endl;
}

////////////////////////////////

produces the following output:

-nan
-nan
-nan

Reason: the CG implementation doesn't check if the initial residual is already zero and then divides by zero. I think at least some other iterative solvers suffer from the same problem.
Comment 1 Desire NUENTSA 2013-03-20 11:24:07 UTC
Thanks for the report. 
I push a fix in the devel branch
https://bitbucket.org/eigen/eigen/commits/a0af9c37220d/
changeset:   a0af9c37220d
user:        dnuentsa
date:        2013-03-20 11:22:45
summary:     Handle zero right hand side in CG and GMRES
affected #:  2 files
Comment 2 Clemens Hofreither 2013-03-20 11:38:03 UTC
Thank you, looks good!

Can you also check BiCGStab, I think it has the same issue?

A more subtle issue is that if a (non-zero) guess is given and it already solves the linear problem, then the same problems will happen.
Comment 3 Desire NUENTSA 2013-03-20 16:19:07 UTC
Thanks, 
The new commit address the subtle issue as well :
https://bitbucket.org/eigen/eigen/commits/38b8e322e272/
changeset:   38b8e322e272
user:        dnuentsa
date:        2013-03-20 16:15:18
summary:     Bug567 : Fix iterative solvers to immediately return when the initial guess is the true solution and for trivial solution
affected #:  4 files

Let us know if it fixes the issue.

Note You need to log in before you can comment on or make changes to this bug.