I'm using Eigen for some matrix task. Say I have matrix A whose size is 4x3, then its transpose A^T is 3x4 size, then A^T * A is 3x3 size, thus the inverse, (A^T * A)^(-1), is also 3x3 size.
I would like to get (A^T * A)^(-1). By using the mentioned formula, and by manually defining A^T * A matrix then do inverse, I got different result for (A^T * A)^(-1) matrix. Curious why this two results mismatch and differs very much.
Eigen version: git commit 972cf0c28a8d2ee0808c1277dea2c5c206591ce6
System: Ubuntu 20.04
Compiler: Clang++, 10.0.0
Code:
#include <iostream>
#include <Eigen/Dense>
int main() {
{
printf("----- A^T * A by calculate\n");
Eigen::MatrixXd A(4, 3);
A << 1.70023e+06, 1303.93, 1,
1.99113e+06, 1411.07, 1,
1.97211e+06, 1404.32, 1,
1.69433e+06, 1301.67, 1;
std::cout << "A:" << std::endl << A << std::endl;
Eigen::MatrixXd AT = A.transpose();
std::cout << "A^T:" << std::endl << AT << std::endl;
Eigen::MatrixXd ATA = AT * A;
std::cout << "A^T * A:" << std::endl << ATA << std::endl;
std::cout << "(A^T * A)^(-1):" << std::endl << ATA.inverse() << std::endl;
}
{
printf("----- A^T * A by define\n");
Eigen::MatrixXd ATA(3, 3);
ATA << 1.36154e+13, 1.00015e+10, 7.3578e+06,
1.00015e+10, 7.3578e+06, 5420.99,
7.3578e+06, 5420.99, 4;
std::cout << "A^T * A:" << std::endl << ATA << std::endl;
std::cout << "(A^T * A)^(-1):" << std::endl << ATA.inverse() << std::endl;
}
return 0;
}
The actual outputs:
----- A^T * A by calculate
A:
1.70023e+06 1303.93 1
1.99113e+06 1411.07 1
1.97211e+06 1404.32 1
1.69433e+06 1301.67 1
A^T:
1.70023e+06 1.99113e+06 1.97211e+06 1.69433e+06
1303.93 1411.07 1404.32 1301.67
1 1 1 1
A^T * A:
1.36154e+13 1.00015e+10 7.3578e+06
1.00015e+10 7.35781e+06 5420.99
7.3578e+06 5420.99 4
(A^T * A)^(-1):
3.49282e-06 -0.00946871 6.40758
-0.00946871 25.6689 -17370.5
6.40758 -17370.5 1.17549e+07
----- A^T * A by define
A^T * A:
1.36154e+13 1.00015e+10 7.3578e+06
1.00015e+10 7.3578e+06 5420.99
7.3578e+06 5420.99 4
(A^T * A)^(-1):
6.1435e-09 -1.66513e-05 0.0112659
-1.66513e-05 0.0452221 -30.658
0.0112659 -30.658 20826.4
This is just a matter of numerical accuracy. As pointed out by @Damien in the comment, the matrix is ill-conditioned, and thus a small difference in the input can lead to a large change in the results. By copying from the output only the first five digits and using them to manually define the second matrix, a significant part is discarded that is handled internally but not displayed with the standard accuracy of std::cout.
Take for instance the entry ATA(0,0). By using is ATA << 1.36154e+13,..., any value of the order of 1.e7 or lower is not considered. This is, however, the order of the other entries in the matrix.
In short, both results are correct, but your manually defined matrix ATA is not the same as the one that is calculated in the first part. (You can take the difference between the two to verify this).