[数据结构与算法]支持向量机SVM最详解推导及代码实现smo算法

带约束条件优化问题

再讲SVM算法之前，先了解下带约束条件的优化问题，为后面svm优化奠定基础。

求一个函数的最(极)值

分三种情况：

1. 无约束条件: 就是我们最常见的情况，直接对变量求导并令其为零，求极值即可；

2. 含有等式约束:引入拉格朗日乘子，拉格朗日乘子的思想是将约束条件函数与原函数联系到一起,构造拉格朗日函数L，求导等于0即为最优解；$L(x,\lambda )=f(x)?{\sum }_{i=1}^n{\lambda }_i{h}_i(x)$

3. 不等式约束 ：在等式约束条件的基础上，再引入一个松弛变量β，也叫KKT乘子。
   $L(x,\lambda ,\beta )=f(x)?{\sum }_{i=1}^n{\lambda }_i{h}_i(x)?{\sum }_{i=1}^n{\beta }_i{g}_i(x)$

SVM支持向量机原理

SVM最开始提出是用来解决二分类问题的，但后来可以应用到多分类，但基础还是多分类。比如，一共有5类，先把1作为1类，把2345作为一类，不属于1分类里面的数据，就在2345里面找，再把2345再分为两类，2和345两类，一直向这样迭代…

最大间隔分类器

SVM的核心思想：我们要找到一个超平面(w,b),使得离超平面最近的样本点距离最大。

如上图，我们有这样两类数据，我们的任务就是要找到这样一条线，把这两类数据分开。
我们以线性分类器为例，我们可以用1表示属于类别C1，-1类别C2，找到一个超平面，使两类数据分开，这个超平面的函数方程是：
$$w^{T}x+b=0$$
可以想象得到，将+1和-1分开的这条线，其实是有很多条的，哪一条最合适呢？比如图上有两个超平面，当然是黑色的合适些，因为蓝色超平面里样本点很近，就很容易分类错误，所以他的泛化误差没有另外一个好。

因此，SVM的核心是：我们要找到一个超平面(w,b),使得离超平面最近的样本点距离最大。

我们的最大间隔分类就是使这个点到直线的距离γ(又称为间隔margin)最大，并且在正确分类(因为这个是一个分类器)情况下，$w^{T}x+b<0，y=?1，w^{T}x+b>0，y=+1$，整合这两个情况，记为$y(w^{T}x+b)>0$。因此，我们的优化目标变为：
$$\begin{gathered} max\gamma \\ s.t{y}_i(w^T{x}_i+b)>0 \end{gathered}$$
又因为点到直线的距离公式：

(第二个等式是因为去绝对值，前面讲过约束条件大于0了)

有n个样本点，就有n个距离，我们要找离超平面最近的，因此，$\gamma =min{\gamma }^{(i)}$

(接下来都是对优化函数的化简，因为其他人讲SVM都是一笔带过，直接把它化成了最后的结果，省去了很多步骤，这里我详细讲一下是如何化简的)

把两个式子带入优化目标中，整理一下得：
$\underset{w,b}{\max }\underset{x^{(i)}}{m}in\frac{1}{||w||y^{(i)}}(w^T{x}^{(i)}+b)$
$s.t.y^i(w^T+b)>0$

上述优化目标也再次说明，SVM的核心思想，求离超平面最近的点的距离最大！
公式中||w||可以提到前面(因为min中和w无关)，因此优化目标变为：

还可以化简…看下面的约束条件𝑦^𝑖 (𝑤^𝑇 𝑥^𝑖+𝑏)>0，所以一定存在?𝛾>0，使得𝑚𝑖𝑛 𝑦^𝑖 (𝑤^𝑇 𝑥^𝑖+𝑏)=𝛾 (因为它的最小值等于𝛾且大于零，所以𝑦^𝑖 (𝑤^𝑇 𝑥^𝑖+𝑏)所有的值都大于0) 。𝑚𝑖𝑛 𝑦^𝑖 (𝑤^𝑇 𝑥^𝑖+𝑏)=𝛾 这个形式带入上面的函数中，整理得：

$\underset{w,b}{m}ax\frac{1}{||w||}\gamma$
$s.t.y_i(w^T{x}_i+b)>0$

