New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 790 - JacobiSVD fails on small rank-deficient float matrices (U is not unitary)
Summary: JacobiSVD fails on small rank-deficient float matrices (U is not unitary)
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: SVD (show other bugs)
Version: 3.2
Hardware: x86 - 64-bit Linux
: Normal major
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks: 3.3
  Show dependency treegraph
 
Reported: 2014-04-10 16:05 UTC by Charles Karney
Modified: 2014-04-14 13:53 UTC (History)
3 users (show)



Attachments

Description Charles Karney 2014-04-10 16:05:01 UTC
(Found during testing of the new alignPoints function.  This bug will
also cause umeyama to fail.)

JacobiSVD can fail spectacularly with small rank-deficient float
matrices.  Here's an example (Eigen 3.2.1, Linux 64-bit, g++ (GCC)
4.7.2 20121109, optimize or debug):

    Matrix<float,4,4> M;
    M << 
      0.030851028,0.275531113,-0.451450854,-0.126374155,
      -0.015783134,-0.140959471,0.230958566,0.064651995,
      0.051904849,0.463563203,-0.759536683,-0.212616309,
      0.037682864,0.336546361,-0.551422834,-0.154359206;
    JacobiSVD< Matrix<float,4,4> > svd(M, ComputeFullU | ComputeFullV);
    std::cout << "U\n" << svd.matrixU() << "\n\n";
    std::cout << "U.determinant() " << svd.matrixU().determinant() << "\n\n";

This prints out

    U
       -0.423173    -0.740408 -1.07744e-11            0
        0.216492     0.378787  5.51207e-12            0
       -0.711961     0.555262  8.08013e-12            0
       -0.516883            0            0            0

    U.determinant() 0

JacobiSVD does OK on M.transpose() and if float is promoted to double.
I realize that there are limits to what I can expect from floats.  But
it would be nice if the basic properties of the SVD were more accurately
met (specifically U should be nearly unitary).
Comment 1 Benoit Jacob 2014-04-10 19:13:40 UTC
That's very embarrassing; I would be very thankful for any debugging here!
Comment 2 Jitse Niesen 2014-04-10 19:19:48 UTC
I can reproduce with Eigen 3.2.1. The current development branch handles this particular example fine because of changeset 9034caa3d917 (Add scaling in JacobiSVD to avoid overflows), which was introduced to handle the example in http://forum.kde.org/viewtopic.php?f=74&t=118446 . However, there does not seem any problem with bad scaling in this example, so there still may be an issue in our current implementation (therefore I'm leaving the bug open pending further investigation).
Comment 3 Charles Karney 2014-04-14 00:23:25 UTC
Yes, there is definitely still an issue!  Consider the following 5x5 matrix

    const int m = 5;
    Matrix<float,m,m> M;
    M << 
      -2,2,4,-1,0,
      0,0,0,0,0,
      2,-2,-4,1,0,
      0,0,0,0,0,
      2,-2,-4,1,0;
    JacobiSVD< Matrix<float,m,m>, FullPivHouseholderQRPreconditioner >
      svd(M, ComputeFullU | ComputeFullV);
    std::cout << "U\n" << svd.matrixU() << "\n\n";
    std::cout << "U.determinant() " << svd.matrixU().determinant() << "\n\n";

With the latest development version 3a35bcba92f4, I get

U
    0.707107     0.707107 -1.58051e-08            0            0
          -0           -0            0            1            0
   -0.707107     0.707107 -2.63418e-08            0            0
-1.77636e-15  5.96046e-08 -1.77636e-15            0            0
          -0           -0            0            0            0

U.determinant() -0

Finally, I recommend that you include tests of rank deficient matrices to
the test suite (especially float ones).
Comment 4 Gael Guennebaud 2014-04-14 13:53:04 UTC
https://bitbucket.org/eigen/eigen/commits/f8b458d46a3a/
Changeset:   f8b458d46a3a
User:        ggael
Date:        2014-04-14 13:52:16
Summary:     Bug 790: fix overflow in real_2x2_jacobi_svd

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