पाइथन त्वरक (साइथन, नुम्बा, f2py) की तुलना में Numpy einsum

Question

की तुलना में मैं पाइथन त्वरक (Numba, Cython, f2py) की तुलना किसी साधारण समस्या के लिए लूप और नम्पी के इन्सम के लिए सरल (नीचे देखें) की तुलना कर रहा हूं। अब तक इस समस्या के लिए Numpy सबसे तेज़ है (कारक 6x तेज), लेकिन अगर मैं कुछ अतिरिक्त अनुकूलन कर रहा हूं, तो मुझे कुछ प्रतिक्रिया चाहिए, या अगर मैं कुछ गलत कर रहा हूं। यह सरल कोड एक बड़े कोड पर आधारित है जिसमें इनमें से कई ईन्सम कॉल हैं, लेकिन लूप के लिए कोई स्पष्ट नहीं है। मैं जांच रहा हूं कि इनमें से कोई भी त्वरक बेहतर कर सकता है या नहीं।

मैक ओएस एक्स योसामेट पर पाइथन 2.7.9 के साथ किए गए समय, जीसीसी-5.3.0 स्थापित (--with-fortran --without-multilib) होमब्रू से। % टाइमिट कॉल भी किया; ये एकल कॉल समय काफी सटीक हैं।

In [1]: %run -i test_numba.py 
test_numpy: 0.0805640220642 
Matches Numpy output: True 

test_dumb: 1.43043899536 
Matches Numpy output: True 

test_numba: 0.464295864105 
Matches Numpy output: True 

test_cython: 0.627640008926 
Matches Numpy output: True 

test_f2py: 5.01890516281 
Matches Numpy output: True 

test_f2py_order: 2.31424307823 
Matches Numpy output: True 

test_f2py_reorder: 0.507861852646 
Matches Numpy output: True 

मुख्य कोड:

import numpy as np 
import numba 
import time 
import test_f2py as tf2py 
import pyximport 
pyximport.install(setup_args={'include_dirs':np.get_include()}) 
import test_cython as tcyth 

def test_dumb(f,b): 
    fnew = np.empty((f.shape[1],f.shape[2])) 
    for i in range(f.shape[0]): 
     for l in range(f.shape[3]): 
      fnew += f[i,:,:,l] * b[i,l] 
    return fnew 


def test_dumber(f,b): 
    fnew = np.empty((f.shape[1],f.shape[2])) 
    for i in range(f.shape[0]): 
     for j in range(f.shape[1]): 
      for k in range(f.shape[2]): 
       for l in range(f.shape[3]): 
        fnew[j,k] += f[i,j,k,l] * b[i,l] 
    return fnew 

@numba.jit(nopython=True) 
def test_numba(f,b): 
    fnew = np.zeros((f.shape[1],f.shape[2])) #NOTE: can't be empty, gives errors 
    for i in range(f.shape[0]): 
     for j in range(f.shape[1]): 
      for k in range(f.shape[2]): 
       for l in range(f.shape[3]): 
        fnew[j,k] += f[i,j,k,l] * b[i,l] 
    return fnew 

def test_numpy(f,b): 
    return np.einsum('i...k,ik->...',f,b) 

def test_f2py(f,b): 
    return tf2py.test_f2py(f,b) 

def test_f2py_order(f,b): 
    return tf2py.test_f2py(f,b) 

def test_f2py_reorder(f,b): 
    return tf2py.test_f2py_reorder(f,b) 

def test_cython(f,b): 
    return tcyth.test_cython(f,b) 

if __name__ == '__main__': 

    #goal is to create: fnew = sum f*b over dim 0 and 3. 
    f = np.random.rand(32,33,2000,64) 
    b = np.random.rand(32,64) 

    f1 = np.asfortranarray(f) 
    b1 = np.asfortranarray(b) 

    f2 = np.asfortranarray(np.transpose(f,[1,2,0,3])) 

    funcs = [test_dumb,test_numba, test_cython, \ 
      test_f2py,test_f2py_order,test_f2py_reorder] 

    tstart = time.time()  
    fnew_numpy= test_numpy(f,b) 
    tstop = time.time() 
    print test_numpy.__name__+': '+str(tstop-tstart) 
    print 'Matches Numpy output: '+str(np.allclose(fnew_numpy,fnew_numpy)) 
    print '' 

    for func in funcs: 
     tstart = time.time() 
     if func.__name__ == 'test_f2py_order': 
      fnew = func(f1,b1) 
     elif func.__name__ == 'test_f2py_reorder': 
      fnew = func(f2,b1) 
     else: 
      fnew = func(f,b) 
     tstop = time.time() 
     print func.__name__+': '+str(tstop-tstart) 
     print 'Matches Numpy output: '+str(np.allclose(fnew,fnew_numpy)) 
     print '' 

