为了得到最精确的结果,应该以哪种顺序添加浮点数?

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

如果我们有一些函数,它会累加浮点数:

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

v是一个std::vector<float> ,例如。

  • 在积累这些数字之前sorting这些数字会更好吗?

  • 哪个命令会给出最准确的答案?

我怀疑,按升序排列数字实际上会使数字错误更less ,但不幸的是我无法certificate自己。

PS我意识到这可能与现实世界的编程无关,只是好奇。

你的本能基本上是正确的,按升序(大小)sorting通常有所改善。 考虑一下我们添加单精度(32位)浮点数的情况,并且有十亿个数值等于1 /(10亿),一个数值等于1.如果第一个数字先出现,那么总和就会来到1,因为1 +(1/10亿)是1,因为精度的损失。 总的来说,每增加一点都没有效果。

如果小数值先到,他们至less会总结一些东西,虽然即使这样我也有2 ^ 30个,而在2 ^ 25左右之后,我又回到了个别不影响总数的情况任何更多。 所以我仍然需要更多的技巧。

这是一个极端的情况,但通常添加两个类似幅度的值比添加两个非常不同的幅度值更精确,因为您以这种方式“丢弃”更小的精度位数。 通过对数字进行sorting,您可以将相似数量的数值组合在一起,然后按照递增顺序将这些数值相加,从而使小数值成为累积地达到更大数字量级的“机会”。

不过,如果涉及负数,很容易“超越”这种方法。 考虑三个值的总和, {1, -1, 1 billionth} 。 算术上正确的总数是1 billionth ,但是如果我的第一个加法涉及微小的值,那么我的最终总和将是0.在六个可能的命令中,只有两个是“正确的” – {1, -1, 1 billionth}{-1, 1, 1 billionth} 。 所有的6个命令都给出了在input的最大值(0.0000001%out)的范围内准确的结果,但是其中4个结果在真实解的范围内(100%输出)是不准确的。 你正在解决的具体问题将告诉你,前者是否够好。

事实上,你可以发挥更多的技巧,而不是按照sorting顺序添加它们。 如果你有很多很小的数值,中等数值的中等数值和less量的大数值,那么把所有的小数值加起来,然后分别总和中等数值,加上这两个总数一起然后添加大的。 find最精确的浮点加法组合并不重要,但为了处理非常糟糕的情况,您可以保持不同大小的所有运行总数,将每个新值添加到与其大小最匹配的总和中,当一个跑步总数开始变得过大时,把它join下一个总起来,并开始一个新的。 采取其逻辑的极端,这个过程相当于执行任意精度types的总和(所以你会这样做)。 但是考虑到按照升序或降序排列的简单select,上升是更好的select。

它与现实世界的编程有一些关系,因为在某些情况下,如果你不小心剔除了一个由大量值组成的“重尾”尾部,那么你的计算就会非常糟糕这个总和,或者如果你从很多小的值中丢掉了太多的精度,而这些小的值只会影响总和的最后几位。 如果尾巴可以忽略不计,你可能不在乎。 例如,如果您只是在一起添加less量的值,而只使用总和的几个有效数字。

还有一种为这种积累操作devise的algorithm,叫做Kahan Summation ,你应该知道。

根据维基百科,

Kahan求和algorithm (也称为补偿求和 )显着地减less了通过添加一系列有限精度浮点数而获得的总数中的数值误差,与明显的方法相比。 这是通过保持一个单独的运行补偿(一个variables来累积小错误)来完成的。

在伪代码中,algorithm是:

 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 

我尝试了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) 

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

如果我在上面的代码中将double s更改为float ,则得到:

 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 summation(with double s)

 #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) 

即使我改变了double float在上面的代码中,我得到:

 Kahan sum = 2.000000000000000 (difference = 0.000000000000000) 

看来卡汉是要走的路!

有一类algorithm可以解决这个确切的问题, 而不需要对数据进行sorting或重新sorting

换句话说,总和可以在数据上一次完成。 这也使得这样的algorithm适用于数据集未知的情况下,例如,如果数据实时到达并且需要保持运行总和。

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

