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

第一个和第二个函数的性能差异并不奇怪。但令我惊讶的是，第三个功能比其他功能慢得多。

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

另外，我什至首先想出第三个函数的原因是因为我担心func1 中向x 分配大量内存的性能。

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

提前致谢。

编辑：

出于好奇，虽然它违背了我创建第 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：

我根据丹尼尔的建议再次修改了代码。

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 与 pure 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），它需要一个包含 16384000 个元素的数组才能使时间翻倍。所以函数 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 绝对是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()

图表是对数对数，我认为这是可视化数据的最佳方式，因为它扩展了几个数量级（我只是希望它仍然可以理解）。

第一个图显示了执行sum 需要多少时间：

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

Answer 2

根据测试（如下所示），您似乎被函数调用开销打败了。除了 NumPy 函数/工具的矢量化功能外，我们还需要为其提供足够的数据来进行处理。对于func3，我们每次调用np.sum 时只给它4 元素。

让我们调查np.sum 的每次调用开销。这是np.sum，从无、一个元素的总和开始 -

In [90]: a = np.array([])

In [91]: %timeit np.sum(a)
1000000 loops, best of 3: 1.6 µs per loop

In [61]: a = np.array([0])

In [62]: %timeit np.sum(a)
1000000 loops, best of 3: 1.66 µs per loop

In [63]: a = np.random.randint(0,9,(100))

In [64]: %timeit np.sum(a)
100000 loops, best of 3: 1.79 µs per loop

In [65]: a = np.random.randint(0,9,(1000))

In [66]: %timeit np.sum(a)
100000 loops, best of 3: 2.25 µs per loop

In [67]: a = np.random.randint(0,9,(10000))

In [68]: %timeit np.sum(a)
100000 loops, best of 3: 7.27 µs per loop

等等。

因此，在这些测试的系统设置中，每次调用 np.sum 至少会产生大约 1.6 u-sec。

让我们看看加法运算符的标量加法是如何执行的-

In [98]: def add_nums(a,b):
    ...:     return a+b
    ...: 

In [99]: %timeit add_nums(2,3)
10000000 loops, best of 3: 71.5 ns per loop

这大约比25x 的每次调用开销快np.sum。

接下来的明显想法是测试func3 在处理更多数据给np.sum 时的表现。

修改func3（使用切片的版本）为每次迭代求和具有可变数据大小：

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

从最初使用的scale_factor = 4 开始 -

In [83]: %timeit func1()
100000 loops, best of 3: 5.39 µs per loop

In [84]: %timeit func2()
10000 loops, best of 3: 39.8 µs per loop

In [85]: %timeit func3(scale_factor = 4)
1000 loops, best of 3: 741 µs per loop

是的，func3 很慢。

现在，让我们在每次调用 np.sum 时提供更多数据，即增加 scale_factor -

In [86]: %timeit func3(scale_factor = 8)
1000 loops, best of 3: 376 µs per loop

In [87]: %timeit func3(scale_factor = 20)
10000 loops, best of 3: 152 µs per loop

In [88]: %timeit func3(scale_factor = 100)
10000 loops, best of 3: 33.5 µs per loop

以此类推，直到我们将整个数据输入到np.sum，以在np.sum 和最小调用开销下实现最大性能限制。

Answer 3

首先，没有人会用 C 编写第三个变体，因为编译器应该进行必要的优化。

所以拿第一个，你有两个 numpy 数组的创建（arange 和 *2）和一个求和。创建像 numpy 数组这样的复杂对象需要一些时间，但是每个向量操作都是用 C 代码编写的，而且速度非常快。

第二个只使用原始的python操作（大约3000，迭代，乘法和求和），用C编写，速度非常快。

您创建了大约 2 * 250 个 numpy 数组的第三个变体（一个相对较慢的操作），与仅创建 2 个 numpy 数组相比，这会导致执行速度慢 100 倍。

如果你对内存使用有顾虑，你应该使用内联操作，它只创建一个数组：

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

如果您仍然需要使用太多内存，请将您的操作分成尽可能大的块。