Eigen中逆矩阵的计算出错

Question

我正在尝试构建一个简单的 input/output 矩阵（如果需求增加，您可以在其中计算简单经济中的乘数效应）。但由于某种原因，最终结果并不相符。

#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中提到的那些"underlying mechanisms"其实是我自己在输入矩阵值时不小心造成的错误。

下面是原始答案，其中包含这些错误。

实际问题不在于 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 检查 仍然没有产生最初定义的 Finald 的准确值，但我认为这与精度或其他一些潜在机制有关。（见注）~~

作为概念证明，这里是使用 MATLAB 获得的一些结果，对 replicated Final Demand 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 相同，但是，如果您手动执行最后一个矩阵乘法（在最终需求向量验证中的乘法），您将看到实际上 IminAinv 的 MATLAB 和案例 2 版本 产生与案例 2 的最终输出相同的结果，即 [88.9219, 125.728, 59.5037, 105.543, 84.5808]。这就是为什么我认为这些差异还涉及一些其他机制。（请参阅 post 之上的编辑）

Eigen中逆矩阵的计算出错

Calculation of inverse matrix in Eigen going wrong

c++

linear-algebra

eigen

案例一

案例二