Hi All,

I'm experiencing a very interesting issue I try to solve the teh generalized eigenvalue problem where all inner products of matrices are positive and matrices are symmetrical. Compared with the other libraries the results are drastically different. Having such a different results made me think that I have to apply additional multiplication, division or any another operation to get the correct results. I use the following code to get the eigenvalues.

Code: Select all

Eigen::GeneralizedSelfAdjointEigenSolver<MatrixXf> es;
es.compute(A,B);

float num;
myfile.open("result.txt");
stringstream ss;

myfile << "\n Eigen alpha values\n";
myfile << es.eigenvalues() ;





What's the output of (A*es.eigenvectors() - B*es.eigenvectors()*es.eigenvalues().asDiagonal()).norm() / (A*es.eigenvectors()).norm() ??

If that's close to zero, then there is an issue and we would need the values of A and B to debug it. Also, note that the matrix B must be positive definite, not only symmetric.

Hi ggael,

Is it possible that there is mis-typing in your expression, I'm not very much familiar with Eigen syntax for the time being, but I already tried every possible combinaation of placing the paranthesis but never worked. AS far as I can deduct from the expression it is kind of back-substitution in the original eigen equation, to testify that the righ hand and left hand side are equal to each other ?

Regards,

works for me with Eigen 3.2:

Code: Select all

#include <Eigen/Dense>
#include <iostream>
using namespace Eigen;
int main()
{
  MatrixXf A(10,10), B(10,10);
  A.setRandom();
  B.setRandom();
  A = A.adjoint()*A;
  B = B.adjoint()*B;
  Eigen::GeneralizedSelfAdjointEigenSolver<MatrixXf> es;
  es.compute(A,B);
  std::cout << (A*es.eigenvectors() - B*es.eigenvectors()*es.eigenvalues().asDiagonal()).norm() / (A*es.eigenvectors()).norm() << std::endl;
}


btw, most numerical tools work on double by default, so it might be fair to compare using double as well.

Sorry my bad, As a result yes it produces a value close to zero 5.25705e-6.


Below is the main appilcation body :

The matrix notation given in IN.txt file is actually stiffness matrix which is widely used in engineering problems, I only store
the upper triangle of matrices and populate the rest in code manually. The other file mass.txt is mass matrix of simplified form
where there is only products on the diagonal of it.

Files needed during the execution are

https://www.dropbox.com/s/uv0xodjo0zopqzx/IN.txt
https://www.dropbox.com/s/d4d99akbehbv98d/mass.txt

Code: Select all

#include <iostream>
#include <Eigen/Eigenvalues>
#include <stdio.h>
#include <fstream>
#include <vector>
#include <Eigen/Eigenvalues>

using namespace std;
using Eigen::MatrixXf;

int main()
{
Eigen::GeneralizedEigenSolver<MatrixXf> ges;
MatrixXf A = MatrixXf::Random(4,4);
MatrixXf B = MatrixXf::Random(4,4);
FILE * file_in;
int max_size, column, row;
double value;

    file_in = fopen(".........................../IN.txt","r");
    if (file_in!=NULL )
    {
        fscanf(file_in,"%d",&max_size);
        //-K= dmatrix(1,max_size,1,max_size);
        A.resize(max_size,max_size);
        while ( fscanf(file_in,"%d",&row) !=EOF )
        {
          fscanf(file_in,"%d",&column);
          fscanf(file_in,"%lf",&value);
          //K[row][column]=K[column][row]=value;
           A(row-1,column-1) = A(column-1,row-1) = value;
        }
    }
    fclose(file_in);

    file_in = fopen("............................./mass.txt","r");
    if (file_in!=NULL )
    {
        fscanf(file_in,"%d",&max_size);
        //-K= dmatrix(1,max_size,1,max_size);
        B.resize(max_size,max_size);
        while ( fscanf(file_in,"%d",&row) !=EOF )
        {
          fscanf(file_in,"%d",&column);
          fscanf(file_in,"%lf",&value);
          //K[row][column]=K[column][row]=value;
           B(row-1,column-1) = B(column-1,row-1) = value;
        }
    }
    fclose(file_in);

//ges.compute(A, B);
Eigen::GeneralizedSelfAdjointEigenSolver<MatrixXf> es;
es.compute(A,B);

fstream myfile;

std::vector< float > vec;
float num;
myfile.open("result.txt");
stringstream ss;

myfile << "\n Eigen alpha values\n";
myfile << es.eigenvalues();
myfile >> num;
myfile.close();
return 0;

}



I don't see anything wrong.

Moreover, I get the same result than with MatLab. However, in the case that not all entries are stored in the files, then make sure that you call A.setZero(); and B.setZero(); after resizing them. By default, Eigen's matrices are not initialized to zero.

If you used the main code, can you share your "result.txt" file output produced from application. I used another independent library Slepc and another source code both are using subspace algorithm to find the lowest eigenvalues and they converged nicely #2) . Although the Eigen and forementioned libraries use different algorithms aren't they supposed to produce the same converging results?

I provide the outputs from both libraries.
#1) Eigen library produced results which start with following values :
2.38269
4.8292
6.82543
8.96906
14.5238
16.0879
16.2533
20.5156
20.5895
27.6132
71.2584
157.663
157.777
164.934
297.029
586.402
586.969
870.542
991.488
1462.9
1463.6
2456.81
2613.91
2614.97
2932.81
3217.96
3238.51
3436.77
3795.85
3814.18
3816.28
3948.37
3974.94
5628.21
6787.2
7017.32
7090.46
7210.12
7267.28
10928.5
16859.8
19314.4
30256.4
32426
44611.2
46697

#2) Slepc
0.168948
1.714001
2.605187
2.688302
3.523136
4.037952
6.982649
7.802687
9.924228

I got:

Code: Select all

0.166705
 1.71996
 2.61373
 2.69565
 3.21322
 4.03804
 4.79298
 6.94157
 7.80227
 7.82321
 9.92228
 10.9546
 14.5828
 33.3783
  34.682
 39.8336
 44.3472
 51.0768
  51.098
 61.2854
 61.6912
 77.3527
 92.6374
 93.1729
 95.1073
 165.708
 167.405
 190.333
 242.753
 249.122
 288.596
 501.648
 520.598
 566.125
 602.389
 627.022
 680.474
 686.067
 686.114


but do not forget to add the setZero()!

OK thanks after your guidance I found also the same result with slight difference. What could be the reason behind of those slight deviations in eigenvalues ? Do you think that this is due to the floating precision ?

Regards,