应该按什么顺序添加浮点数以获得最精确的结果？

Question

这是我在最近的采访中被问到的一个问题，我想知道（我实际上不记得数值分析的理论，所以请帮助我:)

如果我们有一些函数，它会累积浮点数：

std::accumulate(v.begin(), v.end(), 0.0);

例如，v 是 std::vector<float>。

• 在累积这些数字之前对它们进行排序会更好吗？

• 哪个顺序会给出最准确的答案？

我怀疑将数字按升序排序实际上会使数字错误少，但不幸的是我自己无法证明。

附注我确实意识到这可能与现实世界的编程无关，只是好奇。

Answer 1

您的直觉基本上是正确的，按升序（数量级）排序通常会有所改善。考虑我们添加单精度（32 位）浮点数的情况，有 10 亿个值等于 1 /（10 亿），一个值等于 1。如果 1 先出现，那么总和会出现为 1，因为 1 + (1 / 10 亿) 由于精度损失而为 1。每次添加对总数都没有影响。

如果小值首先出现，它们至少会相加，尽管即使那样我也有 2^30 个，而在 2^25 左右之后，我又回到了每个单独不是的情况再影响总量。所以我仍然需要更多技巧。

这是一种极端情况，但通常添加两个幅度相似的值比添加两个幅度非常不同的值更准确，因为您以这种方式“丢弃”较小值的更少精度位。通过对数字进行排序，您可以将大小相似的值组合在一起，并通过将它们按升序添加，您可以为较小的值提供累积达到较大数字大小的“机会”。

不过，如果涉及负数，则很容易“智取”这种方法。考虑三个值的总和，{1, -1, 1 billionth}。算术上正确的总和是1 billionth，但如果我的第一次加法涉及微小的值，那么我的最终总和将为 0。在 6 个可能的订单中，只有 2 个是“正确的” - {1, -1, 1 billionth} 和 {-1, 1, 1 billionth}。所有 6 个阶数给出的结果在输入中的最大数值范围内是准确的（0.0000001% 出），但其中 4 个的结果在真实解的尺度上是不准确的（100% 出）。您正在解决的特定问题会告诉您前者是否足够好。

事实上，除了按排序顺序添加它们之外，您还可以玩更多的技巧。如果您有很多非常小的值、中等数量的中等值和少量的大值，那么首先将所有小值相加，然后分别合计中等值，再将这两个总计相加可能是最准确的一起然后添加大的。找到最准确的浮点加法组合并非易事，但要应对非常糟糕的情况，您可以将整个数组的运行总计保持在不同的大小，将每个新值添加到与其大小最匹配的总数中，当一个连续的总数开始变得太大时，将其添加到下一个总数中并开始一个新的总数。从逻辑上讲，这个过程相当于以任意精度类型执行求和（所以你会这样做）。但考虑到按数量级升序或降序添加的简单选择，升序是更好的选择。

它确实与现实世界的编程有一定的关系，因为在某些情况下，如果你不小心砍掉了由大量值组成的“重”尾，每个值都太小了，你的计算可能会出现非常严重的错误单独影响总和，或者如果您从许多仅影响总和的最后几位的小值中丢弃过多的精度。在无论如何尾巴可以忽略不计的情况下，您可能不在乎。例如，如果您一开始只将少量值相加，并且只使用总和的几个有效数字。

Answer 2

还有一种为这种累加运算设计的算法，称为Kahan Summation，您可能应该知道。

根据维基百科，

Kahan 求和算法（也称为补偿求和）显着降低了通过添加有限精度浮点数序列获得的总数中的数值误差，与显而易见的方法。这是通过保持单独的运行补偿（累积小错误的变量）来完成的。

在伪代码中，算法是：

function kahanSum(input)
 var sum = input[1]
 var c = 0.0          //A running compensation for lost low-order bits.
 for i = 2 to input.length
  y = input[i] - c    //So far, so good: c is zero.
  t = sum + y         //Alas, sum is big, y small, so low-order digits of y are lost.
  c = (t - sum) - y   //(t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
  sum = t             //Algebraically, c should always be zero. Beware eagerly optimising compilers!
 next i               //Next time around, the lost low part will be added to y in a fresh attempt.
return sum

Answer 3

我尝试了 Steve Jessop 提供的答案中的极端示例。

