机器学习实践之支持向量机学习

       本文根据最近学习机器学习书籍网络文章的情况,特将一些学习思路做了归纳整理,详情如下.如有不当之处,请各位大拿多多指点,在此谢过。

       第一部分第3点 寻找最大间隔，部分图片手写，图片没有经过很好的处理，如有不美观之处，请谅解。

一  概述

 1、概念理解

     支持向量机（SupportVectorMachines,SVM），由于理解支持向量机需要一定的相关理论基础，这些理论基础有些复杂，笔者这里只是做些相关概念性的介绍，以后有机会专门介绍深入学习支持向量机相关知识点。简单说，支持向量机就是离分隔超平面最近的点。SVM其实有很多种实现，这里只介绍一种较为流行的实现方法，即序列最小优化（SequentialMinimal Optimization，SMO）算法。后面还会介绍一种称为核函数（kernel）的方式将SVM扩展到更多数据集上。

2、 基于最大间隔分隔数据

    如下图所示，图中的圆形点和方形点是否可以被一条直线（二维空间中是一条线，三维空间是一个平面）分开？由于两者之间分隔距离已经足够多，所以用一条直线把它们分开并不难。在这种情况下，这组数据被称为线性可分（linearlySeparable）数据。

          

      上述将数据集分隔开来的直线被称为分隔超平面（separatinghyperplane）。这里所给的数据集是在二维平面上，所以这里的分隔超平面是一条直线。如果给出的数据集是三维的，则分隔它们的将是一个平面。更高维的数据集以此类推。如果给出的数据集是1024维的数据，则需要一个1023维的某某对象来对数据进行分隔。这个1023维的某某对象就被称为超平面（hyperplane），也就是分类的决策边界。分布在超平面一侧的所有数据属于某一类，分布在另一侧的所有数据属于另一类别。所以，我们希望通过这种方式来建立分类器，即数据点离决策边界越远，最后预测的结果就越可信。

      而这里点到分隔面的距离称为间隔。考虑到我们自身可能犯错或者在有限数据集上训练分类器的情况出现，所以我们希望间隔尽可能要大。

     

 3、寻找最大间隔

     在寻找最大间隔之前，先看下图，有助于理解后面的内容。

               

    

       这里假设苹果和香蕉都是带炸药的水果，且两种水果中间有一条安全警戒线，苹果和香蕉距离警戒线足够近的话，就会引起爆炸伤及对方，但同类别水果内部不会引起爆炸。

      这里可以想一下，将一水果数据分给分类器，就会输出该水果对应的类别标签（苹果或香蕉），这相当于一个Sigmoid函数在作用。这里与Logistic回归有所不同的是，类别标签为-1和+1。

       当计算数据点到分隔面的距离并确定分隔面的位置时，间隔通过

                    

   来计算，这时类别选取-1和+1的优势就体现出来了。所以，我们只需要找出分类器中点w和b就可以了。为了达到这一目的，我们必须找到具有最小间隔的数据点，而这些点就是前面提到的支持向量。我们一旦找到具有最小间隔的数据点，接下来就要对该间隔进行最大化。表达式如下

                   。

       当然直接求解上式难度较大，可以考虑通过拉格朗日乘子对其进行转化优化，这里有个约束条件：

                。

         通过引入拉格朗日乘子，我们就可以基于约束条件来表达原来的问题。由于这里的约束条件都是基于数据点的，所以我们可以将超平面写成数据点的形式。目标函数被优化后变为：

           

        约束条件为：

          

         整体看起来相当不错，但有个假设条件就是数据必须100%线性可分。众所周知，实际工作中数据可以说都不是很“纯净”，解决办法是引入所谓松弛变量（slackvariable），来允许一些数据点可以处于分隔面的错误一侧。结果是，我们一直优化的目标仍然保持不变，但约束条件变为：

      

       这里C作为常数，用于控制“最大化间隔”和“保证大部分点的函数间隔小于1.0”这两个目标的权重。     一旦求出所有的alpha，那么分隔超平面就可以通过这些alpha来表达。

