3

I am trying to determine the eigenvalues and eigenvectors of a sparse array in Eigen. Since I need to compute all the eigenvectors and eigenvalues, and I could not get this done using the unsupported ArpackSupport module working, I chose to convert the system to a dense matrix and compute the eigensystem using SelfAdjointEigenSolver (I know my matrix is real and has real eigenvalues). This works well until I have matrices of size 1024*1024 but then I start getting deviations from the expected results.

In the documentation of this module (https://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html) from what I understood it is possible to change the number of max iterations:

const int m_maxIterations static Maximum number of iterations.

The algorithm terminates if it does not converge within m_maxIterations * n iterations, where n denotes the size of the matrix. This value is currently set to 30 (copied from LAPACK).

However, I do not understand how do you implement this, using their examples:

SelfAdjointEigenSolver<Matrix4f> es;
Matrix4f X = Matrix4f::Random(4,4);
Matrix4f A = X + X.transpose();
es.compute(A);
cout << "The eigenvalues of A are: " <<   es.eigenvalues().transpose() << endl;
es.compute(A + Matrix4f::Identity(4,4)); // re-use es to compute eigenvalues of A+I
cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl

How would you modify it in order to change the maximum number of iterations?

Additionally, will this solve my problem or should I try to find an alternative function or algorithm to solve the eigensystem?

My thanks in advance.

jcarvalho
  • 115
  • 2
  • 11

3 Answers3

2

Increasing the number of iterations is unlikely to help. On the other hand, moving from float to double will help a lot!

If that does not help, please, be more specific on "deviations from the expected results".

ggael
  • 28,425
  • 2
  • 65
  • 71
  • Basically after I obtain the eigenvectors I need to compute an expectation value, that is calculate something like: [a_1.....a_m]*M*[a_1.....a_m]^T, a_m being the coefficients of the eigenvectors. I expect to obtain integer values. For a 256*256 matrix for all the eigenvectors I obtain something like 4.999999... but for a 1024*1024 matrix almost all of them work, but I also obtain something like 0.2 and 0.4. Physically speaking it would be impossible not to obtain non integer values, and 0.2 and 0.4 a a bit too far off. – jcarvalho Oct 03 '16 at 23:45
  • And have you tried with double precision? (i.e., with `MatrixXd`). You can also check the accuracy by testing that the input matrix is properly reconstructed: `norm(A-(eig.eigenvectors()*eig.eigenvalues().asDiagonal()).eval()*eig.eigenvectors().transpose()) / norm(A)`. – ggael Oct 04 '16 at 15:29
  • Sorry I forgot to mention it, all my matrices are implemented using double precision, I wonder if trying quadruple precision would be worth a shot. – jcarvalho Oct 04 '16 at 20:16
  • OK, so what about the precision test I suggested above? If it returns something like 1e-15, there not much we can do on the solver side. Also, it could be interesting to check the eigenvalues for their range and multiplicity. – ggael Oct 04 '16 at 21:22
1

m_maxIterations is a static const int variable, and as such it can be considered an intrinsic property of the type. Changing such a type property usually would be done via a specific template parameter. In this case, however, it is set to the constant number 30, so it's not possible.

Therefore, you're only choice is to change the value in the header file and recompile your program.

However, before doing that, I would try the Singular Value Decomposition. According to the homepage, its accuracy is "Excellent-Proven". Moreover, it can overcome problems due to numerically not completely symmetric matrices.

davidhigh
  • 14,652
  • 2
  • 44
  • 75
  • 1
    Do you mean Jacobi SVD? If so could you elaborate on how do I get Eigevalues and Eigenvectors from this decomposition? – jcarvalho Oct 03 '16 at 23:36
1

I solved the problem by writing the Jacobi algorithm adapted from the Book Numerical Recipes:

void ROTATy(MatrixXd &a, int i, int j, int k, int l, double s, double tau)
{
double g,h;
g=a(i,j);
h=a(k,l);
a(i,j)=g-s*(h+g*tau);
a(k,l)=h+s*(g-h*tau);

}

void jacoby(int n, MatrixXd &a, MatrixXd &v, VectorXd &d )
{
int j,iq,ip,i;
double tresh,theta,tau,t,sm,s,h,g,c;


VectorXd b(n);
VectorXd z(n);

v.setIdentity();    
z.setZero();


for (ip=0;ip<n;ip++)
{   
    d(ip)=a(ip,ip);
    b(ip)=d(ip);
}


for (i=0;i<50;i++) 
{
    sm=0.0;
    for (ip=0;ip<n-1;ip++) 
    {

        for (iq=ip+1;iq<n;iq++)
            sm += fabs(a(ip,iq));
    }

    if (sm == 0.0) {
        break;
    }

    if (i < 3)
    tresh=0.2*sm/(n*n); 
    else
    tresh=0.0;  


    for (ip=0;ip<n-1;ip++)
    {
        for (iq=ip+1;iq<n;iq++)
        {
            g=100.0*fabs(a(ip,iq));
            if (i > 3 && (fabs(d(ip))+g) == fabs(d[ip]) && (fabs(d[iq])+g) == fabs(d[iq]))
            a(ip,iq)=0.0;
            else if (fabs(a(ip,iq)) > tresh)
            {
                h=d(iq)-d(ip);
                if ((fabs(h)+g) == fabs(h))
                {
                    t=(a(ip,iq))/h;
                }   
                else 
                {
                    theta=0.5*h/(a(ip,iq));
                    t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
                    if (theta < 0.0)
                    {
                        t = -t;
                    }
                    c=1.0/sqrt(1+t*t);
                    s=t*c;
                    tau=s/(1.0+c);
                    h=t*a(ip,iq);


                    z(ip)=z(ip)-h;
                    z(iq)=z(iq)+h;
                    d(ip)=d(ip)- h;
                    d(iq)=d(iq) + h;
                    a(ip,iq)=0.0;


                    for (j=0;j<ip;j++)
                        ROTATy(a,j,ip,j,iq,s,tau);
                    for (j=ip+1;j<iq;j++)
                        ROTATy(a,ip,j,j,iq,s,tau);
                    for (j=iq+1;j<n;j++)
                        ROTATy(a,ip,j,iq,j,s,tau);
                    for (j=0;j<n;j++)
                        ROTATy(v,j,ip,j,iq,s,tau);


                }
            } 
        }
    }


}

}

the function jacoby receives the size of of the square matrix n, the matrix we want to calculate the we want to solve (a) and a matrix that will receive the eigenvectors in each column and a vector that is going to receive the eigenvalues. It is a bit slower so I tried to parallelize it with OpenMp (see: Parallelization of Jacobi algorithm using eigen c++ using openmp) but for 4096x4096 sized matrices what I did not mean an improvement in computation time, unfortunately.

Community
  • 1
  • 1
jcarvalho
  • 115
  • 2
  • 11