【问题标题】：numpy ufuncs speed vs for loop speednumpy ufuncs速度与for循环速度
【发布时间】：2016-12-26 00:13:36
【问题描述】：

我读过很多“避免使用 numpy 循环”。所以，我试过了。我正在使用此代码（简化版）。一些辅助数据：

 In[1]: import numpy as np
        resolution = 1000                             # this parameter varies
        tim = np.linspace(-np.pi, np.pi, resolution) 
        prec = np.arange(1, resolution + 1)
        prec = 2 * prec - 1
        values = np.zeros_like(tim)

我的第一个实现是使用for 循环：

 In[2]: for i, ti in enumerate(tim):
           values[i] = np.sum(np.sin(prec * ti))

然后，我摆脱了显式的for 循环，并实现了这一点：

 In[3]: values = np.sum(np.sin(tim[:, np.newaxis] * prec), axis=1)

这个解决方案对于小型阵列来说更快，但是当我扩大规模时，我得到了这样的时间依赖性：

我缺少什么还是正常行为？如果不是，在哪里挖？

编辑：根据 cmets，这里有一些附加信息。时间是用 IPython 的 %timeit 和 %%timeit 测量的，每次运行都是在新内核上执行的。我的笔记本电脑是acer aspire v7-482pg (i7, 8GB)。我正在使用：

python 3.5.2
numpy 1.11.2 + mkl
Windows 10

【问题讨论】：

真的，我正在构建方波，但不是为了污染相关系数，我已经简化了示例。
你有多少内存？如果不够大，tim[:, np.newaxis] * prec 可能需要交换空间，这会导致性能下降。
你是如何对这两个函数进行基准测试的？
@unutbu 我有 8GB。这里尺寸=分辨率*分辨率。是的，我也这么认为，循环次数越多，差异就会越大。但是......这就是我问的原因：）
是的，我在其他 SO 问题中观察到了这一点（我会尝试找到最近的问题）。您的循环案例在大小为resolution 的块上运行。另一个在（大致）resolution**2 的数组上。对于大尺寸，这些大数组的内存管理大致平衡了迭代成本。对于大的resolution，一个块的计算时间相对于迭代开销来说很大。

标签： python performance numpy for-loop numpy-ufunc

【解决方案1】：

这是正常的预期行为。应用 "avoid for loops with numpy" 语句太简单了。如果您正在处理内部循环，那么它（几乎）总是正确的。但是在外部循环的情况下（就像你的情况一样），还有更多的例外。特别是如果替代方案是使用广播，因为这会通过使用更多内存来加快您的操作。

只是为那个 “避免使用 numpy 循环” 语句添加一点背景：

NumPy 数组存储为具有c 类型的连续数组。 Python int 与 C int 不同！因此，每当您遍历数组中的每个项目时，您都需要从数组中插入该项目，将其转换为 Python int，然后对它进行任何您想做的事情，最后您可能需要再次将其转换为 c 整数（称为装箱和拆箱值）。例如，您想使用 Python sum 数组中的项目：

import numpy as np
arr = np.arange(1000)
%%timeit
acc = 0
for item in arr:
    acc += item
# 1000 loops, best of 3: 478 µs per loop

你最好使用 numpy：

%timeit np.sum(arr)
# 10000 loops, best of 3: 24.2 µs per loop

即使你将循环推送到 Python C 代码中，你也离 numpy 的性能还很远：

%timeit sum(arr)
# 1000 loops, best of 3: 387 µs per loop

这条规则可能有例外，但这些情况非常少。至少只要有一些等效的 numpy 功能。因此，如果您要遍历单个元素，那么您应该使用 numpy。

有时一个普通的 python 循环就足够了。它没有被广泛宣传，但与 Python 函数相比，numpy 函数具有巨大的开销。例如考虑一个 3 元素数组：

arr = np.arange(3)
%timeit np.sum(arr)
%timeit sum(arr)

哪个会更快？

解决方案：Python函数比numpy解决方案性能更好：

# 10000 loops, best of 3: 21.9 µs per loop  <- numpy
# 100000 loops, best of 3: 6.27 µs per loop <- python

