Re: [eigen] JacobiSVD with nearly rank-deficient matrix

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Hi,

It is very important to stress that our JacobiSVD uses an extremely
reliable and precise, two-sided-Jacobi algorithm. It is mathematically
proven. By contrast, LAPACK's DGELSD uses a much faster but less
accurate and less reliable bidiagonalizing SVD algorithm.

To show that Eigen's result is the correct one, I've expanded your
test case a little to compute relative error and make it compilable:


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

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;
  std::cout << "A * soln\n" << A*soln << std::endl;
  std::cout << "relative error: " << (A*soln - b).norm() / b.norm() <<
std::endl;
}

int main() { svdFailure(); }



I get this output:


A
         1  1.527e-18 -7.216e-18
        -0          0         -1
         0          0         -1
         1  1.388e-18 -9.021e-18
singular vals
1.41421
1.41421
3.67427e-38
soln
1
-1024
-1
A * soln
1
1
1
1
relative error: 1.31715e-15


A relative error of 1e-15, when the machine epsilon is 1e-16, is very
acceptable. So this solution that we return, is actually very good,
despite the condition number being absolutely awful.

If you want to test the condition number explicitly, this is already
very easy to do with Eigen, from the singularValues().

I understand that LAPACK cuts off very small singular values, but I
don't think that we should do the same in Eigen. The algorithm that
LAPACK uses makes these small singular values very unreliable, but the
same is not true with our algorithm.

There are many different notions of "goodness/badness" tests of
solutions, depending on your use cases, and the path we've chosen in
Eigen is to let the user perform his own tests, like in the example
program above where I compute the relative error on the right-hand
side.

Cheers,
Benoit


2011/2/12 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;
> }
>