机器学习实战【7】（SMO算法实现）

本博客记录《机器学习实战》（MachineLearningInAction）的学习过程，包括算法介绍和python实现。

SMO算法

前两篇文章介绍了SVM的原理，经过一番推导，原始问题转化为：

minα∑i=1n12∑i,j=1nαiαjyiyjK(xi,xj)−αis.t.,0≤αi≤C,i=1,...,n∑i=1naiyi=0

SMO算法就是用来解决这个问题，求解出这些α之后，超平面的参数就可以通过这些α计算出来。

二元优化

SMO算法的核心在于，它把原本n个参数的优化问题拆分为很多个小的子问题，每次只优化其中的两个参数而固定其它参数，两个参数的优化是很快的，从而使得最终算法的效率非常高。
假设在n个参数(α1,...,αn)中选取α1和α1进行优化，其它参数全部视为常数，则原问题化简为：

minΨ(α1,α2)=12K11α21+12K22α22+K12y1y2α1α2−(α1+α2)+y1v1α1+y2v2α2

其中：

vi=∑j=3nαjyjKij=f(xi)−∑j=12αjyjKij−b,i=1,2

由等式约束条件得出：

α1y1+α2y2=−∑i=3nαiyi=ζ

以此消去上式中的α1，得到：

minΨ(α2)=12K11(ζ−α2y2)+12K22α22+K12y2(ζ−α2y2)α2−α2−(ζ−α2y2)y1+v1(ζ−α2y2)+y2v2α2

上式求导：

∂Ψ(α2)∂α2=(K11+K22−2K12)α2+K12y2ζ−K11y2ζ+y1y2−1−v1y2+v2y2=0

带入ζ=αold1y1+αold2y2化简得到：

αnew,unclipped2=αold2+y2(E2−E1)η

其中，E表示误差：

Ei=f(xi)−yi
η=K11+K22−2K12

解的修剪

上式中计算出的α2有可能在边界(0-C)之外，所以需要对结果进行修剪(clip)，可以计算出两种情况下α2的上下界：
y1≠y2时：L=max(0,αold2−αold1)H=min(C,C+αold2−αold1)
y1=y2时：L=max(0,αold2+αold1−C)H=min(C,αold2−αold1)
依此对α2进行修剪之后，就可以进一步计算出α1。

b的更新

因为在每次迭代的过程中需要计算当前模型的预测函数值f(x)，所以需要计算出b以便下次迭代使用，更新时有几种情况。
若0<α1<C，则有y1(wTx1+b)=1，可以推导出：
b1=bold−E1−y1K11(αnew1−αold1)−y2K21(αnew2−αold2)
同样，若0<α1<C，则：
b2=bold−E2−y1K12(αnew1−αold1)−y2K22(αnew2−αold2)
若两者都不满足，则取b=b1+b22。

参数的启发式选择

第一个参数选择那些违反kkt条件的参数，第二个参数选择那些使得更新跨度更大的参数。

第一个参数的确定（外循环）：

kkt条件：

αi=0⇒yiui≥1
αi=C⇒yiui≤1
0<αi<C⇒yiui=1

在外循环中，遍历一次整个样本集，遇到违反kkt条件的α就选作第一个参数进行更新，结束之后多次遍历非边界样本集（满足0<α<C的那些样本），直至边界样本集全部满足kkt条件。此时，再次遍历一次整个样本集，再遍历非边界样本集，以此往复，直至整个样本集都满足kkt条件，此时算法结束。

第二个参数的确定（内循环）：

第二个参数依据更新跨度大小来选择，而由于η计算量较大，所以近似以|E2−E1|来衡量更新跨度，即在确定了第一个参数的情况下，选择第二个参数使得|E2−E1|最大。
第二个参数的选择首先在非边界样本集上进行，若无法找到具有足够更新跨度的参数，则在整体样本集上找，若还是找不到则重新选择第一个参数。

python实现

# 定义结构用于存放参数和中间变量class optStruct:    def __init__(self,dataMatIn, classLabels, C, toler, kTup):        # 输入数据集        self.X = dataMatIn        # 输入类标签        self.labelMat = classLabels        # 参数C，用于调整松弛变量比重        self.C = C        # 作为kkt条件判断时的阈值        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)))        # 核函数计算结果矩阵        self.K = mat(zeros((self.m,self.m)))        for i in range(self.m):            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)# 计算误差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):    maxK = -1; maxDeltaE = 0; Ej = 0    oS.eCache[i] = [1,Ei]    # 误差计算之后会进行缓存，这里先在缓存中遍历误差，选择最优的    validEcacheList = 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# 更新参数，返回参数是否被更新def innerL(i, oS):    Ei = calcEk(oS, i)    # 此处违反kkt条件的判断采用阈值的方式，在边界附近某个小范围内的值都视作在边界上    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.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) #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])        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# smo算法，部分细节与理论不太相同，但影响不大def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup = ('lin', 0)):    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:            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

总结

smo算法一次选择两个参数进行更新，将复杂的优化问题分解为一个个很好解决的小问题进行化解。机器学习实战书中的实现代码相比paper中的理论进行了一些简化，但对整体速度影响不大，但是其中的误差缓存方式貌似有些问题，在更新了参数之后只对对应的数据误差进行了更新，但此时其它的误差应该是失效了，却没有重新计算，可能是因为其它的误差变化不大，所以不会有很大的影响吧。