这里有一篇帖子:https://gist.github.com/JonathanRaiman/f2ce5331750da7b2d4e9,展示了通过调用Fortran库(BLAS / LAPACK / Intel MKL / OpenBLAS / 与NumPy一起安装的任何内容)可以大大提高速度。由于SciPy库已被弃用,我花费了许多时间进行工作,最终成功编译,但没有结果。它比NumPy快2倍。不幸的是,正如另一个用户指出的那样,Fortran例程总是将输出矩阵添加到新计算的结果中,因此它只能在第一次运行时与NumPy匹配。即A:= alpha * x * y.T + A
。因此需要找到一种快速解决方案来解决这个问题。
[更新:对于那些想要使用SciPy接口的人,请访问此处https://github.com/scipy/scipy/blob/master/scipy/linalg/cython_blas.pyx,因为他们已经优化了CPDEF语句中对BLAS / LAPACK的调用,只需将其复制/粘贴到您的Cython脚本中# Python-accessible wrappers for testing:
。此外,在上面链接的cython_lapack.pyx可用,但没有Cython测试脚本]
测试脚本
import numpy as np;
from cyblas import outer_prod;
a=np.random.randint(0,100, 1000);
b=np.random.randint(0,100, 1000);
a=a.astype(np.float64)
b=b.astype(np.float64)
cy_outer=np.zeros((a.shape[0],b.shape[0]));
np_outer=np.zeros((a.shape[0],b.shape[0]));
%timeit outer_prod(a,b,cy_outer)
#%timeit outer_prod(a,b) #use with fixed version instead of above line, results will automatically update cy_outer
%timeit np.outer(a,b, np_outer)
100 loops, best of 3: 2.83 ms per loop
100 loops, best of 3: 6.58 ms per loop
将PYX文件编译成cyblas.pyx(基本上是np.ndarray版本)
# END TEST SCRIPT
import cython
import numpy as np
cimport numpy as np
from cpython cimport PyCapsule_GetPointer
cimport scipy.linalg.cython_blas
cimport scipy.linalg.cython_lapack
import scipy.linalg as LA
REAL = np.float64
ctypedef np.float64_t REAL_t
ctypedef np.uint64_t INT_t
cdef int ONE = 1
cdef REAL_t ONEF = <REAL_t>1.0
ctypedef void (*dger_ptr) (const int *M, const int *N, const double *alpha, const double *X, const int *incX, double *Y, const int *incY, double *A, const int * LDA) nogil
cdef dger_ptr dger=<dger_ptr>PyCapsule_GetPointer(LA.blas.dger._cpointer, NULL) # A := alpha*x*y.T + A
cpdef outer_prod(_x, _y, _output):
#cpdef outer_prod(_x, _y): #comment above line & use this to use the reset output matrix to zeros
cdef REAL_t *x = <REAL_t *>(np.PyArray_DATA(_x))
cdef int M = _y.shape[0]
cdef int N = _x.shape[0]
#cdef np.ndarray[np.float64_t, ndim=2, order='c'] _output = np.zeros((M,N)) #slow fix to uncomment to reset output matrix to zeros
cdef REAL_t *y = <REAL_t *>(np.PyArray_DATA(_y))
cdef REAL_t *output = <REAL_t *>(np.PyArray_DATA(_output))
with nogil:
dger(&M, &N, &ONEF, y, &ONE, x, &ONE, output, &M)
非常感谢。希望这篇文章能够节省其他人的时间(它几乎可以起作用)- 实际上,正如我评论的那样,它在匹配NumPy时仅起作用1次,然后每个后续调用都会再次添加到结果矩阵中。如果我将输出矩阵重置为0并重新运行结果,则匹配NumPy。奇怪…虽然如果取消注释上面的几行,它会起作用,但仅以NumPy的速度。另一种选择是使用
memset
,并将在另一篇文章中介绍…我只是还没有弄清楚如何精确调用它。
dger(M, N, ALPHA, X INCX, Y, INCY, A, LDA)
执行A := alpha*x*y**T + A
。因此,如果要得到X
和Y
的外积,则A
应该全部为零。 - user7138814