Re: [eigen] eigen3 SelfAdjointEigenSolver<Matrix3f> still "buggy"

[ Thread Index | Date Index | More Archives ]


Thanks a lot for the patch! Things look better now, although I wonder why the results on Matrix3f are always the negative of Matrix3d :) Here's an example:
[3double]       0.21281 seconds.
  0.168841  -0.161623  -0.972302
  0.451632   0.889498 -0.0694329
  0.876083    -0.4274   0.223178
[3float]        0.123755 seconds.
-0.168841  0.161623  0.972302
-0.451632 -0.889498 0.0694328
-0.876083    0.4274 -0.223178

For matrices:
covariance_matrix_f << 0.000536227, -1.56178e-05, -9.47391e-05, -1.56178e-05, 0.000297322, -0.000148785, -9.47391e-05, -0.000148785, 9.7827e-05; covariance_matrix_d << 0.000536227, -1.56178e-05, -9.47391e-05, -1.56178e-05, 0.000297322, -0.000148785, -9.47391e-05, -0.000148785, 9.7827e-05;

I also had to apply the following patch to get rid of some compilation warnings about "n defined but not used":
  1 --- Eigen/src/Eigenvalues/Tridiagonalization.h  2010-08-31 09:02:15.000000000 -0700
  2 +++ Eigen/src/Eigenvalues/Tridiagonalization.h  2010-08-31 09:02:30.000000000 -0700
  3 @@ -433,8 +433,8 @@
  4  void ei_tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
  5  {
  6    typedef typename MatrixType::Index Index;
  7 -  Index n = mat.rows();
  8 -  ei_assert(mat.cols()==n && diag.size()==n && subdiag.size()==n-1);
  9 +  //Index n = mat.rows();
 10 +  ei_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
 11    ei_tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
 12  }


On 08/31/2010 01:14 AM, Gael Guennebaud wrote:

I'm glad to see you already made the switch :)

Regarding your issue, it has been fixed after beta1, but you can apply
the following two changesets:


or directly use the attached patch which simply contains these two changesets.

I hope they apply well to beta1 (I did not check).

If you want an even faster 3x3 eigen decomposition, you can have a
look at the file bench/eig33.cpp (in the hg repo), which compares
current Eigen's 3x3 eigensolver against a direct method solving the
underlying degree 3 polynomial. To test it:

g++ -I .. -O2 -lrt eig33.cpp -DNDEBUG&&  ./a.out

Here it is about 4 times faster though the speed of the default method
is not constant since this is an iterative method. I'm still unsure
about how to integrate it into Eigen because I don't want to make it
the default. Indeed, it might suffers from some numerical issues if
the matrix is not well conditioned. It also involves some
trigonometric functions that is not very nice. So I guess we could add
a template option to SelfAdjointEigenSolver or add a special class...
well there are many possibilities!


On Tue, Aug 31, 2010 at 5:10 AM, Radu Bogdan Rusu<rusu@xxxxxxxxxxxxxxxx>  wrote:
We've started adopting Eigen3 in ROS (, and
so far a thumbs up to the developers for doing a great job keeping most of
the API intact! We documented some of the minor changes that we had to do to
our code at the URL mentioned above, but all in all it was a walk in the
park to switch one of our larger code libraries (PCL - to Eigen3 today.

One of the Eigen2 problems however still remains: the usage of
SelfAdjointEigenSolver for Matrix3f matrices still gives a different
("incorrect") solution than for Matrix3d. This was already discussed on this
mailing list, and the patch that we adopted in ROS for Eigen2 was:

  1 --- Eigen/src/QR/Tridiagonalization.h 2010-07-22 14:45:32.000000000 -0700
  2 +++ Eigen/src/QR/Tridiagonalization.h 2010-07-22 14:46:26.000000000 -0700
  3 @@ -391,7 +391,7 @@
  4  {
  5    diag[0] = ei_real(mat(0,0));
  6    RealScalar v1norm2 = ei_abs2(mat(0,2));
  7 -  if (ei_isMuchSmallerThan(v1norm2, RealScalar(1)))
  8 +  if (v1norm2==RealScalar(0))
  9    {
  10      diag[1] = ei_real(mat(1,1));
  11      diag[2] = ei_real(mat(2,2));

This allowed us to perform the eigen decomposition and get the same values
without much precision loss, while keeping everything else as floats.

In Eigen3 however, the structure of QR changed dramatically, and the
3.0-beta1.tar.bz2 snapshot has not fixed this issue, meaning that:

Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d>  ei_symm_d;
Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f>  ei_symm_f;

would give different values on compute (matrix_{f,d}), where matrix_f =
matrix_d (minus the type).

Is there a patch that we can apply against Eigen3 beta1 that would solve
this issue?

Thank you in advance.

PS. One of the reasons why we'd like to stick to
SelfAdjointEigenSolver<Matrix3f>  is the speed. We wrote a few tests in ROS
(under for testing the performance of
different eigenvector estimation methods in Eigen. Here are some results on
an i7 laptop:
[3double]       0.212725 seconds.
[3float]        0.104199 seconds.
[4double]       0.3162 seconds.
[4float]        0.21743 seconds.
[svd_3double]   0.451361 seconds.
[svd_3float]    0.339349 seconds.
[svd_4double]   0.534939 seconds.
[svd_4float]    0.43122 seconds.

where 3double = Matrix3d, 3float = Matrix3f, 4{double,float} = 4x4 matrix
with added 0s for padding (we thought that keeping the matrix SSE aligned
would result in some speedups, but it didn't), and svd_XX = Eigen::SVD<XX>.
Overall, the good old Matrix3f seems to be the fastest.

| Radu Bogdan Rusu |

| Radu Bogdan Rusu |

Mail converted by MHonArc 2.6.19+