Computation of kernel matrix using radial basis kernel in svm

The kernel matrix $K$ or Gram Matrix is a positive-semidefinite symmetric matrix that has the form: $$ K(x_1,x_2,\dots,x_n) = \begin{vmatrix} \langle x_1,x_1 \rangle & \langle x_1,x_2 \rangle & \dots & \langle x_1,x_n \rangle \\ \langle x_2,x_1 \rangle & \langle x_2,x_2 \rangle & \dots & \langle x_2,x_n \rangle \\ \vdots & \vdots & \ddots &\vdots \\ \langle x_n,x_1 \rangle & \langle x_n,x_2 \rangle & \dots & \langle x_n,x_n \rangle \end{vmatrix} $$

Where $\langle a,b \rangle$ detones the inner product/dot product of vectors $a$ and $b$ and $\{x_1,x_2,\dots,x_n\}$ is a set of vectors (in this case, training samples).

In order to calculate the Gramian Matrix you will have to calculate the Inner Product using the Kernel Function. For a RBF kernel function $\kappa_{RBF}$ this can be done by

$$K_{ij} = \kappa_{RBF}(x_i,x_j) = e^{\gamma D_{ist}(x_i,x_j)^2} $$

where $\gamma$ is a function hyperparameter, $K_{ij}$ is the element in row $i$ and column $j$ of the matrix $K$ and $ D_{ist}(x_i,x_j)$ is some distance between two vector measured in some vector space.

Usually, the distance measure used is the $L_2$ norm or euclidean distance. Remember that when computing the Kernel function you will use the square $L_2$ norm, so there is no need to take the square root.

There is a open library named pykernels that has implementations to dozens of kernel functions. This here is the implementation for the RBF (or Gaussian) Kernel:

class RBF(Kernel):
    """
    Radial Basis Function kernel, defined as unnormalized Gaussian PDF
        K(x, y) = e^(-g||x - y||^2)
    where:
        g = gamma
    """

    def __init__(self, gamma=None):
        self._gamma = gamma

    def _compute(self, data_1, data_2):
        if self._gamma is None:
            # libSVM heuristics
            self._gamma = 1./data_1.shape[1]

        dists_sq = euclidean_dist_matrix(data_1, data_2)
        return np.exp(-self._gamma * dists_sq)

    def dim(self):
        return np.inf

it uses the square distance define as:

def euclidean_dist_matrix(data_1, data_2):
    """
    Returns matrix of pairwise, squared Euclidean distances
    """
    norms_1 = (data_1 ** 2).sum(axis=1)
    norms_2 = (data_2 ** 2).sum(axis=1)
    return np.abs(norms_1.reshape(-1, 1) + norms_2 - 2 * np.dot(data_1, data_2.T))

The class is defined as a subclass of the Kernel class and has some methods to deal with kernel combinations.

Note that this is not the fastest way to implement euclidean distance, you can do better but this is a simple form to do it. also that "np.abs" should not be needed

Note 2 that this SO Question has a lot of optimized version of RBF Kernel