将 Matlab eig(A,B)（广义特征值/特征向量）重写为 C/C++答案

【问题标题】：Rewriting Matlab eig(A,B) (Generalized eigenvalues/eigenvectors) to C/C++将 Matlab eig(A,B)（广义特征值/特征向量）重写为 C/C++
【发布时间】：2023-03-21 20:43:01
【问题描述】：

有谁知道如何从用于计算广义特征向量/特征值的 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 那样返回复杂值。此外，.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 方法从特征值返回双精度，而 Matlab 返回复数双精度。我还发现还有其他名为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

第二次尝试是使用 Intel 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

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

任何人有任何其他提议/想法我该如何实施？或者也许指出我的错误。我会很感激的。

编辑：在 cmets 中，我收到一个提示，提示我使用 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 一点都不相似？

【问题讨论】：

随意尝试库是没有效率的。 Eigen 是一个全面且设计精美的库，具有出色的文档。确定不合适？
您好，我已经在 Eigen 上花费了大约 10 个小时。查看他们的文档eigen.tuxfamily.org/dox/group__Eigenvalues__Module.html 我可以看到有一些 Eigen 求解器，其中只有两个用于一般情况。我已经实现了它们并插入了与在 Matlab 中相同的矩阵。在这两种情况下，结果都不同+我对输出数据的结构和类型有疑问。我阅读了整个文档。通常我会试图找到一个线索，然后我会跟着它，直到我没有想法，然后尝试不同的东西。总比坐着不知道下一步该做什么要好。
一些注意事项：您的y 矩阵不是正定矩阵，正如Eigen ("该类解决了广义特征值问题$ Av = \lambda Bv $。在这种情况下，矩阵$ A $ 应该是自伴随的，并且矩阵 $ B $ 应该是正定的。”）。请参阅here 以确定矩阵是否为 PD。同样，x 不是自伴随的，但我没有空间发表此评论。
感谢您的关注。我会调查的。
我研究了正定和自伴条件。我正在对不适合它的矩阵进行测试。我已经从我的新测试中更新了信息，并且特征值正在工作。我不明白为什么特征向量没有给出相同的结果。那一定是一些算法差异，不是吗？

标签： c++ matlab linear-algebra eigen intel-mkl

【解决方案1】：

这里的 Eigen 没有问题。

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

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

编辑

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

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

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

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

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

【讨论】：

请注意：这是一个反复出现的问题，请参阅。 stackoverflow.com/q/17682626/1499402
感谢您的帮助。有什么方法可以将 Matlab 使用的缩放比例与我的库结果之一相匹配。让我们说特征？
请参阅上面的编辑以获取更多问题的答案。

【解决方案2】：

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

【讨论】：