Implementando un producto interno de python más rápido con BLAS

numpy cython

Question

Encontré este tutorial útil sobre el uso de funciones BLAS de bajo nivel (implementadas en Cython) para obtener mejoras a gran velocidad sobre las rutinas estándar de álgebra lineal numpy en python. Ahora, he conseguido que el producto vectorial funcione bien. Primero linalg.pyx siguiente como linalg.pyx :

import cython
import numpy as np
cimport numpy as np

from libc.math cimport exp
from libc.string cimport memset

from scipy.linalg.blas import fblas

REAL = np.float64
ctypedef np.float64_t REAL_t

cdef extern from "/home/jlorince/flda/voidptr.h":
    void* PyCObject_AsVoidPtr(object obj)

ctypedef double (*ddot_ptr) (const int *N, const double *X, const int *incX, const double *Y, const int *incY) nogil
cdef ddot_ptr ddot=<ddot_ptr>PyCObject_AsVoidPtr(fblas.ddot._cpointer)  # vector-vector multiplication 

cdef int ONE = 1
def vec_vec(syn0, syn1, size):
    cdef int lSize = size
    f = <REAL_t>ddot(&lSize, <REAL_t *>(np.PyArray_DATA(syn0)), &ONE, <REAL_t *>(np.PyArray_DATA(syn1)), &ONE)
    return f

(código fuente para voidptr.h disponible here )

Una vez que lo compilo, funciona bien, y es definitivamente más rápido que np.inner :

In [1]: import linalg
In [2]: import numpy as np
In [3]: x = np.random.random(100)
In [4]: %timeit np.inner(x,x)
1000000 loops, best of 3: 1.61 µs per loop
In [5]: %timeit linalg.vec_vec(x,x,100)
1000000 loops, best of 3: 483 ns per loop
In [8]: np.all(np.inner(x,x)==linalg.vec_vec(x,x,100))
Out[8]: True

Ahora, esto es genial, pero solo funciona para el caso de calcular el punto / producto interno (equivalente en este caso) de dos vectores . Lo que necesito hacer ahora, implementar funciones similares (que espero ofrezcan aceleraciones similares) para hacer productos internos de matriz vectorial. Es decir, quiero replicar la funcionalidad de np.inner cuando se pasa una matriz 1D y una matriz 2D:

In [4]: x = np.random.random(5)
In [5]: y = np.random.random((5,5))
In [6]: np.inner(x,y)
Out[6]: array([ 1.42116225,  1.13242989,  1.95690196,  1.87691992,  0.93967486])

Esto es equivalente a calcular los productos internos / de puntos (de nuevo, equivalentes para matrices 1D) de la matriz 1D y cada fila de la matriz:

In [32]: [np.inner(x,row) for row in y]
Out[32]:
[1.4211622497461549, 1.1324298918119025, 1.9569019618096966,1.8769199192990056, 0.93967485730285505]

Por lo que he visto de la documentación de BLAS, creo que debo comenzar con algo como esto (usando dgemv):

ctypedef double (*dgemv_ptr) (const str *TRANS, const int *M, const int *N, const double *ALPHA, const double *A, const int *LDA, const double *X, const int *incX, const double *BETA, const double *Y, const int *incY)
cdef dgemv_ptr dgemv=<dgemv>PyCObject_AsVoidPtr(fblas.dgemv._cpointer)  # matrix vector multiplication

Pero necesito ayuda (a) que define la función real que puedo usar en Python (es decir, una función de vec-matrix análoga a vec_vec anterior), y (b) saber cómo usarla para replicar correctamente el comportamiento de np.inner , que Es lo que necesito para el modelo que estoy implementando.

Edición: Here está el enlace a la documentación relevante de BLAS para dgemv, que necesito usar, que se confirma aquí:

In [13]: np.allclose(scipy.linalg.blas.fblas.dgemv(1.0,y,x), np.inner(x,y))
Out[13]: True

Pero usarlo fuera de la caja como esto es en realidad más lento que puro np.inner, así que todavía necesito ayuda con la implementación de Cython.

Edit2 Aquí está mi último intento de esto, que compila bien pero bloquea Python con una falla de segmentación cada vez que intento ejecutarlo:

cdef int ONE = 1
cdef char tr = ''n''
cdef REAL_t ZEROF = <REAL_t>0.0
cdef REAL_t ONEF = <REAL_t>1.0
def mat_vec(mat,vec,mat_rows,mat_cols):
    cdef int m = mat_rows
    cdef int n = mat_cols
    out = <REAL_t>dgemv(&tr, &m, &n, &ONEF, <REAL_t *>(np.PyArray_DATA(mat)), &m, <REAL_t *>(np.PyArray_DATA(vec)), &ONE, &ZEROF, NULL, &ONE)
    return out

Después de compilar, trato de ejecutar linalg.mat_vec(y,x,5,5) , (usando la misma xey que las anteriores) pero esto simplemente falla. Creo que estoy cerca, pero no sé qué más cambiar ...

Answer 1

Per @Pietro Saccardi:

int dgemv_(char *trans, integer *m, integer *n, doublereal *
           alpha, doublereal *a, integer *lda, doublereal *x, integer *incx, 
           doublereal *beta, doublereal *y, integer *incy)

...

Y      - DOUBLE PRECISION array of DIMENSION at least   
         ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = ''N'' or ''n''   
         and at least   
         ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.   
         Before entry with BETA non-zero, the incremented array Y   
         must contain the vector y. On exit, Y is overwritten by the
         updated vector y.

Dudo que puedas usar NULL para Y en la llamada.