[vtkusers] Jacobi eigenvectors incorrect?

[David Doria]

I am getting quite different eigenvectors from vtkMath::Jacobi than I
am from matlab. I have tried some physical situations (scatter matrix
eigenvalues to indicate the normal direction of the best fit plane of
some points) and the matlab vectors seem to be much better.

Here is a demo:

  double *a[3];
  double a0[3];
  double a1[3];
  double a2[3];
  a[0] = a0; a[1] = a1; a[2] = a2;

  //set the matrix to all zeros
  for(unsigned int i = 0; i < 3; i++)
  {
    a0[i] = a1[i] = a2[i] = 0.0;
  }

  /* Want this matrix
  1 2 3
  1 3 2
  3 4 8
  */
//not sure if this is expecting column or row major order, so I tried
both and neither match.
  a0[0] = 1;
  a0[1] = 2;
  a0[2] = 3;
  a1[0] = 1;
  a1[1] = 3;
  a1[2] = 2;
  a2[0] = 3;
  a2[1] = 4;
  a2[2] = 8;

  /*
  a0[0] = 1;
  a0[1] = 1;
  a0[2] = 3;
  a1[0] = 2;
  a1[1] = 3;
  a1[2] = 4;
  a2[0] = 3;
  a2[1] = 2;
  a2[2] = 8;
  */
  // Extract eigenvectors from covariance matrix
  double *v[3];
  double v0[3];
  double v1[3];
  double v2[3];

  v[0] = v0; v[1] = v1; v[2] = v2;
  double eigval[3];
  vtkMath::Jacobi(a,eigval,v);
  vtkstd::cout << "eigvals: " << eigval[0] << " " << eigval[1] << " "
<< eigval[2] << vtkstd::endl;
  vtkstd::cout << "v0: " << v0[0] << " " << v0[1] << " " << v0[2] <<
vtkstd::endl;
  vtkstd::cout << "v1: " << v1[0] << " " << v1[1] << " " << v1[2] <<
vtkstd::endl;
  vtkstd::cout << "v2: " << v2[0] << " " << v2[1] << " " << v2[2] <<
vtkstd::endl;

  vtkstd::cout << "Evec corresponding to smallest eval: " << v2[0] <<
" " << v2[1] << " " << v2[2] << vtkstd::endl;

Of course these vectors can be scaled etc, but the ratios don't seem
to match. Can anyone confirm these vectors are the same/different with
a different eigen value function from another toolkit and let me know
if they match what Jacobi gives?

Also, this code gets VERY messy/hard to read in a big hurry with all
of the ** type variables that are required to call these functions.
Is there a way to provide a cleaner interface to this? Maybe a
vtkMatrix class to provide GetColumn, a multiply operator, etc so this
type of operation doesn't have to be manually coded each time one
wants to simply use an eigenvalue function? And so the column/row
major order is indicated clearly?

Thanks,

David