在python中加速元素數組乘法

Question

我一直在玩弄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) 

Answer 1

如何使用fortran和ctypes？

elementwise.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' 

予編譯使用

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

輸出的文件的Fortran產生的加速〜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 

Answer 2

編輯：從來沒有這個答案，我錯了（見下面的評論）。

---

恐怕在python中比使用numpy更快的矩陣乘法是非常非常困難的。 NumPy通常使用像ATLAS/LAPACK這樣的內部fortran庫，這些庫非常好的優化。

要檢查您的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的評論）。

Answer 3

你最近在做什麼？

隨機數組的創建佔用了整個計算的一部分，如果將​​其包含在您的計算時間內，您幾乎不會在結果中看到任何實際差異，但是，如果您在前面創建它，實際上比較方法。

這是我的結果，我一直在看你在看什麼。 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的，只是compiled it from source： '0.11.0-3-gea20d11髒'

我Fedora中19使用默認numpy的測試此， '1.7.1' 和 numpy的 '1.6.1' 從源代碼編譯，對鏈接：

Update3 我以前的結果當然是不正確的，我在內循環中返回了D，所以跳過了90％的計算。

這爲ali_m的假設提供了更多的證據，證明它比已經非常優化的c代碼真的很難做得更好。

但是，如果您嘗試do something more complicated，例如，，

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

我可以重現的數字傑克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的迄今最優化很難做得更好。

Answer 4

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

gnumpy