New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Bug 1442 - Product of 3x3 fixed-size matrix and its inverse gives wrong result
Summary: Product of 3x3 fixed-size matrix and its inverse gives wrong result
Status: RESOLVED WONTFIX
Alias: None
Product: Eigen
Classification: Unclassified
Component: LU (show other bugs)
Version: 3.2
Hardware: x86 - 32-bit Windows
: Normal Wrong Result
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2017-06-22 07:14 UTC by ipolok
Modified: 2017-07-18 12:24 UTC (History)
3 users (show)



Attachments
Test for various sizes of matrices. Does not contain main(), (2.42 KB, text/x-csrc)
2017-06-22 07:14 UTC, ipolok
no flags Details

Description ipolok 2017-06-22 07:14:10 UTC
Created attachment 792 [details]
Test for various sizes of matrices. Does not contain main(),

When taking a product of a 3x3 fixed-size matrix (Eigen::Matrix3d) and its inverse, the output will not be an identity matrix. Looking deeper into the issue, the last row of the inverse is not computed correctly. Adding an eval() solves the issue.

    typedef Eigen::Matrix3d Matrix;
    Matrix M = Matrix::Random();
    Matrix I_bad = M * M.inverse(); // buggy
    Matrix I_good = M * M.inverse().eval(); // ok
    Matrix bad_inv = M.inverse().eval() * I_bad; // what M.inverse() evaluated to
    Matrix good_inv = M.inverse(); // what it was supposed to evaluate to

This does not seem to affect dynamic matrices. A more comprehensive test is attached (I've tested up to 100x100 fixed-size matrices, only 3x3 have the issue). I get:

M =
-0.0402539    0.64568   0.717887
 -0.299417    0.49321   0.421003
  0.791925  -0.651784  0.0270699

M * M.inverse() =
  -1.04439    5.07456    1.81471
  -1.19893    3.97596    1.06423
-0.0770894    0.19135    1.06843

the M.inverse() was
 4.19313   -7.073 -1.19833
 4.97643  -8.3002 -2.88525
-5.69557  14.1375   5.0557

real M.inverse() is
 4.19313   -7.073 -1.19833
 4.97643  -8.3002 -2.88525
-2.84779  7.06874  2.52785

M * M.inverse().eval() =
            1 -8.88178e-016  2.22045e-016
 2.22045e-016             1  2.22045e-016
 -9.4369e-016  8.88178e-016             1

This happens in Visual Studio 2008 Professional with SP1 on Windows, with Eigen 3.2.10 and only in x86 (32-bit) debug build. It is somewhat elusive, as it does not happen with Visual Studio 2008 Express on another machine. It does not seem to be an issue in Linux (tried with g++ 4.4.7, with Eigen 3.2.4, 3.2.10 and 3.3.3).
Comment 1 Gael Guennebaud 2017-06-23 12:03:48 UTC
I'm afraid this is msvc 2008 issue that won't be easy to workaround. I can only recommend to move to Eigen 3.3.4 and/or to a more recent version of MSVC. If you really want to investigate this issue, then the only way I see is to generate the assembly for

Matrix3d A = M * M.inverse();

and

Matrix3d A = M * M.inverse().eval();


and to compare them to see where MSVC screw up.
Comment 2 ipolok 2017-06-26 09:56:37 UTC
I don't think 3.3.4 will compile with MSVC 2008, unless there are major compatibility improvements since 3.3.3 (which I tried and failed to build).

Anyway, I realize 2008 is a pretty old piece of software so this is low priority. I just thought you might want to know something bad happens there. I tried to get that assembly and see what's amiss but enabling the "generate assembly" option actually fixes the bug. Duh.

I tried step-by-step debugging and unsurprisingly, the inverse is evaluated into a temp in the CoeffBasedProduct constructor, so the two lines of code should be pretty much identical. In both cases, inverse_impl::evalTo() is called from the matrix copy-constructor (which is invoked from different places but that should make no difference). I'd say this is some compiler glitch related to incremental linking.

Let's close this. Sorry this took your time.
Comment 3 Christoph Hertzberg 2017-07-18 12:24:45 UTC
We officially only support MSVC 2010 and newer --> WONTFIX

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