python - numpy ufuncs 速度与 for 循环速度

Question

我读了很多“避免使用 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)

对于小型阵列，这个解决方案更快，但是当我放大时，我得到了这样的时间依赖性:


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

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

python 3.5.2

numpy 1.11.2 + mkl

Windows 10

Answer 1

这是正常和预期的行为。应用 太简单了“避免使用 numpy 进行循环”声明每一个。如果您正在处理内部循环，它(几乎)总是正确的。但是在外循环的情况下(就像你的情况)有更多的异常(exception)。特别是如果替代方案是使用广播，因为这可以通过使用 加快您的操作速度。更多内存。

只是为了添加一点背景 “避免使用 numpy 进行循环”陈述:

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

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

此规则可能有异常(exception)，但这些将非常稀少。至少只要有一些等效的 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 函数(不是单个元素，甚至不是几个元素)，所以你的 内循环已经使用了优化的功能。这就是为什么两者的性能大致相同(+/- 因子 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

再次使用 line_profiler 和 resolution=100 :

%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 数组，请确保您不会耗尽内存。

但关于优化最重要的一点仍然是:

只有在代码太慢时才优化代码!如果速度太慢，则仅在分析代码后进行优化。

不要盲目相信简化的语句，也不要在没有分析的情况下进行优化。

最后一个想法:

使用 cython 可以轻松实现此类需要循环或广播的功能。 , numba或 numexpr如果 numpy 中没有已经存在的解决方案或 scipy .

例如，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

正如评论中指出的 numexpr还可以非常快速地评估广播计算和 没有 要求 O(n*n)内存:

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