0
Answered

Real generalized eigenvalue problem

Stefan Güttel 4 years ago updated by Pavel Holoborodko 4 years ago 2

I have a generalized eigenvalue problem Ax = lambda*Bx, where A and B are real multiprecision matrices. I used eig(A,B) and expected that non-real eigenvalues came out in perfectly complex conjugate pairs (because I want to group them in pairs), but this is not the case. I'm not sure this should be considered a bug, or rather a question about whether there is a way to enforce that EIG returns perfectly conjugate pairs when A and B are real mp matrices?


Here is a simple example (tested in MATLAB 2015A and MCT version 3.9.4.10443):


H = [ 2.285927555469414e-03 -9.683344516608315e-03 0 0 ;

2.835630264340973e+03 1.995068176067767e+02 0 -3.217604282987038e-01 ;
-1.995068176067769e+02 2.835630264340972e+03 0 1.973947878435759e+00 ;
0 0 9.875839323886169e-01 -1.403672780503961e+01 ];
K = [ -6.767269167134378e-04 -2.108684486049049e-04 0 0 ;
1.424833496089810e+01 0 0 -1.608802141493519e-03 ;
0 1.424833496089809e+01 0 9.869739392178797e-03 ;
0 0 1.002472007707650e+00 0 ];

ee = eig(mp(H),mp(K)); % eigenvalue with real(ee)~100 not a conjugate pair
cplxpair(ee) % hence this fails

% Note: could use cplxpair with a tolerance, but for a real problem the
% eigenvalues should be exactly complex conjugate.
Under review

The cplxpair fails because it has hard-coded checks allowing only 'double'/'single' inputs.


I have changed the cplxpair to accept the mp-parameters as well. Will send it by email in a few sec.


I am not sure if EIG should return perfect conjugate pairs to last digits though, but it returns them grouped:

>> mp.Digits(34);
>> test_ee
ee = 
        0.9851486373638723693809204675873853 +                                       0i    
          99.9999998743453462896070441576835 +       100.00000000000015525899193951503i    
         99.99999987434534628960704415768347 -     100.0000000000001552589919395150301i    
         199.9999999999999523395305332254016 +                                       0i
ans = 
          99.9999998743453462896070441576835 -     100.00000000000015525899193951503i    
          99.9999998743453462896070441576835 +     100.00000000000015525899193951503i    
        0.9851486373638723693809204675873853 +                                     0i    
         199.9999999999999523395305332254016 +                                     0i

Second output is from modified cplxpair - it just finds and truncates pairs by tolerance (100*eps be default).

Answered

Guaranteed bit-to-bit matching has been added to generalized eigen-solver as of 3.9.4.10458 version, hence I am closing the ticket.