Skip to content

Instantly share code, notes, and snippets.

@evancchow
Created December 28, 2013 10:56
Faulty (and also really messy) code for finding max-margin decision boundary given datapoints -- that is, the max-margin line that separates datapoints into 2 classes given a set of 2 dimensional points. Problems occur with w*, which is wopt in the program below.
from cvxopt import solvers, matrix
import numpy as np
import matplotlib.pyplot as plt
# specify # of dimensions
dim = 2
""" Initialize data and plot the max-margin decision boundary.
You already know the data / max-margin decision boundary
(from working it out on paper) so we're going to plot it now
and check our results later from the actual algorithm implementation."""
# x: datapoints; y: labels; b, w: w0, w1, w2
x = np.array([[1,1],[2,2],[0,0]])
y = np.array([1,1,-1])
b = -1
w = np.array([1,1])
length = len(x) # number of data points
# plot points. (1) = green; (-1) = red.
for i in xrange(length):
if y[i] == 1:
plt.scatter(x[i][0],x[i][1], color='g')
else:
plt.scatter(x[i][0],x[i][1], color='r')
# plot title, max-margin line
plt.plot([0,1],[1,0],'m')
plt.title('Green: 1; Red: -1')
# copies of original numpy arrays; convert to cvxopt matrices
xnp = x
ynp = y
x = matrix(x)
y = matrix(y)
""" Begin algorithm implementation """
# Set up a couple matrices for the quadratic program.
# in order: solver.qp(p, q, G, h, A, b)
p = matrix(np.identity(2))
q = matrix(np.zeros([2,1]))
# Algorithm
r = 0
for i in xrange(length):
mag = np.linalg.norm(x[i])
if r < mag:
r = mag # should be 2.0
poslim = 2
neglim = -1 * poslim # for b range -- increase range later
# wopt = None
bopt = None
# construct X
X = np.zeros(xnp.T.shape)
for n in xrange(dim):
for idx, item in enumerate(ynp):
X[n][idx] = xnp[idx][n] * ynp[idx]
X = matrix(X.T)
# local function to test if variable is defined
def varDefined(var):
try:
var
except NameError:
return 0 # not defined
else:
return 1 # defined
# start looping
counter = 0 # track loop number
for b in xrange(neglim, poslim, 1):
# construct c
ones = np.zeros([length,1])
ones[:] = 1
product = (ynp * b).reshape(-1,1)
c = matrix(ones + product)
# solve quadratic problem
w = solvers.qp(p, q, X, c)['x']
print "Loop #%s!" % (counter + 1)
print "W: ", w
if (counter > 1):
print "Wopt: ", wopt
# conditions -- need to refactor this
if ((varDefined(w) == 1 and varDefined(wopt) == 0) or
(varDefined(w) == 1 and np.linalg.norm(w) < np.linalg.norm(wopt))):
wopt = w
bopt = b
counter+=1
if (varDefined(wopt) == 0):
raise Exception('Constraints not satisfiable!')
elif (np.linalg.norm(wopt) > (poslim / float(r))):
raise Exception('Bounding assumption of optimal w violated')
print wopt, bopt
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment