| [eigen] May be a bug in SelfAdjointEigenSolver - EigenSolver |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
- To: eigen <eigen@xxxxxxxxxxxxxxxxxxx>
- Subject: [eigen] May be a bug in SelfAdjointEigenSolver - EigenSolver
- From: Jose Luis Blanco <joseluisblancoc@xxxxxxxxx>
- Date: Sat, 11 Dec 2010 17:26:08 +0100
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:received:date:message-id :subject:from:to:content-type; bh=5LpNEFmx9EFveP5Svl7ivgLBZDrvWYO8DYvcE/AYaZQ=; b=lciH+JcKkvG6nIORD3CBbYYbF2LVU/1I4DnugvZSyu2TWALNfY3GAz3qU0GBmoq4C/ 5it6aVo3aP6HLHN+tpDVTAe6m1VUpzuOugM99w3cSp4ARnh1GrUdc1vtvkU7YUpcPQ+B RGWEY9icopnIMqlpU0XK18132AI/PvDD9yVws=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:date:message-id:subject:from:to:content-type; b=pHAFRQqlMnAGZxCq6WOuTMKt/GmmltBpDsqy9/8BKEOknpNM6WAwBXaMP5dX/LuziH NOnRWPBcb9uoxI8JLuopVUATkUVXNcIBXyM1S86ttll4G40LoU/Ws3hc6Pa4h4ZZgt4i Ri6x8zTX3bUS9+mVnM7Iq4Q+l/RjY0l3Ll9eQ=
Hi again,
Doing unit testing for an Eigen-dependent library, I found what I
think is a bug in either Eigen... or in the virtualization system the
tests run on (Ubuntu PPA servers).
It's about eigenvalues decomposition of a real symmetric 4x4 matrix.
The unit tests (with Google Test gtests) did run on Ubuntu's PPA
building machines for these architectures:
- x86 --> all OK
- amd64 (x86_64) --> Fail
- lpia --> all OK
The problematic test case didn't fail in my machine with GCC+linux
64bit nor with MSVC+Win32. It didn't fail as well in the PPA servers
running upon x86 nor lpia, but *only in the amd64 machines*. It's been
a long debugging since I can't reproduce the fail on my local machine
(which is amd64!): it only fails on Launchpad's amd64 build servers.
The complete log for an amd64 is in [1], but the relevant part is:
======================================================================
void do_test_EigenVal4x4_sym_vs_generic_eigen() [with int ColRowOrder = 0]
SelfAdjointEigenSolver:
eigvecs:
-0.440208 -0.710821 -0.0254769 0.547998
-0.392351 0.674628 0.253235 0.571673
-0.136391 -0.177049 0.92866 -0.296044
0.796037 -0.0909078 0.26984 0.534085
eigvals:
0.975003
3.97658
7.0729
38.4097
EigenSolver:
eigvecs:
(0.509673,0) (0.752303,0) (0.299386,0) (-2.22475,0)
(0.591864,0) (-0.62003,0) (0.505669,0) (1.46388,0)
(-0.302914,0) (0.230792,0) (0.109769,0) (-1.6662,0)
(0.546052,0) (-0.0174474,0) (-0.803457,0) (0.0737898,0)
eigvals:
(37.5398,0)
(5.03294,0)
(1.2362,0)
(6.62521,0)
void do_test_EigenVal4x4_sym_vs_generic_eigen() [with int ColRowOrder = 1]
SelfAdjointEigenSolver:
eigvecs:
0.587602 -0.324628 -0.516793 0.531287
0.0488521 0.776559 0.184898 0.600318
-0.325169 0.38677 -0.835326 -0.216576
-0.739325 -0.376804 -0.0311285 0.557178
eigvals:
0.612734
2.11465
10.7396
36.9672
EigenSolver:
eigvecs:
(0.55922,0) (-0.309765,0) (0.175508,0) (-0.366926,0)
(0.585747,0) (0.432738,0) (-0.839941,0) (-0.40577,0)
(-0.207763,0) (0.80562,0) (-2.89354,0) (0.172306,0)
(0.548642,0) (0.26699,0) (-0.858551,0) (0.877879,0)
eigvals:
(37.5398,0)
(5.03294,0)
(6.62521,0)
(1.2362,0)
====================================================================
The code for generating that is attached. For debugging, the unit test
was modified to just show information to the console so I can show you
these results from the server buildlogs.
My conclusions:
1) I would expect the eigenvalues & eigenvectors of a *real symetric*
matrix to be the same no matter they're computed with
SelfAdjointEigenSolver or EigenSolver (right???). I mean, disregarding
the order of the eigenvalues, that's not the point. You can see in the
logs that the results are different for both solvers!
2) The results for the SelfAdjointEigenSolver aren't the same for
RowMajor vs ColMajor. They are the same for the EigenSolver.
3) From the 4 combinations (SelfAdjointEigenSolver vs EigenSolver,
ColMajor vs RowMajor), the only one giving the same results than
MATLAB/octave is SelfAdjointEigenSolver + ColMajor.
What do you think? Is there anything wrong in my code?
I could just ignore this and remove the failing test cases from my
library if you think it's a bug in the virtualization platforms...
More info:
- The version of Eigen is from 9/DEC/2010.
- I also tested 2x2 and 3x3 eigen-decompositions and all work fine on
all architectures. That's what made me think on a possible bug related
to special SSE optimizations for 4x4 matrices (makes sense??)
Regards,
JL
[1] http://launchpadlibrarian.net/60484561/buildlog_ubuntu-maverick-amd64.mrpt_1:0.9.3svn2312maverick-1~ppa1~maverick_FAILEDTOBUILD.txt.gz
/* +---------------------------------------------------------------------------+
| The Mobile Robot Programming Toolkit (MRPT) C++ library |
| |
| http://mrpt.sourceforge.net/ |
| |
| Copyright (C) 2005-2010 University of Malaga |
| |
| This software was written by the Machine Perception and Intelligent |
| Robotics Lab, University of Malaga (Spain). |
| Contact: Jose-Luis Blanco <jlblanco@xxxxxxxxxxxx> |
| |
| This file is part of the MRPT project. |
| |
| MRPT is free software: you can redistribute it and/or modify |
| it under the terms of the GNU General Public License as published by |
| the Free Software Foundation, either version 3 of the License, or |
| (at your option) any later version. |
| |
| MRPT is distributed in the hope that it will be useful, |
| but WITHOUT ANY WARRANTY; without even the implied warranty of |
| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| GNU General Public License for more details. |
| |
| You should have received a copy of the GNU General Public License |
| along with MRPT. If not, see <http://www.gnu.org/licenses/>. |
| |
+---------------------------------------------------------------------------+ */
#include <Eigen/Dense>
#include <gtest/gtest.h>
using namespace Eigen;
using namespace std;
/** A macro for obtaining the name of the current function: */
#if defined(__BORLANDC__)
#define __CURRENT_FUNCTION_NAME__ __FUNC__
#elif defined(_MSC_VER) && (_MSC_VER>=1300)
#define __CURRENT_FUNCTION_NAME__ __FUNCTION__
#elif defined(_MSC_VER) && (_MSC_VER<1300)
// Visual C++ 6 HAS NOT A __FUNCTION__ equivalent.
#define __CURRENT_FUNCTION_NAME__ ::system::extractFileName(__FILE__).c_str()
#else
#define __CURRENT_FUNCTION_NAME__ __PRETTY_FUNCTION__
#endif
template <int ColRowOrder>
void do_test_EigenVal4x4_sym_vs_generic_eigen()
{
typedef Matrix<double,4,4,ColRowOrder> Mat44;
const double dat_C1[] = {
13.737245,10.248641,-5.839599,11.108320,
10.248641,14.966139,-5.259922,11.662222,
-5.839599,-5.259922,9.608822,-4.342505,
11.108320,11.662222,-4.342505,12.121940 };
const Mat44 C1(dat_C1); // It doesn't mind the row/col major order since data are symetric
// Symetric --------------------
// This solver returns the eigenvectors already sorted.
Eigen::SelfAdjointEigenSolver<Mat44> eigensolver(C1);
// MatrixXd eVecs_s = eigensolver.eigenvectors();
// MatrixXd eVals_s = eigensolver.eigenvalues();
cout << endl << __CURRENT_FUNCTION_NAME__ << endl
<< "SelfAdjointEigenSolver:\n"
<< "eigvecs: " << endl << eigensolver.eigenvectors() << endl
<< "eigvals: " << endl << eigensolver.eigenvalues() << endl;
// Generic ---------------------
Eigen::EigenSolver<Mat44> es(C1, true);
// MatrixXd eVecs_g = es.eigenvectors().real();
// MatrixXd eVals_g = es.eigenvalues().real();
cout << endl
<< "EigenSolver:\n"
<< "eigvecs: " << endl << es.eigenvectors() << endl
<< "eigvals: " << endl << es.eigenvalues() << endl;
}
// Compare the two ways of computing matrix eigenvectors: generic & for symmetric matrices:
TEST(MatricesEigen,EigenVal4x4_sym_vs_generic)
{
do_test_EigenVal4x4_sym_vs_generic_eigen<Eigen::ColMajor>();
do_test_EigenVal4x4_sym_vs_generic_eigen<Eigen::RowMajor>();
}