矩陣的優化NumPy的總和通過每一個元素

Question

我使用numpy的1.9，蟒蛇2.7與OpenCV的工作重複，處理大矩陣，我必須進行以下操作多次

def sumShifted(A): # A: numpy array 1000*1000*10 
    return A[:, 0:-1] + A[:, 1:]  

如果可能，我想優化此操作;我嘗試了Cython，但我沒有得到任何顯着的改進，但我不排除這是因爲我的實施不好。

有沒有辦法讓它更快？

編輯： sumShifted獲取調用在一個循環是這樣的：

for i in xrange(0, 400): 
    # ... Various operations on B 
    A = sumShifted(B) 
    # ... Other operations on B 


#More detailed 
for i in xrange(0, 400): 
    A = sumShifted(a11) 
    B = sumShifted(a12) 
    C = sumShifted(b12) 
    D = sumShifted(b22) 

    v = -upQ12/upQ11 

    W, X, Z = self.function1(input_matrix, v, A, C[:,:,4], D[:,:,4]) 
    S, D, F = self.function2(input_matrix, v, A, C[:,:,5], D[:,:,5]) 
    AA  = self.function3(input_matrix, v, A, C[:,:,6], D[:,:,6]) 
    BB  = self.function4(input_matrix, v, A, C[:,:,7], D[:,:,7]) 

EDIT2：按照你的意見，我創建了這個兩個可運行的基準測試（與用Cython）有關合並的4種sumShifted方法在一個。

A, B, C, D= improvedSumShifted(E, F, G, H) 
#E,F: 1000x1000 matrices 
#G,H: 1000x1000x8 matrices 

#first implementation 
def improvedSumShifted(np.ndarray[dtype_t, ndim=2] a, np.ndarray[dtype_t, ndim=2] b, np.ndarray[dtype_t, ndim=3] c, np.ndarray[dtype_t, ndim=3] d): 
    cdef unsigned int i,j,k; 
    cdef unsigned int w = a.shape[0], h = a.shape[1]-1, z = c.shape[2] 
    cdef np.ndarray[dtype_t, ndim=2] aa = np.empty((w, h)) 
    cdef np.ndarray[dtype_t, ndim=2] bb = np.empty((w, h)) 
    cdef np.ndarray[dtype_t, ndim=3] cc = np.empty((w, h, z)) 
    cdef np.ndarray[dtype_t, ndim=3] dd = np.empty((w, h, z)) 
    with cython.boundscheck(False), cython.wraparound(False), cython.overflowcheck(False), cython.nonecheck(False): 
    for i in range(w): 
     for j in range(h): 
     aa[i,j] = a[i,j] + a[i,j+1] 
     bb[i,j] = b[i,j] + b[i,j+1] 
     for k in range(z): 
      cc[i,j,k] = c[i,j,k] + c[i,j+1,k] 
      dd[i,j,k] = d[i,j,k] + d[i,j+1,k] 
return aa, bb, cc, dd 

#second implementation 
def improvedSumShifted(np.ndarray[dtype_t, ndim=2] a, np.ndarray[dtype_t, ndim=2] b, np.ndarray[dtype_t, ndim=3] c, np.ndarray[dtype_t, ndim=3] d): 
    cdef unsigned int i,j,k; 
    cdef unsigned int w = a.shape[0], h = a.shape[1]-1, z = c.shape[2] 
    cdef np.ndarray[dtype_t, ndim=2] aa = np.copy(a[:, 0:h]) 
    cdef np.ndarray[dtype_t, ndim=2] bb = np.copy(b[:, 0:h]) 
    cdef np.ndarray[dtype_t, ndim=3] cc = np.copy(c[:, 0:h]) 
    cdef np.ndarray[dtype_t, ndim=3] dd = np.copy(d[:, 0:h]) 
    with cython.boundscheck(False), cython.wraparound(False), cython.overflowcheck(False), cython.nonecheck(False): 
    for i in range(w): 
    for j in range(h): 
     aa[i,j] += a[i,j+1] 
     bb[i,j] += b[i,j+1] 
     for k in range(z): 
     cc[i,j,k] += c[i,j+1,k] 
     dd[i,j,k] += d[i,j+1,k] 

return aa, bb, cc, dd 

Answer 1

這是不可能的，這個功能可以加快任何進一步的：它確實在蟒蛇的水平只是四種操作：

1. （2次）對輸入片。這些切片速度非常快，因爲它們只需要少量整數運算來計算新的步幅和大小。
2. 爲輸出分配一個新數組。對於這樣一個簡單的功能，這是一個很大的負擔。
3. 評估兩片上的np.add ufunc，這是一種在numpy中高度優化的操作。

事實上，我的基準測試顯示使用numba或cython都沒有改進。在我的機器上，如果預先分配輸出陣列，則每次調用的持續時間爲30毫秒，如果考慮到內存分配，則持續時間爲〜50毫秒。

純numpy的版本：

import numpy as np 

def ss1(A): 
    return np.add(A[:,:-1,:],A[:,1:,:]) 

def ss2(A,output): 
    return np.add(A[:,:-1,:],A[:,1:,:],output) 

的用Cython版本：

import numpy as np 
cimport numpy as np 
cimport cython 

def ss3(np.float64_t[:,:,::1] A not None): 
    cdef unsigned int i,j,k; 
    cdef np.float64_t[:,:,::1] ret = np.empty((A.shape[0],A.shape[1]-1,A.shape[2]),'f8') 
    with cython.boundscheck(False), cython.wraparound(False): 
     for i in range(A.shape[0]): 
      for j in range(A.shape[1]-1): 
       for k in range(A.shape[2]): 
        ret[i,j,k] = A[i,j,k] + A[i,j+1,k] 
    return ret 

def ss4(np.float64_t[:,:,::1] A not None, np.float64_t[:,:,::1] ret not None): 
    cdef unsigned int i,j,k; 
    assert ret.shape[0]>=A.shape[0] and ret.shape[1]>=A.shape[1]-1 and ret.shape[2]>=A.shape[2] 
    with cython.boundscheck(False), cython.wraparound(False): 
     for i in range(A.shape[0]): 
      for j in range(A.shape[1]-1): 
       for k in range(A.shape[2]): 
        ret[i,j,k] = A[i,j,k] + A[i,j+1,k] 
    return ret 

的numba版本（當前numba 0.14.0不能在優化的函數分配新的數組）：

@numba.njit('f8[:,:,:](f8[:,:,:],f8[:,:,:])') 
def ss5(A,output): 
    for i in range(A.shape[0]): 
     for j in range(A.shape[1]-1): 
      for k in range(A.shape[2]): 
       output[i,j,k] = A[i,j,k] + A[i,j+1,k] 
    return output 

以下是計時：

>>> A = np.random.randn((1000,1000,10)) 
>>> output = np.empty((A.shape[0],A.shape[1]-1,A.shape[2])) 

>>> %timeit ss1(A) 
10 loops, best of 3: 50.2 ms per loop 

>>> %timeit ss2(A,output) 
10 loops, best of 3: 30.8 ms per loop 

>>> %timeit ss3(A) 
10 loops, best of 3: 50.8 ms per loop 

>>> %timeit ss4(A,output) 
10 loops, best of 3: 30.9 ms per loop 

>>> %timeit ss5(A,output) 
10 loops, best of 3: 31 ms per loop