4、SMO高效优化算法

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

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

       SMO算法的工作原理：每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha后，就增大其中一个同时减小另一个。这里所谓的“合适”，是指这些alpha必须符合一些约束条件。一是这些alpha必须在间隔边界之外；二是这些alpha还没有进行区间化优化处理或者不在边界之上。

   

5、SVM应用的一般过程

     (1)收集数据:可以使用任何方法。

    (2)准备数据:需要数值型数据。

    (3)分析数据:有助于可视化分隔超平面。

   (4)训练算法:SVM的大部分时间都来自于训练，训练的过程主要实现两个参数的调优。

    (5)测试算法:十分简单的计算过程就可以实现。

    (6)使用算法:几乎所有的分类问题都可以使用SVM，而且，SVM本身就是一个二分类分类器，对于多类问题应用SVM的话，需要在代码实现上做些微调。

6、SVM的相关特性

      优点：泛化错误率低，计算开销不大，结果容易解释。

      缺点：对参数调节和核函数的选择较为敏感，原始分类器不做任何修改仅适用于处理二分类问题。

      适用数据类型：数值型和标称型。

二 项目一 小规模数据集上的SMO算法

1、项目概述

    对小规模数据点进行分类。

2、数据集格式

      

7.286357        0.251077        1
2.301095        -0.533988       -1
-0.232542       -0.547690       -1
3.457096        -0.082216       -1
3.023938        -0.057392       -1

3、准备数据

     读取相关数据。

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

    

def loadDataSet(fileName):
    """
    对文件进行逐行解析，从而得到第某行的类标签和整个特征矩阵
    Args：
         fileName  文件名
    Returns:
        dataMat 特征矩阵
        labelMat 类标签
    """
    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(i,m):
    """
    随机选择一个整数
    Args：
         i： 第一个alpha的下标
         m： 所有alpha的数目
    Return：
         j： 返回一个不为i的随机数，在0～m之间的整数值。
    """
    j = i
    while(j == i):
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):
    """
    clipAlpha(调整aj的值，使aj处于L<aj<H范围)
    Args： 
         aj： 目标值
         H： 最大值
         L： 最小值
    Return
         aj： 目标值
    """
    
    if aj > H:
        aj=H
    if L > aj:
        aj = L
    return aj

    

       这里如果我们查看数据集中的类别，可以看到如下结果。

            

