Support Vector Machine（一）

最新推荐文章于 2022-12-18 22:14:45 发布

苹果pie

最新推荐文章于 2022-12-18 22:14:45 发布

阅读量212

点赞数

分类专栏：机器学习

机器学习专栏收录该内容

13 篇文章 0 订阅

订阅专栏

参考

1) 崔家华 https://cuijiahua.com/blog/2017/11/ml_8_svm_1.html

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

SVM的英文全称是Support Vector Machines，我们叫它支持向量机。支持向量机是我们用于分类的一种算法。让我们以一个小故事的形式，开启我们的SVM之旅吧。

在很久以前的情人节，一位大侠要去救他的爱人，但天空中的魔鬼和他玩了一个游戏。

魔鬼在桌子上似乎有规律放了两种颜色的球，说："你用一根棍分开它们？要求：尽量在放更多球之后，仍然适用。"

于是大侠这样放，干的不错？

然后魔鬼，又在桌上放了更多的球，似乎有一个球站错了阵营。显然，大侠需要对棍做出调整。

SVM就是试图把棍放在最佳位置，好让在棍的两边有尽可能大的间隙。这个间隙就是球到棍的距离。

现在好了，即使魔鬼放了更多的球，棍仍然是一个好的分界线。

魔鬼看到大侠已经学会了一个trick(方法、招式)，于是魔鬼给了大侠一个新的挑战。

现在，大侠没有棍可以很好帮他分开两种球了，现在怎么办呢？当然像所有武侠片中一样大侠桌子一拍，球飞到空中。然后，凭借大侠的轻功，大侠抓起一张纸，插到了两种球的中间。

æºå¨å¦ä¹ å®ææç¨ï¼å«ï¼ï¼æ¯æåéæºåçç¯ä¹ææçº¿æ§SVM

现在，从空中的魔鬼的角度看这些球，这些球看起来像是被一条曲线分开了。

再之后，无聊的大人们，把这些球叫做data，把棍子叫做classifier, 找到最大间隙的trick叫做optimization，拍桌子叫做kernelling, 那张纸叫做hyperplane。

概述一下：

当一个分类问题，数据是线性可分的，也就是用一根棍就可以将两种小球分开的时候，我们只要将棍的位置放在让小球距离棍的距离最大化的位置即可，寻找这个最大间隔的过程，就叫做最优化。但是，现实往往是很残酷的，一般的数据是线性不可分的，也就是找不到一个棍将两种小球很好的分类。这个时候，我们就需要像大侠一样，将小球拍起，用一张纸代替小棍将小球进行分类。想要让数据飞起，我们需要的东西就是核函数(kernel)，用于切分小球的纸，就是超平面。

线性SVM

æºå¨å¦ä¹ å®ææç¨ï¼å«ï¼ï¼æ¯æåéæºåçç¯ä¹ææçº¿æ§SVM

上图中的(a)是已有的数据，红色和蓝色分别代表两个不同的类别。数据显然是线性可分的，但是将两类数据点分开的直线显然不止一条。上图的(b)和(c)分别给出了B、C两种不同的分类方案，其中黑色实线为分界线，术语称为“决策面”。每个决策面对应了一个线性分类器。虽然从分类结果上看，分类器A和分类器B的效果是相同的。但是他们的性能是有差距的，看下图：

在"决策面"不变的情况下，我又添加了一个红点。可以看到，分类器B依然能很好的分类结果，而分类器C则出现了分类错误。显然分类器B的"决策面"放置的位置优于分类器C的"决策面"放置的位置，SVM算法也是这么认为的，它的依据就是分类器B的分类间隔比分类器C的分类间隔大。

显然每一个可能把数据集正确分开的方向都有一个最优决策面（有些方向无论如何移动决策面的位置也不可能将两类样本完全分开），而不同方向的最优决策面的分类间隔通常是不同的，那个具有“最大间隔”的决策面就是SVM要寻找的最优解。而这个真正的最优解对应的两侧虚线所穿过的样本点，就是SVM中的支持样本点，称为"支持向量"。

数学建模

一个最优化问题通常有两个最基本的因素：

1）目标函数，也就是你希望什么东西的什么指标达到最好；

2）优化对象，你期望通过改变哪些因素来使你的目标函数达到最优。

在线性SVM算法中，目标函数显然就是那个“分类间隔”，而优化对象则是决策面。

决策面

简单的来看二维空间里：

æºå¨å¦ä¹ å®ææç¨ï¼å«ï¼ï¼æ¯æåéæºåçç¯ä¹ææçº¿æ§SVM

直线A、B可表示为： y=ax+b (2.1)

