SVM(Support Vector Machine, 支持向量机)有很多实现,但是本章其中最流行的一种实现,即序列最小优化(Sequence Minimal Optimization, SMO)算法。
最小优化(Sequence Minimal Optimization, SMO)算法 一种求解支持向量机二次规划的算法。

6.1 基于最大间隔分隔数据

支持向量机
优点:泛化错误率低,计算开销不大,结果易解释。
缺点:对参数调节和核函数的选择敏感,原始分类器不加修改仅使用于处理二类问题。
使用数据类型:数值型和标称型数据。

线性可分(linearly separable):
两类点被一条直线完全分开
数学定义:
D 0 D_0 D0 D 1 D_1 D1 n n n维度欧氏空间中的两个点集。如果存在 n n n维向量 w ⃗ \vec w w 和实数 b b b,使得所有属于 D 0 D_0 D0的点 x i x_i xi都有 w x i + b > 0 wx_i + b > 0 wxi+b>0,而对于所有属于 D 1 D_1 D1的点 x j x_j xj则有 w x i + b < 0 wx_i + b < 0 wxi+b<0,则我们称 D 0 D_0 D0 D 1 D_1 D1线性可分。

分割超平面(seqarating hyperplane):
数据集分割开来的直线

超平面(hyperplane):
n维欧氏空间中余维度等于一的线性子空间,也就是必须是(n-1)维度。也就是分类的决策边界。

间隔(margin):
点到分割面的距离。

支持向量(support vector):
离分割超平面最近的那些点。

6.2 寻找最大间隔

分隔超平面的形式可以写成 w T x + b \textbf {w}^T \textbf {x} + b wTx+b
要计算A到分割超平面的距离,就必须给出点到分隔面的法线或垂线的长度,该值为 ∣ w T x + b ∣ / ∥ w ∥ |\textbf {w}^T \textbf {x} + b| / \Vert{\textbf {w}\Vert} wTx+b∣/∥w

6.2.1 分类器求解的优化问题

6.2.2 SVM应用的一般框架

SVM的一般流程:

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

6.3 SMO高效优化算法

6.3.1 Platt的SMO算法

SMO用于训练SVM。SMO表示序列最小优化(Sequence Minimal Optimization)。
SMO算法是将最大优化问题分解为多个小优化问题来求解的。
SMO算法的工作原理是:每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha,那么就增大其中一个同时减小另一个。
这里所谓的‘合适’就是指两个alpha必须要符合一定的条件,条件之一就是这两个alpha必须要在间隔边界之外,而其第二个条件则是这两alpha还没进行过区间化处理或者不在边界上。

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

SMO算法中的辅助函数:

    import random


    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 selectJrand(i, m):
        '''

        :param i: alpha的下标
        :param m: 所有alpha的数目
        :return:
        '''
        j = i
        while(j == i):
            j = int(random.uniform(0, m))
        return j

    def clipAlpha(aj, H, L):
        '''
        用于调整大于H或小于L的alpha值
        :param aj:
        :param H:
        :param L:
        :return:
        '''
        if aj > H:
            aj = H
        if L > aj:
            aj = L
        return aj

    if __name__ == '__main__':
        dataArr, labelArr = loadDataSet('./resource/testSet.txt')
        print(dataArr)
        print(labelArr)

该SMO函数的伪代码大致如下:

创建一个alpha向量并将其初始化为0向量
当迭代次数小于最大迭代次数时(外循环)
    对数据集中的每个数据向量(内循环):
    如果该数据向量可以被优化:
        随机选择另外一个数据向量
        同时优化这两个向量
        如果两个向量都不能被优化,退出内循环
如果所有向量都没被优化,增加迭代数目,继续下一次循环.

简化版SMO算法:

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    '''

    :param dataMatIn: 数据集
    :param classLabels: 类别标签
    :param C: 常数C
    :param toler: 容错率
    :param maxIter: 退出前最大的循环次数
    :return:
    '''
    ## 将输入转换成NumPy矩阵
    dataMatrix = mat(dataMatIn);
    labelMat = mat(classLabels).transpose()
    b = 0;
    m, n = shape(dataMatrix)
    alphas = mat(zeros((m, 1)))
    iter = 0
    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 checks if an example violates KKT conditions
            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])  # update i by the same amount as j
                # the update is in the oppostie direction
                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


