SVM - How can I vectorize a kernalized gram matrix?

Question

I implemented a support vector machine in python using the cvxopt qp solver where I need to compute a gram matrix of two vectors with a kernel function at each element. I implemented it correctly using for loops but this strategy is computationally intensive. I would like to vectorize the code.

Example:

Here is what I have written:

K = np.array( [kernel(X[i], X[j],poly=poly_kernel) 
     for j in range(m)
     for i in range(m)]).reshape((m, m))

How can I vectorize the above code without for loops to achieve the same result faster?

The kernel function computes a gaussian kernel.

Here is a quick explanation of an svm with kernel trick. Second page of this explains the problem.

Here is my full code for context.

EDIT: Here is a quick code snippet that runs what I need to vectorized in an unvectorized form

from sklearn.datasets import make_gaussian_quantiles;
import numpy as np;


X,y = make_gaussian_quantiles(mean=None, cov=1.0, n_samples=100, n_features=2, n_classes=2, shuffle=True, random_state=5);

m = X.shape[0];


def kernel(a,b,d=20,poly=True,sigma=0.5):
    if (poly):
        return np.inner(a,b) ** d;
    else:
        return np.exp(-np.linalg.norm((a - b) ** 2)/sigma**2)

# Need to vectorize these loops

K = np.array([kernel(X[i], X[j],poly=False) 
    for j in range(m)
    for i in range(m)]).reshape((m, m))

Thanks!

Please make it easy for those trying to help you out by including the relevant bits (the kernel function in this case) and some mock data, so they can directly run your example. — Paul Panzer, Dec 22 '17 at 21:57

Paul Panzer · Accepted Answer · 2017-12-23T06:50:33.283

3

Here is a vectorized version. The non poly branch comes in two variants a direct one and a memory saving one in case the number of features is large:

from sklearn.datasets import make_gaussian_quantiles;
import numpy as np;


X,y = make_gaussian_quantiles(mean=None, cov=1.0, n_samples=100, n_features=2, n_classes=2, shuffle=True, random_state=5);
Y,_ = make_gaussian_quantiles(mean=None, cov=1.0, n_samples=200, n_features=2, n_classes=2, shuffle=True, random_state=2);

m = X.shape[0];
n = Y.shape[0]

def kernel(a,b,d=20,poly=True,sigma=0.5):
    if (poly):
        return np.inner(a,b) ** d;
    else:
        return np.exp(-np.linalg.norm((a - b) ** 2)/sigma**2)

# Need to vectorize these loops

POLY = False
LOW_MEM = 0

K = np.array([kernel(X[i], Y[j], poly=POLY) 
              for i in range(m)
              for j in range(n)]).reshape((m, n))

def kernel_v(X, Y=None, d=20, poly=True, sigma=0.5):
    Z = X if Y is None else Y
    if poly:
        return np.einsum('ik,jk', X, Z)**d
    elif X.shape[1] < LOW_MEM:
        return np.exp(-np.sqrt(((X[:, None, :] - Z[None, :, :])**4).sum(axis=-1)) / sigma**2)
    elif Y is None or Y is X:
        X2 = X*X
        H = np.einsum('ij,ij->i', X2, X2) + np.einsum('ik,jk', X2, 3*X2) - np.einsum('ik,jk', X2*X, 4*X)
        return np.exp(-np.sqrt(np.maximum(0, H+H.T)) / sigma**2)
    else:
        X2, Y2 = X*X, Y*Y
        E = np.einsum('ik,jk', X2, 6*Y2) - np.einsum('ik,jk', X2*X, 4*Y) - np.einsum('ik,jk', X, 4*Y2*Y)
        E += np.add.outer(np.einsum('ij,ij->i', X2, X2), np.einsum('ij,ij->i', Y2, Y2))
        return np.exp(-np.sqrt(np.maximum(0, E)) / sigma**2)

print(np.allclose(K, kernel_v(X, Y, poly=POLY)))

edited Dec 23 '17 at 06:50

answered Dec 23 '17 at 01:02

Paul Panzer

51,835
3
54
99

Thanks! Could you explain what the einsums are doing? – random_0620 Dec 23 '17 at 02:24
1

`einsum('ik,jk', X, X)` multiplies elements of `X` with elements of `X` the output will have an axis like the first of the first `X` because `i` in the spec string is unique and an axis like the first of the second 'X' because 'j' is unique. So with respect to these two axes it behaves like an outer product. The second axes of the two factors are not expanded and the individual products are summed along this axis because 'k' occurs in both specs. That is the Einstein convention hence "einsum", so with respect to the second axis it behaves like an inner product. `out_ij = sum_k X_ik X_jk` – Paul Panzer Dec 23 '17 at 02:35
If there is also an output spec as in `np.einsum('ij,ij->i', X2, X2)` in this case it means `i` is not summed over but also not outer expanded, so `out_i = sum_j X2_ij X2_ij` I hope with that you can figure out the others yourself. – Paul Panzer Dec 23 '17 at 02:39
Thanks for the explanation. Why did you multiply 3*X2 and 4*X? – random_0620 Dec 23 '17 at 04:02
Just some trickery `(a-b)^4=a^4 - 4a^3 b + 6 a^2 b^2 - 4 a b^3 + b^4`. because we do this "outer" and the vectors for a and b are both `X` we can calculate the first half of the terms and the second half will then just be the transpose. Since the number of terms is odd we need to split the middle one in half, thus `3`. – Paul Panzer Dec 23 '17 at 04:56
So later in the code I have to use the kernel on two different a and b so would this method not work in that case? – random_0620 Dec 23 '17 at 05:14
Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/161844/discussion-between-s-kirkiles-and-paul-panzer). – random_0620 Dec 23 '17 at 06:33

SVM - How can I vectorize a kernalized gram matrix?

1 Answers1