使用 IDL 的 TOTAL 函数最小化舍入误差
15593 对此文章评分:
1.0
使用 IDL 的 TOTAL 函数最小化舍入误差
匿名作者 2016年10月6日，星期四
我本周博客主题的想法来自于 Ron Kneusel 提出的一个关于使用 IDL 的 TOTAL 函数时精度损失的问题。IDL 函数的实现方式通常是为了使其运行得更快，因为处理速度至关重要。这对 TOTAL 函数来说当然也是如此。如果你没有使用任何关键字和额外参数，那么它基本上是按输入数组的顺序进行累加求和。例外情况是，如果你的数组中有足够多的元素触发了多线程代码路径，对此我暂且不作讨论。
因此，精度损失的问题在以下情况下最容易观察到：如果你从最大的项开始，然后不断添加越来越小的项。在某个时刻，与庞大的和数相比，小项变得如此之小，以至于额外的项会在舍入过程中丢失。
以下是一个我们已知其和应趋近于 128（在第 52 位小数处存在偏差）的求和示例：
IDL> data = (1-1d/128) ^ lindgen(15000)
IDL> 128 - total(data)
5.2580162446247414e-013
所以，在这个例子中，我们在第 13 位小数处出现了误差。这个数组从最大的项开始，以最小的项结束，这对 TOTAL 的实现来说是最坏的情况。如果我们颠倒这些项的顺序，会得到：
IDL> 128 - total(reverse(data))
5.6843418860808015e-014
在这种情况下，误差缩小了 10 倍。我决定将 Ron 建议的算法（称为 Kahan 求和）与我之前用于在大规模并行 GPU（在成千上万个独立线程下运行速度很快）上执行的分治方案进行比较。Kahan 求和算法可以在维基百科上找到，其 IDL 代码相当简短（但速度非常慢）：
; Kahan 算法，源自维基百科
function KahanSum, data
sum = 0d
c = 0d
for i=0, data.length-1 do begin
y = data[i] - c
t = sum + y
c = (t - sum) - y
sum = t
endfor
return, sum
end
大规模并行算法的运行速度将比 Kahan 求和快得多，其 IDL 代码如下：
; 分治求和，
; 该算法适合大规模并行执行
function total_mp, data
compile_opt idl2,logical_predicate
n = ishft(1ull,total(ishft(1ull,lindgen(63)) lt n_elements(data),/integer))
pad = make_array(n, type=size(data,/type))
pad[0] = data[*]
while n gt 1 do begin
pad[0:n/2-1] += pad[n/2:n-1]
n /= 2
endwhile
return, pad[0]
end
以下是对递减幅度项进行求和的结果：
IDL> 128 - KahanSum(data)
5.6843418860808015e-014
IDL> 128 - total_mp(data)
5.6843418860808015e-014
因此，在这种情况下，误差与排好序（最佳情况）后输入给 TOTAL 的结果相似。作为独立测试，我还尝试了随机排序这些项，KahanSum 和 TOTAL_MP 的结果仍然在 5e-14 的数量级上保持一致：
IDL> 128 - KahanSum(data[order])
5.6843418860808015e-014
IDL> 128 - total_mp(data[order])
5.6843418860808015e-014
IDL> 128 - total(data[order])
4.1211478674085811e-013
我的结论是，TOTAL_MP 示例与 Kahan 求和算法一样准确，并且还有一个额外优势，即允许大规模并行执行，而 Kahan 算法则需要高度顺序化的执行。