f2py फ़ाइल (साथ f2py -c -m test_f2py test_f2py.F90 संकलित):

!file: test_f2py 
subroutine test_f2py(f,b,fnew,n1,n2,n3,n4) 

integer :: n1,n2,n3,n4 
real(8), dimension(n1,n2,n3,n4) :: f 
real(8), dimension(n1,n4) :: b 
real(8), dimension(n2,n3) :: fnew 
!f2py intent(in) f 
!f2py intent(in) b 
!f2py intent(out) fnew 
!f2py intent(in) n1 
!f2py intent(in) n2 
!f2py intent(in) n3 
!f2py intent(in) n4 

integer :: i1,i2,i3,i4 

do i1=1,n1 
    do i2=1,n2 
     do i3=1,n3 
      do i4=1,n4 
       fnew(i2,i3) = fnew(i2,i3) + f(i1,i2,i3,i4)*b(i1,i4) 
      enddo 
     enddo 
    enddo 
enddo 

end subroutine test_f2py 

subroutine test_f2py_reorder(f,b,fnew,n1,n2,n3,n4) 

integer :: n1,n2,n3,n4 
real(8), dimension(n1,n2,n3,n4) :: f 
real(8), dimension(n3,n4) :: b 
real(8), dimension(n1,n2) :: fnew 
!f2py intent(in) f 
!f2py intent(in) b 
!f2py intent(out) fnew 
!f2py intent(in) n1 
!f2py intent(in) n2 
!f2py intent(in) n3 
!f2py intent(in) n4 

integer :: i1,i2,i3,i4 

do i3=1,n3 
    do i4=1,n4 
     do i1=1,n1 
      do i2=1,n2 
       fnew(i1,i2) = fnew(i1,i2) + f(i1,i2,i3,i4)*b(i3,i4) 
      enddo 
     enddo 
    enddo 
enddo 

end subroutine test_f2py_reorder 

और Cython .pyx फ़ाइल (मुख्य दिनचर्या में pyximport के साथ संकलित):

#/usr/bin python 
import numpy as np 
cimport numpy as np 

def test_cython(np.ndarray[np.float64_t,ndim=4] f, np.ndarray[np.float64_t,ndim=2] b): 
    # cdef np.ndarray[np.float64_t,ndim=4] f 
    # cdef np.ndarray[np.float64_t,ndim=2] b 
    cdef np.ndarray[np.float64_t,ndim=2] fnew = np.empty((f.shape[1],f.shape[2]),dtype=np.float64) 
    cdef int i,j,k,l 
    cdef int Ni = f.shape[0] 
    cdef int Nj = f.shape[1] 
    cdef int Nk = f.shape[2] 
    cdef int Nl = f.shape[3] 

    for i in range(Ni): 
     for j in range(Nj): 
      for k in range(Nk): 
       for l in range(Nl): 
        fnew[j,k] += f[i,j,k,l] * b[i,l] 
    return fnew 

Answer 1

स्ट्रिंग पैरामीटर को पार्स करने के बाद, einsum सभी अक्षों में समेकित उत्पादों की गणना करने के लिए nditer के संकलित संस्करण का उपयोग करता है। स्रोत कोड आसानी से numpy github पर पाया जाता है।

कुछ समय पहले मैंने एक पैच लिखने के हिस्से के रूप में einsum काम-समान काम किया था। इसके हिस्से के रूप में मैंने cython स्क्रिप्ट लिखी जो कि उत्पाद का योग करता है। आप इस कोड देख सकते हैं:

https://github.com/hpaulj/numpy-einsum

मैं einsum गति से चलाने के लिए मेरी कोड प्राप्त करने की कोशिश नहीं की। मैं बस यह समझने की कोशिश कर रहा था कि यह कैसे काम करता है।

