假设我有三个 32 位浮点值:a
、b
和 c
,这样 (a + b) + c != a + (b + c)
。是否有一种求和算法,可能类似于 Kahan summation,保证这些值可以以任何顺序求和,并且总是得到完全相同(相当准确)的总数?我正在寻找一般情况(即不是只处理 3 个数字的解决方案)。
任意精度算术是唯一的方法吗?我正在处理非常大的数据集,因此我希望尽可能避免使用任意精度算术的开销。
谢谢!
Let's say I have three 32-bit floating point values, a
, b
, and c
, such that (a + b) + c != a + (b + c)
. Is there a summation algorithm, perhaps similar to Kahan summation, that guarantees that these values can be summed in any order and always arrive at the exact same (fairly accurate) total? I'm looking for the general case (i.e. not a solution that only deals with 3 numbers).
Is arbitrary precision arithmetic the only way to go? I'm dealing with very large data sets, so I'd like to avoid the overhead of using arbitrary precision arithmetic if possible.
Thanks!
发布评论
评论(3)
有一个有趣的“全精度求和”算法 这里,它保证最终的总和与被加数的顺序无关(用Python给出的方法;但翻译成其他语言应该不会太困难)。请注意,该链接中给出的方法并不完全正确:主累加循环很好,但在最后一步中,将累加部分和的列表转换为单个浮点结果(< code>msumrecipe),为了得到正确舍入的结果,需要比简单地对部分和求和更加小心。请参阅配方下面的注释以及 Python 的实现(下面链接)以获取解决此问题的方法。
它确实使用某种形式的任意精度算术来保存部分和(中间和表示为双精度数的“非重叠”和),但仍然可能足够快,特别是当所有输入的大小大致相同。它总是给出正确的舍入结果,因此准确性与您所希望的一样好,并且最终的总和与被加数的顺序无关。它基于本文(自适应精度浮点算术和快速稳健几何谓词)作者:Jonathan Shewchuk。
Python 使用该算法来实现 math.fsum,它可以进行正确舍入的与顺序无关的求和;你可以在这里查看Python使用的C实现--- 查找 math_fsum 函数。
There's an interesting 'full-precision-summation' algorithm here, which guarantees that the final sum is independent of the order of the summands (recipe given in Python; but it shouldn't be too difficult to translate to other languages). Note that the recipe as given in that link isn't perfectly correct: the main accumulation loop is fine, but in the final step that converts the list of accumulated partial sums to a single floating-point result (the very last line of the
msum
recipe), one needs to be a little bit more careful than simply summing the partial sums in order to get a correctly-rounded result. See the comments below the recipe, and Python's implementation (linked below) for a way to fix this.It does use a form of arbitrary-precision arithmetic to hold partial sums (the intermediate sums are represented as 'non-overlapping' sums of doubles), but may nevertheless be fast enough, especially when all the inputs are of roughly the same magnitude. And it always gives a correctly rounded result, so accuracy is as good as you could hope for and the final sum is independent of the order of the summands. It's based on this paper (Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates) by Jonathan Shewchuk.
Python uses this algorithm for its implementation of math.fsum, which does correctly-rounded order-independent summation; you can see the C implementation that Python uses here--- look for the math_fsum function.
通过一些有关必须求和的项的附加信息,您可以避免 Shewchuk 算法的开销。
在 IEEE 754 算术中,只要
y/2
x <= 2*y
,xy
就是精确的(Sterbenz 定理,正式证明 此处)因此,如果您可以按顺序排列所有条款,使得每个部分总和为填写上面的表格,然后您就可以免费获得准确的结果。
恐怕在实践中,这种情况肯定会发生的可能性很小。正数和负数随着幅度的增加而交替可能是发生这种情况的一种情况。
注意:最初的问题是关于一种无论求和顺序如何都会给出相同结果的算法。马克的回答引发了向“精确算法”方向的转变,但再次阅读你的问题,我担心当我建议重新排序术语时我把事情推得太远了。您可能无法完成您想要做的事情,而我的回答可能是偏离主题的。好吧,抱歉:)
With some additional information about the terms you have to sum, you can avoid the overhead of Shewchuk's algorithm.
In IEEE 754 arithmetic,
x-y
is exact whenevery/2 <= x <= 2*y
(Sterbenz theorem, formally proved here)So if you can arrange all your terms in an order such that each partial sum is of the form above, then you get the exact result for free.
I am afraid that in practice there is little chance of being in conditions where this is assured to happen. Alternating positive and negatives numbers with increasing magnitudes may be one case where it happens.
Note: the original question was about an algorithm that would give the same result regardless of the summation order. Mark's answer initiated a drift in the direction of "an exact algorithm", but reading again your question, I am afraid that I am pushing things too far when I am suggesting to reorder terms. You probably can't in what you are trying to do, and my answer is probably off-topic. Well, sorry :)
在程序中进行算术运算时,我不太确定 (a + b) + c != a + (b + c) 。
然而,在当今的硬件上使用浮点运算的经验法则是永远不要直接测试相等性。
对于任何应用程序,您都应该选择足够小的 epsilon 并用作
相等性测试。
I am not quite sure that (a + b) + c != a + (b + c) when doing arithmetic in a program.
However the rule of thumb with using floating point arithmetic on present day hardware is to never directly test for equality.
For whatever application you have you should choose an epsilon that is small enough and use
as the equality test.