Eigen 3x3 矩阵逆错误结果

如何解决Eigen 3x3 矩阵逆错误结果

问题描述

我将 Eigen 用于一些矩阵任务。假设我有矩阵 A 的大小是 4x3，那么它的转置 A^T 是 3x4 大小，那么 A^T * A 是 3x3 大小，因此逆矩阵 (A^T * A)^(-1) 也是 3x3 大小.

我想要(A^T * A)^(-1)。 通过使用上述公式，并通过手动定义 A^T * A 矩阵然后进行逆运算，我得到了 (A^T * A)^(-1) 矩阵的不同结果。很好奇为什么这两个结果不匹配并且相差很大。

重现

特征版本：git commit 972cf0c28a8d2ee0808c1277dea2c5c206591ce6

系统：Ubuntu 20.04

编译器：Clang++，10.0.0

代码：

#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.97211e+06,1404.32,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,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;
}

实际输出：

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

解决方法

这只是数值精度的问题。正如@Damien 在评论中指出的那样，矩阵是病态的，因此输入的微小差异会导致结果发生很大变化。通过仅从输出中复制前五位数字并使用它们手动定义第二个矩阵，可以丢弃在内部处理但未以 std::cout 的标准精度显示的重要部分。

以条目 ATA(0,0) 为例。通过使用 is ATA << 1.36154e+13,...，不考虑 1.e7 或更低阶的任何值。然而，这是矩阵中其他条目的顺序。

简而言之，两个结果都是正确的，但是您手动定义的矩阵 ATA 与第一部分中计算的不同。 （你可以拿两者的差异来验证这一点）