但这与您的示例有什么关系？实际上并没有那么多，因为您总是在数组上使用 numpy-functions（不是单个元素，甚至不是几个元素），所以您的 inner loop 已经使用了优化的函数。这就是为什么两者的性能大致相同（+/- 10 倍，元素很少，到 2 倍，大约 500 个元素）。但这并不是真正的循环开销，而是函数调用开销！

您的循环解决方案

使用line-profiler 和resolution = 100：

def fun_func(tim, prec, values):
    for i, ti in enumerate(tim):
        values[i] = np.sum(np.sin(prec * ti))
%lprun -f fun_func fun_func(tim, prec, values)
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2       101          752      7.4      5.7      for i, ti in enumerate(tim):
     3       100        12449    124.5     94.3          values[i] = np.sum(np.sin(prec * ti))

95% 都花在了循环中，我什至将循环体分成几个部分来验证这一点：

def fun_func(tim, prec, values):
    for i, ti in enumerate(tim):
        x = prec * ti
        x = np.sin(x)
        x = np.sum(x)
        values[i] = x
%lprun -f fun_func fun_func(tim, prec, values)
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2       101          609      6.0      3.5      for i, ti in enumerate(tim):
     3       100         4521     45.2     26.3          x = prec * ti
     4       100         4646     46.5     27.0          x = np.sin(x)
     5       100         6731     67.3     39.1          x = np.sum(x)
     6       100          714      7.1      4.1          values[i] = x

消费者在这里的时间是np.multiply、np.sin、np.sum，因为您可以通过比较他们每次通话的时间与他们的开销来轻松检查：

arr = np.ones(1, float)
%timeit np.sum(arr)
# 10000 loops, best of 3: 22.6 µs per loop

因此，与计算运行时相比，只要计算函数调用开销很小，您就会有类似的运行时。即使有 100 个项目，您也非常接近开销时间。诀窍是知道他们在什么时候收支平衡。对于 1000 个项目，调用开销仍然很大：

%lprun -f fun_func fun_func(tim, prec, values)
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2      1001         5864      5.9      2.4      for i, ti in enumerate(tim):
     3      1000        42817     42.8     17.2          x = prec * ti
     4      1000       119327    119.3     48.0          x = np.sin(x)
     5      1000        73313     73.3     29.5          x = np.sum(x)
     6      1000         7287      7.3      2.9          values[i] = x

但是与运行时相比，resolution = 5000 的开销非常低：

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2      5001        29412      5.9      0.9      for i, ti in enumerate(tim):
     3      5000       388827     77.8     11.6          x = prec * ti
     4      5000      2442460    488.5     73.2          x = np.sin(x)
     5      5000       441337     88.3     13.2          x = np.sum(x)
     6      5000        36187      7.2      1.1          values[i] = x

当您在每个 np.sin 通话中花费 500us 时，您不再关心 20us 的开销。

请注意：line_profiler 可能包括每行的一些额外开销，也可能包括每个函数调用，因此函数调用开销变得可以忽略不计的点可能会更低！！！

您的广播解决方案

我从分析第一个解决方案开始，让我们对第二个解决方案做同样的事情：

def fun_func(tim, prec, values):
    x = tim[:, np.newaxis]
    x = x * prec
    x = np.sin(x)
    x = np.sum(x, axis=1)
    return x

再次使用带有resolution=100 的line_profiler：

%lprun -f fun_func fun_func(tim, prec, values)
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2         1           27     27.0      0.5      x = tim[:, np.newaxis]
     3         1          638    638.0     12.9      x = x * prec
     4         1         3963   3963.0     79.9      x = np.sin(x)
     5         1          326    326.0      6.6      x = np.sum(x, axis=1)
     6         1            4      4.0      0.1      return x

这已经大大超过了开销时间，因此我们最终比循环快了 10 倍。

我还为resolution=1000 进行了分析：

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2         1           28     28.0      0.0      x = tim[:, np.newaxis]
     3         1        17716  17716.0     14.6      x = x * prec
     4         1        91174  91174.0     75.3      x = np.sin(x)
     5         1        12140  12140.0     10.0      x = np.sum(x, axis=1)
     6         1           10     10.0      0.0      return x

和precision=5000:

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2         1           34     34.0      0.0      x = tim[:, np.newaxis]
     3         1       333685 333685.0     11.1      x = x * prec
     4         1      2391812 2391812.0    79.6      x = np.sin(x)
     5         1       280832 280832.0      9.3      x = np.sum(x, axis=1)
     6         1           14     14.0      0.0      return x

