-
Notifications
You must be signed in to change notification settings - Fork 21
/
svm.py
125 lines (102 loc) · 3.93 KB
/
svm.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
import numpy as np
import scipy.sparse as spa
import cvxpy
class SVMExample(object):
'''
SVM QP example
'''
def __init__(self, n, seed=1):
'''
Generate problem in QP format and CVXPY format
'''
# Set random seed
np.random.seed(seed)
self.n = int(n) # Number of features
self.m = int(self.n * 100) # Number of data-points
# Generate data
self.N = int(self.m / 2)
self.gamma = 1.0
self.b_svm = np.append(np.ones(self.N), -np.ones(self.N))
A_upp = spa.random(self.N, self.n, density=.15,
data_rvs=np.random.randn)
A_low = spa.random(self.N, self.n, density=.15,
data_rvs=np.random.randn)
self.A_svm = spa.vstack([
A_upp / np.sqrt(self.n) + (A_upp != 0.).astype(float) / self.n,
A_low / np.sqrt(self.n) - (A_low != 0.).astype(float) / self.n
]).tocsc()
self.qp_problem = self._generate_qp_problem()
self.cvxpy_problem, self.cvxpy_variables = \
self._generate_cvxpy_problem()
@staticmethod
def name():
return 'SVM'
def _generate_qp_problem(self):
'''
Generate QP problem
'''
# Construct the problem
# minimize x.T * x + gamma 1.T * t
# subject to t >= diag(b) A x + 1
# t >= 0
P = spa.block_diag((spa.eye(self.n),
spa.csc_matrix((self.m, self.m))), format='csc')
q = np.append(np.zeros(self.n), (self.gamma/2) * np.ones(self.m))
A = spa.vstack([spa.hstack([spa.diags(self.b_svm).dot(self.A_svm),
-spa.eye(self.m)]),
spa.hstack([spa.csc_matrix((self.m, self.n)),
spa.eye(self.m)])
]).tocsc()
l = np.hstack([-np.inf*np.ones(self.m), np.zeros(self.m)])
u = np.hstack([-np.ones(self.m), np.inf*np.ones(self.m)])
# Constraints without bounds
A_nobounds = spa.hstack([spa.diags(self.b_svm).dot(self.A_svm),
-spa.eye(self.m)]).tocsc()
l_nobounds = -np.inf*np.ones(self.m)
u_nobounds = -np.ones(self.m)
bounds_idx = np.arange(self.n, self.n + self.m)
# Separate bounds
lx = np.append(-np.inf*np.ones(self.n), np.zeros(self.m))
ux = np.append(np.inf*np.ones(self.n), np.inf*np.ones(self.m))
problem = {}
problem['P'] = P
problem['q'] = q
problem['A'] = A
problem['l'] = l
problem['u'] = u
problem['m'] = A.shape[0]
problem['n'] = A.shape[1]
problem['A_nobounds'] = A_nobounds
problem['l_nobounds'] = l_nobounds
problem['u_nobounds'] = u_nobounds
problem['bounds_idx'] = bounds_idx
problem['lx'] = lx
problem['ux'] = ux
return problem
def _generate_cvxpy_problem(self):
'''
Generate QP problem
'''
n = self.n
m = self.m
x = cvxpy.Variable(n)
t = cvxpy.Variable(m)
objective = cvxpy.Minimize(.5 * cvxpy.quad_form(x, spa.eye(n))
+ .5 * self.gamma * np.ones(m) * t)
constraints = [t >= spa.diags(self.b_svm).dot(self.A_svm) * x + 1,
t >= 0]
problem = cvxpy.Problem(objective, constraints)
return problem, (x, t)
def revert_cvxpy_solution(self):
'''
Get QP primal and duar variables from cvxpy solution
'''
(x_cvx, t_cvx) = self.cvxpy_variables
constraints = self.cvxpy_problem.constraints
# primal solution
x = np.concatenate((x_cvx.value,
t_cvx.value))
# dual solution
y = np.concatenate((constraints[0].dual_value,
-constraints[1].dual_value))
return x, y