Co powoduje spowolnienie 2x w mojej implementacji Cython w mnożeniu wektorów macierzy?

Question

Obecnie próbuję wdrożyć podstawowy wektor mnożenie macierzy w Cython (jako część większego larger project to reduce computation) i stwierdzeniu, że mój kod jest około 2x wolniej niż Numpy.dot.

Zastanawiam się, czy jest coś, co mi brakuje, że jest w wyniku spowolnienia. Piszę zoptymalizowany kod Cythona, deklarując typy zmiennych, wymagające sąsiadujących tablic i unikając chybionych pamięci podręcznych. Próbowałem nawet mieć Cythona jako opakowanie i wywoływanie natywnego kodu C (zobacz poniżej).

Zastanawiam się: co jeszcze mogę zrobić, aby przyspieszyć moją implementację, tak działa tak szybko jak NumPy dla tej podstawowej operacji?

---

Kod Cython że używam jest beow:

import numpy as np 
cimport numpy as np 
cimport cython 

DTYPE = np.float64; 
ctypedef np.float64_t DTYPE_T 

@cython.boundscheck(False) 
@cython.wraparound(False) 
@cython.nonecheck(False) 
def matrix_vector_multiplication(np.ndarray[DTYPE_T, ndim=2] A, np.ndarray[DTYPE_T, ndim=1] x): 

    cdef Py_ssize_t i, j 
    cdef Py_ssize_t N = A.shape[0] 
    cdef Py_ssize_t D = A.shape[1] 
    cdef np.ndarray[DTYPE_T, ndim=1] y = np.empty(N, dtype = DTYPE) 
    cdef DTYPE_T val 

    for i in range(N): 
     val = 0.0 
     for j in range(D): 
      val += A[i,j] * x[j] 
     y[i] = val 
    return y 

Mam kompilacji tego pliku (seMatrixVectorExample.pyx) stosując następujący skrypt:

from distutils.core import setup 
from distutils.extension import Extension 
from Cython.Distutils import build_ext 
import numpy as np 

ext_modules=[ Extension("seMatrixVectorExample", 
         ["seMatrixVectorExample.pyx"], 
         libraries=["m"], 
         extra_compile_args = ["-ffast-math"])] 

setup(
    name = "seMatrixVectorExample", 
    cmdclass = {"build_ext": build_ext}, 
    include_dirs = [np.get_include()], 
    ext_modules = ext_modules 
) 

i przy użyciu następujących skrypt testowy ocenić wydajność:

import numpy as np 
from seMatrixVectorExample import matrix_vector_multiplication 
import time 

n_rows, n_cols = 1e6, 100 
np.random.seed(seed = 0) 

#initialize data matrix X and label vector Y 
A = np.random.random(size=(n_rows, n_cols)) 
np.require(A, requirements = ['C']) 

x = np.random.random(size=n_cols) 
x = np.require(x, requirements = ['C']) 

start_time = time.time() 
scores = matrix_vector_multiplication(A, x) 
print "cython runtime = %1.5f seconds" % (time.time() - start_time) 

start_time = time.time() 
py_scores = np.exp(A.dot(x)) 
print "numpy runtime = %1.5f seconds" % (time.time() - start_time) 

Na matrycy testowej z n_rows = 10e6 i n_cols = 100 uzyskać:

cython runtime = 0.08852 seconds 
numpy runtime = 0.04372 seconds 

---

Edit: Warto wspomnieć, że spowolnienie nie ustępuje nawet kiedy wdrożyć mnożenia macierzy w natywnym kodzie C, a tylko użycie Cython jako opakowanie.

void c_matrix_vector_multiplication(double* y, double* A, double* x, int N, int D) { 

    int i, j; 
    int index = 0; 
    double val; 

    for (i = 0; i < N; i++) { 
     val = 0.0; 
     for (j = 0; j < D; j++) { 
      val = val + A[index] * x[j]; 
      index++; 
      } 
     y[i] = val; 
     } 
    return; 
} 

i tu jest owinięcie Cython, które po prostu wysyła wskaźnik do pierwszego elementu y, A i x. :

import cython 
import numpy as np 
cimport numpy as np 

DTYPE = np.float64; 
ctypedef np.float64_t DTYPE_T 

# declare the interface to the C code 
cdef extern void c_multiply (double* y, double* A, double* x, int N, int D) 

