Che cosa sta causando il rallentamento 2x nella mia implementazione Cython della moltiplicazione dei vettori di matrici?

Question

Attualmente sto cercando di implementare la moltiplicazione di vettore matrice di base in Cython (come parte di un gran numero di larger project to reduce computation) e di scoprire che il mio codice è di circa 2 volte più lento di Numpy.dot.

Mi chiedo se c'è qualcosa che mi manca che sta causando il rallentamento. Sto scrivendo codice Cython ottimizzato, dichiarando tipi variabili, che richiedono array contigui ed evitando errori di cache. Ho anche provato ad avere Cython come wrapper e chiamare il codice C nativo (vedi sotto).

Mi chiedo: cos'altro posso fare per accelerare la mia implementazione, quindi funziona velocemente come NumPy per questa operazione di base?

---

Il codice Cython che sto usando è 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 

sto compilando questo file di (seMatrixVectorExample.pyx) utilizzando il seguente script:

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 
) 

e utilizzando il seguente script di test per valutare le prestazioni:

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) 

Per una matrice prova con n_rows = 10e6 e n_cols = 100 ottengo:

cython runtime = 0.08852 seconds 
numpy runtime = 0.04372 seconds 

---

Edit: Vale la pena ricordare che il rallentamento persiste anche quando implemento la moltiplicazione matrice codice nativo C, e utilizzare solo Cython come wrapper.

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; 
} 

e qui è l'involucro Cython, che invia solo il puntatore al primo elemento di y, A e 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 finalmente riuscito a ottenere runtime che sono meglio di NumPy!

Ecco cosa (penso) ha causato la differenza: NumPy chiama le funzioni BLAS, che sono codificate in Fortran anziché C, risultando nella differenza di velocità.

Penso che sia importante notare, poiché in precedenza avevo l'impressione che le funzioni BLAS fossero codificate in C e non potevo capire perché sarebbero state eseguite in modo notevolmente più veloce della seconda implementazione C nativa che ho postato nella domanda .

In entrambi i casi, ora posso replicare le prestazioni utilizzando Cython + i puntatori a funzione SciPy Cython BLAS da scipy.linalg.cython_blas.

---

Per completezza, ecco il nuovo codice 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 

Ecco il codice che utilizzo per creare:

!/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, 
) 

E qui è il codice di test (si noti che gli array passati alla funzione BLAS sono F_CONTIGUOUS ora)

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) 

Il risultato di questo test sulla mia macchina è:

Cython runtime = 0.04556 seconds 
Python runtime = 0.05110 seconds