【机器学习实战-6章】SVM支持向量机

最新推荐文章于 2022-01-01 21:18:44 发布

唐生一

最新推荐文章于 2022-01-01 21:18:44 发布

阅读量637

点赞数

分类专栏：机器学习文章标签： python 机器学习算法 svm

本文链接：https://blog.csdn.net/lynn_flame/article/details/106602329

版权

机器学习专栏收录该内容

18 篇文章 1 订阅

订阅专栏

*****SVM有很多实现，本章仅关注其中最流行的一种实现：序列最小优化(Sequential Minimal Optimization,SMO)算法，之后，介绍如何使用核函数(kernel)的方式将SVM扩展到更多数据集上

6.1 基于最大间隔分隔数据

SVM的优缺点

优点：泛化错误率低，计算开销不大，结果易解释
缺点：对参数调节和核函数的选择敏感，原始分类器不加修改仅使用于处理二类问题
使用数据类型：数值型和标称型数据
分隔超平面(separating hyperplane)
如果一组数据集可以用一条直线分成两类，我们称该直线为分隔超平面，更高维的数据则由高维的对象进行分割，称为超平面(hyperplane),也就是分类的决策边界。
我们希望采用这种方式来构建分类器，即如果数据点离决策边界越远，那么其最后的预测结果也就越可信。我们希望找到离分隔超平面最近的点，确保它们离分割面的距离尽可能远。点到分割面的距离被称为间隔(margin).我们希望间隔尽可能地大，这样的分类器会更加健壮。
支持向量(support vector)
支持向量就是离分隔超平面最近的那些点。接下来要试着最大化支持向量到分隔面的距离，需要找到此问题的优化求解方法。

6.2 寻找最大间隔

分隔超平面的形式可以写成:
W.TX+b
计算点A到分隔超平面的距离，即点到分割面的法线或垂线的长度，即|W.TA+b|/|W|——点到直线距离公式
在这里插入图片描述

6.2.1 分类器求解的优化问题

分类器的工作原理：输入数据给分类器会输出一个类别标签，这相当于一个类似于Sigmoid函数在作用
本节使用单位阶跃函数对W.TX+b作用得到f(W.TX+b),其中，当u<0时，f(u)输出-1,反之输出+1
采用-1,+1更利于数学上的处理。当计算数据点到分隔面的距离并确定分隔面的放置位置时，间隔通过label*(W.TX+b)来计算。label为+1，-1,正好将数据分布在分隔面两端
现在的目标是找出分类器定义中的W和b.为此，我们必须找到支持向量，一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化:

label*(W转置X+b)：称为点到分隔面的函数间隔
(label*(W转置X+b))/||W||：称为点到分隔面的几何间隔
引入拉格朗日乘子alpha[i],一番计算之后得到优化问题的最终结果

满足于条件：Σalpha[i]y[i] = 0,alpha>=0,i=1,2,…,n
alpha:拉格朗日乘子
y:label标签：-1，1
x:特征,x1,x2

6.2.2 SVM应用的一般框架

1.收集数据
2.准备数据：需要数值型数据
3.分析数据：有助于可视化分隔超平面
4.训练算法：SVM的大部分时间都源自训练，该过程主要实现两个参数的调优
5.测试算法：十分简单的计算过程就可以实现
6.使用算法：几乎所有分类问题都可以使用SVM，值得一提的是，SVM本身是一个二类分类器，对多类问题应用SVM需要对代码进行修改

6.3 SMO高效优化算法

6.3.1 platt的SMO算法

    SMO:序列最小优化(sequential Minimal Optimization)

SMO算法是将大优化问题分解为多个小优化问题来求解的。这些小优化问题往往很容易求解，并且对它们进行顺序求解的结果与将它们作为整体求解的结果是一致的。
SMO算法的目标是求出一系列alpha和b
SMO算法的工作原理：每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha,那么就增大其中一个同时减少另一个。这里所谓的“合适”就是指两个alpha必须要符合一定的条件，条件之一就是这两个alpha必须要在间隔边界之外，二其第二个条件则是这两个alpha还没有进行过区间化处理或者不在边界上。
解决一个多变量的优化问题是困难的，但如果能够多次迭代，每次选择两个变量优化，同时把其他变量看成是固定的，这样“分而治之”的话，要容易很多。根据Σalpha[i]y[i] = 0(拉格朗日乘法所决定的),故一次更新两个alpha,其他alpha暂定为常数，由等式确定一个alpha后，另一个alpha也能确定。

