Rewriting Matlab eig(A,B) (Generalized eigenvalues/eigenvectors) to C/C++ -

do have idea how can rewrite eig(a,b) matlab used calculate generalized eigenvector/eigenvalues? i've been struggling problem lately. far:

matlab definition of eig function need:

[v,d] = eig(a,b) produces diagonal matrix d of generalized eigenvalues , full matrix v columns corresponding eigenvectors a*v = b*v*d. 

1. so far tried eigen library (http://eigen.tuxfamily.org/dox/classeigen_1_1generalizedselfadjointeigensolver.html)

my implementation looks this:

std::pair<matrix4cd, vector4d> eig(const matrix4cd& a, const matrix4cd& b) {     eigen::generalizedselfadjointeigensolver<matrix4cd> solver(a, b);      matrix4cd v = solver.eigenvectors();     vector4d d = solver.eigenvalues();      return std::make_pair(v, d); } 

but first thing comes mind is, can't use vector4cd .eigenvalues() doesn't return complex values matlab does. furthermore results of .eigenvectors() , .eigenvalues() same matrices not same @ all:

c++:

matrix4cd x; matrix4cd y; pair<matrix4cd, vector4d> result;  (int = 0; < 4; i++) {     (int j = 0; j < 4; j++)     {         x.real()(i,j) = (double)(i+j+1+i*3);         y.real()(i,j) = (double)(17 - (i+j+1+i*3));          x.imag()(i,j) = (double)(i+j+1+i*3);         y.imag()(i,j) = (double)(17 - (i+j+1+i*3));     } } result = eig(x,y); cout << result.first << endl << endl; cout << result.second << endl << endl; 

matlab:

for i=1:1:4     j=1:1:4         x(i,j) = complex((i-1)+(j-1)+1+((i-1)*3), (i-1)+(j-1)+1+((i-1)*3));         y(i,j) = complex(17 - ((i-1)+(j-1)+1+((i-1)*3)), 17 - ((i-1)+(j-1)+1+((i-1)*3)));     end end  [a,b] = eig(x,y) 

so give eig same 4x4 matrices holding values 1-16 ascending (x) , descending (y). receive different results, furthermore eigen method returns double eigenvalues while matlab returns complex dobule. find out there other eigen solver named generalizedeigensolver. 1 in documentation (http://eigen.tuxfamily.org/dox/classeigen_1_1generalizedeigensolver.html) has written solves a*v = b*v*d honest tried , results (matrix sizes) not same size matlab got quite lost how works (examplary results on website i've linked). has .eigenvector method.

c++ results:

(-0.222268,-0.0108754) (0.0803437,-0.0254809) (0.0383264,-0.0233819) (0.0995482,0.00682079) (-0.009275,-0.0182668) (-0.0395551,-0.0582127) (0.0550395,0.03434) (-0.034419,-0.0287563) (-0.112716,-0.0621061) (-0.010788,0.10297) (-0.0820552,0.0294896) (-0.114596,-0.146384) (0.28873,0.257988) (0.0166259,-0.0529934) (0.0351645,-0.0322988) (0.405394,0.424698)  -1.66983 -0.0733194 0.0386832 3.97933 

matlab results:

[a,b] = eig(x,y)   =    columns 1 through 3    -0.9100 + 0.0900i  -0.5506 + 0.4494i   0.3614 + 0.3531i    0.7123 + 0.0734i   0.4928 - 0.2586i  -0.5663 - 0.4337i    0.0899 - 0.4170i  -0.1210 - 0.3087i   0.0484 - 0.1918i    0.1077 + 0.2535i   0.1787 + 0.1179i   0.1565 + 0.2724i    column 4    -0.3237 - 0.3868i    0.2338 + 0.7662i    0.5036 - 0.3720i   -0.4136 - 0.0074i   b =    columns 1 through 3    -1.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i    0.0000 + 0.0000i  -1.0000 - 0.0000i   0.0000 + 0.0000i    0.0000 + 0.0000i   0.0000 + 0.0000i  -4.5745 - 1.8929i    0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i    column 4     0.0000 + 0.0000i    0.0000 + 0.0000i    0.0000 + 0.0000i   -0.3317 + 1.1948i 

2. second try intel ipp seems solves a*v = v*d , support told me it's not supported anymore.

https://software.intel.com/en-us/node/505270 (list of constructors intel ipp)

3. i got suggestion move intel ipp mkl. did , hit wall again. tried check algorithms eigen seems there a*v = v*d problems solved. checking lapack95.lib. list of algorithms used library available there: https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/index.htm#dsyev.htm

somewhere on web find topic on mathworks when said managed solve problem partially usage of mkl:

http://jp.mathworks.com/matlabcentral/answers/40050-generalized-eigenvalue-and-eigenvectors-differences-between-matlab-eig-a-b-and-mkl-lapack-dsygv

person said he/she used dsygv algorithm can't locate on web. maybe it's typo.

anyone has other proposition/idea how can implement it? or maybe point mistake. i'd appreciate that.

---

edit: in comments i've received hint using eigen solver wrong. a matrix wasn't self-adjoint , b matrix wasn't positive-definite. took matrices program want rewrite c++ (from random iteration) , checked if meet requirements. did:

rj =    1.0e+02 *   columns 1 through 3     0.1302 + 0.0000i  -0.0153 + 0.0724i   0.0011 - 0.0042i   -0.0153 - 0.0724i   1.2041 + 0.0000i  -0.0524 + 0.0377i    0.0011 + 0.0042i  -0.0524 - 0.0377i   0.0477 + 0.0000i   -0.0080 - 0.0108i   0.0929 - 0.0115i  -0.0055 + 0.0021i    column 4    -0.0080 + 0.0108i    0.0929 + 0.0115i   -0.0055 - 0.0021i    0.0317 + 0.0000i  rt =    columns 1 through 3     4.8156 + 0.0000i  -0.3397 + 1.3502i  -0.2143 - 0.3593i   -0.3397 - 1.3502i   7.3635 + 0.0000i  -0.5539 - 0.5176i   -0.2143 + 0.3593i  -0.5539 + 0.5176i   1.7801 + 0.0000i    0.5241 + 0.9105i   0.9514 + 0.6572i  -0.7302 + 0.3161i    column 4     0.5241 - 0.9105i    0.9514 - 0.6572i   -0.7302 - 0.3161i    9.6022 + 0.0000i 

as rj a - self-adjoint because rj = rj' , rj = ctranspose(rj). (http://mathworld.wolfram.com/self-adjointmatrix.html)

as rt b - positive-definite checked method linked me. (http://www.mathworks.com/matlabcentral/answers/101132-how-do-i-determine-if-a-matrix-is-positive-definite-using-matlab).

>> [~,p] = chol(rt)  p =       0 

i've rewritten matrices manually c++ , performed eig(a,b) again matrices meeting requirements:

matrix4cd x; matrix4cd y; pair<matrix4cd, vector4d> result;  x.real()(0,0) = 13.0163601949795; x.real()(0,1) = -1.53172561296005; x.real()(0,2) = 0.109594869350436; x.real()(0,3) = -0.804231869422614;  x.real()(1,0) = -1.53172561296005; x.real()(1,1) = 120.406645675346; x.real()(1,2) = -5.23758765476463; x.real()(1,3) = 9.28686785230169;  x.real()(2,0) = 0.109594869350436;  x.real()(2,1) = -5.23758765476463; x.real()(2,2) = 4.76648319080400; x.real()(2,3) = -0.552823839520508;  x.real()(3,0) = -0.804231869422614; x.real()(3,1) = 9.28686785230169; x.real()(3,2) = -0.552823839520508; x.real()(3,3) = 3.16510496622613;  x.imag()(0,0) = -0.00000000000000; x.imag()(0,1) = 7.23946944213164; x.imag()(0,2) = 0.419181335323979; x.imag()(0,3) = 1.08441894337449;  x.imag()(1,0) = -7.23946944213164; x.imag()(1,1) = -0.00000000000000; x.imag()(1,2) = 3.76849276970080; x.imag()(1,3) = 1.14635625342266;  x.imag()(2,0) = 0.419181335323979; x.imag()(2,1) = -3.76849276970080; x.imag()(2,2) = -0.00000000000000; x.imag()(2,3) = 0.205129702522089;  x.imag()(3,0) = -1.08441894337449; x.imag()(3,1) = -1.14635625342266; x.imag()(3,2) = 0.205129702522089; x.imag()(3,3) = -0.00000000000000;  y.real()(0,0) = 4.81562784930907; y.real()(0,1) = -0.339731222392148; y.real()(0,2) = -0.214319720979258; y.real()(0,3) = 0.524107127885349;  y.real()(1,0) = -0.339731222392148; y.real()(1,1) = 7.36354235698375; y.real()(1,2) = -0.553927983436786; y.real()(1,3) = 0.951404408649307;  y.real()(2,0) = -0.214319720979258; y.real()(2,1) = -0.553927983436786; y.real()(2,2) = 1.78008768533745; y.real()(2,3) = -0.730246631850385;  y.real()(3,0) = 0.524107127885349; y.real()(3,1) = 0.951404408649307; y.real()(3,2) = -0.730246631850385; y.real()(3,3) = 9.60215057284395;  y.imag()(0,0) = -0.00000000000000; y.imag()(0,1) = 1.35016928394966; y.imag()(0,2) = -0.359262708214312; y.imag()(0,3) = -0.910512495060186;  y.imag()(1,0) = -1.35016928394966; y.imag()(1,1) = -0.00000000000000; y.imag()(1,2) = -0.517616473138836; y.imag()(1,3) = -0.657235460367660;  y.imag()(2,0) = 0.359262708214312; y.imag()(2,1) = 0.517616473138836; y.imag()(2,2) = -0.00000000000000; y.imag()(2,3) = -0.316090662865005;  y.imag()(3,0) = 0.910512495060186; y.imag()(3,1) = 0.657235460367660; y.imag()(3,2) = 0.316090662865005; y.imag()(3,3) = -0.00000000000000;  result = eig(x,y);  cout << result.first << endl << endl; cout << result.second << endl << endl; 

and results of c++:

(0.0295948,0.00562174) (-0.253532,0.0138373) (-0.395087,-0.0139696) (-0.0918132,-0.0788735) (-0.00994614,-0.0213973) (-0.0118322,-0.0445976) (0.00993512,0.0127006) (0.0590018,-0.387949) (0.0139485,-0.00832193) (0.363694,-0.446652) (-0.319168,0.376483) (-0.234447,-0.0859585) (0.173697,0.268015) (0.0279387,-0.0103741) (0.0273701,0.0937148)  (-0.055169,0.0295393)  0.244233 2.24309 3.24152 18.664 

results of matlab:

>>  [a,b] = eig(rj,rt)  =    columns 1 through 3     0.0208 - 0.0218i   0.2425 + 0.0753i  -0.1242 + 0.3753i   -0.0234 - 0.0033i  -0.0044 + 0.0459i   0.0150 - 0.0060i    0.0006 - 0.0162i  -0.4964 + 0.2921i   0.2719 + 0.4119i    0.3194 + 0.0000i  -0.0298 + 0.0000i   0.0976 + 0.0000i    column 4    -0.0437 - 0.1129i    0.2351 - 0.3142i   -0.1661 - 0.1864i   -0.0626 + 0.0000i  b =     0.2442         0         0         0         0    2.2431         0         0         0         0    3.2415         0         0         0         0   18.6640 

eigenvalues same! nice, why eigenvectors not similar @ all?

there no problem here eigen.

in fact second example run, matlab , eigen produced same result. please remember basic linear algebra eigenvector determined arbitrary scaling factor. (i.e. if v eigenvector same holds alpha*v, alpha non 0 complex scalar.)

it quite common different linear algebra libraries compute different eigenvectors, not mean 1 of 2 codes wrong: means choose different scaling of eigenvectors.

edit

the main problem replicating scaling chosen matlab eig(a,b) driver routine, depending different properties of a , b may call different libraries/routines, , apply steps balancing matrices , on. inspecting example, in case matlab enforcing following condition:

• all(imag(v(end,:))==0) (the last component of each eigenvector real)

but not imposing other constraints. unfortunately means scaling not unique, , depends on intermediate results of generalised eigenvector algorithm used. in case i'm not able give advice on how exactly replicate matlab: knowledge of internal working of matlab required.

as general remark, in linear algebra 1 not care eigenvector scaling, since irrelevant problem solved, when eigenvectors used intermediate results.

the case in scaling has defined exactly, when going give graphic representation of eigenvalues.