机器学习笔记3：支持向量机的SMO高效优化算法

1. SMO算法背景

　　根据上一篇的支持向量机SVM的基础知识介绍，算法的根本目标就是实现对优化目标函数在约束条件下的求解问题。SMO算法表示序列最小优化（Sequential Minimal Optimization），将大优化问题分解为多个小优化问题来求解，小优化问题便于求解，同时对它们进行顺序求解的结果与作为整体求解的结果完全一致，并且SMO算法的求解时间要短很多。
　　SMO算法的目标就是求出一系列的alpha与b，一旦求出这些alpha，就很容易计算出权重向量w并得到分隔超平面。
　　SMO算法的工作原理是：每次循环中选择两个alpha，增大其中的一个同时减少另外一个，所谓的合适就是两个alpha符合一定的条件，条件之一就是两个alpha必须要在间隔边界之外，而第二个条件则是这两个alpha还没有进行过区间化处理或者不在边界上。

2.简化版SMO算法处理小规模数据集

　　简化版的SMO算法虽然执行速度较慢，Platt SMO算法中的外循环确定要优化的最佳alpha对，简化版跳过了这个部分，简化版的SMO算法首先在数据集上遍历每一个alpha，然后在剩下的alpha集合中随机选择另一个alpha，从而建立alpha对。这一点非常重要，同时改变两个alpha，因为我们有一个约束条件：

　　由于改变一个alpha可能会导致该约束条件失效，因此我们同时改变两个alpha。因此，代码整体如下，部分重要部分进行注释。

# coding=utf-8

'''
Date: August 14 2017
Author：Chenhao
Description: simplified SMO Algorithm
'''

from numpy import *
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

#读取文档内数据，数据格式为三列，前两列包含二维数据，第三列为label
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

#用于在区间范围内选择一个整数
#输入函数i代表第一个alpha的下标，m是所有alpha的数目
def selectJrand(i,m):
    j=i             #we want to select any J not equal to i
    while (j==i):
        j = int(random.uniform(0,m))
    return j

#用于在数值太大的情况下进行调整
#用于调整大于H小于L的alpha值
def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

#简化版的SMO算法
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):     #五个输入参数对应表示数据集，类别标签，常数C，容错率，退出前最大循环次数
    dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()    #将输入数据转化为matrix的形式，将类别标签进行转置
    b = 0; m,n = shape(dataMatrix)    #将常数b进行初始化，m，n分别表示数据集的行和列
    alphas = mat(zeros((m,1)))      #将alpha初始化为m行1列的全零矩阵
    iter = 0        #该变量用于存储是在没有任何alpha改变的情况下遍历数据集的次数
    while (iter < maxIter):     #当遍历数据集的次数超过设定的最大次数后退出循环
        alphaPairsChanged = 0       #该参数用于记录alpha是否已进行优化，每次循环先将其设为0
        for i in range(m):      #将数据集中的每一行进行测试
            #multiply是numpy中的运算形式，将alpha与标签相乘的结果的转置与对应行的数据相乘再加b
            #fxi表示将该点带入后进行计算得到的预测的类别，若Ei=0，则该点就在回归线上
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions
            #label[i]*Ei代表误差，如果误差很大可以根据对应的alpha值进行优化，同时检查alpha值使其不能等于0或C
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i,m)       #随机选取另外一个数作为j，也就是选取另外一行数据
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])       #采用与上述相同的方式来计算j行的预测误差
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();     #暂存变量，便于后面对其进行比对
                #调整H与L的数值，便于用来控制alpha的范围
                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     #continue在python中代表本次循环结束直接运行下一次循环
                #eta代表的是alphas[j]的最优修改量，如果eta大于等于0则需要退出for循环的当前迭代过程
                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]赋新值，同时控制alphas[j]的范围在H与L之间
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                #然后检查alphas[j]是否有轻微改变，如改变很小，就进行下一次循环
                if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue
                #将alpha[i]与alpha[j]同时进行改变，改变的数值保持一致，但是改变的方向相反
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                #在完成alpha[i]与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]): b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0
                alphaPairsChanged += 1      #该参数用于记录alpha是否已进行优化
                print "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)      #iter表示遍历的次数
        if (alphaPairsChanged == 0): iter += 1
        else: iter = 0
        print "iteration number: %d" % iter
    return b,alphas

