我用MKL编译了numpy 1.6.2和scipy,希望有更好的性能。
目前我有一个严重依赖 np.einsum() 的代码,并且我被告知 einsum 不适用于 MKL,因为几乎没有矢量化。 =(
所以我想用 np.dot() 和切片重新编写一些代码,只是为了能够获得一些多核加速。
我真的很喜欢 np.einsum() 的简单性并且可读性很好。
无论如何,例如,我有一个以下形式的多维矩阵乘法:
np.einsum('mi,mnijqk->njqk',A,B)
那么我如何在 np.dot() 高效 MKL 操作中转换类似这样的东西,或其他 3,4 和 5 维数组乘法?
我会广告更多信息:
我正在计算这个方程:
为此,我使用以下代码:
np.einsum('mn,mni,nij,nik,mi->njk',a,np.exp(b[:,:,np.newaxis]*U[np.newaxis,:,:]),P,P,X)
这并没有那么快,用 cython 编码的相同内容快了 5 倍:
#STACKOVERFLOW QUESTION:
from __future__ import division
import numpy as np
cimport numpy as np
cimport cython
cdef extern from "math.h":
double exp(double x)
DTYPE = np.float
ctypedef np.float_t DTYPE_t
@cython.boundscheck(False) # turn of bounds-checking for entire function
def cython_DX_h(np.ndarray[DTYPE_t, ndim=3] P, np.ndarray[DTYPE_t, ndim=1] a, np.ndarray[DTYPE_t, ndim=1] b, np.ndarray[DTYPE_t, ndim=2] U, np.ndarray[DTYPE_t, ndim=2] X, int I, int M):
assert P.dtype == DTYPE and a.dtype == DTYPE and b.dtype == DTYPE and U.dtype == DTYPE and X.dtype == DTYPE
cdef np.ndarray[DTYPE_t,ndim=3] DX_h=np.zeros((N,I,I),dtype=DTYPE)
cdef unsigned int j,n,k,m,i
for n in range(N):
for j in range(I):
for k in range(I):
aux=0
for m in range(N):
for i in range(I):
aux+=a[m,n]*exp(b[m,n]*U[n,i])*P[n,i,j]*P[n,i,k]*X[m,i]
DX_h[n,j,k]=aux
return DX_h
有没有办法在纯Python中以cython的性能来做到这一点? (我还没能弄清楚如何张量这个方程)
无法在这个 cython 代码中进行 prange,有很多 gil 和 nogil 错误。