使用高斯消除计算矩阵逆时非常高的误差

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."
    );
}

我已经黔驴技穷了。我绝不是矩阵数学的专家，因为我不得不在工作中学习，但这真的难倒我了。

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

第 404 行是反函数所在的位置。

感谢任何帮助。

Answer 1

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

太好了，您的测试已经发现了代码中的第一个问题 - 这种情况没有得到正确处理，没有引发错误。

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

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

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

对于这种计算，使用众所周知且经过良好测试的库可能是明智的——需要考虑很多问题，以及哪些地方可能会出错。