本征逆矩阵计算出错

Calculation of inverse matrix in Eigen going wrong

本文关键字：计算 出错      更新时间：2023-10-16

我试图建立一个简单的输入/输出矩阵(如果需求增加，您可以在其中计算简单经济中的乘数效应)。但由于某种原因，最终的结果并没有加起来。

#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矩阵的逆，我得到的是与本征不同的东西。

所以在求单位矩阵-系数矩阵的逆时出了问题。但是什么?

谢谢!

EDIT:在深入研究了为什么最终结果会有一些差异之后，我发现案例2中提到的那些"潜在机制"实际上是我在输入矩阵值时的疏忽造成的错误。

以下是原文的答案，并纠正了这些错误。

实际问题不在于矩阵的反转，而在于更微妙的细节。使用此命令执行矩阵:

定义中的除法。
AMatrix = ProdA.array() / (Finald.replicate(1, ProdA.cols())).array();

但是如果你在 final 上检查这个复制操作的结果，你会得到:

...
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

那么，让我们看看在这两种情况下，中间结果和最终结果有什么不同:

案例1

ie当前方法

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

即使用反向的最终需求向量

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获得的一些结果，使用第二种情况(相反)对复制最终需求向量(dom):

>> 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相同，但是，如果您手工执行最后一次矩阵乘法(最终需求向量验证中的一个)，您将看到实际上MATLAB和Case 2版本的IminAinv产生与Case 2的最终输出相同的结果，即[88.9219,125.728,59.5037,105.543,84.5808]。这就是为什么我认为有一些其他机制涉及到这些差异。

相关文章：