我找到了这个有用的教程,介绍如何在Python中使用基于Cython实现的低级BLAS函数来大幅提高标准numpy线性代数库的运算速度。现在,我已经成功地让向量乘积正常工作了。首先,我将以下代码保存为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 的源代码可在 这里 获取)
一旦我编译它,它就可以正常工作,而且肯定比 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
现在,这很好,但仅适用于计算两个向量的点积/内积(在这种情况下等效)。现在我需要做的是实现类似功能的函数(希望能提供类似的加速),用于执行向量-矩阵内积。也就是说,当传递一个1D数组和一个2D矩阵给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])
这相当于计算 1 维数组与矩阵的每一行的内积/点积(对于 1 维数组来说也是等效的):
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
的行为,这是我正在实现的模型所需的。
编辑: 这里是我需要使用的相关BLAS文档的链接,确切地说是dgemv。
In [13]: np.allclose(scipy.linalg.blas.fblas.dgemv(1.0,y,x), np.inner(x,y))
Out[13]: True
但是,像这样直接使用它比纯np.inner要慢,所以我仍然需要Cython实现的帮助。
Edit2 这是我最新的尝试,编译没有问题,但每当我尝试运行它时,Python会崩溃并显示分段错误:
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),但是它会崩溃。我觉得我已经很接近了,但不知道还需要做什么改变...
np.dot
呢? - BeRecursivenp.inner
一样对一个1D数组和一个2D矩阵进行精确计算(即数组和矩阵的每一行的点积/内积),这比迭代矩阵并单独计算每个点积/内积要快得多。 - moustachiodgemv
要求一个非空的y
矩阵来存储结果,而你传递了 NULL? - Pietro Saccardi