In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
%matplotlib inline

In [2]:
# Make an artificial dataset.
Npos = Nneg = 200
N = Npos + Nneg

X = np.vstack([np.random.randn(Npos, 2) + 1.0, np.random.randn(Nneg, 2) - 1.0])
y = np.hstack([+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 SolveLinearSVM1(X, y, C):
Y = np.dot( y.reshape(-1, 1), y.reshape(1, -1) )
K = np.dot(X, X.T)
G = K * Y
alp0 = np.zeros_like(y)
# 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 PlotLinearHeatMap(ax, X, y, w, 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)] )
Zm = Xm*w[0] + Ym*w[1] + b

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 [6]:
%%time

C = 1.0

# Solve a SVM problem.
alp = SolveLinearSVM1(X, y, C)

# Obtain a normal vector and a bias parameter.
w = (alp*y).dot(X)
i = np.argmin( np.abs(C/2 - alp) ) # Find an index i, s.t. 0 < alp[i] < C, ==> min(| C/2 - alp[i] |)
b = y[i] - w.dot(X[i]) # From KKT-condition of SVM problem, the bias parameter is derived the same as this.

# Plot
(fig, ax) = plt.subplots(1,1)
ax.set_title('C = %.1f' % C)
PlotLinearHeatMap(ax, X, y, w, b)

Optimization terminated successfully.    (Exit mode 0)
Current function value: -85.3362325706
Iterations: 25
Function evaluations: 28