Re: [eigen] Issue in SelfAdjointEigenSolver for 3-by-3 matrices

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2010/7/15 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>:
> Hi all,
>
> There seems to be an issue with SelfAdjointEigenSolver when applied to
> 3-by-3 matrices, because the unit test eigensolver_selfadjoint_1 fails
> occasionally on my computer.
>
> I think I boiled it down to this example:
>
> ------------------------------------------------
>
> #include <iostream>
> #include <Eigen/Core>
> #include <Eigen/Eigenvalues>
>
> using namespace Eigen;
> using namespace std;
>
> int main()
> {
>  Matrix3f A = Matrix3f::Identity(3,3);
>  A(0,2) = 0.001;
>  SelfAdjointEigenSolver<Matrix3f> eigOfA(A);
>  cout << eigOfA.eigenvalues() << endl;
> }
>
> ------------------------------------------------
>
> The output of this program is
>
>  1
>  1
>  1
>
> This program computes the eigenvalues of the matrix
>
>  1      0     0.001
>  0      1     0
>  0.001  0     1
>
> They are easily computed by hand; the exact eigenvalues are 0.999, 1, and
> 1.001. I think this is rather too big an error.
>
> The issue seems to be in
>
>  ei_tridiagonalization_inplace_selector<MatrixType,3>::run()
>
> in Tridiagonalization.h, where we have
>
>    RealScalar v1norm2 = ei_abs2(mat(2,0));
>    if (ei_isMuchSmallerThan(v1norm2, RealScalar(1)))
>    {

I dont have time now, but I can already tell you that this code is out
of place in a decomposition. It uses the default "dummy" precision
which for float is 1e-5f, so it considers v1 to be "zero" if v1norm2 <
1e-5f, which means it considers v1 to be "zero" if its norm is less
than sqrt(1e-5) which is 3.1e-3, which is big enough to account for
the inaccuracy you are reporting.

So this is doubly wrong:
 - don't use the default dummy precision when doing fuzzy compares in
a decomposition.
 - check for formulas that amount to replacing a precision value by
its sqrt, which is much bigger.

Benoit

>        ...
>
> If the condition in the if statement is true, then mat(2,0) is simply
> ignored. This does not seem the right condition, and I'm not even sure why
> we need two branches here. Can somebody please help me here? I'd rather not
> touch code I don't understand ...
>
> This is all old code and the same issue is present in Eigen 2, so perhaps we
> want to backport the fix before the next version is released?
>
> Cheers,
> Jitse
>
>
>