Answer 2

आम तौर पर इन त्वरकों का उपयोग पायथन लूप या बहुत से मध्यवर्ती परिणामों के साथ कोड को गति देने के लिए किया जाता है, जबकि einsum पहले से ही बहुत अच्छी तरह अनुकूलित (see source) है। आपको उम्मीद नहीं करनी चाहिए कि वे आसानी से einsum को हराते हैं, लेकिन आप प्रदर्शन में इसके करीब आ सकते हैं।

नुंबा के लिए बेंचमार्क से संकलन समय को बाहर करना महत्वपूर्ण है। यह केवल दो बार फ़ंक्शन को चलाकर पूरा किया जा सकता है (उसी प्रकार के इनपुट के साथ)। जैसे IPython साथ मैं:

f = np.random.rand(32,33,500,64) 
b = np.random.rand(32,64) 

%time _ = test_numba(f,b) # First invocation 
# Wall time: 466 ms 
%time _ = test_numba(f,b) 
# Wall time: 73 ms 
%timeit test_numba(f, b) 
# 10 loops, best of 3: 72.7 ms per loop 
%timeit test_numpy(f, b) 
# 10 loops, best of 3: 62.8 ms per loop 

अपने Cython कोड सुधार के एक नंबर बनाया जा सकता है के लिए: सरणी सीमा और wraparound के लिए

1. अक्षम चेक, compiler directives देखते हैं।
2. निर्दिष्ट करें कि सरणी संगत हैं।
3. typed memoryviews का उपयोग करें।

कुछ की तरह:

cimport cython 
import numpy as np 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def test_cython(double[:,:,:,::1] f, double[:,::1] b): 
    cdef int i, j, k, l, Ni, Nj, Nk, Nl 
    Ni = f.shape[0] 
    Nj = f.shape[1] 
    Nk = f.shape[2] 
    Nl = f.shape[3] 

    fnew = np.empty((Nj, Nk)) 
    cdef double[:,::1] fnew_v = fnew 

    for i in range(Ni): 
     for j in range(Nj): 
      for k in range(Nk): 
       for l in range(Nl): 
        fnew_v[j,k] += f[i,j,k,l] * b[i,l] 
    return fnew 

एक अप-टू-डेट उबंटू 15।10 (x86) यह मुझे einsum के समान गति देता है। हालांकि, एनाकोंडा वितरण के साथ एक ही पीसी पर विंडोज (x86) पर यह साइथन कोड लगभग einsum की आधा गति है। मुझे लगता है कि इसे जीसीसी संस्करणों (5.2.1 बनाम 4.7.0) और एसएसई निर्देशों को सम्मिलित करने की क्षमता (einsum एसएसई 2 इंट्रिनिक्स के साथ कोडित) के साथ करना पड़ सकता है। शायद विभिन्न कंपाइलर विकल्पों की आपूर्ति करने में मदद मिलेगी, लेकिन मुझे यकीन नहीं है।

मुझे शायद ही कोई फोरट्रान पता है, इसलिए मैं उस पर टिप्पणी नहीं कर सकता।

चूंकि आपका लक्ष्य einsum को हरा देना है, मुझे लगता है कि स्पष्ट अगला कदम बढ़ती समांतरता को देखना है। cython.parallel के साथ कुछ धागे को फैलाना काफी आसान होना चाहिए। यदि यह अभी तक आपके सिस्टम मेमोरी बैंडविड्थ को संतृप्त नहीं करता है, तो आप AVX2 और फ़्यूज्ड मल्टीप्ली-एड जैसे नवीनतम CPU निर्देशों को स्पष्ट रूप से शामिल करने का प्रयास कर सकते हैं।

एक और चीज जिसे आप कोशिश कर सकते हैं f को पुन: व्यवस्थित करना और दोबारा बनाना है और np.dot के साथ अपना ऑपरेशन करना है। यदि आपका न्यूम्पी एक अच्छी बीएलएएस लाइब्रेरी के साथ आता है, तो यह सामान्यता के नुकसान की लागत और शायद f सरणी की एक बहुत ही महंगी प्रतिलिपि के बारे में सोचने के लिए हर अनुकूलन को सक्षम कर सकता है।