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

[ Thread Index | Date Index | More Archives ]

Btw, I noticed the same negation in your eigen33.cpp test program too:

for(int k=0;k<3;++k)
    evecs.col(k) = -evecs.col(k);

Again, this seems to almost never happen with Matrix3d. Just curious if this is a bug or a feature.

PS. The compilation warnings that I mentioned only happen when -DEIGEN_NO_DEBUG is defined, which makes perfect sense.


On 08/31/2010 09:17 AM, Radu Bogdan Rusu wrote:

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,
covariance_matrix_d << 0.000536227, -1.56178e-05, -9.47391e-05,
-1.56178e-05, 0.000297322, -0.000148785, -9.47391e-05, -0.000148785,

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() &&
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

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
2 +++ Eigen/src/QR/Tridiagonalization.h 2010-07-22 14:46:26.000000000
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
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
with added 0s for padding (we thought that keeping the matrix SSE
would result in some speedups, but it didn't), and 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+