“银行家四舍五入”真的在数值上更稳定吗？

Question

我的意思是银行家的四舍五入

1. “四舍五入到最近，平成平”

作为recommended by IEEE 754:

四舍五入到最接近的值；如果该数字在中途下降，则将其四舍五入到具有偶数（零）最低有效位的最接近的值。这是二进制浮点的默认值和十进制的推荐默认值。

据说这种方法优于

1. “四舍五入，远离零”

on the grounds that 它“在对四舍五入的数字求和时最小化预期误差”。显然这是because“在最合理的分布上，它不会像从零接近一半的方法那样受到负面或正面偏差的影响”。

我不明白为什么会这样。直观地说，如果0.0 朝零舍入，0.5“应该”从零舍入（如方法 2 中所示）。这样一来，等量的数字就会朝零舍入和远离零。简而言之，如果浮点数用 1 个十进制数字表示，则在 10 个数字中，0.0、...、0.9 将向下舍入 5 个，使用方法 2 将 5 个向上舍入。@987654328 类似@, ..., 1.9 等等。

当然浮点数用二进制尾数表示，但我认为上述推理仍然适用。请注意，对于 IEEE 754 双精度，整数和整数加半数都可以精确表示为近似高达 2^52 的绝对值，因此这些精确值实际上会在实践中显示。

那么方法1如何更好？

Answer 1

是的！它确实在数值上更稳定。

对于您正在查看的情况，数字 [0.0, 0.1, ..., 0.9]，请注意，在round-ties-to-away 下，这些数字中只有 四个 向下舍入（0.1 到0.4)，5 个向上舍入，1 个 (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 版本无关的舍入函数：

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 算法计算浮点数列表的总和。该算法将要计算的总和分成两个大致相等的部分，递归地对这两个部分求和，然后将结果相加。它比简单的求和要准确得多，但通常不如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 的普通加法算法，该算法在典型机器上将解析为标准的 IEEE 754 加法，使用 round-ties-to-even 舍入模式。

我们想要查看来自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 总和并返回两个版本的错误，最后以单位测量。

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

我得到以下情节：

我本来打算更进一步，对两个样本的偏差进行统计检验，但是从零关系方法产生的偏差非常明显，看起来没有必要。有趣的是，虽然 tie-off-from-zero 方法给出的结果较差，但它确实给出了更小的错误分布。