https://en.wikipedia.org/wiki/Kahan_summation_algorithm

目的:用于降低有限精度浮点数组成的数列的求和误差

普通算法:最坏情形误差正比于 \(n\) 增长,随机输入情形均方误差正比于 \(\sqrt{n}\)。

补偿求和:最坏情形误差随 \(n\) 的增长要慢许多。

Python 代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def kahanSum(vals):
    s = 0.0
    c = 0.0
    for v in vals:
        # y = v - err_of_previous_step
        y = v - c
        # t = s + v - err_of_previous_step + err_of_this_step
        t = s + y
        # c = err_of_this_step
        c = (t - s) - y
        # s = v0 + v1 + ε1 + v2 + ε2 - ε1 + v3 + ε3 - ε2 + ...
        #   = v0 + v1 + v2 + v3 + ... + vn + εn
        s = t
    return s

试验

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
vals = [100000, 2.8, 2.7]
# 在只有六位有效数字的情况下,简单求和的过程是
s0 = 100000
s1 = 100002
s2 = 100004
# 使用补偿算法后的过程是
y0 = 100000
t0 = 100000
c0 = 0
s0 = 100000
y1 = 2.8
t1 = 100002
c1 = -0.8
s1 = 100002
y2 = 3.5
t2 = 100005
c2 = -0.5
s2 = 100005
# 这个当然更接近精确结果

误差界

对于求和 \[ S_n = \sum_{i=1}^n x_i, \] 其误差为 \[ |E_n| \le \left[2\varepsilon + O(n\varepsilon^2)\right] \sum_{i=1}^n|x_i|. \]

Condition number

定义为 \[ \frac{\sum_{i=1}^n \left|x_i\right|}{\left|\sum_{i=1}^n x_i\right|}. \]

给定条件数的情形,由于 \(n\) 通常远小于 \(1/\varepsilon\),误差实质上与 \(n\) 无关。

计算数列方差时也应该避免这种舍入误差。

其它方法: 成对求和

https://en.wikipedia.org/wiki/Pairwise_summation

优点:加法次数与普通的简单求和一样; 应用于快速傅立叶变换

误差: \(O(\varepsilon\sqrt{\log n})\) (假定舍入误差的正负号随机); 最坏情形 \(O(\varepsilon{\log n})\)

直觉理解: 为什么分治求和可以降低误差?

1
2
3
4
5
6
7
8
9
def pairwiseSum(vals):
    if len(vals) <= N:
        s = 0.0
        for v in vals:
            s += v
    else:
        m = floor(len(vals) / 2)
        s = pairwiseSum(vals[0:m]) + pairwiseSum(vals[m:])
    return s