c++ - 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 矩阵的逆矩阵，我会得到与特征值不同的结果。

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

谢谢!

Answer 1

编辑: 进一步深究为什么最终结果会有一些差异后，我发现案例2中提到的那些“底层机制”实际上是我在输入矩阵值时由于不慎而导致的错误。

以下是处理了这些错误的原始答案。

---

实际问题不在于 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

您可以像这样转置复制的最终需求 vector :

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 的行列式

案例2

即使用反向的最终需求 vector

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 检查 仍然没有产生最初定义的 Finald 的准确值，但我相信这与此有关具有精度或其他一些潜在机制。 (见注意)

作为概念证明，这里是使用 MATLAB 获得的一些结果，使用第二种情况(相反)作为复制最终需求 vector (denom):

>> 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 相同，但是，如果您执行最后一个矩阵乘法(在Final Demand vector 的验证)手动，您会发现实际上 IminAinv 的 MATLAB 和案例 2 版本产生与案例 2 的最终输出相同的结果, 即 [88.9219, 125.728, 59.5037, 105.543, 84.5808]。 这就是为什么我认为这些差异还涉及一些其他机制。(请参阅帖子顶部的编辑)