现在我们做个小小的改变，让原来的轴变成 x_1 轴，变成 x_2 轴，于是公式(2.1)中的直线方程会变成下面的样子

x_2=ax_1+b (2.2)

ax_1+(-1)x_2+b=0 (2.3)

公式（2.3）的向量形式可以写成

$[a,-1]\left[ \begin{array}{c}x_1\\x_2\end{array} \right] +b=0$ (2.4)

令 $\boldsymbol{\omega}=[\omega_1,\omega_2]^T$ ， $\boldsymbol{x}=[x_1,x_2]^T$ 则2.4式可表示为 $\boldsymbol{\omega}^T\boldsymbol{x}+\gamma=0$ (2.5)（一般说W，X均指列向量）

æºå¨å¦ä¹ å®ææç¨ï¼å«ï¼ï¼æ¯æåéæºåçç¯ä¹ææçº¿æ§SVM

向量 $\boldsymbol{\omega}=[\omega_1,\omega_2]^T$ 跟直线 $\boldsymbol{\omega}^T\boldsymbol{x}+\gamma=0$ 是相互垂直的，也就是说 $\boldsymbol{\omega}=[\omega_1,\omega_2]^T$ 控制了直线的方向。 $\gamma$ 就是截距，它控制了直线的位置。

分类间隔

间隔的大小实际上就是支持向量对应的样本点到决策面的距离的二倍

æºå¨å¦ä¹ å®ææç¨ï¼å«ï¼ï¼æ¯æåéæºåçç¯ä¹ææçº¿æ§SVM

$d = \frac{|\boldsymbol{\omega}^T\boldsymbol{x}+\gamma|}{||\boldsymbol{\omega}||}$ (2.6) , $||\boldsymbol{\omega}||$ 是向量 $\boldsymbol{\omega}$ 的模，表示在空间中向量的长度，。追求的最大化也就是寻找的最大化。

约束条件

（1）并不是所有的方向都存在能够实现100%正确分类的决策面，我们如何判断一条直线是否能够将所有的样本点都正确分类？

对于上图的例子，为每个样本点 $\boldsymbol{x}_i$ 加上一个类别标签 y_i ：

如果我们的决策面方程能够完全正确地对图2中的样本点进行分类，就会满足下面的公式（参见逻辑回归）

