python - 为什么 numpy ma.average 比 arr.mean 慢 24 倍？

Question

我在 Python 的 numpy 中发现了一些有趣的东西。 ma.average 比 arr.mean 慢很多(arr 是一个数组)

>>> arr = np.full((3, 3), -9999, dtype=float)
array([[-9999., -9999., -9999.],
       [-9999., -9999., -9999.],
       [-9999., -9999., -9999.]])

%timeit np.ma.average(arr, axis=0)
The slowest run took 49.32 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 191 µs per loop

%timeit arr.mean(axis=0)
The slowest run took 6.63 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 7.41 µs per loop

随机数

arr = np.random.random((3,3))

%timeit arr.mean(axis=0)
The slowest run took 6.17 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 7.78 µs per loop

%timeit np.ma.average(arr, axis=0)
1000 loops, best of 3: 186 µs per loop

--> 慢了将近 24 倍。

文档

numpy.ma.average(a, axis=None, weights=None, returned=False)

Return the weighted average of array over the given axis.

numpy.mean(a, axis=None, dtype=None, out=None, keepdims)

Compute the arithmetic mean along the specified axis.

---

为什么 ma.average 比 arr.mean 慢那么多？在数学上它们是相同的(如果我错了请纠正我)。我的猜测是它与 ma.average 上的加权选项有关，但如果没有传递权重，不应该有回退吗？

Answer 1

找出为什么速度变慢的一个好方法是分析它。我将使用第三方库 line_profiler和 IPython 命令 %lprun(参见示例 this blog):

%load_ext line_profiler

import numpy as np
arr = np.full((3, 3), -9999, dtype=float)

%lprun -f np.ma.average np.ma.average(arr, axis=0)

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   519                                           def average(a, axis=None, weights=None, returned=False):
   ...
   570         1         1810   1810.0     30.5      a = asarray(a)
   571         1           15     15.0      0.3      m = getmask(a)
   572                                           
   573                                               # inspired by 'average' in numpy/lib/function_base.py
   574                                           
   575         1            5      5.0      0.1      if weights is None:
   576         1         3500   3500.0     59.0          avg = a.mean(axis)
   577         1          591    591.0     10.0          scl = avg.dtype.type(a.count(axis))
   578                                               else: 
   ...
   608                                           
   609         1            7      7.0      0.1      if returned:
   610                                                   if scl.shape != avg.shape:
   611                                                       scl = np.broadcast_to(scl, avg.shape).copy()
   612                                                   return avg, scl
   613                                               else:
   614         1            5      5.0      0.1          return avg

我删除了一些不相关的行。

所以实际上有 30% 的时间花在 np.ma.asarray 上(arr.mean 不需要这样做!)。

然而，如果您使用更大的数组，相对时间会发生巨大变化:

arr = np.full((1000, 1000), -9999, dtype=float)

%lprun -f np.ma.average np.ma.average(arr, axis=0)
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   519                                           def average(a, axis=None, weights=None, returned=False):
   ...
   570         1          609    609.0      7.6      a = asarray(a)
   571         1           14     14.0      0.2      m = getmask(a)
   572                                           
   573                                               # inspired by 'average' in numpy/lib/function_base.py
   574                                           
   575         1            7      7.0      0.1      if weights is None:
   576         1         6924   6924.0     86.9          avg = a.mean(axis)
   577         1          404    404.0      5.1          scl = avg.dtype.type(a.count(axis))
   578                                               else:
   ...
   609         1            6      6.0      0.1      if returned:
   610                                                   if scl.shape != avg.shape:
   611                                                       scl = np.broadcast_to(scl, avg.shape).copy()
   612                                                   return avg, scl
   613                                               else:
   614         1            6      6.0      0.1          return avg

这次 np.ma.MaskedArray.mean 函数几乎占用了 90% 的时间。

注意:您还可以更深入地研究 np.ma.asarray 或 np.ma.MaskedArray.count 或 np.ma.MaskedArray。 mean 并检查他们的线路配置文件。但我只是想表明有很多被调用的函数会增加开销。

所以下一个问题是:np.ndarray.mean 和 np.ma.average 之间的相对时间是否也发生了变化？至少在我的电脑上，差异现在要小得多:

%timeit np.ma.average(arr, axis=0)
# 2.96 ms ± 91 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit arr.mean(axis=0)
# 1.84 ms ± 23.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

这次它甚至没有慢 2 倍。我假设对于更大的阵列，差异会变得更小。

---

这也是 NumPy 实际上很常见的东西:

即使对于普通的 numpy 函数，常数因子也非常高(例如，参见我对问题 "Performance in different vectorization method in numpy" 的回答)。对于 np.ma，这些常数因子甚至更大，尤其是当您不使用 np.ma.MaskedArray 作为输入时。但是，尽管常数因子可能很高，但这些函数在处理大型数组时表现出色。