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

（现场演示）

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

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

您可以np.mean使用dtype关键字参数进行调用，该参数指定累加器的类型（默认为与浮点数组的数组相同的类型）。

因此，调用a.mean(dtype=np.float64)将解决您的玩具示例，也许还有更大数组的问题。

于 2013-07-04T06:29:20.930 回答
3

您可以通过使用内置的 来部分解决此问题math.fsum，它可以追踪部分总和（文档包含指向 AS 配方原型的链接）：

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

据我所知，numpy没有模拟。

于 2013-07-04T08:14:44.553 回答
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 回答