python - 错误的 numpy 平均值？

Question

我通常使用大型模拟。有时，我需要计算一组粒子的质心。我注意到在很多情况下，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 时不会发生这种情况，因此解决单精度问题应该很好。

Answer 1

这不是 NumPy 问题，而是 float 问题。同样的情况发生在 C:

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

( Live demo )

问题是 float 的精度有限；随着累加器值相对于添加到其中的元素的增长，相对精度下降。

一种解决方案是通过构建加法器树来限制相对增长。这是一个 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

( Live demo )