[Insight-users] calculating eigen values and vector of matrix

[Miller, James V (Research)]

Leila, 
 
Here is an example of using the vnl routines.  The code a eigenvalue decomposition of a matrix M and extracts the eigenvector corresponding to the smallest eigenvalue.  I use the vnl_real_eigensystem to compute the eigensystem for the real matrix M.  The results (eigenvalues and eigenvectors) are always
"returned" as complex numbers matrices of complex numbers (eigenvalues stored in a vnl_diag_matrix).
 
To start this process with an itk::Matrix, you should be able to just the result of GetVnlMatrix() to the eigensolver.
 
#include <vnl/vnl_matrix.h>
#include <vnl/vnl_math.h>
#include <vnl/vnl_vector.h>
#include <vnl/algo/vnl_real_eigensystem.h>

{
  vnl_matrix<double> M(3, 3);
  vnl_vector<double> params;
  double lambda;

  M = .......;  // some matrix M
 
  // construct the eigenvector problem
  vnl_real_eigensystem eig(M);

  // extract eigenvector corresponding to the smallest positive real eigenvalue
  vnl_diag_matrix<vcl_complex<double> >::iterator it;
  vnl_diag_matrix<vcl_complex<double> >::iterator small_it = eig.D.end();
  double small = vcl_numeric_limits<double>::max();
  const double eigenvalue_zero_tol = 1.0e-11;
  
  for (it = eig.D.begin(); it != eig.D.end(); ++it)
    {
    if ((*it).real() > eigenvalue_zero_tol
        && vcl_fabs((*it).imag()) < eigenvalue_zero_tol
        && (*it).real() < small)
      {
      small_it = it;
      }
    }

  params.set_size( eig.V.rows() ); // destroys previous values
  for (i=0; i < eig.V.rows(); i++)
    {
    params(i) = eig.V(i, small_it - eig.D.begin()).real();
    }
  lambda = (*small_it).real();
}

	-----Original Message----- 
	From: Leila baghdadi [mailto:baghdadi at sickkids.ca] 
	Sent: Wed 4/6/2005 4:44 PM 
	To: Miller, James V (Research) 
	Cc: 
	Subject: RE: [Insight-users] calculating eigen values and vector of matrix
	
	

	Hello Jim,
	
	I see what you mean but I am having difficulty using the vnl methods
	from my class,
	
	can you please give me a simple example,
	
	
	Thanks
	
	Leila
	
	On Wed, 2005-04-06 at 16:03, Miller, James V (Research) wrote:
	> Leila,
	> 
	> You'll have to call GetVnlMatrix() on your itk::Matrix and pass that to a vnl routine. There are several vnl eigensystem methods to choose from depending on the structure of your matrix (real, symmetric, ...)
	> 
	> Jim
	>
	>       -----Original Message-----
	>       From: insight-users-bounces at itk.org on behalf of Leila baghdadi
	>       Sent: Wed 4/6/2005 3:52 PM
	>       To: Insight-users at itk.org
	>       Cc:
	>       Subject: [Insight-users] calculating eigen values and vector of matrix
	>      
	>      
	>
	>       Hi everyone,
	>      
	>       I know that the vnl utilities should have methods for the above
	>      
	>       but can you just use itkMatrix,
	>      
	>       I do not see any public methods for calculating eigen values or vector
	>      
	>      
	>       any hints
	>      
	>       thanks
	>      
	>       Leila
	>      
	>       _______________________________________________
	>       Insight-users mailing list
	>       Insight-users at itk.org
	>       http://www.itk.org/mailman/listinfo/insight-users
	>      
	>