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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi,

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:

a332a2b37265
1f979f8b9d99

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!

gael


On Tue, Aug 31, 2010 at 5:10 AM, Radu Bogdan Rusu <rusu@xxxxxxxxxxxxxxxx> wrote:
> We've started adopting Eigen3 in ROS (http://www.ros.org/wiki/eigen3), 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 -
> http://www.ros.org/wiki/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:
> (http://bitbucket.org/eigen/eigen/issue/149/selfadjointeigensolver-matrix3f-vs)
>
>  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));
> (https://code.ros.org/svn/ros-pkg/stacks/geometry/trunk/eigen/matrix.patch)
>
> 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;
> and
> 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 http://www.ros.org/wiki/eigen_tests) 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.
>
> Cheers,
> Radu.
> --
> | Radu Bogdan Rusu | http://rbrusu.com/
>
>
>

# HG changeset patch
# User Gael Guennebaud <g.gael@xxxxxxx>
# Date 1279291138 -7200
# Node ID a332a2b372655990e9f23e32ad0cf6c2e5b41780
# Parent  7f593359c5dc7234605feb84ae6c6c603b4fb68b
fix bad fuzzy comparison in 3x3 tridiagonalization

diff -r 7f593359c5dc -r a332a2b37265 Eigen/src/Eigenvalues/Tridiagonalization.h
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h	Fri Jul 16 14:02:20 2010 +0200
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h	Fri Jul 16 16:38:58 2010 +0200
@@ -458,7 +458,7 @@
 };
 
 /** \internal
-  * Specialization for 3x3 matrices.
+  * Specialization for 3x3 real matrices.
   * Especially useful for plane fitting.
   */
 template<typename MatrixType>
@@ -472,7 +472,7 @@
   {
     diag[0] = ei_real(mat(0,0));
     RealScalar v1norm2 = ei_abs2(mat(2,0));
-    if (ei_isMuchSmallerThan(v1norm2, RealScalar(1)))
+    if(v1norm2 == RealScalar(0) && ei_imag(mat(1,0))==RealScalar(0))
     {
       diag[1] = ei_real(mat(1,1));
       diag[2] = ei_real(mat(2,2));
@@ -495,8 +495,8 @@
       if (extractQ)
       {
         mat << 1,   0,    0,
-              0, m01,  m02,
-              0, m02, -m01;
+               0, m01,  m02,
+               0, m02, -m01;
       }
     }
   }
# HG changeset patch
# User Gael Guennebaud <g.gael@xxxxxxx>
# Date 1279297320 -7200
# Node ID 1f979f8b9d99cc9ce5550f282420e471bed4ce3c
# Parent  a332a2b372655990e9f23e32ad0cf6c2e5b41780
disable the optimized 3x3 path for complexes which was not working at all

diff -r a332a2b37265 -r 1f979f8b9d99 Eigen/src/Eigenvalues/Tridiagonalization.h
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h	Fri Jul 16 16:38:58 2010 +0200
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h	Fri Jul 16 18:22:00 2010 +0200
@@ -384,7 +384,9 @@
 }
 
 // forward declaration, implementation at the end of this file
-template<typename MatrixType, int Size=MatrixType::ColsAtCompileTime>
+template<typename MatrixType,
+         int Size=MatrixType::ColsAtCompileTime,
+         bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
 struct ei_tridiagonalization_inplace_selector;
 
 /** \brief Performs a full tridiagonalization in place
@@ -439,7 +441,7 @@
 /** \internal
   * General full tridiagonalization
   */
-template<typename MatrixType, int Size>
+template<typename MatrixType, int Size, bool IsComplex>
 struct ei_tridiagonalization_inplace_selector
 {
   typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
@@ -462,7 +464,7 @@
   * Especially useful for plane fitting.
   */
 template<typename MatrixType>
-struct ei_tridiagonalization_inplace_selector<MatrixType,3>
+struct ei_tridiagonalization_inplace_selector<MatrixType,3,false>
 {
   typedef typename MatrixType::Scalar Scalar;
   typedef typename MatrixType::RealScalar RealScalar;
@@ -470,14 +472,14 @@
   template<typename DiagonalType, typename SubDiagonalType>
   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
   {
-    diag[0] = ei_real(mat(0,0));
+    diag[0] = mat(0,0);
     RealScalar v1norm2 = ei_abs2(mat(2,0));
-    if(v1norm2 == RealScalar(0) && ei_imag(mat(1,0))==RealScalar(0))
+    if(v1norm2 == RealScalar(0))
     {
-      diag[1] = ei_real(mat(1,1));
-      diag[2] = ei_real(mat(2,2));
-      subdiag[0] = ei_real(mat(1,0));
-      subdiag[1] = ei_real(mat(2,1));
+      diag[1] = mat(1,1);
+      diag[2] = mat(2,2);
+      subdiag[0] = mat(1,0);
+      subdiag[1] = mat(2,1);
       if (extractQ)
         mat.setIdentity();
     }
@@ -485,13 +487,13 @@
     {
       RealScalar beta = ei_sqrt(ei_abs2(mat(1,0)) + v1norm2);
       RealScalar invBeta = RealScalar(1)/beta;
-      Scalar m01 = ei_conj(mat(1,0)) * invBeta;
-      Scalar m02 = ei_conj(mat(2,0)) * invBeta;
-      Scalar q = RealScalar(2)*m01*ei_conj(mat(2,1)) + m02*(mat(2,2) - mat(1,1));
-      diag[1] = ei_real(mat(1,1) + m02*q);
-      diag[2] = ei_real(mat(2,2) - m02*q);
+      Scalar m01 = mat(1,0) * invBeta;
+      Scalar m02 = mat(2,0) * invBeta;
+      Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
+      diag[1] = mat(1,1) + m02*q;
+      diag[2] = mat(2,2) - m02*q;
       subdiag[0] = beta;
-      subdiag[1] = ei_real(ei_conj(mat(2,1)) - m01 * q);
+      subdiag[1] = mat(2,1) - m01 * q;
       if (extractQ)
       {
         mat << 1,   0,    0,