numpy中不同矢量化方法的表现

Question

我想在python中测试矢量化代码的性能:

import timeit
import numpy as np

def func1():
  x = np.arange(1000)
  sum = np.sum(x*2)
  return sum

def func2():
  sum = 0
  for i in xrange(1000):
    sum += i*2
  return sum

def func3():
  sum = 0
  for i in xrange(0,1000,4):
    x = np.arange(i,i+4,1)
    sum += np.sum(x*2)
  return sum

print timeit.timeit(func1, number = 1000)
print timeit.timeit(func2, number = 1000)
print timeit.timeit(func3, number = 1000)

代码提供以下输出:

0.0105729103088
0.069864988327
0.983253955841

第一和第二功能的性能差异并不令人惊讶.但我很惊讶第3个功能明显慢于其他功能.

我在C中的代码中比在Python中更熟悉,第三个函数更像C - 运行for循环并在每个循环中的一条指令中处理4个数字.根据我的理解,numpy调用C函数,然后在C中对代码进行矢量化.因此,如果是这种情况,我的代码也会一次传递4个数字到numpy.当我一次传递更多数字时,代码应该不会更好.那为什么它要慢得多呢？是因为调用numpy函数的开销？

再说,我即使摆在首位的第三个功能上来的原因是因为我担心大量内存分配的性能x中func1.

我的担心有效吗？为什么以及如何改进它或为什么不改进？

提前致谢.

编辑:

出于好奇的缘故,虽然它打破了我创建第3版的最初目的,但我已经研究了roganjosh的建议,并尝试了以下编辑.

def func3():
  sum = 0
  x = np.arange(0,1000)
  for i in xrange(0,1000,4):
    sum += np.sum(x[i:i+4]*2)
  return sum

输出:

0.0104308128357
0.0630609989166
0.748773813248

虽然有所改进,但与其他功能相比仍有很大差距.

是因为x[i:i+4]还是会创建一个新阵列吗？

编辑2:

我根据Daniel的建议再次修改了代码.

def func1():
  x = np.arange(1000)
  x *= 2
  return x.sum()

def func3():
  sum = 0
  x = np.arange(0,1000)
  for i in xrange(0,1000,4):
    x[i:i+4] *= 2
    sum += x[i:i+4].sum()
  return sum

输出:

0.00824999809265
0.0660569667816
0.598328828812

还有另一种加速.所以numpy数组的声明肯定是个问题.现在func3中应该只有一个数组声明,但是时间仍然慢一些.是因为调用numpy数组的开销吗？

Answer 1

看起来你最感兴趣的是你的函数3与纯 NumPy(函数1)和Python(函数2)方法之间的区别.答案很简单(特别是如果你看一下函数4):

• NumPy函数具有"巨大的"常数因子.

您通常需要数千个元素才能进入运行状态,其中运行时np.sum实际上取决于数组中元素的数量.使用IPython和matplotlib(情节在答案的最后),您可以轻松检查运行时依赖性:

import numpy as np

n = []
timing_sum1 = []
timing_sum2 = []
for i in range(1, 25):
    num = 2**i
    arr = np.arange(num)
    print(num)
    time1 = %timeit -o arr.sum()    # calling the method
    time2 = %timeit -o np.sum(arr)  # calling the function
    n.append(num)
    timing_sum1.append(time1)
    timing_sum2.append(time2)

np.sum(缩短)的结果非常有趣:

4
22.6 µs ± 297 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
16
25.1 µs ± 1.08 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
64
25.3 µs ± 1.58 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
256
24.1 µs ± 1.48 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
1024
24.6 µs ± 221 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
4096
27.6 µs ± 147 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
16384
40.6 µs ± 1.29 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
65536
91.2 µs ± 1.03 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
262144
394 µs ± 8.09 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
1048576
1.24 ms ± 4.38 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
4194304
4.71 ms ± 22.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
16777216
18.6 ms ± 280 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

似乎常数因素大致20µs在我的计算机上)并且它需要一个具有16384千个元素的数组来加倍.因此,函数3和4的时序主要是常数因子的时间乘法.

在功能3中,您包括常数因子2次,一次使用np.sum,一次使用np.arange.在这种情况下arange非常便宜,因为每个数组大小相同,因此NumPy和Python以及您的操作系统可能会重用上一次迭代的数组内存.然而,即使这需要时间(大致2µs对于我的计算机上的非常小的阵列).

更一般地说:要识别瓶颈,您应该始终分析功能!

