SVM支持向量机

SVM简介

前几次的文章中提到过，机器学习的主要任务分为分类问题和回归问题；而SVM是一个典型的解决分类问题的算法。现在假定如下图的一个分类任务的情况

如图a所示，在二维空间内分布着一些共有两种类别的分类样本，对于分类问题而言，最基本的解决思路，就是找到一个划分的超平面，将不同的样本分开；其中图b中的线A和图c中的线B都是这样的超平面。

但正如上图所示，事实上这样的超平面可以找到无数个；那么我们应该以哪个作为最终选定的超平面呢？从直观上来讲，我们会认为图b中的线A是一个更加优秀的选择，因为对于这个超平面而言，其具备更强的鲁棒性；也就是实际应用中，由于各种原因产生的样本扰动并不会对图b产生太大的影响，换言之，实际使用的对象（也包括测试集）中会出现大量的样本更逼近这个超平面。但A的超平面由于宏观上与不同样本分类间的“间隔”最大，所以应对这个问题的鲁棒性也就更强，对未见示例的泛化性能也更出色。

总而言之，SVM的算法就是试图找出那个泛化性能更强的划分超平面，依据则是那些距离超平面最近的支持向量与超平面的间隔。(后面会对支持向量的概念进行补充)

支持向量与间隔

那么在这里就可以给出SVM中划分超平面的定义式

其中需要说明的是：

• w={w1;w2;w3…wd}为法向量，决定超平面的方向
• b为位移项，据欸的那个了超平面与坐标系原点的距离
• x为训练样本
• 上面的例子是二维，但实际应用中维度可以更高，所以是”超平面“

而如上图所示，样本空间中任意点到超平面的距离都可以用上面的式子表示；这里用到的是类似与中学的点到直线距离的公式，推导过程简单，在此不再赘述。

这里继续推进，假设超平面已经能将样本正确分类，且yi的取值为{-1，+1}（这里要注意，默认在超平面上方的样本为+1，反之为-1，这里仅能取两个值，而不是一个范围）
那么存在

也可以简化为如下方程

显然，是距离平面最近的几个训练样本使得上式成立，它们也被称为支持向量；而两个不同类的支持向量到超平面的距离之和被称为间隔。

综上，为了找到具有最大间隔的划分超平面，也就是要找到满足上式中的约束参数w和b使得γ最大，即

当然，为了最大化间隔的方便，这里可以将最大化问题等价转化为最小化问题，上式可以重写为下

这里需要补充，虽然式子上看起来间隔只与w有关，但事实上b通过约束隐式地影响着w的取值，进而对间隔产生影响

范例——简单求解超平面

上面的数学推导过程可能有些简略，一时间也许有些难以理解，那么我们在这里引入一个简单的范例，给定三个样本，找出一个超平面将其分隔。

将上面的式子代入上面提到的优化式如下

可以得到一个如下的方程组

3w1+3w2+b≥1
4w1+3w2+b≥1
-w1-w2-b≥1

联立求解后得到最终解应该是

对偶问题

这一部分会涉及到一部分数学知识，不甚了解的可以先查看下面的这个博客进行补足
https://blog.youkuaiyun.com/pipisorry/article/details/52135854

事实上，实际求解问题原比上面的例子复杂的多，因为实际的问题的维度往往很高，且一般的带约束优化问题无法直接通过梯度下降等迭代方法求解，所以这里引入一个高效的数学工具：拉格朗日乘子法

具体来说，首先对前面提到的优化式添加拉格朗日乘子αi≥0；则可写为

接着，令L（w，b，a） 对w和b的偏导为零可得

两式联立，消去w和b，并且将约束条件带入；可得到优化式的对偶问题

解得α后，再继续求出w和b即可得到模型如下

到了这一步我们可以注意到上面式子得结果需要满足KKT条件，也就是要求

SMO算法

为了求解上面的对偶问题，我们可以简单的使用二次规划算法进行求解；但显而易见的，正比于训练样本规模的二次规划算法会在实际任务中产生大量的算法开销；因此我们根据问题特性，选择一个更有效的方法SMO（要注意这只是方法之一，还有其他的方法有同样效果。）