还是看$min{y}^i(w^T{x}^i+b)=\gamma$，可以推出$y^i(w^T{x}^i+b)\ge \gamma$(最小值等于𝛾，所以它的所有值都大于等于𝛾)，因此我们的约束条件变为：
$$\begin{gathered} \underset{w,b}{m}ax\frac{1}{||w||}\gamma \\ s.t.y_i(w^T{x}_i+b)\ge \gamma \end{gathered}$$

(接下来要讲的是比较重要的，为什么𝛾=1，很多书都没写)
因为我们的超平面方程为$w^{T}x+b=0$ ，同比例放大缩小w和b，我们的超平面不变，为什么？因为等式两边都乘以2之后：$2w^{T}x+2b=2?0$， 这个等式没变，所以超平面也没变。

我是这样理解的，我觉得这种理解方法更容易理解，举个例子：比如在三维空间内，我们同比例缩放w和b，就相当于把我们的超平面的面积增大或者缩小，但其实它的所在位置没有变化，因此这个超平面没有变化！

所以我们将𝛾除以一个适当的值，我们的𝛾是可以等于1的，因此，最后化简为：
$\underset{w,b}{m}ax\frac{1}{||w||}$
$s.t.y_i(w^T{x}_i+b)\ge 1$

再化简：
$\underset{w,b}{m}in\frac{1}{2}||w|{|}^2$
$s.t.y_i(w^T{x}_i+b)\ge 1$

（这的二分之一和平方是因为后面要优化这个函数，对他求导，使计算更加简便，并且对我们找这个超平面无影响）

所以问题就变成了一个带约束条件的优化问题，前面提到过，引入拉格朗日乘子，构造拉格朗日函数L，求它的最大值：

因此，优化目标函数变为
$$\underset{w,b}{m}in\underset{\lambda }{m}axL(w,b,\lambda )$$

为什么可以这样化简呢？可以这样理解，以下不是证明，只是逻辑上的理解：

再经过强对偶化简为：
$$\underset{\lambda }{m}ax\underset{w,b}{m}inL(w,b,\lambda )$$

强对偶是可以进行证明的，但是这个特别复杂，有兴趣的可以单独了解。为什么要这样化简，因为优化拉格朗日函数L，是一个凸优化问题，求它的最大值不方便，求最小值方便，所以要用强对偶。因此，$\underset{w,b}{m}inL(w,b,\lambda )$对w和b求偏导令其为0，$w^{?}={\sum }_{i=1}^n{\lambda }_i{y}_i{x}_i，{\sum }_{i=1}^n{\lambda }_i{y}_i=0$

把w带入拉格朗日函数L中，化简整理得：

和前面一样，凸优化问题，求最大值比较困难，我们取相反数求最小值，最终将问题转化为：

为什么要用对偶问题

1. 首先是我们有不等式约束方程，这就需要我们写成min max的形式来得到最优解。而这种写成这种形式对w，b不能求导，所以我们需要转换成max min的形式。而为了满足这种对偶变换成立，就需要满足KKT条件（KKT条件是原问题与对偶问题等价的必要条件，当原问题是凸优化问题时，变为充要条件）；
2. 对偶问题将原始问题中的不等式约束条件转为了对偶问题中的等式约束条件，方便引入拉格朗日函数；
3. 方便核函数的引入kernel。

核函数kernel function

核函数$K(x,z)=\text{?}(x{)}^T\text{?}(z)=<\text{?}(x),\text{?}(z)>=(x^{T}z{)}^2$, 𝜙是从n维到m维映射。举个例子：

$$x=\underset{3维}{?\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}?},z=\underset{3维}{?\begin{array}{c}{z}_1\\ z_2\\ z_3\end{array}?},\text{?}(x)=\underset{9维}{?\begin{array}{c}{x}_1{x}_1\\ x_1{x}_2\\ x_1{x}_3\\ x_2{x}_1\\ x_2{x}_2\\ x_2{x}_3\\ x_3{x}_1\\ x_3{x}_2\\ x_3{x}_3\end{array}?},\text{?}(z)=\underset{9维}{?\begin{array}{c}{z}_1{z}_1\\ z_1{z}_2\\ z_1{z}_3\\ z_2{z}_1\\ z_2{z}_2\\ z_2{z}_3\\ z_3{z}_1\\ z_3{z}_2\\ z_3{z}_3\end{array}?},K(x,z)=(x^T,z{)}^2$$
Kernel函数将原本9维的计算变成了求x和z内积的平方，变成3维的问题。因此kernel函数就是帮我们省去了在高位空间里进行繁琐计算的简便运算法.

