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

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

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


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


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?


Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/