Kahan's summation Formula

保持精度的小trick:Kahan's summation Formula

由于最近用GPU编程,涉及到了float数组,就不得不涉及精度问题。在 CPU 上进行计算时,我们使用 double(即 64 bits 浮点数)来累进计算过程,而在 GPU 上则只能用 float(32 bits 浮点数)。在累加大量数字的时候,由于累加结果很快会变大,因此后面的数字很容易被舍去过多的位数。对于双精度如C中double以及Fortran中real(kind = 8),一般运算的精度足以保持,但是单精度数组,在大量操作后极易出现“大数吃小数”等不稳定现象。在不能使用更高精度数组的前提下,可以用一个小技巧来保持精度:Kahan求和。

function KahanSum(input)
    var sum = 0.0
    var c = 0.0             //A running compensation for lost low-order bits.
    for i = 1 to input.length do
        y = input[i] - c    //So far, so good: c is zero.
        t = sum + y         //Alas, sum is big, y small, so low-order digits of y are lost.
        c = (t - sum) - y   //(t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
        sum = t             //Algebraically, c should always be zero. Beware eagerly optimising compilers!
        //Next time around, the lost low part will be added to y in a fresh attempt.
    return sum

 

见下面一段Fortran代码:

program main
implicit none
    integer, parameter:: N = 1000000
    integer i
    real(kind = 4), parameter:: ELEMENT = 0.001
    real(kind = 4) s, eps, y, t

    write (*, "('Theoretical value: ', F10.5)") N*ELEMENT

    s = 0.0
    do i = 1, N
        s = s+ELEMENT
    enddo
    write (*, "('Naive method value: ', F10.5)") s

   s = 0.0
    eps = 0.0
    do i = 1, N
        y = ELEMENT-eps
        t = s+y
        eps = (t-s)-y
        s = t
    enddo
    write (*, "('Kahan method value: ',F10.5)") s

stop
end

运行结果为:

Theoretical value: 1000.00006
Naive method value:  991.14154
Kahan method value: 1000.00006

对于N个0.001,普通方法累加到991左右就已经丢失精度了。可以看到用“Kahan method”能够得到近乎于理论的精度数值。分析一下他的原理。我们发现,如果没有精度损失,eps永远为0,y就是ELEMENT=0.001。一旦在 i 到了某个数值出现了大数吃小数 的情形时,不妨激进的设小数部分全部被截断,则如s = 991.0000时,由于eps之前为0,则y=0.0010.之后t=s+y,得到的就是“吃掉”的结果,如991.0000,绝对误差达0.001.此时:eps=(t-s)-y=(991.0000-991.0000)-0.0010=-0.001,可见eps起了保存“损失位”的作用。此时s=t=991.0000.下个循环:y = 0.001--0.001=0.002,t = s+y=991.0000,eps=-0.002,如此反复,这样足够多循环后,eps足可以复现大的校正值,从而保证结果的高精度。当eps足够大时候,(t-s)-y=0,从而使eps重新为0,继续起保存损失的作用。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值