SVM and SMO（Sequencial Minimal Optimization）

最新推荐文章于 2022-11-29 14:54:37 发布

转载最新推荐文章于 2022-11-29 14:54:37 发布 · 575 阅读

文章标签：

#SMO #SVM

ML 专栏收录该内容

14 篇文章

订阅专栏

本文深入解析支持向量机(SVM)原理，涵盖线性与非线性分类，介绍SMO算法优化过程，包括参数选择、更新及Python实现，适合机器学习初学者。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

主要参考

1.SVM_SMO-Python代码实现 https://blog.youkuaiyun.com/zouxy09/article/details/17292011

2.通俗讲解SVM SVM---https://blog.youkuaiyun.com/v_JULY_v/article/details/7624837

3.SMO论文部分翻译解析 http://www.cnblogs.com/jerrylead/archive/2011/03/18/1988419.html

4.SMO论文分享：http://ishare.iask.sina.com.cn/f/13399644.html?from

5.机器学习实战[美]Peter Harrington.图灵程序设计丛书

6.数据集是源自github：https://github.com/pbharrin/machinelearninginaction 第六章-Ch06，当然上面也有源码

主要参考

线性可分支持向量机

最大间隔分类器：即最大化几何间隔

Lagrange dual function

线性可分支持向量机

假设数据是线性可分的，则存在一个支撑超平面 $f(X)=W^TX+b$ ,将数据分开。当 $f(X)=0$ 时，X便是位于超平面上的点。

而当 $f(X)>0$ 时，对应y=1的数据点， $f(X)<0$ ，对应于y=-1的数据点。y为数据类别标记

那么SVM做的就是找到一条直线，使得这条直线距离两边的数据间隔最大。

在超平面确定的情况下， $|W^TX+b|$ 能够表示数据点x到超平面距离的远近；另外通过观察 $W^TX+b$ 的符号与标记y是否一致判断分类是否正确，故用 $y*(W^TX+b)$ 的正负性判定或表示分类是否正确，这也称作函数间隔，即

$\hat{\gamma}=y(W^TX+b)=yf(X)$

而超平面（W,b）关于所有训练集样本点 $(x_i,y_i)$ 的函数间隔最小值便为超平面关于训练数据集的函数间隔，即

$\hat{\gamma}=min \hat{\gamma}_i,(i=1,...,n)$ （1）

但是，如果成比例地改变w和b，函数间隔也会等比例改变（超平面并没有改变）。所以需要对法向量W添加一些约束条件，从而引出真正定义点到超平面的距离----几何距离

，如图所示，由几何知识得 $x-x_0=\gamma\frac{w}{||w||}$ , $\frac{w}{||w||}$ 为单位向量

两边同乘以 $w^T$ ,得 $w^Tx-w^Tx_0=\gamma\frac{w^Tw}{||w||}$

又有 $w^Tx_0+b=0,w^Tw=||w||^2$ $\Rightarrow \gamma=\frac{w^Tx+b}{||w||}=\frac{f(x)}{||w||}$

同理，为了得到 $\gamma$ 的绝对值，令其乘以对应的类别y，即得到几何间隔

$\tilde{\gamma}=y \gamma=\frac{\hat{\gamma}}{||w||}$ (2)

从上述函数间隔和几何间隔的定义可以看出：几何间隔就是函数间隔除以||w||，而且函数间隔y*(wx+b) = y*f(x)实际上就是|f(x)|，只是人为定义的一个间隔度量，而几何间隔|f(x)|/||w||才是直观上的点到超平面的距离.

最大间隔分类器：即最大化几何间隔

目标函数： $max\tilde{\gamma}$ , 即 $arg_{w,b} max\left \{ \frac{1}{||w||}min_i[y_i(w^Tx_i+b)] \right \}$ ,前面(2)已经提到 $\hat\gamma=min_i[y_i(w^Tx_i+b)]$

（可以理解为求满足最小函数间隔的最大几何间隔）

约束条件： $y_i(w^Tx_i+b)=\hat{\gamma}_i\geq \hat{\gamma},i=1,...,n$

为了方便推导和优化，令 $\hat\gamma$ =1

则目标函数转化为 $max \frac{1}{||w||}$ ， $s.t. ,y_i(w^Tx_i+b)\geq 1,i=1,...,n$ (3)

那么，满足 $\mathbf{y(w^Tx+b)}=1$ ,即在间隔边界上的点是支持向量，＞1的则为不是支持向量的点。

(3)式可以转化为

$min\frac{1}{2}||w||^2$

$s.t. ,-y_i(w^Tx_i+b)+1\leq 0,i=1,...,n$

这里||w||的转换的做法，一方面是等价的，另一方面也是为了便于推导

这样就转化为凸优化（二次规划）问题

Lagrange dual function

$\begin{align*} g(\lambda,v) &= inf_{x \in D}L(x,\lambda,v)\\ &= inf_{x\in D}\left \{ f_0(x)+\sum_{i=1}^m\lambda_if_i(x)+\sum_{i=1}^pv_ih_i(x) \right \} \end{align*}$

若没有下确界，定义：

$g(\lambda,v)=-\infty$
根据定义，显然有：对∀λ>0，∀v，若原优化题有最优值p*，则

