Gael,
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 }
Thanks,
Radu.
On 08/31/2010 01:14 AM, Gael Guennebaud wrote:
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/