How to calculate a Gaussian kernel effectively in numpy [closed]

Question

I have a numpy array with m columns and n rows, the columns being dimensions and the rows datapoints.

I now need to calculate kernel values for each combination of data points.

For a linear kernel $K(\mathbf{x}_i,\mathbf{x}_j) = \langle \mathbf{x}_i,\mathbf{x}_j \rangle$ I can simply do dot(X,X.T)

How can I effectively calculate all values for the Gaussian Kernel $K(\mathbf{x}_i,\mathbf{x}_j) = \exp{-\frac{\|\mathbf{x}_i-\mathbf{x}_j\|_2^2}{s^2}}$ with a given s?

Answer 1

I think the main problem is to get the pairwise distances efficiently. Once you have that the rest is element wise.

To do this, you probably want to use scipy. The function scipy.spatial.distance.pdist does what you need, and scipy.spatial.distance.squareform will possibly ease your life.

So if you want the kernel matrix you do

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)

Documentation can be found here. 

Answer 2

As a small addendum to bayerj's answer, scipy's pdist function can directly compute squared euclidean norms by calling it as pdist(X, 'sqeuclidean'). The full code can then be written more efficiently as

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)

Answer 3

You can also write square form by hand:

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 but this works 30% slower

Answer 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)

which is equal to

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

You can effectively calculate the RBF from the above code note that the gamma value is 1, since it is a constant the s you requested is also the same constant.

Answer 5

I think this will help:

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