python - una - ¿Qué está causando la ralentización de 2x en mi implementación Cython de la multiplicación vectorial de matriz?

que es un array en python

Question

Actualmente estoy tratando de implementar la multiplicación de vectores de matriz básica en Cython (como parte de un proyecto mucho más grande para reducir el cálculo ) y descubriendo que mi código es aproximadamente 2 Numpy.dot más lento que Numpy.dot .

Me pregunto si hay algo que me falta y que está provocando la desaceleración. Estoy escribiendo un código de Cython optimizado, declarando tipos de variables, requiriendo matrices contiguas y evitando errores de caché. Incluso traté de tener Cython como envoltorio y llamar al código C nativo (ver a continuación).

Me pregunto: ¿qué más podría hacer para acelerar mi implementación, así que corre tan rápido como NumPy para esta operación básica?

El código de Cython que estoy usando es 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

Estoy compilando este archivo ( seMatrixVectorExample.pyx ) con el siguiente 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
)

y utilizando el siguiente script de prueba para evaluar el rendimiento:

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)

Para una matriz de prueba con n_rows = 10e6 y n_cols = 100 obtengo:

cython runtime = 0.08852 seconds
numpy runtime = 0.04372 seconds

Editar: Vale la pena mencionar que la desaceleración persiste incluso cuando implemento la multiplicación de la matriz en el código C nativo, y solo uso Cython como envoltorio.

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

y aquí está el contenedor de Cython, que simplemente envía el puntero al primer elemento de y , A y 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 logró obtener tiempos de ejecución que son mejores que NumPy!

Esto es lo que (creo) causó la diferencia: NumPy está llamando funciones BLAS, que están codificadas en Fortran en lugar de C, lo que resulta en la diferencia de velocidad.

Creo que esto es importante tener en cuenta, ya que antes tenía la impresión de que las funciones BLAS estaban codificadas en C y no podía ver por qué funcionarían notablemente más rápido que la segunda implementación C nativa que publiqué en la pregunta.

En cualquier caso, ahora puedo replicar el rendimiento usando Cython + los punteros de función SciPy Cython BLAS de scipy.linalg.cython_blas.

Para completar, aquí está el nuevo código de 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

Aquí está el código que utilizo para compilar:

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

Y aquí está el código de prueba (tenga en cuenta que las matrices pasadas a la función BLAS son F_CONTIGUOUS ahora)

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)

El resultado de esta prueba en mi máquina es:

Cython runtime = 0.04556 seconds
Python runtime = 0.05110 seconds