[eigen] JacobiSVD with nearly rank-deficient matrix

[Bill Greene <w.h.greene@xxxxxxxxx>]

I've  been using JacobiSVD from a recent Eigen3 build to
solve a system of equations that may be rank deficient.

In most cases this works fine but, in some cases, like the
example matrix below-- that is very nearly a rank-2 matrix--
it does not. The singular values for this matrix are
1.41421
1.41421
3.67427e-038
I like the fact that JacobiSVD returns all singular values
including the very small one. However, I don't like
the fact that the solve() method returns
1
-7.19424e+018
-1
I would like to be able to specify a zero-tolerance value
to solve() to exclude the third singular value but this does
not seem to be possible-- at least I didn't see how to do this.

Lapack method DGELSD() includes the RCOND tolerance
and if this is set to, say, 1.0e-8,
returns the expected solution to this system.

Thanks.

Bill Greene
----------------------------------------------------------------------------------------------------
void svdFailure()
{
  Eigen::MatrixXd A(4, 3);
  A <<  1.000e+000,   1.527e-018,  -7.216e-018,
       -0.000e+000,   0.000e+000,  -1.000e+000,
        0.000e+000,   0.000e+000,  -1.000e+000,
        1.000e+000,   1.388e-018,  -9.021e-018;
 
  std::cout << "A\n" << A << std::endl;
  Eigen::VectorXd b = Eigen::VectorXd::Ones(4);
  Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
  std::cout << "singular vals\n" << svd.singularValues() << std::endl;
  Eigen::VectorXd soln = svd.solve(b);
  std::cout << "soln\n" << soln << std::endl;
}