本文已参与「新人创作礼」活动，一起开启掘金创作之路

在上一篇文章SVM支持向量机与SMO学习笔记(一）中，已经将待求解的未知量由w变成了α，这一节主要介绍求解α的常用算法-SMO算法，序列最小化(Sequence Minimal Optimization）算法，由John Platt于1996年发布。

1. 任务提取

在上一篇中，我们需要求解的目标函数与约束条件为： 通过构造拉格朗日函数，利用KKT条件进行求解出的极值条件为：

$1. \frac{\partial L(w,b,\alpha)}{\partial w}=w-\sum_{i=1}^{m}\alpha_iy_ix_i=0$

$2. \frac{\partial L(w,b,\alpha)}{\partial b}=\sum_{i=1}^m y_i\alpha_i=0$

Tip:这里是累加求和为零,并不一定每个式子为0

约束条件为： $1. \alpha_i*(1-y_i(w^Tx_i + b))=0$ $2.\alpha_i \ge0$ $3.y_i(w^Tx_i + b)-1\ge0$

推导得出其他条件:

$L(w,b,\alpha)=\frac{1}{2}||w||^2+\sum_{i=1}^{m} \alpha_i *(1-y_i(w^Tx_i + b))$

$1.\alpha_i *(1-y_i(w^Tx_i + b))=0$

$2. min_w max_\alpha L(w,b,\alpha)=max_\alpha min_w L(w,b,\alpha) = min_wf(w)=f(w^*)$

原待求解式为$min_w max_\alpha L(w,b,\alpha)$,经对偶转换可以先进行求解$min_w L(w,b,\alpha)$，将极值条件代入到$max_\alpha min_w L(w,b,\alpha)$中，可以得到：

$min_w max_\alpha L(w,b,\alpha)=max_\alpha min_w L(w,b,\alpha) = max_\alpha min_w\frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{m}\alpha_i\alpha_jy_iy_jx_ix_j^T + \sum_{i=1}^{m}(\alpha_i-\alpha_iy_i\sum_{j=1}^{m}\alpha_jy_jx_j^Tx_i-b*\alpha_iy_i) = max_\alpha \sum_{i=1}^{m}\alpha_i - \frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{m}\alpha_i\alpha_jy_iy_jx_ix_j^T$

SMO算法的目标是求解一系列$\alpha$和b来求得上述式子中的最优值。

SMO算法工作原理为，每次选择一对α进行优化处理，一旦找到一对合适的α，就增大一个α同时减小另一个α(这是由限制条件$\sum_{i=1}^m y_i\alpha_i=0$的原因)。这里的合适的α需要满足一定条件，条件之一是两个α必须在间隔边界外（限制条件$0\le\alpha_i\le C$），而第二个条件则是两个α没有进行过区间化处理或者不在边界上。

选一对α进行优化：

$max_\alpha (\alpha_1+\alpha_2) - \frac{1}{2}(y_1^2x_1x_1^T*\alpha_1^2+y_2^2x_2x_2^T * \alpha_2^2 + y_1y_2x_1x_2^T*\alpha_1\alpha_2 + y_1y_2x_1^Tx_2*\alpha_1\alpha_2+y_1x_1\sum_{i=1}^{m-2}y_ix_i^T\alpha_i *\alpha_1 + y_1x_1^T\sum_{i=1}^{m-2}y_ix_i\alpha_i*\alpha_1 + y_2x_2\sum_{i=1}^{m-2}y_ix_i^T\alpha_i * \alpha_2 + y_2x_2^T\sum_{i=1}^{m-2}y_i\alpha_ix_i *\alpha_2)$

由$y_i^2=1$；简化方程表示令$K_{ij} = x_i*x_j^T=x_i^T*x_j$;便于最小化求解，对整个式子取相反数可得:

原式$=min_{\alpha_1,\alpha_2} \frac{1}{2}K_{11}\alpha_1^2+\frac{1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1\alpha_2 + y_1\sum_{i=1}^{m-2}K_{i1}y_i\alpha_i *\alpha_1+y_2\sum_{i=1}^{m-2}K_{i2}y_i\alpha_i * \alpha_2-(\alpha_1+\alpha_2)$

约束条件 $1.\alpha_1y_1 + \alpha2y_2=-\sum_{i=1}^{m-2}\alpha_iy_i = \zeta$

$2.0\le\alpha_i\le C$

现在上述问题转变成为了一个类似于二次函数求极值的问题。

2.问题求解

首先，根据$\alpha_1y_1 + \alpha_2y_2=-\sum_{i=1}^{m-2}\alpha_iy_i = \zeta$得到$\alpha_1=(\zeta-\alpha_2y_2)/y_1$,将该式代入到求解式中进行消元。 消元后，原式$=\frac{1}{2}K_{11}(\zeta-\alpha_2y_2)^2/y_1^2+\frac{1}{2}K_{22}\alpha_2^2+K_{12}y_2\alpha_2(\zeta-\alpha_2y_2)+(\zeta-\alpha_2y_2)\sum_{i=1}^{m-2}K_{i1}y_i\alpha_i+y_2\alpha_2\sum_{i=1}^{m-2}K_{i2}y_i\alpha_i-((\zeta-\alpha_2y_2)/y_1+\alpha_2)$

2.1.极值点求解

首先对$\alpha_2$进行求导，求出其极值点：

$f^{'}(\alpha_2)=K_{11}(\zeta-\alpha_2y_2)+K_{22}\alpha_2+K_{12}y_2\zeta-2K_{12}-y_2\sum_{i=1}^{m-2}K_{i2}y_i\alpha_i-(1-y_2/y_1)$

令$f^{'}(\alpha_2)=0$，可以推出$(K_{11}+K_{22}-2K_{12})\alpha_2=y_2(y_2-y_1+(K_{11}-K_{12})\zeta+\sum_{i=3}^{m}\alpha_iy_i(K_{i1}-K_{i2}))$

由上一篇提到的核函数可以得出$g(x)=\sum_{i=1}^m\alpha_iy_iK(x_i,x)+b$

则$v_i=\sum_{j=3}^m\alpha_jy_jK(x_i,x_j)=g(x_i)-\sum_{i=1}^2\alpha_jy_jK(x_i,x_j)-b$

则上式可变为$(K_{11}+K_{22}-2K_{12})\alpha_2=y_2(y_2-y_1+(K_{11}-K_{12})\zeta+v_1-v_2)$

$=y_2(y_2-y_1+(K_{11}-K_{12})\zeta+g(x_1)-g(x_2)+\alpha_1y_1(K_{12}-K_{11})+\alpha_2y_2(K_{22}-K_{21}))$

再将$\alpha_2^{old}y_2+\alpha_1^{old}y_1=\zeta$代入

$(K_{11}+K_{22}-2K_{12})\alpha_2=\alpha_2^{old}(K_{11}+K_{22}-2K_{12})+y_2(E_1-E_2)$

其中$E_i=g(x_i)-y_i，g(x_i)为估计值,y_i为实际值,E_i为误差$

$综上可得:\alpha_2^{new}=\alpha_2^{old}+y_2(E_1-E_2)/\eta;\eta=(K_{11}+K_{22}-2K_{12})$

2.2.定义域考察

在高中我们就知道，求二次函数极值时，还需要考察极值点是否在定义域内,高中数学老师经常说的就是千万不要忘记定义域！

首先我们的约束条件为:

$1.\alpha_1y_1+\alpha_2y_2=\zeta$

$2.0\le\alpha_i\le C$

其中约束范围可通过画图看出来，图为SMO论文内的图。 图中方形区域就是指$0\le\alpha_1\le C$和$0\le\alpha_2\le C$共同围起来的区域，几条直线就是对应$\alpha_1y_1+\alpha_2y_2=\zeta$的几种不同的情况.

在$\alpha_1y_1+\alpha_2y_2=\zeta$式子中，$y_1$与$y_2$分别可以去-1，1，所以一共可分为4种情况：

（1）$y_1=1，y_2=1$，式子变为$\alpha_1+\alpha_2=\zeta$

图像变为图像中右图所示，其截距为$\zeta$,我们针对$\zeta$不同取值范围来讨论.

①$\zeta\lt0 or \zeta>2C$

直线未落到矩形可行范围内，即可行域无解，为空集。

②$0\le\zeta\le C$

图形对应右图中左下角的直线

$\alpha_2$可行范围为$[0,\zeta]$

③$C\lt\zeta\le2C$

图形对应右图中右上角的直线

$\alpha_2$可行范围为$[\zeta-C,C]$

综上考虑$\zeta$有可行域时，$\zeta$范围可表示为$[max(0,\zeta-C),min(\zeta,C)]即[max(0,\alpha_1+\alpha_2-C),min(\alpha_1+\alpha_2,C)]$

（2）$y_1=-1，y_2=-1$，式子变为$\alpha_1+\alpha_2=-\zeta$ 图像仍为图像中右图所示，其截距为$-\zeta$,我们针对$-\zeta$不同取值范围来讨论. 其情况同(1)类似,所以可以看出图中截距k不一定为$\zeta$ ①$-\zeta\lt0 or -\zeta>2C$

直线未落到矩形可行范围内，即可行域无解，为空集。

②$0\le-\zeta\le C$

图形对应右图中左下角的直线

$\alpha_2$可行范围为$[0,-\zeta]$

③$C\lt-\zeta\le2C$

图形对应右图中右上角的直线

$\alpha_2$可行范围为$[-\zeta-C,C]$

综上考虑$-\zeta$有可行域时，$-\zeta$范围可表示为$[max(0,-\zeta-C),min(-\zeta,C)]即[max(0,\alpha_1+\alpha_2-C),min(\alpha_1+\alpha_2,C)]$

（3）$y_1=1，y_2=-1$，式子变为$\alpha_1-\alpha_2=\zeta$

图形变为左图中的情况，截距为$\zeta$,我们针对$\zeta$不同取值范围来讨论.

①$\zeta\lt-Cor\zeta\gt C$

直线未落到矩形可行范围内，即可行域无解，为空集。

②$0\lt\zeta\le C$

图形对应左图中右下角的直线

$\alpha_2$可行范围为$[0,C-\zeta]$

③$-C\le \zeta \le0$

图形对应左图中左上角的直线

$\alpha_2$可行范围为$[-\zeta,C]$

综上考虑$\zeta$有可行域时，$\zeta$范围可表示为$[max(0,-\zeta),min(C-\zeta,C)]即[max(0,\alpha_2-\alpha_1),min(C+\alpha_2-\alpha_1,C)]$

（4）$y_1=-1，y_2=1$，式子变为$\alpha_1-\alpha_2=-\zeta$

情况同前面类似

综合以上四种情况，我们可以得到， （1）$y_1=y_2时$

$\alpha_2可行范围为[max(0,\alpha_1+\alpha_2-C),min(\alpha_1+\alpha_2,C)]$

$区间端点L=max(0,\alpha_1+\alpha_2-C),H=min(\alpha_1+\alpha_2,C)$

（2）$y_1!=y_2时$

$\alpha_2可行范围为[max(0,\alpha_2-\alpha_1),min(C+\alpha_2-\alpha_1,C)]$

$区间端点L=max(0,\alpha_2-\alpha_1)，H=min(C+\alpha_2-\alpha_1,C)$

2.3.开口方向

我们的目标是上文的$min_{\alpha_1,\alpha_2} \frac{1}{2}K_{11}\alpha_1^2+\frac{1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1\alpha_2 + y_1\sum_{i=1}^{m-2}K_{i1}y_i\alpha_i *\alpha_1+y_2\sum_{i=1}^{m-2}K_{i2}y_i\alpha_i * \alpha_2-(\alpha_1+\alpha_2)$，要求其最小值，还需要判断极值点是否为极小值点。

（1）开口向上：$K_{11}+K_{22}-2K_{12}>0$

判断极值点是否在定义域内，进行参数更新，可以想象一个二次函数

①如果这个极值点在可行域左边，那么我们可以得到这个可行域内二次函数一定在单增，所以此时L应该是那个取最小值的地方。就如大括号的第三种情况

②如果这个极值点在可行域右边，那么此时可行域内一定单减，所以此时H就是那个取最小值的地方，就是大括号里的第一种情况。 （2）开口向下或者退化成为一次函数：$K_{11}+K_{22}-2K_{12}\le0$

这种情况下极小值在端点处取得，可以通过计算两端端点值取最小的。

3.参数更新总结

①求出的可行域

②根据这一项$K_{11}+K_{22}-2K_{12}$和0的关系：

如果$K_{11}+K_{22}-2K_{12}>0$，那么就直接利用之前求好的公式求出新的极值点处的$\alpha_2$，然后看这个$\alpha_2$和可行域的关系，得到$\alpha_2$的新值

如果$K_{11}+K_{22}-2K_{12}\le0$，那么就直接求出H和L对应的边界值，哪个小，就让$\alpha_2$取哪个值

③根据新的$\alpha_2$求出新的$\alpha_1$

④更新ω很简单,利用$w-\sum_{i=1}^{m}\alpha_iy_ix_i=0$求解，现在还需要说一下b的更新：

利用$1-y_i(w^Tx_i + b)=0$进行求解

上式进一步推导,$b=1/y_i-\sum_{j=1}^my_i\alpha_iK(x_i,x_j)$

$b_1^{new}=y_1-y_1\alpha_1^{new}K_{11}-y_2\alpha_2^{new}K_{12}-\sum_{i=3}^my_i\alpha_iK_{i1}$

$E_1=g(x_1)-y_1=\alpha_1^{old}y_1K_{11}+\alpha_2^{old}y_2K_{12}+\sum_{i=3}^my_i\alpha_iK_{i1}+b_1^{old} -y_1$

消元后可得:

$b_1^{new}=b_1^{old}+(\alpha_1^{old}-\alpha_1^{new})y_1K_{11}+(\alpha_2^{old}-\alpha_2^{new})y_2K_{12}-E_1$

4.更新参数$\alpha$的选取

4.1.寻找第一个待更新$\alpha$

①遍历一遍整个数据集，对每个不满足KKT条件的参数，选作第一个待修改参数

②在上面对整个数据集遍历一遍后，选择那些参数满足$0\lt\alpha\lt C$的子集，开始遍历，如果发现一个不满足KKT条件的，作为第一个待修改参数，然后找到第二个待修改的参数并修改，修改完后，重新开始遍历这个子集

③遍历完子集后，重新开始①②，直到在执行①和②时没有任何修改就结束

4.2.寻找第二个待更新$\alpha$

①启发式找，找能让$|E_1-E_2|$最大的

②寻找一个随机位置的满足$0\lt\alpha\lt C$的可以优化的参数进行修改

③在整个数据集上寻找一个随机位置的可以优化的参数进行修改

④都不行那就找下一个第一个参数

ri = Eiyi 关于寻找不满足KKT条件第一个α认识（有疑问），伪码中 if ((r2<-tol) and (α<C))or((r2>tol)and(α>0)) tol为容错率

KKT条件为：

$1.w-\sum_{i=1}^{m}\alpha_iy_ix_i=0$

$2.\sum_{i=1}^m y_i\alpha_i=0$

$3.y_i(w^Tx_i + b)-1\ge0$

$4.0\le\alpha\le C$

$5.\alpha_i*(1-y_i(w^Tx_i + b))=0$

$r_i=E_iy_i=(g(x_i)-y_i)y_i=y_ig(x_i)-1$

$r_i<-tol=>y_ig(x_i)<1-tol<1$又KKT条件中$y_i(w^Tx_i + b)\ge1$

$r_i>tol=>y_2g(x_i)>1+tol>1$

以下为部分在机器学习实战中的代码，加深算法理解

import numpy as np
import random

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

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

def svmSimple(dataMatIn,classLabels,C,toler,maxIter):
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()
    b = 0
    m,n = np.shape(dataMatIn)
    alphas = np.zeros((m,1))
    iter = 0
    while iter < maxIter:
        alphaPairsChanged = 0
        for i in range(m):
            fXi = float(np.multiply(alphas,labelMat).T * dataMatrix * dataMatrix[i,:].T) + b
            Ei = fXi - float(labelMat[i])
            if ((labelMat[i]*Ei > toler) and (alphas[i] > 0)) or ((labelMat[i] * Ei < -toler) and (alphas[i] < C)): #误差超过容错率
                j = selectJrand(i,m)
                fXj = float(np.multiply(alphas,labelMat).T * dataMatrix * dataMatrix[j,:].T) + b
                Ej = fXj - 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 + alpha[j] - alpha[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 #减号因为eta符号与计算是相反的
                alphas[j] = clipAlpha(alphas[j],H,L)
                if abs(alphas[j] - alphaJold) < 0.00001:
                    print("j not moving enough")
                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 (alphas[i] > 0) and (alphas[i] < C):
                    b = b1
                elif (alphas[j] > 0) and (alphas[j] < C):
                    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

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
        self.eCache = np.mat(np.zeros((self.m,2)))

def calcEk(os, k):
    fXk = float(np.multiply(os.alphas,os.labelMat[k,:]).T * os.X * os.X[k,:].T) + 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 = np.nonzero(os.eCache[:,0].A)[0] #.A将mat转成array,返回第一维为0的可用eCache
    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 updateEk(os,k):
    Ek = calcEk(os,k)
    os.eCache[k] = [1,Ek]

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

def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin',0)):
    os = optStruct(np.mat(dataMatIn),np.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))
        else:
            nonBounfIs = np.nonzero((os.alphas.A > 0)*(os.alphas.A < C))[0]
            for i in nonBounfIs:
                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

    return os.b,os.alphas

参考链接：SMO算法的个人理解