SVM算法原理与实现

一、线性可分SVM最大硬间隔

1.1 SVM简单介绍

  支持向量机（support vector machines， SVM) 是一种二分类模型。它的基本模型是定义在特征空间上的间隔最大的线性分类器。对于线性不可分数据，可以运用核技巧来分类，下面我们分别介绍线性可分模型和线性不可分模型。

1.2 符号约定

  假设给定一个特征空间上的训练数据集：
T = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x m , y m ) } T=\{(x_1,y_1),(x_2,y_2),...,(x_m,y_m)\} T={(x1​,y1​),(x2​,y2​),...,(xm​,ym​)}
其中，$x_i\in \chi =R^d$， y i ∈ Y = { + 1 , − 1 } y_i \in Y=\{+1,-1\} yi​∈Y={+1,−1}，$x_i$表示第i个特征向量。

图1 SVM分类示意图

1.3 函数间隔

  对于给定的训练数据集T和超平面$(w,b)$，定义超平面$(w,b)$关于T样本点$(x_i,y_i)$的函数间隔$\widehat{{\gamma }_i}$与几何间隔${\gamma }_i$分别为：
$${\hat{\gamma }}_i=y_i(w^T{x}_i+b)\text{\hspace{1em}(1.1)}$$
$${\gamma }_i=\frac{y_i(w^T{x}_i+b)}{||w||}\text{\hspace{1em}(1.2)}$$
  定义超平面$(w,b)$关于训练数据集T 的函数间隔为超平面$(w,b)$关于T中所有样本点$(x_i,y_i)$的函数间隔之最小值，即：
$$\hat{\gamma }=\underset{i=1,...,m}{\min }{\hat{\gamma }}_i\text{\hspace{1em}(1.3)}$$
  同理，定义超平面$(w,b)$关于训练数据集T 的函数间隔为：
$$\gamma =\underset{i=1,...,m}{\min }{\gamma }_i\text{\hspace{1em}(1.4)}$$
  通过上面的定义可以得到函数间隔和几何间隔存在如下关系：
$$\hat{\gamma }=\gamma ||w||\text{\hspace{1em}(1.5)}$$

1.4 硬间隔最大化

  支持向量机学习的基本想法是：求解能够正确划分训练数据集并且几何间隔最大的分离超平面。具体地， 这个问题可以表示为下面的约束最优化问题：
$$\begin{gathered} \underset{w,b}{\max }\gamma \\ s.t{y}_i\frac{w^T{x}_i+b}{||w||}\ge \gamma i=1,2,...,m\text{\hspace{1em}(1.6)} \end{gathered}$$
  考虑几何间隔和函数间的关系,公式$(1.6)$可以改写为：
$$\begin{gathered} \underset{w,b}{\max }\frac{\hat{\gamma }}{||w||} \\ s.t{y}_i(w^T{x}_i+b)\ge \hat{\gamma }i=1,2,...,m\text{\hspace{1em}(1.7)} \end{gathered}$$
  显然，函数间隔$\hat{\gamma }$的取值并不影响最优化问题的解，令$\hat{\gamma }=1$，注意到最大化$\frac{1}{||w||}$和最小化$\frac{1}{2}||w|{|}^2$是等价的，所以(1.7)式可写为：
$$\begin{gathered} \underset{w,b}{\min }\frac{1}{2}||w|{|}^2 \\ s.t{y}_i(w^T{x}_i+b)\ge 1i=1,2,...,m\text{\hspace{1em}(1.8)} \end{gathered}$$
  式(1.8)就是我们要求解的SVM分类问题，同时这也是一个二次凸规划问题，并将(1.8)式叫做SVM模型的原问题。这个问题的直接求解可以使用动态规划求解，但我们一般不这么做，因为样本数量过大时，时间复杂度太大。但是将这个问题转换为拉格朗日对偶问题，可以得到许多有趣的性质。我们先将这个原始问题转化为等价的拉格朗日无约束问题。
  首先我们写出拉格朗日函数：
$$Lagrange(w,b,\lambda )=\frac{1}{2}{w}^{T}w+\sum _{i=1}^m{\lambda }_i[1-y_i(w^T{x}_i+b)]{\lambda }_i\ge 0\text{\hspace{1em}(1.9)}$$
  则(1.8)式的优化就等价于下面的公式的优化
$$\begin{gathered} \underset{w,b}{\min }\underset{\lambda }{\max }Lagrange(w,b,\lambda ) \\ s.t{\lambda }_i\ge 0i=1,2,...,m\text{\hspace{1em}(1.10)} \end{gathered}$$
这里简单说明下(1.8)和(1.10)为什么等价，首先固定$(w,b)$来看$\max Lagrange(w,b,\lambda )$,这是一个关于$\lambda$的函数。