SMO的基本思路实现固定αi以外的所有参数，然后求αi上的极值，由于存在约束；固定了αi后其可由其他变量导出；于是，SMO每次选择两个变量αi和αj且固定其他参数。这样，在参数初始化后，SMO不断执行如下两个步骤直至收敛

• 选取一对需要更新的变量αi和αj
• 固定其他参数，求解更新后的αi和αj

这样，在仅考虑αi和αj时，对偶问题的约束就变为

上面的式子反代入前面的对偶问题的式子后可以得到关于αi的单变量二次规划问题，且具有闭式解。这时就可以不必调用数值优化算法即可高性地计算更新后的αi和αj。

而偏移项b的确定则由支持向量完成，也即是使用所有支持向量求解的平均值。

这里需要注意

当然，本部分的描叙比较简略，适合对高数和SMO算法有一定了解的读者，如果不甚理解同样推荐阅读下面这篇博客
https://blog.youkuaiyun.com/qq_39521554/article/details/80723770

核函数

前面的问题解决后，任然存在一个问题：我们默认训练样本是线性可分的。可实际任务中，存在大量线性不可分的样本空间，就如同下面的例子

这里引入一个被证实的数学原理：假使原始空间是有限维，那么一定存在一个更高维的空间可以使得原本不可分的样本空间可分。

以上面的图为示意，当样本组在左边的二维空间中线性不可分时，以某种映射规则投射到三维空间后，它就变得线性可分了。而中间的那个映射关系的内积就是我们这节内容的核函数。映射函数本身仅仅是一种映射关系，并没有增加维度的特性，不过可以利用核函数的特性，构造可以增加维度的核函数，这通常是我们希望的。

要注意，核函数和映射没有关系。核函数只是用来计算映射到高维空间之后的内积的一种简便方法。

一般英文文献对Kernel有两种提法，一是Kernel Function，二是Kernel Trick。从Trick一词中就可以看出，这只是一种运算技巧而已，不涉及什么高深莫测的东西。

具体巧在哪里呢？我们如果想进行原本就线性不可分的数据集进行分割，那么选项一是容忍错误分类，即引入Soft Margin；选项二是我们可以对Input Space做Feature Expansion，把数据集映射到高维中去，形成了Feature Space。我们几乎可以认为（引用Caltech的课堂用语“We are safe but not certain”）原本在低维中线性不可分的数据集在足够高的维度中存在线性可分的超平面。

常用核函数

• Linear Kernel——线性核是最简单的核函数，核函数的数学公式如下：
  

• Polynomial Kernel——多项式核实一种非标准核函数，它非常适合于正交归一化后的数据
  

• Gaussian Kernel——高斯核函数，鲁棒径向基核对于数据中的噪音有着较好的抗干扰能力，其参数决定了函数作用范围，超过了这个范围，数据的作用就“基本消失”。
  

• Exponential Kernel——指数核函数就是高斯核函数的变种，它仅仅是将向量之间的L2距离调整为L1距离，这样改动会对参数的依赖性降低，但是适用范围相对狭窄。
  

• Laplacian Kernel——拉普拉斯核完全等价于指数核，唯一的区别在于前者对参数的敏感性降低，也是一种径向基核函数。
  

• ANOVA Kernel——ANOVA 核也属于径向基核函数一族，其适用于多维回归问题
  

• Sigmoid Kernel——Sigmoid 核来源于神经网络，具体作用在前几篇中已有论述
  

软间隔与正则化

在前面的讨论中，我们一直假定训练样本在样本空间中是线性可分的，存在一个超平面能够把不同类的样本完全划分开。然而，在现实任务中往往是很难找到合适的核函数将训练样本“完全划分”。即使找到了某个核函数刚好将训练集在特征空间中线性可分，但是这种情况也有可能是过拟合造成的。

对于上述情况我们必须引入“软间隔”概念。左图要求所有的样本必须划分正确，这称之为“硬间隔”。而“软间隔”允许SVM在一些样本上出错，即允许某些样本不满足约束:

