应该按啥顺序添加浮点数以获得最精确的结果？
Posted
技术标签:
【中文标题】应该按啥顺序添加浮点数以获得最精确的结果？【英文标题】：In which order should floats be added to get the most precise result?应该按什么顺序添加浮点数以获得最精确的结果？ 【发布时间】：2011-10-05 15:47:24 【问题描述】：这是我在最近的采访中被问到的一个问题，我想知道（我实际上不记得数值分析的理论，所以请帮助我:)
如果我们有一些函数，它会累积浮点数：
std::accumulate(v.begin(), v.end(), 0.0);
例如，v 是 std::vector<float>。
在累积这些数字之前对它们进行排序会更好吗？
哪个顺序会给出最准确的答案？
我怀疑将数字按升序排序实际上会使数字错误少，但不幸的是我自己无法证明。
附注我确实意识到这可能与现实世界的编程无关，只是好奇。
【问题讨论】：
这实际上与现实世界的编程有关。然而，许多应用程序并不真正关心计算的绝对最佳精度，只要它“非常接近”。工程应用？极其重要。医疗应用？极其重要。大规模统计？稍微低一点的准确性是可以接受的。 请不要回答，除非您确实知道并且可以指向详细解释您的推理的页面。已经有太多关于浮点数的废话了，我们不想再添加了。如果你认为你知道。停止。因为如果你只认为你知道，那么你可能错了。 @Zéychin “工程应用？非常重要。医疗应用？非常重要。”？？？如果你知道真相，我想你会感到惊讶:) @Zeychin 绝对错误无关紧要。重要的是相对误差。如果百分之几的弧度是 0.001%，那么谁在乎呢？ 我非常推荐阅读这篇文章：“每个计算机科学家需要了解的浮点知识”perso.ens-lyon.fr/jean-michel.muller/goldberg.pdf 【参考方案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% 出）。您正在解决的特定问题会告诉您前者是否足够好。
事实上，除了按排序顺序添加它们之外，您还可以玩更多的技巧。如果您有很多非常小的值、中等数量的中等值和少量的大值，那么首先将所有小值相加，然后分别合计中等值，再将这两个总计相加可能是最准确的一起然后添加大的。找到最准确的浮点加法组合并非易事，但要应对非常糟糕的情况，您可以将整个数组的运行总计保持在不同的大小，将每个新值添加到与其大小最匹配的总数中，当一个连续的总数开始变得太大时，将其添加到下一个总数中并开始一个新的总数。从逻辑上讲，这个过程相当于以任意精度类型执行求和（所以你会这样做）。但考虑到按数量级升序或降序添加的简单选择，升序是更好的选择。
它确实与现实世界的编程有一定的关系，因为在某些情况下，如果你不小心砍掉了由大量值组成的“重”尾，每个值都太小了，你的计算可能会出现非常严重的错误单独影响总和，或者如果您从许多仅影响总和的最后几位的小值中丢弃过多的精度。在无论如何尾巴可以忽略不计的情况下，您可能不在乎。例如，如果您一开始只将少量值相加，并且只使用总和的几个有效数字。
【讨论】：
+1 用于解释。这有点违反直觉，因为加法通常在数值上是稳定的（与减法和除法不同）。 @Konrad，它可能在数值上是稳定的，但考虑到操作数的大小不同，它并不精确:) @6502：它们按数量级排序，所以 -1 排在最后。如果总和的真实值是 1，那很好。如果您将三个值相加：1 / 十亿、1 和 -1，那么您将得到 0，此时您必须回答一个有趣的实际问题——您是否需要一个准确的答案？真正的总和，还是您只需要一个在最大值范围内准确的答案？对于一些实际应用，后者已经足够好，但如果不是，则需要更复杂的方法。量子物理学使用重整化。 如果你要坚持这个简单的方案，我总是将两个最小的数字相加，然后将总和重新插入集合中。 （好吧，在这里合并排序可能效果最好。您可以使用数组中包含先前求和的数字的部分作为部分和的工作区域。） @Kevin Panko：简单的版本是单精度浮点数有 24 个二进制数字，其中最大的是数字中最大的设置位。因此，如果您将两个幅度相差超过 2^24 的数字加在一起，您将遭受较小值的全部损失，如果它们的幅度相差较小，那么您将损失相应数量的较小的精度号码。【参考方案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
【讨论】：
+1 可爱的添加到这个线程。任何“热切优化”这些语句的编译器都应该被禁止。 通过使用两个不同大小的求和变量sum 和c，这是一种几乎将精度提高一倍的简单方法。它可以简单地扩展到 N 个变量。
@ChrisA。好吧，您可以在所有重要的编译器上显式控制它（例如，通过 GCC 上的-ffast-math）。
@Konrad Rudolph 感谢您指出这是-ffast-math 的可能优化。我从这次讨论和this link 中学到的是，如果您关心数值精度，您可能应该避免使用-ffast-math，但在许多应用程序中，您可能受 CPU 限制但不关心精确的数值计算，（例如游戏编程），-ffast-math 使用是合理的。因此，我想修改我措辞强硬的“禁止”评论。
对sum, c, t, y 使用双精度变量会有所帮助。您还需要在return sum 之前添加sum -= c。【参考方案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是要走的路！
【讨论】：
我的“大”值等于 1，而不是 1e9。您的第二个答案，按大小顺序添加，在数学上是正确的（10 亿，加上十亿分之一，是 10 亿和 1），尽管更幸运的是该方法的任何一般合理性:-) 请注意double 没有将十亿分之一加在一起时不会遭受严重的精度损失，因为它有 52 个有效位，而 IEEE float 只有 24 个并且会。
@Steve，我的错误，抱歉。我已将示例代码更新为您想要的。
Kahan 的精度仍然有限，但要构造一个杀手级案例，您需要主和和误差累加器c 包含比下一个和数大得多的值。这意味着总和比主要总和要小得多，因此它们中的数量必须非常多才能加起来。尤其是 double 算术。【参考方案4】：
基于史蒂夫的答案，首先按升序对数字进行排序，我将介绍另外两个想法：
确定两个数字的指数差异，如果超过这个差异，您可能会认为您会损失太多的精度。
然后按顺序将数字相加，直到累加器的指数对于下一个数字来说太大，然后将累加器放到一个临时队列中，并以下一个数字启动累加器。继续，直到用完原始列表。
您使用临时队列（已对其进行排序）重复该过程，并且指数差异可能更大。
如果你必须一直计算指数，我认为这会很慢。
我快速运行了一个程序，结果是 1.99903
【讨论】：
【参考方案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 产生最大 数字，但我们想要最小。我本可以为队列提供更多模板参数，但这种方法似乎更简单。
【讨论】：
【参考方案6】：有一类算法可以解决这个确切的问题，无需对数据进行排序或重新排序。
换句话说，求和可以在数据上一次完成。这也使得此类算法适用于事先不知道数据集的情况，例如如果数据是实时到达的，需要维护运行总和。
这是最近一篇论文的摘要：
我们提出了一种新颖的在线算法，用于对流进行精确求和 的浮点数。 “在线”是指算法 一次只需要查看一个输入，并且可以取任意一个 此类输入的长度输入流，同时只需要常量 记忆。 “精确”是指我们的内部数组的总和 算法完全等于所有输入的总和，并且 返回的结果是正确舍入的总和。正确性证明 适用于所有输入（包括非规范化数字，但模 中间溢出），并且与被加数的数量无关 或求和的条件数。该算法渐近需要 每个加法运算只有 5 次 FLOP，并且由于指令级并行性 运行速度只比明显的、快速但愚蠢的慢 2--3 倍 当和数为时的“普通递归求和”循环 大于 10,000。因此，据我们所知，它是最快、最 在已知算法中准确且内存效率最高。的确，它 很难看出一个更快的算法或一个需要 如果没有硬件改进，可能存在的 FLOP 会显着减少。 提供大量求和的应用程序。
来源：Algorithm 908: Online Exact Summation of Floating-Point Streams。
【讨论】：
@Inverse：周围仍然有实体库。或者，在线购买 PDF 的费用为 5-15 美元（取决于您是否是 ACM 会员）。最后，DeepDyve 似乎愿意以 2.99 美元的价格将论文借出 24 小时（如果您是 DeepDyve 的新手，您甚至可以在他们的免费试用中免费获得它）：deepdyve.com/lp/acm/…【参考方案7】：这并不能完全回答您的问题，但明智的做法是计算两次总和，一次使用 rounding mode“向上取整”，一次使用“向下取整”。比较这两个答案，您就会知道/如何/您的结果不准确，因此您需要使用更聪明的求和策略。不幸的是，大多数语言并没有让浮点舍入模式变得如此简单，因为人们不知道它在日常计算中实际上是有用的。
看看Interval arithmetic，你在哪里做所有这样的数学运算，随时保持最高和最低值。它带来了一些有趣的结果和优化。
【讨论】：
【参考方案8】：提高准确性的最简单排序是按绝对值升序排序。这让最小的震级值有机会在与可能导致精度损失的较大震级值交互之前累积或取消。
也就是说，您可以通过跟踪多个不重叠的部分和来做得更好。这是描述该技术并提供准确性证明的论文：www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
该算法和其他精确浮点求和的方法是在简单的 Python 中实现的：http://code.activestate.com/recipes/393090/ 其中至少两个可以简单地转换为 C++。
【讨论】：
【参考方案9】：您的浮点数应该以双精度添加。这将为您提供比任何其他技术都更高的精度。为了获得更高的精度和更快的速度，您可以创建四个总和，并在最后将它们相加。
如果您要添加双精度数字，请使用 long double 作为总和 - 但是，这只会在 long double 实际上比 double 精度更高的实现中产生积极影响（通常是 x86，PowerPC 取决于编译器设置）。
【讨论】：
“这将为您提供比任何其他技术都更高的精度”您是否意识到您的答案是在较早的一个描述如何使用精确求和的较晚答案之后一年多？ “long double”类型太可怕了，你不应该使用它。【参考方案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;
【讨论】：
这听起来有点像Malcolm 1971 的方法，或者更像是它的变体，它使用Demmel and Hida 的指数（“算法3”）。还有另一种算法可以像您一样执行基于进位的循环，但我目前找不到。 @ZachB - 这个概念类似于bottom up merge sort for linked list，它也使用一个小数组，其中array[i] 指向具有2^i 个节点的列表。我不知道这可以追溯到多远。就我而言，这是 1970 年代的自我发现。【参考方案11】：关于排序，在我看来，如果您希望取消，那么数字应该以 降序 的数量级添加，而不是升序。例如：
((-1 + 1) + 1e-20) 将给出 1e-20
但是
((1e-20 + 1) - 1) 将给出 0
在第一个等式中，两个大数被抵消，而在第二个等式中，1e-20 项在添加到 1 时会丢失，因为没有足够的精度来保留它。
另外，pairwise summation 非常适合对大量数字求和。
【讨论】：
以上是关于应该按啥顺序添加浮点数以获得最精确的结果？的主要内容，如果未能解决你的问题，请参考以下文章