15

我读过很多“避免使用 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)。我在用着:

  • 蟒蛇3.5.2
  • numpy 1.11.2 + mkl
  • 视窗 10
4

1 回答 1

15

这是正常和预期的行为。应用“避免使用 numpy 的循环”语句太简单了。如果您正在处理内部循环,那么它(几乎)总是正确的。但是在外部循环的情况下(就像你的情况一样),还有更多的例外。特别是如果替代方案是使用广播,因为这会通过使用更多内存来加速您的操作。

只是为“避免使用 numpy 循环”语句添加一些背景知识:

NumPy 数组存储为具有类型的连续数组。Pythonint与 C 不同int!因此,每当您遍历数组中的每个项目时,您都需要从数组中插入该项目,将其转换为 Python int,然后对它进行任何您想做的事情,最后您可能需要再次将其转换为 ac 整数(称为装箱并拆箱值)。例如,您想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

此规则可能有例外,但这些将非常稀少。至少只要有一些等效的 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和 a 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

当您在每次通话中花费 500 美元np.sin时,您不再关心 20 美元的开销。

需要注意的是: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-array 的子数组,请确保每个循环中的函数调用开销与函数中花费的时间相比可以忽略不计。
  • 如果您广播 numpy 数组,请确保您不会耗尽内存。

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

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

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


最后一个想法:

中没有现有的解决方案,则可以使用轻松实现此类需要循环或广播的功能。

例如,将循环解决方案的内存效率与低广播解决方案的速度相结合的 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)')
于 2016-12-28T03:05:17.943 回答