我们提出了一种新颖的在线algorithm,用于精确求和浮点数字stream。 通过“在线”,我们的意思是algorithm一次只需要看到一个input,并且可以采取任意长度的这种input的inputstream,而只需要固定的存储器。 “精确”是指我们algorithm的内部数组的总和正好等于所有input的总和,返回的结果是正确的四舍五入总和。 所有input(包括非标准化数字,但模中间溢出)的正确性certificate都是有效的,并且与加数的数量或总和的条件数无关。 该algorithm渐近地每个加法器只需要5个FLOP,并且由于指令级的并行运行,当加数数目大于10,000时,比明显的,快速但笨的“普通recursion求和”循环慢2-3倍左右。 因此,据我们所知,它是已知algorithm中最快,最准确和最有效的存储器。 事实上,很难看出如果没有硬件改进,更快的algorithm或需要显着更less的FLOP的algorithm是可以存在的。 提供了大量的加数的应用程序。

来源: algorithm908:浮点stream在线精确求和 。

在史蒂夫首先按照升序对数字进行sorting的答案的基础上,我会再介绍两点:

  1. 决定两个数字的指数之间的差异,你可能会认为你会失去太多的精度。

  2. 然后按顺序添加数字,直到累加器的指数对于下一个数字太大,然后将累加器放到一个临时队列中,然后用下一个数字启动累加器。 继续,直到你用尽了原始列表。

你用临时队列(已经sorting)和指数可能更大的差异来重复这个过程。

如果你必须一直计算指数,我认为这将是相当缓慢的。

我有一个快速的程序,结果是1.99903

在积累之前,我认为你可以比sorting更好,因为在积累的过程中,累加器变得越来越大。 如果你有大量的相似的数字,你会开始很快失去精度。 这是我会build议,而不是:

 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 

当然这个algorithm对于优先队列而不是列表来说是最有效的。 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数是最大的数字,但是我们希望最小 。 我可以给队列提供更多的模板参数,但是这种方法看起来更简单。

这并不完全回答你的问题,但一个聪明的事情是要两次运行总和,一次是舍入模式 “四舍五入”,一次是“ 舍入 ”。 比较两个答案,你知道/如何/不准确的结果,如果你因此需要使用聪明的总结策略。 不幸的是,大多数语言都不会像改变浮点舍入模式一样简单,因为人们并不知道它在日常计算中实际上是有用的。

看一下Intervalalgorithm ,你可以像这样做所有的math运算 ,随时保持最高和最低值。 它导致一些有趣的结果和优化。

提高精度的最简单的sorting是按升序绝对值sorting。 这使得最小的幅度值有机会累积或取消之前与更大的幅度值相互作用,会导致精度的损失。

也就是说,你可以通过跟踪多个不重叠的部分总和来做得更好。 下面是一篇论文,描述了这种技术并提出了一个准确性certificate:www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

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

对于IEEE 754单精度或双精度或已知格式的数字,另一种方法是使用由指数索引的一组数字(由调用者传递,或在C ++的类中)。 将数字添加到数组中时,只添加指数相同的数字(直到find一个空插槽并存储数字)。 当要求求和时,数组从最小到最大以减less截断。 单精度示例:

 /* 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; } 

你的花车应该以双精度join。 这会给你比其他任何技术更高的精度。 为了更精确和速度更快,您可以创build四个总和,并在最后添加它们。

如果要添加双精度数字,则使用long double作为总和 – 但是,这对于long double实际上比double更精确(通常是x86,PowerPC取决于编译器设置)的实现只会起到积极作用。

关于sorting,在我看来,如果你期望取消,那么数字应该按照大小递增的顺序添加,而不是递增的。 例如:

((-1 + 1)+ 1e-20)会得到1e-20

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

在第一个方程中,两个大数字被抵消,而第二个1e-20的术语在加到1时丢失,因为没有足够的精度来保留它。

而且, 两两相加对于大量的数字来说是相当不错的。