机器学习：支持向量机2

本文来自同步博客。

P.S. 不知道怎么显示数学公式以及更好的排版内容。所以如果觉得文章下面格式乱的话请自行跳转到上述链接。后续我将不再对数学公式进行截图，毕竟行内公式截图的话排版会很乱。看原博客地址会有更好的体验。

上一篇文章介绍了机器学习中支持向量机的基本原理，并且在文章末尾介绍了一种利用Python求解二项规划问题极值的方法。这篇文章我将利用这种方法一步步求解上文中提及的$-vec{alpha}-$、$-vec{w}-$、$-b-$，借此复习和验证支持向量机的知识点。

数据

下面看一组测试数据：

data = {
    '+': [
        [1, 7],
        [2, 8],
        [3, 8],
        [2, 6.5]
    ],
    '-': [
        [5, 1],
        [6, -1],
        [7, 3]
    ]
}

数据data是拥有两种已经分好类的数据，每种类型的数据的元素都是二维向量，可以在笛卡尔坐标系中表示。

依照上一篇文章讲述的原理，我们需要利用这些数据求解一个$-vec{alpha}-$向量。也就是我们需要求解使得二项规划方程值最小时的$-vec{alpha}-$向量：
$$
F(vec{alpha}) = frac{1}{2}vec{alpha}_{T}Hvec{alpha} + vec{c}vec{alpha} + c_0， vec{y}^{T}vec{alpha} = 0， vec{alpha} ge 0
$$

很明显，在支持向量机中，$-c_0 = 0-$。

参数求解

首先利用输入的测试data准备上述方程中出现的变量$-H，c，c_0-$。参考下面代码：

def parseXYC(d):
    X = []
    y = []
    c = []
    for _, v in enumerate(d['+']):
        X.append(np.array(v))
        y.append(1)
        c.append(-1)
    for _, v in enumerate(d['-']):
        X.append(np.array(v))
        y.append(-1)
        c.append(-1)
    return X, y, c, 0

X, y, c, c0 = parseXYC(data)

parseXYC函数把data格式化成$-X, y, c, c_0-$。

然后计算$-H-$矩阵的值。比较简单，一行代码就可以得到：

H = np.array([y[i] * y[j] * np.dot(X[i], X[j]) for i in range(len(X)) for j in range(len(X))]).reshape(len(X), len(X))

求解$-vec{alpha}-$

所有数据都准备好了，接下来就是带入optimize.minimize函数中计算结果。

这里有几个超出本文描述范围的难点需要简单提及一下：

1. optimize.minimize函数求解二项规划使用的SLSQP方法既需要用到二项方程的雅各比导函数，也需要用到约束条件函数的雅各比导函数。不清楚这点导致我在测试过程中一直无法求解到正确的值。
2. 不等式约束条件$-vec{alpha} ge 0-$无法在作为约束条件参数constraints传递给optimize.minimize函数。我猜测是因为我构造的不等式参数是错误的，因此无法让不等式约束条件生效。我尚无法解决这个问题，希望了解该问题的同学能留言赐教。作为一种补救方法，我利用边界约束参数bounds描述$-vec{alpha} ge 0-$这个不等式。
3. 求解出来的$-vec{alpha}-$向量中，部分应该为0的元素无法完整精确到0。我观察测试结果总结出来的精确度应该在1e-16，因此在负16次方这个精确度下的值我都假设它就是0。经过绘图验证，我发现这个假设是合理的。

下面开代码实现：

# 定义二项规划方程fun及其雅各比方程jac
def fun(x, sign=1.):
    return sign * (0.5 * np.dot(x.T, np.dot(H, x))+ np.dot(c, x) + c0)
def jac(x, sign=1.):
    return sign * (np.dot(x.T, H) + c)

# 定义等式约束条件方程feq及其雅各比方程jeq
def feq(x):
    return np.dot(y, x)
def jeq(x):
    return np.array(y)

# 生成相关参数
diff = 1e-16
bounds = [(0, None) for _ in range(len(y))] # x >= 0
constraints = [{ 'type': 'eq', 'fun': feq, 'jac': jeq }]# y*x = 0
options = { 'ftol': diff, 'disp': True }
guess = np.array([0 for _ in range(len(X))])

# 计算结果
res_cons = optimize.minimize(fun, guess, method='SLSQP', jac=jac, bounds=bounds, constraints=constraints, options=options)
alpha = [ 0 if abs(x - 0) <= diff else x for x in res_cons.x ]

# 输出结果与校验y*alpha的值是否为0
print('raw alpha: ', res_cons.x)
print('fmt alpha: ', alpha)
print('check y*alpha: ', 'is 0'if (abs(np.dot(y, res_cons.x) - 0) < diff ) else 'is not 0')

求解$-vec{w}-$和$-b-$

# 计算w = sum(xi*yi*Xi)
w = np.sum([ np.array([0, 0]) if alpha[i] == 0 else (alpha[i] * y[i] * X[i]) for i in range(len(alpha))], axis=0)
print('w: ', w)

# 计算b，对support vector有：yi(w*xi + b) = 1，既有：b = 1/yi - w*xi
B = [( 0 if alpha[i] == 0 else ( 1 / y[i] - np.dot(w, X[i]) ) ) for i in range(len(alpha))]
B = list(filter(lambda x: x != 0, B))
b = 0 if len(B) <= 0 else B[0]
print('b: ', b)

至此支持向量机的参数求解过程完毕。

运行结果如下图所示：

 

绘图

最后把数据绘制成图像。

limit = 11
plt.xlim(-2, limit)
plt.ylim(-2, limit)
# 绘制数据点
[plt.scatter(X[i][0],X[i][1], s=100, color=('r' if y[i] > 0 else 'y')) for i in range(len(X))]
# 绘制分割超平面L： wx + b = 0
plt.plot([i for i in range(limit)], [(-b - w[0]*i)/w[1] for i in range(limit)])
# 绘制上下边： wx + b = 1/-1
plt.plot([i for i in range(limit)], [(1-b - w[0]*i)/w[1] for i in range(limit)])
plt.plot([i for i in range(limit)], [(-1-b - w[0]*i)/w[1] for i in range(limit)])
plt.show()

效果如下图。其中红点为'+'样本，绿点为'-'样本。中间的蓝色线为分类的标准线。边界线，即红色线和绿色线分别穿过各自类别中最靠近分类标准线的点。这些点就是支持向量，只有这些向量所对应的$-vec{alpha}-$分量才为非零值。

 

 

本文源码