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

[Manoj Rajagopalan <rmanoj@xxxxxxxxx>]

Just an observation: you probably meant A(2,0) = 0.001;
The documentation says only the lower triangular part is referenced.

Anyway, I can confirm what you have noticed with this correction.

-- Manoj


On Thursday 15 July 2010 03:59:09 pm Jitse Niesen wrote:
> 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)))
>      {
>  	...
>
> 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