1. 对于满足约束条件的$(x_i,y_i)$ ，即$y_i(w^T{x}_i+b)\ge 1$，那么$1-y_i(w^T{x}_i+b)\le 0$，又因为${\lambda }_i\ge 0$ ，此时只有$\lambda =0$取最大，则$\max Lagrange(w,b,\lambda )=0$；
2. 对于不满足约束条件的$(x_i,y_i)$ ，即$y_i(w^T{x}_i+b)<1$，那么$1-y_i(w^T{x}_i+b)>0$，时只需要令$\lambda =+\infty$，则$\max Lagrange(w,b,\lambda )=+\infty$。
   综上，$\max Lagrange(w,b,\lambda )$的取值：
   $$\max Lagrange(w,b,\lambda )=\{\begin{array}{l}\frac{1}{2}{w}^{T}w{y}_i(w^T{x}_i+b)\ge 1\\ +\infty y_i(w^T{x}_i+b)<1\end{array}\text{\hspace{1em}(1.11)}$$
   $$\underset{w,b}{\min }\underset{\lambda }{\max }Lagrange(w,b,\lambda )=\underset{w,b}{\min }(\frac{1}{2}{w}^{T}w,+\infty )=\underset{w,b}{\min }\frac{1}{2}{w}^{T}w\text{\hspace{1em}(1.12)}$$
   现在回到公式(1.10)，对(1.10)使用拉格朗日对偶法，则对公式(1.10)的优化问题等价于：
   $$\begin{gathered} \underset{\lambda }{\max }\underset{w,b}{\min }Lagrange(w,b,\lambda ) \\ s.t{\lambda }_i\ge 0i=1,2,...,m\text{\hspace{1em}(1.13)} \end{gathered}$$
   对(1.13)的求解，需要先对拉格朗日函数对$w,b$求极小，在对$\lambda$求极大(关于什么条件下对偶问题有解，最优解和原始问题最优解相同的问题暂不说明)。

1. 求解${\min }_{w,b}L(w,b,\lambda )$
   将拉格朗日函数$L(w,b,\lambda )$分别对$w,b$求偏导数并令其等于0。
   $$\begin{gathered} {▽}_{w}L(w,b,\lambda )=w-\sum _{i=1}^m{\lambda }_i{y}_i{x}_i=0\Rightarrow w=\sum _{i=1}^m{\lambda }_i{y}_i{x}_i \\ {▽}_{b}L(w,b,\lambda )=\sum _{i=1}^m{\lambda }_i{y}_i=0\text{\hspace{1em}(1.14)} \end{gathered}$$
   将(1.14) 带入(1.13)得：
   $$\begin{gathered} \underset{\lambda }{\max }-\frac{1}{2}\sum _{i=1}^m\sum _{i=1}^m{\lambda }_i{\lambda }_j{y}_i{y}_j{x}_i^T{x}_j+\sum _{i=1}^m{\lambda }_i \\ s.t\sum _{i=1}^m{\lambda }_i{y}_i=0{\lambda }_i\ge 0i=1,2,...,m\text{\hspace{1em}(1.15)} \end{gathered}$$
   将式(1.16)的目标函数由求极大转换成求极小， 就得到下面与之等价的对偶最优化问题。
   $$\begin{gathered} \underset{\lambda }{\min }\frac{1}{2}\sum _{i=1}^m\sum _{i=1}^m{\lambda }_i{\lambda }_j{y}_i{y}_j{x}_i^T{x}_j-\sum _{i=1}^m{\lambda }_i \\ s.t\sum _{i=1}^m{\lambda }_i{y}_i=0{\lambda }_i\ge 0i=1,2,...,m\text{\hspace{1em}(1.16)} \end{gathered}$$
   如果${\lambda }^{\ast }=({\lambda }_1^{\ast },{\lambda }_2^{\ast },...,{\lambda }_m^{\ast })$ 是最优化问题(1.16)的解，则存在下标j，使得${\lambda }_j>0$，可按下式求得原始问题的最优解。
   $$\begin{gathered} w^{\ast }=\sum _{i=1}^m{\lambda }_i^{\ast }{y}_i{x}_i \\ {b}^{\ast }=y_j-\sum _{i=1}^m{\lambda }_i^{\ast }{y}_i{x}_i^T{x}_j\text{\hspace{1em}(1.17)} \end{gathered}$$

1.4 硬间隔最大化SMO

  目前，我们已经知道通过求解对偶问题的最优解${\lambda }^{\ast }$可以获得原始问题的最优解$w^{\ast },b^{\ast }$，那么接来下的问题就是：怎么求得对偶问题的最优解呢，我们可以通过最小序列优化(Sequential Minimal Optimization)求解。

SMO 算法是一种启发式算法， 其基本思路是： 如果所有变量的解都满足此最优化问题的 KK丁条件（ Karush-Kuhn-Tucker conditions)， 那么这个最优化问题的解就得到了. 因为 KKT 条件是该最优化问题的充分必要条件。否则， 选择两个变量， 固定其他变量，针对这两个变量构建一个二次规划问题。

注意，子问题的两个变最中只有一个是自由变量，假设选择${\lambda }_1,{\lambda }_2$这两个变量，${\lambda }_3,{\lambda }_4,..,{\lambda }_m$固定由等式(1.14)可知
$$\sum _{i=1}^m{\lambda }_i{y}_i=0$$
如果 ${\lambda }_1$确定，那么${\lambda }_2$也随之确定 ，所以子问题中同时更新两个变量。

1. 启发式选择第一个变量${\lambda }_1$
   SMO 称选择第 1 个变量的过程为外层循环，外层循环首先在整个训练集上遍历，计算每个样本是否违反KKT条件，如果违反就可以用来优化。单次扫描整个训练集之后，外层循环遍历所有拉格朗日乘数不为0的样本(非边界样本)，且检查非边界样本是否违反KKT条件，违反的将被用来做优化。外层循环扫描那些非边界样本，直到所有的非边界样本在误差ε内满足KKT条件，然后外层循环重新遍历训练集。外层循环交替进行在整个训练集上单次扫描和在非边界子集上多次扫描，直到整个训练集在误差ε内符合KKT条件，此时算法终止。(注：误差ε内很重要，一般设置为0.001)
2. 启发式选择第二个变量${\lambda }_2$
   SMO 称选择第 2 个变量的过程为内层循环，在要在内层循环中找第2个变量${\lambda }_2$。 第2个变量选择的标准是希望
   能使${\lambda }_2$有足够大的变化，选择使$|E_1-E_2|$最大。因为${\lambda }_1$已定，$E_1$也确定了，如果$E_1$是正的，那么选择最小的$E_i$作为$E_2$；如果$E_1$是负的，那么选择最大的$E_i$作为$E_2$。在特殊情况下， 如果内层循环通过以上方法选择的${\lambda }_2$不能使目标函数有足够的下降， 那么采用以下启发式规则继续选择${\lambda }_2$。遍历在间隔边界上的支持向量点，依次将其对应的变量作为${\lambda }_2$用，直到目标函数有足够的下降。若找不到合适的${\lambda }_2$那么遍历训练数据集；若仍找不到合适的${\lambda }_2$，则放弃第1个${\lambda }_1$，再通过外层循环寻求另外的${\lambda }_1$。
