如何最小化求和运算的数值误差？C++代码优化求助

我在以下代码中遇到了数值误差问题：尝试了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，但理论上std::abs(sum1 - sum2)应该为0。如何优化代码，让这个差值尽可能小？

背景说明 #

实际应用中已获知average_of_q，且无需遍历所有i——因为大部分std::abs(f[i] - q[i])的值极小，因此必须使用sum2的计算公式。

已尝试的优化 #

• 放大求和项：

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

补充场景信息 #

实际应用中，多数f[i]远小于q[i]。可简化为：所有q[i] = 1，多数f[i]约为1e-10，少数接近1。

优化方案 #

1. 利用场景特性重构计算逻辑 #

结合q[i]=1、多数f[i]<<1的场景，拆分std::abs(f[i]-q[i])的计算：

• 当f[i] <= q[i]（绝大多数情况），std::abs(f[i]-q[i]) = q[i] - f[i]，代入求和项得：std::abs(f[i]-q[i]) - q[i] = -f[i]
• 当f[i] > q[i]（少数情况），求和项为：(f[i]-q[i]) - q[i] = f[i] - 2q[i]

这种拆分避免了q[i]的抵消操作，减少浮点数相减带来的精度损失，代码示例：

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

2. 使用更高精度的累加类型 #

将sum2的类型从double改为long double，利用更高精度存储累加中间值，最后再转换回double计算结果，减少中间过程的精度丢失：

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

3. 分组求和+误差补偿 #

对于大规模数据，将求和项分成若干小组，每组内用Kahan求和法，再将各组结果相加。这种方法既减少了大数“吞噬”小数的问题，又结合了Kahan的误差补偿机制：

{
    constexpr size_t group_size = 100;
    std::vector<double> group_sums(n / group_size + 1, 0.0);
    std::vector<double> group_corrections(n / group_size + 1, 0.0);

    for (size_t i = 0; i < n; ++i)
    {
        size_t group_idx = i / group_size;
        double x = f[i] > q[i] ? (f[i] - 2 * q[i]) : (-f[i]);

        double temp = x - group_corrections[group_idx];
        double new_sum = group_sums[group_idx] + temp;
        group_corrections[group_idx] = (new_sum - group_sums[group_idx]) - temp;
        group_sums[group_idx] = new_sum;
    }

    double sum2 = 0.0;
    double c = 0.0;
    for (double s : group_sums)
    {
        double temp = s - c;
        double new_sum = sum2 + temp;
        c = (new_sum - sum2) - temp;
        sum2 = new_sum;
    }

    sum2 = sum2 / n + average_of_q;
    std::cout << "grouped kahan: " << std::abs(sum1 - sum2) << std::endl;
}

4. 优化average_of_q的精度匹配 #

如果average_of_q是预先计算的，用long double类型存储它，或者在计算sum2/n + average_of_q时，先将sum2/n转换为long double再与average_of_q相加，避免低精度下的加法误差。

内容的提问来源于stack exchange，提问作者0xbadf00d