title: Rewriting Matlab eig(A,B) (Generalized eigenvalues/eigenvectors) to C/C++
body: <p>Do anyone have any idea how can I rewrite <code>eig(A,B)</code> from Matlab used to calculate generalized eigenvector/eigenvalues? I've been struggling with this problem lately. So far:</p>
  
  <p>Matlab definition of <code>eig</code> function I need:</p>
  
  <pre><code>[V,D] = eig(A,B) produces a diagonal matrix D of generalized
  eigenvalues and a full matrix V whose columns are the corresponding
  eigenvectors so that A*V = B*V*D.
  </code></pre>
  
  <ol>
  <li>So far I tried the <code>Eigen</code> library (<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralizedSelfAdjointEigenSolver.html" rel="nofollow">http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralizedSelfAdjointEigenSolver.html</a>)</li>
  </ol>
  
  <p>My implementation looks like this:</p>
  
  <pre><code>std::pair&lt;Matrix4cd, Vector4d&gt; eig(const Matrix4cd&amp; A, const Matrix4cd&amp; B)
  {
      Eigen::GeneralizedSelfAdjointEigenSolver&lt;Matrix4cd&gt; solver(A, B);
  
      Matrix4cd V = solver.eigenvectors();
      Vector4d D = solver.eigenvalues();
  
      return std::make_pair(V, D);
  }
  </code></pre>
  
  <p>But first thing that comes to my mind is, that I can't use <code>Vector4cd</code> as <code>.eigenvalues()</code> doesn't return complex values where Matlab does. Furthermore results of <code>.eigenvectors()</code> and <code>.eigenvalues()</code> for the same matrices are not the same at all:</p>
  
  <p>C++:</p>
  
  <pre><code>Matrix4cd x;
  Matrix4cd y;
  pair&lt;Matrix4cd, Vector4d&gt; result;
  
  for (int i = 0; i &lt; 4; i++)
  {
      for (int j = 0; j &lt; 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 &lt;&lt; result.first &lt;&lt; endl &lt;&lt; endl;
  cout &lt;&lt; result.second &lt;&lt; endl &lt;&lt; endl;
  </code></pre>
  
  <p>Matlab:</p>
  
  <pre><code>for i=1:1:4
      for 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)
  </code></pre>
  
  <p>So I give <code>eig</code> the same 4x4 matrices holding values 1-16 ascending (x) and descending (y). But I receive different results, furthermore <code>Eigen</code> method returns double from eigenvalues while Matlab returns complex dobule. I also find out that there is other <code>Eigen</code> solver named <code>GeneralizedEigenSolver</code>. That one in the documentation (<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralizedEigenSolver.html" rel="nofollow">http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralizedEigenSolver.html</a>) has written that it solves <code>A*V = B*V*D</code> but to be honest I tried it and results (matrix sizes) are not the same size as Matlab so I got quite lost how it works (examplary results are on the website I've linked). It also has only .eigenvector method.</p>
  
  <p>C++ results:</p>
  
  <pre><code>(-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
  </code></pre>
  
  <p>Matlab results:</p>
  
  <pre><code>[A,B] = eig(x,y)
  
  
  A =
  
    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
  </code></pre>
  
  <ol start="2">
  <li>Second try was with Intel IPP but it seems that it solves only <code>A*V = V*D</code> and support told me that it's not supported anymore.</li>
  </ol>
  
  <p><a href="https://software.intel.com/en-us/node/505270" rel="nofollow">https://software.intel.com/en-us/node/505270</a> (list of constructors for Intel IPP)</p>
  
  <ol start="3">
  <li>I got suggestion to move from Intel IPP to MKL. I did it and hit the wall again. I tried to check all algorithms for <code>Eigen</code> but it seems that there are only <code>A*V = V*D</code> problems solved. I was checking <code>lapack95.lib</code>. The list of algorithms used by this library is available there: 
  <a href="https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/index.htm#dsyev.htm" rel="nofollow">https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/index.htm#dsyev.htm</a></li>
  </ol>
  
  <p>Somewhere on the web I could find topic on Mathworks when someone said that managed to solve my problem partially with usage of MKL:</p>
  
  <p><a href="http://jp.mathworks.com/matlabcentral/answers/40050-generalized-eigenvalue-and-eigenvectors-differences-between-matlab-eig-a-b-and-mkl-lapack-dsygv" rel="nofollow">http://jp.mathworks.com/matlabcentral/answers/40050-generalized-eigenvalue-and-eigenvectors-differences-between-matlab-eig-a-b-and-mkl-lapack-dsygv</a></p>
  
  <p>Person said that he/she used <code>dsygv</code> algorithm but I can't locate anything like that on the web. Maybe it's a typo.</p>
  
  <p>Anyone has any other proposition/idea how can I implement it? Or maybe point my mistake. I'd appreciate that.</p>
  
  <hr>
  
  <p>EDIT:
  In comments I've received a hint that I was using <code>Eigen</code> solver wrong. My <code>A</code> matrix wasn't self-adjoint and my <code>B</code> matrix wasn't positive-definite. I took matrices from program I want to rewrite to C++ (from random iteration) and checked if they meet the requirements. They did:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>As for <code>Rj</code> which is now my <code>A</code> - it is self-adjoint because <code>Rj = Rj'</code> and <code>Rj = ctranspose(Rj)</code>. (<a href="http://mathworld.wolfram.com/Self-AdjointMatrix.html" rel="nofollow">http://mathworld.wolfram.com/Self-AdjointMatrix.html</a>)</p>
  
  <p>As for <code>Rt</code> which is now my <code>B</code> - it is Positive-Definite what is checked with method linked to me. (<a href="http://www.mathworks.com/matlabcentral/answers/101132-how-do-i-determine-if-a-matrix-is-positive-definite-using-matlab" rel="nofollow">http://www.mathworks.com/matlabcentral/answers/101132-how-do-i-determine-if-a-matrix-is-positive-definite-using-matlab</a>). So </p>
  
  <pre><code>&gt;&gt; [~,p] = chol(Rt)
  
  p =
  
       0
  </code></pre>
  
  <p>I've rewritten matrices manually to C++ and performed <code>eig(A,B)</code> again with matrices meeting requirements:</p>
  
  <pre><code>Matrix4cd x;
  Matrix4cd y;
  pair&lt;Matrix4cd, Vector4d&gt; 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 &lt;&lt; result.first &lt;&lt; endl &lt;&lt; endl;
  cout &lt;&lt; result.second &lt;&lt; endl &lt;&lt; endl;
  </code></pre>
  
  <p>And the results of C++:</p>
  
  <pre><code>(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
  </code></pre>
  
  <p>Results of MATLAB:</p>
  
  <pre><code>&gt;&gt;  [A,B] = eig(Rj,Rt)
  
  A =
  
    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
  </code></pre>
  
  <p><code>Eigenvalues</code> are the same! Nice, but why <code>Eigenvectors</code> are not similar at all?</p>