Kahan's Summation Formula原理—它是如何处理大数吃小数的

Kahan's Summation Formula原理—它是如何避免大数吃小数的

Kahan求和公式原理：

       首先，这个算法就是用来求和的，求a1+a2+a3+...为什么不直接相加呢，而要用Kahan求和公式呢，这个算法的用武之地在哪呢，一一道来

       kahan求和算法能避免大数吃小数的情况。

       大数吃小数是什么意思呢？举个例子，我们用两个float相加，float是32位，它的精度是小数点后6-7位(详见http://blog.csdn.net/zhangpinghao/article/details/8138732），设有a=123456;b=2.189；a+b应该是123458.189但是由于float的精度只有小数点后6-7位，所以必然得不到123458.189，后面的89可能会截掉，8不一定，9是必然会截掉的。好的，才做一个加法就产生至少了0.009的误差，做1000个这样的加法，误差就是9了，这显然不是我们想要的。

       kahan求和算法可以避免这种情况，它有一个数用来记住那个被截断的小数，同样做下面的计算，设有a=123456;b=2.189；计算a+b。kahan求和算法是这样做的：sum=a+b（不准确）; temp= (a+b)-a-b;temp等于多少呢，初看这不就是0吗？不是的，计算机此时算的可不是0，而是等于-0.009，就是被截断的那个小数。通过一个临时变量我们就记住了这个误差，当计算下一个加法的时候，可以把这个误差补上，并且更新误差到sum。

      其实也可以这样理解，sum不是由于数太大，占用了小数的精度吗，而这个小数在当前一步看似是可以忽略的，但是由于，迭代的次数旁道，小数会累积成大误差，那么我们另外用的float专门记住这个误差小数不就得了吗。

详细的可以看维基百科，不过是英文的

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

Kahan summation algorithm 的具体算法如下：

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

 
 

  
  #include "stdafx.h"
 
 
 
 

  
  int _tmain(int argc, _TCHAR* argv[])
 
 
 
 

  
  {
 
 
 
 

  
  int i;
 
 
 
 

  
  float x=0.001;
 
 
 
 

  
  float y;
 
 
 
 

  
  float t;
 
 
 
 

  
  float sum;
 
 
 
 

  
  float eps=0;
 
 
 
 

  
  

 
 
 
 

  
  printf("\n理论值：%f\n",x*1000000);
 
 
 
 

  
  

 
 
 
 

  
  sum=0;
 
 
 
 

  
  for(i=0;i<1000000;i++)
 
 
 
 

  
  {
 
 
 
 

  
  sum+=x;
 
 
 
 

  
  }
 
 
 
 

  
  printf("\n累加值：%f\n",sum);
 
 
 
 

  
  

 
 
 
 

  
  sum=0;
 
 
 
 

  
  for(i=0;i<1000000;i++)
 
 
 
 

  
  {
 
 
 
 

  
  y=x-eps;
 
 
 
 

  
  t=sum+y;
 
 
 
 

  
  eps=(t-sum)-y;
 
 
 
 

  
  sum=t;
 
 
 
 

  
  }
 
 
 
 

  
  printf("\nKahan累加值：%f\n",sum);
 
 
 
 

  
  

 
 
 
 

  
  getchar();
 
 
 
 

  
  return 0;
 
 
 
 

  
  }
 
 
 
 

  
  运行结果：
 
 
 
 

  
  Summation Formula" name=image_operate_671312371788031 alt="Kahan's Summation Formula" src="http://s5.sinaimg.cn/middle/59154d37ga9985d8fd894&690" width=677 height=442 real_src="http://s5.sinaimg.cn/middle/59154d37ga9985d8fd894&690">
 
 
 
 

  
   
 
 
 
 

  
  保持精度的小trick：Kahan求和

  
  
由于最近用GPU编程，涉及到了float数组，就不得不涉及精度问题。对于双精度如C中double以及Fortran中real(kind = 8)，一般运算的精度足以保持，但是单精度数组，在大量操作后极易出现“大数吃小数”等不稳定现象。在不能使用更高精度数组的前提下，可以用一个小技巧来保持精度：Kahan求和。
  
  
1 见下面一段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，继续起保存损失的作用。
例如，在GPU计算大型矩阵或向量时，如果涉及到reduction操作，可以在加和中使用这种技巧。当然GPU计算能力在2.x是可以有双精度运算，但是要比单精度慢20倍左右，一般不经常使用。
  
  
 
  
  
注意：有人称一些优化激进的编译器可能会做这样的化简：eps=(t-s)-y=(s+y-s)-s=0。我试验了intel fortran和GNU fortran，包括开启了-O3，-funsafe-math-optimizations等方法，都没有问题。原来ANIS C标准规定不允许做表达式重排序优化。
 
 

做项目终于遇到精度产生的问题了，由于免费的CULA库只能用float精度，所以程序里有些数组是float的，然后在做累加运算的时候问题就出来了——计算小节点系统就收敛，大节点系统发散，明显是精度问题。于是研究了一下Kahan 求和公式，这货果然厉害，能将累加的精度保持在float级别。

以下是我用VC写的演示程序：

#include "stdafx.h"

int _tmain(int argc, _TCHAR* argv[])

{

int i;

float x=0.001;

float y;

float t;

float sum;

float eps=0;

printf("\n理论值：%f\n",x*1000000);

sum=0;

for(i=0;i<1000000;i++)

{

sum+=x;

}

printf("\n累加值：%f\n",sum);

sum=0;

for(i=0;i<1000000;i++)

{

y=x-eps;

t=sum+y;

eps=(t-sum)-y;

sum=t;

}

printf("\nKahan累加值：%f\n",sum);

getchar();

return 0;

}

运行结果：