我怎样才能最小化这个求和中的数值误差？

Question

我在以下代码示例中遇到数字错误（我在下面添加了 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-16

，而我们理论上期望它应该是

0

。我怎样才能优化这段代码，使

std::abs(sum1 - sum2)

尽可能小？

为了激发这一点：在我的实际应用中，我已经知道

average_of_q

并且我不需要遍历每个

i

，因为我知道

std::abs(f[i] - q[i])

对于大多数

i

来说非常小，这是为什么我需要使用

sum2

.

的公式

编辑：我也在 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：

这可能是一个有用的信息：在我的实际应用中，许多

f[i]

与

q[i]

相比是小的。为简单起见，您可以假设所有

q[i] = 1

，许多

f[i]

都在1e-10左右，但少数接近

1

。

Answer 1

那应该如何工作？

for (auto const& q : q)
    average_of_q += q;