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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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