[vtkusers] Jacobi eigenvectors incorrect?

[Bryn Lloyd]

Bill,

Sorry, I missed that in the documentation. I tested it again for 
symmetric some nxn matrices:

--> The result of vtkMath::Jacobi is the same as in Matlab.





Bill Lorensen wrote:
> 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
> 
> On Tue, Oct 27, 2009 at 7:56 AM, David Doria <daviddoria+vtk at gmail.com> wrote:
>> On Tue, Oct 27, 2009 at 4:18 AM, Bryn Lloyd <blloyd at vision.ee.ethz.ch> wrote:
>>> Hi David,
>>>
>>> I agree. The result is different to what I get in Matlab for a different
>>> matrix. vtkMath::JacobiN returns the same result as vtkMath::Jacobi.
>>>
>>>
>>> My test programm:
>>>
>>>
>>>
>>>
>>> #include <vtkMath.h>
>>>
>>>
>>> template<class TReal>
>>> TReal **create_matrix(long nrow, long ncol) {
>>>  typedef TReal* TRealPointer;
>>>  TReal **m = new TRealPointer[nrow];
>>>
>>>  TReal* block = (TReal*)calloc(nrow*ncol, sizeof(TReal));
>>>  m[0] = block;
>>>  for ( int row = 1; row < nrow; ++row ) {
>>>    m[ row ] = &block[ row * ncol ];
>>>  }
>>>  return m;
>>> }
>>> /* free a TReal matrix allocated with matrix() */
>>> template<class TReal>
>>> void free_matrix(TReal **m) {
>>>  free (m[0]);
>>>  delete[] m;
>>> }
>>>
>>>
>>> int main(int argc, char ** argv)
>>> {
>>>  // Matrix
>>>  // 1  3  2
>>>  // 1 -2  3
>>>  // 1  2 -3
>>>
>>>  double **a = create_matrix<double>(3,3);
>>>  a[0][0] = 1; a[0][1] = 3;  a[0][2] = 2;
>>>  a[1][0] = 1; a[1][1] = -2; a[1][2] = 3;
>>>  a[2][0] = 1; a[2][1] = 2;  a[2][2] = -3;
>>>
>>>  double **w = create_matrix<double>(3,3);
>>>  double v[3];
>>>
>>>  printf("Jacobi\n");
>>>  vtkMath::Jacobi(a,v,w);
>>>  for(int i=0; i<3; i++)
>>>    printf("%g \t%g %g %g\n",v[i],w[i][0],w[i][1],w[i][2]);
>>>
>>>
>>>  a[0][0] = 1; a[0][1] = 3;  a[0][2] = 2;
>>>  a[1][0] = 1; a[1][1] = -2; a[1][2] = 3;
>>>  a[2][0] = 1; a[2][1] = 2;  a[2][2] = -3;
>>>
>>>  printf("\nJacobiN\n");
>>>  vtkMath::JacobiN(a,3,v,w);
>>>  for(int i=0; i<3; i++)
>>>    printf("%g \t%g %g %g\n",v[i],w[i][0],w[i][1],w[i][2]);
>>>
>>> /*
>>> Output from this program:
>>> Jacobi
>>> 4.369           0.729826 -0.677467 0.0916155
>>> -2.78956        0.540772 0.490118 -0.68363
>>> -5.57944        0.418234 0.548474 0.724056
>>>
>>> JacobiN
>>> 4.369           0.729826 -0.677467 0.0916155
>>> -2.78956        0.540772 0.490118 -0.68363
>>> -5.57944        0.418234 0.548474 0.724056
>>> */
>>>
>>>
>>> /*
>>> % Matlab result:
>>>
>>> a = [[1,3,2];[1,-2,3];[1,2,-3]];
>>> [v,w] = eig(a)
>>>
>>> v =
>>>
>>>    0.8924    0.7816    0.1413
>>>    0.3570   -0.5541   -0.7279
>>>    0.2761   -0.2864    0.6709
>>>
>>>
>>> w =
>>>
>>>    2.8190         0         0
>>>         0   -1.8598         0
>>>         0         0   -4.9593
>>> */
>>>  return 0;
>>> }
>> Thanks for taking a look Bryn. JacobiN seems like it is just a
>> particular iteration of Jacobi()? Just curious why you would choose to
>> compare Jacobi to JacobiN?
>>
>> This seems like a mildly big problem... as any algorithms that are
>> using eigenvalues/vectors could be producing bogus results silently.
>> Is anyone familiar with the internals of the Jacobi functions?
>>
>> 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
>>
> _______________________________________________
> 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
-------------------------------------------------