Thực hiện sản phẩm bên trong python nhanh hơn với BLAS

Question

Tôi đã tìm thấy this useful tutorial về sử dụng chức năng BLAS mức thấp (được triển khai trong Cython) để có được cải thiện tốc độ lớn so với quy trình đại số tuyến tính tiêu chuẩn trong python. Bây giờ, tôi đã mua thành công sản phẩm vector hoạt động tốt. Trước tiên tôi lưu sau khi 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 

(mã nguồn cho voidptr.h sẵn here)

Khi tôi biên dịch nó, nó hoạt động tốt, và chắc chắn là nhanh hơn so với 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 

Bây giờ, điều này thật tuyệt vời, nhưng chỉ hoạt động đối với trường hợp tính toán dấu chấm/sản phẩm bên trong (tương đương trong trường hợp này) của hai vector . Những gì tôi cần làm bây giờ, thực hiện các chức năng tương tự (mà tôi hy vọng sẽ cung cấp tăng tốc tương tự) để làm vector-matrix bên trong sản phẩm. Đó là, tôi muốn tái tạo các chức năng của np.inner khi đã thông qua một mảng 1D và một ma trận 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]) 

này tương đương với việc tính toán các sản phẩm bên trong/chấm (một lần nữa, tương đương cho mảng 1D) của mảng 1D và mỗi hàng của ma trận:

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

Từ những gì tôi đã nhìn thấy các tài liệu BLAS, tôi nghĩ tôi cần phải bắt đầu với một cái gì đó như thế này (sử dụng 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 

Nhưng tôi cần trợ giúp (a) xác định hàm thực tế mà tôi có thể sử dụng trong Python (tức là một hàm vec-matrix tương tự như vec_vec ở trên), và (b) biết cách sử dụng nó để sao chép đúng hành vi của np.inner, đó là những gì tôi cần cho mô hình mà tôi đang triển khai.

Edit:Here là liên kết đến tài liệu BLAS liên quan cho dgemv, mà tôi cần phải được sử dụng, được khẳng định ở đây:

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

Nhưng sử dụng nó ra khỏi hộp như thế này thực sự là chậm so với np.inner thuần túy, vì vậy tôi vẫn cần trợ giúp với việc triển khai Cython.

Edit2 Đây là nỗ lực mới nhất của tôi lúc này, trong đó biên dịch tốt nhưng bị treo python với một lỗi segmentation bất cứ khi nào tôi cố gắng chạy nó:

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 

Sau khi biên soạn, tôi cố gắng chạy , (sử dụng cùng x và y như trên) nhưng điều này chỉ bị treo. Tôi nghĩ mình là gần, nhưng không biết làm gì khác để thay đổi ...

Answer 1

mỗi @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. 

tôi nghi ngờ bạn có thể sử dụng NULL cho Y trong cuộc gọi.