Q:

RBF kernel algorithm Python

I have this algorithm to compute the RBF kernel and it seems to work just fine. But I would like to understand what kind of operations are involved, for example:

What are the trnorms vectors? What are they for?
What is the meaning of creating the matrices k1 and k2?
Is this algorithm any different from the sklearn implementation?

def compute_RBF(mat1, mat2, sigma):

    trnorms1 = np.mat([(v * v.T)[0, 0] for v in mat1]).T
    trnorms2 = np.mat([(v * v.T)[0, 0] for v in mat2]).T

    k1 = trnorms1 * np.mat(np.ones((mat2.shape[0], 1), dtype=np.float64)).T

    k2 = np.mat(np.ones((mat1.shape[0], 1), dtype=np.float64)) * trnorms2.T

    k = k1 + k2

    k -= 2 * np.mat(mat1 * mat2.T)

    k *= - 1./(2 * np.power(sigma, 2))

    return np.exp(k)

A:

Say that mat1 is $n \times d$ and mat2 is $m \times d$.
Recall that the Gaussian RBF kernel is defined as $k(x, y) = \exp\left( - \frac{1}{2 \sigma^2} \lVert x - y \rVert^2 \right)$.
But we can write $\lVert x - y \rVert^2$ as $(x - y)^T (x - y) = x^T x + y^T y - 2 x^T y$. The code uses this decomposition.
First, the trnorms1 vector stores $x^T x$ for each input $x$ in mat1, and trnorms2 stores $y^T y$ for each $y$ in mat2.
Then, the k1 matrix is obtained by multiplying the $n \times 1$ matrix of $x^T x$ entries by a $1 \times m$ matrix of ones, getting an $n \times m$ matrix with $x^T x$ entries repeated across the rows, so that k1[i, j] is $x_i^T x_i$.
The next line does basically the same thing for the $y$ norms repeated across columns, getting an $n \times m$ matrix with k2[i, j] of $y_j^T y_j$.
k is then their sum, so that k[i, j] is $x_i^T x_i + y_j^T y_j$. The next line then subtracts twice the product of the data matrices, so that k[i, j] becomes $x_i^T x_i + y_j^T y_j - 2 x_i^T y_j = \lVert x_i - y_j \rVert^2$.
Then, the code multiplies by $\frac{-1}{2 \sigma^2}$ and finally takes the elementwise $\exp$, getting out the Gaussian kernel.
If you dig into the scikit-learn implementation, it's exactly the same, except:

It's parameterized instead with $\gamma = \frac{1}{2 \sigma^2}$.
It's written in much better Python, not wasting memory all over the place and doing computations in a needlessly slow way.
It's broken up into helper functions.

But, algorithmically, it's doing the same basic operations.