#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
    long billion = 1000000000;
    double big = 1.0;
    double small = 1e-9;
    double expected = 2.0;

    double sum = big;
    for (long i = 0; i < billion; ++i)
        sum += small;
    std::cout << std::scientific << std::setprecision(1) << big << " + " << billion << " * " << small << " = " <<
        std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    sum = 0;
    for (long i = 0; i < billion; ++i)
        sum += small;
    sum += big;
    std::cout  << std::scientific << std::setprecision(1) << billion << " * " << small << " + " << big << " = " <<
        std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    return 0;
}

我得到了以下结果：

1.0e+00 + 1000000000 * 1.0e-09 = 2.000000082740371    (difference = 0.000000082740371)
1000000000 * 1.0e-09 + 1.0e+00 = 1.999999992539933    (difference = 0.000000007460067)

第一行的错误是第二行的十倍以上。

如果我在上面的代码中将doubles 更改为floats，我会得到：

1.0e+00 + 1000000000 * 1.0e-09 = 1.000000000000000    (difference = 1.000000000000000)
1000000000 * 1.0e-09 + 1.0e+00 = 1.031250000000000    (difference = 0.968750000000000)

两个答案都没有接近 2.0（但第二个稍微接近）。

如 Daniel Pryden 所述，使用 Kahan 求和（doubles）：

#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
    long billion = 1000000000;
    double big = 1.0;
    double small = 1e-9;
    double expected = 2.0;

    double sum = big;
    double c = 0.0;
    for (long i = 0; i < billion; ++i) {
        double y = small - c;
        double t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }

    std::cout << "Kahan sum  = " << std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    return 0;
}

我得到的正好是 2.0：

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

即使我在上面的代码中将doubles 更改为floats，我也会得到：

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

看来Kahan是要走的路！

Answer 4

基于史蒂夫的答案，首先按升序对数字进行排序，我将介绍另外两个想法：

1. 确定两个数字的指数差异，如果超过这个差异，您可能会认为您会损失太多的精度。

2. 然后按顺序将数字相加，直到累加器的指数对于下一个数字来说太大，然后将累加器放到一个临时队列中，并以下一个数字启动累加器。继续，直到用完原始列表。

您使用临时队列（已对其进行排序）重复该过程，并且指数差异可能更大。

如果你必须一直计算指数，我认为这会很慢。

我快速运行了一个程序，结果是 1.99903

Answer 5

我认为你可以做得更好，而不是在累积之前对数字进行排序，因为在累积的过程中，累加器会越来越大。如果您有大量相似的数字，您将很快开始失去精度。这是我的建议：

while the list has multiple elements
    remove the two smallest elements from the list
    add them and put the result back in
the single element in the list is the result

当然，这个算法使用优先级队列而不是列表会最有效。 C++代码：

template <typename Queue>
void reduce(Queue& queue)
{
    typedef typename Queue::value_type vt;
    while (queue.size() > 1)
    {
        vt x = queue.top();
        queue.pop();
        vt y = queue.top();
        queue.pop();
        queue.push(x + y);
    }
}

司机：

#include <iterator>
#include <queue>

template <typename Iterator>
typename std::iterator_traits<Iterator>::value_type
reduce(Iterator begin, Iterator end)
{
    typedef typename std::iterator_traits<Iterator>::value_type vt;
    std::priority_queue<vt> positive_queue;
    positive_queue.push(0);
    std::priority_queue<vt> negative_queue;
    negative_queue.push(0);
    for (; begin != end; ++begin)
    {
        vt x = *begin;
        if (x < 0)
        {
            negative_queue.push(x);
        }
        else
        {
            positive_queue.push(-x);
        }
    }
    reduce(positive_queue);
    reduce(negative_queue);
    return negative_queue.top() - positive_queue.top();
}

队列中的数字是负数，因为top 产生最大 数字，但我们想要最小。我本可以为队列提供更多模板参数，但这种方法似乎更简单。

Answer 6

有一类算法可以解决这个确切的问题，无需对数据进行排序或重新排序。

换句话说，求和可以在数据上一次完成。这也使得此类算法适用于事先不知道数据集的情况，例如如果数据是实时到达的，需要维护运行总和。

这是最近一篇论文的摘要：

我们提出了一种新颖的在线算法，用于对流进行精确求和 的浮点数。 “在线”是指算法 一次只需要查看一个输入，并且可以取任意一个 此类输入的长度输入流，同时只需要常量 记忆。 “精确”是指我们的内部数组的总和 算法完全等于所有输入的总和，并且 返回的结果是正确舍入的总和。正确性证明 适用于所有输入（包括非规范化数字，但模 中间溢出），并且与被加数的数量无关 或求和的条件数。该算法渐近需要 每个加法运算只有 5 次 FLOP，并且由于指令级并行性 运行速度只比明显的、快速但愚蠢的慢 2--3 倍 当和数为时的“普通递归求和”循环 大于 10,000。因此，据我们所知，它是最快、最 在已知算法中准确且内存效率最高。的确，它 很难看出一个更快的算法或一个需要 如果没有硬件改进，可能存在的 FLOP 会显着减少。 提供大量求和的应用程序。