def plot():
    dataplot = array(dataArr)
    n = shape(dataplot)[0]
    plotx1 = []; ploty1 = []
    plotx2 = []; ploty2 = []
    for i in range(n):
        if int(labelArr[i]) == 1:
            plotx1.append(dataplot[i,0])
            ploty1.append(dataplot[i,1])
        else:
            plotx2.append(dataplot[i, 0])
            ploty2.append(dataplot[i, 1])
    fig = plt.figure()
    ax =fig.add_subplot(111)
    ax.scatter(plotx1,ploty1,c='red')
    ax.scatter(plotx2,ploty2,c='blue')
    plt.show()


dataArr,labelArr = loadDataSet('testSet.txt')

b,alphas = smoSimple(dataArr,labelArr,0.6,0.001,40)

for i in range(100):
    if alphas[i]>0.0:
        print dataArr[i],labelArr[i]

plot()

　　代码运行结果如下：
在alpha大于零的前提下的数据及其label：
[4.658191, 3.507396] -1.0
[3.457096, -0.082216] -1.0
[2.893743, -1.643468] -1.0
[6.080573, 0.418886] 1.0

3.利用完整的Platt SMO算法加速优化

　　上述的简化版SMO算法在几百个点的小规模数据集上运行是没有问题的，但是在更大的数据集上运行就会导致运行速度变慢。完整版的Platt SMO算法在实现alpha的更改和代数运算的优化环节与简化版的SMO算法一模一样，唯一不同的就是选择alpha方式。
　　Platt SMO算法是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间交替：一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界alpha中实现单遍扫描，所谓非边界alpha指的是那些不等于边界0或C的alpha值。对整个数据集的扫描相当容易，二实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后对这个表进行遍历。同时，该步骤会跳过那些已知不变的alpha值。
　　在选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha值，在优化过程中，会通过最大化步长的方式来获得第二个alpha值，在简化版的SMO算法中，我们会在选择j之后计算错误率Ej。但是在这里，我们会建立一个全局的缓存用于保存误差值，并从中选择步长或者说是Ei-Ej最大的alpha值。
　　下面是在使用Platt SMO算法前需要准备的用于清理代码的数据结构和3个用于对于E进行缓存的辅助函数。前面的数据结构以及辅助函数如下：

# coding=utf-8
#建立一个全局的数据结构用来保存所有重要的数值
#将数据转移到一个结构来实现，省去手工输入麻烦
class optStruct:
    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      #最大循环次数
        #m表示数据的行数
        self.m = shape(dataMatIn)[0]
        #初始化为m行1列的全零矩阵
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        #eCache的第一列表示是否有效的标志位，第二列是实际E值
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag

#该函数用来计算E值并返回
#k可以表示i或者j的数值
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T * (oS.X*oS.X[k,:].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])    #数据点代入alpha和b后的值与真实标签的差值表示Ek
    return Ek


#该函数用于选择第二个alpha值或者说是内循环的alpha值，保证每次循环采用最大步长
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]  # 将eCache的第一列数据也就是标志位设置为有效，接下来选择合适的alpha使得步长也就是Ei-Ej最大
    #构建出一个非零表
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]     #将eCache的标志位到处存入到validECacheList中
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:  #通过eCache的有效标志位作为判断依据进行循环找到满足步长最大化的alpha
            if k == i: continue     #如果k等于i，跳出循环，避免浪费时间
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k;
                maxDeltaE = deltaE;
                Ej = Ek
        return maxK, Ej
    else:   #这种情况是在第一次循环过程中eCache没有有效值，就采用遍历数据集的方式
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

#将计算误差值缓存到alpha值进行优化之后会用到
def updateEk(oS, k):  # after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]

　　各个函数及其解释如上代码片，上面的全局变量的定义及其三个辅助函数的构建目的是为了便于Platt SMO算法的高效执行，下面是完整的Platt SMO算法及其注释。

