5

我一直在做一个矩阵类（作为一个学习练习），我在测试我的反函数时遇到了问题。

我这样输入一个任意矩阵：

2 1 1
1 2 1
1 1 2

并得到它来计算逆，我得到了正确的结果：

0.75 -0.25 -0.25
-0.25 0.75 -0.25
-0.25 -0.25 0.75

但是当我尝试将两者相乘以确保得到单位矩阵时，我得到：

1 5.5111512e-017 0
0 1 0
-1.11022302e-0.16 0 1

为什么我会得到这些结果？我会理解如果我将奇怪的数字相乘，我可以理解一些舍入错误，但它所做的总和是：

2 * -0.25 + 1 * 0.75 + 1 * -0.25

这显然是 0，而不是 5.111512e-017

如果我手动让它进行计算；例如：

std::cout << (2 * -0.25 + 1 * 0.75 + 1 * -0.25) << "\n";

我按预期得到0？

所有数字都表示为双精度数。这是我的乘法重载：

Matrix operator*(const Matrix& A, const Matrix& B)
{
    if(A.get_cols() == B.get_rows())
    {
        Matrix temp(A.get_rows(), B.get_cols());
        for(unsigned m = 0; m < temp.get_rows(); ++m)
        {
            for(unsigned n = 0; n < temp.get_cols(); ++n)
            {
                for(unsigned i = 0; i < temp.get_cols(); ++i)
                {
                    temp(m, n) += A(m, i) * B(i, n);
                }
            }
        }

        return temp;
    }

    throw std::runtime_error("Bad Matrix Multiplication");
}

和访问功能：

double& Matrix::operator()(unsigned r, unsigned c)
{
    return data[cols * r + c];
}

double Matrix::operator()(unsigned r, unsigned c) const
{
    return data[cols * r + c];
}

下面是求逆函数：

Matrix Inverse(Matrix& M)
{
    if(M.rows != M.cols)
    {
        throw std::runtime_error("Matrix is not square");
    }

    int r = 0;
    int c = 0;
    Matrix augment(M.rows, M.cols*2);
    augment.copy(M);

    for(r = 0; r < M.rows; ++r)
    {
        for(c = M.cols; c < M.cols * 2; ++c)
        {
            augment(r, c) = (r == (c - M.cols) ? 1.0 : 0.0);
        }
    }

    for(int R = 0; R < augment.rows; ++R)
    {
        double n = augment(R, R);
        for(c = 0; c < augment.cols; ++c)
        {
            augment(R, c) /= n;
        }

        for(r = 0; r < augment.rows; ++r)
        {
            if(r == R) { continue; }
            double a = augment(r, R);

            for(c = 0; c < augment.cols; ++c)
            {
                augment(r, c) -= a * augment(R, c);
            }
        }
    }

    Matrix inverse(M.rows, M.cols);
    for(r = 0; r < M.rows; ++r)
    {
        for(c = M.cols; c < M.cols * 2; ++c)
        {
            inverse(r, c - M.cols) = augment(r, c);
        }
    }

    return inverse;
}
4

5 回答 5

11

请阅读这篇论文：每个计算机科学家都应该知道的关于浮点运算的知识

于 2011-06-03T18:04:51.627 回答
9

你的倒置矩阵中有像 0.250000000000000005 这样的数字，它们只是为了显示而四舍五入，所以你会看到像 0.25 这样的小整数。

于 2011-06-03T18:08:51.810 回答
2

你总是会遇到这样的浮点舍入错误，尤其是在处理没有精确二进制表示的数字时（即，你的数字不等于 2^(N) 或 1/(2^N)，其中 N 是某个整数值）。

话虽如此，有多种方法可以提高结果的精度，您可能需要在谷歌上搜索使用固定精度浮点值的数值稳定的高斯消除算法。

如果您愿意加快速度，您也可以合并一个使用有理数的无限精度数学库，如果您选择，请避免使用会产生无理数的根。有许多库可以帮助您使用有理数，例如GMP。您也可以rational自己创建一个类，但请注意，如果您只使用无符号 64 位值以及有理数分量的额外符号标志变量，则相对容易溢出多个数学运算的结果。这就是 GMP 及其无限长度的整数字符串对象派上用场的地方。

于 2011-06-03T18:29:55.707 回答
2

你不应该对这些数字有任何问题，因为对于这个特定的矩阵，逆是 2 的所有幂，并且可以准确地表示。一般来说，浮点数运算会引入可能累积的小错误，结果可能令人惊讶。

在你的情况下，我很确定逆是不准确的，你只是显示前几位数字。即，它不完全是 0.25 (=1/4)、0.75 (=3/4) 等。

于 2011-06-03T18:10:05.223 回答
1

这只是简单的浮点错误。即使double是计算机上的 s 也不是 100% 准确的。没有办法 100% 准确地表示具有有限位数的二进制以 10 为底的十进制数。

于 2011-06-03T18:04:40.100 回答