18

私は同様の質問を認識していますが、実用的なコストで浮動小数点数を可能な限り正確に合計するアルゴリズムについて意見を求めたいと思います。

これが私の最初の解決策です：

put all numbers into a min-absolute-heap. // EDIT as told by comments below
pop the 2 smallest ones.
add them.
put the result back into the heap.
continue until there is only 1 number in the heap.

これは、通常の O(n) の代わりに O(n*logn) を取ります。それは本当に価値がありますか？

2番目の解決策は、私が取り組んでいるデータの特性から来ています。これは、同様の大きさの正数の膨大なリストです。

a[size]; // contains numbers, start at index 0
for(step = 1; step < size; step<<=1)
    for(i = step-1; i+step<size; i+=2*step)
        a[i+step] += a[i];
    if(i < size-1)
        a[size-1] += a[i];

基本的な考え方は、「二分木」方式で合計を行うことです。

注: これは疑似 C コードです。step<<=1ステップを2倍することを意味します。これはO（n）を取ります。もっと良いアプローチがありそうな気がします。推奨/批判できますか?

4

4 に答える 4

21

Kahan の合計アルゴリズムは単純な合計よりもはるかに正確であり、O(n) で実行されます (データ アクセスと比較した浮動小数点の速度に応じて、単純な合計よりも 1 倍から 4 倍遅くなります。デスクトップでは明らかに 4 倍未満の速度です)ハードウェア、およびデータのシャッフルなし)。

または、通常の x86 ハードウェアを使用していて、コンパイラが 80 ビットlong double型へのアクセスを許可している場合は、 type のアキュムレータを使用して単純な合計アルゴリズムを使用しlong doubleます。結果を最後にのみ変換しdoubleます。

本当に多くの精度が必要な場合は、Kahan の総和アルゴリズムで変数 、 、 を使用して、上記の 2 つのソリューションを組み合わせることがlong doubleできcますy。tsum

于 2012-11-16T13:40:32.667 に答える
9

合計の数値誤差を減らすことに関心がある場合は、カハンのアルゴリズムに興味があるかもしれません。

于 2012-11-16T13:40:07.237 に答える
2

私の推測では、バイナリ分解はカハンの合計とほぼ同じように機能すると思います。

これを説明する例を次に示します。

#include <stdio.h>
#include <stdlib.h>
#include <algorithm>

void sumpair( float *a, float *b)
{
    volatile float sum = *a + *b;
    volatile float small = sum - std::max(*a,*b);
    volatile float residue = std::min(*a,*b) - small;
    *a = sum;
    *b = residue;
}

void sumpairs( float *a,size_t size, size_t stride)
{
    if (size <= stride*2 ) {
        if( stride<size )
            sumpair(a+i,a+i+stride);
    } else {
        size_t half = 1;
        while(half*2 < size) half*=2;;
        sumpairs( a , half , stride );
        sumpairs( a+half , size-half , stride );
    }
}

void sumpairwise( float *a,size_t size )
{
    for(size_t stride=1;stride<size;stride*=2)
        sumpairs(a,size,stride);
}

int main()
{
    float data[10000000];
    size_t size= sizeof data/sizeof data[0];
    for(size_t i=0;i<size;i++) data[i]=((1<<30)*-1.0+random())/(1.0+random());

    float naive=0;
    for(size_t i=0;i<size;i++) naive+=data[i];
    printf("naive      sum=%.8g\n",naive);

    double dprec=0;
    for(size_t i=0;i<size;i++) dprec+=data[i];
    printf("dble prec  sum=%.8g\n",(float)dprec);

    sumpairwise( data , size );
    printf("1st approx sum=%.8g\n",data[0]);
    sumpairwise( data+1 , size-1);
    sumpairwise( data , 2 );
    printf("2nd approx sum=%.8g\n",data[0]);
    sumpairwise( data+2 , size-2);
    sumpairwise( data+1 , 2 );
    sumpairwise( data , 2 );
    printf("3rd approx sum=%.8g\n",data[0]);
    return 0;
}

オペランドを volatile と宣言し、x86 アーキテクチャで余分な精度を避けるために -float-store でコンパイルしました

g++  -ffloat-store  -Wl,-stack_size,0x20000000 test_sum.c

そして取得：（0.03125は1ULPです）

naive      sum=-373226.25
dble prec  sum=-373223.03
1st approx sum=-373223
2nd approx sum=-373223.06
3rd approx sum=-373223.06

これには少し説明が必要です。

  • 最初に単純な合計を表示します
  • 次に、倍精度の合計（カハンはそれとほぼ同等です）
  • 最初の近似は、バイナリ分解と同じです。合計を data[0] に保存し、剰余を保存することを除いて。このようにして、合計の前後のデータの正確な合計は変更されません
  • これにより、1 回目の反復を修正するために 2 回目の反復で剰余を合計することで誤差を概算できます (バイナリ和に Kahan を適用するのと同じです)。
  • さらに反復することで、結果をさらに絞り込むことができ、収束が見られます
于 2012-11-18T01:16:41.463 に答える
1

要素は昇順でヒープに配置されるため、代わりに 2 つのキューを使用できます。数値が事前にソートされている場合、これにより O(n) が生成されます。

この疑似コードは、アルゴリズムと同じ結果を生成しO(n)、入力が事前に並べ替えられており、並べ替えアルゴリズムがそれを検出した場合に実行されます。

Queue<float> leaves = sort(arguments[0]).toQueue();
Queue<float> nodes = new Queue();

popAny = #(){
       if(leaves.length == 0) return nodes.pop();
  else if(nodes.length == 0) return leaves.pop();
  else if(leaves.top() > nodes.top()) return nodes.pop();
  else return leaves.pop();
}

while(leaves.length>0 || nodes.length>1) nodes.push(popAny()+popAny());

return nodes.pop();
于 2012-11-16T14:06:42.823 に答える