numpy ufuncs速度vs循环速度

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）。我使用的是：

• 蟒蛇3.5.2
• numpy的1.11.2 + MKL
• 的Windows 10

Answer 1

这是正常的预期行为。这太简单了，不适用“避免与numpy循环”声明everywere。如果你正在处理内部循环，它几乎总是如此。但是在外部循环的情况下（比如你的情况）有更多的例外。特别是如果替代方案是使用广播，因为这通过使用很多内存加速您的操作。

只是为了一点背景知识添加到“避免与numpy的循环”声明：

NumPy的数组存储为连续阵列，c类型。 Python int与C int不一样！因此，无论何时迭代数组中的每个项目，都需要从数组中插入项目，将其转换为Python int，然后执行任何您想要的操作，最后您可能需要再次将其转换为ac整数（称为装箱和拆箱的价值）。例如，你想用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函数（而不是单个元素，甚至几个元素），所以内部循环已经使用了优化函数。这就是为什么两者执行大致相同（+/-约10倍，只有很少的元素在约500个元素的因子2）。但这不是真正的循环开销，而是函数调用开销！

你的循环液

使用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在这里，你可以轻松地CK通过每次呼叫与他们的开销比较它们的时间：

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

所以只要比你有类似的运行时间计算运行时COMULATIVE函数调用的开销很小。即使有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的阵列做在单品蟒蛇循环。
• 如果循环遍历numpy数组的子阵列，请确保每个循环中的函数调用开销与函数耗用的时间相比可忽略不计。
• 如果您广播numpy阵列，请确保您没有耗尽内存。

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

• 只有优化的代码，如果是太慢了！如果速度太慢，那么只有在分析代码后才进行优化。

• 不要盲目信任简化语句，并再次从不优化没有分析。

---

最终的一个想法：

，要么需要一个回路或广播这样的功能可使用cython，numba或numexpr容易地实现，如果没有已经在numpy或scipy现有的解决方案。

例如，结合了来自在低resolutions广播解决方案的速度循环溶液中的存储器效率是这样的一个numba功能：

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)')