python - 使用 BLAS 实现更快的 python 内积

Question

我找到了 this useful tutorial关于使用低级 BLAS 函数(在 Cython 中实现)在 python 中获得比标准 numpy 线性代数例程更大的速度改进。现在，我已经成功地使 vector 产品工作正常。首先我将以下内容保存为 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

(voidptr.h的源代码可用here)

一旦我编译它，它就可以正常工作，而且肯定比 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

现在，这很好，但仅适用于计算两个 向量 的点/内积(在这种情况下等效)的情况。我现在需要做的是，实现类似的功能(我希望能提供类似的加速)来做向量矩阵 inner 产品。也就是说，我想在传递一维数组和二维矩阵时复制 np.inner 的功能:

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])

这相当于计算一维数组和矩阵每一行的内积/点积(同样，相当于一维数组):

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

根据我对 BLAS 文档的了解，我认为我需要从这样的内容开始(使用 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

但我需要帮助 (a) 定义我可以在 Python 中使用的实际函数(即类似于上述 vec_vec 的 vec-matrix 函数)，以及 (b ) 知道如何使用它来正确复制 np.inner 的行为，这是我正在实现的模型所需要的。

编辑: Here是我需要使用的 dgemv 相关 BLAS 文档的链接，在此确认:

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

但是像这样直接使用它实际上比纯 np.inner 慢，所以我仍然需要 Cython 实现方面的帮助。

Edit2这是我最近的尝试，它编译得很好，但每当我尝试运行它时，它都会因段错误而崩溃:

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

编译后，我尝试运行 linalg.mat_vec(y,x,5,5)，(使用与上面相同的 x 和 y)但这只是崩溃。我想我已经接近了，但不知道还有什么要改变的......

Answer 1

根据@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.

我怀疑你可以在调用中将 NULL 用于 Y。