s_kirkiles - 6 months ago 85

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

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!

Recommended for you: Get network issues from **WhatsUp Gold**. **Not end users.**

Answer Source

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

Recommended from our users: **Dynamic Network Monitoring from WhatsUp Gold from IPSwitch**. ** Free Download**