c++ - 将 Matlab eig(A,B)(广义特征值/特征向量)重写为 C/C++

Question

有谁知道如何重写 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.

1. 到目前为止，我尝试了 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 那样返回复数值。此外，对于相同的矩阵，.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 方法从特征值返回 double 值，而 Matlab 返回复数 double 值。我还发现还有一个名为 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

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

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

1. 有人建议我从英特尔 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

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

任何人有任何其他建议/想法我该如何实现？或者也许指出我的错误。我将不胜感激。

---

编辑:在评论中，我收到一条提示，提示我错误地使用了 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 )

至于现在是我的 B 的 Rt - 使用与我相关的方法检查的内容是正定的。 ( 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

特征值 是一样的!很好，但为什么 特征向量 一点都不相似？

Answer 1

Eigen这里没有问题。

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

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

编辑

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

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

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

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

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