$g(\lambda,v)\leq p^*$
进一步：Lagrange对偶函数为凹函数。

可以使用Lagrange对偶函数求解

$\begin{align*} L(w,b,\alpha)&=\frac{1}{2}||w||^2+\sum_{i=1}^n\alpha_i(-y_i(w^Tx_i+b)+1) \\ &= \frac{1}{2}||w||^2-\sum_{i=1}^n\alpha_i(y_i(w^Tx_i+b)-1) \end{align*}$ （4）

然后令 $\Theta (w)=max_{\alpha_i\geq0}L(w,b,\alpha)$

容易验证，当某个约束条件不满足时，例如 $y_i(w^Tx_i+b)<1$ ，那么显然有 $\Theta (w)=\infty$ （只要令 $\alpha_i=\infty$ 即可）。而当所有约束条件都满足时，则最优值为 $\Theta (w)=\frac{1}{2}||w||^2$ ，亦即最初要最小化的量。

因此，在要求约束条件得到满足的情况下最小化 $\frac{1}{2}||w||^2$ ，实际上等价于直接最小化 $\Theta (w)$ （当然，这里也有约束条件，就是 $\alpha _i$ ≥0,i=1,…,n），因为如果约束条件没有得到满足， $\Theta (w)$ 会等于无穷大，自然不会是我们所要求的最小值。

故目标函数变成极小极大问题

$\Theta (w)=min_{w,b}max_{\alpha_i\geq0}L(w,b,\alpha)=p^*$

直接求解，会直接面对W,b两个参数，且 $\alpha _i$ 又是不等式约束，

故交换位置极大极小问题

$\Theta (w)=max_{\alpha_i\geq0}min_{w,b}L(w,b,\alpha)=d^*$

当然，要想使得 $d^*\leq p^*$ 需要满足KKT条件,

求解步骤：

① 固定 $\alpha$ ，关于w,b最小化L(对w,b求偏导)

$\frac{\partial L}{\partial w}=0\Rightarrow w=\sum_{i=1}^n\alpha_iy_ix_i$ （5）

$\frac{\partial L}{\partial b}=0\Rightarrow \sum_{i=1}^n\alpha_iy_i=0$ （6）

②回代L

$\begin{align*} L &=\frac{1}{2}w^Tw-\sum_{i=1}^n\alpha_i(y_i(w^Tx_i+b)-1)\\ &=\frac{1}{2}w^T(\sum_{i=1}^n\alpha_iy_ix_i) -\sum_{i=1}^n\alpha_iy_iw^Tx_i-\sum_{i=1}^n\alpha_iy_ib+\sum_{i=1}^n\alpha_i\\ &=-\frac{1} {2}w^T(\sum_{i=1}^n\alpha_iy_ix_i)-b\cdot 0+\sum_{i=1}^n\alpha_i\\ &= -\frac{1}{2}\sum_{i=1}^n\sum_{j=1}^n\alpha_i\alpha_jy_iy_jx_i^Tx_j +\sum_{i=1}^n\alpha_i \end{align*}$ （7）

---------------------------------------------------------------------------------------------------------------------------------------------------------------

非线性支持向量机（核函数）

附加的说明：正如PRML中线性回归章节所说，对于 $f(X)=W^TX+b$ 中的X，并非一定是线性映射，即

$f(X)=W^T\phi (X)+b$ （8）

$\mathbf{\phi }$ :表示X——>F是从原输入空间(低维)到某个特征空间(高维)的映射

将（5）式带入上式（8）得

$f(X)=\sum_{i=1}^n\alpha_iy_i<\phi (x_i),\phi(X)>+b$ ,其中 $<\phi (x_i),\phi(X)>$ 即 $\phi (x_i)^T\phi(X)$ 表示向量内积

但这存在一个问题：维度爆炸（直接在特征空间（高维空间）计算会导致维度呈指数级爆炸增长，会导致映射计算十分困难，甚至无法计算（维度过高））

因此，引出直接在原来的低维空间进行（内积）计算，而不需要显式地写出映射后的结果，即核函数。

比如， $x_1=\binom{\eta_1}{\eta_2},x_2=\binom{\xi_1}{\xi_2}$ , $\phi(\cdot )$ 为五维空间的映射，

用多项式核函数 $(<x_1,x_2>+1)^2=2\eta _1\xi _1+\eta_1^2\xi_1^2+2\eta_2\xi_2+\eta_2^2\xi_2^2+2\eta_1\eta_2\xi_1\xi_2+1$

代替 $\varphi (X_1,X_2)=(\sqrt{2}X_1,X_1^2,\sqrt{2}X_2,X_2^2,\sqrt{2}X_1X_2,1)^T$ 之后的内积 $<\varphi(x_1),\varphi(x_2)>$ 的结果。