3. 更新${\lambda }_1,{\lambda }_2$
   先定义特征到结果的输出函数为$u_i$、残差$E_i$：
   $$u_i=w^T{x}_i+b\text{\hspace{1em}(1.18)}$$
   $$E_i=w^T{x}_i+b-y_i=u_i-y_i\text{\hspace{1em}(1.19)}$$
   SMO 的最优化问题的子问题可以写成
   $$\begin{gathered} \underset{{\lambda }_1,{\lambda }_2}{\min }W({\lambda }_1,{\lambda }_2)=\frac{1}{2}\left[\left(2y_1\sum _{i=3}^m{\lambda }_i{y}_i{x}_i{x}_1^T\right){\lambda }_1+x_1{x}_1^T{\lambda }_1^2+\left(2y_2\sum _{i=3}^m{\lambda }_i{y}_i{x}_i{x}_2^T\right){\lambda }_2+x_2{x}_2^T{\lambda }_2^2+2y_1{y}_2{x}_1{x}_2^T{\lambda }_1{\lambda }_2\right]-{\lambda }_1-{\lambda }_2+const \\ s.t{\lambda }_1{y}_1+{\lambda }_2{y}_2=-\sum _{i=3}^m{y}_i{\lambda }_i=\zeta \text{\hspace{1em}(1.20)} \end{gathered}$$
   令$K_{ij}$表示內积$x_i^T{x}_j$，将${\lambda }_2=(\zeta -{\lambda }_1{y}_1)y_2,y_i^2=1$带入$W({\lambda }_1,{\lambda }_2)$得到只含有${\lambda }_1$的目标函数：
   $$\begin{array}{rl}W({\lambda }_1)& =W({\lambda }_1,(\zeta -{\lambda }_1{y}_1)y_2)\\ & ={\lambda }_1{y}_1\sum _{i=3}^m{\lambda }_i{y}_i{K}_{i1}+\frac{1}{2}{K}_{11}{\lambda }_1^2+(\zeta -{\lambda }_1{y}_1)\sum _{i=3}^m{\lambda }_i{y}_i{K}_{i2}+\frac{1}{2}{K}_{22}(\zeta -{\lambda }_1{y}_1{)}^2+y_1{K}_{12}{\lambda }_1(\zeta -{\lambda }_1{y}_1)-{\lambda }_1-(\zeta -{\lambda }_1{y}_1)y_2+const\\ & =\frac{1}{2}(K_{11}+K_{22}-2K_{12}){\lambda }_1^2+\left[y_1\sum _{i=3}^m{\lambda }_i{y}_i{K}_{i1}-y_1\sum _{i=3}^m{\lambda }_i{y}_i{K}_{i2}-y_1\zeta K_{22}+y_1\zeta K_{12}-1+y_1{y}_2\right]{\lambda }_1+const\end{array}\text{\hspace{1em}(1.21)}$$
   对${\lambda }_1$求偏导，令其等于0，并带入$y_1{\lambda }_1^{old}+y_2{\lambda }_2^{old}=\zeta$ ：
   $$\frac{\partial W({\lambda }_1)}{\partial {\lambda }_1}=(K_{11}+K_{22}-2K_{12}){\lambda }_1+\left[y_1\sum _{i=3}^m{\lambda }_i{y}_i{K}_{i1}-y_1\sum _{i=3}^m{\lambda }_i{y}_i{K}_{i2}-y_1\zeta K_{22}+y_1\zeta K_{12}-1+y_1{y}_2\right]=0\text{\hspace{1em}(1.22)}$$
   $$\begin{array}{rl}(2K_{12}-K_{11}-K_{22}){\lambda }_1& =y_1\left[y_2-y_1+\zeta K_{12}-\zeta K_{22}+\sum _{i=3}^m{\lambda }_i{y}_i{K}_{i1}-\sum _{i=3}^m{\lambda }_i{y}_i{K}_{i2}\right]\\ & =y_1\left[y_2-y_1+\zeta K_{12}-\zeta K_{22}+\sum _{i=1}^m{\lambda }_i{y}_i{K}_{i1}-{\lambda }_1{y}_1{K}_{11}-{\lambda }_2{y}_2{K}_{21}-\sum _{i=0}^m{\lambda }_i{y}_i{K}_{i2}+{\lambda }_1{y}_1{K}_{12}+{\lambda }_2{y}_2{K}_{22}\right]\\ & =y_1\left[y_2-y_1+\zeta K_{12}-\zeta K_{22}+w^T{x}_1-{\lambda }_1^{old}{y}_1{K}_{11}-{\lambda }_2^{old}{y}_2{K}_{21}-w^T{x}_2+{\lambda }_1^{old}{y}_1{K}_{12}+{\lambda }_2^{old}{y}_2{K}_{22}\right]\\ & =y_1\left[E_1-E_2+2y_1{K}_{12}{\lambda }_1^{old}-y_1{K}_{11}{\lambda }_1^{old}-y_1{K}_{22}{\lambda }_1^{old}\right]\\ & =(2K_{12}-K_{11}-K_{22}){\lambda }_1^{old}+y_1(E_1-E_2)\end{array}\text{\hspace{1em}(1.23)}$$
   令$\eta =2K_{12}-K_{11}-K_{22}$，由(1.23)得到${\lambda }_1$：
   $${\lambda }_1^{new}={\lambda }_1^{old}+\frac{y_1(E_1-E_2)}{\eta }\text{\hspace{1em}(1.24)}$$
   $${\lambda }_2^{new}={\lambda }_2^{old}+y_1{y}_2({\lambda }_1^{old}-{\lambda }_1^{new})\text{\hspace{1em}(1.25)}$$
   此时得到是未经修剪的${\lambda }_1$，需满足约束${\lambda }_i\ge 0$，所以需要确定${\lambda }_1$的上下界。
   $$y_1{\lambda }_1^{old}+y_2{\lambda }_2^{old}=y_1{\lambda }_1^{new}+y_2{\lambda }_2^{new}=const\text{\hspace{1em}(1.26)}$$
   由式(1.25)得到${\lambda }_1$上下界的取值范围：
   $$\begin{gathered} y_1=y_2,L=0;H={\lambda }_1^{old}+{\lambda }_2^{old} \\ {y}_1!=y_2,L=max(0,{\lambda }_1^{old}-{\lambda }_2^{old});H=+\infty \text{\hspace{1em}(1.27)} \end{gathered}$$
