使用高斯消元法计算矩阵的逆时非常不准确

Question

我现在正在开发一个 c++ 代码库，它使用矩阵库来计算各种东西。其中之一就是计算矩阵的逆。它使用高斯消除来实现这一点。但结果非常不准确。如此之多以至于将逆矩阵与原始矩阵相乘甚至不能接近单位矩阵。

下面是用于计算逆的代码，矩阵以数字类型和行和列为模板：

/// \brief Take the inverse of the matrix.
/// \return A new matrix which is the inverse of the current one.
matrix<T, M, M> inverse() const
{
    static_assert(M == N, "Inverse matrix is only defined for square matrices.");

    // augmented the current matrix with the identiy matrix.
    auto augmented = this->augment(matrix<T, M, M>::get_identity());

    for (std::size_t i = 0; i < M; i++)
    {
        // divide the current row by the diagonal element.
        auto divisor = augmented[i][i];
        for (std::size_t j = 0; j < 2 * M; j++)
        {
            augmented[i][j] /= divisor;
        }

        // For each element in the column of the diagonal element that is currently selected
        // set all element in that column to 0 except the diagonal element by using the currently selected row diagonal element.

        for (std::size_t j = 0; j < M; j++)
        {
            if (i == j)
            {
                continue;
            }

            auto multiplier = augmented[j][i];
            for (std::size_t k = 0; k < 2 * M; k++)
            {
                augmented[j][k] -= multiplier * augmented[i][k];
            }
        }
    }

    // Slice of the the new identity matrix on the left side.
    return augmented.template slice<0, M, M, M>();
}

现在我做了一个单元测试，使用预先计算的值测试逆运算是否正确。我尝试了两个矩阵，一个 3x3 和一个 4x4。我用这个网站来计算倒数：https://matrix.reshish.com/，它们确实在一定程度上匹配。因为单元测试确实成功了。但是一旦我计算出原始矩阵 * 逆矩阵，甚至没有实现类似于单位矩阵的东西。请参阅下面代码中的注释。

BOOST_AUTO_TEST_CASE(matrix_inverse)
{
    auto m1 = matrix<double, 3, 3>({
        {7, 8, 9},
        {10, 11, 12},
        {13, 14, 15}
    });

    auto inverse_result1 = matrix<double,3, 3>({
        {264917625139441.28, -529835250278885.3, 264917625139443.47},
        {-529835250278883.75, 1059670500557768, -529835250278884.1},
        {264917625139442.4, -529835250278882.94, 264917625139440.94}
    });

    auto m2 = matrix<double, 4, 4>({
        {7, 8, 9, 23},
        {10, 11, 12, 81},
        {13, 14, 15, 11},
        {1, 73, 42, 65}
    });

    auto inverse_result2 = matrix<double, 4, 4>({
        {-0.928094660194201, 0.21541262135922956, 0.4117111650485529, -0.009708737864078209},
        {-0.9641231796116679, 0.20979975728155775, 0.3562651699029188, 0.019417475728154842},
        {1.7099261731391882, -0.39396237864078376, -0.6169346682848 , -0.009708737864076772 },
        {-0.007812499999999244, 0.01562499999999983, -0.007812500000000278, 0}
    });


    // std::cout << (m1.inverse() * m1) << std::endl;
    // results in
    // 0.500000000     1.000000000     -0.500000000
    // 1.000000000     0.000000000     0.500000000
    // 0.500000000     -1.000000000    1.000000000
    // std::cout << (m2.inverse() * m2) << std::endl; 
    // results in
    // 0.396541262     -0.646237864    -0.689016990    -2.162317961
    // 1.206917476     2.292475728     1.378033981     3.324635922
    // -0.884708738    -0.958737864    -0.032766990    -3.756067961
    // -0.000000000    -0.000000000    -0.000000000    1.000000000

    BOOST_REQUIRE_MESSAGE(
        m1.inverse().fuzzy_equal(inverse_result1, 0.1) == true,
        "3x3 inverse is not the expected result."
    );

    BOOST_REQUIRE_MESSAGE(
        m2.inverse().fuzzy_equal(inverse_result2, 0.1) == true,
        "4x4 inverse is not the expected result."
    );
}

我束手无策。我绝不是矩阵数学方面的专家，因为我必须在工作中学习它，但这真的让我很困惑。

完整的代码矩阵 class 可在以下位置获得： https://codeshare.io/johnsmith

404行是反函数所在的地方

感谢任何帮助。

Answer 1

正如评论中已经确定的那样，感兴趣的矩阵是奇异的，因此没有逆矩阵。

很好，您的测试已经发现代码中的第一个问题 - 这种情况处理不当，没有引发错误。

更大的问题是，这不容易检测到：如果没有由于舍入错误而导致的错误，那将是小菜一碟 - 只需测试除数不为 0 即可！但是浮点运算存在舍入误差，所以除数会是一个非常小的非零数。

并且无法判断这个非零值是由于舍入误差还是由于矩阵接近奇异（但不是奇异）这一事实。但是，如果矩阵接近奇异，则它的条件很差，因此无论如何结果都不可信。

因此，理想情况下，算法不仅应该计算逆矩阵，还应该（估计）原始矩阵的条件，以便调用者可以对糟糕的条件做出反应。

对于这种计算，使用众所周知且经过充分测试的库可能是明智的 - 有很多需要考虑的地方，也有可能做错的地方。

使用高斯消元法计算矩阵的逆时非常不准确

Very high inaccuray when calculating inverse of matrix using gauss elimination

c++

math

floating-point

matrix

matrix-inverse