[vtkusers] Jacobi eigenvectors incorrect?

[David Doria]

On Wed, Oct 28, 2009 at 5:57 AM, Fred Fred <stan1313 at hotmail.fr> wrote:
> BTW, which is the easiest way to compute the real eigenvalues of a
> non-symetric matrix in VTK? Of course such eigenvalues may not exist but in
> case they exist, I wonder which is the better way to use.
>
>> Date: Tue, 27 Oct 2009 08:39:34 -0400
>> From: bill.lorensen at gmail.com
>> To: daviddoria+vtk at gmail.com
>> CC: vtkusers at vtk.org
>> Subject: Re: [vtkusers] Jacobi eigenvectors incorrect?
>>
>> According to the documentation in vtkMath.cxx, the code
>> vtkJacobiN finds the solution of an nxn real SYMMETRIC matrix.
>> Can you try your tests with a symmetric matrix and compare the results
>> between vtk and matlab.
>>
>> Bill

I am still not getting the same evecs even with a symmetric matrix:

#include "vtkMath.h"

int main (int argc, char *argv[])
{
  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;
  }

  /*
  1 2 3
  2 5 9
  3 9 8
  */
  a0[0] = 1;
  a0[1] = 2;
  a0[2] = 3;
  a1[0] = 2;
  a1[1] = 5;
  a1[2] = 9;
  a2[0] = 3;
  a2[1] = 9;
  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;
  return 0;
}

Here is the matlab:


a=[1 2 3; 2 5 9; 3 9 8];
[v d] = eig(a);

The columns of V are the eigenvectors of a:

v =

    0.1372    0.9644    0.2259
    0.7382   -0.2516    0.6259
   -0.6605   -0.0809    0.7464

The output of my vtk c++:
Evec corresponding to smallest eval: 0.746449 0.0808535 -0.660513

I see this vector reversed in the bottom row of the matlab eigenvector
matrix - is this a coincidence or am I misinterpreting the output from
the Jacobi function and just getting these things in the wrong order?
This is exactly what I was talking about when I suggested that we
either reimplement or make a nice wrapper so that users can be sure
what is going on - setting (r,c) type entries in a "matrix" and then
getting "vectors" back instead of ** style arrays.

Thanks,

David