我们原始的问题是

在引入软间隔理论后，我们的优化目标就变成了下面的形式：

这里 C>0为 惩罚参数，可以理解为我们一般回归和分类问题正则化时候的参数， C越大，对于误分类的惩罚越大，反之越小。也就是说我们希望目标函数尽量小，误分类的点尽可能的少。C是协调两者关系的正则化惩罚系数。在实际应用中，需要调参来选择。这个目标函数的优化和SVM线性可分的优化方式类似，下面我们来看看如何对线性分类SVM的软间隔最大化进行学习。

和线性可分SVM的优化方式类似，首先我们将软间隔最大化的约束问题用拉格朗日函数转化为无约束问题，如下：

也就是说，我们现在优化的目标函数是：

这个优化目标也满足KKT条件，也就是说，我们可以通过拉格朗日对偶将我们的优化问题转化为等价的对偶问题来求解，如下：

剩下的部分和上面的SVM部分一样，最终 我们得到的优化目标如下：

总的来说，软间隔相当于正则化的过程，可以有效地防止训练过程中模型过拟合，从而更好地进行分类工作。

简易的SMO算法demo

需要注意，本部分使用的代码来自人名邮电出版社《机器学习实战》一书的源码，备注为笔者添加

首先查看一下源码中自带数据集的样子

可以看到，是一个二维的数据集，图中只截取了一部分；接下来使用书中的代码进行测试

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)              #转换为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('testSet.txt')
    showDataSet(dataMat, labelMat)

上面的代码很简单，就是直接将样本集读入并显示，结果如下

接下来采用书上的代码进行测试；这部分就是对smo算法的实现

def selectJrand(i, m):
    j = i                                 #选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return j


def clipAlpha(aj,H,L):  #调整alpha值
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    
    dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose() #转换为numpy的mat存储
    b = 0; m,n = np.shape(dataMatrix)   #初始化b并统计dataMatrix的维度
    alphas = np.mat(np.zeros((m,1)))    #初始化alpha
    iter_num = 0
    while (iter_num < maxIter):
        alphaPairsChanged = 0
        for i in range(m):      #计算误差
            fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):   #优化alpha
                j = selectJrand(i,m)
                fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
                if (labelMat[i] != labelMat[j]):    #计算上下界L和H
                    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
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T     #计算eta
                if eta >= 0: print("eta>=0"); continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta  #更新alpha_j并调整
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小"); continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#更新alpha_i
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T#更新b_1和b_2
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
                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 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 get_w(dataMat, labelMat, alphas):  #更新W
    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()

可以看到，确实漂亮的完成了任务。并且也在图中标注了支持向量

实例——手写体数字的分类

在本学期的第一次实验中，我们曾经完成过手写体数字的分类；但当时采用的是KNN算法，这里简单的判断一下两者的不同和优缺点。

knn没有训练过程，他的基本原理就是找到训练数据集里面离需要预测的样本点距离最近的k个值（距离可以使用比如欧式距离，k的值需要自己调参），然后把这k个点的label做个投票，选出一个label做为预测。对于KNN，没有训练过程。只是将训练数据与训练数据进行距离度量来实现分类。

svm需要超平面wx+b来分割数据集（此处以线性可分为例），因此会有一个模型训练过程来找到w和b的值。训练完成之后就可以拿去预测了，根据函数y=wx+b的值来确定样本点x的label，不需要再考虑训练集。对于SVM，是先在训练集上训练一个模型，然后用这个模型直接对测试集进行分类。

而对应的knn没有训练过程，但是预测过程需要挨个计算每个训练样本和测试样本的距离，当训练集和测试集很大时，预测效率感人。svm有一个训练过程，训练完直接得到超平面函数，根据超平面函数直接判定预测点的label，预测效率很高（一般我们更关心预测效率）。

https://blog.youkuaiyun.com/Keiji582/article/details/120515323
上面贴的是笔者的KNN算法实验博客，有一部分的本节代码是从这里修改而来，所以不会进行过多的介绍

再实验前再次回顾一下本次实验的数据集格式

