Was verursacht die 2x Verlangsamung in meiner Cython-Implementierung der Matrix-Vektor-Multiplikation?

Question

Ich versuche derzeit grundlegende Matrix Vektor Multiplikation in Cython zu implementieren (als Teil einer viel larger project to reduce computation) und feststellen, dass mein Code etwa 2x langsamer als Numpy.dot ist.

Ich frage mich, ob es etwas gibt, das ich vermisse, das in der Verlangsamung resultiert. Ich schreibe optimierten Cython-Code, deklariere Variablentypen, benötige zusammenhängende Arrays und vermeide Cache-Misses. Ich habe sogar versucht, Cython als Wrapper zu verwenden und nativen C-Code aufzurufen (siehe unten).

Ich frage mich: was sonst könnte ich tun, um meine Implementierung zu beschleunigen, so läuft so schnell wie NumPy für diese grundlegende Operation?

---

Der Cython Code, ich verwende ist 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 

ich diese Datei bin Kompilieren (seMatrixVectorExample.pyx) das folgende Skript:

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 
) 

und mit dem folgenden Testskript zur Beurteilung der Leistung:

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) 

Für eine Testmatrix mit n_rows = 10e6 und n_cols = 100 ich:

cython runtime = 0.08852 seconds 
numpy runtime = 0.04372 seconds 

---

Edit: Es ist erwähnenswert, dass die Verlangsamung bleibt auch, wenn ich die Matrixmultiplikation in nativen C-Code implementieren, und nur Gebrauch Cython als 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; 
} 

und hier ist die Cython Wrapper, der nur den Zeiger auf das erste Element sendet von y, A und 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 endlich geschafft, Laufzeiten zu bekommen, die besser als NumPy sind!

Hier ist, was (ich denke) den Unterschied verursacht: NumPy ruft BLAS-Funktionen, die in Fortran statt C codiert sind, in der Geschwindigkeitsdifferenz.

Ich denke, dass dies wichtig zu beachten ist, da ich zuvor den Eindruck hatte, dass die BLAS-Funktionen in C codiert waren und nicht sehen konnten, warum sie merklich schneller als die zweite native C-Implementierung, die ich in der Frage veröffentlicht, laufen würden .

In jedem Fall habe ich jetzt die Leistung durch die Verwendung Cython + die SciPy Cython BLAS Funktionszeiger von scipy.linalg.cython_blas.

---

Der Vollständigkeit halber hier replizieren kann, ist der neue Cython Code 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 

Hier ist der Code, den ich verwende:

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

Und hier ist der Testcode (beachten Sie, dass Arrays die BLAS-Funktion übergeben sind F_CONTIGUOUS jetzt) ​​

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) 

Das Ergebnis aus diesem Test auf meinem Rechner ist:

Cython runtime = 0.04556 seconds 
Python runtime = 0.05110 seconds