用Eigen在C++中实现Matlab矩阵逆函数

Posted 2023-02-16

【中文标题】用Eigen在C++中实现Matlab矩阵逆函数

【英文标题】：Implementation of Matlab matrix inverse function in C++ with Eigen

【发布时间】：2015-08-11 08:28:10

【问题描述】：

所以我需要将矩阵右手除法从 Matlab 重写为 C++：

At = (xPow*yPow')/(yPow*yPow'); 

我模拟了一些矩阵：

>> xPow*yPow'

ans =

0.0004    0.0040    0.0004    0.0004
0.0014    0.0263    0.0014    0.0014
0.0004    0.0012    0.0004    0.0004
0.0012    0.0053    0.0012    0.0012

和

>> yPow*yPow'

ans =

0.0001    0.0004    0.0001    0.0001
0.0004    0.0256    0.0004    0.0004
0.0001    0.0004    0.0001    0.0001
0.0001    0.0004    0.0001    0.0001

Matlab 为 (xPow*yPow')/(yPow*yPow') 和 xPow*yPow' * inv(yPow*yPow') 返回相同的结果。

>> xPow*yPow' * inv(yPow*yPow')

ans =

36.1259    0.1127  -30.3163   -2.6999
40.6472    0.8810  -19.7529  -11.8430
-1.5578   -0.0182   12.1397   -7.0087
124.4466    0.0594 -130.0163   16.6710

>> At = (xPow*yPow')/(yPow*yPow')

At =

 36.1259    0.1127  -30.3163   -2.6999
 40.6472    0.8810  -19.7529  -11.8430
 -1.5578   -0.0182   12.1397   -7.0087
124.4466    0.0594 -130.0163   16.6710

Eigen 库有函数 .inverse() 所以我想我可以用它来实现这些矩阵的除法：

xyPowMult(0,0) = 0.0004;
xyPowMult(0,1) = 0.0040;
xyPowMult(0,2) = 0.0004;
xyPowMult(0,3) = 0.0004;
xyPowMult(1,0) = 0.0014;
xyPowMult(1,1) = 0.0263;
xyPowMult(1,2) = 0.0014;
xyPowMult(1,3) = 0.0014;
xyPowMult(2,0) = 0.0004;
xyPowMult(2,1) = 0.0012;
xyPowMult(2,2) = 0.0004;
xyPowMult(2,3) = 0.0004;
xyPowMult(3,0) = 0.0012;
xyPowMult(3,1) = 0.0053;
xyPowMult(3,2) = 0.0012;
xyPowMult(3,3) = 0.0012;

yyPowMult(0,0) = 0.0001;
yyPowMult(0,1) = 0.0004;
yyPowMult(0,2) = 0.0001;
yyPowMult(0,3) = 0.0001;
yyPowMult(1,0) = 0.0004;
yyPowMult(1,1) = 0.0256;
yyPowMult(1,2) = 0.0004;
yyPowMult(1,3) = 0.0004;
yyPowMult(2,0) = 0.0001;
yyPowMult(2,1) = 0.0004;
yyPowMult(2,2) = 0.0001;
yyPowMult(2,3) = 0.0001;
yyPowMult(3,0) = 0.0001;
yyPowMult(3,1) = 0.0004;
yyPowMult(3,2) = 0.0001;
yyPowMult(3,3) = 0.0001;

AtTemp = xyPowMult * yyPowMult.inverse();

cout << " x*y' " << endl;
cout << xyPowMult << endl << endl;

cout << " y*y' " << endl;
cout << yyPowMult << endl << endl;

cout << " x*y' * inv(y*y') " << endl;
cout << xyPowMult * yyPowMult.inverse() << endl << endl;

控制台结果显示不同的结果：

 x*y'
0.0004  0.004 0.0004 0.0004
0.0014 0.0263 0.0014 0.0014
0.0004 0.0012 0.0004 0.0004
0.0012 0.0053 0.0012 0.0012

 y*y'
0.0001 0.0004 0.0001 0.0001
0.0004 0.0256 0.0004 0.0004
0.0001 0.0004 0.0001 0.0001
0.0001 0.0004 0.0001 0.0001

 x*y' * inv(y*y')
-1.#IND -1.#IND -1.#IND -1.#IND
-1.#IND -1.#IND -1.#IND -1.#IND
-1.#IND -1.#IND -1.#IND -1.#IND
-1.#IND -1.#IND -1.#IND -1.#IND

所以我的问题是：

为什么结果不同？ 如何使用 Eigen 实现矩阵“斜线”除法？ 如果不是 Eigen，您能否提出任何简单的方法？

【参考方案1】：

你有两个问题。首先，正如 Ilya Popov 指出的，y*y' 是单数的。实际上，x*y' 也是如此。其次，matlab的\operator实际上求解了一个线性方程组（Ax=B -->求解x）。使用简单的方法（例如inverse()）求解奇异（或接近奇异）矩阵通常不是一个好主意。

因此，要对 Eigen 做同样的事情，您需要建立方程来求解并使用解 (intro)。但是，您将根据先验知识选择特定算法。例如，

A.fullPivLu().solve(b);

会给你A\b使用LU。

【参考方案2】：

如您所见，y*y' 的最后两行是相同的。矩阵是奇异的（就像v*v' 这样的任何矩阵一样）。这样的矩阵没有逆矩阵。所以你不能在这里期待任何有意义的结果。奇怪的是，Matlab 产生任何东西。你确定你的公式是正确的吗？

【讨论】：

可能是因为我自己将数字插入到 C++ 中，而在 Matlab 中它们是从正确的乘法中获得的。 Matlab 显示的只是我重写为 C++ 的近似值。在 Matlab 中，第 4 个数字之后的数字可能不同，但我只重写了 4 个，所以我创建了奇异矩阵？

是的，可以。请注意，如果矩阵接近奇异矩阵（又名病态，具有较大的条件数），那么即使其系数中的最小误差也会导致逆矩阵出现较大误差。