14

我一直在玩 numba 和 numexpr 试图加速一个简单的元素矩阵乘法。我一直没能得到更好的结果，它们基本上（速度方面）都相当于 numpys 乘法函数。有没有人在这方面有运气？我是否使用了 numba 和 numexpr 错误（我对此很陌生），或者这完全是一种尝试加快速度的坏方法。这是一个可重现的代码，在此先感谢您：

import numpy as np
from numba import autojit
import numexpr as ne

a=np.random.rand(10,5000000)

# numpy
multiplication1 = np.multiply(a,a)

# numba
def multiplix(X,Y):
    M = X.shape[0]
    N = X.shape[1]
    D = np.empty((M, N), dtype=np.float)
    for i in range(M):
        for j in range(N):
            D[i,j] = X[i, j] * Y[i, j]
    return D

mul = autojit(multiplix)
multiplication2 = mul(a,a)

# numexpr
def numexprmult(X,Y):
    M = X.shape[0]
    N = X.shape[1]
    return ne.evaluate("X * Y")

multiplication3 = numexprmult(a,a)
4

5 回答 5

12

使用和怎么样？

逐元素.F90： 

subroutine elementwise( a, b, c, M, N ) bind(c, name='elementwise')
  use iso_c_binding, only: c_float, c_int

  integer(c_int),intent(in) :: M, N
  real(c_float), intent(in) :: a(M, N), b(M, N)
  real(c_float), intent(out):: c(M, N)

  integer :: i,j

  forall (i=1:M,j=1:N)
    c(i,j) = a(i,j) * b(i,j)
  end forall

end subroutine

elementwise.py：

from ctypes import CDLL, POINTER, c_int, c_float
import numpy as np
import time

fortran = CDLL('./elementwise.so')
fortran.elementwise.argtypes = [ POINTER(c_float), 
                                 POINTER(c_float), 
                                 POINTER(c_float),
                                 POINTER(c_int),
                                 POINTER(c_int) ]

# Setup    
M=10
N=5000000

a = np.empty((M,N), dtype=c_float)
b = np.empty((M,N), dtype=c_float)
c = np.empty((M,N), dtype=c_float)

a[:] = np.random.rand(M,N)
b[:] = np.random.rand(M,N)


# Fortran call
start = time.time()
fortran.elementwise( a.ctypes.data_as(POINTER(c_float)), 
                     b.ctypes.data_as(POINTER(c_float)), 
                     c.ctypes.data_as(POINTER(c_float)), 
                     c_int(M), c_int(N) )
stop = time.time()
print 'Fortran took ',stop - start,'seconds'

# Numpy
start = time.time()
c = np.multiply(a,b)
stop = time.time()
print 'Numpy took ',stop - start,'seconds'

我使用编译了 Fortran 文件

gfortran -O3 -funroll-loops -ffast-math -floop-strip-mine -shared -fPIC \
         -o elementwise.so elementwise.F90

输出产生约 10% 的加速：

 $ python elementwise.py 
Fortran took  0.213667869568 seconds
Numpy took  0.230120897293 seconds
 $ python elementwise.py 
Fortran took  0.209784984589 seconds
Numpy took  0.231616973877 seconds
 $ python elementwise.py 
Fortran took  0.214708089828 seconds
Numpy took  0.25369310379 seconds
于 2013-10-18T20:33:28.040 回答
6

你的时间安排如何？

随机数组的创建占用了计算的全部部分，如果将​​其包含在计时中，您将几乎看不到结果的任何实际差异，但是，如果您预先创建它，您实际上可以比较这些方法。

这是我的结果，我一直在看到你所看到的。numpy 和 numba 给出大致相同的结果（numba 快一点。）

（我没有可用的 numexpr）

In [1]: import numpy as np
In [2]: from numba import autojit
In [3]: a=np.random.rand(10,5000000)

In [4]: %timeit multiplication1 = np.multiply(a,a)
10 loops, best of 3: 90 ms per loop

In [5]: # numba

In [6]: def multiplix(X,Y):
   ...:         M = X.shape[0]
   ...:         N = X.shape[1]
   ...:         D = np.empty((M, N), dtype=np.float)
   ...:         for i in range(M):
   ...:                 for j in range(N):
   ...:                         D[i,j] = X[i, j] * Y[i, j]
   ...:         return D
   ...:         

In [7]: mul = autojit(multiplix)

In [26]: %timeit multiplication1 = np.multiply(a,a)
10 loops, best of 3: 182 ms per loop

In [27]: %timeit multiplication1 = np.multiply(a,a)
10 loops, best of 3: 185 ms per loop

In [28]: %timeit multiplication1 = np.multiply(a,a)
10 loops, best of 3: 181 ms per loop

In [29]: %timeit multiplication2 = mul(a,a)
10 loops, best of 3: 179 ms per loop

In [30]: %timeit multiplication2 = mul(a,a)
10 loops, best of 3: 180 ms per loop

In [31]: %timeit multiplication2 = mul(a,a)
10 loops, best of 3: 178 ms per loop

更新：我使用了最新版本的 numba，只是从源代码编译它：'0.11.0-3-gea20d11-dirty'

我使用从源代码编译的 Fedora 19 中的默认 numpy '1.7.1' 和numpy '1.6.1' 对此进行了测试，链接到：

Update3 我之前的结果当然是不正确的，我在内循环中返回了 D，所以跳过了 90% 的计算。

这为 ali_m 的假设提供了更多证据，即确实很难比已经非常优化的 c 代码做得更好。

但是，如果您尝试做一些更复杂的事情，例如，

np.sqrt(((X[:, None, :] - X) ** 2).sum(-1))

我可以重现 Jake Vanderplas 得到的数字：

In [14]: %timeit pairwise_numba(X)
10000 loops, best of 3: 92.6 us per loop

In [15]: %timeit pairwise_numpy(X)
1000 loops, best of 3: 662 us per loop

所以看起来你正在做一些到目前为止已经被 numpy 优化过的东西，很难做得更好。

于 2013-10-16T14:23:48.267 回答
4

编辑：不要介意这个答案，我错了（见下面的评论）。

恐怕在 python 中比使用 numpy 更快的矩阵乘法会非常非常困难。NumPy 通常使用内部的 fortran 库，例如 ATLAS/LAPACK，它们经过了非常好的优化。

要检查您的 NumPy 版本是否支持 LAPACK：打开终端，转到 Python 安装目录并键入：

for f in `find lib/python2.7/site-packages/numpy/* -name \*.so`; do echo $f; ldd $f;echo "\n";done | grep lapack

请注意，路径可能因您的 python 版本而异。如果你打印了一些行，你肯定有 LAPACK 支持......所以在单核上实现更快的矩阵乘法将很难实现。

现在我不知道使用多核来执行矩阵乘法，所以你可能想研究一下（见 ali_m 的评论）。

于 2013-10-16T09:33:55.450 回答
2

使用 GPU。使用以下软件包。

粗鲁的

于 2013-10-19T09:38:54.930 回答
2

速度np.multiply很大程度上依赖于大小完全相同的数组。

a = np.random.rand(80000,1)
b = np.random.rand(80000,1)

c = np.multiply(a, b)

速度快得要命，而下面的代码需要一分钟多的时间并用完我所有的 16 GB 内存：

a = np.squeeze(np.random.rand(80000,1))
b = np.random.rand(80000,1)

c = np.multiply(a, b)

所以我的建议是使用完全相同尺寸的数组。希望这对寻找如何加速元素乘法的人有用。

于 2020-08-05T14:24:04.923 回答