In [2]:

    

%matplotlib inline 

from __future__ import division

import sys
if '../' not in sys.path: sys.path.append("../")

import numpy as np
import amnet
import cvxpy

from numpy.random import seed, randn
import matplotlib.pyplot as plt


    

In [3]:

    

# Replicate BV Fig 8.10 using CVXPY
# data generation
n = 2
N = 10
M = 10
seed(0)
Y = np.vstack((
    np.hstack((1.5 + 0.9*randn(1, int(0.6*N)), 1.5 + 0.7*randn(1, int(0.4*N)))),
    np.hstack((2*(randn(1, int(0.6*N)) + 1), 2*(randn(1, int(0.4*N)) - 1)))
))
X = np.vstack((
    np.hstack((-1.5 + 0.9*randn(1,int(0.6*M)),  -1.5 + 0.7*randn(1, int(0.4*M)))),
    np.hstack((2*(randn(1, int(0.6*M)) - 1), 2*(randn(1, int(0.4*M)) + 1))),
))
T = np.array([[-1, 1], [1, 1]])
Y = np.dot(T, Y)
X = np.dot(T, Y)


    

In [6]:

    

# solution via CVXPY
a = cvxpy.Variable(n)
b = cvxpy.Variable(1)
u = cvxpy.Variable(N)
v = cvxpy.Variable(M)
obj = cvxpy.Minimize(sum(u) + sum(v))
cons = [X.T * a - b >= 1 - u,
        Y.T * a - b <= -(1 - v),
        u >= 0,
        v >= 0,
        cvxpy.norm1(u) >=1,
        cvxpy.norm1(v) >=1]
prob = cvxpy.Problem(obj, cons)

result = prob.solve()
print result


    

    




---------------------------------------------------------------------------
DCPError                                  Traceback (most recent call last)
<ipython-input-6-47816946dcf7> in <module>()
     13 prob = cvxpy.Problem(obj, cons)
     14 
---> 15 result = prob.solve()
     16 print result

/home/spc923/.local/lib/python2.7/site-packages/cvxpy/problems/problem.pyc in solve(self, *args, **kwargs)
    207             return func(self, *args, **kwargs)
    208         else:
--> 209             return self._solve(*args, **kwargs)
    210 
    211     @classmethod

/home/spc923/.local/lib/python2.7/site-packages/cvxpy/problems/problem.pyc in _solve(self, solver, ignore_dcp, warm_start, verbose, parallel, **kwargs)
    281                       "Solving a convex relaxation.")
    282             else:
--> 283                 raise DCPError("Problem does not follow DCP rules.")
    284 
    285         if solver == s.LS:

DCPError: Problem does not follow DCP rules.

In [4]:

    

# graph CVXPY solution
t_min = min(np.hstack((X[0,:], Y[0,:])))
t_max = max(np.hstack((X[0,:], Y[0,:])))
tt = np.linspace(t_min-1, t_max+1, 100)
p = np.ravel(-a.value[0]/a.value[1]*tt + b.value/a.value[1])
p1 = np.ravel(-a.value[0]*tt/a.value[1] + (b.value+1)/a.value[1])
p2 = np.ravel(-a.value[0]*tt/a.value[1] + (b.value-1)/a.value[1])
plt.plot(X[0,:], X[1,:], 'o', Y[0,:], Y[1,:], 'o')
plt.plot(tt, p, '-r', tt, p1, '--r', tt, p2, '--r')
plt.title('Appriximate Linear Discrimination')
plt.axis('equal')
plt.xlim(-10, 10)
plt.ylim(-10, 10)


    

    
Out[4]:





(-10, 10)






    

In [7]:

    

# solution via AMNET
x = amnet.Variable(n + 1 + N + M, name='x')
a = x[0:n]
b = x[n:n+1]
u = x[n+1:n+1+N]
v = x[n+1+N:n+1+N+M]
Em = np.ones((M,1))
En = np.ones((N,1))
assert len(a) == n and len(b) == 1 and len(u) == N and len(v) == M

obj = amnet.opt.Minimize(amnet.atoms.add_all(u) + amnet.atoms.add_all(v))
cons = [X.T * a - Em * b >= 1 - u,
        Y.T * a - En * b <= -(1 - v),
        u >= 0,
        v >= 0,
       amnet.atoms.norm1(u)>=1,
       amnet.atoms.norm1(v)>=1]
prob = amnet.opt.Problem(obj, cons)

result = prob.solve()


    

    




itr | lo          | hi          | gam         | res   | obj         | feas        
===============================================================================
  1 | -1.0486e+06 |  1.0486e+06 |           0 | unsat | None        | False 
  2 |           0 |  1.0486e+06 |  5.2429e+05 | sat   |      10.396 | False 
  3 |           0 |  5.2429e+05 |  2.6214e+05 | sat   |      10.396 | False 
  4 |           0 |  2.6214e+05 |  1.3107e+05 | sat   |      10.396 | False 
  5 |           0 |  1.3107e+05 |       65536 | sat   |      10.396 | False 
  6 |           0 |       65536 |       32768 | sat   |      10.396 | False 
  7 |           0 |       32768 |       16384 | sat   |      10.396 | False 
  8 |           0 |       16384 |        8192 | sat   |      10.396 | False 
  9 |           0 |        8192 |        4096 | sat   |      10.396 | False 
 10 |           0 |        4096 |        2048 | sat   |      10.396 | False 
 11 |           0 |        2048 |        1024 | sat   |      10.396 | False 
 12 |           0 |        1024 |         512 | sat   |      10.396 | False 
 13 |           0 |         512 |         256 | sat   |      10.396 | False 
 14 |           0 |         256 |         128 | sat   |      10.396 | False 
 15 |           0 |         128 |          64 | sat   |      10.396 | False 
 16 |           0 |          64 |          32 | sat   |      10.396 | False 
 17 |           0 |          32 |          16 | sat   |      10.396 | False 
 18 |           0 |          16 |           8 | sat   |           8 | False 
 19 |           0 |           8 |           4 | unsat |         8.0 | False 
 20 |           4 |           8 |           6 | sat   |           6 | False 
 21 |           4 |           6 |           5 | sat   |           5 | False 
 22 |           4 |           5 |         4.5 | sat   |         4.5 | False 
Solution found.
  objval: 4.500000
   point: [ 0.86592261 -0.05959819  1.63639614  0.          0.          0.          0.
  0.          1.87981122  0.          0.          0.          0.          0.
  1.59988557  0.          0.          0.          0.94011293  0.          0.
  0.          0.08019029]

In [8]:

    

# graph AMNET solution
t_min = min(np.hstack((X[0,:], Y[0,:])))
t_max = max(np.hstack((X[0,:], Y[0,:])))
tt = np.linspace(t_min-1, t_max+1, 100)
av = a.eval(result.value)
bv = b.eval(result.value)
p = -av[0]/av[1]*tt + bv/av[1]
p1 = -av[0]*tt/av[1] + (bv+1)/av[1]
p2 = -av[0]*tt/av[1] + (bv-1)/av[1]
plt.plot(X[0,:], X[1,:], 'o', Y[0,:], Y[1,:], 'o')
plt.plot(tt, p, '-r', tt, p1, '--r', tt, p2, '--r')
plt.title('Appriximate Linear Discrimination')
plt.axis('equal')
plt.xlim(-10, 10)
plt.ylim(-10, 10)


    

    
Out[8]:





(-10, 10)