如何最小化此求和中的数值误差？

How can I minimize the numerical error in this summation?

提问人：0xbadf00d 提问时间：4/29/2023 最后编辑：0xbadf00d 更新时间：4/30/2023 访问量：193

问：

我在下面的代码示例中遇到了一个数字错误（我在下面添加了 Kahan 求和尝试和一个更聪明的幼稚但仍然幼稚的版本;但不幸的是，情况更糟）：

#include <algorithm>
#include <iostream>
#include <random>
#include <vector>


int main()
{
    std::random_device rd;
    std::mt19937 g{ rd() };
    std::uniform_real_distribution<> u;

    static std::size_t constexpr n = 1000;

    std::vector<double> q(n);
    std::generate_n(q.begin(), q.size(), [&]() { return u(g); });

    double average_of_q{};
    for (auto const& q : q)
        average_of_q += q;
    average_of_q /= n;

    std::vector<double> f(n);
    std::generate_n(f.begin(), n, [&]() { return u(g); });

    double sum1{};
    for (std::size_t i = 0; i < n; ++i)
        sum1 += std::abs(f[i] - q[i]);
    sum1 /= n;

    {
        double sum2{};
        for (std::size_t i = 0; i < n; ++i)
            sum2 += std::abs(f[i] - q[i]) - q[i];
        sum2 = sum2 / n + average_of_q;

        std::cout << "naive: " << std::abs(sum1 - sum2) << std::endl;
    }
    {
        double sum2{},
            c{};
        for (std::size_t i = 0; i < n; ++i)
        {
            double const x = std::abs(f[i] - q[i]) - q[i] - c,
                s = sum2 + x;
            c = (s - sum2) - x;

            sum2 = s;
        }
        sum2 = sum2 / n + average_of_q;

        std::cout << "kahan: " << std::abs(sum1 - sum2) << std::endl;
    }
    {
        double sum2{};
        for (std::size_t i = 0; i < n; ++i)
        {
            if (f[i] - q[i] >= 0)
                sum2 += f[i] - 2 * q[i];
            else
                sum2 -= f[i];
        }
        sum2 = sum2 / n + average_of_q;

        std::cout << "more clever, but still naive: " << std::abs(sum1 - sum2) << std::endl;
    }

    return 0;
}

输出是 ，而我们理论上期望它应该是 。如何优化此代码以使其尽可能小？1.11022e-160std::abs(sum1 - sum2)

为了激励这一点：在我的实际应用程序中，我已经知道并且我不需要遍历每个 ，因为我知道对于大多数 来说，这是非常小的，这就是为什么我需要将公式用于 。average_of_qistd::abs(f[i] - q[i])isum2

编辑：我也在MSE上询问了这个问题的理论部分（但略有不同;我不想在这里把事情搞得太复杂）：https://math.stackexchange.com/q/4688917/47771。

编辑2：我还试图通过将它们乘以一个因子来“提升”总和的项：

{
    double sum2{};
    for (std::size_t i = 0; i < n; ++i)
        sum2 += 1000 * (std::abs(f[i] - q[i]) - q[i]);
    sum2 = sum2 / (1000 * n) + average_of_q;

    std::cout << "boosted: " << std::abs(sum1 - sum2) << std::endl;
}

编辑3：

这可能是一个有用的信息：在我的实际应用中，与 .为简单起见，您可以假设所有 ，其中许多都在 1e-10 左右，但有一些接近 。f[i]q[i]q[i] = 1f[i]1

C++ C++20 精度 数值 方法

评论

0赞 Samuel Liew 4/30/2023
评论已移至聊天室;请不要在这里继续讨论。在在此下方发表评论之前，请查看评论的目的。不要求澄清或提出改进建议的评论通常属于 Meta Stack Overflow 或 Stack Overflow Chat 中的答案。继续讨论的评论可能会被删除。

答：

1赞 njuffa 4/30/2023 #1

这种情况基本上可以归结为比较计算相同数量的多种数学等效方法之间的累积舍入误差。正如在评论中指出的那样，为了尽量减少使用有限精度浮点计算时的数值差异，中间计算必须以高于目标精度的方式进行。具体而言，在此代码中，这需要应用于 、 和 的计算。average_of_qsum1sum2

此处的目标精度为 ，很可能映射到 IEEE-754 二进制浮点格式。各种编译器提供了某种形式的四精度浮点类型，这些类型可能映射到 IEEE-754 的 .例如，使用英特尔编译器 （） 时，将提供该类型，并且适用于此代码。然而，更便携的解决方案可以使用 Kahan 求和来累积准四倍精度。下面演示了这一点。doublebinary64binary128icl_Quad

