Eigen 3x3 matrix inverse 错误结果

Question

问题描述

我正在使用 Eigen 执行一些矩阵任务。假设我有大小为 4x3 的矩阵 A，那么它的转置 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

系统：Ubuntu20.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,
             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;
}

实际输出：

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

Answer 1

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

标准精度显示的重要部分

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

总之，两个结果都是正确的，但是你手动定义的矩阵ATA和第一部分计算的不一样。（你可以拿两者的区别来验证一下）

Answer 2

我今天对逆矩阵 3x3 的练习是使用

template<typename Derived>
template<typename ResultType>
inline void MatrixBase<Derived>::computeInverseAndDetWithCheck(
    ResultType& inverse,
    typename ResultType::Scalar& determinant,
    bool& invertible,
    const RealScalar& absDeterminantThreshold
  ) const

设置手动 absDeterminantThreshold 为默认值： 1) 对于浮点数 1e-5 2) 对于双 1e-12 3）长双1e-15 例如您的代码：

Eigen::Matrix3d ATA_inv;
bool invertible = false;
double determinant = 0.;
double Threshold = abs(ATA.diagonal().prod());//maybe any test
ATA.computeInverseAndDetWithCheck(ATA_inv,determinant,invertible,Threshold<1e-12?Threshold*0.1:1e-12);

Eigen 3x3 matrix inverse 错误结果

Eigen 3x3 matrix inverse wrong result

c++

matrix

eigen

问题描述

复制