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

Question

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

4

有人知道如何从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);
}

但我首先想到的是，由于 .eigenvalues() 不返回复数值，所以我无法使用 Vector4cd，而 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

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

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

有人建议我从Intel 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

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

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

关于现在我的A是Rj - 它是自伴的，因为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

"特征值是相同的！非常好，但为什么特征向量完全不相似呢？"

- F1sher

1

随意尝试库并不是一种有效的方法。Eigen是一个全面而精美设计的库，文档也非常出色。你确定它不适合吗？ - David Heffernan

你好，我已经花费了约10个小时在Eigen上。从他们的文档http://eigen.tuxfamily.org/dox/group__Eigenvalues__Module.html中可以看出，有几个Eigen求解器，只有两个用于广义情况。我已经实现了它们并插入了与Matlab相同的矩阵。但是在两种情况下结果都不同+我遇到了输出数据结构和类型问题。我已经阅读了整个文档。通常我会尝试找到线索，然后一直跟随它直到我没有想法，然后尝试不同的方法。这比无所事事地坐着不知道下一步该做什么要好得多。 - F1sher

2

一些注意事项：您的 y 矩阵不是 Eigen 所需的正定矩阵（“该类解决广义特征值问题 $ Av = \lambda Bv $。在这种情况下，矩阵 $ A $ 应为自伴随矩阵，矩阵 $ B $ 应为正定矩阵。”）。请参见此处以确定矩阵是否为 PD。同样，x 不是自伴随矩阵，但我已经没有足够的空间来进行评论。 - Avi Ginsburg

谢谢注意到这个问题。我会调查一下。 - F1sher

我已经研究了正定和自伴随条件。我对不适合进行测试的矩阵进行了测试。我从新的测试中更新了信息，特征值正在工作。但是我不明白为什么特征向量没有给出相同的结果。这一定是算法差异导致的，不是吗？ - F1sher

为什么不尝试一些更简单的东西，比如SVD，并调试您的代码呢？ - μολὼν.λαβέ

2个回答

0

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

- Schulzey

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Stefano M · Accepted Answer

这里Eigen没有问题。实际上，在第二个示例运行中，Matlab和Eigen产生了完全相同的结果。请从基本线性代数中记住特征向量是由任意缩放因子决定的（即如果v是一个特征向量，那么对于非零复标量α*v也成立）。不同的线性代数库计算不同的特征向量是很常见的，但这并不意味着两个代码中有一个是错误的：这只是意味着它们选择了不同的特征向量缩放方式。

编辑

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

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

但没有强制施加其他限制。不幸的是，这意味着比例不是唯一的，并且可能取决于通用特征向量算法的中间结果。在这种情况下，我无法为您提供有关如何精确复制Matlab的建议：需要了解Matlab的内部工作原理。

总的来说，在线性代数中，人们通常不太关心特征向量缩放，因为这通常对所解决的问题完全无关，当特征向量只是用作中间结果时。唯一需要精确定义比例的情况是当您要给出特征值的图形表示时。