如何通过使用 numpy einsum 和 numexpr 来计算内核函数来加快性能?

Question

我正在尝试为sklearn库中的方法定义一些著名的内核，如RBF，双曲正切，傅里叶等svm.SVR。我开始研究rbf(我知道 svm 中有一个用于 rbf 的默认内核，但我需要定义一个以便以后能够自定义它)，并在这里找到了一些有用的链接并选择了这个：

def my_kernel(X,Y):
K = np.zeros((X.shape[0],Y.shape[0]))
for i,x in enumerate(X):
for j,y in enumerate(Y):
K[i,j] = np.exp(-1*np.linalg.norm(x-y)**2)
return K
clf=SVR(kernel=my_kernel)

我使用这个是因为我可以将它用于具有不同形状的火车(形状为 [3850,4])和测试数据(形状为 [1200,4])。但问题是它太慢了，我必须等待很长时间才能得到结果。我什至在cython中使用了静态类型和内存视图，但它的性能不如默认svmrbf内核。我还发现了这个链接，它大约是同样的问题，但与numpy.einsum和numexpr.evaluate一起工作对我来说有点混乱。事实证明，就速度性能而言，这是最好的代码：

从 scipy.linalg.blas import sgemm

def app2(X, gamma, var):
X_norm = -gamma*np.einsum('ij,ij->i',X,X)
return ne.evaluate('v * exp(A + B + C)', {
'A' : X_norm[:,None],
'B' : X_norm[None,:],
'C' : sgemm(alpha=2.0*gamma, a=X, b=X, trans_b=True),
'g' : gamma,
'v' : var
})

这段代码仅适用于一个输入 (X)，我找不到一种方法来针对我的情况修改它(两个具有两种不同大小的输入 - 内核函数根据 svm 文档获取形状为 (m，n) 和 (l，n) 的矩阵和输出 (m，l) )。我想我只需要在第二个代码中替换第一个代码中的K[i,j] = np.exp(-1*np.linalg.norm(x-y)**2)即可加快速度。任何帮助将不胜感激。

Answer 1

你可以通过 pythran 管道传输你的原始内核

kernel.py：

import numpy as np
#pythran export my_kernel(float64[:,:], float64[:,:])
def my_kernel(X,Y):
K = np.zeros((X.shape[0],Y.shape[0]))
for i,x in enumerate(X):
for j,y in enumerate(Y):
K[i,j] = np.exp(-1*np.linalg.norm(x-y)**2)
return K

编译步骤：

> pythran kernel.py

没有重写步骤(不过你需要将内核放在一个单独的文件中)，而且加速是显着的：在我的笔记本电脑上快 19 倍，使用

> python -m timeit -s 'from numpy.random import random; x = random((100,100)); y = random((100,100)); from svr_kernel import my_kernel as k' 'k(x,y)'

以收集时间。

Answer 2

三种可能的变体

变体 1 和 3 利用

(a-b)^2 = a^2 + b^2 - 2ab

如此处或此处所述。但对于像小的二维这样的特殊情况，变体 2 也是可以的。

import numpy as np
import numba as nb
import numexpr as ne
from scipy.linalg.blas import sgemm
def vers_1(X,Y, gamma, var):
X_norm = -gamma*np.einsum('ij,ij->i',X,X)
Y_norm = -gamma*np.einsum('ij,ij->i',Y,Y)
return ne.evaluate('v * exp(A + B + C)', {
'A' : X_norm[:,None],
'B' : Y_norm[None,:],
'C' : sgemm(alpha=2.0*gamma, a=X, b=Y, trans_b=True),
'g' : gamma,
'v' : var
})
#Maybe easier to read but slow, if X.shape[1] gets bigger
@nb.njit(fastmath=True,parallel=True)
def vers_2(X,Y):
K = np.empty((X.shape[0],Y.shape[0]),dtype=X.dtype)
for i in nb.prange(X.shape[0]):
for j in range(Y.shape[0]):
sum=0.
for k in range(X.shape[1]):
sum+=(X[i,k]-Y[j,k])**2
K[i,j] = np.exp(-1*sum)
return K
@nb.njit(fastmath=True,parallel=True)
def vers_3(A,B):
dist=np.dot(A,B.T)
TMP_A=np.empty(A.shape[0],dtype=A.dtype)
for i in nb.prange(A.shape[0]):
sum=0.
for j in range(A.shape[1]):
sum+=A[i,j]**2
TMP_A[i]=sum
TMP_B=np.empty(B.shape[0],dtype=A.dtype)
for i in nb.prange(B.shape[0]):
sum=0.
for j in range(B.shape[1]):
sum+=B[i,j]**2
TMP_B[i]=sum
for i in nb.prange(A.shape[0]):
for j in range(B.shape[0]):
dist[i,j]=np.exp((-2.*dist[i,j]+TMP_A[i]+TMP_B[j])*-1)
return dist

计时

gamma = 1.
var = 1.
X=np.random.rand(3850,4)
Y=np.random.rand(1200,4)
res_1=vers_1(X,Y, gamma, var)
res_2=vers_2(X,Y)
res_3=vers_3(X,Y)
np.allclose(res_1,res_2)
np.allclose(res_1,res_3)

%timeit res_1=vers_1(X,Y, gamma, var)
19.8 ms ± 615 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit res_2=vers_2(X,Y)
16.1 ms ± 938 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit res_3=vers_3(X,Y)
13.5 ms ± 162 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)