6.3.2 应用简化版SMO算法处理小规模数据集

platt SMO算法中的外循环确定要优化的最佳alpha对，而简化版跳过这一部分，首先在数据集上遍历每一个alpha,然后在剩下的alpha集合中随机选择另一个alpha,从而构建alpha对。这里有一点非常重要，就是我们要同时改变两个alpha,之所以这样做是因为我们有一个约束条件：
Σαlabel[i] = 0
伪代码
创建一个alpha向量并将其初始化为0向量
当迭代次数小于最大迭代次数时(外循环)
对数据集中的每个数据向量(内循环)：
如果该数据向量可以被优化：
随机选择另外一个数据向量
同时优化这两个向量
如果两个向量都不能被优化，退出内循环
如果所有向量都没被优化，增加迭代数目，继续下一次循环

#简化的SMO算法
class SimplfiedSMO():
    """简化的SMO算法"""
    #SMO算法中的辅助函数
    def loadDataSet(self,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 selectJrand(self,i,m):
        """在区间范围内随机选择一个整数
        i:是第一个alpha的下标\n
        m:是所有alpha的数目"""
        j = i
        #只要函数值不等于输入值i,函数就会进行随机选择
        while(j==i):
            j = int(random.uniform(0,m))
        return j
    def clipAlpha(self,aj,H,L):
        """调整α值的大小"""
        if aj > H:
            aj = H
        if L > aj:
            aj = L
        return aj
    def smoSimple(self,dataMatIn,classLabels,C,toler,maxIter):
        """简化的SMO算法
        dataMatIn:数据集
        classLabels：类别标签
        C：常数,用于控制“最大化间隔”和“保证大部分点的函数间隔小于1.0-function margin”这两个目标的权重
        toler：容错率
        maxIter：退出前最大的循环次数"""
        dataMatrix = np.mat(dataMatIn)
        labelMat = np.mat(classLabels).transpose()  #转换为列向量
        b = 0
        m,n = np.shape(dataMatrix)  #m为实例个数，n为特征个数
        alphas = np.mat(np.zeros((m,1)))    #alpha列矩阵
        iter = 0    #存储的是在没有任何alpha改变的情况下遍历数据集的次数，当该变量达到输入值maIter时，函数结束运行
        while(iter<maxIter):
            alphaPairsChanged = 0   #记录是否已经优化
            for i in range(m):
                """fXi是我们预测的类别，然后，基于这个实例的预测结果和真实结果的比对，就可以计算误差Ei.如果误差很大，那么可以对该数据实例所对应的alpha值进行优化"""
                fXi = float(np.multiply(alphas,labelMat)).T*(dataMatrix*dataMatrix[i,:].T) +b
                Ei = fXi - float(labelMat[i])
                #如果alpha可以更改进入优化过程
                """在if语句中，不管是正间隔还是负间隔都会被测试，并且在该if语句中，也要同时检查alpha值，以保证其不能等于0或C。由于后面alpha小于0或大于C时将被调整为0或C，所以一旦在该if语句中它们等于这两个值的话，那么它们就已经在‘边界’上了，因而不饿能再继续减小或增大，这时，利用selectJrand()函数随机选择第二个alpha值"""
                if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                    #随机选择第二个alpha
                    j = self.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()
                    #保证alpha在0与C之间
                    if (labelMat[i] != labelMat[j]):
                        L = max(0,alphas[j]+alphas[i]-C)
                        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[j,:]*dataMatrix[j,:].T   #是alpha[j]的最优修改量，如果eta为0，则退出for循环
                    if eta >= 0 :
                        print("eta>=0")
                        continue
                    alphas[i] -= labelMat[j]*(Ei-Ej)/eta
                    alphas[j] = self.clipAlpha(alphas[j],H,L)
                    if abs(alphas[j]-alphaJold) < 0.00001:
                        print("j not moving enough")
                        continue
                    #对alphas[i]进行修改，修改量与j相同，但方向相反
                    alphas[i] += labelMat[j]*labelMat[i]*(alphaJold-alphas[j])  #alphaJold-alphas[j]为alpha[j]的修改量
                    #优化过后，为两个alpha值设置一个常数项b.
                    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
                    if 0<alphas[i] and C>alphas[i]: #只有alpha[i]还在优化区间的时候才加上常数值
                        b = b1
                    elif 0<alphas[j] and C>alphas[j]:
                        b = b2
                    else:
                        b = (b1+b2)/2.0
                    alphaPairsChanged += 1
                    print('iter:%d i:%d,pairs changed %d' % (iter,i,alphaPairsChanged))
            #如果程序执行到for循环的最后一行都不执行continue语句，那么就已经成功地改变了一对alpha,同时可以增加alphaPairsChanged的值。
            if alphaPairsChanged == 0:  #只有在所有数据集上遍历maxIter次，且不再发生任何alpha修改之后，程序才会停止并退出while循环
                iter += 1
            else:
                iter = 0
            print('iteration number:%d' % iter)
        return b,alphas
if __name__ == '__main__':
    load_data = SimplfiedSMO()
    dataArr,labelArr = load_data.loadDataSet('./M_06_SVM/testSet.txt')
    b,alphas = smoP(dataArr,labelArr,0.6,0.001,40)
    weights = calcWs(alphas,dataArr,labelArr)
    print("============开始测试训练结果============")
    #测试数据
    testX = np.mat(dataArr[0])
    print(testX)
    result = testX*weights + b  #y=WX+b
    if result < 0 and labelArr[0] == -1:
        print("该数据的标签是",labelArr[0])
        print(-1)
    else:
        print(1)

在这里插入图片描述

6.4 利用完整Platt SMO算法加速优化

两个版本的alpha更改和代数运算的优化环节一模一样

Platt SMO算法是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间进行交替：
一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界alpha中实现单遍扫描。而所谓的非边界alpha指的就是那些不等于边界0或C的alpha值。实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后再对这个表进行遍历。同时，该步骤会跳过那些已知的不会改变的alpha值。
选择第一个alpha值后，算法会通过一个内循环来选择第二个
alpha值：
在优化过程中，会通过最大化步长的方式来获得第二个alpha值。在简化版SMO算法中，我们会选择j之后计算错误率Ej.但在完整的SMO算法中，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说Ei-Ej最大的alpha值。
C值
常数C给出的是不同优化问题的权重。常数C一方面要保障所有样例的间隔不小于1.0，另一方面又要使得分类间隔要尽可能大，并且要在这两方面之间平衡。如果C很大，那么分类器将力图通过分隔超平面对所有的样例都正确分类。

#完整的SMO算法
#完整版platt SMO的支持函数
class optStruct():
    """求解器。建立一个数据结构——对象保存所有的重要值。可以通过pyhon字典完成值的转递
    dataMatIn:np.array数据集,classLabels：数据标签,
    C：惩罚参数,常数,用于控制“最大化间隔”和“保证大部分点的函数间隔小于1.0-function margin”这两个目标的权重
    toler：容错率,alpha调整的阈值"""
    def __init__(self,dataMatIn:np.array,classLabels:np.array,C,toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = np.shape(dataMatIn)[0]
        #初始alphas是全0的矩阵，alphas[i]从0开始更新，alphas[j]根据误差进行更新
        #oS.alphas[j] -= oS.labelMat[j]*(Ei-Ej)/eta
        #oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        self.alphas = np.mat(np.zeros((self.m,1)))
        self.b = 0
        self.eCache = np.mat(np.zeros((self.m,2)))  #误差缓存
def calcEk(oS,k):
    """对于给定的alpha值，返回对应的误差"""
    # latter = oS.X*oS.X[k,:].T
    # print(latter)
    # former = np.multiply(oS.alphas,oS.labelMat).T
    # print(former)
    # fXk = former*latter
    # print(fXk)
    # fXk = float(fXk)
    # fXk = fXk + oS.b
    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 clipAlpha(aj,H,L):
        """调整α值的大小"""
        if aj > H:
            aj = H
        if L > aj:
            aj = L
        return aj
def selectJrand(i,m):
        """在区间范围内随机选择一个整数
        i:是第一个alpha的下标\n
        m:是所有alpha的数目"""
        j = i
        #只要函数值不等于输入值i,函数就会进行随机选择
        while(j==i):
            j = int(random.uniform(0,m))
        return j
def selectJ(i,oS,Ei):
    """用于选择第二个alpha或者说内循环的alpha值。这里的目标是选择合适的第二个alpha值以保证在每次优化中采用最大步长。该函数的误差值与第一个alpha值Ei和下标i有关"""
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1,Ei]
    validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]   #返回非零E值所对应的alpha本身
    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):
    """计算误差值并存入缓存当中，在对alpha值进行优化之后会用到这个值"""
    Ek = calcEk(oS,k)
    oS.eCache[k] = [1,Ek]
#内循环函数
def innerL(i,oS):
    """内循环：寻找一次迭代中所有数据的最小geometry margin.完整Platt SMO算法中的优化例程:更新alpha[i],alpha[j]的一种优化方法"""
    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() #copy()复杂的object，如list中套着list的情况，shallow copy中的子list，并未从原object真的「独立」出来。也就是说，如果你改变原object的子list中的一个元素，你的copy就会跟着一起变。这意味着之前优化了的alphas[i]随着优化一直在跟着变化
        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  #从中选择使得步长或者说Ei-Ej最大的alpha值
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS,j)
        if abs(oS.alphas[j] - alphaJold) < 0.00001:
            print('j not moving enough')
            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,kTup=('lin',0)):
    """完整的外循环：外循环寻找所有迭代中最大的geometry margin\n
    return oS.b,oS.alphas"""
    print(np.shape(np.mat(dataMatIn)))
    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):   #内循环
                """调用InnerL()来选择第二个alpha[j],并在可能时对其进行优化处理。如果有任意一对alpha值发生改变，那么返回1"""
                alphaPairsChanged += innerL(i,oS)
            print('fullSet,iter:%d i:%d,pairs changed %d' % (iter,i,alphaPairsChanged))
            iter += 1
        else:
            nonBoundIs = np.nonzero((oS.alphas.A > 0)* (oS.alphas.A < C))[0]  #将矩阵转换为数组.A.#nonzero()返回非零值所对应的alpha本身
            print(nonBoundIs)
            for i in nonBoundIs:
                """遍历所有的非边界alpha值，也就是不在边界0或C上的值"""
                alphaPairsChanged += innerL(i,oS)
                print('non-bound,iter: %d i:%d,pairs changed %d' % (iter,i,alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False   #切换整个循环到NonboundIs中
        elif alphaPairsChanged == 0:
            entireSet = True    #如果alpha不再变化，完成一次迭代
        print('iteration number:%d' % iter)
    return oS.b,oS.alphas
def calcWs(alphas,dataArr,classLabels):
    """基于alpha值得到超平面
    return weights"""
    X = np.mat(dataArr)
    labelMat = np.mat(classLabels).transpose()
    m,n = np.shape(X)
    weights = np.zeros((n,1)) #同特征数量相同
    for i in range(m):
        weights += np.multiply(alphas[i]*labelMat[i],X[i,:].T)#大部分呢alpha值为0，而非零alpha所对应的也就是支持向量
    return weights

在这里插入图片描述

6.5 在复杂数据上应用核函数(kernel)

在这里插入图片描述
前面我们讨论了线性可分的情况，对于非线性可分的数据又怎么分类呢？

6.5.1 利用核函数将数据映射到高维空间

我们将数据从一个特征空间转换到另一个特征空间
通过核函数，我们能实现这样的空间转换，可以把核函数想象成一个包装器或接口，它能把数据从某个很难处理的形式转换成为另一个较容易处理的形式。映射到高维空间后，我们可以在高维空间中解决线性问题，等价于低维空间中解决非线性问题
SVM优化中一个特别好的地方是，所有的运算都可以写成内积(inner product,也称点积)的形式。向量的内积指的是两个向量相乘之后得到单个标量。我们可以把内积运算替换成核函数，而不必做简化处理。将内积替换成核函数的方式被称为核技巧(kernel trick),很多能写成内积形式的算法都可以使用核技巧。

6.5.2 径向基核函数

径向基函数是一个采用向量作为自变量的函数，能够基于向量距离运算输出一个标量。这个距离可以是从<0,0>向量或者其他向量开始计算的距离。接下来，我们将会使用到径向基函数的高斯版本，具体公式是
在这里插入图片描述
其中,σ是用户定义的用于确定到达率(reach)或者说函数值跌落到0的速度参数
上述高斯核函数将数据从其特征空间映射到更高维的空间，具体来说这里是映射到一个无穷维空间。
K[j] = deltaRow * deltaRow.T #将特征映射到高维x[i],x[j],K[j].类比y=ax+b,包含两个变量(x,y)二维空间，(x[i],x[j],k[j])三维空间，在我们的径向基函数中，第三维就是两个向量的距离
如果使用exp()函数，最终会将结果映射到无穷维

def kernelTrans(X:np.mat,A:np.mat,Ktup:tuple):
    """ 核转换函数:将特征从低维空间映射到高维空间

    X:数据集
    A:数据集的第i行，即某个传入的example实例
    Ktup:是一个包含核函数信息的元组。核函数类型：lin:线性，rbf:高斯核，σ-Ktup[1]：由用户定义的用于确定到达率(reach)或者说函数值跌落到0的速度参数

    返回(m,1)经过转换后得到的数组K"""
    m,n = np.shape(X)   #(m,n)
    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    #将特征映射到高维x[i],x[j],K[j].类比y=ax+b,包含两个变量(x,y)二维空间，(x[i],x[j],k[j])三维空间，在我们的径向基函数中，第三维就是两个向量的距离
        K = exp(K/(-1*Ktup[1]**2))  #高斯核转换，对矩阵元素展开计算
    else:
        raise NameError('Houston We have a Problem,the kernel is not recognized')
    return K
class optStruct():
    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)))
        self.b = 0
        self.eCache = np.mat(np.zeros((self.m,2)))  
        self.K = np.mat(np.zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X,self.X[i,:],Ktup)

