12

私は通常、巨大なシミュレーションを扱っています。場合によっては、一連の粒子の重心を計算する必要があります。多くの状況で、numpy.mean() によって返される平均値が間違っていることに注意しました。アキュムレータの飽和が原因であることがわかります。この問題を回避するために、すべての粒子の合計を小さな粒子のセットに分割することができますが、それは不快です。この問題をエレガントな方法で解決する方法を知っている人はいますか?

好奇心を刺激するためだけに、次の例は、私のシミュレーションで観察したものと同様のものを生成します。

import numpy as np
a = np.ones((1024,1024), dtype=np.float32)*30504.00005

最大値と最小値を確認すると、次のようになります。

a.max() 
30504.0
a.min() 
30504.0

ただし、平均値は次のとおりです。

a.mean()
30687.236328125

ここで何かがおかしいことがわかります。これは dtype=np.float64 を使用する場合には発生しないため、単精度の問題を解決するとよいでしょう。

4

4 に答える 4

7

これは NumPy の問題ではなく、浮動小数点の問題です。Cでも同じことが起こります：

float acc = 0;
for (int i = 0; i < 1024*1024; i++) {
    acc += 30504.00005f;
}
acc /= (1024*1024);
printf("%f\n", acc);  // 30687.304688

(ライブデモ)

問題は、浮動小数点の精度が限られていることです。アキュムレータの値が、追加される要素に対して相対的に大きくなると、相対的な精度が低下します。

解決策の 1 つは、加算器ツリーを構築して相対的な成長を制限することです。これは C の例です (私の Python は十分ではありません...):

float sum(float *p, int n) {
    if (n == 1) return *p;
    for (int i = 0; i < n/2; i++) {
        p[i] += p[i+n/2];
    }
    return sum(p, n/2);
}

float x[1024*1024];
for (int i = 0; i < 1024*1024; i++) {
    x[i] = 30504.00005f;
}

float acc = sum(x, 1024*1024);

acc /= (1024*1024);
printf("%f\n", acc);   // 30504.000000

(ライブデモ)

于 2013-07-04T06:24:23.893 に答える
3

math.fsumこれは、部分和を追跡する組み込みの を使用して部分的に修正できます(ドキュメントには AS レシピ プロトタイプへのリンクが含まれています)。

>>> fsum(a.ravel())/(1024*1024)
30504.0

私の知る限りnumpy、アナログはありません。

于 2013-07-04T08:14:44.553 に答える
3

np.meanアキュムレータの型を指定するキーワード引数を使用して呼び出すことができますdtype(デフォルトは浮動小数点配列の配列と同じ型になります)。

したがって、呼び出すa.mean(dtype=np.float64)ことでおもちゃの例が解決され、おそらく大きな配列の問題が解決されます。

于 2013-07-04T06:29:20.930 に答える
0

迅速で汚い答え

assert a.ndim == 2
a.mean(axis=-1).mean()

これにより、1024*1024 マトリックスの期待される結果が得られますが、もちろん、これはより大きな配列には当てはまりません...

平均の計算がコードのボトルネックにならない場合は、Python でアドホック アルゴリズムを実装します。ただし、詳細はデータ構造によって異なります。

平均の計算がボトルネックである場合は、特殊な (並列) 削減アルゴリズムで問題を解決できます。

編集

このアプローチはばかげているように見えるかもしれませんが、確実に問題を軽減し、.mean()それ自体とほぼ同じくらい効率的です。

In [65]: a = np.ones((1024,1024), dtype=np.float32)*30504.00005

In [66]: a.mean()
Out[66]: 30687.236328125

In [67]: a.mean(axis=-1).mean()
Out[67]: 30504.0

In [68]: %timeit a.mean()
1000 loops, best of 3: 894 us per loop

In [69]: %timeit a.mean(axis=-1).mean()
1000 loops, best of 3: 906 us per loop

より賢明な答えを出すには、データ構造、そのサイズ、およびターゲット アーキテクチャに関する詳細情報が必要です。

于 2013-07-04T06:51:12.977 に答える