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

12

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.(xja,xjot)=xja,xjotdot(X,X.T)

Jak mogę skutecznie obliczyć wszystkie wartości dla jądra Gaussa z podanym s ?K.(xja,xjot)=exp-xja-xjot2)2)s2)

Peter Smit
źródło
1
Cóż, jeśli nie przejmujesz się zbytnio współczynnikiem dwóch obliczeń, zawsze możesz po prostu zrobić a następnie gdzie oczywiście to ty element . Prawdopodobnie nie jest to jednak najbardziej stabilna liczbowo. K ( x i , x j ) = exp ( - ( S i i + S j j - 2 S i j ) / s 2 ) S i j ( i , j ) SS.=XXT.K.(xja,xjot)=exp(-(S.jaja+S.jotjot-2)S.jajot)/s2))S.jajot(ja,jot)S.
kardynał
2
(Lata później) dla dużych rzadkich tablic, zobacz sklearn.metrics.pairwise.pairwise_distances.html w scikit-learn.
den

Odpowiedzi:

26

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

bayerj
źródło
3
Wydaje mi się, że odpowiedź Bayerja wymaga drobnych modyfikacji w celu dopasowania do formuły, na wypadek, gdyby ktoś inny jej potrzebował:K = scipy.exp(-pairwise_dists**2 / s**2)
chloe
Jeśli ktoś jest ciekawy, zastosowany algorytm pdistjest bardzo prosty: jest to po prostu pętla zaimplementowana w języku C, która bezpośrednio oblicza odległości w oczywisty sposób , przy czym pętla jest wykonywana tutaj ; żadna fantazyjna wektoryzacja ani nic poza tym, co kompilator może wykonać automatycznie.
Dougal,
11

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)
tenedor
źródło
1
Lub po prostu pairwise_sq_dists = cdist(X, X, 'sqeuclidean')co daje to samo.
user1721713
5

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

spetz911
źródło
Jest to metoda sugerowana przez kardynała w komentarzach, którą można nieco przyspieszyć za pomocą operacji inplace. To jak scikit-learn to robi , ze na einsumwezwanie Dla Państwa X2.
Dougal,
4
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.

Jan
źródło
Witamy na naszej stronie! Kładziemy nieco inny nacisk na przepełnienie stosu, ponieważ generalnie mniej koncentrujemy się na kodzie, a bardziej na podstawowych pomysłach, więc warto zanotować kod lub podać krótki pomysł, jakie są jego kluczowe pomysły, ponieważ niektóre zrobiły to inne odpowiedzi. Pomogłoby to wyjaśnić, jak Twoja odpowiedź różni się od innych.
Silverfish
Będzie to znacznie wolniejsze niż inne odpowiedzi, ponieważ używa pętli Python zamiast wektoryzacji.
Dougal,
-1

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

def GaussianKernel(v1, v2, sigma):
    return exp(-norm(v1-v2, 2)**2/(2.*sigma**2))
Jądro
źródło
3
Witamy na stronie @Kernel. Możesz wyświetlić matematykę, umieszczając wyrażenie między znakami $ i używając składni podobnej do LateX. Możesz wyświetlać kod (z podświetlaniem składni) wcinając linie o 4 spacje. Zobacz pomoc dotyczącą edycji przeceny w celu uzyskania wskazówek na temat formatowania, a często zadawane pytania w celu uzyskania bardziej ogólnych wskazówek .
Antoine Vernet,
1
Czy to nie odbija się echem w pytaniu?
whuber