6.5.3 在测试中使用核函数

    通过调整σ来调整训练和测试的速度，得到不同的训练错误率和测试错误率，找到最优的σ值。

支持向量的数目存在一个最优值。SVM的优点在于它能对数据进行高效分类。如果支持向量太少，就可能得到一个很差的决策边界；如果支持向量太多，也就相当于每次都利用整个数据集进行分类，这种分类方法称为K近邻。

def testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('./M_06_SVM/testSetRBF.txt')
    b,alphas = smoP(dataArr,labelArr,200,0.00001,10000,('rbf',k1))
    dataMat = np.mat(dataArr)
    labelMat = np.mat(labelArr).transpose()
    #构建支持向量矩阵
    svInd = np.nonzero(alphas.A > 0)[0] #返回非零alpha值所在的行索引
    print(svInd)
    sVs = dataMat[svInd]    #只使用支持向量进行训练
    labelSV = labelMat[svInd]
    print(labelSV)
    print('there are %d Support Vectors' % np.shape(sVs)[0])
    m,n = np.shape(dataMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,dataMat[i,:],('rbf',k1))
        print(kernelEval.shape)   #(支持向量个数,1)行数同支持向量个数相同
        print(type(kernelEval)) #<class 'numpy.ma.core.MaskedArray'>
        #weights += np.multiply(alphas[i]*labelMat[i],X[i,:].T)#大部分呢alpha值为0，而非零alpha所对应的也就是支持向量
        predict = np.mat(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))
    testDataArr,testLabelArr = loadDataSet('./M_06_SVM/testSetRBF2.txt')
    errorCount = 0
    testDataMat = np.mat(dataArr)
    testLabelMat = np.mat(labelArr).transpose()
    m,n = np.shape(testDataMat)
    for i in range(m):
        kernelEval - kernelTrans(sVs,testDataMat[i,:],('rbf',k1))
        predictTset = kernelEval.T*np.multiply(labelSV,alphas[svInd]) + b
        if np.sign(predict) != np.sign(testLabelArr[i]):
            errorCount += 1
    print("the test error rate is :%f" % (float(errorCount)/m))