来源：Algorithm 908: Online Exact Summation of Floating-Point Streams。

Answer 7

这并不能完全回答您的问题，但明智的做法是计算两次总和，一次使用 rounding mode“向上取整”，一次使用“向下取整”。比较这两个答案，您就会知道/如何/您的结果不准确，因此您需要使用更聪明的求和策略。不幸的是，大多数语言并没有让浮点舍入模式变得如此简单，因为人们不知道它在日常计算中实际上是有用的。

看看Interval arithmetic，你在哪里做所有这样的数学运算，随时保持最高和最低值。它带来了一些有趣的结果和优化。

Answer 8

提高准确性的最简单排序是按绝对值升序排序。这让最小的震级值有机会在与可能导致精度损失的较大震级值交互之前累积或取消。

也就是说，您可以通过跟踪多个不重叠的部分和来做得更好。这是描述该技术并提供准确性证明的论文：www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

该算法和其他精确浮点求和的方法是在简单的 Python 中实现的：http://code.activestate.com/recipes/393090/ 其中至少两个可以简单地转换为 C++。

Answer 9

您的浮点数应该以双精度添加。这将为您提供比任何其他技术都更高的精度。为了获得更高的精度和更快的速度，您可以创建四个总和，并在最后将它们相加。

如果您要添加双精度数字，请使用 long double 作为总和 - 但是，这只会在 long double 实际上比 double 精度更高的实现中产生积极影响（通常是 x86，PowerPC 取决于编译器设置）。

Answer 10

对于 IEEE 754 单精度或双精度或已知格式的数字，另一种选择是使用由指数索引的数字数组（由调用者传递，或在 C++ 的类中）。将数字添加到数组中时，仅添加具有相同指数的数字（直到找到一个空槽并存储该数字）。当调用求和时，数组从最小到最大求和以最小化截断。单精度示例：

/* clear array */
void clearsum(float asum[256])
{
size_t i;
    for(i = 0; i < 256; i++)
        asum[i] = 0.f;
}

/* add a number into array */
void addtosum(float f, float asum[256])
{
size_t i;
    while(1){
        /* i = exponent of f */
        i = ((size_t)((*(unsigned int *)&f)>>23))&0xff;
        if(i == 0xff){          /* max exponent, could be overflow */
            asum[i] += f;
            return;
        }
        if(asum[i] == 0.f){     /* if empty slot store f */
            asum[i] = f;
            return;
        }
        f += asum[i];           /* else add slot to f, clear slot */
        asum[i] = 0.f;          /* and continue until empty slot */
    }
}

/* return sum from array */
float returnsum(float asum[256])
{
float sum = 0.f;
size_t i;
    for(i = 0; i < 256; i++)
        sum += asum[i];
    return sum;
}

双精度示例：

/* clear array */
void clearsum(double asum[2048])
{
size_t i;
    for(i = 0; i < 2048; i++)
        asum[i] = 0.;
}

/* add a number into array */
void addtosum(double d, double asum[2048])
{
size_t i;
    while(1){
        /* i = exponent of d */
        i = ((size_t)((*(unsigned long long *)&d)>>52))&0x7ff;
        if(i == 0x7ff){         /* max exponent, could be overflow */
            asum[i] += d;
            return;
        }
        if(asum[i] == 0.){      /* if empty slot store d */
            asum[i] = d;
            return;
        }
        d += asum[i];           /* else add slot to d, clear slot */
        asum[i] = 0.;           /* and continue until empty slot */
    }
}

/* return sum from array */
double returnsum(double asum[2048])
{
double sum = 0.;
size_t i;
    for(i = 0; i < 2048; i++)
        sum += asum[i];
    return sum;
}

Answer 11

关于排序，在我看来，如果您希望取消，那么数字应该以 降序 的数量级添加，而不是升序。例如：

((-1 + 1) + 1e-20) 将给出 1e-20

但是

((1e-20 + 1) - 1) 将给出 0

在第一个等式中，两个大数被抵消，而在第二个等式中，1e-20 项在添加到 1 时会丢失，因为没有足够的精度来保留它。

另外，pairwise summation 非常适合对大量数字求和。