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
            Gradient evaluations: 24
CPU times: user 12.9 s, sys: 168 ms, total: 13.1 s
Wall time: 8.19 s