为什么要用核函数？

举个例子：假设把a和b之间红色部分里的所有点定为正类，两边的黑色部分里的点定为负类。试问能找到一个线性函数把两类正确分开吗？
答案是不能的，因为二维空间里的线性函数就是指直线，显然找不到符合条件的直线。但我们可以找到一条二次曲线将它们分开，如下图。

因此，这个二次曲线的函数表达式$g(x)=c_0+c_{1}x+c_2{x}^2=<a,y>$,新建两个向量：

$$a=?\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}?=?\begin{array}{c}{c}_0\\ c_1\\ c_2\end{array}?,y=?\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}?=?\begin{array}{c}1\\ 2\\ x^2\end{array}?$$

所以，原来在二维空间中一个线性不可分的问题，映射到高维空间后，变成了线性可分的。
看下图，左边的数据本来是线性不可分的，经过了$\text{?}$函数后，有一个超平面可以将原本线性不可分的数据变成了线性可分的。

软间隔支持向量机 (soft margin)

我们之前讲是hard margin，引入soft的原因是因为：在hard margin的时候是认为这个数据是最理想的情况下是线性可分的，但是实际情况是比较复杂的，要么数据是不可分的，或者这个数据是可分的但是这个数据是存在噪声的。
在这种情况下，我们可以把噪声点剔除，可以利用软间隔支持向量机，通俗的来说就是可以容错，如下图：


软误差：${\sum }_{i=1}^m{\xi }_i$ ，所有数据元组松弛变量的和，C>0为惩罚参数， C的取值用来在数据拟合之间找到平衡。C越大，对目标函数的损失也越大，此时就暗示着你非常不愿意放弃这些离群点。

SVM涉及到的内容基本都讲到了，但是这里面还是很多细节，有兴趣的大家可以单独去了解，去学习。接下来讲一下SVM代码实现，最简答的算法smo算法。

SMO算法实现

SMO算法思路

