机器学习实战之四——支持向量机：完整Platt SMO算法

上一篇文章中介绍了简化版的SMO算法及其应用，这篇文章将会介绍完整的Platt SMO算法。完整版与简化版最大的区别在于，简化版中选择要改变的第二个参数alpha[j]时，是随机选择的；而在完整版中，alpha[j]的选取原则是使得对第j个数据的预测误差Ej和第i个数据预测误差Ei二者之间的差值最大的alpha[j]。下面作者将《机器学习实战》一书的主要代码抄写过来，并在必要的地方作了注释，以期各位读者能够看懂。

import numpy as np
import copy

def selectJrand(i,m):
    j = i
    while j == i:
        j = int(np.random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

#通过创建一个对象来建立一个数据结构，以保存重要的变量，其目的是为了在一定程度省去写代码的麻烦。
class OptStruct:
    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)))
        self.b = 0
        #eCache的第二列是每条数据的预测值与实际值的误差；第一列是一个标识列，当该列的值为1时， 表明已经已经根据其对应的alpha计算出了Ej，可以直接拿来使用，因此节省了大量时间。
        self.eCache = np.mat(np.zeros((self.m,2)))

#这里的函数就是用来计算预测误差的，由于简化版SMO中很多地方都用到了这个计算过程，因此此处将其封装起来，可以直接调用。        
def calcEk(oS,k):
    w = np.sum(np.multiply(np.multiply(oS.alphas,oS.X),oS.labelMat),axis=0)
    fXk = w*oS.X[k].T + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek

#这个函数是通过前面说的启发式方法来选择j；注意到当第一次使用这个函数时，所有的alpha都为0，即都是无效的，因此此时就使用selectJrand()方法来随机选择。
def selectJ(i,oS,Ei):
    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 = 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

#对于修改了的第k个alpha[k]，要在缓冲区中将其变为有效的
def updateEk(oS,k):
    Ek = calcEk(oS,k)
    oS.eCache[k] = [1,Ek]

#这个函数是内循环，用来选择alpha[j]，以及计算和修改参数    
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 = copy.deepcopy(oS.alphas[i])
        alphaJold = copy.deepcopy(oS.alphas[j])
        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
        K12 = oS.X[i,:]*oS.X[j,:].T
        K11 = oS.X[i,:]*oS.X[i,:].T
        K22 = oS.X[j,:]*oS.X[j,:].T
        eta = 2.0 * K12 - K11 - K22
        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,i)
        b1 = oS.b - Ei - oS.labelMat[i]*(oS.alphas[i]-alphaIold)*K11 - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*K12
        b2 = oS.b - Ej - oS.labelMat[i]*(oS.alphas[i]-alphaIold)*K12 - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*K22
        if (oS.alphas[i] > 0) and (oS.C > oS.alphas[i]):
            oS.b = b1
        elif (oS.alphas[j] > 0) and (oS.C > oS.alphas[j]):
            oS.b = b2
        else:
            oS.b = (b1 + b2)/2.0
        return 1
    else:
        return 0

#完整算法的主体程序，该程序在外循环中选择第一个alpha，也就是alpha[i]时，会在两种方式之间交替进行：一种方式是在所有数据集上进行单边扫描，即下面的entireSet；另一种是在非边界alpha中实现单边扫描，即下面的nonBoundIs
def smoP(dataMatIn,classLabels,C,toler,maxIter,kTup=('lin',0)):
    oS = OptStruct(np.mat(dataMatIn),np.mat(classLabels).tranpose(),C,toler)
    iter_ = 0
    entireSet = True
    alphaPairsChanged = 0
    #alphaPairsChanged用来标识是否有alpha发生了改变。显然，程序停止运行的条件除了达到最大迭代次数外，还有一种情况：遍历了所有的数据集也没有alpha发生改变。
    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用来寻找不满足KKT条件的alpha
            nonBoundIs = np.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

该算法用在点数据集上效果如图：

测试成功率为0.9.

在乳腺癌数据集上测试成功率为0.7