使用numexpr提升numpy代码的运行时:分析

你的假设是,加速只来自于/大部分来自于并行化是错误的。正如@brenla已经指出的那样,numexpr加速的最大份额通常来自对缓存的更好利用。但也有一些其他的原因。

首先,numpy和numexpr不计算相同的表达式:

• numpy计算 x**3 和 x**2 作为 pow(x,3) 和 pow(x,2) .
• numexpr有权将其评估为 x**3=x*x*x 和 x**2=x*x .

pow 比一个或两个乘法更复杂,因此更慢,比较:

ne.set_num_threads(1)
%timeit ne.evaluate("0.25*x**3 + 0.75*x**2 + 1.5*x - 2")
# 60.7 ms ± 1.2 ms, base line on my machine

%timeit 0.25*x**3 + 0.75*x**2 + 1.5*x - 2
# 766 ms ± 4.02 ms
%timeit 0.25*x*x*x + 0.75*x*x + 1.5*x - 2 
# 130 ms ± 692 µs 

现在,numexpr的速度只有原来的两倍。我的猜测是 战俘 -版本受CPU限制,而乘法版本受内存限制更大。

当数据较大时,numexpr通常会发光-大于l3缓存(例如,我的计算机上的15MB),在您的示例中给出了,如 x 大约76MB:

• numexp按块计算——即对一个块的所有操作进行计算,每个块至少适合三级缓存,从而最大限度地提高缓存的利用率。只有完成一个块后,才能计算另一个块。
• numpy对整个数据逐一评估一个操作,从而在可以重用数据之前将数据从缓存中移出。

例如,我们可以使用 valgrind (见本帖附录中的脚本):

>>> valgrind --tool=cachegrind python np_version.py
...
...
==5676== D   refs:      1,144,572,370  (754,717,376 rd   + 389,854,994 wr)
==5676== D1  misses:      220,844,716  (181,436,970 rd   +  39,407,746 wr)
==5676== LLd misses:      217,056,340  (178,062,890 rd   +  38,993,450 wr)
==5676== D1  miss rate:          19.3% (       24.0%     +        10.1%  )
==5676== LLd miss rate:          19.0% (       23.6%     +        10.0%  )
....

我们感兴趣的部分是 LLd-misses (即三级未命中,见 here 有关输出解释的信息),大约25%的读取访问是未命中的。

对numexpr的相同分析显示:

>>> valgrind --tool=cachegrind python ne_version.py 
...
==5145== D   refs:      2,612,495,487  (1,737,673,018 rd   + 874,822,469 wr)
==5145== D1  misses:      110,971,378  (   86,949,951 rd   +  24,021,427 wr)
==5145== LLd misses:       29,574,847  (   15,579,163 rd   +  13,995,684 wr)
==5145== D1  miss rate:           4.2% (          5.0%     +         2.7%  )
==5145== LLd miss rate:           1.1% (          0.9%     +         1.6%  )
...

只有5%的读取未命中!

然而,numpy也有一些优势:在引擎盖下,numpy使用mkl例程(至少在我的机器上),而numexpr不使用。因此,numpy最终使用压缩的SSE操作。( movups + mulpd + addpd ,而numexpr最终使用的是标量版本( movsd + mulsd )

这就解释了numpy版本的25%未命中率:一次读取是128位( 莫夫普 )这意味着在4次读取之后,将处理缓存线(64字节),并产生一个未命中。它可以在配置文件中看到(例如 perf 在Linux上):

 32,93 │       movups 0x10(%r15,%rcx,8),%xmm4                                                                               
  1,33 │       movups 0x20(%r15,%rcx,8),%xmm5                                                                               
  1,71 │       movups 0x30(%r15,%rcx,8),%xmm6                                                                               
  0,76 │       movups 0x40(%r15,%rcx,8),%xmm7                                                                               
 24,68 │       movups 0x50(%r15,%rcx,8),%xmm8                                                                               
  1,21 │       movups 0x60(%r15,%rcx,8),%xmm9                                                                               
  2,54 │       movups 0x70(%r15,%rcx,8),%xmm10 

每四分之一 莫夫普 需要更多的时间,因为它等待内存访问。

numpy适合较小的数组大小,适合一级缓存(但随后最大的份额是开销,而不是计算本身,numpy的计算速度更快,但这不起很大作用):

x = np.linspace(-1, 1, 10**3)
%timeit ne.evaluate("0.25*x*x*x + 0.75*x*x + 1.5*x - 2")
# 20.1 µs ± 306 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit 0.25*x*x*x + 0.75*x*x + 1.5*x - 2
# 13.1 µs ± 125 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

附带说明:将函数评估为 ((0.25*x + 0.75)*x + 1.5)*x - 2 .

这都是因为CPU使用率较低:

# small x - CPU bound
x = np.linspace(-1, 1, 10**3)
%timeit ((0.25*x + 0.75)*x + 1.5)*x - 2
#  9.02 µs ± 204 ns 

内存访问也更少:

# large x - memory bound
x = np.linspace(-1, 1, 10**7)
%timeit ((0.25*x + 0.75)*x + 1.5)*x - 2
#  73.8 ms ± 3.71 ms

列表:

一 np_version.py :

import numpy as np

x = np.linspace(-1, 1, 10**7)
for _ in range(10):
    cubic_poly = 0.25*x*x*x + 0.75*x*x + 1.5*x - 2

乙 ne_version.py :

import numpy as np
import numexpr as ne

x = np.linspace(-1, 1, 10**7)
ne.set_num_threads(1)
for _ in range(10):
    ne.evaluate("0.25*x**3 + 0.75*x**2 + 1.5*x - 2")