Hello,

I'm searching eigen values (lambda vector) and eigen vectors (x) for the A.x = lambda.B.x system.
A and B are 6x6 squared symetric non-inversible matrixs (see below values).

Code: Select all

A = [2.926E+10 1.854E+10 1.220E+10 1.461E+08 92694180.0 739772.0
1.854E+10 1.220E+10 8.323E+09 92694180.0 61350634.0 469726.0
1.220E+10 8.323E+09 5.884E+09 61350634.0 42305670.0 312844.0
1.461E+08 92694180.0 61350634.0 739772.0 469726.0 3810.0
92694180.0 61350634.0 42305670.0 469726.0 312844.0 2418.0
739772.0 469726.0 312844.0 3810.0 2418.0 20.0];

B = [0 0 2 0 0 0
0 -1 0 0 0 0
2 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0 ]

With Eigen 3, I've try the "GeneralizedSelfAdjointEigenSolver<MatrixXd>" class and the below program :

Code: Select all

   // Lecture du système à résoudre A.x = lambda B.x
   MatrixXd matA = CLinAlg::getMatrice( "matA2.txt" );
   cout << "Matrix A : " << endl << matA << endl;

   MatrixXd matB = CLinAlg::getMatrice( "matB2.txt" );
   cout << "Matrix B : " << endl << matB << endl;;

   // Résolution de A.x = lambda B.x avec Eigen
   cout << "Resolution A.x = lambda B.x" << endl;
   MatrixXd vectValeurs_propres;
   MatrixXd matVecteurs_propres;
   ComputationInfo eInfo;
   vector<string> vInfos_retour;
   vInfos_retour.push_back( "Succes ! (= 0)" );
   vInfos_retour.push_back( "The provided data did not satisfy the prerequisites. (= 1)" );
   vInfos_retour.push_back( "Iterative procedure did not converge. NoConvergence (= 2) ");
   vInfos_retour.push_back( "The inputs are invalid, or the algorithm has been improperly called. InvalidInput (= 3)");

   GeneralizedSelfAdjointEigenSolver<MatrixXd> es;
   es.compute( matA, matB, ComputeEigenvectors | Ax_lBx );//, Ax_lBx EigenvaluesOnly );
   eInfo = es.info();
   cout << "Computation result : " << vInfos_retour[ eInfo ] << endl;
   vectValeurs_propres = es.eigenvalues();
   cout << "Eigen values = " << endl << vectValeurs_propres << endl;;

   matVecteurs_propres = es.eigenvectors();
   cout << "Eigen vectors = " << endl << matVecteurs_propres << endl;

After the execution I obtain the following results :

Code: Select all

Matrix A :
  2.926e+010   1.854e+010    1.22e+010   1.461e+008 9.26942e+007       739772
  1.854e+010    1.22e+010   8.323e+009 9.26942e+007 6.13506e+007       469726
   1.22e+010   8.323e+009   5.884e+009 6.13506e+007 4.23057e+007       312844
  1.461e+008 9.26942e+007 6.13506e+007       739772       469726         3810
9.26942e+007 6.13506e+007 4.23057e+007       469726       312844         2418
      739772       469726       312844         3810         2418           20
Matrix B :
 0  0  2  0  0  0
 0 -1  0  0  0  0
 2  0  0  0  0  0
 0  0  0  0  0  0
 0  0  0  0  0  0
 0  0  0  0  0  0
Resolution A.x = lambda B.x
Computation result : Iterative procedure did not converge. NoConvergence (= 2)
Eigen values =
nan
nan
nan
nan
nan
nan
Eigen vectors =
nan nan nan nan nan nan
nan nan nan nan nan nan
nan nan nan nan nan nan
nan nan nan nan nan nan
nan nan nan nan nan nan
nan nan nan nan nan nan


With Eigen, I don't have a result but with the "spec" function of Scilab, I find :

Code: Select all

Computation with the spec( A, B) function of Scilab :
- lambdas (eigen values) ( we find the same eigne values with the "schur(A, B)" function of Scilab)
 
 10^13 *
   - 0.0000006 
  - 0.0000001 
    8.334D-09 
    7.793D+10 
  - 1793612.7 
  - 4.278D+08 