从软件的角度来看，重要的是指示编译器不要重新关联浮点表达式，以便保留 Kahan 求和的数值属性。要强制执行的命令行标志因编译器而异，在我的情况下是 ，并且以下代码的整个编译器调用是/fp:stricticl /W4 /Ox /QxHOST /fp:strict array_sum_issue.cpp

有区别 ||通常为 0，但偶尔为 2-54 对应双精度 epsilon。n=1000000sum1 - sum2

#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <random>
#include <vector>

int main()
{
    std::random_device rd;
    std::mt19937 g{ rd() };
    std::uniform_real_distribution<> u;

    static std::size_t constexpr n = 1000000;

    std::vector<double> q(n);
    std::generate_n(q.begin(), q.size(), [&]() { return u(g); });

    double average_of_q{};
    {
        double sum = 0, c = 0;
        for (std::size_t i = 0; i < n; ++i) {
            double y = q[i] - c;
            double t = sum + y;
            c = (t - sum) - y;
            sum = t;
        }
        average_of_q = sum / n;
    }

    std::vector<double> f(n);
    std::generate_n(f.begin(), n, [&]() { return u(g); });

    double sum1{};
    {
        double sum = 0, c = 0;
        for (std::size_t i = 0; i < n; ++i) {
            double y = std::abs(f[i] - q[i]) - c;
            double t = sum + y;
            c = (t - sum) - y;
            sum = t;
        }
        sum1 = sum / n;
    }

    double sum2{};
    {
        double sum = 0, c = 0;
        for (std::size_t i = 0; i < n; ++i) {
            double y = (std::abs(f[i] - q[i]) - q[i]) - c;
            double t = sum + y;
            c = (t - sum) - y;
            sum = t;
        }
        sum2 = std::fma (sum, 1.0 / n, average_of_q);
    }
    double diff = std::abs (sum1 - sum2);
    printf ("average_of_q = % 23.16e (% 23.13a)\n", average_of_q, average_of_q);
    printf ("sum1         = % 23.16e (% 23.13a)\n", sum1, sum1);
    printf ("sum2         = % 23.16e (% 23.13a)\n", sum2, sum2);
    printf ("|sum1-sum2|  = % 23.16e (% 23.13a)\n", diff, diff);

    return EXIT_SUCCESS;
}

上述程序的示例输出（根据生成的随机数，数值会有所不同）：

average_of_q =  5.0031728235599426e-01 (   0x1.0029963aaf686p-1)
sum1         =  3.3347262871877092e-01 (   0x1.5579d949d5452p-2)
sum2         =  3.3347262871877092e-01 (   0x1.5579d949d5452p-2)
|sum1-sum2|  =  0.0000000000000000e+00 (   0x0.0000000000000p+0)

评论

0赞 0xbadf00d 4/30/2023
谢谢你的回答。当我替换 时，则 的输出是 。为什么？std::mt19937 g{ rd() };std::mt19937 g{ 0 };diff5.55112e-17
1赞 njuffa 4/30/2023
在通过有限精度浮点运算实现时，在数学上返回零的计算通常不会这样做。不同的输入序列将导致具有不同数量或误差的不同输出。使用卡汉求和可以减少累积误差，它不能完全消除它并增加新的误差。5.55e-17 表示 2**（-54）。IEEE-754 双精度格式有 53 个有效位（其中 52 个被存储，1 个隐式）。我们正处于精度极限。您应该看到右侧十六进制浮点数输出的最低有效数字相差 1。sum/n
0赞 njuffa 4/30/2023
@0xbadf00d 强烈推荐阅读：Nicholas J. Higham，Accuracy and Stability of Numerical Algorithms，第 2 版，SIAM 2002
0赞 0xbadf00d 5/3/2023
当我真的需要使用 Kahan 的算法并且普通求和就足够时，您能给出建议吗？另外：在什么情况下我需要使用？我们是否使用 or 是否重要（代数显然不是，但我可以想象当 or 为负或整数时，一个优先于另一个）。std::fmastd::fma(a, b, c)std::fma(b, a, c)ab
1赞 njuffa 5/3/2023
@0xbadf00d 在以下情况下使用 Kahan 求和：（1） 您非常关心在数学上应该为零的结果，但在浮点算术中计算时接近浮点 epsilon（这在实践中通常不是问题）或 （2） 简单求和和和 Kahan 求和之间存在显着的数值差异。.如果您的产品 a * b + c，其中 （a * b） 和 c 之间存在减法取消的风险，请使用 FMA。一本关于实用数值分析的好书会让你跟上速度;我已经推荐了一个。fma(a,b,c) == fma(b,a,c)

上一个：如何截断计算列中的小数点？

下一个：是否可以使此表达式更准确并且不返回负值？