错误的numpy平均值？答案

【问题标题】：Wrong numpy mean value?错误的numpy平均值？
【发布时间】：2013-07-04 06:18:38
【问题描述】：

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

【问题讨论】：

如果这些答案中的任何一个解决了您的问题，您应该接受它。

标签： python numpy

【解决方案1】：

这不是 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

(Live demo)

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

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

【讨论】：

谢谢 Oli，我知道这不是 numpy 的问题。我认为有一个函数可以自己拆分累加器以避免这个问题（在 numpy 中实现）
谢谢奥利，我喜欢你的方法。很有用
这将覆盖输入向量p ...不确定这是否可以接受。
@StefanoM：确实，这只是 PoC 代码。将其重写为不合适的工作是微不足道的。
@OliCharlesworth，很抱歉在 PoC 代码上争论，这很愚蠢。但根据我的经验，实现高效且稳定的缩减绝不是微不足道的......

【解决方案2】：

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

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

【讨论】：

是的，问题中已说明。 np.float64 正如你所说的那样解决了这个问题。但是在不改变 dtype 的情况下手动计算平均值时可以解决这个问题。如果你取数据的一小部分并计算部分求和，即使是单精度，你也会得到更好的结果
正确的做法是使用（Welford 的方法）[stackoverflow.com/questions/895929/…，或类似的变体，但在 numpy.制作np.float64 数组的下一个最佳方法是使用dtype 关键字告诉np.mean 使用np.float64 累加器。

【解决方案3】：

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

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

据我所知，numpy 没有模拟。

【讨论】：

+1 的准确性，但在我的机器上比 a.mean() 或 a.mean(axis=-1).mean() 慢 100 倍以上。
确定是的，它是纯python。即使这种事情进入 numpy，与仅仅总结事情相比，还有很多工作要做。但问题当然是这样做是否会在您的真实代码中造成瓶颈——您在原始帖子中提到“有时”:-)。
math.fsum 是用 C 实现的，AS 配方只是一个参考。可能 AS python 代码慢了数千倍......由于 OP 谈到了huge 问题，我虽然速度是一个问题，但在这里我一个人。为了速度和小内存占用而交易准确性并没有错...
当然，你是对的。我的意思是“fsum 会导致 python 开销”。

【解决方案4】：

快速而肮脏的答案

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

要给出更明智的答案，需要更多关于数据结构、数据大小和目标架构的信息。

【讨论】：