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

Question

因此我需要将矩阵右手除法从 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 returns (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 实现矩阵 "slash" 除法？
如果不是 Eigen，那么你能推荐任何简单的方法吗？

Answer 1

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

Answer 2

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

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

A.fullPivLu().solve(b);

使用 LU 会给你 A\b。

Answer 3

这里：

#include <Eigen/Dense>

Eigen::MatrixXd inv_eigen (Eigen::MatrixXd m)
{
    // Matlab: inv(X) = X^(-1) = X\speye(size(X))
    return m.fullPivLu().solve(Eigen::MatrixXd::Identity(m.rows(), m.cols()));
}

结果：

X = 3×3

  1     0     2
 -1     5     0
  0     3    -9

Matlab inv(X)

  0.8824   -0.1176    0.1961
  0.1765    0.1765    0.0392
  0.0588    0.0588   -0.0980

Eigen inv(X)

  0.882353  -0.117647   0.196078
  0.176471   0.176471   0.0392157
  0.0588235  0.0588235 -0.0980392

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

Implementation of Matlab matrix inverse function in C++ with Eigen

c++

matlab

matrix

eigen