Numpy原地操作的性能表现

Question

我正在比较numpy数组的原地操作与常规操作。 这是我所做的(使用Python版本3.7.3)：

    a1, a2 = np.random.random((10,10)), np.random.random((10,10))

为了进行比较：

    def func1(a1, a2):
        a1 = a1 + a2

    def func2(a1, a2):
        a1 += a2

%timeit func1(a1, a2)
%timeit func2(a1, a2)

因为原地操作避免了每次循环的内存分配，所以我本来期望 func1 比 func2 更慢。 但是实际结果是：

In [10]: %timeit func1(a1, a2)
595 ns ± 14.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

In [11]: %timeit func2(a1, a2)
1.38 µs ± 7.87 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

In [12]: np.__version__
Out[12]: '1.16.2'

这表明func1所需的时间仅为func2所需时间的一半。有人可以帮忙解释一下为什么会这样吗？

- jli256

Answer 1

我发现这非常有趣，决定自己测试一下。但是我不仅仅检查 10x10 的数组，还使用了 NumPy 1.16.2 测试了许多不同的数组大小：

这清楚地表明，对于小型数组，普通加法更快，只有在中等大小的数组上，原地操作才更快。在大约 100000 个元素周围还有一个奇怪的突起，我无法解释（它接近我的电脑页面大小，也许使用了不同的分配方案）。

分配临时数组比预期要慢，因为：

• 必须分配那段内存
• 必须迭代三个数组才能执行操作，而不是两个。

特别是第一点（分配内存）可能没有在基准测试中考虑到（不用 %timeit 也不用 simple_benchmark.run）。这是因为反复请求相同大小的内存将是一种可能得到极度优化的东西。这会使具有额外数组的加法似乎比实际上快一些。

还有一点要提到的是，原地加法可能具有更高的常数因子。如果您进行原地加法，则在执行操作之前必须进行更多代码检查，例如对于重叠输入。这可能会给原地加法带来更高的常数因子。

作为更一般的建议：微基准测试可能有所帮助，但它们并不总是真正准确的。您还应该对调用其代码进行基准测试，以更明智地评估代码的实际性能。通常，这样的微基准测试会命中某些高度优化的情况（例如反复分配相同数量的内存并释放它），这在实际使用代码时不会经常发生。

这里也是我用于绘制图形的代码，使用了我的库 simple_benchmark：

from simple_benchmark import BenchmarkBuilder, MultiArgument
import numpy as np

b = BenchmarkBuilder()

@b.add_function()
def func1(a1, a2):
    a1 = a1 + a2

@b.add_function()
def func2(a1, a2):
    a1 += a2

@b.add_arguments('array size')
def argument_provider():
    for exp in range(3, 28):
        dim_size = int(1.4**exp)
        a1 = np.random.random([dim_size, dim_size])
        a2 = np.random.random([dim_size, dim_size])
        yield dim_size ** 2, MultiArgument([a1, a2])

r = b.run()
r.plot()

- MSeifert

---

非常有趣，感谢提供数据和图表！我从未想到原地操作会更慢。但你的观点是原地操作可能需要额外的检查，这很有道理。 - gebbissimo

Answer 2

因为你忽略了向量化操作和预读对小矩阵的影响。

注意到你的矩阵大小只有10x10，因此分配临时存储所需的时间并不重要(目前)，对于具有大缓存大小的处理器，这些小矩阵可能仍然完全适合L1缓存，因此从执行向量化操作等方面获得的速度提升将弥补从分配临时矩阵中损失的时间以及直接添加到已分配内存位置的速度提升。

但是当你增加矩阵的大小时，情况就会变得不同。

In [41]: k = 100

In [42]: a1, a2 = np.random.random((k, k)), np.random.random((k, k))

In [43]: %timeit func2(a1, a2)
4.41 µs ± 3.01 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [44]: %timeit func1(a1, a2)
6.36 µs ± 4.18 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [45]: k = 1000

In [46]: a1, a2 = np.random.random((k, k)), np.random.random((k, k))

In [47]: %timeit func2(a1, a2)
1.13 ms ± 1.49 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [48]: %timeit func1(a1, a2)
1.59 ms ± 2.06 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [49]: k = 5000
In [50]: a1, a2 = np.random.random((k, k)), np.random.random((k, k))

In [51]: %timeit func2(a1, a2)
30.3 ms ± 122 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [52]: %timeit func1(a1, a2)
94.4 ms ± 58.3 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

编辑：这是为了 k = 10 而展示，你观察到的小矩阵的特性也同样适用于我的机器。

In [56]: k = 10

In [57]: a1, a2 = np.random.random((k, k)), np.random.random((k, k))

In [58]: %timeit func2(a1, a2)
1.06 µs ± 10.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

In [59]: %timeit func1(a1, a2)
500 ns ± 0.149 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

- lightalchemist

---

这两个操作对我来说看起来都是矢量化的。就地操作符没有这样做的原因吗？ - kabanus

我并没有说原地操作没有向量化。实际上，如果它们被编译器向量化了，func1和func2中的操作很可能是不同的。我的意思是，人们不应该低估某些向量化操作在数据（包括临时存储）完全适合顶层缓存的性能。当每个级别的缓存大小达到极限时，性能通常会下降，这会导致缓存未命中并需要在缓存之间移动数据，这就是当矩阵的大小变大时我们会遇到的情况。 - lightalchemist

谢谢，这是一个很好的解释。所以基本上你是说额外的分配不是问题。但我仍然不清楚为什么一个会比另一个好（而且好很多！）。你的答案似乎表明时间应该是可比的。 - kabanus

时间不可比较，因为特定的向量化操作可能非常不同。这必须与预取和缓存未命中的影响一起考虑。此外，数组大小的差异也可能导致运行不同的代码（除非我检查Numpy的代码，否则无法确定）。我的观点是，我们不应该仅仅基于我们天真地认为所涉及的操作来期望func1比func2执行得更慢/更快。例如，当Numpy使用像OpenBLAS等多线程BLAS时，情况可能会变得更加复杂。 - lightalchemist

补充一下，这种现象在科学计算领域非常出名。人们可以查阅相关的科学计算、数值算法和分布式计算教材来了解这种效应。在Stackoverflow上也有类似/相关的问题。 - lightalchemist