#Platt SMO算法的内循环部分，确定第二个alpha的值
def innerL(i, oS):
    Ei = calcEk(oS, i)      #计算第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)):
        #selectJ函数用来选择第二个alpha值来保证步长得到最大
        j,Ej = selectJ(i, oS, Ei)
        #将alphai与alphaj的旧值存储
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        # 调整H与L的数值，便于用来控制alpha的范围
        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     #如果L=H时，则返回0
        # eta代表的是alphas[j]的最优修改量，如果eta大于等于0则需要退出for循环的当前迭代过程
        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
        # 计算alphaj的数值，给alphas[j]赋新值，同时控制alphas[j]的范围在H与L之间
        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
        # 然后检查alphas[j]是否有轻微改变，如改变很小，就进行下一次循环
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; return 0
        # 将alpha[i]与alpha[j]同时进行改变，改变的数值保持一致，但是改变的方向相反
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        #将Ek存入缓存
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        # 在完成alpha[i]与alpha[j]的优化之后，给这两个alpha设置一个常数项b
        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

#选择第一个Alpha值的外循环
#五个输入参数对应表示数据集，类别标签，常数C，容错率，退出前最大循环次数
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
    iter = 0        #该变量用于存储是在没有任何alpha改变的情况下遍历数据集的次数
    entireSet = True
    # 该参数用于记录alpha是否已进行优化，每次循环先将其设为0
    alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   #go over all
            for i in range(oS.m):       #开始的for循环在数据集上遍历任意可能的alpha
                alphaPairsChanged += innerL(i,oS)       #用innerL来选择第二个alpha
                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)       #遍历非边界值来确定第二个alpha
                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

　　部分运行结果如下：
fullSet, iter: 0 i:0, pairs changed 1
fullSet, iter: 0 i:1, pairs changed 1
fullSet, iter: 0 i:2, pairs changed 1
fullSet, iter: 0 i:3, pairs changed 1
fullSet, iter: 0 i:4, pairs changed 1
…
…（省略一部分）
fullSet, iter: 2 i:94, pairs changed 0
fullSet, iter: 2 i:95, pairs changed 0
fullSet, iter: 2 i:96, pairs changed 0
fullSet, iter: 2 i:97, pairs changed 0
fullSet, iter: 2 i:98, pairs changed 0
fullSet, iter: 2 i:99, pairs changed 0
iteration number: 3

　　在完成了alpha的计算过程之后，需要通过alpha来计算得到w才可以得到分隔超平面。下面的函数就实现了w的计算。

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

　　我们可以添加如下代码，来用第一行数据作为验证。

dataArr,labelArr = loadDataSet('testSet.txt')
b,alphas = smoP(dataArr, labelArr , 0.6 , 0.001 , 40 )
ws = calcWs(alphas,dataArr,labelArr)
print ws
datMat = mat(dataArr)
test = datMat[0]*mat(ws)+b
print test

　　结果如下所示：

ws = [[ 0.65307162]
* [-0.17196128]]*
test = [[-0.92555695]]
　　
　　第一个label是-1，因此符合验证的结果。进一步的，可以添加代码，将数据及其分隔平面展示出来。

def plot():
    dataplot = array(dataArr)
    n = shape(dataplot)[0]
    plotx1 = []; ploty1 = []
    plotx2 = []; ploty2 = []
    for i in range(n):
        if int(labelArr[i]) == 1:
            plotx1.append(dataplot[i,0])
            ploty1.append(dataplot[i,1])
        else:
            plotx2.append(dataplot[i, 0])
            ploty2.append(dataplot[i, 1])
    fig = plt.figure()
    ax =fig.add_subplot(111)
    ax.scatter(plotx1,ploty1,c='red')
    ax.scatter(plotx2,ploty2,c='blue')
    plotx = arange(0.0 , 10.0 , 0.1)
    ploty = -(ws_float0*plotx+ b_float)/ws_float1
    ax.plot(plotx,ploty)
    plt.show()

dataArr,labelArr = loadDataSet('testSet.txt')
b,alphas = smoP(dataArr, labelArr , 0.6 , 0.001 , 40 )
ws = calcWs(alphas,dataArr,labelArr)
print ws
datMat = mat(dataArr)
b_float = float(b[0])
ws_float0 = float(ws[0])
ws_float1 = float(ws[1])
test = datMat[1]*mat(ws)+b
print test
plot()

　　最终分隔的结果如下图所示：

　　项目的完整代码以及数据源已经上传到我的github上，其中SMO.py是简化版的SMO算法，Platt SMO.py是SMO完整版的算法，就实验结果来看，二者算法的运算速度差距很大，结果在这个实验中差距不大，后面将会用SVM中的算法与核函数原理进行数据预测来验证这个算法的有效性。
　　
　　项目数据与源代码

August 17 , 2017