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

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?

• Well if you don't care too much about a factor of two increase in computations, you can always just do $\newcommand{\m}{\mathbf} \m S = \m X \m X^T$ and then $K(\m x_i, \m x_j ) = \exp( - (S_{ii} + S_{jj} - 2 S_{ij})/s^2 )$ where, of course, $S_{ij}$ is the $(i,j)$th element of $\m S$. This is probably not the most numerically stable, either, though. Commented Sep 20, 2011 at 13:46
• (Years later) for large sparse arrays, see sklearn.metrics.pairwise.pairwise_distances.html in scikit-learn . Commented Dec 20, 2015 at 14:27

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
pairwise_dists = squareform(pdist(X, 'euclidean'))
K = scip.exp(-pairwise_dists ** 2 / s ** 2)


Documentation can be found here

• If anyone is curious, the algorithm used by pdist is very simple: it's just a C-implemented loop that directly computes distances in the obvious way, the looping being done here; no fancy vectorization or anything beyond whatever the compiler can accomplish automatically. Commented Mar 25, 2017 at 2:27

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
pairwise_sq_dists = squareform(pdist(X, 'sqeuclidean'))
K = scip.exp(-pairwise_sq_dists / s**2)

• Or simply pairwise_sq_dists = cdist(X, X, 'sqeuclidean') which gives the same. Commented Jan 6, 2018 at 2:09

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

• This, which is the method suggested by cardinal in the comments, could be sped up a bit by using inplace operations. It's how scikit-learn does it, with an einsum call for your X2. Commented Mar 25, 2017 at 2:18
• You wrote: K0 = X2 + X2.T - 2 * X * X.T - how does it can work with X and X.T having different dimensions? Commented Jan 16, 2022 at 0:17
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.

• Welcome to our site! We have a slightly different emphasis to Stack Overflow, in that we generally have less focus on code and more on underlying ideas, so it might be worth annotating your code or giving a brief idea what the key ideas to it are, as some of the other answers have done. That would help explain how your answer differs to the others. Commented Mar 24, 2017 at 1:41
• This will be much slower than the other answers because it uses Python loops rather than vectorization. Commented Mar 25, 2017 at 2:13

I think this will help:

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

• Welcome to the site @Kernel. You can display mathematic by putting the expression between \$ signs and using LateX like syntax. And you can display code (with syntax highlighting) by indenting the lines by 4 spaces. See the markdown editing help for formatting guidelines, and the faq for more general ones. Commented Jan 11, 2013 at 15:24
• Doesn't this just echo what is in the question?
– whuber
Commented Jan 11, 2013 at 16:09