【367】通过 python 实现 SVM 硬边界 算法

参考: 支持向量机整理

SVM 硬边界的结果如下:

$$min \quad \frac{1}{2} \sum_{i=1}^m\sum_{j=1}^m \alpha_i\alpha_jy_iy_j \vec x_i \vec x_j - \sum_{i=1}^m\alpha_i \\ s.t. \quad \alpha_i\ge0 \quad i=1...m \\ \quad \sum_{i=1}^m \alpha_i y_i=0$$

---

一. 数据准备

测试数据如下所示, 前两个为 -1, 后面三个为 1, 如下图可以看到分割线即为:

$$y = x + 1$$

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

import numpy as np import matplotlib.pyplot as plt   X = np.array([[1,3],              [0,2],              [0,0],              [2,0],              [2,2]])   x = np.linspace(-2, 3, 100)   y = np.array([-1,-1,1,1,1])   plt.figure() plt.scatter(X[:2,0],X[:2,1]) plt.scatter(X[2:,0],X[2:,1]) plt.plot(x, x+1) plt.show()

二. 获取 QP 的参数并计算α

将下面的结果带入到二次规划问题中分别求得 P/p/G/h/A/b 的值.

这个过程不是很容易, 虽然数据量这么少, 我反反复复弄了好几遍最终才做对.

$$min \quad \frac{1}{2} \sum_{i=1}^m\sum_{j=1}^m \alpha_i\alpha_jy_iy_j \vec x_i \vec x_j - \sum_{i=1}^m\alpha_i \\ s.t. \quad \alpha_i\ge0 \quad i=1...m \\ \quad \sum_{i=1}^m \alpha_i y_i=0$$

按照下面的形式进行获取参数.

$$min \quad \frac{1}{2}x^TPx + q^Tx\\ s.t. \quad Gx \le h\\ \quad\quad Ax = b$$

?

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

# 需要将 X 中的数据彼此相乘, 得到一个 5*5 的矩阵 # 同时需要注意 y 的符号会影响 P = matrix([[10.0,6.0,0.0,-2.0,-8.0],             [6.0,4.0,0.0,0.0,-4.0],             [0.0,0.0,0.0,0.0,0.0],             [-2.0,0.0,0.0,4.0,4.0],             [-8.0,-4.0,0.0,4.0,8.0]])   # 为了得到一个常数, q 为 5*1 的矩阵, 转置后正好可以用 q = matrix(-1.0, (5,1))   # 首先将 ≥ 调整为 ≤, 然后按照向量的形式表示 # 结果 h 为 5*1 的矩阵 # 因此 G 为 5*5 的矩阵(α 是 5*1 矩阵) G = matrix([[-1.0,0.0,0.0,0.0,0.0],             [0.0,-1.0,0.0,0.0,0.0],             [0.0,0.0,-1.0,0.0,0.0],             [0.0,0.0,0.0,-1.0,0.0],             [0.0,0.0,0.0,0.0,-1.0]])   h = matrix(0.0, (5,1))   # 结果为常数的形式, 因此 A 是一个 1*5 矩阵 A = matrix([1.0,1.0,-1.0,-1.0,-1.0]).T   b = matrix(0.0, (1,1))

将上面的内容带入到二次规划的函数中进行求解.　　

?

1 2 3

sol = solvers.qp(P,q,G,h,A,b) alpha = sol['x'] print(alpha)

     pcost       dcost       gap    pres   dres
 0: -1.4151e+00 -3.0463e+00  1e+01  3e+00  2e+00
 1: -3.4780e-01 -2.4147e+00  2e+00  7e-16  8e-16
 2: -9.2882e-01 -1.0856e+00  2e-01  3e-16  5e-16
 3: -9.9882e-01 -1.0010e+00  2e-03  2e-16  3e-16
 4: -9.9999e-01 -1.0000e+00  2e-05  1e-16  3e-16
 5: -1.0000e+00 -1.0000e+00  2e-07  2e-16  2e-16
Optimal solution found.
[ 4.31e-01]
[ 5.69e-01]
[ 2.84e-01]
[ 5.88e-08]
[ 7.16e-01]

三. 根据α来计算w

目前已经求出了所有的$\alpha$, 根据下面的公式将所有的样本点数据带入求得$\vec w$. 根据$\alpha$的结果可以判断哪些是支持向量, 包括 index = 0, 1, 2, 4 都满足.

$$\vec w=\sum_{i=1}^m \alpha_i y_i \vec x_i$$

?

1 2 3 4 5 6 7 8 9 10

X0 = X[:,0].flatten() X1 = X[:,1].flatten()   w1 = (w*y*X0).sum() w2 = (w*y*X1).sum()   W = np.array([w1,w2])   print("w1=", w1) print("w2=", w2)

w1= 1.0000000446896518
w2= -1.000000054585139

四. 根据w来求b

$\vec w$已经求出了, 这时候只要带入任何一个支持向量里面即可, 公式如下:

$$y_i(\vec w ^T\vec x_i+b) = 1$$

化简后得到:

$$b = y_i - \vec w^T\vec x_i$$

由上面计算可知, 第一个点在支持向量上面, 因此可以计算获得b值.

$$b = y_1 - \vec w^T\vec x_1$$

?

1 2 3 4 5 6

W = np.mat(W) xx = np.mat(X[1,:].T) xx = xx.T   b = int(-1 - W*xx) print("b=", b)

b= 1

所以最终的结果就是:

$$x_1 - x_2 + 1 = 0$$

将$x_1$换成$x$, 将$x_2$换成$y$, 则得到:

$$y = x + 1$$