Производительность ядра Гаусса

Question

Следующий метод вычисляет ядро ​​Гаусса:

import numpy as np
def gaussian_kernel(X, X2, sigma):
    """
    Calculate the Gaussian kernel matrix

        k_ij = exp(-||x_i - x_j||^2 / (2 * sigma^2))

    :param X: array-like, shape=(n_samples_1, n_features), feature-matrix
    :param X2: array-like, shape=(n_samples_2, n_features), feature-matrix
    :param sigma: scalar, bandwidth parameter

    :return: array-like, shape=(n_samples_1, n_samples_2), kernel matrix
    """

    norm = np.square(np.linalg.norm(X[None,:,:] - X2[:,None,:], axis=2).T)    
    return np.exp(-norm/(2*np.square(sigma)))

# Usage example
%timeit gaussian_kernel(np.random.rand(5000, 10), np.random.rand(5000, 10), 1)

1,43 с ± 39,3 мс на л oop (среднее ± стандартное отклонение из 7 прогонов, 1 л oop каждое)

Мой вопрос: есть ли способы повысить производительность, используя numpy?

Answer 1

Для довольно маленьких массивов вы можете написать простую реализацию l oop и скомпилировать ее с помощью Numba. Для больших массивов переформулировка алгебры c с использованием np.dot () будет быстрее.

Пример

#from version 0.43 until 0.47 this has to be set before importing numba
#Bug: https://github.com/numba/numba/issues/4689
from llvmlite import binding
binding.set_option('SVML', '-vector-library=SVML')
import numba as nb
import numpy as np

@nb.njit(fastmath=True,error_model="numpy",parallel=True)
def gaussian_kernel_2(X, X2, sigma):
    res=np.empty((X.shape[0],X2.shape[0]),dtype=X.dtype)
    for i in nb.prange(X.shape[0]):
        for j in range(X2.shape[0]):
            acc=0.
            for k in range(X.shape[1]):
                acc+=(X[i,k]-X2[j,k])**2/(2*sigma**2)
            res[i,j]=np.exp(-1*acc)
    return res

Время

X1=np.random.rand(5000, 10)
X2=np.random.rand(5000, 10)

#Your solution
%timeit gaussian_kernel(X1,X2, 1)
#511 ms ± 10.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit gaussian_kernel_2(X1,X2, 1)
#90.1 ms ± 9.14 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Answer 2

Этот пост: { ссылка } дал ответ.

Вкратце, для копирования части numpy:

import numpy as np
def gaussian_kernel(X, X2, sigma):
    """
    Calculate the Gaussian kernel matrix

        k_ij = exp(-||x_i - x_j||^2 / (2 * sigma^2))

    :param X: array-like, shape=(n_samples_1, n_features), feature-matrix
    :param X2: array-like, shape=(n_samples_2, n_features), feature-matrix
    :param sigma: scalar, bandwidth parameter

    :return: array-like, shape=(n_samples_1, n_samples_2), kernel matrix
    """
    X_norm = np.sum(X ** 2, axis = -1)
    X2_norm = np.sum(X2 ** 2, axis = -1)
    norm = X_norm[:,None] + X2_norm[None,:] - 2 * np.dot(X, X2.T)
    return np.exp(-norm/(2*np.square(sigma)))

# Timing
%timeit gaussian_kernel(np.random.rand(5000, 10), np.random.rand(5000, 10), 1)

976 мс ± 73,5 мс на л oop (среднее ± стандартное отклонение из 7 прогонов, 1 л oop каждый)