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

In [2]:
def GenerateDummyData(nPos=100, nNeg=100):
N = nPos + nNeg
Xpos = np.random.normal(+6, 1, (nPos, 2))
Xneg = np.random.normal(+3, 1, (nNeg, 2))
Ypos =  np.ones(nPos)
Yneg = -np.ones(nNeg)
X = np.r_[Xpos, Xneg]
Y = np.r_[Ypos, Yneg]
perm = np.random.permutation(N)
X = X[perm,:]
Y = Y[perm]
return X, Y

In [3]:
def Primal(X, Y, lam, w):
r = lam * w.dot(w) / 2
l = 1 - Y * X.dot(w)
l = ( l * (l > 0) ).sum() / l.size
return r + l

def Dual(X, Y, lam, w, alp):
r = -lam * w.dot(w) / 2
l = (Y * alp).sum() / alp.size
return r + l

def SolveSVM(X, Y, lam, T=10):
N, d = X.shape

alpha = np.zeros(N)
w = np.zeros(d) # in general cace s.t. alpha != 0 , w = alpha.dot(X) / (lam*N)
p, d, sv = [ np.zeros(T) for _ in range(3) ] # to save learning history.

for t in range(T):
perm = np.random.permutation(N)
for i in perm:
x, y, alp_old = X[i], Y[i], alpha[i]
alp_new = y * np.clip( y*alp_old + lam*N*( 1-y*x.dot(w) ) / x.dot(x), 0, 1 )
alpha[i] = alp_new
w += (alp_new - alp_old) * x / (lam*N)

p[t] = Primal(X,Y,lam,w)
d[t] = Dual(X,Y,lam,w,alpha)
sv[t] = (alpha != 0).sum()

return w, alpha, p, d, sv

In [4]:
%%time
np.random.seed(2018)
lam = 1e-3
I = int(3e3) # the number of iteration.

X, Y = GenerateDummyData()
N = Y.size
X = np.c_[X, np.ones(X.shape[0]) ] # Add bias

w,alpha,p,d,sv = SolveSVM(X, Y, lam, I)

CPU times: user 9.82 s, sys: 4 ms, total: 9.82 s
Wall time: 9.82 s

In [5]:
# Compute separating line.
x = X[:,0]
x = np.linspace(x.min(), x.max())
y = -( w[0]*x + w[2] ) / w[1]

# Plot
fig = plt.figure(figsize=(10, 4))
ax1 = plt.subplot2grid((2,3), (0,0), rowspan=2)
ax2 = plt.subplot2grid((2,3), (0,1))
ax3 = plt.subplot2grid((2,3), (1,1))
ax4 = plt.subplot2grid((2,3), (0,2), rowspan=2)

# Plot
ax = ax1
ax.set_title('Data distribution and SVs.')
ax.plot(x, y, '--')
ax.plot(*X[alpha != 0,:2].T, 'go', alpha=0.3)
ax.plot(*X[Y==+1,:2].T, 'r+')
ax.plot(*X[Y==-1,:2].T, 'b+')

ax = ax2
ax.set_title('Objective value')
ax.plot(p, label='primal')
ax.plot(d, label='dual')
ax.set_xscale('log')
ax.set_yscale('log')
ax.legend()

ax = ax3
ax.set_title('Duality gap')
ax.plot(p-d)
ax.set_xscale('log')
ax.set_yscale('log')

ax = ax4
ax.set_title('# of Supporting vectors')
ax.plot(sv)

for ax in [ax1,ax2,ax3,ax4] : ax.grid()
fig.tight_layout()