用微小的数字代替零

Tiny numbers in place of zero?

本文关键字:数字      更新时间:2023-10-16

我一直在制作一个矩阵类(作为学习练习),在测试逆函数时遇到了问题。

我输入一个任意矩阵,如下所示:

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

请阅读这篇论文:每一位计算机科学家都应该知道浮点算术

你的倒矩阵中有0.250000000000000005这样的数字,它们只是四舍五入以显示,所以你可以看到0.25这样的漂亮的小整数。

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

在你的情况下,我确信倒数是不准确的,你只是显示了前几个数字。也就是说,它不是精确地0.25(=1/4)、0.75(=3/4)等。

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

也就是说,有很多方法可以提高结果的精度,你可能想在谷歌上搜索使用固定精度浮点值的数值稳定的高斯消除算法。

如果你愿意快速命中,你也可以加入一个使用有理数的无限精度数学库,如果你选择了这个选项,只需避免使用会产生无理数的根。有很多库可以帮助你使用有理数,比如GMP。您也可以自己创建rational类,不过要注意,如果您只使用无符号64位值和一个额外的符号标志变量作为有理数的组成部分,则相对容易溢出多个数学运算的结果。这就是GMP的用武之地,它有无限长度的整数字符串对象。

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