而SVM里需要计算的数据向量正是以内积的形式出现（如公式(7)），因此针对非线性数据，SVM的处理方法是选择一个核函数 $k\mathbf{(x,x')}$ ，通过将数据映射到高维空间，来解决原始空间中线性不可分的问题。

常用的核函数有，多项式核： $k\mathbf{(x1,x2)}=(<\mathbf{x_1,x_2}>+R)^d$ ;

高斯核（径向基函数(radial basis function，rbf)）： $k(\mathbf{x1,x2})=exp\left ( -\frac{||\mathbf{x_1-x_2}||^2}{2\sigma ^2} \right )$

-----------------------------------------------------------------------------------------------------------------------------------------------------------------------

③此时L只包含一个变量 $\alpha_i$ ,然后对 $\alpha_i$ 求极大，便能求出w，b

即目标函数变为： $max\Psi (\alpha)=max_\alpha\sum_{i=1}^n\alpha_i -\frac{1}{2}\sum_{i,j=1}^n\alpha_i\alpha_jy_iy_jx_i^Tx_j$

求出 $\alpha_i$ 后，根据 $w=\sum_{i=1}^n\alpha_iy_ix_i$ ，即可求出w，然后通过 $b=-w^Tx$ 即可求出b

④求 $\alpha_i$ 的极大，主要利用SMO算法，如下所示，先将目标函数转化为求最小最优化问题，然后利用启发式算法选择其中两个Lagrange乘子 $\alpha_1,\alpha_2$ ,固定其余的乘子，逐步求解。

------------------------------------------------------------------------------------------------------------------------------------------------------------------------

线性支持向量机（松弛变量）

关于推导过程中参数C的说明：

前面的一系列推导都是假设数据是线性可分或者高维线性可分，但是对于有噪声点的数据（偏离正常位置很远的数据点，outliers），可能会导致超平面的位置变化（超平面本身就是由少数几个支持向量组成的，如果这些 support vector 里又存在 outlier 的话，其影响就很大了），如下图所示。

用黑圈圈起来的那个蓝点是一个 outlier，可以看出它对分隔超平面的影响

为了处理这种情况，SVM 允许数据点在一定程度上偏离一下超平面。将原来的约束条件公式（3）

$s.t. ,y_i(w^Tx_i+b)\geq 1,i=1,...,n$

改为

$s.t. ,y_i(w^Tx_i+b)\geq 1-\xi_i,i=1,...,n$ , $\xi_i\geq 0$ 称为松弛变量，

但是如果使得 $\xi_i$ 任意大的话，那么任意的超平面都是符合条件的，因此在原来的目标函数后面加上一个权重C，使得 $\xi_i$ (的总和)的影响不至于过大。[即是一个参数，用于控制目标函数中两项（“寻找 margin 最大的超平面”和“保证数据点偏差量最小”）之间的权重]

因此，原始目标函数变成

$min \frac{1}{2}||w||^2+C\sum_{i=1}^n\xi_i$ (9)

$s.t.,y_i(w^Tx_i+b)\geq1-\xi_i,i=1,...,n$

$\xi_i\geq0,i=1,...,n$

与之前相同，

求其Lagrange对偶函数:

$L(w,b,\xi,\alpha,r)=\frac{1}{2}||w||^2+C\sum_{i=1}^n\xi_i+\sum_{i=1}^n\alpha_i\left ( y_i(w^Tx_i+b) -1+\xi_i\right )-\sum_{i=1}^nr_I\xi_i$

分别对 $w,b,\xi_i$ 求偏导，

$\frac{\partial L}{\partial w}=0\Rightarrow w =\sum_{i=1}^n\alpha_iy_ix_i$

$\frac{\partial L}{\partial b}=0\Rightarrow \sum_{i=1}^n\alpha_iy_i=0$

$\frac{\partial L}{\partial \xi_i}=0\Rightarrow C-\alpha_i-r_i=0$

$i=1,...,n$

将 $w$ 回带到L并化简得到与公式（7）一样的目标函数

只是约束条件有所变化，

$C-\alpha_i-r_i=0\ and\ r_i\geq0\Rightarrow \alpha_i\leq C$

即对偶问题变为

$max_\alpha\sum_{i=1}^n\alpha_i-\frac{1}{2}\sum_{i,j=1}^n\alpha_i\alpha_jy_iy_j<x_i,x_j>$ （10）

$s.t.,0\leq \alpha_i \leq C,i=1,...,n$

$\sum_{i=1}^n\alpha_iy_i=0$

注意，C是已知的

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------

SMO推导过程及参数更新

那么 $\alpha_1^{new}=\alpha_1^{old}+s(\alpha_2^{old}-\alpha_2^{new,clipped})$

α1和α2的选取--启发式算法

坐标下降法通常每次迭代只调整一个参数的值，其他参数固定不变，但这里由于约束条件 $\sum_{i=1}^n\alpha_iy_i=0$ ,只改变一个α可能会导致约束条件失效，因此这里需要同时改变两个

对于，即第一个乘子，可以通过刚刚说的那3种不满足KKT的条件来找；
而对于第二个乘子可以寻找满足条件：的乘子。

第一个变量 $\alpha_i$ 的选择：

前面说，SVM的学习问题可以转化为对偶问题（公式(10)），需要满足的KKT条件：

$\alpha_i=0 \Leftrightarrow y_i(w^Tx_i+b)\geq1$

$0 < \alpha_i<C \Leftrightarrow y_i(w^Tx_i+b)= 1$

$\alpha_i=C \Leftrightarrow y_i(w^Tx_i+b)\leq 1$

这里的 $\alpha_i$ 还是拉格朗日乘子：

            1.对于第1种情况，表明 $\alpha_i$ 是正常分类，在间隔边界内部（我们知道正确分类的点）；
      2.对于第2种情况，表明 $\alpha_i$ 是支持向量，在间隔边界上；
            3.对于第3种情况，表明 $\alpha_i$ 是在两条间隔边界之间；

而最优解需要满足KKT条件，即上述3个条件都得满足，以下几种情况出现将会出现不满足：

            $y_i(w^Tx_i+b)\leq1$ 但是 $\alpha_i$ <C则是不满足的，而原本 $\alpha_i$ =C
            $y_i(w^Tx_i+b)\geq1$ 但是 $\alpha_i$ >0则是不满足的，而原本 $\alpha_i$ =0
            $y_i(w^Tx_i+b)=1$ 但是 $\alpha_i$ =0或者 $\alpha_i$ =C则是不满足的，而原本应该是 $0<\alpha_i<C$

SMO称选择第一个变量的过程为外层循环。外层训练在训练样本中选取违法KKT条件最严重的样本点（即上面的不满足的情况）。并将其对应的变量作为第一个变量。

该检验是在ε范围内进行的。在检验过程中，外层循环首先遍历所有满足条件 $0<\alpha_i<C$ 的样本点，即在间隔边界上的支持向量点，检验他们是否满足KKT条件，然后选择违反KKT条件最严重的 $\alpha_i$ 。如果这些样本点都满足KKT条件，那么遍历整个训练集，检验他们是否满足KKT条件，然后选择违反KKT条件最严重的 $\alpha_i$ 。

优先选择遍历非边界数据样本，因为非边界数据样本更有可能需要调整，边界数据样本常常不能得到进一步调整而留在边界上。由于大部分数据样本都很明显不可能是支持向量，因此对应的α乘子一旦取得零值就无需再调整。遍历非边界数据样本并选出他们当中违反KKT 条件为止。当某一次遍历发现没有非边界数据样本得到调整时，遍历所有数据样本，以检验是否整个集合都满足KKT条件。如果整个集合的检验中又有数据样本被进一步进化，则有必要再遍历非边界数据样本。这样，不停地在遍历所有数据样本和遍历非边界数据样本之间切换，直到整个样本集合都满足KKT条件为止。以上用KKT条件对数据样本所做的检验都以达到一定精度ε就可以停止为条件。如果要求十分精确的输出算法，则往往不能很快收敛。
对整个数据集的遍历扫描相当容易，而实现对非边界 $\alpha_i$ 的扫描时，首先需要将所有非边界样本的 $\alpha_i$ 值（也就是满足 $0<\alpha_i<C$ ）保存到新的一个列表中，然后再对其进行遍历。同时，该步骤跳过那些已知的不会改变的 $\alpha_i$ 值

第二个变量 $\alpha_j$ 的选择

计算阈值b：

优化 $\alpha_i$ 和 $\alpha_j$ 后，我们就可以更新阈值b，使得对两个样本i和j都满足KKT条件。如果优化后 $\alpha_i$ 不在边界上（也就是满足 $0<\alpha_i<C$ ，这时候根据KKT条件，可以得到 $y_iu_i=1,u_i=w^Tx_i+b=\sum_{j=1}^n\alpha_jy_j<x_i,x_j>+b$ ，这样我们才可以计算b），那下面的阈值b1是有效的，因为当输入 $x_i$ 时它迫使SVM输出 $y_i$

$b_1^{new}=b^{old}-E_1-y_1(\alpha_1^{new}-\alpha_1^{old})k(x_1,x_1)-y_2(\alpha_2^{new,clipped}-\alpha_2^{old})k(x_1,x_2)$

同样，如果 $0<\alpha_j<C$ ，那么下面的b2也是有效的：

$b_2^{new}=b^{old}-E_2-y_1(\alpha_1^{new}-\alpha_1^{old})k(x_1,x_2)-y_2(\alpha_2^{new,clipped}-\alpha_2^{old})k(x_2,x_2)$

如果 $0<\alpha_i<C$ 和 $0<\alpha_j<C$ 都满足，那么b1和b2都有效，而且他们是相等的。如果他们两个都处于边界上（也就是 $\alpha_i=0 \ \ or\ \ \alpha_i=C$ ，同时 $\alpha_j=0 \ \ or\ \ \alpha_j=C$ ），那么在b1和b2之间的阈值都满足KKT条件，一般我们取他们的平均值b=(b1+b2)/2。所以，总的来说对b的更新如下：

$b=\left\{\begin{matrix} b1& if \ \ 0<\alpha_i^{new}<C\\ b2& if \ \ 0<\alpha_j^{new,clipped}<C \\ (b_1+b_2)/2&otherwise \end{matrix}\right.$

每做完一次最小优化，必须更新每个数据样本的误差 $E_i$ ，以便用修正过的分类面对其他数据样本再做检验，在选择第二个配对优化数据样本时用来估计步长。

SMO算法的基本思路是：

如果说有变量的解都满足此最优化问题的KKT条件，那么这个最优化问题的解就得到了。因为KKT条件是该最优化问题的充分必要条件（证明请参考文献）。所以我们可以监视原问题的KKT条件，所以所有的样本都满足KKT条件，那么就表示迭代结束了。但是由于KKT条件本身是比较苛刻的，所以也需要设定一个容忍值，即所有样本在容忍值范围内满足KKT条件则认为训练可以结束

----------------------------------------------------------------------------------------------------------------------------------------------------------------

SMO算法的Python实现：

以下为个人根据机器学习实战这本书上的代码结合前面的推导以及理解写的代码（大部分相同，更多的是作为个人记录之用）

简版：

SMO_utils.py

import numpy as np

import pandas as pd
import matplotlib.pyplot as plt
import random
import pylab
from pandas import DataFrame, Series


plt.rcParams['font.sans-serif'] = ['SimHei']  #指定默认字体
plt.rcParams['axes.unicode_minus'] = False  #解决保存图像是负号'-'显示为方块的问题

def loadDataSet1(fileName):#读取数据
    df=pd.read_csv(fileName)
    data=df.ix[:,:-2]
    label=df.ix[:,-1]
    return data,label
def loadDataSet(fileName):
    dataMat=[];labelMat=[]
    f=open(fileName)
    for line in f.readlines():
        lineArr=line.strip().split('\t')
        dataMat.append([float(lineArr[0]),float(lineArr[1])])
        labelMat.append(float(lineArr[-1]))
    return dataMat,labelMat

def selectJrand(i,m):#随机选择一个不与αi重复的αj，
    j=i
    while j==i:
        j=int(random.uniform(0,m))
    return j
def clipAlphaJ(aj,H,L):#裁剪αj使得满足约束条件
    if aj>H:
        aj=H
    if aj<L:
        aj=L
    return aj

data,label=loadDataSet('datasets/testSet.txt')

simple_SMO.py

import numpy as np
import operator
import pandas as pd
import matplotlib.pyplot as plt
import pylab
from pandas import DataFrame, Series
from SMO_utils import selectJrand,loadDataSet,clipAlphaJ
plt.rcParams['font.sans-serif'] = ['SimHei']  #指定默认字体
plt.rcParams['axes.unicode_minus'] = False  #解决保存图像是负号'-'显示为方块的问题

class SimpleSMO(object):
    def __init__(self):
        '''
        该SMO的伪代码大致如下：
        创建一个alpha向量并将其初始化为0向量
        当迭代次数小于最大迭代次数时（外循环）：
            对数据集中的每个数据向量（内循环）：
                如果该数据向量（对应alpha_i）可以被优化：
                    随机选择另外一个数据向量（对应alpha_j）
                    同时优化这两个向量
                    如果两个向量都不能被优化，退出内循环
            如果所有向量都没被优化，增加迭代次数，继续下一次循环
        '''
        pass

    def smoSimple(self,data,labels,C,toler,maxIter):#toler容错率
        dataMatrix=np.mat(data);labelMat=np.mat(labels).transpose()#1维列向量
        b=0;m,n=np.shape(dataMatrix)
        alphas=np.mat(np.zeros((m,1)))#1维列向量
        iter=0
        while(iter<maxIter):#最大迭代次数
            alpha_pairs_changed=0#选择的alpha对是否已经进行优化
            for i in range(m):#更新αi 和αj（遍历整个alphas）
                f_Xi=np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)+b#u_i
                Ei=f_Xi-labelMat[i]#计算误差E_i=u_i-y_i
                if (labelMat[i]*Ei<-toler and alphas[i]<C) or (labelMat[i]*Ei>toler and alphas[i]>C) or (labelMat[i]*Ei==0 and (alphas[i]==C or alphas[i]==0)):#找到不满足KKT条件的αi，即y_i*u_i-1<0 and αi<C; y_i*u_i-1>0 and αi>C
                    j=selectJrand(i,m)#该简版直接随机选择j对应的αj
                    f_Xj=np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)+b#u_j
                    Ej=f_Xj-labelMat[j]
                    alpha_i_old=alphas[i].copy()
                    alpha_j_old=alphas[j].copy()
                    #将αj调整到(0,C)
                    if labelMat[i]!=labelMat[j]:#y_i,y_j异号
                        L=max(0,alphas[j]-alphas[i])
                        H=min(C,C+alphas[j]-alphas[i])
                    else:#同号
                        L=max(0,alphas[i]+alphas[j]-C)
                        H=min(C,alphas[j]+alphas[i])
                    if L==H:
                        print("L=H")
                        continue
                    eta=dataMatrix[i,:]*dataMatrix[i,:].T+dataMatrix[j,:]*dataMatrix[j,:].T-2.0*dataMatrix[i,:]*dataMatrix[j,:].T#η=K11+K22-2*K12
                    if eta<=0:#通常情况下目标函数是正定的，也就是说，能够在直线约束方向上求得最小值，并且η>0
                        print("eta<=0")
                        continue
                    alphas[j]+=labelMat[j]*(Ei-Ej)/eta#更新αj
                    alphas[j]=clipAlphaJ(alphas[j],H,L)#裁剪到(L,H)区间内
                    if abs((alphas[j]-alpha_j_old)<1e-7):#当更新变化不大时，停止更新--容忍值
                        print("j not moving enough")
                        continue
                    alphas[i]+=labelMat[j]*labelMat[i]*(alpha_j_old-alphas[j])#更新αi
                    b1=b-Ei-labelMat[i]*(alphas[i]-alpha_i_old)*dataMatrix[i,:]*dataMatrix[i,:].T-labelMat[j]*(alphas[j]-alpha_j_old)*dataMatrix[i,:]*dataMatrix[j,:].T#
                    b2=b-Ej-labelMat[i]*(alphas[i]-alpha_i_old)*dataMatrix[i,:]*dataMatrix[j,:].T-labelMat[j]*(alphas[j]-alpha_j_old)*dataMatrix[j,:]*dataMatrix[j,:].T
                    if alphas[i]>0 and alphas[i]<C:#更新满足条件后的b值
                        b=b1
                    elif alphas[j]>0 and alphas[j]<C:
                        b=b2
                    else:
                        b=(b1+b2)/2.
                    alpha_pairs_changed+=1
                    print("iter:%d,i:%d,pair changed:%d"%(iter,i,alpha_pairs_changed))
            if alpha_pairs_changed==0:
                iter+=1
            else:
                iter=0
            print("iteration number:%d"%iter)
        return b,alphas

smo_simple=SimpleSMO()
data,label=loadDataSet('datasets/testSet.txt')
b,alphas=smo_simple.smoSimple(data,label,0.6,0.01,100)
print(b)
idx=np.argwhere(alphas>0)
print(idx)
idx_dict={k:alphas[k,0] for k in idx[:,0]}
print(idx_dict)
sorted_idx_dict=sorted(idx_dict.items(),key=operator.itemgetter(1),reverse=True)
print(sorted_idx_dict,'-----')
fig=plt.figure()
ax=fig.add_subplot(1,1,1)
data=np.array(data)
sv_idx=sorted_idx_dict[:5]#取最大的几个α值作为支撑向量
x=[data[sv_idx[i][0],0] for i in range(len(sv_idx))]
y=[data[sv_idx[j][0],1] for j in range(len(sv_idx))]

for a,b in zip(x,y):
    # ax.annotate('*', xy=(a,b),verticalalignment='center',horizontalalignment='center')
    plt.text(a, b-0.1 , '*' , ha='center', va='center', fontsize=18,color='r')
ax.scatter(data[:,0],data[:,1])
plt.show()

plattSMO.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import operator
import random

plt.rcParams['font.sans-serif'] = ['SimHei']  #指定默认字体
plt.rcParams['axes.unicode_minus'] = False  #解决保存图像是负号'-'显示为方块的问题
'''
Platt SMO算法
①通过一个外循环来选择第一个alpha值，并且其选择过程会在两种方式之间进行交替：
    1、在所有数据集上进行单遍扫描；
    2、在非边界alpha中实现单遍扫描。[即(0,C)的alpha值]
对整个数据集扫描相当容易，而实现非边界alpha值得扫描时，首先需要建立这些alpha值的列表，然后再对这个表进行遍历。同时跳过已知的不会改变的alpha值
②选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha的值。在优化过程中，会通过最大化步长的方式来获得第二个alpha值。[ 建立一个全局的缓存用于保存误差值Ei，并从中选择max|Ei-Ej|,即最大化步长 ]
    
'''

class PlattSmoUtils(object):
    def __init__(self,data,label,C,toler,maxIter,kernel_tuple=('lin',0)):
        self.data=data#仅用于plot图像
        self.dataMatrix=np.mat(data)#数据矩阵化
        self.labelMatrix=np.mat(label).transpose()#一维列向量
        self.C=C
        self.toler=toler#容错率
        self.m=np.shape(np.mat(data))[0]#训练数据样本数量
        self.alphas=np.mat(np.zeros((self.m,1)))#初始化α为0矩阵，一维列向量
        self.b=0
        self.ECache=np.mat(np.zeros((self.m,2)))#全局变量，缓存误差
        self.maxIter=maxIter#不收敛时的最大迭代次数
        self.K=np.mat(np.zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i]=self.kernelTransfer(self.dataMatrix,self.dataMatrix[i,:],kernel_tuple)

    def kernelTransfer(self,X, A, kernel_tuple):
        m, n = X.shape
        K = np.mat(np.zeros((m, 1)))
        if kernel_tuple[0] == 'lin':#线性核函数
            K = X * A.T
        elif kernel_tuple[0] == 'rbf':#高斯核函数
            for j in range(m):
                delta_row = X[j, :] - A
                K[j] = delta_row * delta_row.T
            K = np.exp(K / (-1 * kernel_tuple[1] ** 2))
        else:
            raise NameError('Here we use only line kernel function and radial basis function')
        return K
    def cal_Ek(self,k):#计算误差Ei=u_i-y_i
        # f_Xk=np.multiply(self.alphas,self.labelMatrix).T*(self.dataMatrix*self.dataMatrix[k,:].T)+self.b#u_i
        f_Xk=np.multiply(self.alphas,self.labelMatrix).T*self.K[:,k]+self.b#u_i
        Ek=f_Xk-self.labelMatrix[k]
        return Ek

    def clipAlphaJ(self,aj, H, L):  # 裁剪αj使得满足约束条件
        if aj > H:
            aj = H
        if aj < L:
            aj = L
        return aj

    def selectJrand(self,i, m):  # 随机选择一个不与αi重复的αj，
        j = i
        while j == i:
            j = int(random.uniform(0, m))
        return j
    def selectJ(self,i,Ei):#根据max|Ei-Ej|选取αj
        max_k=-1
        max_delta_E=0#增量，保存Ei-Ej
        Ej=0
        self.ECache[i]=Ei
        valid_ECache_list=np.nonzero(self.ECache[:,0].A)[0]
        if len(valid_ECache_list)>1:
            for k in valid_ECache_list:
                if k==i:
                    continue
                Ek=self.cal_Ek(k)
                delta_E=abs(Ei-Ek)
                if delta_E>max_delta_E:
                    max_delta_E=delta_E
                    max_k=k
                    Ej=Ek
            return max_k,Ej
        else:
            j=self.selectJrand(i,self.m)
            Ej=self.cal_Ek(j)
        return j,Ej

    def updateEk(self,k):#通过更新后的α和b更新误差Ei
        Ek=self.cal_Ek(k)
        self.ECache[k]=[1,Ek]

    def innerLoop(self,i):
        Ei=self.cal_Ek(i)
        # print(Ei.shape)
        #or (self.labelMatrix[i] * Ei == 0 and (self.alphas[i] == self.C or self.alphas[i] == 0))#y_i*u_i-1=0 and αi=0或C
        if (self.labelMatrix[i]*Ei<-self.toler and self.alphas[i]<self.C) or (self.labelMatrix[i]*Ei>self.toler and self.alphas[i]>self.C):#找到不满足KKT条件的αi，即y_i*u_i-1<0 and αi<C; y_i*u_i-1>0 and αi>C;
            j,Ej=self.selectJ(i,Ei)
            alpha_i_old=self.alphas[i].copy()
            alpha_j_old=self.alphas[j].copy()
            if self.labelMatrix[i] != self.labelMatrix[j]:  # y_i,y_j异号
                L = max(0, self.alphas[j] - self.alphas[i])
                H = min(self.C, self.C + self.alphas[j] - self.alphas[i])
            else:  # 同号
                L = max(0, self.alphas[i] + self.alphas[j] - self.C)
                H = min(self.C, self.alphas[j] + self.alphas[i])
            if L == H:
                print("L=H")
                return 0
            # eta = self.dataMatrix[i, :] * self.dataMatrix[i, :].T + self.dataMatrix[j, :] * self.dataMatrix[j, :].T - 2.0 * self.dataMatrix[i,:] * self.dataMatrix[j,:].T  # η=K11+K22-2*K12
            eta=self.K[i,i]+self.K[j,j]-2.0*self.K[i,j]
            if eta <= 0:  # 通常情况下目标函数是正定的，也就是说，能够在直线约束方向上求得最小值，并且η>0
                print("eta<=0")
                return 0
            self.alphas[j] += self.labelMatrix[j] * (Ei - Ej) / eta  # 更新αj
            self.alphas[j] = self.clipAlphaJ(self.alphas[j], H, L)  # 裁剪到(L,H)区间内
            # print(self.alphas[j])
            self.updateEk(j)#更新Ej
            if abs((self.alphas[j] - alpha_j_old) < 1e-7):  # 当更新变化不大时，停止更新--容忍值
                print("j not moving enough")
                return 0
            self.alphas[i] += self.labelMatrix[j] * self.labelMatrix[i] * (alpha_j_old - self.alphas[j])  # 更新αi
            self.updateEk(i)#更新Ei
            # b1 = self.b - Ei - self.labelMatrix[i] * (self.alphas[i] - alpha_i_old) * self.dataMatrix[i, :] * self.dataMatrix[i, :].T - self.labelMatrix[j] * (self.alphas[j] - alpha_j_old) * self.dataMatrix[i, :] * self.dataMatrix[j, :].T
            b1=self.b-Ei-self.labelMatrix[i]*(self.alphas[i] - alpha_i_old)*self.K[i,i]-self.labelMatrix[j] * (self.alphas[j] - alpha_j_old) *self.K[i,j]
            # b2 = self.b - Ej - self.labelMatrix[i] * (self.alphas[i] - alpha_i_old) * self.dataMatrix[i, :] * self.dataMatrix[j, :].T - self.labelMatrix[j] * (self.alphas[j] - alpha_j_old) * self.dataMatrix[j, :] * self.dataMatrix[j, :].T  #
            b2 = self.b - Ej - self.labelMatrix[i] * (self.alphas[i] - alpha_i_old) * self.K[i,j] - self.labelMatrix[j] * (self.alphas[j] - alpha_j_old) * self.K[j,j]
            if self.alphas[i] > 0 and self.alphas[i] < self.C:  # 更新满足条件后的b值
                self.b = b1
            elif self.alphas[j] > 0 and self.alphas[j] < self.C:
                self.b = b2
            else:
                self.b = (b1 + b2) / 2.
            return 1
        else:
            return 0
    def start(self,sv_num=0):
        iter=0
        entire_set=True;alpha_pair_changed=0
        while(iter<self.maxIter) and (alpha_pair_changed>0 or entire_set):#外循环
            alpha_pair_changed=0#选择的alpha对是否已经进行优化
            if entire_set:#遍历全集
                for i in range(self.m):
                    alpha_pair_changed+=self.innerLoop(i)
                    print("fullSet--iter:%d i:%d, pairs changed %d"%(iter,i,alpha_pair_changed))
                iter+=1
            else:#遍历非边界α值
                non_bound_i=np.nonzero((self.alphas.A>0)*self.alphas.A<self.C)[0]#选取非边界αi
                for i in non_bound_i:
                    alpha_pair_changed+=self.innerLoop(i)
                    print("non bound,iter:%d i:%d,pairs changed %d"%(iter,i,alpha_pair_changed))
                iter+=1
            if entire_set:
                entire_set=False
            elif alpha_pair_changed==0:#跳过不再改变的值，迭代次数+1
                entire_set=True
            print("iteration number: %d"%iter)
        print(self.b)
        self.plot(self.data,self.alphas,sv_num)
        return self.b,self.alphas

    def plot(self,data, alphas,sv_num):#绘图
        idx = np.argwhere(alphas > 0)
        print(idx)
        idx_dict = {k: alphas[k, 0] for k in idx[:, 0]}
        # print(idx_dict)
        sorted_idx_dict = sorted(idx_dict.items(), key=operator.itemgetter(1), reverse=True)
        print(sorted_idx_dict, '-----','\n',len(sorted_idx_dict))

        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1)
        data = np.array(data)
        if sv_num==0:
            sv_idx=sorted_idx_dict[:]
        else:
            sv_idx = sorted_idx_dict[:sv_num]  # 取最大的几个α值作为支撑向量
        x = [data[sv_idx[i][0], 0] for i in range(len(sv_idx))]
        y = [data[sv_idx[j][0], 1] for j in range(len(sv_idx))]
        for a, b in zip(x, y):
            plt.text(a, b, '*', ha='center', va='center', fontsize=18, color='r')
        ax.scatter(data[:, 0], data[:, 1])
        plt.show()