6.6 示例：手写体识别问题回顾

因为可以仅用支持向量样本进行测试，所以相比KNN算法，SVM算法占用内存更低

import random
import numpy as np
import copy
from numpy.ma import exp
from M_0601_SVM_SMO import selectJ,clipAlpha,updateEk,calcWs
import os
import matplotlib.pyplot  as plt
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 kernelTrans(X:np.mat,A:np.mat,Ktup:tuple):
    """ 核转换函数

    X:数据集
    A:数据集的第i行，即某个传入的example实例
    Ktup:是一个包含核函数信息的元组。核函数类型：lin:线性，rbf:高斯核，σ-Ktup[1]：由用户定义的用于确定到达率(reach)或者说函数值跌落到0的速度参数

    返回(m,1)经过转换后得到的数组K"""
    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 = exp(K/(-1*Ktup[1]**2))  #高斯核转换，对矩阵元素展开计算
    else:
        raise NameError('Houston We have a Problem,the kernel is not recognized')
    return K
class optStruct():
    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)))
        self.b = 0
        self.eCache = np.mat(np.zeros((self.m,2)))  
        self.K = np.mat(np.zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X,self.X[i,:],Ktup)

#使用核函数需要重写innerL()及calcEk()函数
def calcEk(oS,k):
    """对于给定的alpha值，返回对应的误差"""
    #fXk = float(np.multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b 
    fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek
def innerL(i,oS):
    """内循环：寻找一次迭代中所有数据的最小geometry margin.完整Platt SMO算法中的优化例程:更新alpha[i],alpha[j]的一种优化方法"""
    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   
        #决策函数
        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] 
        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('j not moving enough')
            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.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]
        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,Ktup=('lin',0)):
    """完整的外循环：外循环寻找所有迭代中最大的geometry margin\n
    return oS.b,oS.alphas"""
    print(np.shape(np.mat(dataMatIn)))
    oS = optStruct(np.mat(dataMatIn),np.mat(classLabels).transpose(),C,toler,Ktup=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):   #内循环
                """调用InnerL()来选择第二个alpha[j],并在可能时对其进行优化处理。如果有任意一对alpha值发生改变，那么返回1"""
                alphaPairsChanged += innerL(i,oS)
            print('fullSet,iter:%d i:%d,pairs changed %d' % (iter,i,alphaPairsChanged))
            iter += 1
        else:
            nonBoundIs = np.nonzero((oS.alphas.A > 0)* (oS.alphas.A < C))[0]  #将矩阵转换为数组.A.#nonzero()返回非零值所对应的alpha本身
            print(nonBoundIs)
            for i in nonBoundIs:
                """遍历所有的非边界alpha值，也就是不在边界0或C上的值"""
                alphaPairsChanged += innerL(i,oS)
                print('non-bound,iter: %d i:%d,pairs changed %d' % (iter,i,alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False   #切换整个循环到NonboundIs中
        elif alphaPairsChanged == 0:
            entireSet = True    #如果alpha不再变化，完成一次迭代
        print('iteration number:%d' % iter)
    return oS.b,oS.alphas
def testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('./M_06_SVM/testSetRBF.txt')
    b,alphas = smoP(dataArr,labelArr,200,0.00001,10000,('rbf',k1))
    dataMat = np.mat(dataArr)
    labelMat = np.mat(labelArr).transpose()
    #构建支持向量矩阵
    svInd = np.nonzero(alphas.A > 0)[0] #返回非零alpha值所在的行索引
    print(svInd)
    sVs = dataMat[svInd]    #只使用支持向量进行训练
    labelSV = labelMat[svInd]
    print(labelSV)
    print('there are %d Support Vectors' % np.shape(sVs)[0])
    m,n = np.shape(dataMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,dataMat[i,:],('rbf',k1))
        print(kernelEval.shape)   #(支持向量个数,1)行数同支持向量个数相同
        predict = np.mat(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))
    testDataArr,testLabelArr = loadDataSet('./M_06_SVM/testSetRBF2.txt')
    errorCount = 0
    testDataMat = np.mat(dataArr)
    tettLabelMat = np.mat(labelArr).transpose()
    m,n = np.shape(testDataMat)
    for i in range(m):
        kernelEval - kernelTrans(sVs,testDataMat[i,:],('rbf',k1))
        predictTset = kernelEval.T*np.multiply(labelSV,alphas[svInd]) + b
        if np.sign(predict) != np.sign(testLabelArr[i]):
            errorCount += 1
    print("the test error rate is :%f" % (float(errorCount)/m))
def img2vector(filename):
    """准备数据：将图片转换为向量，循环读出文件的前32行，并将每行的头32个字符值存储在Numpy数组中"""
    return_vect = np.zeros((1,1024))
    fr = open(filename)
    for i in range(32):
        line_str = fr.readline()
        for j in range(32):
            return_vect[0,32*i+j] = int(line_str[j])
    return return_vect
def loadImages(dirname):
    """"同KNN中的数字识别不同，在KNN中直接使用类别标签，而同支持向量机一起使用时，类别标签为-1或者+1.因此一旦碰到数字9，则输出类别标签-1，否则输出+1.由于这里我们只做二类分类，因此除了1和9意外的数字都被去掉了。"""
    hwLabels = []
    trainingFileList = os.listdir(dirname)
    m = len(trainingFileList)
    trainingMat = np.zeros((m,1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('_')[0]
        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)):
    """训练集得到SVM的aplha权重以及b,然后"""
    dataArr,labelArr = loadImages('./M_06_SVM/digits/trainingDigits')
    #训练得到b,alpha
    b,alphas = smoP(dataArr,labelArr,200,0.0001,1000,kTup)
    dataMat = np.mat(dataArr)
    labelMat = np.mat(labelArr).transpose()
    #支持向量矩阵
    svInd = np.nonzero(alphas.A > 0)[0]
    sVs = dataMat[svInd]
    labelSV = labelMat[svInd]
    print("there are %d support Vectors" % np.shape(sVs)[0])
    m,n = np.shape(dataMat)
    trainErrorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,dataMat[i,:],kTup)
        predict = np.mat(kernelEval.T)*np.multiply(labelSV,alphas[svInd]) + b
        if np.sign(predict) != np.sign(labelArr[i]):
            trainErrorCount += 1
    print("the training error rate is:%f" % (float(trainErrorCount)/m))
    dataArr,labelArr = loadImages('./M_06_SVM/digits/testDigits')
    testErrorCount = 0
    dataMat = np.mat(dataArr)
    labelMat = np.mat(labelArr).transpose()
    m,n = np.shape(dataMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,dataMat[i,:],kTup)
        predict = np.mat(kernelEval.T)*np.multiply(labelSV,alphas[svInd]) + b
        if np.sign(predict) != np.sign(labelArr[i]):
            testErrorCount += 1
    print("the test error rate is:%f" % (float(testErrorCount)/m))
    return b,alphas,trainErrorCount,testErrorCount
if __name__ == '__main__':
    #testRbf()
    #尝试不同的σ-径向基核函数中的速度参数，观察不同σ对准确率的影响
    trainErrorCounts = []
    testErrorCounts = []
    speedParams = [0.1,1,5,10,20,50,100]
    for xigema in speedParams:
        b,alphas,trainErrorCount,testErrorCount = testDigits() #保存模型训练结果
        trainErrorCounts.append(trainErrorCount)
        testErrorCounts.append(testErrorCount)
    plt.plot(speedParams,trainErrorCounts)
    plt.show()

6.7 本章小结

核方法或者说核技巧会将数据从一个低维数据映射到一个高维空间，可以将一个在低维空间中的非线性问题转换成高维空间下的线性问题来求解。核方法不止在SVM中使用，还可以用于其他算法。而其中的径向基函数是一个常用的度量两个向量距离的核函数。
支持向量机是一个二类分类器。当用其解决多类问题时，则需要额外的方法对其进行扩展。SVM的效果也对优化参数和所用核函数中的参数敏感。