我用line-profiler显示函数的结果.因此我稍微改变了函数,因此它们每行只执行一次操作:

import numpy as np

def func1():
    x = np.arange(1000)
    x = x*2
    return np.sum(x)

def func2():
    sum_ = 0
    for i in range(1000):
        tmp = i*2
        sum_ += tmp
    return sum_

def func3():
    sum_ = 0
    for i in range(0, 1000, 4):  # I'm using python3, so "range" is like "xrange"!
        x = np.arange(i, i + 4, 1)
        x = x * 2
        tmp = np.sum(x)
        sum_ += tmp
    return sum_

def func4():
    sum_ = 0
    x = np.arange(1000)
    for i in range(0, 1000, 4):
        y = x[i:i + 4]
        y = y * 2
        tmp = np.sum(y)
        sum_ += tmp
    return sum_

结果:

%load_ext line_profiler

%lprun -f func1 func1()
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     4                                           def func1():
     5         1           62     62.0     23.8      x = np.arange(1000)
     6         1           65     65.0     24.9      x = x*2
     7         1          134    134.0     51.3      return np.sum(x)

%lprun -f func2 func2()
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     9                                           def func2():
    10         1            7      7.0      0.1      sum_ = 0
    11      1001         2523      2.5     30.9      for i in range(1000):
    12      1000         2819      2.8     34.5          tmp = i*2
    13      1000         2819      2.8     34.5          sum_ += tmp
    14         1            3      3.0      0.0      return sum_

%lprun -f func3 func3()
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    16                                           def func3():
    17         1            7      7.0      0.0      sum_ = 0
    18       251          909      3.6      2.9      for i in range(0, 1000, 4):
    19       250         6527     26.1     21.2          x = np.arange(i, i + 4, 1)
    20       250         5615     22.5     18.2          x = x * 2
    21       250        16053     64.2     52.1          tmp = np.sum(x)
    22       250         1720      6.9      5.6          sum_ += tmp
    23         1            3      3.0      0.0      return sum_

%lprun -f func4 func4()
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    25                                           def func4():
    26         1            7      7.0      0.0      sum_ = 0
    27         1           49     49.0      0.2      x = np.arange(1000)
    28       251          892      3.6      3.4      for i in range(0, 1000, 4):
    29       250         2177      8.7      8.3          y = x[i:i + 4]
    30       250         5431     21.7     20.7          y = y * 2
    31       250        15990     64.0     60.9          tmp = np.sum(y)
    32       250         1686      6.7      6.4          sum_ += tmp
    33         1            3      3.0      0.0      return sum_

我不会去到结果的细节,但你可以看到np.sum是definetly中的瓶颈func3和func4.我已经猜到np.sum是瓶颈之前我写的答案,但这些线成型部实际验证它的瓶颈.

这在使用NumPy时会产生一个非常重要的事实:

• 知道什么时候使用它!小阵列不值得(大多数情况下).
• 了解NumPy功能并使用它们.他们已经使用(如果可用的话)编译器优化标志来展开循环.

如果你真的相信某些部分太慢,那么你可以使用:

• NumPy的C API并使用C处理数组(使用Cython可以非常简单,但您也可以手动执行)
• Numba(基于LLVM).

但一般来说,你可能无法击败NumPy的数量级(数千个条目和更多)阵列.

可视化的时间:

%matplotlib notebook

import matplotlib.pyplot as plt

# Average time per sum-call
fig = plt.figure(1)
ax = plt.subplot(111)
ax.plot(n, [time.average for time in timing_sum1], label='arr.sum()', c='red')
ax.plot(n, [time.average for time in timing_sum2], label='np.sum(arr)', c='blue')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('elements')
ax.set_ylabel('time it takes to sum them [seconds]')
ax.grid(which='both')
ax.legend()

# Average time per element
fig = plt.figure(1)
ax = plt.subplot(111)
ax.plot(n, [time.average / num for num, time in zip(n, timing_sum1)], label='arr.sum()', c='red')
ax.plot(n, [time.average / num for num, time in zip(n, timing_sum2)], label='np.sum(arr)', c='blue')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('elements')
ax.set_ylabel('time per element [seconds / element]')
ax.grid(which='both')
ax.legend()

这些图是log-log,我认为这是可视化数据的最佳方式,因为它扩展了几个数量级(我只希望它仍然可以理解).

第一个图表显示了执行以下操作所需的时间sum:

第二个图显示了执行此操作所需的平均时间sum除以数组中的元素数.这只是解释数据的另一种方式: