numpy 比 numba 和 cython 更快，如何改进 numba 代码
numpy faster than numba and cython , how to improve numba code
我这里有一个简单的例子来帮助我理解使用 numba 和 cython。我是 numba 和 cython 的新手。我已经尽我最大的努力结合所有技巧来使 numba 更快，并且在某种程度上，对于 cython 也是如此，但我的 numpy 代码几乎比 numba 快 2 倍（对于 float64），如果使用 float32 则快 2 倍以上。不确定我在这里遗漏了什么。
我在想问题可能不再是编码，而是更多关于编译器和我不太熟悉的东西。
我浏览了很多关于 numpy、numba 和 cython 的 Whosebug post，但没有找到直接的答案。
numpy 版本：
def py_expsum(x):
return np.sum( np.exp(x) )
numba 版本：
@numba.jit( nopython=True)
def nb_expsum(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in range(ny):
val += np.exp(x[ix, iy])
return val
Cython 版本：
import numpy as np
import cython
from libc.math cimport exp
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef double cy_expsum2 ( double[:,:] x, int nx, int ny ):
cdef:
double val = 0.0
int ix, iy
for ix in range(nx):
for iy in range(ny):
val += exp(x[ix, iy])
return val
使用大小为 2000 x 1000 的数组并循环超过 100 次。对于numba，第一次激活不计入循环。
使用python 3（蟒蛇分布），window 10
float64 / float32
1. numpy : 0.56 sec / 0.23 sec
2. numba : 0.93 sec / 0.74 sec
3. cython: 0.83 sec
cython 接近 numba。所以对我来说最大的问题是为什么 numba 不能打败 numpy 的运行时间？我在这里做错了什么或遗漏了什么？其他因素如何起作用，我如何找出答案？
添加并行化。在 Numba 中，只涉及制作外循环 prange 并将 parallel=True 添加到 jit 选项：
@numba.jit( nopython=True,parallel=True)
def nb_expsum2(x):
nx, ny = x.shape
val = 0.0
for ix in numba.prange(nx):
for iy in range(ny):
val += np.exp( x[ix, iy] )
return val
在我的电脑上，比 non-parallel 版本提速 3.2 倍。也就是说，在我的 PC 上，Numba 和 Cython 都击败了 Numpy。
您也可以 parallelization in Cython - 我没有在这里测试过，但我希望它在性能上与 Numba 相似。 （另请注意，对于 Cython，您可以从 x.shape[0] 和 x.shape[1] 获得 nx 和 ny，因此您不必关闭 bounds-checking 然后完全依赖用户输入保持在范围内）。
正如我们将看到的，行为取决于使用的 numpy-distribution。
此答案将重点关注 Anacoda-distribution 与 Intel 的 VML（矢量数学库），millage 可能因其他硬件和 numpy-version 而异。
还将展示如何通过 Cython 或 numexpr 使用 VML，以防不使用 Anacoda-distribution，其中 plugs-in VML 对某些人来说numpy-operations.
我可以针对以下维度重现您的结果
N,M=2*10**4, 10**3
a=np.random.rand(N, M)
我得到：
%timeit py_expsum(a) # 87ms
%timeit nb_expsum(a) # 672ms
%timeit nb_expsum2(a) # 412ms
calculation-time 的最大份额（约 90%）用于评估 exp- 函数，正如我们将看到的，它是一个 CPU-intensive 任务。
快速浏览一下 top 统计数据显示，numpy 的版本是并行执行的，但 numba 并非如此。但是，在我只有两个处理器的 VM 上，单独的并行化无法解释因子 7 的巨大差异（如 DavidW 的版本 nb_expsum2 所示）。
通过 perf 分析两个版本的代码显示如下：
nb_expsum
Overhead Command Shared Object Symbol
62,56% python libm-2.23.so [.] __ieee754_exp_avx
16,16% python libm-2.23.so [.] __GI___exp
5,25% python perf-28936.map [.] 0x00007f1658d53213
2,21% python mtrand.cpython-37m-x86_64-linux-gnu.so [.] rk_random
py_expsum
31,84% python libmkl_vml_avx.so [.] mkl_vml_kernel_dExp_E9HAynn ▒
9,47% python libiomp5.so [.] _INTERNAL_25_______src_kmp_barrier_cpp_38a91946::__kmp_wait_te▒
6,21% python [unknown] [k] 0xffffffff8140290c ▒
5,27% python mtrand.cpython-37m-x86_64-linux-gnu.so [.] rk_random
正如你所看到的：numpy 在底层使用英特尔的并行矢量化 mkl/vml-version，它很容易胜过 numba 使用的 gnu-math-library (lm.so) 的版本（或并行numba 的版本或 cython 的版本）。可以通过使用并行化稍微平整地面，但 mkl 的矢量化版本仍然优于 numba 和 cython。
但是，仅查看一种尺寸的性能并不是很有启发性，在 exp 的情况下（对于其他先验函数），有两个维度需要考虑：
- 数组中的元素数量 - 缓存效果和不同大小的不同算法（在 numpy 中并非闻所未闻）会导致不同的性能。
- 根据
x值，计算exp(x)需要不同的时间。通常有三种不同类型的输入导致不同的计算时间：非常小、正常和非常大（non-finite 结果）
我正在使用 perfplot 来可视化结果（请参阅附录中的代码）。对于 "normal" 范围，我们得到以下性能：
虽然 0.0 的性能相似，但我们可以看到，一旦结果变为无穷大，Intel 的 VML 就会产生相当大的负面影响：
但是还有其他需要注意的地方：
- 对于向量大小
<= 8192 = 2^13，numpy 使用 non-parallelized glibc-version of exp（同样使用 numba 和 cython）。
我使用 - Anaconda-distribution，overrides numpy's functionality and plugs Intel's VML-library 用于大小 > 8192，它是矢量化和并行化的 - 这解释了大小约为 10^4 时 运行 时间的下降。
- numba 对于较小的尺寸来说很容易击败通常的 glibc-version（对于 numpy 来说开销太大），但是对于更大的数组来说（如果 numpy 不切换到 VML）差别不大。
- 这似乎是一项 CPU-bound 任务 - 我们在任何地方都看不到 cache-boundaries。
- Parallized numba-version 只有在元素超过 500 个时才有意义。
那么后果是什么？
- 如果不超过8192个元素，应该使用numba-version。
- 否则numpy版本（即使没有VML-plugin也不会损失多少）
注意：numba 无法自动使用英特尔 VML 中的 vdExp（如评论中部分建议的那样），因为它单独计算 exp(x)，而 VML 在整个数组上运行。
可以减少写入和加载数据时的缓存未命中，这是由 numpy-version 使用以下算法执行的：
- 对适合缓存但也不太小（开销）的部分数据执行 VML
vdExp。
- 总结生成的工作数组。
- 执行1.+2。对于下一部分数据，直到处理完整个数据。
但是，与 numpy 的版本相比，我不希望获得超过 10% 的收益（但也许我错了），因为无论如何 90% 的计算时间都花在了 MVL 上。
尽管如此，这里是 Cython 中可能的快速实现：
%%cython -L=<path_mkl_libs> --link-args=-Wl,-rpath=<path_mkl_libs> --link-args=-Wl,--no-as-needed -l=mkl_intel_ilp64 -l=mkl_core -l=mkl_gnu_thread -l=iomp5
# path to mkl can be found via np.show_config()
# which libraries needed: https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor
# another option would be to wrap mkl.h:
cdef extern from *:
"""
// MKL_INT is 64bit integer for mkl-ilp64
// see https://software.intel.com/en-us/mkl-developer-reference-c-c-datatypes-specific-to-intel-mkl
#define MKL_INT long long int
void vdExp(MKL_INT n, const double *x, double *y);
"""
void vdExp(long long int n, const double *x, double *y)
def cy_expsum(const double[:,:] v):
cdef:
double[1024] w;
int n = v.size
int current = 0;
double res = 0.0
int size = 0
int i = 0
while current<n:
size = n-current
if size>1024:
size = 1024
vdExp(size, &v[0,0]+current, w)
for i in range(size):
res+=w[i]
current+=size
return res
然而，正是 numexpr 会做的，它也使用 Intel 的 vml 作为后端：
import numexpr as ne
def ne_expsum(x):
return ne.evaluate("sum(exp(x))")
至于时间我们可以看到如下：
具有以下值得注意的细节：
- numpy、numexpr 和 cython 版本对于更大的数组具有几乎相同的性能——这并不奇怪，因为它们使用相同的 vml-functionality.
- 在这三个中，cython-version 开销最小，numexpr 开销最大
- numexpr-version 可能是最容易编写的（假设不是每个 numpy 分发插件 mvl-functionality）。
清单：
地块：
import numpy as np
def py_expsum(x):
return np.sum(np.exp(x))
import numba as nb
@nb.jit( nopython=True)
def nb_expsum(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in range(ny):
val += np.exp( x[ix, iy] )
return val
@nb.jit( nopython=True, parallel=True)
def nb_expsum2(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in nb.prange(ny):
val += np.exp( x[ix, iy] )
return val
import perfplot
factor = 1.0 # 0.0 or 1e4
perfplot.show(
setup=lambda n: factor*np.random.rand(1,n),
n_range=[2**k for k in range(0,27)],
kernels=[
py_expsum,
nb_expsum,
nb_expsum2,
],
logx=True,
logy=True,
xlabel='len(x)'
)
这取决于 exp 实现和并行化
如果您在 Numpy 中使用英特尔 SVML，那么也可以在 Numba、Numexpr 或 Cython 等其他软件包中使用它。 Numba performance tips
如果 Numpy 命令被并行化，也尝试在 Numba 或 Cython 中对其进行并行化。
代码
import os
#Have to be before importing numpy
#Test with 1 Thread against a single thread Numba/Cython Version and
#at least with number of physical cores against parallel versions
os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
#from version 0.43 until 0.47 this has to be set before importing numba
#Bug: https://github.com/numba/numba/issues/4689
from llvmlite import binding
binding.set_option('SVML', '-vector-library=SVML')
import numba as nb
def py_expsum(x):
return np.sum( np.exp(x) )
@nb.njit(parallel=False,fastmath=True) #set it to True for a parallel version
def nb_expsum(x):
val = nb.float32(0.)#change this to float64 on the float64 version
for ix in nb.prange(x.shape[0]):
for iy in range(x.shape[1]):
val += np.exp(x[ix,iy])
return val
N,M=2000, 1000
#a=np.random.rand(N*M).reshape((N,M)).astype(np.float32)
a=np.random.rand(N*M).reshape((N,M))
基准测试
#float64
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "1"
#7.44 ms ± 86.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "6"
#4.83 ms ± 139 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=false
#2.49 ms ± 25.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) ##parallel=true
#568 µs ± 45.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
#float32
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "1"
#3.44 ms ± 66.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "6"
#2.59 ms ± 35.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=false
#1 ms ± 12.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit nb_expsum(a) #parallel=true
#252 µs ± 19.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
使用 SVML 的 Perfplot
import numpy as np
#from version 0.43 until 0.47 this has to be set before importing numba
#Bug: https://github.com/numba/numba/issues/4689
from llvmlite import binding
binding.set_option('SVML', '-vector-library=SVML')
import numba as nb
def py_expsum(x):
return np.sum(np.exp(x))
@nb.jit( nopython=True,parallel=False,fastmath=False)
def nb_expsum_single_thread(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in range(ny):
val += np.exp( x[ix, iy] )
return val
#fastmath makes SIMD-vectorization possible
#val+=some_value is not vectorizable (scalar depends on scalar)
#This would also prevents the usage of SVML
@nb.jit( nopython=True,parallel=False,fastmath=True)
def nb_expsum_single_thread_vec(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in range(ny):
val += np.exp( x[ix, iy] )
return val
@nb.jit(nopython=True,parallel=True,fastmath=False)
def nb_expsum_parallel(x):
nx, ny = x.shape
val = 0.0
#parallelization over the outer loop is almost every time faster
#except for rare cases like this (x.shape -> (1,n))
for ix in range(nx):
for iy in nb.prange(ny):
val += np.exp( x[ix, iy] )
return val
#fastmath makes SIMD-vectorization possible
#val+=some_value is not vectorizable (scalar depends on scalar)
#This would also prevents the usage of SVML
@nb.jit(nopython=True,parallel=True,fastmath=True)
def nb_expsum_parallel_vec(x):
nx, ny = x.shape
val = 0.0
#parallelization over the outer loop is almost every time faster
#except for rare cases like this (x.shape -> (1,n))
for ix in range(nx):
for iy in nb.prange(ny):
val += np.exp( x[ix, iy] )
return val
import perfplot
factor = 1.0 # 0.0 or 1e4
perfplot.show(
setup=lambda n: factor*np.random.rand(1,n),
n_range=[2**k for k in range(0,27)],
kernels=[
py_expsum,
nb_expsum_single_thread,
nb_expsum_single_thread_vec,
nb_expsum_parallel,
nb_expsum_parallel_vec,
cy_expsum
],
logx=True,
logy=True,
xlabel='len(x)'
)
检查是否使用了SVML
可用于检查一切是否按预期工作。
def check_SVML(func):
if 'intel_svmlcc' in func.inspect_llvm(func.signatures[0]):
print("found")
else:
print("not found")
check_SVML(nb_expsum_parallel_vec)
#found
我这里有一个简单的例子来帮助我理解使用 numba 和 cython。我是 numba 和 cython 的新手。我已经尽我最大的努力结合所有技巧来使 numba 更快，并且在某种程度上，对于 cython 也是如此，但我的 numpy 代码几乎比 numba 快 2 倍（对于 float64），如果使用 float32 则快 2 倍以上。不确定我在这里遗漏了什么。
我在想问题可能不再是编码，而是更多关于编译器和我不太熟悉的东西。
我浏览了很多关于 numpy、numba 和 cython 的 Whosebug post，但没有找到直接的答案。
numpy 版本：
def py_expsum(x):
return np.sum( np.exp(x) )
numba 版本：
@numba.jit( nopython=True)
def nb_expsum(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in range(ny):
val += np.exp(x[ix, iy])
return val
Cython 版本：
import numpy as np
import cython
from libc.math cimport exp
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef double cy_expsum2 ( double[:,:] x, int nx, int ny ):
cdef:
double val = 0.0
int ix, iy
for ix in range(nx):
for iy in range(ny):
val += exp(x[ix, iy])
return val
使用大小为 2000 x 1000 的数组并循环超过 100 次。对于numba，第一次激活不计入循环。
使用python 3（蟒蛇分布），window 10
float64 / float32
1. numpy : 0.56 sec / 0.23 sec
2. numba : 0.93 sec / 0.74 sec
3. cython: 0.83 sec
cython 接近 numba。所以对我来说最大的问题是为什么 numba 不能打败 numpy 的运行时间？我在这里做错了什么或遗漏了什么？其他因素如何起作用，我如何找出答案？
添加并行化。在 Numba 中，只涉及制作外循环 prange 并将 parallel=True 添加到 jit 选项：
@numba.jit( nopython=True,parallel=True)
def nb_expsum2(x):
nx, ny = x.shape
val = 0.0
for ix in numba.prange(nx):
for iy in range(ny):
val += np.exp( x[ix, iy] )
return val
在我的电脑上，比 non-parallel 版本提速 3.2 倍。也就是说，在我的 PC 上，Numba 和 Cython 都击败了 Numpy。
您也可以 parallelization in Cython - 我没有在这里测试过，但我希望它在性能上与 Numba 相似。 （另请注意，对于 Cython，您可以从 x.shape[0] 和 x.shape[1] 获得 nx 和 ny，因此您不必关闭 bounds-checking 然后完全依赖用户输入保持在范围内）。
正如我们将看到的，行为取决于使用的 numpy-distribution。
此答案将重点关注 Anacoda-distribution 与 Intel 的 VML（矢量数学库），millage 可能因其他硬件和 numpy-version 而异。
还将展示如何通过 Cython 或 numexpr 使用 VML，以防不使用 Anacoda-distribution，其中 plugs-in VML 对某些人来说numpy-operations.
我可以针对以下维度重现您的结果
N,M=2*10**4, 10**3
a=np.random.rand(N, M)
我得到：
%timeit py_expsum(a) # 87ms
%timeit nb_expsum(a) # 672ms
%timeit nb_expsum2(a) # 412ms
calculation-time 的最大份额（约 90%）用于评估 exp- 函数，正如我们将看到的，它是一个 CPU-intensive 任务。
快速浏览一下 top 统计数据显示，numpy 的版本是并行执行的，但 numba 并非如此。但是，在我只有两个处理器的 VM 上，单独的并行化无法解释因子 7 的巨大差异（如 DavidW 的版本 nb_expsum2 所示）。
通过 perf 分析两个版本的代码显示如下：
nb_expsum
Overhead Command Shared Object Symbol
62,56% python libm-2.23.so [.] __ieee754_exp_avx
16,16% python libm-2.23.so [.] __GI___exp
5,25% python perf-28936.map [.] 0x00007f1658d53213
2,21% python mtrand.cpython-37m-x86_64-linux-gnu.so [.] rk_random
py_expsum
31,84% python libmkl_vml_avx.so [.] mkl_vml_kernel_dExp_E9HAynn ▒
9,47% python libiomp5.so [.] _INTERNAL_25_______src_kmp_barrier_cpp_38a91946::__kmp_wait_te▒
6,21% python [unknown] [k] 0xffffffff8140290c ▒
5,27% python mtrand.cpython-37m-x86_64-linux-gnu.so [.] rk_random
正如你所看到的：numpy 在底层使用英特尔的并行矢量化 mkl/vml-version，它很容易胜过 numba 使用的 gnu-math-library (lm.so) 的版本（或并行numba 的版本或 cython 的版本）。可以通过使用并行化稍微平整地面，但 mkl 的矢量化版本仍然优于 numba 和 cython。
但是，仅查看一种尺寸的性能并不是很有启发性，在 exp 的情况下（对于其他先验函数），有两个维度需要考虑：
- 数组中的元素数量 - 缓存效果和不同大小的不同算法（在 numpy 中并非闻所未闻）会导致不同的性能。
- 根据
x值，计算exp(x)需要不同的时间。通常有三种不同类型的输入导致不同的计算时间：非常小、正常和非常大（non-finite 结果）
我正在使用 perfplot 来可视化结果（请参阅附录中的代码）。对于 "normal" 范围，我们得到以下性能：
虽然 0.0 的性能相似，但我们可以看到，一旦结果变为无穷大，Intel 的 VML 就会产生相当大的负面影响：
但是还有其他需要注意的地方：
- 对于向量大小
<= 8192 = 2^13，numpy 使用 non-parallelized glibc-version of exp（同样使用 numba 和 cython）。
我使用 - Anaconda-distribution，overrides numpy's functionality and plugs Intel's VML-library 用于大小 > 8192，它是矢量化和并行化的 - 这解释了大小约为 10^4 时 运行 时间的下降。
- numba 对于较小的尺寸来说很容易击败通常的 glibc-version（对于 numpy 来说开销太大），但是对于更大的数组来说（如果 numpy 不切换到 VML）差别不大。
- 这似乎是一项 CPU-bound 任务 - 我们在任何地方都看不到 cache-boundaries。
- Parallized numba-version 只有在元素超过 500 个时才有意义。
那么后果是什么？
- 如果不超过8192个元素，应该使用numba-version。
- 否则numpy版本（即使没有VML-plugin也不会损失多少）
注意：numba 无法自动使用英特尔 VML 中的 vdExp（如评论中部分建议的那样），因为它单独计算 exp(x)，而 VML 在整个数组上运行。
可以减少写入和加载数据时的缓存未命中，这是由 numpy-version 使用以下算法执行的：
- 对适合缓存但也不太小（开销）的部分数据执行 VML
vdExp。 - 总结生成的工作数组。
- 执行1.+2。对于下一部分数据，直到处理完整个数据。
但是，与 numpy 的版本相比，我不希望获得超过 10% 的收益（但也许我错了），因为无论如何 90% 的计算时间都花在了 MVL 上。
尽管如此，这里是 Cython 中可能的快速实现：
%%cython -L=<path_mkl_libs> --link-args=-Wl,-rpath=<path_mkl_libs> --link-args=-Wl,--no-as-needed -l=mkl_intel_ilp64 -l=mkl_core -l=mkl_gnu_thread -l=iomp5
# path to mkl can be found via np.show_config()
# which libraries needed: https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor
# another option would be to wrap mkl.h:
cdef extern from *:
"""
// MKL_INT is 64bit integer for mkl-ilp64
// see https://software.intel.com/en-us/mkl-developer-reference-c-c-datatypes-specific-to-intel-mkl
#define MKL_INT long long int
void vdExp(MKL_INT n, const double *x, double *y);
"""
void vdExp(long long int n, const double *x, double *y)
def cy_expsum(const double[:,:] v):
cdef:
double[1024] w;
int n = v.size
int current = 0;
double res = 0.0
int size = 0
int i = 0
while current<n:
size = n-current
if size>1024:
size = 1024
vdExp(size, &v[0,0]+current, w)
for i in range(size):
res+=w[i]
current+=size
return res
然而，正是 numexpr 会做的，它也使用 Intel 的 vml 作为后端：
import numexpr as ne
def ne_expsum(x):
return ne.evaluate("sum(exp(x))")
至于时间我们可以看到如下：
具有以下值得注意的细节：
- numpy、numexpr 和 cython 版本对于更大的数组具有几乎相同的性能——这并不奇怪，因为它们使用相同的 vml-functionality.
- 在这三个中，cython-version 开销最小，numexpr 开销最大
- numexpr-version 可能是最容易编写的（假设不是每个 numpy 分发插件 mvl-functionality）。
清单：
地块：
import numpy as np
def py_expsum(x):
return np.sum(np.exp(x))
import numba as nb
@nb.jit( nopython=True)
def nb_expsum(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in range(ny):
val += np.exp( x[ix, iy] )
return val
@nb.jit( nopython=True, parallel=True)
def nb_expsum2(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in nb.prange(ny):
val += np.exp( x[ix, iy] )
return val
import perfplot
factor = 1.0 # 0.0 or 1e4
perfplot.show(
setup=lambda n: factor*np.random.rand(1,n),
n_range=[2**k for k in range(0,27)],
kernels=[
py_expsum,
nb_expsum,
nb_expsum2,
],
logx=True,
logy=True,
xlabel='len(x)'
)
这取决于 exp 实现和并行化
如果您在 Numpy 中使用英特尔 SVML，那么也可以在 Numba、Numexpr 或 Cython 等其他软件包中使用它。 Numba performance tips
如果 Numpy 命令被并行化，也尝试在 Numba 或 Cython 中对其进行并行化。
代码
import os
#Have to be before importing numpy
#Test with 1 Thread against a single thread Numba/Cython Version and
#at least with number of physical cores against parallel versions
os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
#from version 0.43 until 0.47 this has to be set before importing numba
#Bug: https://github.com/numba/numba/issues/4689
from llvmlite import binding
binding.set_option('SVML', '-vector-library=SVML')
import numba as nb
def py_expsum(x):
return np.sum( np.exp(x) )
@nb.njit(parallel=False,fastmath=True) #set it to True for a parallel version
def nb_expsum(x):
val = nb.float32(0.)#change this to float64 on the float64 version
for ix in nb.prange(x.shape[0]):
for iy in range(x.shape[1]):
val += np.exp(x[ix,iy])
return val
N,M=2000, 1000
#a=np.random.rand(N*M).reshape((N,M)).astype(np.float32)
a=np.random.rand(N*M).reshape((N,M))
基准测试
#float64
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "1"
#7.44 ms ± 86.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "6"
#4.83 ms ± 139 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=false
#2.49 ms ± 25.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) ##parallel=true
#568 µs ± 45.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
#float32
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "1"
#3.44 ms ± 66.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "6"
#2.59 ms ± 35.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=false
#1 ms ± 12.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit nb_expsum(a) #parallel=true
#252 µs ± 19.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
使用 SVML 的 Perfplot
import numpy as np
#from version 0.43 until 0.47 this has to be set before importing numba
#Bug: https://github.com/numba/numba/issues/4689
from llvmlite import binding
binding.set_option('SVML', '-vector-library=SVML')
import numba as nb
def py_expsum(x):
return np.sum(np.exp(x))
@nb.jit( nopython=True,parallel=False,fastmath=False)
def nb_expsum_single_thread(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in range(ny):
val += np.exp( x[ix, iy] )
return val
#fastmath makes SIMD-vectorization possible
#val+=some_value is not vectorizable (scalar depends on scalar)
#This would also prevents the usage of SVML
@nb.jit( nopython=True,parallel=False,fastmath=True)
def nb_expsum_single_thread_vec(x):
nx, ny = x.shape
val = 0.0
for ix in range(nx):
for iy in range(ny):
val += np.exp( x[ix, iy] )
return val
@nb.jit(nopython=True,parallel=True,fastmath=False)
def nb_expsum_parallel(x):
nx, ny = x.shape
val = 0.0
#parallelization over the outer loop is almost every time faster
#except for rare cases like this (x.shape -> (1,n))
for ix in range(nx):
for iy in nb.prange(ny):
val += np.exp( x[ix, iy] )
return val
#fastmath makes SIMD-vectorization possible
#val+=some_value is not vectorizable (scalar depends on scalar)
#This would also prevents the usage of SVML
@nb.jit(nopython=True,parallel=True,fastmath=True)
def nb_expsum_parallel_vec(x):
nx, ny = x.shape
val = 0.0
#parallelization over the outer loop is almost every time faster
#except for rare cases like this (x.shape -> (1,n))
for ix in range(nx):
for iy in nb.prange(ny):
val += np.exp( x[ix, iy] )
return val
import perfplot
factor = 1.0 # 0.0 or 1e4
perfplot.show(
setup=lambda n: factor*np.random.rand(1,n),
n_range=[2**k for k in range(0,27)],
kernels=[
py_expsum,
nb_expsum_single_thread,
nb_expsum_single_thread_vec,
nb_expsum_parallel,
nb_expsum_parallel_vec,
cy_expsum
],
logx=True,
logy=True,
xlabel='len(x)'
)
检查是否使用了SVML
可用于检查一切是否按预期工作。
def check_SVML(func):
if 'intel_svmlcc' in func.inspect_llvm(func.signatures[0]):
print("found")
else:
print("not found")
check_SVML(nb_expsum_parallel_vec)
#found