接下来，在上面的代码的基础上，补充一些针对数字分类这个任务的函数

def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin',0)):

    oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler, kTup)#初始化
    iter = 0#初始化当前迭代次数
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):#遍历整个数据集都alpha也没有更新或者超过最大迭代次数,则退出循环
        alphaPairsChanged = 0
        if entireSet:#遍历整个数据集
            for i in range(oS.m):
                alphaPairsChanged += innerL(i,oS)#使用优化的SMO算法
            iter += 1
        else:#遍历非边界值
            nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]         
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
            iter += 1
        if entireSet: #改为非边界遍历
            entireSet = False
        elif (alphaPairsChanged == 0):#如果alpha没有更新,计算全样本遍历
            entireSet = True
    return oS.b,oS.alphas#返回SMO算法计算的b和alphas

def img2vector(filename):
    returnVect = np.zeros((1, 1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0, 32 * i + j] = int(lineStr[j])
    return returnVect
def loadImages(dirName):
    from os import listdir
    hwLabels = []
    trainingFileList = listdir(dirName)  # load the training set
    m = len(trainingFileList)
    trainingMat = np.zeros((m, 1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]  # take off .txt
        classNumStr = int(fileStr.split('_')[0])
        if classNumStr == 9:
            hwLabels.append(-1)
        else:
            hwLabels.append(1)
        trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))
    return trainingMat, hwLabels

def testDigits(kTup=('rbf', 10)):	#核函数超参数设置
    dataArr, labelArr = loadImages('trainingDigits')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    datMat = np.mat(dataArr);
    labelMat = np.mat(labelArr).transpose()
    svInd = np.nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]
    labelSV = labelMat[svInd];
    print("there are %d Support Vectors" % np.shape(sVs)[0])
    m, n = np.shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
        if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1
    print("the training error rate is: %f" % (float(errorCount) / m))
    dataArr, labelArr = loadImages('testDigits')
    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, :], kTup)
        predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
        if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1
    print("the test error rate is: %f" % (float(errorCount) / m))

上面的函数有部分因为是改写自KNN算法中的函数，所以没有特别注明！
接下来看看实际运行的效果

但我们试着修改def testDigits(kTup=('rbf', 10)): #核函数超参数设置中超参数的设置后看看结果如何

当参数为（Linear）时

而当我们把α的值设置为50时

这里我们可以看到一个明显的现象：SVM本身的性能和表现极其受制于参数调节和核函数的选择。这也是其明显的缺点和不足。

总结

需要补充的是，上面的数字分类实际上是一个二分类任务，样本只包含1和9两种类别。

但我们要怎么解决多分类问题呢？比如我们想要将上面的代码修改为区分0-9的数字要怎么做呢？

一个简单且行之有效的思路是：对每一个待分类类别训练一个SVM将其与其他的类别区分开，只需要将上面的代码小小修改就可以做到。（因笔者临近期末时间较为紧张，会在后续有时间时补上）

SVM的优点

• SVM 是一种有坚实理论基础的新颖的小样本学习方法。它基本上不涉及概率测度及大数定律等,因此不同于现有的统计方法。从本质上看,它避开了从归纳到演绎的传统过程,实现了高效的从训练样本到预报样本的转导推理,大大简化了通常的分类和回归等问题。
• SVM 的最终决策函数只由少数的支持向量所确定,计算的复杂性取决于支持向量的数目,而不是样本空间的维数,这在某种意义上避免了维数灾难。
• 少数支持向量决定了最终结果,这不但可以帮助我们抓住关键样本、“剔除”大量冗余样本,而且注定了该方法不但算法简单,而且具有较好的鲁棒性。

SVM的缺点

• SVM算法对大规模训练样本难以实施
• SVM本身的性能和表现极其受制于参数调节和核函数的选择

当然，还有一个不算缺点的缺点，就是支持向量机的数学推导过程不再像前几个实验一样浅显易懂，受制于笔者捉襟见肘的表达能力，对于数学部分不是很理解的读者可以移步他处寻找更加完整的推导过程；本篇博客仅针对于已经了解部分原理和过程的读者。