Python手撸机器学习系列（十）：软间隔SVM

软间隔SVM

对于硬间隔SVM，其目标函数为：
$$\begin{gathered} \underset{w,b}{\min }\frac{1}{2}{w}^{T}w \\ s.t.y_i(w^T{x}_i+b)\ge 1 \end{gathered}$$
具体细节可参考https://blog.csdn.net/qq_43601378/article/details/122027005
这个式子是建立在数据集线性可分的情况下得出的

那么问题来了，如果是线性不可分的情况，或者是有噪点的线性可分，该怎么办

在这种情况下，我们假设其是线性可分的，并允许其犯一点点错误，在最优化时给其加上一点点损失，这就是软间隔SVM（soft-margin SVM）

1.1 软间隔SVM模型原理

由上述我们可得软间隔SVM是在硬间隔SVM的基础上加上一点损失，即：
$$\begin{gathered} \underset{w,b}{\min }\frac{1}{2}{w}^{T}w+loss \\ s.t.y_i(w^T{x}_i+b)\ge 1 \end{gathered}$$
现在需要确定这个loss

从感知机的知识我们可以知道，有两种思路：一种是统计所有误分类的点，一种是计算误分类点到超平面的距离。

第一种不利于我们求导，一般选择第二种(此处可以恶补感知机)。

按距离设置loss，可以写作：
l o s s = m a x { 0 , 1 − y i ( w T x i + b ) } loss = max\{0,1-y_i(w^Tx_i+b)\} loss=max{0,1−yi​(wTxi​+b)}
即当$y_i(w^T{x}_i+b)\ge 1$时分类正确，此时损失为0，反之分类有误，损失为$1-y_i(w^T{x}_i+b)$

这个损失也被称为合页损失（hinge loss，亦称为铰链损失）

将所有点的损失相加，就是所有误分类点的损失

当确定loss后，软间隔SVM的优化目标也确定了：
min ⁡ w , b 1 2 w T w + C ∑ i = 1 N m a x { 0 , 1 − y i ( w T x i + b ) } s . t .    y i ( w T x i + b ) ≥ 1 \min_{w,b}\frac{1}{2}w^Tw+C\displaystyle\sum_{i=1}^Nmax\{0,1-y_i(w^Tx_i+b)\} \\ s.t. \ \ y_i(w^Tx_i+b)\geq1 w,bmin​21​wTw+Ci=1∑N​max{0,1−yi​(wTxi​+b)}s.t.  yi​(wTxi​+b)≥1
引入$\xi =1-y_i(w^T{x}_i+b),\xi \ge 0$，则原问题重新表述为：
$$\begin{gathered} \underset{w,b,\xi }{\min }\frac{1}{2}{w}^{T}w+C\sum _{i=1}^N{\xi }_i \\ s.t.y_i(w^T{x}_i+b)\ge 1-{\xi }_i \\ \xi \ge 0 \end{gathered}$$

• 这里的$\xi$被称为松弛变量，即适当放松原来问题的要求，且每个$x_i$都有一个对应的${\xi }_i$
• 这里的$C$被称为惩罚参数，且$C>0$。$C$越大，表示对误分类点的惩罚越大，迫使所有样本均满足约束，即尽可能将其两组点分开，反之允许更多的点跨过分界线。

1.2 软间隔SVM模型求解

根据前面硬间隔的方式，依然将其转化为对偶形式进行求解。

首先引入拉格朗日乘子（注意拉格朗日乘子法中限定不等式必须都是小于等于0的形式，所以有一点转化）：
$$L(w,b,\xi ,\alpha ,\beta )=\frac{1}{2}{w}^{T}w+C\sum _{i=1}^N{\xi }_i+\sum _{i=1}^N{\alpha }_i(1-{\xi }_i-y_i(w^T{x}_i+b))+\sum _{i=1}^N\beta (-\xi )$$
原问题的对偶问题为（假设已满足KKT条件，交换了min和max的位置）：
$$\begin{gathered} \underset{\alpha ,\beta }{\max }\underset{w,b,\xi }{\min }L(w,b,\xi ,\alpha ,\beta ) \\ s.t.{\alpha }_i\ge 0,i=1,2,...,N \\ {\beta }_i\ge 0,i=1,2,...,N \end{gathered}$$
求解$\underset{w,b,\xi }{\min }L(w,b,\xi ,\alpha ,\beta )$