4. 更新$b$
   在每次完成两个变量的优化后， 都要重新计算阈值b，根据KKT条件
   $$KKT条件\{\begin{array}{l}{\lambda }_i\ge 0\\ y_i(w^T{x}_i+b)-1\ge 0\\ {\lambda }_i[y_i(w^T{x}_i+b)-1]=0\end{array}\text{\hspace{1em}(1.28)}$$
   $$违反KKT条件\{\begin{array}{l}1、y_i(w^T{x}_i+b)-1=y_i{E}_i<0\\ 2、y_i(w^T{x}_i+b)-1=y_i{E}_i>0,{\lambda }_i>0,\end{array}$$
   当${\lambda }_1^{new}>0$时：
   $$\begin{gathered} w^T{x}_1+b_1^{new}=y_1\Rightarrow b_1^{new}=y_1-\sum _{i=0}^m{\lambda }_i{y}_i{x}_i{x}_1^T \\ {b}_1^{new}=y_1-\sum _{j=3}^m{\lambda }_i{y}_i{x}_i{x}_1^T-{\lambda }_1^{new}{y}_1{K}_{11}-{\lambda }_2^{new}{y}_2{K}_{21} \\ {E}_1=\sum _{j=3}^m{\lambda }_i{y}_i{x}_i{x}_1^T+{\lambda }_1^{old}{y}_1{K}_{11}+{\lambda }_2^{old}{y}_2{K}_{21}+b^{old}-y_1 \\ {b}_1^{new}=-E_1+{\lambda }_1^{old}{y}_1{K}_{11}+{\lambda }_2^{old}{y}_2{K}_{21}+b^{old}-{\lambda }_1^{new}{y}_1{K}_{11}-{\lambda }_2^{new}{y}_2{K}_{21} \\ {b}_1^{new}=-E_1-y_1{K}_{11}({\lambda }_1^{new}-{\lambda }_1^{old})-y_2{K}_{21}({\lambda }_2^{new}-{\lambda }_2^{old})+b^{old}\text{\hspace{1em}(1.29)} \end{gathered}$$
   同理，当${\lambda }_2^{new}>0$时：
   $$b_2^{new}=-E_2-y_1{K}_{21}({\lambda }_1^{new}-{\lambda }_1^{old})-y_2{K}_{22}({\lambda }_2^{new}-{\lambda }_2^{old})+b^{old}\text{\hspace{1em}(1.30)}$$
   如果 ${\lambda }_1^{new},{\lambda }_2^{new}$，同时满足条件${\lambda }_i^{new}>0$，那么$b_1^{new}=b_2^{new}$；如果 ${\lambda }_1^{new},{\lambda }_2^{new}$是0，那么$b_1^{new}$和$b_2^{new}$以及它们之间的数都是符合 KKT 条件的阈值， 这时选择它们的中点作为 $b^{new}$。
5. 更新$E_i$
   $$E_j^{new}=\sum _S{y}_i{\lambda }_{i}Kij+b^{new}-y_j\text{\hspace{1em}(1.31)}$$
   其中S是所有支持向量$x_i$的集合，根据KKT条件，支持向量的点${\lambda }_i\ge 0$，不是支持向量的点${\lambda }_i=0$，所以式（1.31）可以写成（1.32），只是计算量增加，其结果不变。
   $$E_j^{new}=\sum _{i=1}^m{y}_i{\lambda }_{i}Kij+b^{new}-y_j\text{\hspace{1em}(1.32)}$$

1.5 硬间隔SMO只代码实现

#外层循环 选择alpha1, KKT条件的检查在internal_loop进行
    def outer_loop(self):
        numChanged = 0
        examineAll = 1
        while numChanged > 0 or examineAll:
            numChanged = 0
            if examineAll:
            	# 1 首次在整个训练集上遍历单次扫描，2 所有的非边界样本满足KKT，然后重新遍历整个样本集
                for i in range(self.alphas.shape[0]):
                    numChanged += self.internal_loop(i)
            else:
                # 外层循环遍历所有拉格朗日乘数不为0的样本(非边界样本)，在非边界子集上多次扫描
                for i in np.where(self.alphas != 0)[0]:
                    numChanged += self.internal_loop(i)
            if examineAll:
                examineAll = 0
            elif numChanged == 0:
                #表示外层循环交替进行在整个训练集上单次扫描    和   在非边界子集上多次扫描
                examineAll = 1

    # 内层循环
    def internal_loop(self, i1):
        yi = self.y[i1]
        ei = self.errors[i1]
        alphai = self.alphas[i1]
        #在误差范围内0.001，检查第一个点alpha1是否违反KKT条件 
        if (yi * ei < -self.tol) or (yi * ei > self.tol and alphai > 0):
           # 在边界上的点选择|E1-E2|最大
            if len(self.alphas[(self.alphas != 0)]) > 1:
                if ei > 0:
                    i2 = np.argmin(self.errors)
                else:
                    i2 = np.argmax(self.errors)
                step_result = self.take_step(i1, i2)
                if step_result:
                    return 1
            # 若上一步不能更新，就从 0<alpha_1 随机选择一个j，能够更新alpha就更新。
            for i2 in np.roll(np.where(self.alphas != 0)[0], np.random.choice(self.m)):
                step_result = self.take_step(i1, i2)
                if step_result:
                    return 1
            # 若上一步还是不能更新，不考虑范围，从整个数据集中选择一个样本j
            for i2 in np.roll(np.arange(self.m), np.random.choice(np.arange(self.m))):
                step_result = self.take_step(i1, i2)
                if step_result:
                    return 1
        return 0

    # 更新参数alpha1，alpha2，对应文中的lambda1、lambda2
    def take_step(self, i1, i2):
        if i1 == i2:
            return 0
        alph1 = self.alphas[i1]
        alph2 = self.alphas[i2]
        y1 = self.y[i1]
        y2 = self.y[i2]
        E1 = self.errors[i1]
        E2 = self.errors[i2]
        s = y1 * y2
        # 计算 alpha1的上下界
        if y1 != y2:
            L = max(0, alph1 - alph2)
        elif y1 == y2:
            H = alph1 + alph2

        k11 = self.X[i1] @ self.X[i1].T
        k12 = self.X[i1] @ self.X[i2].T
        k22 = self.X[i2] @ self.X[i2].T
        eta = 2 * k12 - k11 - k22
        #计算新的alpha1
        if eta < 0:
            a1 = alph1 + y1 * (E1 - E2) / eta
            if y1 != y2:
                a1 = max(L, a1)
            else:
                a1 = min(H, a1)
        else:
            return 0

        if a1 < 1e-5:
            a1 = 0.0

        #更新过小，则放弃更新
        if np.abs(a1 - alph1) < self.eps*0.01:
            return 0
        a1 = np.round(a1, decimals=5)
        a2 = alph2 + s * (alph1 - a1)
        a2 = np.round(a2, decimals=5)
       
        b1 = -E1 - y1 * (a1 - alph1) * k11 - y2 * (a2 - alph2) * k12 + self.b
        b2 = -E2 - y1 * (a1 - alph1) * k12 - y2 * (a2 - alph2) * k22 + self.b
        b_new = (b1 + b2) * 0.5

        self.alphas[i1] = a1
        self.alphas[i2] = a2
        self.b = b_new
        self.errors = self.decision_function(self.X) - self.y
        return 1

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs


class linersvm:
    def __init__(self, X_train, y, alphas, errors, tol=0.001, eps=0.001, b=0.0):
        self.X = X_train
        self.y = y
        self.alphas = alphas
        self.tol = tol
        self.eps = eps
        self.b = b
        self.errors = errors
        self.m = len(self.X)

    def decision_function(self, x_test):
        return (self.alphas * self.y) @ (self.X @ x_test.T) + self.b
    
    def showClassifer(self):
        plt.scatter(self.X[:, 0], self.X[:, 1], c=self.y, s=30, cmap=plt.cm.Paired)
        xrange = np.linspace(self.X[:, 0].min(), self.X[:, 0].max(), 1000)
        yrange = np.linspace(self.X[:, 1].min(), self.X[:, 1].max(), 1000)
        XX, YY = np.meshgrid(xrange, yrange)
        xy = np.vstack([XX.ravel(), YY.ravel()]).T
        Z = self.decision_function(xy).reshape(XX.shape)
        #画决策边界
        plt.contour(XX, YY, Z, colors='k', levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
        #画支持向量的点
        mask = np.round(self.alphas, decimals=3) != 0.0
        plt.scatter(self.X[mask, 0], self.X[mask, 1], s=100, linewidth=1, facecolors='none', edgecolors='k')
        plt.show()

if __name__ == '__main__':
    sample_numbers = 500
    samples, labels = make_blobs(n_samples=sample_numbers, centers=2, n_features=2, random_state=1)
    labels[labels == 0] = -1
    model = linersvm(samples, labels, np.zeros(sample_numbers), np.zeros(sample_numbers))
    model.errors = model.decision_function(samples) - labels
    model.outer_loop()
    model.showClassifer()

实验结果：

二、线性SVM最大软间隔

2.1 最大软间隔算法介绍

  对于数据集线性不可分的时候，不等式约束并不能都成立，怎么才能将它扩展到线性不可分问题呢？ 这就需要修改硬间隔最大化， 使其成为软间隔最大化。
  线性不可分意味着某些样本点$(x_i,y_i)$不能满足函数间隔大于等于1的约束条件。为了解决这个问题，可以对每个样本点$(x_i,y_i)$引进一个松弛变量${\xi }_i\ge 0$使函数间隔加上松弛变量大于等于1。约束条件变为：

$$y_i(w^T{x}_i+b)\ge 1-{\xi }_i\text{\hspace{1em}(2.1)}$$
同时， 对每个松弛变量${\xi }_i$，支付一个代价${\xi }_i$函数由原来的$\frac{1}{2}||w|{|}^2$变成：
$$\frac{1}{2}||w|{|}^2+C\sum _{i=1}^m{\xi }_i\text{\hspace{1em}(2.2)}$$
这里，C>0称为惩罚参数，C值大时对误分类的惩罚增大，C值小时对误分类的惩罚减小。最小化目标函数(2.2)包含两层含义：使$\frac{1}{2}||w|{|}^2$尽量小即间隔尽量大，同时使误分类点的个数尽量小，C是调和二者的系数。软间隔最大化的原始问题：
$$\begin{gathered} \underset{w,b,\xi }{\min }\frac{1}{2}||w|{|}^2+C\sum _{i=1}^m{\xi }_i \\ s.t{y}_i(w^T{x}_i+b)\ge 1-{\xi }_{i}i=1,2,...,m \\ {\xi }_i\ge 0i=1,2,...,m\text{\hspace{1em}(2.3)} \end{gathered}$$
原始问题的拉格朗日函数是：
$$L(w,b,\xi ,\lambda ,\mu )=\frac{1}{2}||w|{|}^2+C\sum _{i=1}^m{\xi }_i+\sum _{i=1}^m{\lambda }_i[1-{\xi }_i-y_i(w^T{x}_i+b)]-\sum _{i=1}^m{\mu }_i{\xi }_i\text{\hspace{1em}(2.4)}$$
原始问题等价于对拉格朗日函数的极小极大：
$$\begin{gathered} \underset{w,b,\xi }{\min }\underset{\lambda ,\mu }{\max }L(w,b,\xi ,\lambda ,\mu ) \\ s.t{\lambda }_i\ge 0,{\mu }_i\ge 0i=1,2,...,m\text{\hspace{1em}(2.5)} \end{gathered}$$
拉格朗日函数的极小极大的对偶问题是极大极小：
$$\begin{gathered} \underset{\lambda ,\mu }{\max }\underset{w,b,\xi }{\min }L(w,b,\xi ,\lambda ,\mu ) \\ s.t{\lambda }_i\ge 0,{\mu }_i\ge 0i=1,2,...,m\text{\hspace{1em}(2.6)} \end{gathered}$$
求$L(w,b,\xi ,\lambda ,\mu )$对$w,b,\xi$的极小
$$\begin{gathered} \frac{\partial L}{\partial w}=w-\sum _{i=1}^m{\lambda }_i{y}_i{x}_i=0\Rightarrow w=\sum _{i=1}^m{\lambda }_i{y}_i{x}_i \\ \frac{\partial L}{\partial b}=-\sum _{i=1}^m{\lambda }_i{y}_i=0\Rightarrow \sum _{i=1}^m{\lambda }_i{y}_i=0 \\ \frac{\partial L}{\partial {\xi }_i}=C-{\lambda }_i-{\mu }_i=0\Rightarrow C={\lambda }_i+{\mu }_i\text{\hspace{1em}(2.7)} \end{gathered}$$
将(2.7)带入(2.6)得：
$$\underset{w,b,\xi }{\max }L(w,b,\xi ,\lambda ,\mu )=-\frac{1}{2}\sum _{i=1}^m\sum _{i=1}^m{\lambda }_i{\lambda }_j{y}_i{y}_j(x_i\bullet x_j)+\sum _{i=1}^m{\lambda }_i\text{\hspace{1em}(2.8)}$$
再对$\max {\min }_{w,b,\xi }L(w,b,\xi ,\lambda ,\mu )$求$\lambda$极大,即得对偶问题：
$$\begin{gathered} \underset{\lambda }{\max }-\frac{1}{2}\sum _{i=1}^m\sum _{i=1}^m{\lambda }_i{\lambda }_j{y}_i{y}_j(x_i\bullet x_j)+\sum _{i=1}^m{\lambda }_i \\ s.t\sum _{i=1}^m{\lambda }_i{y}_i=0{\lambda }_i\ge 0,{\mu }_i\ge 0i=1,2,...,m \\ C={\lambda }_i+{\mu }_i\text{\hspace{1em}(2.9)} \end{gathered}$$
将式(2.7)的目标函数由求极大转换成求极小， 就得到下面与之等价的对偶最优化问题。
$$\begin{gathered} \underset{\lambda }{\min }\frac{1}{2}\sum _{i=1}^m\sum _{i=1}^m{\lambda }_i{\lambda }_j{y}_i{y}_j(x_i\bullet x_j)-\sum _{i=1}^m{\lambda }_i \\ s.t\sum _{i=1}^m{\lambda }_i{y}_i=0{\lambda }_i\ge 0,{\mu }_i\ge 0i=1,2,...,m \\ C={\lambda }_i+{\mu }_i\text{\hspace{1em}(2.10)} \end{gathered}$$
我们可以看到对于(2.10) 式目标函数的优化函数与最大硬间隔一样，多了一个约束条件。所以对$\lambda ,b,E_i$的更新与1.4节稍微不同，即$\lambda$的上下界改变了，b的更新条件改变了，KKT条件也稍有不同。
$$KKT条件\{\begin{array}{l}{\lambda }_i\ge 0,{\mu }_i\ge 0,{\xi }_i\ge 0\\ {\mu }_i{\xi }_i=0\\ C={\lambda }_i+{\mu }_i\\ y_i(w^T{x}_i+b)-1+{\xi }_i\ge 0\\ {\lambda }_i[y_i(w^T{x}_i+b)-1+{\xi }_i]=0\end{array}\text{\hspace{1em}(2.11)}$$
由(2.10)可以得出下面的结论：
$$\{\begin{array}{l}1、{\lambda }_i=0\Rightarrow {\mu }_i=C,{\xi }_i=0,y_i(w^T{x}_i+b)-1\ge 0\\ 2、0<{\lambda }_i<C\Rightarrow {\mu }_i>0,{\xi }_i=0,y_i(w^T{x}_i+b)-1=0\\ 3、{\lambda }_i=C\Rightarrow {\mu }_i=0,{\xi }_i\ge 0,y_i(w^T{x}_i+b)-1\le 0\end{array}\text{\hspace{1em}(2.12)}$$
则违反KKT条件：
$$违反KKT条件\{\begin{array}{l}1、y_i(w^T{x}_i+b)-1=y_i{E}_i\le 0,{\lambda }_i<C\\ 2、y_i(w^T{x}_i+b)-1=y_i{E}_i\ge 0,{\lambda }_i>0\\ 3、y_i(w^T{x}_i+b)-1=y_i{E}_i=0,{\lambda }_i=0或{\lambda }_i=C\end{array}\text{\hspace{1em}(2.13)}$$
确定${\lambda }_1$的上下界，根据下面的几个等式：
$$\begin{gathered} y_1{\lambda }_1^{old}+y_2{\lambda }_2^{old}=y_1{\lambda }_1^{new}+y_2{\lambda }_2^{new}=const \\ 0\le {\lambda }_1\le C,0\le {\lambda }_2\le C\text{\hspace{1em}(2.14)} \end{gathered}$$
如果$y_1=y_2$:
$$\begin{gathered} {\lambda }_1^{old}+{\lambda }_2^{old}={\lambda }_1^{new}+{\lambda }_2^{new}=const \\ {\lambda }_2^{new}={\lambda }_1^{old}+{\lambda }_2^{old}-{\lambda }_1^{new} \\ \Rightarrow 0\le {\lambda }_1^{old}+{\lambda }_2^{old}-{\lambda }_1^{new}\le C \\ \Rightarrow {\lambda }_1^{old}+{\lambda }_2^{old}-C\le {\lambda }_1^{new}\le {\lambda }_1^{old}+{\lambda }_2^{old} \\ L:max(0,{\lambda }_1^{old}+{\lambda }_2^{old}-C)H:min(C,{\lambda }_1^{old}+{\lambda }_2^{old}) \end{gathered}$$
如果$y_1!=y_2$:
$$\begin{gathered} {\lambda }_1^{old}-{\lambda }_2^{old}={\lambda }_1^{new}-{\lambda }_2^{new}=const \\ {\lambda }_2^{new}={\lambda }_1^{new}-{\lambda }_1^{old}+{\lambda }_2^{old} \\ \Rightarrow 0\le {\lambda }_1^{new}-{\lambda }_1^{old}+{\lambda }_2^{old}\le C \\ \Rightarrow {\lambda }_1^{old}-{\lambda }_2^{old}\le {\lambda }_1^{new}\le C+{\lambda }_1^{old}-{\lambda }_2^{old} \\ L:max(0,{\lambda }_1^{old}-{\lambda }_2^{old})H:min(C,C+{\lambda }_1^{old}-{\lambda }_2^{old}) \end{gathered}$$
综上${\lambda }_1^{new}$的上下界如下：
$$\{\begin{array}{l}{y}_1=y_{2}L:max(0,{\lambda }_1^{old}+{\lambda }_2^{old}-C)H:min(C,{\lambda }_1^{old}+{\lambda }_2^{old})\\ y_1!=y_{2}L:max(0,{\lambda }_1^{old}-{\lambda }_2^{old})H:min(C,C+{\lambda }_1^{old}-{\lambda }_2^{old})\end{array}\text{\hspace{1em}(2.15)}$$
${\lambda }^{new}$参数更新：
$$\begin{gathered} {\lambda }_1^{new}={\lambda }_1^{old}+\frac{y_1(E_1-E_2)}{\eta } \\ {\lambda }_2^{new}={\lambda }_2^{old}+y_1{y}_2({\lambda }_1^{old}-{\lambda }_1^{new}) \\ \text{\hspace{1em}(2.16)} \end{gathered}$$
更新b，当$0<{\lambda }_1<C$时，由KKT 条件(2.11)可知：
$$\sum _{i=1}^m{\lambda }_i{y}_i{K}_{i1}+b=y_1$$
参考(1.29)解得$b_1^{new}$：
$$b_1^{new}=-E_1-y_1{K}_{11}({\lambda }_1^{new}-{\lambda }_1^{old})-y_2{K}_{21}({\lambda }_2^{new}-{\lambda }_2^{old})+b^{old}\text{\hspace{1em}(2.17)}$$
同理，如果$0<{\lambda }_2<C$，那么：
$$b_2^{new}=-E_2-y_1{K}_{21}({\lambda }_1^{new}-{\lambda }_1^{old})-y_2{K}_{22}({\lambda }_2^{new}-{\lambda }_2^{old})+b^{old}\text{\hspace{1em}(2.18)}$$
$E_i$的更新与(1.32)一样：
$$E_j^{new}=\sum _{i=1}^m{y}_i{\lambda }_{i}Kij+b^{new}-y_j\text{\hspace{1em}(2.19)}$$

2.2 最大软间隔代码实现

    # 内层循环和最大硬间隔一样，这个函数是最大软间隔的外层循环
    def internal_loop(self, i1):
        yi = self.y[i1]
        ei = self.errors[i1]
        alphai = self.alphas[i1]
        if (yi * ei < -self.tol and alphai < self.C) or (yi * ei > self.tol and alphai > 0):
            if len(self.alphas[(self.alphas != 0) & (self.alphas != self.C)]) > 1:
                if ei > 0:
                    i2 = np.argmin(self.errors)
                else:
                    i2 = np.argmax(self.errors)
                step_result = self.take_step(i1, i2)
                if step_result:
                    return 1
            # 若上一步不能更新，就从 0<alpha_i<C 随机选择一个j，能够更新alpha就更新。
            for i2 in np.roll(np.where((self.alphas != 0) & (self.alphas != self.C))[0],
                              np.random.choice(self.m)):
                step_result = self.take_step(i1, i2)
                if step_result:
                    return 1
            # 若上一步还是不能更新，不考虑范围，从整个数据集中选择一个样本j
            for i2 in np.roll(np.arange(self.m), np.random.choice(np.arange(self.m))):
                step_result = self.take_step(i1, i2)
                if step_result:
                    return 1
        return 0

   # 更新参数alpha1，alpha2，对应文中的lambda1、lambda2
    def take_step(self, i1, i2):
        if i1 == i2:
            return 0
        alph1 = self.alphas[i1]
        alph2 = self.alphas[i2]
        y1 = self.y[i1]
        y2 = self.y[i2]
        E1 = self.errors[i1]
        E2 = self.errors[i2]
        s = y1 * y2
        
        if (y1 != y2):
            L = max(0, alph1 - alph2)
            H = min(self.C, self.C + alph1 - alph2)
        elif (y1 == y2):
            L = max(0, alph1 + alph2 - self.C)
            H = min(self.C, alph1 + alph2)
        if (L == H):
            return 0

        k11 = self.X[i1] @ self.X[i1].T
        k12 = self.X[i1] @ self.X[i2].T
        k22 = self.X[i2] @ self.X[i2].T
        eta = 2 * k12 - k11 - k22
        
        if eta < 0:
            a1 = alph1 + y1 * (E1 - E2) / eta
            # Clip a2 based on bounds L & H
            if L < a1 < H:
                a1 = a1
            elif a1 <= L:
                a1 = L
            elif a1 >= H:
                a1 = H
        else:
            return 0
            
        if a1 < 1e-4:
            a1 = 0.0
        elif a1 > (self.C - 1e-4):
            a1 = self.C
         
        if (np.abs(a1 - alph1) < self.eps * (a1 + alph1 + self.eps)):
            return 0
        a1 = np.round(a1, decimals=5)
        a2 = alph2 + s * (alph1 - a1)
        a2 = np.round(a2, decimals=5)

        b1 = -E1 - y1 * (a1 - alph1) * k11 - y2 * (a2 - alph2) * k12 + self.b
        b2 = -E2 - y1 * (a1 - alph1) * k12 - y2 * (a2 - alph2) * k22 + self.b

        if 0 < a1 < self.C:
            b_new = b1
        elif 0 < a2 < self.C:
            b_new = b2
        else:
            b_new = (b1 + b2) * 0.5

        self.alphas[i1] = a1
        self.alphas[i2] = a2
        
        self.b = b_new
        self.errors = self.decision_function(self.X) - self.y
        return 1

if __name__ == '__main__':
	mean = [4, 4]
    cov = [[1, 0], [0, 2]]
    smpcount = 15
    samples1 = np.random.multivariate_normal(mean, cov, smpcount, "raise")
    samples1 = np.c_[samples1, np.ones([smpcount, 1])]
    mean = [7, 7]
    cov = [[2, 0], [0, 1]]
    samples2 = np.random.multivariate_normal(mean, cov, smpcount, "raise")
    samples2 = np.c_[samples2, np.ones([smpcount, 1])*(-1)]
    samples = np.vstack([samples1, samples2])
    np.random.shuffle(samples)
    model = linersvm(samples[:, :-1], samples[:, -1], 0.1, np.zeros(smpcount*2), np.zeros(smpcount*2))
    model.errors = model.decision_function(samples[:, :-1]) - samples[:, -1]
    model.outer_loop()
    model.showClassifer()

2.3 实验结果

3 非线性SVM

  前面介绍的线性SVM，是假设样本集是线性可分或大部分数据线性可分，而实际情况往往是线性不可分。即在样本空间不存在一个超平面能够正确划分样本类别。面对这样的情况，我们可以将样本投影到高维特征空间，是的样本在这个高维特征空间是线性可分的。则划分超平面可以表示为：
$$f(x)=w^T\varphi (x)+b\text{\hspace{1em}(3.1)}$$
关于$w,b$的原始问题为：
$$\begin{gathered} \underset{w,b}{\min }\frac{1}{2}||w|{|}^2 \\ s.t{y}_i(w^T\varphi (x_i)+b)\ge 1i=1,2,...,m\text{\hspace{1em}(3.2)} \end{gathered}$$
对偶问题变为：
$$\begin{gathered} \underset{\lambda }{\min }\frac{1}{2}\sum _{i=1}^m\sum _{i=1}^m{\lambda }_i{\lambda }_j{y}_i{y}_j\varphi (x_i{)}^T\varphi (x_j)-\sum _{i=1}^m{\lambda }_i \\ s.t\sum _{i=1}^m{\lambda }_i{y}_i=0{\lambda }_i\ge 0i=1,2,...,m\text{\hspace{1em}(3.3)} \end{gathered}$$
其中，$\varphi (x_i{)}^T\varphi (x_j)$表示映射在特征里的內积，直接计算$\varphi (x_i{)}^T\varphi (x_j)$是比较困难。所以我们引入了核函数：
$$\kappa (x_i,x_j)=\varphi (x_i{)}^T\varphi (x_j)\text{\hspace{1em}(3.4)}$$
即，$x_i,x_j$在特征空间的內积等于在样本的原空间中通过$\kappa (.,.)$计算得到，称$\kappa (x,z)$为核函数。下面是集中常用的核函数：

名称	表达式	参数
线性核	$\kappa (x_i,x_j)=x_i^T{x}_j$	–
多项式核	$\kappa (x_i,x_j)=(x_i^T{x}_j{)}^T$	$d\ge 1$为多项式的次数
高斯核	$\kappa (x_i,x_j)=exp(-\frac{\Vert x_i-x_j{\Vert }^2}{2{\sigma }^2})$	$\sigma >0$
拉普拉斯核	$\kappa (x_i,x_j)=exp(-\frac{\Vert x_i-x_j\Vert }{\sigma })$	$\sigma >0$
Sigmoid核	$\kappa (x_i,x_j)=tanh(\beta x_i^T{x}_j+\theta )$	tanh为双曲正切函数 $\beta >0,\theta <0$

核函数的性质：

• 若${\kappa }_1$和${\kappa }_2$为核函数，对于任意${\lambda }_1,{\lambda }_2$的线性组合${\lambda }_1{\kappa }_1+{\lambda }_2{\kappa }_2$也是核函数。
• 若${\kappa }_1$和${\kappa }_2$为核函数，核函数的乘积${\kappa }_1(x,z){\kappa }_2(x,z)$也是核函数。
• 若${\kappa }_1$为核函数，对任意函数$g(x)$，有$g(x){\kappa }_1(x,z)g(z)$也是核函数。

参考

【1】《统计学习方法》 李航
【2】《机器学习》周志华
【3】Python3：《机器学习实战》之支持向量机（3）完整版SMO