Как я могу повысить производительность, используя numpy einsum и numexpr при вычислении функций ядра?

Я пытаюсь определить несколько известных ядер, таких как RBF, гиперболический тангенс, Фурье и т. Д., Для метода svm.SVR в библиотеке sklearn. Я начал работать над 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, но его производительность не так хороша, как ядро ​​svm rbf по умолчанию. Я также нашел это ссылка, которая посвящена той же проблеме, но работа с numpy.einsum и numexpr.evaluate меня немного сбивает с толку. Получается, что это был лучший код с точки зрения быстродействия:

из scipy.linalg.blas импортировать 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), и я не смог найти способ изменить его для моего случая (два входа с двумя разными размерами. Функция ядра получает матрицы с формой (m, n) и (l, n) и выводит (m,l) в соответствии с svm документы ). Думаю, мне нужно только заменить K[i,j] = np.exp(-1*np.linalg.norm(x-y)**2) из первого кода во втором, чтобы ускорить его. Любая помощь будет оценена по достоинству.


python performance numpy-einsum numexpr
person Masoud Masoumi Moghadam    schedule 23.04.2019    source источник
comment
Эта проблема связана с евклидовыми расстояниями. см.: stackoverflow.com/a/42994680/4045774 или stackoverflow.com/a/53380192/4045774   -  person max9111    schedule 24.04.2019
comment
Вы правы. Но я не знаю, как использовать одно и то же решение для двух матриц разной формы!   -  person Masoud Masoumi Moghadam    schedule 24.04.2019

Ответы (2)


arrow_upward
1
arrow_downward

Три возможных варианта

Варианты 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)
person max9111    schedule 24.04.2019

arrow_upward
1
arrow_downward

Вы можете просто передать исходное ядро ​​через pythran

ядро.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)'

собирать тайминги.

person serge-sans-paille    schedule 27.04.2019