分别对$w,b,\xi$求导并令其为0：
$$\{\begin{array}{l}\frac{\partial L}{\partial w}=w-\sum _{i=1}^N{\alpha }_i{y}_i{x}_i=0\\ \frac{\partial L}{\partial b}=-\sum _{i=1}^N{\alpha }_i{y}_i=0\\ \frac{\partial L}{\partial \xi }=C-{\alpha }_i-{\beta }_i=0\end{array}$$
可得：
$$\{\begin{array}{l}w=\sum _{i=1}^N{\alpha }_i{y}_i{x}_i\\ \sum _{i=1}^N{\alpha }_i{y}_i=0\\ {\beta }_i=C-{\alpha }_i\end{array}$$
带入$L(w,b,\xi ,\alpha ,\beta )$，得：
$$L(w,b,\xi ,\alpha ,\beta )=-\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j{x}_i{x}_j+\sum _{i=1}^N{\alpha }_i$$
原问题简化为：
$$\begin{gathered} \underset{\alpha }{\max }-\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j{x}_i{x}_j+\sum _{i=1}^N{\alpha }_i \\ s.t.\sum _{i=1}^N{\alpha }_i{y}_i=0 \\ 0\le {\alpha }_i\le C \end{gathered}$$
转化一下，可得到：

$$\begin{gathered} \underset{\alpha }{\min }\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j{x}_i{x}_j-\sum _{i=1}^N{\alpha }_i \\ s.t.\sum _{i=1}^N{\alpha }_i{y}_i=0 \\ 0\le {\alpha }_i\le C \end{gathered}$$
对于全部的样本，存在$(x_k,y_k)$，使得$1-{\xi }_k-y_k(w^T{x}_k+b)=0$

此时可得$b=y_k-y_k{\xi }_k-w^T{x}_k$

如果${\xi }_k>0$，两边同时乘上$y_k$，则可得${\xi }_k=1-y_k(w^T{x}_k+b)$，无解

如果${\xi }_k=0$，则可求出b，此时$0<{\alpha }_k<C$，

求导时我们以及得到
$$w^{\ast }=\sum _{i=1}^N{\alpha }_i{y}_i{x}_i$$
求得
$$b^{\ast }=y_k-w^{\ast }{}^T{x}_k=y_k-\sum _{i=1}^N{y}_i{\alpha }_i(x_i,x_k)$$
即可求出超平面$w^{\ast }{}^{T}x+b$

后面就是和硬间隔SVM一样的使用SMO算法求得$\alpha$，再求出$w$和$b$，SMO算法的实现推导可以参考我的上一篇博客，点击跳转

与硬间隔SVM的对偶形式相比，仅多了$0\le {\alpha }_i\le C$这一个限制条件，而这个参数一般由人为设定，即超参数

1.3 代码实现

import numpy as np
import random
import matplotlib.pyplot as plt

def simple_smo(dataset, labels, C, max_iter):
    ''' 简化版SMO算法实现，未使用启发式方法对alpha对进行选择.
    :param dataset: 所有特征数据向量
    :param labels: 所有的数据标签
    :param C: 软间隔常数, 0 <= alpha_i <= C
    :param max_iter: 外层循环最大迭代次数
    '''
    dataset = np.array(dataset)
    m, n = dataset.shape #样本数量，特征数量
    labels = np.array(labels)

    # 初始化参数λ，b为0
    lambds = np.zeros(m) #每个样本都有一个λ乘子
    b = 0
    it = 0


    while it < max_iter:
        pair_changed = 0 #选取的一对值相较于之前是否有变化，没有的话可以停止迭代
        for i in range(m):
            λ_i, x_i, y_i = lambds[i], dataset[i], labels[i] #选取一组λ
            fx_i = SVM_predict(x_i,lambds,dataset,labels,b)
            E_i = fx_i - y_i
            j = select_j(i, m) #选取另一个λ
            λ_j, x_j, y_j = lambds[j], dataset[j], labels[j]
            fx_j = SVM_predict(x_j,lambds,dataset,labels,b)
            E_j = fx_j - y_j
            K_ii, K_jj, K_ij = np.dot(x_i, x_i), np.dot(x_j, x_j), np.dot(x_i, x_j)
            eta = K_ii + K_jj - 2*K_ij #
            if eta <= 0:
                print('WARNING  eta <= 0')
                continue
            # 获取更新的alpha对
            λ_i_old, λ_j_old = λ_i, λ_j #未更新前的参数
            λ_j_new = λ_j_old + y_j*(E_i - E_j)/eta
            # 对alpha进行修剪
            if y_i != y_j:
                L = max(0, λ_j_old - λ_i_old)
                H = min(C, C + λ_j_old - λ_i_old)
            else:
                L = max(0, λ_i_old + λ_j_old - C)
                H = min(C, λ_j_old + λ_i_old)
            λ_j_new = clip(λ_j_new, L, H) #根据上下界修剪
            λ_i_new = λ_i_old + y_i*y_j*(λ_j_old - λ_j_new) #根据公式反推另一个参数
            if abs(λ_j_new - λ_j_old) < 0.00001: #这个参数已经优化到最佳，换下一个
                #print('WARNING   alpha_j not moving enough')
                continue

            #更新b
            lambds[i], lambds[j] = λ_i_new, λ_j_new
            b_i = -E_i - y_i*K_ii*(λ_i_new - λ_i_old) - y_j*K_ij*(λ_j_new - λ_j_old) + b
            b_j = -E_j - y_i*K_ij*(λ_i_new - λ_i_old) - y_j*K_jj*(λ_j_new - λ_j_old) + b
            if 0 < λ_i_new < C:
                b = b_i
            elif 0 < λ_j_new < C:
                b = b_j
            else:
                b = (b_i + b_j)/2
            pair_changed += 1
            print('INFO   iteration:{}  i:{}  pair_changed:{}'.format(it, i, pair_changed))
        if pair_changed == 0: #参数优化完成，下一轮
            it += 1
        else: #参数没有优化完成，继续迭代
            it = 0
        print('iteration number: {}'.format(it))
    return lambds, b

def SVM_predict(x,lambds,data,label,b):
    "SVM分类器函数 y = w^Tx + b,即文中的f(x)"
    res = 0
    for i in range(data.shape[0]):
        res += lambds[i]*label[i]*(data[i].dot(x.T))
    return res + b

def get_w(lambdas, dataset, labels):
    #通过λ求w
    w = 0
    for i in range(len(dataset)):
        w += lambdas[i]*y[i]*dataset[i]
    return w

def clip(alpha, L, H):
    #修建λ的值到L和H之间.

    if alpha < L:
        return L
    elif alpha > H:
        return H
    else:
        return alpha

def select_j(i, m):
    #在m中随机选择除了i之外剩余的
    l = list(range(m))
    seq = l[: i] + l[i+1:]
    return random.choice(seq)

def get_point():
    x_true =  [[1,1.5],[1,3],[4,5],[2,4]]
    x_false = [[1,0.5],[4,2],[5,1],[4,1]]
    x_true =  [[1,2],[1,3],[4,5],[2,4],[2,1.5]]
    x_false = [[1,0.5],[4,2],[5,1],[4,1],[3,2]]
    x_all = np.array(x_true+x_false)
    y = [1]*len(x_true) + [-1]*len(x_false)
    return x_all,y,x_true,x_false

def plot(x_true,x_false,w,b):
    plot_x = np.arange(0,7,0.1)
    plot_y = -(w[0]*plot_x+b)/w[1]
    plt.scatter([x[0] for x in x_true],[x[1] for x in x_true] , c='r' , label='+1')
    plt.scatter([x[0] for x in x_false],[x[1] for x in x_false] , c='b',label='-1')
    plt.plot(plot_x,plot_y,c = 'green')
    plt.xlim(0,7)
    plt.ylim(0,7)
    plt.legend()
    plt.plot()
    plt.show()

if __name__ == '__main__':
    x,y,x_true,x_false = get_point()
    MAX_EPOCHES = 20 #最大迭代次数
    C = 50 #λ边界
    #调节C的值会有不同的改变
    lambdas, b = simple_smo(x, y, C, MAX_EPOCHES)
    w = get_w(lambdas, x, y)
    print('-'*40+'result'+'-'*40)
    print('lambdas:{}\nw:{}\nb:{}'.format(lambdas,w,b))
    plot(x_true,x_false,w,b)

结果：

当C分别取50，10，1，0.1时：

1.4 参考文献

1. https://zhuanlan.zhihu.com/p/309944707
2. https://zhuanlan.zhihu.com/p/227896319
3. https://zhuanlan.zhihu.com/p/29212107