if __name__ == '__main__':
    dataArr, labelArr = loadDataSet('./resource/testSet.txt')
    b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
    print(b)
    print(alphas[alphas > 0])
    print(shape(alphas[alphas > 0]))

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

完整版Platt SMO的支持函数:

class optStructK:
    def __init__(self,dataMatIn, classLabels, C, toler):  # Initialize the structure with the parameters 
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag


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

def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            if k == i: continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:   #in this case (first time around) we don't have any valid eCache values
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

完整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) #this has been changed from selectJrand
        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.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
        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) #added this for the Ecache
        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])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        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

完整版Platt SMO的外循环代码:

def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):  # full Platt SMO
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
    iter = 0
    entireSet = True;
    alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:  # go over all
            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:  # go over non-bound (railed) alphas
            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  # toggle entire set loop
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas

if __name__ == '__main__':
    dataArr, labelArr = loadDataSet('./resource/testSet.txt')
    b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 0.4)
    print(b)
    print(alphas[alphas > 0])
    print(shape(alphas[alphas > 0]))

实现分类:

def calcWs(alphas, dataArr, classLabels):
    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


if __name__ == '__main__':
    dataArr, labelArr = loadDataSet('./resource/testSet.txt')
    b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 0.4)
    ws = calcWs(alphas, dataArr, labelArr)
    dataMat = mat(dataArr)
    print(ws)
    print(dataMat[0] * mat(ws) + b)

6.5 在复杂数据上应用核函数

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

从某个特征空间到另外一个特征空间的映射是通过核函数来实现的。
核函数想象成一个包装器(wrapper)或者接口(interface),它能把数据从某个很难处理的形式转换成为另一个较容易的处理的形式。

6.5.2 径向基核函数

径向基函数:一个采用向量作为自变量的函数,能够基于向量距离运算输出一个标量。
径向基函数的高斯版本:
k ( x , y ) = e x p ( − ∥ x − y ∥ 2 2 α 2 ) k(x, y) = exp(\frac {-\Vert x- y \Vert^2} {2 \alpha ^2}) k(x,y)=exp(2α2xy2)

其中, α \alpha α是用户定义的用于确定到达率(reach)或者说函数值跌落到0的速度参数。

核转换函数:

def kernelTrans(X, A, kTup):  # calc the kernel or transform data to a higher dimensional space
    m, n = shape(X)
    K = mat(zeros((m, 1)))
    if kTup[0] == 'lin':
        K = X * A.T  # linear kernel
    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))  # divide in NumPy is element-wise not matrix like Matlab
    else:
        raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')
    return K

class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters 
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

6.5.3在测试中使用核函数:

利用核函数进行分类的径向基测试函数:

def testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('./resource/testSetRBF.txt')
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd] #get matrix of only support vectors
    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('./resource/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))

6.6 示例:手写识别问题回顾

示例:基于SVM的数字识别:

  1. 收集数据:提供的文本。
  2. 准备数据:基于二值图像构造向量。
  3. 分析数据:对图像向量进行目测。
  4. 训练算法:采用两种不同的核函数,并对径向基核函数采用不同的的设置来运行SMO算法。
  5. 测试算法:编写一个测试函数来测试不同的核函数并计算错误率。
  6. 使用算法:一个图像识别的完整应用含需要一些图像处理的知识,这里并不打算深入介绍。

基于SVM的手写数字识别:

def loadImages(dirName):
    from os import listdir
    hwLabels = []
    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)):
    dataArr,labelArr = loadImages('./resource/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))
    dataArr,labelArr = loadImages('./resource/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))

if __name__ == '__main__':
    testDigits(('rbf', 20))

6.7 本章小结

点赞(0) 打赏

评论列表 共有 0 条评论

暂无评论

微信公众账号

微信扫一扫加关注

发表
评论
返回
顶部