使用高斯消除计算矩阵的逆时非常不准确答案

【问题标题】：Very high inaccuray when calculating inverse of matrix using gauss elimination使用高斯消除计算矩阵的逆时非常不准确
【发布时间】：2017-09-25 02:11:42
【问题描述】：

我现在正在开发一个 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 行是反函数所在的位置。

感谢任何帮助。

【问题讨论】：

有些矩阵没有逆矩阵。就像第一个 3x3 示例一样。对他们来说，程序会产生无意义的结果。
确保矩阵的行列式非零（即矩阵非奇异且可以反转）。即使行列式它不是零但非常接近于零，你也必须小心（阅读条件数），参见例如math.stackexchange.com/questions/1622610/…
@vsoftco 在这两种情况下，行列式都是非零的。
@anatolyg 我怎么知道是否有逆？
还有一个问题：使用这种方法对矩阵求逆时，行的对角元素（将用作除数）可能为零，在这种情况下，该行需要与该列中没有零的后面的行交换。需要更新代码以处理这种情况。

标签： c++ math matrix floating-point matrix-inverse

【解决方案1】：

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

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

更大的问题是，这不容易检测：如果由于舍入错误而没有错误，那将是一块蛋糕 - 只需测试除数不为 0！但是浮点运算存在舍入误差，因此除数将是一个非常小的非零数。

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

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

使用知名且经过良好测试的库进行此类计算可能是明智之举 - 需要考虑的因素很多，而且可能会出错。

【讨论】：