$\left\{\begin{array}{ll} \boldsymbol{\omega}^T\boldsymbol{x}_i+\gamma>0 & \textrm{for~~} y_i=1\\\boldsymbol{\omega}^T\boldsymbol{x}_i+\gamma<0 & \textrm{for~~} y_i=-1\end{array}\right.$ (2.8)

?????????????????????????????????????????

$\left\{\begin{array}{ll} (\boldsymbol{\omega}^T\boldsymbol{x}_i+\gamma)/||\boldsymbol{\omega}||\geq d & \forall~ y_i=1\\(\boldsymbol{\omega}^T\boldsymbol{x}_i+\gamma)/||\boldsymbol{\omega}||\leq -d & \forall~y_i=-1\end{array}\right.$ （2.9）

左右两边除上d，就可得到：

$\left\{\begin{array}{ll} \boldsymbol{\omega}_d^T\boldsymbol{x}_i+\gamma_d\geq 1 & \textrm{for~~} y_i=1\\\boldsymbol{\omega}_d^T\boldsymbol{x}_i+\gamma_d\leq -1 & \textrm{for~~} y_i=-1\end{array}\right.$ (2.10)

其中 $\boldsymbol{\omega}_d = \frac{\boldsymbol{\omega}}{||\boldsymbol{\omega}||d},~~ \gamma_d = \frac{\gamma}{||\boldsymbol{\omega}||d}$

对wd和γd重新起个名字，就叫它们w和γ。因此，就可以说："对于存在分类间隔的两类样本点，我们一定可以找到一些超平面，使其对于所有的样本点均满足下面的条件："

(2.11)

(2.12) 距离超平面距离最近的点称为支持向量，去除除了支持向量的点，超平面不变

$d = \frac{|\boldsymbol{\omega}^T\boldsymbol{x}_i+\gamma|}{||\boldsymbol{\omega}||}=\frac{1}{||\boldsymbol{\omega}||},~~\textrm{if} ~\boldsymbol{x}_i \textrm {is a support vector}$ (2.13)

可能有人会问了，为什么标记为1和-1呢？因为这样标记方便我们将上述方程变成如下形式：

(2.14)

正是因为标签为1和-1，才方便我们将约束条件变成一个约束方程，从而方便我们的计算。

原来的任务是找到一组参数 $\boldsymbol{\omega},~\gamma$ 使得分类间隔 W=2d 最大化，根据公式(2.13)就可以转变为 $||\boldsymbol{\omega}||$ 的最小化问题，也等效于 $\frac{1}{2}||\boldsymbol{\omega}||^2$ 的最小化问题。我们之所以要在 $||\boldsymbol{\omega}||$ 上加上平方和1/2的系数，是为了以后进行最优化的过程中对目标函数求导时比较方便，但这绝不影响最优化问题最后的解。

将最终的目标函数和约束条件放在一起进行描述：

æºå¨å¦ä¹ å®ææç¨ï¼å«ï¼ï¼æ¯æåéæºåçç¯ä¹ææçº¿æ§SVM

。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。

接着干

我们已经得到支持向量机的基本数学模型，接下来的问题就是如何根据数学模型，求得我们想要的最优解。

1）有约束优化问题的几何意象：闭上眼睛你看到什么？

2）拉格朗日乘子法：约束条件怎么跑到目标函数里面去了？

3）KKT条件：约束条件是不等式该怎么办？

4）拉格朗日对偶：最小化问题怎么变成了最大化问题？

5）实例演示：拉格朗日对偶函数到底啥样子？

6）SVM优化算法的实现：数学讲了辣么多，到底要怎么用啊？

首先

æºå¨å¦ä¹ å®ææç¨ï¼å«ï¼ï¼æ¯æåéæºåçç¯ä¹ææçº¿æ§SVM

一个非凸函数，我们无法获得全局最优解的，只能获得局部最优解。对于我们的目标函数：它是一个凸函数。

对于Question 1 ：

无约束优化问题，可以写为：

尝试使用的方法就是费马大定理(Fermat)，即使用求取函数f(x)的导数，然后令其为零，可以求得候选最优值，再在这些候选值中验证；如果是凸函数，可以保证是最优解。

有等式约束的优化问题，可以写为：

常常使用的方法就是拉格朗日乘子法（Lagrange Multiplier) ，即把等式约束h_i(x)用一个系数与f(x)写为一个式子，称为拉格朗日函数，而系数称为拉格朗日乘子。通过拉格朗日函数对各个变量求导，令其为零，可以求得候选值集合，然后验证求得最优值。

有不等式约束的优化问题，可以写为：

常常使用的方法就是KKT条件。同样地，我们把所有的等式、不等式约束与f(x)写为一个式子，也叫拉格朗日函数，系数也称拉格朗日乘子，通过一些条件，可以求出最优值的必要条件，这个条件称为KKT条件。

对于Question2 ：

将有约束的原始目标函数转换为无约束的新构造的拉格朗日目标函数

æºå¨å¦ä¹ å®ææç¨ï¼å«ï¼ï¼æ¯æåéæºåçç¯ä¹ææçº¿æ§SVM 其中αi是拉格朗日乘子，αi大于等于0，是构造新目标函数时引入的系数变量(自己设置)。

令

当样本点不满足约束条件时，即在可行解区域外：

此时，将αi设置为正无穷，此时θ(w)显然也是正无穷。

当样本点满足约束条件时，即在可行解区域内：

此时，显然θ(w)为原目标函数本身。将上述两种情况结合一下，就得到了新的目标函数：

此时，就建立了一个在可行解区域内与原目标函数相同，在可行解区域外函数值趋近于无穷大的新函数。

现在，问题变成了求新目标函数的最小值，即：

这里用p*表示这个问题的最优值，且和最初的问题是等价的。

先求最大值，再求最小值。这样的话，首先就要面对带有需要求解的参数w和b的方程，而αi又是不等式约束，这个求解过程不好做。所以，需要使用拉格朗日函数对偶性，将最小和最大的位置交换一下，这样就变成了：

交换以后的新问题是原始问题的对偶问题，这个新问题的最优值用d*来表示。而且d*<=p*。

Question3、4 ：

需要什么条件才能让d=p呢？

首先必须满足这个优化问题是凸优化问题。
其次，需要满足KKT条件。

凸优化问题和KKT都满足了，问题转换成了对偶问题。而求解这个对偶学习问题，可以分为三个步骤：首先要让L(w,b,α)关于w和b最小化，然后求对α的极大，最后利用SMO算法求解对偶问题中的拉格朗日乘子。

第一步

根据上述推导已知：

首先固定α，要让L(w,b,α)关于w和b最小化，我们分别对w和b偏导数，令其等于0，即：

将上述结果带回L(w,b,α)得到：

从上面的最后一个式子，我们可以看出，此时的L(w,b,α)函数只含有一个变量，即αi。

第二步

现在内侧的最小值求解完成，我们求解外侧的最大值，从上面的式子得到

通过smo优化算法能得到α，再根据α，我们就可以求解出w和b，进而求得我们最初的目的：找到超平面，即"决策平面"。

将目标函数变形，在前面增加一个符号，将最大值问题转换成最小值问题：

通过引入所谓的松弛变量ξ(slack variable)和惩罚参数C，来允许有些数据点可以处于超平面的错误的一侧。此时约束条件有所改变：

同时，考虑到松弛变量和松弛变量ξ和惩罚参数C，目标函数变为：

原始问题的拉格朗日函数变为：

对偶问题拉格朗日函数的极大极小问题，得到以下等价优化问题：

???????????????????????????????????

此处省略很多，实在看不懂了。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。

。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。

SMO算法

1996年，John Platt发布了一个称为SMO的强大算法，用于训练SVM。Platt的SMO算法是将大优化问题分解为多个小优化问题来求解的。这些小优化问题往往很容易求解，并且对它们进行顺序求解的结果与将它们作为整体来求解的结果完全一致的。在结果完全相同的同时，SMO算法的求解时间短很多。

SMO算法的工作原理是：每次循环中选择两个alpha进行优化处理。一旦找到了一对合适的alpha，那么就增大其中一个同时减小另一个。这里所谓的"合适"就是指两个alpha必须符合以下两个条件，条件之一就是两个alpha必须要在间隔边界之外，而且第二个条件则是这两个alpha还没有进行过区间化处理或者不在边界上。

SMO算法的步骤：

步骤1：计算误差：

步骤2：计算上下界L和H：

步骤3：计算η：

步骤4：更新αj：

步骤5：根据取值范围修剪αj：

步骤6：更新αi：

步骤7：更新b1和b2：

步骤8：根据b1和b2更新b：

动手吧

1、数据集下载

2、查看数据集

import matplotlib.pyplot as plt
import numpy as np


def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():  # 逐行读取，滤除空格等
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])  # 添加数据
        labelMat.append(float(lineArr[2]))  # 添加标签
    return dataMat, labelMat


