将 Matlab eig(A,B)（广义 eigenvalues/eigenvectors）重写为 C/C++

Question

有谁知道如何从用于计算广义 eigenvector/eigenvalues 的 Matlab 重写 eig(A,B)？我最近一直在努力解决这个问题。到目前为止：

我需要的 eig 函数的 Matlab 定义：

[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.

到目前为止，我尝试了 Eigen 库 (http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralizedSelfAdjointEigenSolver.html)

我的实现是这样的：

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);
}

但我首先想到的是，我不能使用 Vector4cd，因为 .eigenvalues() 不像 Matlab 那样 return 复杂值。此外，相同矩阵的 .eigenvectors() 和 .eigenvalues() 的结果完全不同：

C++:

Matrix4cd x;
Matrix4cd y;
pair<Matrix4cd, Vector4d> result;

for (int i = 0; i < 4; i++)
{
    for (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
    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)

所以我给出 eig 相同的 4x4 矩阵，其中包含值 1-16 升序 (x) 和降序 (y)。但是我收到了不同的结果，此外 Eigen 方法 returns 从特征值加倍，而 Matlab returns 复数加倍。我还发现还有另一个名为 GeneralizedEigenSolver 的 Eigen 求解器。文档中的那个 (http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralizedEigenSolver.html) 写道它解决了 A*V = B*V*D 但老实说我试过了结果（矩阵大小）与 Matlab 的大小不一样所以我完全不知道它是如何工作的（示例结果在我链接的网站上）。它也只有 .eigenvector 方法。

C++ 结果：

(-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 结果：

[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

第二次尝试使用英特尔 IPP，但它似乎只能解决 A*V = V*D 并且支持人员告诉我它不再受支持。

https://software.intel.com/en-us/node/505270（英特尔 IPP 的构造函数列表）

我收到了从英特尔 IPP 转移到 MKL 的建议。我做到了，又撞墙了。我试图检查 Eigen 的所有算法，但似乎只解决了 A*V = V*D 个问题。我正在检查 lapack95.lib。该库使用的算法列表可在此处获得： https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/index.htm#dsyev.htm

当有人说使用 MKL 成功地部分解决了我的问题时，我在网上的某个地方找到了关于 Mathworks 的主题：

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

有人说 he/she 使用了 dsygv 算法，但我在网上找不到类似的东西。可能打错了。

谁有其他的 proposition/idea 我该如何实施？或者也许指出我的错误。非常感谢。

编辑：在评论中，我收到一条提示，提示我错误地使用了 Eigen 求解器。我的 A 矩阵不是自伴矩阵，我的 B 矩阵不是正定矩阵。我从我想重写为 C++ 的程序中取出矩阵（从随机迭代）并检查它们是否满足要求。他们做到了：

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

至于 Rj 现在是我的 A - 它是自伴的，因为 Rj = Rj' 和 Rj = ctranspose(Rj)。 (http://mathworld.wolfram.com/Self-AdjointMatrix.html)

至于 Rt 现在是我的 B - 它是正定的，用与我相关的方法检查。 (http://www.mathworks.com/matlabcentral/answers/101132-how-do-i-determine-if-a-matrix-is-positive-definite-using-matlab)。所以

>> [~,p] = chol(Rt)

p =

     0

我已将矩阵手动重写为 C++ 并再次执行 eig(A,B)，矩阵满足要求：

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;

而 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

MATLAB 的结果：

>>  [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

Eigenvalues都一样！不错，但为什么 Eigenvectors 一点都不相似？

Answer 1

这里Eigen没问题

事实上，对于第二个例子运行，Matlab 和 Eigen 产生了完全相同的结果。请记住基本线性代数中的特征向量是由任意比例因子决定的。（即，如果 v 是特征向量，则同样适用于 alpha*v，其中 alpha 是非零复数标量。）

不同的线性代数库计算不同的特征向量是很常见的，但这并不意味着两个代码之一是错误的：它只是意味着他们选择了不同的特征向量缩放比例。

编辑

精确复制 matlab 选择的缩放比例的主要问题是 eig(A,B) 是一个 驱动程序 例程，这取决于 A 的不同属性B 可能会调用不同的 libraries/routines，并应用额外的步骤，例如平衡矩阵等。通过快速检查您的示例，我会说在这种情况下 matlab 强制执行以下条件：

all(imag(V(end,:))==0)（每个特征向量的最后一个分量是实数）

但不施加其他限制。不幸的是，这意味着缩放不是唯一的，并且可能取决于所使用的广义特征向量算法的中间结果。在这种情况下，我无法就如何完全复制 matlab 向您提供建议：需要了解 matlab 的内部工作原理。

一般来说，在线性代数中，人们通常不太关心特征向量缩放，因为当特征向量仅用作中间结果时，这通常与所解决的问题完全无关。

唯一需要精确定义缩放比例的情况是当您要给出特征值的图形表示时。

Answer 2

Matlab 中的特征向量缩放似乎是基于将它们归一化为 1.0（即每个向量中最大项的绝对值为 1.0）。在应用程序中，我还使用了 returns 左特征向量，而不是更常用的右特征向量。这可以解释 Matlab 和 Lapack MKL 中的特征求解器之间的差异。

将 Matlab eig(A,B)（广义 eigenvalues/eigenvectors）重写为 C/C++

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

c++

matlab

linear-algebra

intel-mkl

eigen