[-1.0,
 -1.0,
 1.0,
 -1.0,
 1.0,
 1.0,
 1.0,
 -1.0,
 -1.0,
 -1.0,
 -1.0,
 -1.0,
 -1.0,
 1.0,
 -1.0,
 1.0,
 1.0,
 -1.0,
 1.0,
 -1.0,
 -1.0,
...

4、分析数据

            这里暂无相关内容。

5、训练算法

   （1） 伪代码如下：

    创建一个alpha向量并将其初始化为0向量

    当迭代次数小于最大迭代次数时（外循环）

                对数据集中的每个向量（内循环）：

                        如果该数据向量可以被优化：

                              随机选择另外一个数据向量

                              同时优化这两个向量

                              如果这两个向量都不能被优化，则退出内循环

                如果所有数据向量都没有被优化，增加迭代次数，进行下一次循环

      (2）实现代码（简化版）如下：

             

#简化版 SMO算法
def smoSimple(dataMatIn,classLabels,C, toler, maxIter):
    
    """
       smoSimple
    Args:
         dataMatIn   特征集合
         classLabels  类别标签
         C      松弛变量常量，允许有些数据点处于分隔面的错误一侧；控制最大化间隔和保证大部分点的函数间隔小于1.0这两个目标的权重。
         toler  容错率（指在某个体系中能减少一些因素或选择对某个系统产生不稳定的概率）
         maxIter 退出前最大的循环次数
    Returns:
         b  模型的常量值
         alphas  拉格朗日乘子
    """
    dataMatrix = mat(dataMatIn) #矩阵转换，和 .T一样的功能。
    labelMat = mat(classLabels).transpose()
    b = 0; m,n = shape(dataMatrix) #初始化b和alpha值，alpha有点类似于权重值。
    alphas = mat(zeros((m,1)))
    iter = 0  #没有任何alpha改变的情况下，遍历数据的次数。
    while(iter < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            fXi = float(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)):
                
                j = selectJrand(i,m)
                              
                fXj = float(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 = 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
                if (eta>=0):print("eta>=0");continue
                alphas[j] -= labelMat[j]*(Ei-Ej)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j]-alphaJold)<0.00001):print("j not moving enough "); continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold-alphas[j])
                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]):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))
        if(alphaPairsChanged == 0): iter +=1
        else: iter = 0
        print ("iteration number : %d" % iter)
    return b,alphas

        这里smoSimple（）函数有五个输入参数：数据集、类别标签、常数C、容错率、退出前最大的循环次数。上面函数将多个列表和输入参数转换成Numpy矩阵，这样一来，就可以简化很多数学处理操作。每次循环中，将alphaPairsChanged先设为0，然后对整个数据集进行遍历。alphaPairsChanged用于记录alpha是否已经被优化。

         
    在原始数据集上对这些支持向量进行画圈之后的效果如下：

       

      (3）实现代码（完整版）如下：

        

#完整版Platt SMO算法中的优化例程

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("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
    
    
# 完整版Platt SMO的外循环

def smoP(dataMatIn, classLabels, C, toler, maxIter):
    oS = optStruct(mat(dataMatIn), 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("fullSet, iter : %d i: %d, pairs changed %d" % (iter, i, alphaPairsChanged))
            iter += 1
            
        else:
            nonBoundIs = nonzero((oS.alphas.A >0 )* (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
                print ("non-bound, iter: %d i: %d, pairs changed %d" % (iter ,i ,alphaPairsChanged))
                
            iter += 1
            
        if entireSet: entireSet = False
        elif(alphaPairsChanged == 0): entireSet = True
        
        print("iteration number : %d" % iter)
    return oS.b,oS.alphas

      在原始数据集上对这些支持向量进行画圈之后的效果如下：

        

         

      （4）对比简化版和完整版

•     在实现alpha的更改和代数运算优化环节都是一样的。在优化过程中，唯一不同的就是alpha的选择方式。完整版里应用里一些能够提速的启发方法。

•     常数C一方面要保障所有样例的间隔不小于1.0，另一方面又要使得分类间隔尽可能大，且要在这两个之间进行平衡。如果C值很大，则分类器将力图通过分隔超平面对所有样例都正确分类。这种优化运行的结果，如完整版运行后的效果图所示，显然比简化版运行后的效果图拥有更多的支持向量。究其原因，是因为简化版算法中是通过随机方式选择alpha对的，不如完整版覆盖整个数据集效果好。

• 如果数据集是非线性可分的，就会发现支持向量会在超平面附近聚集成团。

三   项目二 手写识别（有核）

1、项目概述

    结合笔者第一篇关于KNN算法中的手写识别示例。后来上级说：你写的那个手写识别程序不错，但却非常占内存。如何做到在保持其性能的同时，减少所占内存呢？

2、 数据集格式

       

00000000000001100000000000000000
00000000000011111100000000000000
00000000000111111111000000000000
00000000011111111111000000000000
00000001111111111111100000000000
00000000111111100011110000000000
00000001111110000001110000000000
00000001111110000001110000000000
00000011111100000001110000000000
00000011111100000001111000000000
00000011111100000000011100000000
00000011111100000000011100000000
00000011111000000000001110000000
00000011111000000000001110000000
00000001111100000000000111000000
00000001111100000000000111000000
00000001111100000000000111000000
00000011111000000000000111000000
00000011111000000000000111000000
00000000111100000000000011100000
00000000111100000000000111100000
00000000111100000000000111100000
00000000111100000000001111100000
00000000011110000000000111110000
00000000011111000000001111100000
00000000011111000000011111100000
00000000011111000000111111000000
00000000011111100011111111000000
00000000000111111111111110000000
00000000000111111111111100000000
00000000000011111111110000000000
00000000000000111110000000000000

3、开发流程

 

      (1)收集数据:提供的文本文件。

    (2)准备数据:基于二值图像构造向量。

    (3)分析数据:对图像向量进行目测。

   (4)训练算法:  采用两种不同的核函数，并对径向基核函数采用不同的设置来运行SMO算法。

    (5)测试算法： 编写一个函数来测试不同的核函数并计算错误率。

    (6)使用算法: 一个图像识别的完整应用还需要一些图像处理的知识，不做过深说明。

4、代码实现

            代码实现过程中，涉及与SMO完整版代码一样，读者可以自行填充，这里就不再贴出。

      

def kernelTrans(X, A, kTup):  # calc the kernel or transform data to a higher dimensional space
    """
    核转换函数
    Args:
        X     dataMatIn数据集
        A     dataMatIn数据集的第i行的数据
        kTup  核函数的信息
    Returns:
    """
    m, n = shape(X)
    K = mat(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 -- That Kernel is not recognized')
    return K


def calcWs(alphas, dataArr, classLabels):
    """
    基于alpha计算w值
    Args:
        alphas        拉格朗日乘子
        dataArr       feature数据集
        classLabels   目标变量数据集
    Returns:
        wc  回归系数
    """
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m, n = shape(X)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w


def testRbf(k1=1.3):
    dataArr, labelArr = loadDataSet('testSetRBF.txt')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))  
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd] 
    labelSV = labelMat[svInd]
    print("there are %d Support Vectors" % shape(sVs)[0])
    m, n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))

        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the training error rate is: %f" % (float(errorCount) / m))

    dataArr, labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]):
            errorCount += 1
    print("the test error rate is: %f" % (float(errorCount) / m))