def showDataSet(dataMat, labelMat):
    data_plus = []
    data_minus = []
    for i in range(len(dataMat)):
        if labelMat[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = np.array(data_plus)
    data_minus_np = np.array(data_minus)
    plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1])
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1])
    plt.show()


if __name__ == '__main__':
    dataMat, labelMat = loadDataSet('testSet.txt')
    showDataSet(dataMat, labelMat)

3、简单的smo算法

from time import sleep
import matplotlib.pyplot as plt
import numpy as np
import random
import types


def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():  # 逐行读取，滤除空格等
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])  # 添加数据
        labelMat.append(float(lineArr[2]))  # 添加标签
    return dataMat, labelMat


# 在0到m的范围内随机选一个与i不同的j
def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j


def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    # 转换为numpy的mat存储
    dataMatrix = np.mat(dataMatIn);
    labelMat = np.mat(classLabels).transpose()
    # 初始化b参数，统计dataMatrix的维度
    b = 0
    m, n = np.shape(dataMatrix)  # m行n列
    # 初始化alpha参数，设为0
    alphas = np.mat(np.zeros((m, 1)))  # 生成m行1列元素为0的矩阵
    # 初始化迭代次数
    iter_num = 0
    # 最多迭代matIter次
    while (iter_num < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            # 步骤1：计算误差Ei
            # **********************************************************************
            fXi = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
            Ei = fXi - float(labelMat[i])
            # 优化alpha，更设定一定的容错率。
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
                # 随机选择另一个与alpha_i成对优化的alpha_j
                j = selectJrand(i, m)
                # 步骤1：计算误差Ej
                fXj = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])
                # 保存更新前的aplpha值，使用深拷贝
                alphaIold = alphas[i].copy();
                alphaJold = alphas[j].copy();
                # 步骤2：计算上下界L和H
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L == H: print("L==H"); continue
                # 步骤3：计算eta
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - \
                      dataMatrix[i, :] * dataMatrix[i, :].T - \
                      dataMatrix[j, :] * dataMatrix[j, :].T
                if eta >= 0: print("eta>=0"); continue
                # 步骤4：更新alpha_j
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                # 步骤5：修剪alpha_j
                alphas[j] = clipAlpha(alphas[j], H, L)
                if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小"); continue
                # 步骤6：更新alpha_i
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])
                # 步骤7：更新b_1和b_2
                b1 = b - Ei - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[i, :].T - \
                     labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[i, :] * dataMatrix[j, :].T
                b2 = b - Ej - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[j, :].T - \
                     labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[j, :] * dataMatrix[j, :].T
                # 步骤8：根据b_1和b_2更新b
                if (0 < alphas[i]) and (C > alphas[i]):
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                else:
                    b = (b1 + b2) / 2.0
                # 统计优化次数
                alphaPairsChanged += 1
                # 打印统计信息
                print("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num, i, alphaPairsChanged))
        # 更新迭代次数
        if (alphaPairsChanged == 0):
            iter_num += 1
        else:
            iter_num = 0
        print("迭代次数: %d" % iter_num)
    return b, alphas


def showClassifer(dataMat, w, b):
    # 绘制样本点
    data_plus = []  # 正样本
    data_minus = []  # 负样本
    for i in range(len(dataMat)):
        if labelMat[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = np.array(data_plus)  # 转换为numpy矩阵
    data_minus_np = np.array(data_minus)  # 转换为numpy矩阵
    plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)  # 正样本散点图
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7)  # 负样本散点图
    # 绘制直线
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b - a1 * x1) / a2, (-b - a1 * x2) / a2
    plt.plot([x1, x2], [y1, y2])
    # 找出支持向量点
    for i, alpha in enumerate(alphas):
        if alpha > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()


def get_w(dataMat, labelMat, alphas):
    alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)
    w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)
    return w.tolist()


if __name__ == '__main__':
    dataMat, labelMat = loadDataSet('testSet.txt')
    b, alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
    w = get_w(dataMat, labelMat, alphas)
    showClassifer(dataMat, w, b)

smo函数太长了。。。

完整版SMO

优化：

最外层循环，首先在样本中选择违反KKT条件的一个乘子作为最外层循环，然后用"启发式选择"选择另外一个乘子并进行这两个乘子的优化
在非边界乘子中寻找使得 |Ei - Ej| 最大的样本
如果没有找到，则从整个样本中随机选择一个样本

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


class optStruct:

    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = np.shape(dataMatIn)[0]
        self.alphas = np.mat(np.zeros((self.m, 1)))
        self.b = 0
        self.eCache = np.mat(np.zeros((self.m, 2)))


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


def calcEk(oS, k):
    fXk = float(np.multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k, :].T) + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek


def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j


def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1, Ei]
    validEcacheList = np.nonzero(oS.eCache[:, 0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i: continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej


def updateEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]


def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or (
            (oS.labelMat[i] * Ei > oS.tol) and (oS.alphas[i] > 0)):
        j, Ej = selectJ(i, oS, Ei)
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L == H:
            print('L==H')
            return 0
        eta = 2.0 * oS.X[i, :] * oS.X[j, :].T - oS.X[i, :] * oS.X[i, :].T - oS.X[j, :] * oS.X[j, :].T
        if eta >= 0:
            print('eta>=0')
            return 0
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print('aplha_j变化太小')
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        updateEk(oS, i)
        b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[i, :].T - \
             oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.X[i, :] * oS.X[j, :].T
        b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[j, :].T - \
             oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.X[j, :] * oS.X[j, :].T
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
            oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
            oS.b = b2
        else:
            oS.b = (b1 + b2) / 2.0
        return 1
    else:
        return 0


def smoP(dataMatIn, classLabels, C, toler, maxIter):
    oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:
            for i in range(oS.m):
                alphaPairsChanged += innerL(i, oS)
                print('全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d' % (iter, i, alphaPairsChanged))
            iter += 1
        else:
            nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i, oS)
                print('非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d' % (iter, i, alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False
        elif (alphaPairsChanged == 0):
            entireSet = True
        print('迭代次数: %d' % iter)
    return oS.b, oS.alphas


def showClassifer(dataMat, classLabels, w, b):
    data_plus = []
    data_minus = []
    for i in range(len(dataMat)):
        if classLabels[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = np.array(data_plus)
    data_minus_np = np.array(data_minus)
    plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7)
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b - a1 * x1) / a2, (-b - a1 * x2) / a2
    plt.plot([x1, x2], [y1, y2])
    for i, alpha in enumerate(alphas):
        if abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()


def calcWs(alphas, dataArr, classLabels):
    X = np.mat(dataArr)
    labelMat = np.mat(classLabels).transpose()
    m, n = np.shape(X)
    w = np.zeros((n, 1))
    for i in range(m):
        w += np.multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w


if __name__ == '__main__':
    dataArr, classLabels = loadDataSet('testSet.txt')
    b, alphas = smoP(dataArr, classLabels, 0.6, 0.001, 40)
    w = calcWs(alphas, dataArr, classLabels)
    showClassifer(dataArr, classLabels, w, b)