In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
%matplotlib inline
In [2]:
# Make an artificial dataset.
Npos = Nneg = 50
N = Npos + Nneg
X = np.vstack([np.random.randn(Npos, 2) - 3.0, np.random.randn(Nneg, 2) - 1.0, np.random.randn(Npos, 2) + 1.0, np.random.randn(Nneg, 2) + 3.0])
y = np.hstack([+np.ones(Npos), -np.ones(Nneg), +np.ones(Npos), -np.ones(Nneg)])
In [3]:
# Visualize
(fig, axs) = plt.subplots(1,1)
ax = axs
ax.scatter(X[y==+1, 0], X[y==+1, 1], color='red', marker='.')
ax.scatter(X[y==-1, 0], X[y==-1, 1], color='blue', marker='.')
ax.axis('equal')
ax.grid()
In [4]:
def SolveKernelSVM(K, y, C):
y1 = y.reshape(-1, 1)
Y = np.dot(y1, y1.T)
G = K * Y
alp0 = np.zeros(y.size)
# Objective function
func = lambda alpha : 0.5 * alpha.dot(G).dot(alpha) - alpha.sum()
func_deriv = lambda alpha : alpha.dot(G) - np.ones_like(alpha)
# Constraints
cons = [
{'type':'eq', 'fun':lambda alpha : alpha.dot(y), 'jac':lambda alpha:y },
{'type':'ineq', 'fun':lambda alpha : alpha, 'jac':lambda alpha:np.eye(alpha.size) },
{'type':'ineq', 'fun':lambda alpha : C-alpha, 'jac':lambda alpha:-np.eye(alpha.size) },
]
# Solve using SLSQP optimization algorithm.
res = scipy.optimize.minimize(func, alp0, method='SLSQP', constraints=cons, jac=func_deriv, options={'disp':True, 'maxiter':int(1e4)})
alp = res.x
return alp
In [5]:
def DistanceMatrix(X, Y):
''' Compute distance matrix among data matrices X, Y, whose row-vector is representating a sample.'''
assert( X.shape[1] == Y.shape[1] ) # Check if the datasets have the same dimension.
nx, ny = [ (M**2).sum(axis=1).reshape(-1,1) for M in [X, Y] ]
In, Im = [ np.ones(M.shape[0]).reshape(-1,1) for M in [X, Y] ]
D = np.dot(nx, Im.T) + np.dot(In, ny.T) - 2 * np.dot(X, Y.T)
return D
In [6]:
def KernelMatrixLinear(X, Y):
K = np.dot(X, Y.T)
return K
In [7]:
def KernelMatrixPolynomial(X, Y, p):
K = np.power(np.dot(X, Y.T) + 1, p)
return K
In [8]:
def KernelMatrixRBF(X, Y, gamma):
D = DistanceMatrix(X, Y)
K = np.exp(-gamma * D)
return K
In [9]:
def KernelMatrix(X, Y):
gamma = 0.1
return KernelMatrixRBF(X, Y, gamma)
In [10]:
def PlotHeatMap(ax, KernelMatrix, alp, y, b):
smin, smax = -5, 5 # definition of range of samples.
nMesh = 101 # The number of meshes to divide.
Xm, Ym = np.meshgrid( *[np.linspace(smin, smax, nMesh) for _ in range(2)] )
XX = np.zeros((2, nMesh, nMesh))
XX[0], XX[1] = Xm, Ym
X2 = XX.reshape((2,-1)).T # Feature vectors.
K2 = KernelMatrix(X, X2) # Kernel Matrix
Y2 = np.dot(alp*y, K2) + b # Prediction
Zm = Y2.reshape((nMesh, nMesh)) # Reshape into the mesh-grid style.
ax.pcolor(Xm, Ym, Zm, cmap='jet')
cont = ax.contour(Xm, Ym, Zm, levels=[-1.0, 0.0, 1.0], colors='k')
cont.clabel(fmt='%1.1f', fontsize=14, inline=1)
ax.scatter(X[y==+1, 0], X[y==+1, 1], color='red', marker='.')
ax.scatter(X[y==-1, 0], X[y==-1, 1], color='blue', marker='.')
ax.grid()
ax.set_xlim(smin, smax)
ax.set_ylim(smin, smax)
In [11]:
%%time
Cs = [1e-1, 1e0, 1e1]
gammas = [1e-1, 1e0, 1e1]
nCols = len(Cs)
nRows = len(gammas)
figSize = 5
(fig, axs) = plt.subplots(nRows, nCols, figsize=(figSize*nCols, figSize*nRows) )
for i in range(nRows):
for j in range(nCols):
ax = axs[i, j]
gamma = gammas[i]
C = Cs[j]
# Definition of a kernel function.
KernelMatrix = lambda X, Y : KernelMatrixRBF(X, Y, gamma)
# Train SVM and compute parameters.
K = KernelMatrix(X, X)
alp = SolveKernelSVM(K, y, C)
k = np.argmin(np.abs(C/2 - alp))
b = y[k] - (y*alp).dot(K[k])
# Plot
ax.set_title("$\gamma$ = %.1f, C = %.1f" % (gamma, C) )
PlotHeatMap(ax, KernelMatrix, alp, y, b)