- Eigen vectors (1 vector by column )
     0.0000701    0.0000140  - 0.0000833  - 1.674D-20  - 2.472D-20  - 1.622D-20 
  - 0.0001784    0.0000567    0.0001497  - 7.836D-20  - 4.949D-20  - 1.526D-20 
    0.0000269  - 0.0000775  - 0.0001585  - 3.269D-20  - 9.897D-21    8.395D-20 
  - 0.0041215  - 0.0138240    0.0093706    0.0112268    0.0097287    0.0010433 
    0.0244922    0.0082455    0.0104291  - 0.0001927  - 0.0001927    0.0196257 
  - 1.           1.         - 1.         - 1.         - 1.         - 1.   

I'm interrested by the smallest positive value of lambda.
The last 3 values are not exact because it must be divided by zero (I've remplaced zero by %eps), but its is not my soluce.

It was just for information for the Eigen developpers... Thank's for this library !

Do you know wich algorithms are used in Eigen and Scilab for this problem ?

As explain in the doc, to use GeneralizedSelfAdjointEigenSolver, B must be selfadjoint positive definite because its Cholesky factorization LL^T=B has to be formed to transform the initial problem to a standard selfadjoint eigen-problem. If you only need the eigenvalues, you can use the RealQZ factorization of Eigen.

Hello,

Thank's for this answer, ggael !

ggael wrote:As explain in the doc, to use GeneralizedSelfAdjointEigenSolver, B must be selfadjoint positive definite

Yes, it's true, I have not enough read the Eigen manual...
My matrix B is not positive definite.

ggael wrote:If you only need the eigenvalues, you can use the RealQZ factorization of Eigen.

I need eigenvalues and eigenvectors, but eigenvalues is a good start, thank's !
I have not find the "RealQZ.h" file in my Eigen v3.1.1 source files...
So I have copy it from the Eigen website to my hard disk.
Here is my source file test :

Code: Select all

#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h>
#include <Eigen/src/Eigenvalues/RealQZ.h>
using namespace Eigen;// v3.0

int main(int argc, char *argv[]) {
   ...
   RealQZ< MatrixXd > qz( 6 );
   qz.setMaxIterations( 4000 );// default is 1000
   qz.compute( matA, matB ); // A = Q S Z, B = Q T Z
   // print original matrices and result of decomposition
   cout << "A:\n" << matA << "\n\n" << "B:\n" << matB << "\n";
   MatrixXd matQ = qz.matrixQ();
   MatrixXd matS = qz.matrixS();
   MatrixXd matT = qz.matrixT();
   MatrixXd matZ = qz.matrixZ();
   cout << "S:" << endl << matS << endl;
   cout << "T:" << endl << matT << endl;
   cout << "Q:" << endl << matQ << endl;
   cout << "Z:" << endl << matZ << endl;

   MatrixXd matVerif;
   cout << "Verification Q.S.Z == A ?" << endl;
   matVerif = matQ * matS * matZ;
   cout << matVerif << endl;

   cout << "Verification Q.T.Z == B ?" << endl;
   matVerif = matQ * matT * matZ;
   cout << matVerif << endl;

   // Eigen values computing
   vectValeurs_propres.resize( 6, 1 );
   for ( int i = 0; i < 6; ++i ) {
      double nNum = matS( i, i );
      double nDen = matT( i, i );
      if ( nDen != 0.0 ) {
         vectValeurs_propres( i ) = nNum / nDen;
      }
      else {
         vectValeurs_propres( i ) = -9999;
      }
   }
   cout << "Eigen values = " << endl << vectValeurs_propres << endl;
   ...
}// main


And here is the execution results :

Code: Select all

Resolution A.x = lambda B.x
A:
   2.93e+010    1.85e+010    1.22e+010    1.46e+008    9.27e+007     7.4e+005
   1.85e+010    1.22e+010    8.32e+009    9.27e+007    6.14e+007     4.7e+005
   1.22e+010    8.32e+009    5.88e+009    6.14e+007    4.23e+007    3.13e+005
   1.46e+008    9.27e+007    6.14e+007     7.4e+005     4.7e+005    3.81e+003
   9.27e+007    6.14e+007    4.23e+007     4.7e+005    3.13e+005    2.42e+003
    7.4e+005     4.7e+005    3.13e+005    3.81e+003    2.42e+003           20

B:
 0  0  2  0  0  0
 0 -1  0  0  0  0
 2  0  0  0  0  0
 0  0  0  0  0  0
 0  0  0  0  0  0
 0  0  0  0  0  0
S:
          432     3.99e+004     9.75e+004    -3.45e+010    -8.76e+009     -2.5e+010
          395     4.73e+003     1.48e+004    -4.62e+009    -1.17e+009    -4.23e+009
            0             0    -4.44e+005     1.19e+010     3.02e+009      9.2e+009
            0             0             0     1.73e+008     4.38e+007     1.26e+008
            0             0             0             0     5.13e+003     -3.5e+006
            0             0             0             0             0     7.13e+005
T:
     0.000317        0.0121        -0.104        -0.439          1.31          -1.4
            0       -0.0023       -0.0228         -1.95        -0.301         0.308
           -0             0        0.0354        0.0149        -0.518        -0.934
            0             0             0             0     -2.6e-018    -1.11e-017
            0             0             0             0             0     1.62e-019
            0             0             0             0             0     0
Q:
       -0.953         0.214        -0.217    -6.64e-018     1.93e-019    -1.14e-020
        -0.21        0.0553         0.976    -9.94e-018      1.6e-019     6.53e-020
       -0.221        -0.975       0.00784             0             0  0
   -4.94e-018     1.17e-018     6.83e-018         0.844         0.526        -0.102
   -4.83e-018     1.08e-018    -1.44e-018         0.536        -0.827         0.169
   -3.49e-020     9.28e-021     1.71e-019       0.00428        -0.197         -0.98
Z:
   -3.5e-005    6.64e-005    -0.000151      -0.0887       -0.151       -0.985
   -0.000218      0.00267     -0.00602        0.168        0.972       -0.164
      0.0227       -0.055       0.0431       -0.979         0.18       0.0607
           1      0.00125     -0.00098       0.0223     -0.00388     -0.00145
           0        0.797         -0.6      -0.0696      0.00701      0.00534
           0        0.602        0.799      0.00182      0.00286    -0.000684

Verification Q.S.Z == A ?
    2.93e+010     1.85e+010     1.22e+010     1.72e+008    -1.46e+006     -1.7e+007
    1.85e+010     1.22e+010     8.32e+009     1.11e+008     1.18e+006    -1.13e+007
    1.22e+010     8.32e+009     5.88e+009     7.41e+007     2.19e+006    -7.78e+006
    1.46e+008     9.27e+007     6.14e+007     8.72e+005    -7.06e+003    -8.58e+004
    9.27e+007     6.14e+007     4.23e+007     5.61e+005     7.59e+003    -5.74e+004
     7.4e+005      4.7e+005     3.13e+005     4.49e+003         -37.3          -442

Verification Q.T.Z == B ?
   -5.84e-016     8.03e-017             2    -7.09e-017     1.75e-017     4.25e-018
   -5.49e-017            -1     1.14e-016    -2.29e-017     6.88e-018     1.42e-018
            2    -2.27e-016    -1.13e-016    -2.08e-017     7.29e-018     1.91e-018
   -2.31e-033    -1.51e-017     8.72e-019     1.36e-019    -4.19e-020    -5.39e-021
   -2.06e-033    -4.43e-018     6.28e-018     8.58e-020    -2.71e-020    -3.29e-021
   -2.19e-035    -2.32e-019    -6.05e-020     6.29e-022    -3.05e-022     -5.1e-024

Eigen values =
1.36e+006
-2.06e+006
-1.25e+007
-1e+004
-1e+004
-1e+004


It's look work fine, thank's.
But if you look at "Verification Q.S.Z == A ?" line, the A( 6, 6 ) value (=20) and the (Q.S.Z)( 6, 6 ) value (=-442), there is a large difference...

Like for Scilab Schur decomposition ([As,Es,Q,Z] = schur(A,E) returns in addition two unitary matrices Q and Z such that As=Q'*A*Z and Es=Q'*E*Z. ), I try to compute the eigenvalues with a As(i,i)/Es(i,i) division ( diag( As )./ diag( Es ) ).
The Scilab eigenvalues are very different... (I am only interrested by the smaller positive one ( 83340 ), and it's associate eigenvector.

Code: Select all

Scilab Eigen values :
 10^13 *
   - 0.0000006 
  - 0.0000001 
    8.334D-09 
    7.793D+10 
  - 1793612.7 
  - 4.278D+08


I've look at the RealSchur Eigen function, but it seems doesn't works with two matrix for argument like Scilab [As,Bs,Q,Z] = schur(A,B)...

Do you have an idea to improve my Eigen-C++ numerical results ?