@cython.boundscheck(False) 
@cython.wraparound(False) 
@cython.nonecheck(False) 
def multiply(np.ndarray[DTYPE_T, ndim=2, mode="c"] A, np.ndarray[DTYPE_T, ndim=1, mode="c"] x): 

    cdef int N = A.shape[0] 
    cdef int D = A.shape[1] 
    cdef np.ndarray[DTYPE_T, ndim=1, mode = "c"] y = np.empty(N, dtype = DTYPE) 

    c_multiply (&y[0], &A[0,0], &x[0], N, D) 

    return y 

Answer 1

OK udało się wreszcie uzyskać czasy działania lepsze niż NumPy!

Oto co (chyba) spowodował różnicę: NumPy dzwoni BLAS funkcje, które są zakodowane w Fortran zamiast C, w wyniku różnicy prędkości.

Uważam, że należy to odnotować, ponieważ wcześniej miałem wrażenie, że funkcje BLAS zostały zakodowane w C i nie widziałem powodu, dla którego działałyby zauważalnie szybciej niż druga natywna implementacja C, którą napisałem w pytaniu. .

W każdym przypadku, mogę teraz replikować wydajności za pomocą Cython + kursory funkcyjnych scipy Cython BLAS z scipy.linalg.cython_blas.

---

pod względem kompletności, tutaj jest nowy kod Cython blas_multiply.pyx:

import cython 
import numpy as np 
cimport numpy as np 
cimport scipy.linalg.cython_blas as blas 

DTYPE = np.float64 
ctypedef np.float64_t DTYPE_T 

@cython.boundscheck(False) 
@cython.wraparound(False) 
@cython.nonecheck(False) 

def blas_multiply(np.ndarray[DTYPE_T, ndim=2, mode="fortran"] A, np.ndarray[DTYPE_T, ndim=1, mode="fortran"] x): 
    #calls dgemv from BLAS which computes y = alpha * trans(A) + beta * y 
    #see: http://www.nag.com/numeric/fl/nagdoc_fl22/xhtml/F06/f06paf.xml 

    cdef int N = A.shape[0] 
    cdef int D = A.shape[1] 
    cdef int lda = N 
    cdef int incx = 1 #increments of x 
    cdef int incy = 1 #increments of y 
    cdef double alpha = 1.0 
    cdef double beta = 0.0 
    cdef np.ndarray[DTYPE_T, ndim=1, mode = "fortran"] y = np.empty(N, dtype = DTYPE) 

    blas.dgemv("N", &N, &D, &alpha, &A[0,0], &lda, &x[0], &incx, &beta, &y[0], &incy) 

    return y 

Oto kod, którego używam do zbudowania:

!/usr/bin/env python 

from distutils.core import setup 
from distutils.extension import Extension 
from Cython.Distutils import build_ext 

import numpy 
import scipy 

ext_modules=[ Extension("blas_multiply", 
         sources=["blas_multiply.pyx"], 
         include_dirs=[numpy.get_include(), scipy.get_include()], 
         libraries=["m"], 
         extra_compile_args = ["-ffast-math"])] 

setup(
    cmdclass = {'build_ext': build_ext}, 
    include_dirs = [numpy.get_include(), scipy.get_include()], 
    ext_modules = ext_modules, 
) 

A oto kod testowanie (zauważ, że tablice przekazywany do funkcji BLAS są F_CONTIGUOUS teraz)

import numpy as np 
from blas_multiply import blas_multiply 
import time 

#np.__config__.show() 
n_rows, n_cols = 1e6, 100 
np.random.seed(seed = 0) 

#initialize data matrix X and label vector Y 
X = np.random.random(size=(n_rows, n_cols)) 
Y = np.random.randint(low=0, high=2, size=(n_rows, 1)) 
Y[Y==0] = -1 
Z = X*Y 
Z.flags 
Z = np.require(Z, requirements = ['F']) 

rho_test = np.random.randint(low=-10, high=10, size= n_cols) 
set_to_zero = np.random.choice(range(0, n_cols), size =(np.floor(n_cols/2), 1), replace=False) 
rho_test[set_to_zero] = 0.0 
rho_test = np.require(rho_test, dtype=Z.dtype, requirements = ['F']) 

start_time = time.time() 
scores = blas_multiply(Z, rho_test) 
print "Cython runtime = %1.5f seconds" % (time.time() - start_time) 


Z = np.require(Z, requirements = ['C']) 
rho_test = np.require(rho_test, requirements = ['C']) 
start_time = time.time() 
py_scores = np.exp(Z.dot(rho_test)) 
print "Python runtime = %1.5f seconds" % (time.time() - start_time) 

Wynik z tego testu na moim komputerze jest:

Cython runtime = 0.04556 seconds 
Python runtime = 0.05110 seconds