def img2vector(filename):
    returnVect = 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 = []
    print(dirName)
    trainingFileList = listdir(dirName)  # load the training set
    m = len(trainingFileList)
    trainingMat = 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)):

    # 1. 导入训练数据
    dataArr, labelArr = loadImages('trainingDigits')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]
    labelSV = labelMat[svInd]
    # print("there are %d Support Vectors" % shape(sVs)[0])
    m, n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]): errorCount += 1
    print("the training error rate is: %f" % (float(errorCount) / m))

    # 2. 导入测试数据
    dataArr, labelArr = loadImages('testDigits')
    errorCount = 0
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]): errorCount += 1
    print("the test error rate is: %f" % (float(errorCount) / m))

   

四   小结

   

     支持向量机是一种分类器，因为它会产生一种二值决策结果，即它是一种决策“机”，所以它被称为“机”就可以理解。支持向量机的泛化错误率低，换句话讲，它具有良好的学习能力，并且学习到的结果具有良好的推广性。这些优点使得支持向量机应用非常流行，甚至有人认为它是监督学习中最好的定式算法。

      支持向量机试图通过求解一个二次优化问题来最大化分类间隔问题。在JohnPlatt引入SMO之前，训练支持向量常采用非常复杂且低效的二次规划求解，但Platt引入SMO算法之后，每次只需要优化两个alpha值来加快SVM训练速度即可。

     核方法或核技巧可以将数据（有时是非线性数据）从一个低维空间映射到一个高维空间，可以将一个在低维空间中的非线性问题转换成高维空间下的线性问题以便求解。当然，核方法不止适用于SVM，其他算法也开展良好。这里特别说明一句，其中的径向基函数是一个常用的度量两个向量距离的核函数。

     支持向量机是一个二类分类器，面对多类问题时，需要对其进行相应的扩展。SVM的效果对优化参数和所用核函数中的参数也敏感。

       由于笔者知识水平有限，尤其面对支持向量机这类需要大量理论功底的算法，暂且先梳理至此，后续有机会再深入SVM各个环节单独深入研究。若有不足之处，还望各位多多指点。