Support Vector Machine (SVM)

The maximum margin classifier as the solution of a quadratic programming problem

First cook up some data to test it on:

In [37]:
from numpy import *
d = 2
n = 30
x = random.rand(d,n)
X = empty((d+1,n))
X[0,:] = 1
X[1:,:] = x
W = array([-1,-2,3.])
y = sign(dot(W,X))
# let's open up a little space between the classes
halfgap = 0.05
X[-1,y>0] += halfgap
X[-1,y<0] -= halfgap
xp = X[1:,y>0]
xm = X[1:,y<0] 
In [38]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [39]:
subplot(111,aspect=1)
plot(xp[0],xp[1],'ro')
plot(xm[0],xm[1],'bo');
# We are given X and y ONLY!!

def drawline(W,color='g',alpha=1.):
    plot([0,1],[-W[0]/W[2],-(W[0]+W[1])/W[2]],color=color,alpha=alpha)
    
In [40]:
X
Out[40]:
array([[  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00],
       [  6.85486133e-01,   3.88492985e-01,   1.26410428e-01,
          1.50498124e-01,   2.10074578e-01,   4.92595297e-01,
          1.23332336e-04,   3.33778422e-01,   3.95830835e-01,
          8.52578740e-01,   6.88937643e-01,   5.15304348e-02,
          3.23396808e-01,   1.26657689e-01,   4.30538852e-01,
          5.49957649e-02,   6.32539138e-01,   9.06908651e-01,
          6.01296336e-01,   6.30430822e-01,   2.18189820e-01,
          3.64665393e-01,   1.83755065e-01,   9.04610689e-01,
          6.97742641e-01,   2.64643028e-01,   1.90625224e-01,
          7.08726761e-01,   5.17495578e-01,   1.39016345e-01],
       [  9.17538250e-01,   9.98623036e-01,   7.45653046e-02,
          8.61462459e-01,   1.75747473e-01,   5.63798470e-01,
          4.26698374e-01,   1.94344018e-01,   6.52665420e-01,
          8.24195611e-01,   4.87537967e-01,   6.46644365e-01,
          7.04684980e-01,   7.58255225e-01,   3.51230897e-01,
          8.46642311e-01,   8.89428090e-02,   5.98895068e-01,
          1.61744127e-01,   4.15594137e-01,   7.40471110e-01,
          8.93579322e-01,   9.68452316e-01,   4.04312191e-01,
          5.97680161e-01,   2.38221831e-01,   4.38072827e-02,
          2.43206359e-01,   8.82048662e-01,   9.41318340e-01]])
In [47]:
y
Out[47]:
array([ 1.,  1., -1.,  1., -1., -1.,  1., -1.,  1., -1., -1.,  1.,  1.,
        1., -1.,  1., -1., -1., -1., -1.,  1.,  1.,  1., -1., -1., -1.,
       -1., -1.,  1.,  1.])

Prep the parameters for cvxopt solver

In [41]:
from cvxopt import matrix,solvers

P = eye(d+1)
P[0,0] = 0
P = matrix(P)

q = zeros(d+1)
q = matrix(q)

G = (-X*y).T
G = matrix(G)

h = -ones(n)
h = matrix(h)

Solve the problem

In [42]:
sol=solvers.qp(P,q,G,h)
     pcost       dcost       gap    pres   dres
 0:  2.4029e+00  2.4172e+01  7e+01  2e+00  9e+00
 1:  1.0059e+01  1.6840e+01  3e+01  7e-01  4e+00
 2:  2.1080e+01  7.0171e+01  3e+01  7e-01  3e+00
 3:  6.8257e+01  1.0140e+02  4e+01  4e-01  2e+00
 4:  1.1457e+02  1.1601e+02  3e+00  2e-02  1e-01
 5:  1.1673e+02  1.1675e+02  4e-02  2e-04  1e-03
 6:  1.1676e+02  1.1676e+02  4e-04  2e-06  1e-05
 7:  1.1676e+02  1.1676e+02  4e-06  2e-08  1e-07
 8:  1.1676e+02  1.1676e+02  4e-08  2e-10  1e-09
Optimal solution found.
In [43]:
type(sol)
Out[43]:
dict
In [44]:
sol.keys()
Out[44]:
dict_keys(['s', 'status', 'relative gap', 'primal slack', 'primal objective', 'dual slack', 'dual objective', 'iterations', 'z', 'gap', 'dual infeasibility', 'primal infeasibility', 'y', 'x'])
In [45]:
W = array(sol['x']).reshape(d+1)

Plot the solution

In [46]:
subplot(111,aspect=1)
plot(xp[0],xp[1],'ro')
plot(xm[0],xm[1],'bo');
# We are given X and y ONLY!!

def drawline(W,color='g',alpha=1.):
    plot([0,1],[-W[0]/W[2],-(W[0]+W[1])/W[2]],color=color,alpha=alpha)

drawline(W,color='b')

Nice!

In [ ]: