“机器学习实战”刻意练习——分类问题：支持向量机

灰太狼 2023-06-17 15:54 0阅读 0赞

参考：

1.  [Python3《机器学习实战》学习笔记（八）：支持向量机原理篇之手撕线性SVM - Jack-Cui - CSDN博客][Python3_SVM - Jack-Cui - CSDN]
2.  [深入理解拉格朗日乘子法（Lagrange Multiplier) 和KKT条件 - Xianling Mao的专栏 - CSDN博客][Lagrange Multiplier_ _KKT_ - Xianling Mao_ - CSDN]
3.  [Python3《机器学习实战》学习笔记（九）：支持向量机实战篇之再撕非线性SVM - Jack-Cui - CSDN博客][Python3_SVM - Jack-Cui - CSDN 1]

# 一、SVM原理概述 #

支持向量机，即Support Vector Machines，简称SVM。

有些人认为，SVM是最好的现成的分类器，这里说的“现成”指的是分类器不加修改即可直接使用。同时，这就意味着在数据上应用基本形式的SVM分类器就可以得到低错误率的结果。SVM能够对训练集之外的数据点做出很好的分类决策。

SVM有很多实现，但是我们只关注其中最流行的一种实现，即**序列最小优化（Sequential Minimal Optimization，SMO）算法**。以及如何使用一种称为**核函数**（kernel）的方式将SVM扩展到更多数据集上。

*  优缺点：  
    优点：**泛化错误率低**，**计算开销不大**，结果易解释。  
    缺点：**对参数调节和核函数的选择敏感**，原始分类器不加修改仅适用于处理二类问题。
 *  适用数据类型：  
    **数值型**和**标称型数据**。
 *  SVM的一般流程  
    (1) 收集数据：可以使用任意方法。  
    (2) 准备数据：需要数值型数据。  
    (3) 分析数据：有助于**可视化分隔超平面**。  
    (4) 训练算法：SVM的大部分时间都源自训练，该过程主要实现**两个参数的调优**。  
    (5) 测试算法：十分简单的计算过程就可以实现。  
    (6) 使用算法：**几乎所有分类问题**都可以使用SVM，值得一提的是，SVM本身是一个**二类分类器**，对多类问题应用SVM需要对代码做一些修改。

## 1.基于最大间隔分隔数据 ##

在介绍SVM这个主题之前，先解释几个概念：

两组数据点，若它们之间已经分隔得足够开，很容易就可以在图中画出一条直线将两组数据点分开。  
在这种情况下，这组数据被称为**线性可分（linearly separable）数据**。  
而将数据集分隔开来的直线称为**分隔超平面**（separating hyperplane）。

若数据点都在二维平面上，则此时分隔超平面就只是一条直线。但是，如果所给的数据集是三维的，那么此时用来分隔数据的就是一个平面。显而易见，更高维的情况可以依此类推。如果数据集是1024维的，那么就需要一个1023维的某某对象来对数据进行分隔。这个1023维的某某对象到底应该叫什么？N-1维呢？  
该对象被称为**超平面**（hyperplane），也就是**分类的决策边界**。

分布在超平面一侧的所有数据都属于某个类别，而分布在另一侧的所有数据则属于另一个类别。我们希望能采用这种方式来构建分类器，**即如果数据点离决策边界越远，那么其最后的预测结果也就越可信**。

这里点到分隔面的距离被称为**间隔**（margin）。我们希望**间隔尽可能地大**，这是因为如果我们犯错或者在有限数据上训练分类器的话，我们希望分类器尽可能健壮。

**支持向量**（support vector）就是离分隔超平面最近的那些点。  
接下来要试着**最大化支持向量到分隔面的距离**，需要找到此问题的优化求解方法。

## 2.寻找最大间隔 ##

如何求解数据集的最佳分隔直线？

**分隔超平面的形式**可以写成  w T x + b w^\{T\}x+b wTx\+b

这里的常数 b b b类似于Logistic回归中的截距 w 0 w\_\{0\} w0

要计算点A到分隔超平面的距离，就必须给出点到分隔面的**法线或垂线的长度**，该值为  d = ∣ w T A + b ∣ ∣ ∣ w ∣ ∣ d=\\frac\{|w^\{T\}A+b|\}\{||w||\} d=∣∣w∣∣∣wTA\+b∣  
(可以类比初中学过的点到直线距离公式 d = ∣ A x + B y + C A 2 + B 2 ∣ d=\\left | \{\\frac\{Ax+By+C\}\{\\sqrt\{A^2+B^2\}\}\} \\right | d=∣∣∣A2\+B2Ax\+By\+C∣∣∣来理解)

接下来我们讨论分类器的工作原理：  
输入数据给分类器会**输出一个类别标**，这相当于一个类似于Sigmoid的函数在作用。  
下面将使用类似**海维赛德阶跃函数**（即单位阶跃函数）的函数对  w T x + b w^\{T\}x+b wTx\+b 作用得到  f ( w T x + b ) f(w^\{T\}x+b) f(wTx\+b)  
其中,当 u < 0 u<0 u<0时 f ( u ) f(u) f(u)输出-1，反之则输出+1。

**这里的类别标签采用-1和+1，而不是0和1主要有以下两个好处：**  
第一，由于-1和+1仅仅相差一个符号，方便数学上的处理。我们**可以通过一个统一公式来表示间隔或者数据点到分隔超平面的距离**，同时不必担心数据到底是属于-1还是+1类。  
第二，当计算数据点到分隔面的距离并确定分隔面的放置位置时，间隔通过 l a b e l ∗ ( w T x + b ) label \* (w^\{T\}x+b) label∗(wTx\+b)来计算。如果数据点处于正方向（即+1类）并且离分隔超平面很远的位置时， w T x + b w^\{T\}x+b wTx\+b会是一个很大的正数，同时 l a b e l ∗ ( w T x + b ) label \* (w^\{T\}x+b) label∗(wTx\+b)也会是一个**很大的正数**。而如果数据点处于负方向（-1类）并且离分隔超平面很远的位置时，此时由于类别标签为-1，则 l a b e l ∗ ( w T x + b ) label \* (w^\{T\}x+b) label∗(wTx\+b)仍然是一个**很大的正数**。

现在的目标就是找出分类器定义中的 w w w和 b b b。  
为此，我们必须找到具有最小间隔的数据点，而这些数据点也就是前面提到的支持向量。**一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化**。  
这就可以写作：  
![在这里插入图片描述][20191130183413868.png]  
直接求解上述问题相当困难，所以我们将它转换成为另一种更容易求解的形式如下：  
![在这里插入图片描述][20191130184438145.png]  
转换过程主要有两步（详细推导过程看上方参考文章，不赘述了）：  
一是令所有**支持向量**的 l a b e l ∗ ( w T x + b ) label \* (w^\{T\}x+b) label∗(wTx\+b)都为1，那么就可以通过求 ∣ ∣ w ∣ ∣ − 1 ||w||^\{-1\} ∣∣w∣∣−1的最大值来得到最终解。  
二是利用**拉格朗日乘子法及KKT条件**求解带约束条件的优化问题。  
其中，KKT条件是说最优值必须满足以下条件：

1.  L(a, b, x)对x求导为零；
2.  h(x) =0;
3.  a\*g(x) = 0;

转换后式子约束条件为：  
![在这里插入图片描述][20191130190245928.png]  
其中，常数C为引入的**松弛变量**（slack variable），来允许有些数据点可以处于分隔面的错误一侧。C用于控制“最大化间隔”和“保证大部分点的函数间隔小于1.0”这两个目标的权重。在优化算法的实现代码中，常数C是一个参数，因此我们就可以**通过调节该参数得到不同的结果**。

## 3.SMO 高效优化算法 ##

SMO表示序列最小优化（Sequential Minimal Optimization）。

SMO算法的目标是求出一系列  α \\alpha α 和  b b b ，一旦求出了这些  α \\alpha α ，就很容易计算出权重向量  w w w 并得到分隔超平面。

SMO算法的**工作原理**是：  
每次循环中选择**两个**  α \\alpha α 进行优化处理。  
一旦找到一对**合适**的  α \\alpha α ，那么就**增大其中一个同时减小另一个**。

*  这里所谓的“合适”就是指两个  α \\alpha α 必须要符合一定的条件  
    条件一是这两个  α \\alpha α 必须要在**间隔边界之外**  
    条件二是这两个  α \\alpha α **还没有进行过区间化处理或者不在边界上**

## 4.在复杂数据上应用核函数 ##

现在想象这样一个情形：A类数据点处于一个圆中，而B类数据点将将其环形包围。  
人类的大脑能够意识到这一点并找到其分割线。然而，对于分类器而言，它只能识别分类器的结果是大于0还是小于0。如果只在x和y轴构成的坐标系中插入直线进行分类的话，我们并不会得到理想的结果。

所以，我们可以将数据从一个特征空间转换到**另一个特征空间**。  
在新空间下，我们可以很容易利用已有的工具对数据进行处理，得到大于0或者小于0的测试结果。

数学家们喜欢将这个过程称之为从一个特征空间到另一个特征空间的**映射**。在通常情况下，这种映射会**将低维特征空间映射到高维空间**。  
经过空间转换之后，我们可以在高维空间中解决线性问题，这也就等价于在低维空间中解决非线性问题。

这种从某个特征空间到另一个特征空间的映射是通过**核函数**来实现的。  
读者可以把核函数想象成一个包装器（wrapper）或者是接口（interface），它能把数据从某个很难处理的形式转换成为另一个较容易处理的形式。

SVM优化中一个特别好的地方就是，所有的运算都可以写成内积（inner product，也称点积）的形式。向量的内积指的是两个向量相乘，之后得到单个标量或者数值。我们可以**把内积运算替换成核函数**，而不必做简化处理。  
将内积替换成核函数的方式被称为**核技巧（kernel trick**）或者核“变电”（kernel substation）。

核函数并不仅仅应用于支持向量机，很多其他的机器学习算法也都用到核函数。

**径向基函数**是SVM中常用的一个核函数。径向基函数是一个采用向量作为自变量的函数，能够基于向量距离运算输出一个标量。这个距离可以是从<0,0>向量或者其他向量开始计算的距离。  
我们经常使用到径向基函数的高斯版本，其具体公式为：  
 k ( x , y ) = e x p ( − ∥ x − y ∥ 2 2 σ 2 ) k(x,y)=exp(\\frac\{-\\left \\| x-y \\right \\|^2\}\{2\\sigma ^2\}) k(x,y)=exp(2σ2−∥x−y∥2)  
上述高斯核函数将数据从其特征空间映射到更高维的空间，具体来说这里是映射到一个**无穷维**的空间。

# 二、代码实现（python3） #

## 1.可视化数据集 ##

这一部分与之前的Logistic回归相同

import matplotlib.pyplot as plt
    import numpy as np
    
    def loadDataSet(fileName):
        """ 建立数据集 - - - - 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):
        """ 数据可视化 - - - - 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)              #转换为numpy矩阵
        data_minus_np = np.array(data_minus)            #转换为numpy矩阵
        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('5-SVM/testSet.txt')
        showDataSet(dataMat, labelMat)

结果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg_size_16_color_FFFFFF_t_70]  
可以看出使用的二维数据集显然线性可分。

## 2.简化SMO算法 ##

原Platt SMO算法的完整实现需要大量代码。我们将会对算法进行简化处理，以便了解算法的基本工作思路。当然，简化版代码虽然代码量少但会使执行速度变慢。

Platt SMO算法中的外循环确定要优化的最佳  α \\alpha α 对。而简化版却会跳过这一部分：首先在数据集上遍历每一个  α \\alpha α ，然后在剩下的  α \\alpha α 集合中**随机选择**另一个  α \\alpha α ，从而构建  α \\alpha α 对。

之所以这样做是因为我们有一个约束条件：  
![在这里插入图片描述][20191130202128204.png]  
由于改变一个  α \\alpha α 可能会导致该约束条件失效，因此我们总是同时改变两个  α \\alpha α 。

为此，我们将构建一个辅助函数，用于在某个区间范围内随机选择一个整数。同时，我们也需要另一个辅助函数，用于在数值太大时对其进行调整。

import random
    
    def selectJrand(i, m):
        """ 随机选择alpha - - - - i - alpha m - alpha参数个数 """
        j = i
        #随机选择一个j，使其与i不相等
        while (j == i):
            j = int(random.uniform(0, m))
        return j
    
    def clipAlpha(aj,H,L):
        """ 调整alpha值 - - - - aj - alpha值 H - alpha上限 L - alpha下限 """
        if aj > H:
            aj = H
        if L > aj:
            aj = L
        return aj
    
    def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
        """ 简化SMO算法 - - - - dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 maxIter - 最大迭代次数 """
        dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
        #初始化b参数，统计dataMatrix的维度
        b = 0; m,n = np.shape(dataMatrix)
        #初始化alpha参数，设为0
        alphas = np.mat(np.zeros((m,1)))
        #初始化迭代次数
        iter_num = 0
        #最多迭代matIter次
        while (iter_num < maxIter):
            #记录alpha是否已经进行优化
            alphaPairsChanged = 0
            for i in range(m):
                #步骤1：计算误差Ei，其中预测的分类f(Xi)=w^Tx[i]+b
                fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
                Ei = fXi - float(labelMat[i])
                #对满足优化条件的点，比较优化alpha
                #labelMat[i]*Ei=fXi*labelMat[i]-labelMat[i]*labelMat[i]
                # =fXi*yi-1
                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])
                    #保存更新前的alpha值，使用深拷贝
                    alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy()
                    #步骤2：判断是否同侧，计算用于调整alpha_j的上下界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：计算alpha_j最优修改量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 get_w(dataMat, labelMat, alphas):
        """ 计算直线法向量w - - - - dataMat - 数据矩阵 labelMat - 数据标签 alphas - 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()
    
    def showClassifer(dataMat, w, b):
        """ 分类结果可视化 - - - - 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)  
        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()
    
    if __name__ == '__main__':
        dataMat, labelMat = loadDataSet('5-SVM/testSet.txt')
        b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
        w = get_w(dataMat, labelMat, alphas)
        showClassifer(dataMat, w, b)

结果：

L==H
    L==H
    L==H
    第0次迭代 样本:3, alpha优化次数:1
    L==H
    第0次迭代 样本:5, alpha优化次数:2
    L==H
    L==H
    第0次迭代 样本:10, alpha优化次数:3
    第0次迭代 样本:13, alpha优化次数:4
    第0次迭代 样本:17, alpha优化次数:5
    alpha_j变化太小
    alpha_j变化太小
    alpha_j变化太小
    L==H
    alpha_j变化太小
    L==H
    L==H
    L==H
    alpha_j变化太小
    L==H
    alpha_j变化太小
    alpha_j变化太小
    alpha_j变化太小
    L==H
    L==H
    alpha_j变化太小
    alpha_j变化太小
    L==H
    alpha_j变化太小
    第0次迭代 样本:69, alpha优化次数:6
    ......
    alpha_j变化太小
    alpha_j变化太小
    alpha_j变化太小
    迭代次数: 40

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg_size_16_color_FFFFFF_t_70 1]  
虽然程序运行的时间比较长，但是可以看出分类结果还是很不错的

## 3.利用完整 Platt SMO 算法加速优化 ##

在几百个点组成的小规模数据集上，简化版SMO算法的运行是没有什么问题的，但是在更大的数据集上的运行速度就会变慢。

完整 Platt SMO 算法实现alpha的更改和代数运算的优化环节和简化版的一模一样。在优化过程中，唯一的不同就是**选择  α \\alpha α 的方式**。

完整版的Platt SMO算法应用了一些能够提速的启发方法。  
Platt SMO算法是通过一个**外循环**来选择第一个  α \\alpha α 值的，并且其选择过程会在两种方式之间进行交替：  
一种方式是在所有数据集上进行单遍扫描  
另一种方式则是在非边界  α \\alpha α 中实现单遍扫描。  
而所谓**非边界**  α \\alpha α 指的就是那些不等于边界0或C的  α \\alpha α 值。

在选择第一个  α \\alpha α 值后，算法会通过一个内循环来选择第二个  α \\alpha α 值。在优化过程中，会通过**最大化步长**的方式来获得第二个  α \\alpha α 值。

每个部分实现的功能如下：

*  **class optStruct:**  
    首要的事情就是建立一个**数据结构**来保存所有的重要值，而这个过程可以通过一个**对象**来完成。这样在将值传给函数使用时，可以省掉一些手工输入的麻烦。
 *  **calcEk(oS, k):**  
    由于**计算误差E值**在此版本SMO算法中**出现频繁**，这里必须要将其单独拎出来。
 *  **selectJ(i, oS, Ei)：**  
    选择合适的第二个alpha值以保证在每次优化中采用最大步长。  
    首先将输入值Ei在缓存中设置成为**有效**的。这里的有效（valid）意味着  
    它已经计算好了。接着程序会在所有的值上进行循环并选择其中使得改变最大的那个值。如果这是第一次循环的话，那么就随机选择一个alpha值。
 *  **updateEk()**  
    它会计算误差值并**存入缓存**当中。在对alpha值进行优化之后会用到这个值。
 *  **updateEk(oS, k):**  
    smoSimple()函数一模一样，但是这里的代码已经使用了自己的**数据结构**。另外，使用**selectJ()而不是selectJrand()来选择第二个alpha的值 。最后，在alpha值改变时更新Ecache**。
 *  **smoP(dataMatIn, classLabels, C, toler, maxIter):**  
    完整版的Platt SMO（外循环）算法。这里的一次迭代定义为一次循环过程，而不管该循环具体做了什么事。此时，如果在优化过程中**存在波动就会停止**，因此这里的做法优于smoSimple()函数中的计数方法。
 *  **calcWs(alphas,dataArr,classLabels):**  
    计算直线法向量w

class optStruct:
        """ 数据结构，维护所有需要操作的值 - - - - dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 """
        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)))  #根据矩阵行数初始化alpha参数为0 
            self.b = 0                                  #初始化b参数为0
            self.eCache = np.mat(np.zeros((self.m,2)))  #根据矩阵行数初始化误差缓存，第一列为是否有效的标志位，第二列为实际的误差E的值。
    
    def calcEk(oS, k):
        """ 计算误差 - - - - oS - 数据结构 k - 标号为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 selectJ(i, oS, Ei):
        """ 选择内循环alpha值 - - - - i - 标号为i的数据的索引值 oS - 数据结构 Ei - 标号为i的数据误差 """
        #初始化
        maxK = -1; maxDeltaE = 0; Ej = 0 
        #根据Ei更新误差缓存
        oS.eCache[i] = [1,Ei]      
        #返回误差不为0的数据的索引值
        validEcacheList = np.nonzero(oS.eCache[:,0].A)[0] 
        #存在不为0的误差 
        if (len(validEcacheList)) > 1:
            #遍历,找到相差最大的Ek 
            for k in validEcacheList:
                #跳过i
                if k == i: continue
                Ek = calcEk(oS, k)
                deltaE = abs(Ei - Ek)
                #找到maxDeltaE
                if (deltaE > maxDeltaE):
                    maxK = k; maxDeltaE = deltaE; Ej = Ek
            return maxK, Ej
        #没有不为0的误差
        else:
            #随机选择alpha_j的索引值
            j = selectJrand(i, oS.m)
            Ej = calcEk(oS, j)
        return j, Ej
    
    def updateEk(oS, k):
        """ 计算Ek,并更新误差缓存 - - - - oS - 数据结构 k - 标号为k的数据的索引值 """
        Ek = calcEk(oS, k)
        oS.eCache[k] = [1,Ek]
    
    def innerL(i, oS):
        """ 优化的SMO算法 - - - - i - 标号为i的数据的索引值 oS - 数据结构 """
        #步骤1：计算误差Ei
        Ei = calcEk(oS, i)
        #优化alpha,设定一定的容错率。
        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)):
            #使用内循环启发方式2选择alpha_j,并计算Ej
            j,Ej = selectJ(i, oS, Ei)
            #保存更新前的aplpha值，使用深拷贝
            alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
            #步骤2：计算上下界L和H
            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
            #步骤3：计算eta
            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
            #步骤4：更新alpha_j
            oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta
            #步骤5：修剪alpha_j
            oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
            #更新Ej至误差缓存
            updateEk(oS, j)
            if (abs(oS.alphas[j] - alphaJold) < 0.00001):
                print("alpha_j变化太小")
                return 0
            #步骤6：更新alpha_i
            oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
            #更新Ei至误差缓存
            updateEk(oS, i)
            #步骤7：更新b_1和b_2
            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
            #步骤8：根据b_1和b_2更新b
            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):
        """ 完整的线性SMO算法 - - - - dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 maxIter - 最大迭代次数 """
        #初始化
        oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler)
        iter = 0
        entireSet = True; alphaPairsChanged = 0
        #超过最大迭代次数或遍历整个数据集alpha也没有更新,则退出循环
        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:
                #遍历不在边界0和C的alpha
                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
            #如果alpha没有更新,计算全样本遍历
            elif (alphaPairsChanged == 0):
                entireSet = True 
            print("迭代次数: %d" % iter)
        return oS.b,oS.alphas
    
    def calcWs(alphas,dataArr,classLabels):
        """ 计算直线法向量w - - - - dataArr - 数据矩阵 classLabels - 数据标签 alphas - alphas值 """
        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__':
        #1.简化SMO
        """ dataMat, labelMat = loadDataSet('5-SVM/testSet.txt') b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40) w = get_w(dataMat, labelMat, alphas) """
        #2.完整SMO
        dataArr, classLabels = loadDataSet('5-SVM/testSet.txt')
        b, alphas = smoP(dataArr, classLabels, 0.6, 0.001, 40)
        w = calcWs(alphas,dataArr, classLabels)
        showClassifer(dataArr, classLabels, w, b)

结果：

全样本遍历:第0次迭代 样本:0, alpha优化次数:1
    全样本遍历:第0次迭代 样本:1, alpha优化次数:1
    全样本遍历:第0次迭代 样本:2, alpha优化次数:1
    全样本遍历:第0次迭代 样本:3, alpha优化次数:2
    全样本遍历:第0次迭代 样本:4, alpha优化次数:3
    全样本遍历:第0次迭代 样本:5, alpha优化次数:3
    ......
    全样本遍历:第2次迭代 样本:96, alpha优化次数:0
    全样本遍历:第2次迭代 样本:97, alpha优化次数:0
    全样本遍历:第2次迭代 样本:98, alpha优化次数:0
    全样本遍历:第2次迭代 样本:99, alpha优化次数:0
    迭代次数: 3

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg_size_16_color_FFFFFF_t_70 2]  
可以看出，完整版SMO算法选出的支持向量样点更多，更接近理想的分隔超平面。

## 4.非线性SVM ##

使用数据集分布如下图：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg_size_16_color_FFFFFF_t_70 3]

class optStruct2:
        """ 数据结构，维护所有需要操作的值 - - - - dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 kTup - 包含核函数信息的元组,第一个参数存放核函数类别，第二个参数存放必要的核函数需要用到的参数 """
        def __init__(self, dataMatIn, classLabels, C, toler, kTup):
            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)))  #根据矩阵行数初始化alpha参数为0 
            self.b = 0                                  #初始化b参数为0
            self.eCache = np.mat(np.zeros((self.m,2)))  #根据矩阵行数初始化误差缓存，第一列为是否有效的标志位，第二列为实际的误差E的值。
            self.K = np.mat(np.zeros((self.m,self.m)))  #初始化核K
            for i in range(self.m):                     #计算所有数据的核K
                self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
    
    def kernelTrans(X, A, kTup):
        """ 通过核函数将数据转换更高维的空间 - - - - X - 数据矩阵 A - 单个数据的向量 kTup - 包含核函数信息的元组 """
        m,n = np.shape(X)
        K = np.mat(np.zeros((m,1)))
        #线性核函数,只进行内积。
        if kTup[0] == 'lin': K = X * A.T 
        #高斯核函数,根据高斯核函数公式进行计算
        elif kTup[0] == 'rbf':
            for j in range(m):
                deltaRow = X[j,:] - A
                K[j] = deltaRow*deltaRow.T
            #计算高斯核K
            K = np.exp(K/(-1*kTup[1]**2))
        else: raise NameError('核函数无法识别')
        return K
    
    def calcEk2(oS, k):
        """ 计算误差 - - - - oS - 数据结构 k - 标号为k的数据 """
        fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
        Ek = fXk - float(oS.labelMat[k])
        return Ek
    
    def selectJ2(i, oS, Ei):
        """ 内循环启发方式2 Parameters： i - 标号为i的数据的索引值 oS - 数据结构 Ei - 标号为i的数据误差 Returns: j, maxK - 标号为j或maxK的数据的索引值 Ej - 标号为j的数据误差 """
        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 = calcEk2(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 = calcEk2(oS, j) 
        return j, Ej 
    
    def updateEk2(oS, k):
        """ 计算Ek,并更新误差缓存 - - - - oS - 数据结构 k - 标号为k的数据的索引值 """
        Ek = calcEk2(oS, k)
        oS.eCache[k] = [1,Ek]
    
    def innerL2(i, oS):
        """ 优化的SMO算法 - - - - i - 标号为i的数据的索引值 oS - 数据结构 """
        #步骤1：计算误差Ei
        Ei = calcEk2(oS, i)
        #优化alpha,设定一定的容错率。
        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)):
            #使用内循环启发方式2选择alpha_j,并计算Ej
            j,Ej = selectJ2(i, oS, Ei)
            #保存更新前的aplpha值，使用深拷贝
            alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy()
            #步骤2：计算上下界L和H
            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
            #步骤3：计算eta
            eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]
            if eta >= 0:
                print("eta>=0")
                return 0
            #步骤4：更新alpha_j
            oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta
            #步骤5：修剪alpha_j
            oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
            #更新Ej至误差缓存
            updateEk2(oS, j)
            if (abs(oS.alphas[j] - alphaJold) < 0.00001):
                print("alpha_j变化太小")
                return 0
            #步骤6：更新alpha_i
            oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
            #更新Ei至误差缓存
            updateEk2(oS, i)
            #步骤7：更新b_1和b_2
            b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]
            b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]
            #步骤8：根据b_1和b_2更新b
            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 smoP2(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin',0)):
        """ 完整的线性SMO算法 - - - - dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 maxIter - 最大迭代次数 kTup - 包含核函数信息的元组 """
        oS = optStruct2(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler, kTup) 
        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 += innerL2(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 += innerL2(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 testRbf(k1 = 1.3):
        """ 非线形Rbf核函数测试函数 - - - - k1 - 使用高斯核函数的时候表示到达率 """
        dataArr,labelArr = loadDataSet('5-SVM/testSetRBF.txt')
        #根据训练集计算b和alphas
        b,alphas = smoP2(dataArr, labelArr, 200, 0.0001, 100, ('rbf', k1))
        datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()
        svInd = np.nonzero(alphas.A > 0)[0]
        sVs = datMat[svInd]                                                     
        labelSV = labelMat[svInd]
        print("支持向量个数:%d" % np.shape(sVs)[0])
        m,n = np.shape(datMat)
        errorCount = 0
        for i in range(m):
            #计算各个点的核
            kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))             
            #根据支持向量的点，计算超平面，返回预测结果 
            predict = kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b     
            #返回数组中各元素的正负符号，用1和-1表示，并统计错误个数
            if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1      
        #打印错误率 
        print("训练集错误率: %.2f%%" % ((float(errorCount)/m)*100))
        dataArr,labelArr = loadDataSet('5-SVM/testSetRBF2.txt')
        errorCount = 0
        datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()         
        m,n = np.shape(datMat)
        for i in range(m):
            kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))      
            predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b 
            if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1 
        print("测试集错误率: %.2f%%" % ((float(errorCount)/m)*100)) 
    
    if __name__ == '__main__':
    	#3.非线性SVM
        testRbf()

结果：

全样本遍历:第0次迭代 样本:0, alpha优化次数:1
    全样本遍历:第0次迭代 样本:1, alpha优化次数:1
    全样本遍历:第0次迭代 样本:2, alpha优化次数:1
    全样本遍历:第0次迭代 样本:3, alpha优化次数:2
    全样本遍历:第0次迭代 样本:4, alpha优化次数:2
    全样本遍历:第0次迭代 样本:5, alpha优化次数:3
    ......
    全样本遍历:第3次迭代 样本:98, alpha优化次数:0
    全样本遍历:第3次迭代 样本:99, alpha优化次数:0
    迭代次数: 4
    支持向量个数:21
    训练集错误率: 2.00%
    测试集错误率: 4.00%

# 三、小结 #

支持向量机是一种分类器。之所以称为“机”是因为它会产生一个二值决策结果，即它是一种**决策“机”**。

支持向量机的**泛化错误率较低**，也就是说它具有**良好的学习能力**，且学到的结果具有**很好的推广性**。这些优点使得支持向量机十分流行，有些人认为它是监督学习中最好的定式算法。

支持向量机试图**通过求解一个二次优化问题来最大化分类间隔**。  
John Platt引入了**SMO算法**，此算法可以通过每次**只优化2个alpha值**来加快SVM的训练速度。

**核方法**或者说核技巧会将数据（有时是非线性数据）从一个低维空间**映射**到一个高维空间，可以将一个在低维空间中的非线性问题转换成高维空间下的线性问题来求解。  
核方法不止在SVM中适用，还可以用于其他算法中。  
而其中的**径向基函数**是一个常用的度量两个向量距离的核函数。

支持向量机是一个二类分类器。当用其解决多类问题时，则需要额外的方法对其进行扩展。SVM的效果也**对优化参数和所用核函数中的参数敏感**。

[Python3_SVM - Jack-Cui - CSDN]: https://blog.csdn.net/c406495762/article/details/78072313
[Lagrange Multiplier_ _KKT_ - Xianling Mao_ - CSDN]: https://blog.csdn.net/xianlingmao/article/details/7919597
[Python3_SVM - Jack-Cui - CSDN 1]: https://blog.csdn.net/c406495762/article/details/78158354
[20191130183413868.png]: https://img-blog.csdnimg.cn/20191130183413868.png
[20191130184438145.png]: https://img-blog.csdnimg.cn/20191130184438145.png
[20191130190245928.png]: https://img-blog.csdnimg.cn/20191130190245928.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/20191130210038487.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg==,size_16,color_FFFFFF,t_70
[20191130202128204.png]: https://img-blog.csdnimg.cn/20191130202128204.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg_size_16_color_FFFFFF_t_70 1]: https://img-blog.csdnimg.cn/20191201115040941.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg==,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg_size_16_color_FFFFFF_t_70 2]: https://img-blog.csdnimg.cn/20191202162726720.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg==,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg_size_16_color_FFFFFF_t_70 3]: https://img-blog.csdnimg.cn/20191202164432210.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L25hbmFzaGlfRg==,size_16,color_FFFFFF,t_70