[vtkusers] Jacobi eigenvectors incorrect?

[Bryn Lloyd]

I would expect the matrix to be in row-major order, since this is 
standard for C/C++.


David Doria wrote:
> 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
> _______________________________________________
> Powered by www.kitware.com
> 
> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
> 
> Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
> 
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
> 


-- 
-------------------------------------------------
Bryn Lloyd
Computer Vision Laboratory
ETH Zürich, Sternwartstrasse 7, ETF C110
CH - 8092 Zürich, Switzerland
Tel: +41 44 63 26668
Fax: +41 44 63 21199
-------------------------------------------------