"banker's rounding" 真的在数值上更稳定吗？
Is "banker's rounding" really more numerically stable?
我的意思是银行家四舍五入
- "round to nearest, ties to even"
Rounds to the nearest value; if the number falls midway it is rounded to the nearest value with an even (zero) least significant bit. This is the default for binary floating-point and the recommended default for decimal.
据说这种方法优于
- "round to nearest, ties away from zero"
on the grounds that it "minimizes the expected error when summing over rounded figures". Apparently this is because "it does not suffer from negative or positive bias as much as the round half away from zero method over most reasonable distributions".
我不明白为什么会这样。直观地说，如果 0.0 向零四舍五入，0.5 "should" 将从零四舍五入（如方法 2 中所示）。这样，等量的数字将四舍五入到零或远离零。简单来说，如果浮点数用 1 位小数表示，那么十个数字 0.0, ..., 0.9 中有五个会向下舍入，五个会用方法 2 向上舍入。 1.0, ..., 1.9 等
同样
当然浮点数是用二进制尾数来表示的，但我认为上面的推理仍然适用。请注意，对于 IEEE 754 双精度，整数和加半整数都可以精确地表示为大约 2^52 的绝对值，因此这些精确值实际上会出现在实践中。
那么方法 1 有何优点？
是的！它确实在数值上更稳定。
对于您正在查看的情况，数字 [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 的快速演示。为了避免由于 Python 2 / Python 3 内置 round 函数的差异，我们给出了两个 Python-version-agnostic 舍入函数：
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 具有正均值：平均误差偏离零。
一个更现实的例子
这是一个半现实的例子，它演示了数值算法中远离零的舍入关系的偏差。我们将使用 pairwise summation algorithm. This algorithm breaks the sum to be computed into two roughly equal parts, recursively sums those two parts, then adds the results. It's substantially more accurate than a naive sum, but typically not as good as more sophisticated algorithms like Kahan summation 计算浮点数列表的总和。这是 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”算法将两个浮点数相加并获得正确舍入的总和以及该总和中的 exact 错误：
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] 中的浮点数列表，使用正常的圆结对偶数加法和我们的自定义圆结远离零版本添加它们，与 exact 总和和 returns 两个版本的误差进行比较，在最后位置以单位测量。
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 ulp 的误差，以及从零舍入时的 9.6 ulp 的误差。它肯定 看起来 好像从零结的方法更糟糕，但是单个 运行 并不是特别有说服力。让我们这样做 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
我得到以下情节：
我本来打算更进一步，对两个样本的偏差进行统计检验，但是零值结方法的偏差非常明显，看起来没有必要。有趣的是，虽然零结结点法给出的结果较差，但它 确实 给出了更小的误差分布。
我的意思是银行家四舍五入
- "round to nearest, ties to even"
Rounds to the nearest value; if the number falls midway it is rounded to the nearest value with an even (zero) least significant bit. This is the default for binary floating-point and the recommended default for decimal.
据说这种方法优于
- "round to nearest, ties away from zero"
on the grounds that it "minimizes the expected error when summing over rounded figures". Apparently this is because "it does not suffer from negative or positive bias as much as the round half away from zero method over most reasonable distributions".
我不明白为什么会这样。直观地说，如果 0.0 向零四舍五入，0.5 "should" 将从零四舍五入（如方法 2 中所示）。这样，等量的数字将四舍五入到零或远离零。简单来说，如果浮点数用 1 位小数表示，那么十个数字 0.0, ..., 0.9 中有五个会向下舍入，五个会用方法 2 向上舍入。 1.0, ..., 1.9 等
当然浮点数是用二进制尾数来表示的，但我认为上面的推理仍然适用。请注意，对于 IEEE 754 双精度，整数和加半整数都可以精确地表示为大约 2^52 的绝对值，因此这些精确值实际上会出现在实践中。
那么方法 1 有何优点？
是的！它确实在数值上更稳定。
对于您正在查看的情况，数字 [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 的快速演示。为了避免由于 Python 2 / Python 3 内置 round 函数的差异，我们给出了两个 Python-version-agnostic 舍入函数：
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 具有正均值：平均误差偏离零。
一个更现实的例子
这是一个半现实的例子，它演示了数值算法中远离零的舍入关系的偏差。我们将使用 pairwise summation algorithm. This algorithm breaks the sum to be computed into two roughly equal parts, recursively sums those two parts, then adds the results. It's substantially more accurate than a naive sum, but typically not as good as more sophisticated algorithms like Kahan summation 计算浮点数列表的总和。这是 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”算法将两个浮点数相加并获得正确舍入的总和以及该总和中的 exact 错误：
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] 中的浮点数列表，使用正常的圆结对偶数加法和我们的自定义圆结远离零版本添加它们，与 exact 总和和 returns 两个版本的误差进行比较，在最后位置以单位测量。
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 ulp 的误差，以及从零舍入时的 9.6 ulp 的误差。它肯定 看起来 好像从零结的方法更糟糕，但是单个 运行 并不是特别有说服力。让我们这样做 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
我得到以下情节：
我本来打算更进一步，对两个样本的偏差进行统计检验，但是零值结方法的偏差非常明显，看起来没有必要。有趣的是，虽然零结结点法给出的结果较差，但它 确实 给出了更小的误差分布。