Monday, October 1, 2012

Jacobi eigenvalue method implementation in C++


The Jacobi eigenvalue problem is an algorithm to compute the eigenvalues of a matrix by canceling out the off-diagonal elements by multiplying the matrices with rotation matrices.

I started with the code in Wikipedia, which is essentially wrong.

Here I offer an implementation that uses Boost uBLAS matrices. It uses a trick to avoid computing sines and cosines. It needs to be improved to avoid multiplication of full matrices when doing the rotations.

This code is copylefted, use it as you wish.

#include "stdafx.h"

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace boost::numeric::ublas;
using namespace boost;

#define EPS 10e-08


void abs(matrix<double> &M)
{
 int n = M.size1();

 for (int k=0; k<n; k++)
 {
  for (int i=0; i<n;i++)
  {
   M(k,i)=abs(M(k,i));
  }
 }
}


void findMax(matrix<double> &M, int &row, int &col)
{
 double m = M(0,0);
 row = col = 0;
 int n = M.size1();

 for (int k=0; k<n; k++)
 {
  for (int l=0; l<n;l++)
  {
   if(M(k,l) > m)
   {
    row = k;
    col = l;
    m = M(k,l);
   }
  }
 }
}


void zeroRowCol(matrix<double> &S, matrix<double> &U, int row, int col)
{
 double t, c, s, theta;
 int n = S.size1();
 //matrix<double> S(n,n);
 //S = S0;

 if (row == col) return;

 theta=(S(row,row)-S(col,col))/(2*S(row,col));
 if (theta < EPS)
 {
  t = 1/(abs(theta)+sqrt(theta*theta+1));
 }
 else {
  t = 1/abs(2*theta);
 }

 if (theta<0) t = -t;

 c = 1/sqrt(t*t + 1);
 s = c*t;
 matrix<double> R(n,n);
 R = identity_matrix<double>(n,n);
 R(row,row) = R(col,col) = c; R(row,col) = s; R(col,row) = -s;
 S = prod(S, trans(R));
 S = prod(R, S);
 U = prod(U, R);
}


int jacobi2(matrix<double> &S, vector<double> &e, matrix<double>  &U)
{
 int col, row;
 bool iterating = true;
 int n = S.size1();
 if (S.size2() != n)
 {
  return -1;
 }
 matrix<double> M(n,n);

 U = identity_matrix<double>(n,n);
 
 while(iterating)
 {
  M = S;
  abs(M);

  for (int k=0; k<n; k++)
  {
   M(k,k)=0;
  }

  findMax(M, row, col);
  if (row == col)
  {
   for (int i=0; i<n; i++) e(i) = S(i,i);
   return 0;
  }
  double Smax = S(row,col);
  zeroRowCol (S, U, row, col);
  if (Smax < EPS * norm_frobenius(S)) iterating = false;
 }

 for (int i=0; i<n; i++) e(i) = S(i,i);

 return 0;  
} 


int _tmain(int argc, _TCHAR* argv[])
{
 int n = 3;

 mmatrix<double> M(n,n);
 M = boost::numeric::ublas::identity_matrix<double>(3,3);
 M(0,1)=.5;
 M(1,0)=.5;

 std::cout << "The matrix is:\n" << M << std::endl;

 matrix<double> M(n,n);
 matrix<double> U(n,n);
 vector<double> e(n); 

 jacobi2(M, e, U);

 std::cout << M << std::endl;
 std::cout << e << std::endl;
 std::cout << U << std::endl;
}

3 comments:

  1. Thank you, this looks very valuable for my scientific project. Basically, my task is to investigate eigenvalues in a sequence of matrices, which are symmetrical and consist of three unit elements in each row and column, and zero elements elsewhere. The maximal size I need to tackle in 1024*1024. I'm not an expert in programming, despite my attempts to improve my skills by doing the project, and it might take me some time to comprehend the algorithm. Here are some basic questions I have regarding the code:
    1) how fast is it supposed to process a 1024*1024 matrix as described above?
    2) does it require any libraries not included into the standard C++ package?
    3) where are the descriptions of matrix and similar structures? is there another part of the code that I need to make it work?

    I'm looking forward to hearing from you!

    Max.

    ReplyDelete
  2. I think your code may have some issues. For instance, it does not correctly solve the example on the wikipedia page (http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm). The element in the example matrix with greatest magnitude is actually negative. When you compare Smax to EPS it will always be false on the first step. Since computation stops after one iteration, the answer is obviously incorrect.

    Is there a simple fix to this issue for your implementation?

    ReplyDelete
  3. Hi there! I am sorry for the delay but I have been essentially away for a couple of months.

    I admit I did not test it too thoroughly but I believe it did OK for some examples. I have some other priorities right now but I will have a look back to it.

    ReplyDelete