Jak skutecznie obliczyć jądro Gaussa w numpy [zamknięte]

Mam tablicy numpy z m kolumn i n wierszy, kolumn i wymiarach będących rzędów punktów danych.

Teraz muszę obliczyć wartości jądra dla każdej kombinacji punktów danych.

Dla liniowego jądra mogę po prostu zrobić$K(\mathbf{x}_i,\mathbf{x}_j) = \langle \mathbf{x}_i,\mathbf{x}_j \rangle$dot(X,X.T)

Jak mogę skutecznie obliczyć wszystkie wartości dla jądra Gaussa z podanym s ?$K(\mathbf{x}_i,\mathbf{x}_j) = \exp{-\frac{\|\mathbf{x}_i-\mathbf{x}_j\|_2^2}{s^2}}$

Myślę, że głównym problemem jest efektywne uzyskanie odległości parami. Gdy już to zrobisz, reszta jest mądra.

Aby to zrobić, prawdopodobnie chcesz użyć scipy. Ta funkcja scipy.spatial.distance.pdistspełnia scipy.spatial.distance.squareformTwoje oczekiwania i prawdopodobnie ułatwi Ci życie.

Więc jeśli chcesz macierz jądra, to robisz

from scipy.spatial.distance import pdist, squareform
  # this is an NxD matrix, where N is number of items and D its dimensionalites
X = loaddata() 
pairwise_dists = squareform(pdist(X, 'euclidean'))
K = scip.exp(-pairwise_dists ** 2 / s ** 2)

Dokumentacja znajduje się tutaj . 

Jako niewielki dodatek do odpowiedzi Bayerja, pdistfunkcja Scipy'ego może bezpośrednio obliczyć kwadratowe normy euklidesowe, nazywając to pdist(X, 'sqeuclidean'). Pełny kod można następnie zapisać bardziej efektywnie jako

from scipy.spatial.distance import pdist, squareform
  # this is an NxD matrix, where N is number of items and D its dimensionalites
X = loaddata() 
pairwise_sq_dists = squareform(pdist(X, 'sqeuclidean'))
K = scip.exp(-pairwise_sq_dists / s**2)

Możesz również ręcznie napisać kwadratowy formularz:

import numpy as np
def vectorized_RBF_kernel(X, sigma):
    # % This is equivalent to computing the kernel on every pair of examples
    X2 = np.sum(np.multiply(X, X), 1) # sum colums of the matrix
    K0 = X2 + X2.T - 2 * X * X.T
    K = np.power(np.exp(-1.0 / sigma**2), K0)
    return K

PS, ale działa to o 30% wolniej

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)

co jest równe

clf=SVR(kernel="rbf",gamma=1)

Możesz skutecznie obliczyć RBF na podstawie powyższego kodu, że wartość gamma wynosi 1, ponieważ jest to stała, o którą prosiłeś, jest również taka sama.

Myślę, że to pomoże:

def GaussianKernel(v1, v2, sigma):
    return exp(-norm(v1-v2, 2)**2/(2.*sigma**2))