1000 大小仍然更快，但正如我们在那里看到的那样，循环解决方案中的调用开销仍然不可忽略。但是对于resolution = 5000，每一步花费的时间几乎相同（有些慢一些，有些快一些，但总体上非常相似）

另一个影响是，当您进行乘法运算时，实际的广播会变得很重要。即使使用非常智能的 numpy 解决方案，这仍然包括一些额外的计算。对于resolution=10000，您会看到广播乘法开始占用与循环解决方案相关的更多“% 时间”：

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def broadcast_solution(tim, prec, values):
     2         1           37     37.0      0.0      x = tim[:, np.newaxis]
     3         1      1783345 1783345.0    13.9      x = x * prec
     4         1      9879333 9879333.0    77.1      x = np.sin(x)
     5         1      1153789 1153789.0     9.0      x = np.sum(x, axis=1)
     6         1           11     11.0      0.0      return x


Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     8                                           def loop_solution(tim, prec, values):
     9     10001        62502      6.2      0.5      for i, ti in enumerate(tim):
    10     10000      1287698    128.8     10.5          x = prec * ti
    11     10000      9758633    975.9     79.7          x = np.sin(x)
    12     10000      1058995    105.9      8.6          x = np.sum(x)
    13     10000        75760      7.6      0.6          values[i] = x

但是除了实际花费的时间之外还有另一件事：内存消耗。您的循环解决方案需要O(n) 内存，因为您总是处理n 元素。然而，广播解决方案需要O(n*n) 内存。如果您在循环中使用resolution=20000，您可能需要等待一段时间，但它仍然只需要8bytes/element * 20000 element ~= 160kB，但使用广播您将需要~3GB。这忽略了常数因素（如临时数组又名中间数组）！假设你走得更远，你会很快耗尽内存！

是时候再次总结一下要点了：

如果您对 numpy 数组中的单个项目执行 python 循环，那么您做错了。
如果您循环遍历 numpy 数组的子数组，请确保每个循环中的函数调用开销与在函数中花费的时间相比可以忽略不计。
如果您广播 numpy 数组，请确保您不会耗尽内存。

但是关于优化最重要的一点还是：

只有在代码太慢时才优化！如果速度太慢，则仅在分析代码后进行优化。
不要盲目相信简化的语句，并且再次永远不要在没有分析的情况下进行优化。

最后一个想法：

如果numpy 或scipy 中尚无解决方案，则可以使用cython、numba 或numexpr 轻松实现此类需要循环或广播的功能。

例如，一个 numba 函数将循环解决方案的内存效率与低resolutions 的广播解决方案的速度相结合，如下所示：

from numba import njit

import math

@njit
def numba_solution(tim, prec, values):
    size = tim.size
    for i in range(size):
        ti = tim[i]
        x = 0
        for j in range(size):
            x += math.sin(prec[j] * ti)
        values[i] = x

正如 cmets 中指出的，numexpr 也可以非常快速地评估广播的计算并且无需需要O(n*n) 内存：

>>> import numexpr
>>> tim_2d = tim[:, np.newaxis]
>>> numexpr.evaluate('sum(sin(tim_2d * prec), axis=1)')

【讨论】：

这个类比似乎没有用。类比中的情况机制如何与实际问题的机制相对应并不明显，而且很容易产生错误的印象，认为问题是试图在一次调用中完成大量工作，而不是试图使用一个巨大的工作集。
@user2357112 感谢您的反馈！再次删除会更好吗？它是广播与 for 循环的隐喻，如果低估了大小，广播可能会导致内存错误或运行时间过长。这篇文章比我预期的要长得多，所以删除一些次优的部分可能会很好。 :-)
我会说删除它，而是先快速总结一下 Python 外循环如何将问题分块减少内存消耗，同时随着 resolution 的增加而按比例减少开销（因为它是外循环）。
@MSeifert 谢谢你如此详细的回复，有些地方对我来说并不明显，Numba 很棒。我认为我的误解与我的背景有关 - Mathematica。 p.s.：你用什么分析器？
@MSeifert 我发现numexpr.evaluate('sum(sin(prec*ti), axis=1)') 是最简单的加速方法。如果您可以将此添加到您的答案中，那就太好了:)