Calculation of inverse matrix in Eigen going wrong

Question

Im trying to build a simple input/output matrix (where you can calculate the multiplier effect in a simple economy if demand increases). But for some reason the final result is not adding up.

#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;

When I verify if the (I-A)inverse matrix (or IminAinv) is properly calculated it doesn't add up. By multiplying the IminAinv by Internal demand (int), I should get the same Intd. That is if Intd isn't changed. Instead I get a bigger number. Also if I calculate the inverse of the IminA matrix myself, I get something different then eigen.

So something goes wrong in getting the inverse of Identity matrix - Coefficient matrix. But what?

Thanks!

Answer 1

EDIT : After some more digging into why there are some differences in the final results, I discovered that those "underlying mechanisms" mentioned in Case 2 were, in fact, my own mistake out of inadvertence while entering the matrices' values.

Here follows the original answer with these mistakes taken care of.

---

The actual problem does not lie in the inversion of the AMatrix but in a much more subtle detail. You are using this command to perform the division in the definition of AMatrix:

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

But if you check the results of this replicate operation on Finald you get:

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

whereas the correct one should be:

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

You can transpose the replicated final demand vector like this:

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

and then of course:

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

which yields:

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

So, let's see what are the differences in your intermediate and final results in both these cases:

Case 1

ie Your current approach

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

I have also added the determinant of IminA

Case 2

ie using the reversed final demand 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

Now, I understand that the Finald check still does not produce the exact values of the originally defined Finald, but I believe that this has something to do with the precision or some other underlying mechanism. (see NOTE)

As a proof of concept here are some results obtained with MATLAB, using the second case (reversed) for the 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

It is quite clear that the second case results are (almost) identical to the MATLAB ones.

NOTE: Here you can see that the MATLAB output is identical to the original Finald, however, if you perform the last matrix multiplication (the one in the validation of the Final Demand vector) by hand, you will see that in fact both MATLAB and Case 2 versions of IminAinv yield the same result as the final output of Case 2, ie [88.9219, 125.728, 59.5037, 105.543, 84.5808]. This is why I think that there is some other mechanism involved in these differences. (See EDIT on top of post)