def loadDataSet(fileName):
    dataMat=[];labelMat=[]
    f=open(fileName)
    for line in f.readlines():
        lineArr=line.strip().split('\t')
        dataMat.append([float(lineArr[0]),float(lineArr[1])])
        labelMat.append(float(lineArr[-1]))
    return dataMat,labelMat


def validate_line_kernel(t,sv_num=0):
    data, label = loadDataSet('datasets/testSet.txt')
    pSMO = PlattSmoUtils(data, label, 0.6, 0.001, 40)
    b, alphas = pSMO.start(sv_num)  # 绘图时绘出的支持向量的个数
    X=np.mat(data);labelMat=np.mat(label).transpose()
    m,n=X.shape
    w=np.zeros((n,1))
    for i in range(m):
        w+=np.multiply(alphas[i]*labelMat[i],X[i,:].T)

    f_X=X[t]*np.mat(w)+b
    print(f_X, label[t])
    if f_X*label[t]>0:
        return "True"
    else:
        return "False"
# bool_=validate_line_kernel(0,6)
# print(bool_)

def validate_rbf_kernel(sv_num=0,k1=0.1):#k1高斯核带宽
    data, label = loadDataSet('datasets/testSetRBF.txt')
    pSMO = PlattSmoUtils(data, label, 200, 0.0001, 1000,('rbf',k1))
    b, alphas = pSMO.start(sv_num)  # 绘图时绘出的支持向量的个数
    X = np.mat(data)
    labelMat = np.mat(label).transpose()
    sv_idx=np.nonzero(alphas.A>0)[0]
    # print(len(sv_idx))
    sv=X[sv_idx]#这几行代码效果等同于前面plot()函数中的代码[获取α>0的值对应的索引；然后根据索引得到对应的数据；数据标签]
    sv_label=labelMat[sv_idx]
    m, n = X.shape
    error_cnt=0
    for i in range(m):
        kernel_eval=pSMO.kernelTransfer(sv,X[i,:],kernel_tuple=('rbf',k1))#得到高斯核函数的内积映射
        predict=kernel_eval.T*np.multiply(sv_label,alphas[sv_idx])+b#计算f_x
        # print(np.sign(predict))
        if np.sign(predict)!=np.sign(label[i]):#决策函数sign()，判断准确性
            error_cnt+=1
    print("the training error rate is: %2.f"%(error_cnt/m))

    #测试集
    data, label = loadDataSet('datasets/testSetRBF2.txt')
    error_cnt=0
    X = np.mat(data);labelMat = np.mat(label).transpose()
    m, n = X.shape
    for i in range(m):
        kernel_eval = pSMO.kernelTransfer(sv, X[i, :], kernel_tuple=('rbf', k1))
        predict = kernel_eval.T * np.multiply(sv_label, alphas[sv_idx]) + b
        if np.sign(predict) != np.sign(label[i]):
            error_cnt += 1
    print("the test error rate is: %2.f" % (error_cnt / m))

validate_rbf_kernel(0,k1=1.3)

线性核函数对应图：validate_line_kernel()执行后的图（*为支持向量）：