其中$\lambda =[{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4\dots \dots {\lambda }_n]$

smo算法思路：选取两个优化变量${\lambda }_i,{\lambda }_j$ ，把剩余变量看做是固定的常数，进行优化目标函数，再迭代${\lambda }_i,{\lambda }_j$，再优化，直到所有的$\lambda$被优化完，迭代终止。

为什么要选取两个变量进行优化，因为有限制条件${\sum }_{i=1}^n{\lambda }_i{y}_i=0$，只优化一个的话，会不满足这个约束，所以我们选择两个，其中一个变化，另外一个也随着变化，然后把剩下的变量看成一个常数。

为了方便讲解，先举一对变量${\lambda }_1,{\lambda }_2$优化，令$K_{ij}=x_i^T{x}_j$，把含有${\lambda }_1,{\lambda }_2$变量的写出来，剩下的变量看做成常数用R表示:

$f({\lambda }_1,{\lambda }_2)=\frac{1}{2}{\lambda }_1^2{K}_{11}+\frac{1}{2}{\lambda }_2^2{K}_{22}+{\lambda }_1{\lambda }_2{y}_1{y}_2{K}_{12}+y_1{\lambda }_1{\sum }_{i=3}^n{y}_i{\lambda }_i{K}_{i1}+y_2{\lambda }_2{\sum }_{i=3}^n{y}_i{\lambda }_i{K}_{i2}?{\lambda }_1?{\lambda }_2+R$

于是就是一个二元函数优化：
$argma{x}_{{\lambda }_1,{\lambda }_2}f({\lambda }_1,{\lambda }_2)$

约束条件同样写成：

${\lambda }_1{y}_1+{\lambda }_2{y}_2+{\sum }_{i=3}^n{\lambda }_i{y}_i=0$，将${\sum }_{i=3}^n{\lambda }_i{y}_i$看成常数C

所以可以写成：${\lambda }_1{y}_1+{\lambda }_2{y}_2=?{\sum }_{i=3}^n{\lambda }_i{y}_i=C$，等式两边同时乘以$y_1$

因为$y_i=\pm 1$，所以$y_1^2=1$，因此${\lambda }_1=y_1(C?{\lambda }_2{y}_2)$，带入优化函数f中，就只有含有一个变量${\lambda }_2$，直接优化即可，令$v_1={\sum }_{i=3}^N{\lambda }_i{y}_i{K}_{i,1},v_2={\sum }_{i=3}^N{\lambda }_i{y}_i{K}_{i,2}$：

$$f({\lambda }_2)=\frac{1}{2}{\lambda }_1^2{K}_{11}+\frac{1}{2}{\lambda }_2^2{K}_{22}+{\lambda }_1{\lambda }_2{y}_1{y}_2{K}_{12}+y_1{\lambda }_1{v}_1+y_2{\lambda }_2{v}_2?{\lambda }_1?{\lambda }_2+C$$
对f函数求导得：
$\frac{?f({\lambda }_2)}{?{\lambda }_2}=?(K_{11}+K_{22}?2K_{12}){\lambda }_2+K_{11}C{y}_2?K_{12}C{y}_2+v_1{y}_2?v_2{y}_2?y_1{y}_2+y_2^2=0$

下面我们稍微对上式进行下变形，因为参数$\lambda$是不断更新迭代的，所以会有更新前${\lambda }_{old}$更新后${\lambda }_{new}$的值。
因为SVM对数据点的预测值为：$f(x)=w^T{x}_i+b={\sum }_{i=1}^N{\lambda }_i{y}_{i}K(x_i,x)+b$
则用$v_1,v_2$表示成：

$v_1={\sum }_{i=3}^N{\lambda }_i{y}_i{K}_{1i}=f(x_1)?{\lambda }_1{y}_1{K}_{11}?{\lambda }_2{y}_2{K}_{12}?b$

$v_2={\sum }_{i=3}^N{\lambda }_i{y}_i{K}_{2i}=f(x_2)?{\lambda }_1{y}_1{K}_{12}?{\lambda }_2{y}_2{K}_{22}?b$

已知${\lambda }_1=y_2(C?{\lambda }_2{y}_2)$，可得到：

$v_1?v_2=f(x_1)?f(x_2)?K_{11}C+k_{12}C+(K_{11}+K_{22}?2K_{12}){\lambda }_2{y}_2$

将$v_1?v_2$带入表达式$\frac{?f({\lambda }_2)}{?{\lambda }_2}$得到：

$\frac{?f({\lambda }_2)}{?{\lambda }_2}=?(K_{11}+K_{22}?2K_{12}){\lambda }_2^{new}+(K_{11}+K_{22}?2K_{12}){\lambda }_2^{old}+y_2[y_2?y_1+f(x_1)?f(x_2)]$

我们记$E_i$为svm预测值与真实值的误差：$E_i=f(x_i)?y_i$
令$\eta =K_{11}+K_{22}?2K_{12}$

$\frac{?f({\lambda }_2)}{?{\lambda }_2}=?\eta {\lambda }_2^{new}+\eta {\lambda }_2^{old}+y_2(E_1?E_2)=0$，得到：

${\lambda }_2^{new}={\lambda }_2^{old}+\frac{y_2(E_1?E_2)}{\eta }$

这样我们就可以通过不断迭代更新参数进行优化了。
以上是没有讨论了带约束条件的，但是在SVM中是有约束条件的，即：

${\lambda }_1{y}_1+{\lambda }_2{y}_2=?{\sum }_{i=3}^n{\lambda }_i{y}_i=\gamma$
$0\le {\lambda }_i\le C$

因为$y_i=\pm 1$，所以还要分情况讨论，$y_i=1$或者-1。以y1 y2 为例，${\lambda }_1=y_1(\gamma ?{\lambda }_2{y}_2)=y_1\gamma ?{\lambda }_2{y}_1{y}_2$

分情况讨论

第一种情况：当$y_1和y_2$异号时，${\lambda }_1=y_{1}C+{\lambda }_2$，这里面又分为两种情况，当y1等于1或者-1时，因此有：$(1){\lambda }_1=C+{\lambda }_2$ 和$(2){\lambda }_1=?C+{\lambda }_2$，这就变成了一个线性规划问题，把图画出来：

${\lambda }_1和{\lambda }_2$既要在矩形方框内，又要在直线上，因此可以写为：
下界：$L=max(0,{\lambda }_2^{old}?{\lambda }_1^{old})$
上界：$H=min(C,C+{\lambda }_2^{old}?{\lambda }_1^{old})$

第二种情况：当$y_1和y_2$同号时，同理，上下界：
下界：$L=max(0,{\lambda }_2+{\lambda }_1?C)$
上界：$H=min(C,{\lambda }_2+{\lambda }_1)$

根据${\lambda }_1^{old}{y}_1+{\lambda }_2^{old}{y}_2={\lambda }_1^{new}{y}_1+{\lambda }_2^{new}{y}_2$得到

${\lambda }_1^{new}={\lambda }_1^{old}+y_1{y}_2({\lambda }_2^{old}?{\lambda }_2^{new})$

这样，我们的一对${\lambda }_i,{\lambda }_j$就能优化更新了。

b阈值更新

更新了一对${\lambda }_i,{\lambda }_j$之后我们要计算b，因为有$f(x)=w^T{x}_i+b={\sum }_{i=1}^N{\lambda }_i{y}_{i}K(x_i,x)+b$

且优化样本都满足kkt条件

当${\lambda }_1^{new}$不在边界上，即在$0<{\lambda }_1^{new}<C$，根据kkt条件相应的数据点为支持向量(支持向量就是离超平面最近的那个点)，满足$y_1(w^{T}x+b)=1$,两边同乘以y1得${\sum }_{i=1}^N{\lambda }_i{y}_i{K}_{i1}+b=y_1$，因此得到

$b_1^{new}=y_1?{\sum }_{i=3}^N{\lambda }_i{y}_i{K}_{i1}?{\lambda }_1^{new}{y}_1{K}_{11}?{\lambda }_2^{new}{y}_2{K}_{21}$

根据${\sum }_{i=1}^N{\lambda }_i{y}_i{K}_{i1}+b=y_1$，前b中的两项等于：

$y_1?{\sum }_{i=3}^N{\lambda }_i{y}_i{K}_{i1}=?E_1+{\lambda }_1^{old}{y}_1{K}_{11}+{\lambda }_2^{old}{y}_2{K}_{21}+b^{old}$

当$0<{\lambda }_2^{new}<C$，可以得到：

$b_2^{new}=?E_2?y_1{K}_{12}({\lambda }_1^{new}?{\lambda }_1^{old})?y_2{K}_{22}({\lambda }_2^{new}?{\lambda }_2^{old})+b^{old}$

当b1和b2都有效的时候既是他们相等的时候，即$b^{new}=b_1^{new}=b_2^{new}$

当${\lambda }_1,{\lambda }_2$都在边界上时，SMO选择他们的重点作为新的值：

$b^{new}=\frac{b_1^{new}+b_2^{new}}{2}$

SMO代码实现

首先定义辅助函数，后面会用到

# 辅助函数
# 修剪alpha的值 L和H之间
def clip(alpha, L, H):
    if alpha < L:
        return L
    if alpha > H:
        return H
    else:
        return alpha

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

# 通过已知数据点和拉格朗日乘子获得分割超平面参数w
def get_w(alphas, dataset, labels):
    alphas, dataset, labels = np.array(alphas), np.array(dataset), np.array(labels)
    yx = labels.reshape(1,-1).T*np.array([1,1])*dataset  # reshape(1,-1)变成一行
    w = np.dot(yx.T, alphas)
    return w.tolist()


# smo代码
def simple_smo(dataset, labels, C, max_iter):
    # 初始化参数
    dataset = np.array(dataset)
    m, n = dataset.shape
    labels = np.array(labels)
    alphas = np.zeros(m)
    b = 0
    it = 0
    def f(x):
        "SVM分类器函数 y = w^Tx + b"
        # Kernel function vector.
        x = np.matrix(x).T
        data = np.matrix(dataset)
        ks = data*x
        # Predictive value.
        wx = np.matrix(alphas*labels)*ks
        fx = wx + b
        return fx[0, 0]
    # 迭代开始,选一对ai和aj进行迭代
    while it < max_iter:
        pair_changed = 0 # 迭代一次计数一次
        for i in range(m):
            a_i, x_i, y_i = alphas[i], dataset[i], labels[i]
            fx_i = f(x_i) # 分类器函数y=wTx+b
            E_i = fx_i - y_i # 预测值和真实值误差
            j = select_j(i, m)
            a_j, x_j, y_j = alphas[j], dataset[j], labels[j]
            fx_j = f(x_j)
            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对
            a_i_old, a_j_old = a_i, a_j
            a_j_new = a_j_old + y_j*(E_i - E_j)/eta
            # 对alpha进行修剪
            if y_i != y_j:
                L = max(0, a_j_old - a_i_old)
                H = min(C, C + a_j_old - a_i_old)
            else:
                L = max(0, a_i_old + a_j_old - C)
                H = min(C, a_j_old + a_i_old)
            # 一对aiaj迭代了之后，范围也要更新
            a_j_new = clip(a_j_new, L, H)
            a_i_new = a_i_old + y_i*y_j*(a_j_old - a_j_new)
            if abs(a_j_new - a_j_old) < 0.00001:
                #print('WARNING   alpha_j not moving enough')
                continue
            alphas[i], alphas[j] = a_i_new, a_j_new

            # 更新阈值b
            b_i = -E_i - y_i*K_ii*(a_i_new - a_i_old) - y_j*K_ij*(a_j_new - a_j_old) + b
            b_j = -E_j - y_i*K_ij*(a_i_new - a_i_old) - y_j*K_jj*(a_j_new - a_j_old) + b
            if 0 < a_i_new < C:
                b = b_i
            elif 0 < a_j_new < C:
                b = b_j
            else:
                b = (b_i + b_j)/2
            pair_changed += 1
            print('第{}次迭代  pair_changed:{}'.format(it, pair_changed))
        if pair_changed == 0:
            it += 1
        else:
            it = 0
        print('迭代次数: {}'.format(it))
    return alphas, b




if '__main__' == __name__:
    # 生成测试数据
    a = 200  # 样本数
    dataset, labels = make_blobs(n_samples=a, centers=2, n_features=2, random_state=2)
    labels[labels == 0] = -1

    # 使用简化版SMO算法优化SVM
    alphas, b = simple_smo(dataset,labels,0.6,20)
    # print('参数alphas: {}'.format(alphas))
    # 分类数据点
    classified_pts = {'+1':[], '-1':[]}
    for point, label in zip(dataset, labels):
        if label == 1.0:
            classified_pts['+1'].append(point)
        else:
            classified_pts['-1'].append(point)

    # print(classified_pts.items())

    fig = plt.figure() # 创建一个子图plt.figure(figsize = (5,5))
    ax = fig.add_subplot(111)
    # 绘制数据点
    for label, pts in classified_pts.items():
        pts = np.array(pts)
        ax.scatter(pts[:,0],pts[:,1], label = label)
    # 绘制分割线
    w = get_w(alphas,dataset,labels)
    x1, _ = max(dataset, key=lambda x: x[0]) #在样本点第1列数的最大值[1,3] =3
    x2, _ = min(dataset, key=lambda x: x[0])
    w1, w2 = w # w=(3,4)
    y1, y2 = (-b - w1*x1)/w2, (-b - w1*x2)/w2
    ax.plot([x1,x2],[y1,y2],c='black')
    # 绘制支持向量
    for i, alpha in enumerate(alphas): # 遍历alphas
        if abs(alpha) > 1e-3:
            x, y = dataset[i]
            ax.scatter([x],[y],s=50,c='r',marker="v" )
    plt.show()

下面运行代码，分别是迭代了50次和80次的效果，其中红色三角符就是支持向量，即离超平面最近的点：

有些图是借鉴，侵则删！