Eigen中求逆矩阵的计算出现错误

Question

我正在尝试构建一个简单的输入/输出矩阵（可以计算需求增加时简单经济中的乘数效应）。但是出于某种原因，最终结果没有加起来。

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;  

void InputOutput(){
MatrixXf ProdA(5, 5);;
VectorXf Intd(5);
VectorXf Finald(5);
ProdA <<
    10, 20, 0, 0, 5,
    20, 30, 20, 10, 10,
    10, 10, 0, 10, 10,
    10, 40, 20, 5, 5,
    20, 20, 30, 5, 5;

Intd << 55, 40, 20, 30, 10;

Finald << 0, 0, 0, 0, 0;

VectorXf ones(5);
ones << 1, 1, 1, 1, 1;

Finald = ProdA * ones + Intd;

MatrixXf AMatrix = MatrixXf::Zero(ProdA.rows(), ProdA.cols());

AMatrix = ProdA.array() / (Finald.replicate(1, ProdA.cols())).array();
cout << "Here is the Coefficient vector production needed:\n" << AMatrix << endl;

MatrixXf IminA(5, 5);;

IminA = MatrixXf::Identity(AMatrix.rows(), AMatrix.cols()) - AMatrix;

cout << "Here is the matrix of production:\n" << ProdA << endl;
cout << "Here is the vector Internal demand:\n" << Intd << endl;
cout << "Here is the vector Final demand:\n" << Finald << endl;
cout << "Here is the Coefficient vector production needed:\n" << AMatrix << endl;

MatrixXf IminAinv(5, 5);;
IminAinv = IminA.inverse();

cout << "The inverse of CoMatrix - Imatrix is:\n" << IminAinv << endl;

cout << "To check, final demand is:\n" << (IminAinv * Intd) << endl;

当我验证(I-A)逆矩阵(IminAinv)是否正确计算时，结果不符合预期。通过将IminAinv乘以内部需求(int)，我应该得到相同的Intd，前提是Intd没有改变。但实际上我得到的数字更大。另外，如果我自己计算IminA矩阵的逆矩阵，得到的结果与特征值不同。

因此，在获取单位矩阵-系数矩阵的逆矩阵时出现了问题。但是具体问题是什么呢？

谢谢！

- TheProgramMAN123

---

请注意，如果 IminA 接近奇异（病态）状态，则使用 IminA.inverse(); 是危险的。有其他可以安全使用的方法，请参阅文档以了解详情。 - Avi Ginsburg

Answer 1

编辑： 经过进一步的研究，我发现在第二种情况下提到的那些“潜在机制”实际上是我自己在输入矩阵值时不小心犯的错误。以下是已经修正了这些错误的原始答案。

---

实际问题并不在于反转AMatrix，而在于一个更加微妙的细节。 您正在使用此命令执行AMatrix定义中的除法操作：

AMatrix = ProdA.array() / (Finald.replicate(1, ProdA.cols())).array();

但是，如果您检查对Finald进行的此复制操作的结果，则会得到：

...
cout << "Here is the replicated final demand vector:\n" << (Finald.replicate(1, ProdA.cols())).array() << endl;    
...
>>
Here is the replicated final demand vector:
     90  90  90  90  90
    130 130 130 130 130
     60  60  60  60  60
    110 110 110 110 110
     90  90  90  90  90

而正确的应该是：

90   130    60   110    90
90   130    60   110    90
90   130    60   110    90
90   130    60   110    90
90   130    60   110    90

你可以这样转置复制的最终需求向量：

MatrixXf Finaldrep(5,5);
Finaldrep = (Finald.replicate(1, ProdA.cols())).array().transpose();

当然还有：

AMatrix = ProdA.array() / Finaldrep.array();

得到：

cout << "Here is the transposed replicated final demand vector:\n" << Finaldrep << endl;
...
>>
Here is the transposed replicated final demand vector:
 90 130  60 110  90
 90 130  60 110  90
 90 130  60 110  90
 90 130  60 110  90
 90 130  60 110  90

所以，让我们看看这两种情况下你的中间和最终结果有什么区别：

情况一

即您目前的方法

Here is the Coefficient vector production needed:
 0.111111  0.222222         0         0 0.0555556
 0.153846  0.230769  0.153846 0.0769231 0.0769231
 0.166667  0.166667         0  0.166667  0.166667
0.0909091  0.363636  0.181818 0.0454545 0.0454545
 0.222222  0.222222  0.333333 0.0555556 0.0555556
The determinant of IminA is: 0.420962
The inverse of CoMatrix - Imatrix is:
  1.27266  0.468904  0.131153 0.0688064   0.13951
 0.443909   1.68132  0.377871  0.215443  0.240105
 0.451292  0.628205   1.25318  0.287633  0.312705
 0.404225  0.841827  0.423093   1.20242  0.224877
 0.586957  0.777174  0.586957   0.23913   1.27174
To check, final demand is:
94.8349
 108.09
86.7689
102.689
     95

我也添加了 IminA 的行列式。

第二种情况，即使用反转的最终需求向量。

Here is the Coefficient vector production needed:
 0.111111  0.153846         0         0 0.0555556
 0.222222  0.230769  0.333333 0.0909091  0.111111
 0.111111 0.0769231         0 0.0909091  0.111111
 0.111111  0.307692  0.333333 0.0454545 0.0555556
 0.222222  0.153846       0.5 0.0454545 0.0555556
The determinant of IminA is: 0.420962
The inverse of CoMatrix - Imatrix is:
  1.27266  0.324626  0.196729 0.0562962   0.13951
 0.641202   1.68132  0.818721  0.254615  0.346818
 0.300861  0.289941   1.25318  0.156891   0.20847
 0.494053  0.712316   0.77567   1.20242   0.27485
 0.586957  0.538044  0.880435  0.195652   1.27174
To check, final demand is:
 90
130
 60
110
 90

现在，我理解Finald check仍然不能产生最初定义的Finald的确切值，但我相信这与精度或其他基本机制有关（见注释）。

作为概念证明，以下是使用MATLAB获得的一些结果，使用第二种情况（反转）用于replicated Final Demand Vector（分母）：

>> AMatrixcm = ProdA ./ Finaldfullcm

AMatrixcm =

    0.1111    0.1538         0         0    0.0556
    0.2222    0.2308    0.3333    0.0909    0.1111
    0.1111    0.0769         0    0.0909    0.1111
    0.1111    0.3077    0.3333    0.0455    0.0556
    0.2222    0.1538    0.5000    0.0455    0.0556

>> IminAcm = eye(5) - AMatrixcm

IminAcm =

    0.8889   -0.1538         0         0   -0.0556
   -0.2222    0.7692   -0.3333   -0.0909   -0.1111
   -0.1111   -0.0769    1.0000   -0.0909   -0.1111
   -0.1111   -0.3077   -0.3333    0.9545   -0.0556
   -0.2222   -0.1538   -0.5000   -0.0455    0.9444

>> det(IminAcm)

ans =

    0.4210

>> IminAinvcm = inv(IminAcm)

IminAinvcm =

    1.2727    0.3246    0.1967    0.0563    0.1395
    0.6412    1.6813    0.8187    0.2546    0.3468
    0.3009    0.2899    1.2532    0.1569    0.2085
    0.4941    0.7123    0.7757    1.2024    0.2748
    0.5870    0.5380    0.8804    0.1957    1.2717

>> Finaldcheckcm = IminAinvcm * Intdc

Finaldcheckcm =

   90.0000
  130.0000
   60.0000
  110.0000
   90.0000

很明显，第二种情况的结果与MATLAB的结果（几乎）完全相同。请注意：在这里，您可以看到MATLAB输出与原始Finald完全相同，但是，如果您手动执行最后一个矩阵乘法（验证最终需求向量中的那个），您将看到实际上Case 2版本的IminAinv与MATLAB产生的相同结果，即Case 2的最终输出[88.9219，125.728，59.5037，105.543，84.5808]。这就是我认为存在某些差异的其他机制的原因。（请参见帖子顶部的编辑）

- sokin

---

嘿，非常感谢！编辑：算了，我自己解决了，再次感谢。 - TheProgramMAN123

不用谢。我刚刚编辑了我的答案，加入了转置。 - sokin