Skip to content

Commit 54961a6

Browse files
author
louwill
authored
Add files via upload
1 parent 1797e16 commit 54961a6

File tree

3 files changed

+559
-0
lines changed

3 files changed

+559
-0
lines changed

soft_margin_svm/Qp-cvxopt.pdf

139 KB
Binary file not shown.
Lines changed: 291 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,291 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": 1,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"import numpy\n",
10+
"from cvxopt import matrix\n",
11+
"from cvxopt import solvers"
12+
]
13+
},
14+
{
15+
"cell_type": "code",
16+
"execution_count": 2,
17+
"metadata": {},
18+
"outputs": [],
19+
"source": [
20+
"# 定义二次规划参数\n",
21+
"P = matrix([[1.0,0.0],[0.0,0.0]])\n",
22+
"q = matrix([3.0,4.0])\n",
23+
"G = matrix([[-1.0,0.0,-1.0,2.0,3.0],[0.0,-1.0,-3.0,5.0,4.0]])\n",
24+
"h = matrix([0.0,0.0,-15.0,100.0,80.0])"
25+
]
26+
},
27+
{
28+
"cell_type": "code",
29+
"execution_count": 3,
30+
"metadata": {},
31+
"outputs": [
32+
{
33+
"name": "stdout",
34+
"output_type": "stream",
35+
"text": [
36+
" pcost dcost gap pres dres\n",
37+
" 0: 1.0780e+02 -7.6366e+02 9e+02 4e-17 4e+01\n",
38+
" 1: 9.3245e+01 9.7637e+00 8e+01 8e-17 3e+00\n",
39+
" 2: 6.7311e+01 3.2553e+01 3e+01 8e-17 1e+00\n",
40+
" 3: 2.6071e+01 1.5068e+01 1e+01 7e-17 7e-01\n",
41+
" 4: 3.7092e+01 2.3152e+01 1e+01 1e-16 4e-01\n",
42+
" 5: 2.5352e+01 1.8652e+01 7e+00 9e-17 4e-16\n",
43+
" 6: 2.0062e+01 1.9974e+01 9e-02 7e-17 2e-16\n",
44+
" 7: 2.0001e+01 2.0000e+01 9e-04 8e-17 2e-16\n",
45+
" 8: 2.0000e+01 2.0000e+01 9e-06 1e-16 2e-16\n",
46+
"Optimal solution found.\n"
47+
]
48+
}
49+
],
50+
"source": [
51+
"# 构建求解\n",
52+
"sol = solvers.qp(P,q,G,h)"
53+
]
54+
},
55+
{
56+
"cell_type": "code",
57+
"execution_count": 8,
58+
"metadata": {},
59+
"outputs": [
60+
{
61+
"name": "stdout",
62+
"output_type": "stream",
63+
"text": [
64+
"[ 7.13e-07]\n",
65+
"[ 5.00e+00]\n",
66+
" 20.00000617311241\n"
67+
]
68+
}
69+
],
70+
"source": [
71+
"# 获取最优值\n",
72+
"print(sol['x'],sol['primal objective'])"
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": 9,
78+
"metadata": {},
79+
"outputs": [],
80+
"source": [
81+
"import numpy as np\n",
82+
"from numpy import linalg\n",
83+
"import cvxopt\n",
84+
"import cvxopt.solvers\n",
85+
"import pylab as pl"
86+
]
87+
},
88+
{
89+
"cell_type": "code",
90+
"execution_count": 10,
91+
"metadata": {},
92+
"outputs": [],
93+
"source": [
94+
"# 定义一个线性核\n",
95+
"def linear_kernel(x1, x2):\n",
96+
" return np.dot(x1, x2)"
97+
]
98+
},
99+
{
100+
"cell_type": "code",
101+
"execution_count": 11,
102+
"metadata": {},
103+
"outputs": [],
104+
"source": [
105+
"def gen_non_lin_separable_data():\n",
106+
" mean1 = [-1, 2]\n",
107+
" mean2 = [1, -1]\n",
108+
" mean3 = [4, -4]\n",
109+
" mean4 = [-4, 4]\n",
110+
" cov = [[1.0, 0.8], [0.8, 1.0]]\n",
111+
" X1 = np.random.multivariate_normal(mean1, cov, 50)\n",
112+
" X1 = np.vstack((X1, np.random.multivariate_normal(mean3, cov, 50)))\n",
113+
" y1 = np.ones(len(X1))\n",
114+
" X2 = np.random.multivariate_normal(mean2, cov, 50)\n",
115+
" X2 = np.vstack((X2, np.random.multivariate_normal(mean4, cov, 50)))\n",
116+
" y2 = np.ones(len(X2)) * -1\n",
117+
" return X1, y1, X2, y2\n",
118+
"\n",
119+
"X1, y1, X2, y2 = gen_non_lin_separable_data()"
120+
]
121+
},
122+
{
123+
"cell_type": "code",
124+
"execution_count": 12,
125+
"metadata": {},
126+
"outputs": [
127+
{
128+
"name": "stdout",
129+
"output_type": "stream",
130+
"text": [
131+
"(180, 2) (180,) (20, 2) (20,)\n"
132+
]
133+
}
134+
],
135+
"source": [
136+
"def split_train(X1, y1, X2, y2):\n",
137+
" X1_train = X1[:90]\n",
138+
" y1_train = y1[:90]\n",
139+
" X2_train = X2[:90]\n",
140+
" y2_train = y2[:90]\n",
141+
" X_train = np.vstack((X1_train, X2_train))\n",
142+
" y_train = np.hstack((y1_train, y2_train))\n",
143+
" return X_train, y_train\n",
144+
"\n",
145+
"\n",
146+
"def split_test(X1, y1, X2, y2):\n",
147+
" X1_test = X1[90:]\n",
148+
" y1_test = y1[90:]\n",
149+
" X2_test = X2[90:]\n",
150+
" y2_test = y2[90:]\n",
151+
" X_test = np.vstack((X1_test, X2_test))\n",
152+
" y_test = np.hstack((y1_test, y2_test))\n",
153+
" return X_test, y_test\n",
154+
"\n",
155+
"X_train, y_train = split_train(X1, y1, X2, y2)\n",
156+
"X_test, y_test = split_test(X1, y1, X2, y2)\n",
157+
"print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)"
158+
]
159+
},
160+
{
161+
"cell_type": "code",
162+
"execution_count": 27,
163+
"metadata": {},
164+
"outputs": [],
165+
"source": [
166+
"def fit(X, y, C):\n",
167+
" n_samples, n_features = X.shape\n",
168+
"\n",
169+
" # Gram matrix\n",
170+
" K = np.zeros((n_samples, n_samples))\n",
171+
" for i in range(n_samples):\n",
172+
" for j in range(n_samples):\n",
173+
" K[i, j] = linear_kernel(X[i], X[j])\n",
174+
"\n",
175+
" P = cvxopt.matrix(np.outer(y, y) * K)\n",
176+
" q = cvxopt.matrix(np.ones(n_samples) * -1)\n",
177+
" A = cvxopt.matrix(y, (1, n_samples))\n",
178+
" b = cvxopt.matrix(0.0)\n",
179+
"\n",
180+
" if C is None:\n",
181+
" G = cvxopt.matrix(np.diag(np.ones(n_samples) * -1))\n",
182+
" h = cvxopt.matrix(np.zeros(n_samples))\n",
183+
" else:\n",
184+
" tmp1 = np.diag(np.ones(n_samples) * -1)\n",
185+
" tmp2 = np.identity(n_samples)\n",
186+
" G = cvxopt.matrix(np.vstack((tmp1, tmp2)))\n",
187+
" tmp1 = np.zeros(n_samples)\n",
188+
" tmp2 = np.ones(n_samples) * C\n",
189+
" h = cvxopt.matrix(np.hstack((tmp1, tmp2)))\n",
190+
"\n",
191+
" # solve QP problem\n",
192+
" solution = cvxopt.solvers.qp(P, q, G, h, A, b)\n",
193+
"\n",
194+
" # Lagrange multipliers\n",
195+
" a = np.ravel(solution['x'])\n",
196+
" # Support vectors have non zero lagrange multipliers\n",
197+
" sv = a > 1e-5\n",
198+
" ind = np.arange(len(a))[sv]\n",
199+
" a = a[sv]\n",
200+
" sv_x = X[sv]\n",
201+
" sv_y = y[sv]\n",
202+
" print(\"%d support vectors out of %d points\" % (len(a), n_samples))\n",
203+
"\n",
204+
" # Intercept\n",
205+
" b = 0\n",
206+
" for n in range(len(a)):\n",
207+
" b += sv_y[n]\n",
208+
" b -= np.sum(a * sv_y * K[ind[n], sv])\n",
209+
" b /= len(a)\n",
210+
"\n",
211+
" # Weight vector\n",
212+
" w = np.zeros(n_features)\n",
213+
" for n in range(len(a)):\n",
214+
" w += a[n] * sv_y[n] * sv[n]\n",
215+
" else:\n",
216+
" w = None"
217+
]
218+
},
219+
{
220+
"cell_type": "code",
221+
"execution_count": 28,
222+
"metadata": {},
223+
"outputs": [
224+
{
225+
"name": "stdout",
226+
"output_type": "stream",
227+
"text": [
228+
" pcost dcost gap pres dres\n",
229+
" 0: -5.5569e+04 -7.4952e+07 8e+07 1e-02 5e-11\n",
230+
" 1: -7.3687e+04 -9.6455e+05 9e+05 1e-04 5e-11\n",
231+
" 2: -8.2531e+04 -1.6849e+05 9e+04 1e-05 6e-11\n",
232+
" 3: -1.2558e+05 -1.4468e+05 2e+04 1e-06 8e-11\n",
233+
" 4: -1.2939e+05 -1.3440e+05 5e+03 3e-07 8e-11\n",
234+
" 5: -1.3084e+05 -1.3145e+05 6e+02 2e-08 9e-11\n",
235+
" 6: -1.3103e+05 -1.3115e+05 1e+02 4e-09 9e-11\n",
236+
" 7: -1.3107e+05 -1.3109e+05 2e+01 6e-10 1e-10\n",
237+
" 8: -1.3108e+05 -1.3108e+05 2e+00 4e-11 1e-10\n",
238+
" 9: -1.3108e+05 -1.3108e+05 2e-01 3e-11 1e-10\n",
239+
"10: -1.3108e+05 -1.3108e+05 2e-03 4e-11 1e-10\n",
240+
"Optimal solution found.\n",
241+
"158 support vectors out of 180 points\n"
242+
]
243+
}
244+
],
245+
"source": [
246+
"w, b = fit(X_train, y_train, C=1000.1)"
247+
]
248+
},
249+
{
250+
"cell_type": "code",
251+
"execution_count": null,
252+
"metadata": {},
253+
"outputs": [],
254+
"source": []
255+
}
256+
],
257+
"metadata": {
258+
"kernelspec": {
259+
"display_name": "Python 3",
260+
"language": "python",
261+
"name": "python3"
262+
},
263+
"language_info": {
264+
"codemirror_mode": {
265+
"name": "ipython",
266+
"version": 3
267+
},
268+
"file_extension": ".py",
269+
"mimetype": "text/x-python",
270+
"name": "python",
271+
"nbconvert_exporter": "python",
272+
"pygments_lexer": "ipython3",
273+
"version": "3.6.5"
274+
},
275+
"toc": {
276+
"base_numbering": 1,
277+
"nav_menu": {},
278+
"number_sections": true,
279+
"sideBar": true,
280+
"skip_h1_title": false,
281+
"title_cell": "Table of Contents",
282+
"title_sidebar": "Contents",
283+
"toc_cell": false,
284+
"toc_position": {},
285+
"toc_section_display": true,
286+
"toc_window_display": false
287+
}
288+
},
289+
"nbformat": 4,
290+
"nbformat_minor": 2
291+
}

0 commit comments

Comments
 (0)