"银行家的四舍五入"真的在数值上更稳定吗？
Lui*_*ndo 9 language-agnostic floating-point rounding floating-accuracy
通过银行家的舍入我的意思是
- "绕到最近,连接到均匀"
按照IEEE 754的建议:
舍入到最接近的值; 如果数字中间下降,则使用偶数(零)最低有效位舍入到最接近的值.这是二进制浮点的默认值,也是十进制的推荐默认值.
据说该方法优于
- "圆到最近,与零相关"
理由是它"在对四舍五入的数字求和时最小化预期误差".显然,这是因为 "它不会受到负面或正偏差的影响,就像在大多数合理分布上从零方法的一半那样".
我不明白为什么会这样.直观地,如果0.0向零舍入,则0.5"应该"从零舍入(如方法2中所示).这样,相同数量的数字将朝向零并且远离零.换句话说,如果浮动数字用1个十进制数字表示,十个数字中的数字0.0,......,0.9五个将向下舍入,五个用方法2进行四舍五入.类似地1.0,... 1.9等等.
当然,浮点数用二进制尾数表示,但我认为上述推理仍然适用.请注意,对于IEEE 754双精度,整数和整数加半精度可以精确表示绝对值到2^52大约,因此这些精确值实际上会显示出来.
那么方法1如何更好？
Mar*_*son 16
是! 它确实在数值上更稳定.
对于您正在查看的情况,数字[0.0, 0.1, ..., 0.9],请注意,在圆形关系下,这些数字中只有四个向下舍入(0.1通过0.4),五个向上舍入,一个(0.0)由舍入不变操作,然后当然,这种模式重复1.0贯穿1.9,2.0通过2.9等等.因此,平均而言,更多的值从零到舍入而不是朝向它.但在圆形关系下,我们会得到:
- 五个值四舍五入,四个四舍五入
[0.0, 0.9] - 四个值四舍五入,五个四舍五入
[1.0, 1.9]
等等.平均而言,我们得到的数值相同,即四舍五入.更重要的是,舍入引入的预期误差(在对输入分布的适当假设下)接近于零.
这是使用Python的快速演示.为了避免因内置round函数中的Python 2/Python 3差异而导致的困难,我们提供了两个与Python版本无关的舍入函数:
def round_ties_to_even(x):
"""
Round a float x to the nearest integer, rounding ties to even.
"""
if x < 0:
return -round_ties_to_even(-x) # use symmetry
int_part, frac_part = divmod(x, 1)
return int(int_part) + (
frac_part > 0.5
or (frac_part == 0.5 and int_part % 2.0 == 1.0))
def round_ties_away_from_zero(x):
"""
Round a float x to the nearest integer, rounding ties away from zero.
"""
if x < 0:
return -round_ties_away_from_zero(-x) # use symmetry
int_part, frac_part = divmod(x, 1)
return int(int_part) + (frac_part >= 0.5)
现在我们来看看通过在范围内的一位数后点小数值上应用这两个函数而引入的平均误差[50.0, 100.0]:
>>> test_values = [n / 10.0 for n in range(500, 1001)]
>>> errors_even = [round_ties_to_even(value) - value for value in test_values]
>>> errors_away = [round_ties_away_from_zero(value) - value for value in test_values]
我们使用最近添加的statistics标准库模块来计算这些错误的平均值和标准差:
>>> import statistics
>>> statistics.mean(errors_even), statistics.stdev(errors_even)
(0.0, 0.2915475947422656)
>>> statistics.mean(errors_away), statistics.stdev(errors_away)
(0.0499001996007984, 0.28723681870533313)
这里的关键点是errors_even零均值:平均误差为零.但errors_away有正面意思:平均误差偏离零.
一个更现实的例子
这是一个半现实的例子,它展示了数值算法中从零到零的偏差.我们将使用成对求和算法计算浮点数列表的总和.该算法将要计算的总和分成两个大致相等的部分,递归地对这两个部分求和,然后将结果相加.它比天真的总和更准确,但通常不如Kahan求和更复杂的算法.这是NumPy sum函数使用的算法.这是一个简单的Python实现.
import operator
def pairwise_sum(xs, i, j, add=operator.add):
"""
Return the sum of floats xs[i:j] (0 <= i <= j <= len(xs)),
using pairwise summation.
"""
count = j - i
if count >= 2:
k = (i + j) // 2
return add(pairwise_sum(xs, i, k, add),
pairwise_sum(xs, k, j, add))
elif count == 1:
return xs[i]
else: # count == 0
return 0.0
我们在add上面的函数中包含了一个参数,表示要用于添加的操作.默认情况下,它使用Python的常规加法算法,在典型的机器上将使用round-ties-to-even舍入模式解析为标准的IEEE 754加法.
我们希望从pairwise_sum函数中查看预期的误差,使用标准加法和使用从零添加的圆形连接版本.我们的第一个问题是我们没有一种简单易用的方法来从Python中改变硬件的舍入模式,而二进制浮点的软件实现会很大而且很慢.幸运的是,我们可以使用一种技巧在仍然使用硬件浮点的同时获得零距离.对于该技巧的第一部分,我们可以使用Knuth的"2Sum"算法来添加两个浮点数并获得正确舍入的和以及该总和中的确切误差:
def exact_add(a, b):
"""
Add floats a and b, giving a correctly rounded sum and exact error.
Mathematically, a + b is exactly equal to sum + error.
"""
# This is Knuth's 2Sum algorithm. See section 4.3.2 of the Handbook
# of Floating-Point Arithmetic for exposition and proof.
sum = a + b
bv = sum - a
error = (a - (sum - bv)) + (b - bv)
return sum, error
有了这个,我们可以很容易地使用误差项来确定精确总和何时达到平局.我们有一个平局,当且仅当它error是非零并且sum + 2*error是完全可表示的,并且在那种情况下sum并且sum + 2*error是最接近该平局的两个浮点数.使用这种想法,下面是两个数相加,并给出了正确的舍入结果,但查房关系的功能远从零.
def add_ties_away(a, b):
"""
Return the sum of a and b. Ties are rounded away from zero.
"""
sum, error = exact_add(a, b)
sum2, error2 = exact_add(sum, 2.0*error)
if error2 or not error:
# Not a tie.
return sum
else:
# Tie. Choose the larger of sum and sum2 in absolute value.
return max([sum, sum2], key=abs)
现在我们可以比较结果.sample_sum_errors是一个函数,它生成[1,2]范围内的浮点数列表,使用正常的round-ties-to-even加法和我们的自定义round-ties-from-zero版本添加它们,与精确的总和进行比较并返回两个版本的错误,以最后一个单位为单位进行测量.
import fractions
import random
def sample_sum_errors(sample_size=1024):
"""
Generate `sample_size` floats in the range [1.0, 2.0], sum
using both addition methods, and return the two errors in ulps.
"""
xs = [random.uniform(1.0, 2.0) for _ in range(sample_size)]
to_even_sum = pairwise_sum(xs, 0, len(xs))
to_away_sum = pairwise_sum(xs, 0, len(xs), add=add_ties_away)
# Assuming IEEE 754, each value in xs becomes an integer when
# scaled by 2**52; use this to compute an exact sum as a Fraction.
common_denominator = 2**52
exact_sum = fractions.Fraction(
sum(int(m*common_denominator) for m in xs),
common_denominator)
# Result will be in [1024, 2048]; 1 ulp in this range is 2**-44.
ulp = 2**-44
to_even_error = (fractions.Fraction(to_even_sum) - exact_sum) / ulp
to_away_error = (fractions.Fraction(to_away_sum) - exact_sum) / ulp
return to_even_error, to_away_error
这是一个示例运行:
>>> sample_sum_errors()
(1.6015625, 9.6015625)
所以这是使用标准加法的1.6 ulps的误差,并且当将线条从零处舍入时误差为9.6 ulps.这当然看起来好像关系外之调零法是差,但单次运行是不是特别有说服力.让我们这样做10000次,每次使用不同的随机样本,并绘制我们得到的错误.这是代码:
import statistics
import numpy as np
import matplotlib.pyplot as plt
def show_error_distributions():
errors = [sample_sum_errors() for _ in range(10000)]
to_even_errors, to_away_errors = zip(*errors)
print("Errors from ties-to-even: "
"mean {:.2f} ulps, stdev {:.2f} ulps".format(
statistics.mean(to_even_errors),
statistics.stdev(to_even_errors)))
print("Errors from ties-away-from-zero: "
"mean {:.2f} ulps, stdev {:.2f} ulps".format(
statistics.mean(to_away_errors),
statistics.stdev(to_away_errors)))
ax1 = plt.subplot(2, 1, 1)
plt.hist(to_even_errors, bins=np.arange(-7, 17, 0.5))
ax2 = plt.subplot(2, 1, 2)
plt.hist(to_away_errors, bins=np.arange(-7, 17, 0.5))
ax1.set_title("Errors from ties-to-even (ulps)")
ax2.set_title("Errors from ties-away-from-zero (ulps)")
ax1.xaxis.set_visible(False)
plt.show()
当我在我的机器上运行上述功能时,我看到:
Errors from ties-to-even: mean 0.00 ulps, stdev 1.81 ulps
Errors from ties-away-from-zero: mean 9.76 ulps, stdev 1.40 ulps
我得到以下情节:
我打算更进一步,对两个样本的偏差进行统计检验,但是零关系方法的偏差是如此明显,以至于看起来没必要.有趣的是,虽然从零开始的方法给出了较差的结果,但它确实给出了较小的误差扩散.