Efficient dominant eigenvalue algorithm or the magic of Matlab

Question

Given a symmetric positive semidefinite matrix $A$, find an eigenpair $(\lambda,v)$ such that $A v = \lambda v$ with $\lambda$ maximal in the shortest amount of time.

My approaches:

1. power iteration

pseudocode:

x:=random vector of length numrows(A) x=x/||x|| 
while l changes (e.g. abs(l-l_last)<tolerance)
    x=Av
    l=||x||
    x=x/||x||

very slow but it works. I generated 100 times a 100x100 random A (sym.,pos.def.) and i get an average runtime of 0.015s in matlab. (i also implemented this in c++ with little improvements)

2. inverse iteration:

pseudocode:

x:=random vector of length numrows(A)
x=x/||x||
l=2*max(A(i,j)) % start with a value above the maximal lambda
while l changes
    x=(A-l*I)^-1 * v
    x=x/||x||
    l=x'*A*x %rayleigh quotient

i get an average runtime of 0.0136s. (i compute (A-l*I)^-1 with lu decomposition)

Now the mindtwisting part:

Matlab has a build in function called "eig" that computes all(!) eigenvalues about 10 times faster (average ~3e-4s with the above A)... To top it all off: I did compare the eigenvalues of 100 reruns (for the same A) and the matlab solution does change exactly 0. So they use a symbolic computation???... I know matlab is very secret about everything they do, but is there maybe a pseudo code or even a hint how i can do better? Seems like the power iteration/inverse power iteration/rayleigh... is very out-of-date?

(btw: matlab does not use svd=single value decomposition, since svd takes much longer than eig.)

Answer 1

For anyone that has the same problem, here is a solution/template using lapack (thanks MrYouMath) in c++:

//compile: g++ test_lapack.cpp -o test_lapack -llapack
//run: ./test_lapack matrix.txt

#include <iostream>
#include <fstream>
#include <vector>


using namespace std;


extern "C" {

//dspevx 
//-      double
// --    symmetric matrix
//   --- eigenvalue expert (-> first k eigenvalues...)
extern int dspevx_(char*, //jobz 'V':  Compute eigenvalues and eigenvectors.
                  char*, //RANGE 'I': the IL-th through IU-th eigenvalues will be found.
                  char*, //?? UPLO 'U'/'L'
                  int*, //The order of the matrix A.  N >= 0.
                  double*, // AP 
                  double*, //if range V
                  double*, //if range V
                  int*, //IL smallest eigenvalue to be returned.1 <= IL <= IU <= N
                  int*, //IU largest eigenvalue to be returned.1 <= IL <= IU <= N
                  double*, // ABSTOL The absolute error tolerance for the eigenvalues.
                  int*, //M The total number of eigenvalues found. 'I', M = IU-IL+1
                  double*, //W DOUBLE PRECISION array, dimension (N)
                  double*, //Z DOUBLE PRECISION array, dimension (LDZ, max(1,M))
                  int*, //LDZ
                  double*, //WORK
                  int*, //IWORK
                  int*, //IFAIL If INFO > 0, then IFAIL contains the indices of the eigenvectors that failed to converge.If JOBZ = 'N', then IFAIL is not referenced.
                  int*); //INFO is
}

int main(int argc, char** argv){

  // check for an argument
  if (argc<2){
    cout << "Usage: " << argv[0] << " " << " filename" << endl;
    return -1;
  }

  int n,m;
  double *data;

  // read in a text file that contains a real matrix stored in column major format
  // but read it into row major format
  ifstream fin(argv[1]);
  if (!fin.is_open()){
    cout << "Failed to open " << argv[1] << endl;
    return -1;
  }
  fin >> n >> m;  // n is the number of rows, m the number of columns
  data = new double[n*m];
  for (int i=0;i<n;i++){
    for (int j=0;j<m;j++){
      fin >> data[j*n+i];
    }
  }
  if (fin.fail() || fin.eof()){
    cout << "Error while reading " << argv[1] << endl;
    return -1;
  }
  fin.close();

  // check that matrix is square
  if (n != m){
    cout << "Matrix is not square" <<endl;
    return -1;
  }

  n=3;
  // allocate data
  char Nchar='N';
  char Vchar='V';
  char Ichar='I';
  char Uchar='U';
  char Achar='A';
  double *eigenvalues=new double[1];
  double *vl,*vr;
  int IL=n-1;
  int UL=n;
  int M=UL-IL+1;
  int lwork=6*n;
  double *work=new double[8*n];
  int info;
  int *iwork=new int[5*n];
  int *ifail=new int[n];
  int ldz=n;
  vector<double> eigenvectors(ldz*n);
  eigenvectors[0]=1;


  double abstol = 0; //= default value
  vector<double> laplacian(n*n);

laplacian[0] = 0;           
laplacian[1] = 1;            
laplacian[2] = 1;            
laplacian[3] = 1;          
laplacian[4] = 1;            
laplacian[5] = 1;          
laplacian[6] = 1; 

  dspevx_(&Vchar, //jobz 'V':  Compute eigenvalues and eigenvectors.
       &Ichar, //RANGE 'I': the IL-th through IU-th eigenvalues will be found.
       &Uchar, //?? UPLO 'U'/'L'
       &n, //The order of the matrix A.  N >= 0.
       &laplacian[0], // AP 
       vl, //if range V
       vr, //if range V
       &IL, //IL smallest eigenvalue to be returned.1 <= IL <= IU <= N
       &UL, //IU largest eigenvalue to be returned.1 <= IL <= IU <= N
       &abstol, // ABSTOL The absolute error tolerance for the eigenvalues.
       &M, //M The total number of eigenvalues found. 'I', M = IU-IL+1
       eigenvalues, //W DOUBLE PRECISION array, dimension (N)
       &eigenvectors[0], //Z DOUBLE PRECISION array, dimension (LDZ, max(1,M))
       &ldz, //LDZ
       work, //WORK
       iwork, //IWORK
       ifail, //IFAIL
       &info); //INFO is


  // check for errors
  if (info!=0){
    cout << "Error: dgeev returned error code " << info << endl;
    return -1;
  }

  // output eigenvalues to stdout
  cout << "--- Eigenvalues ---" << endl;
  for (int i=0;i<M;i++){
    cout << "( " << eigenvalues[i] << " )\n";
    for (int j=0;j<n;j++){
      cout << eigenvectors[i*n+j] << " ";
    }
    cerr << endl;
  }
  cout << endl;

  // deallocate
  delete [] data;
  delete [] eigenvalues;
  delete [] work;
  delete [] iwork;
  delete [] ifail;

  return 0;
}