参考: 支持向量机整理
SVM 硬边界的结果如下:
$$
min quad frac{1}{2} sum_{i=1}^msum_{j=1}^m alpha_ialpha_jy_iy_j vec x_i vec x_j - sum_{i=1}^malpha_i
\
s.t. quad alpha_ige0 quad i=1...m
\
quad sum_{i=1}^m alpha_i y_i=0
$$
一. 数据准备
测试数据如下所示, 前两个为 -1, 后面三个为 1, 如下图可以看到分割线即为:
$$
y = x + 1
$$
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}^msum_{j=1}^m alpha_ialpha_jy_iy_j vec x_i vec x_j - sum_{i=1}^malpha_i
\
s.t. quad alpha_ige0 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\
quadquad Ax = b
$$
# 需要将 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))
将上面的内容带入到二次规划的函数中进行求解.
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
$$
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 ^Tvec x_i+b) = 1
$$
化简后得到:
$$
b = y_i - vec w^Tvec x_i
$$
由上面计算可知, 第一个点在支持向量上面, 因此可以计算获得b值.
$$
b = y_1 - vec w^Tvec x_1
$$
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
$$