Platt SMO算法

Platt SMO算法
1996年，John Platt 发布了一个称为SMO的强大算法，用于训练SVM
SMO表示序列最小优化（SequentialMinimal Optimization）。Platt 的 SMO 算法是将大优化问题分解为多个小优化问题来求解的。这些小优化问题往往容易求解，并且对它们进行顺序求解的结果与它们作为整体来求解的结果是完全一致的。在结果完全相同的同时，SMO算法的求解时间短很多。SMO算法的目标是求出一系列的alpha和b，一旦求出了这些alpha，就很容易计算出权重w病得到分隔超平面。
SMO算法的工作原理是：每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha，那么就增大其中一个同时减小另一个。这里所谓的“合适”就是指两个alpha必须要符合一定的条件，条件之一就是这两个alpha必须要在间隔边界之外，而其第二个条件则是这两个alpha还没有进行过区间化处理或者不在边界上。

接下来进行SMO算法简化版本描述
Platt SMO 算法中的外循环确定要优化的最佳alpha对，而简化版会跳过这一部分。
首先，遍历数据集上每一个alpha，然后在剩下的alpha集合中随机选择另一个alpha，从而构建alpha对。这里有一点相当重要，就是我们要同时改变两个alpha。之所以这样做是因为我没有一个约束条件：

由于改变一个alpha可能会导致该约束条件失效，因此我们总是改变两个alpha。

相关全部源代码下载地址：
代码

SMO函数的伪代码如下：

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

简化版SMO算法，下面用python代码呈现

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    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]):#保证alpha在0与C之间 
                    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

这个函数比较大，该函数有5个输入参数，分别是：数据集、类别标签、常数C、容错率和取消前最大的循环次数。
每次循环中，将alphaPairsChanged先设为0，然后在对整个集合顺序遍历。变量alphaPairsChanged用于记录alpha是否已经进行优化。当然，在循环结束时就会得知这一点。
fXi能够计算出来，就是我们要预测的类别。然后，基于这个实例的预测结果和真是结果的比对，就可以计算误差Ei。如果误差很大，那么可以对该数据实例所对应的alpha值进行优化。
在if语句中，不管是正间隔还是负间隔都会被测试。并且在该if语句中，也要同时检查alpha值，以保证其不能等于0或者C。由于后面alpha小于0或者大于C时将被调整为O或C，所以一旦在该if语句中它们等于这两个值的话，那么表明它们就在“边界”上了，因此不再对它们进行优化。

接下来需要用辅助函数来随机选择第二个alpha值，即alpha[j]。计算方法同alpha[i]

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

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

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

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

有一点需要注意：Python会通过引用的方式传递所有列表，所以必须明确地告知Python要为alphaIold和alphaJold分配新的内存；否则的话，对新值和旧值进行比较时，我们就看不到新旧值的变化。之后我们开始计算L和H，他们用于将alpha[j]调整到O到C之间。

Eta是alpha[j]的最优修改量，如果eta为0，那就是说需要退出for循环的当前迭代过程。该过程对真实SMO算法进行了简化处理。如果eta为O，那么计算新的alpha[j]就比较麻烦了，这里我们就不对此进行详细的介绍。有需要的渎职可以阅读Platt的原文了解更多细节。现实中，这种情况并不常发生。

然后，就需要检查alpha[j]是否有轻微改变。如果是的话，就退出for循环。然后，alpha[i]和alpha[j]同样进行改变，虽然改变的大小一样，但是改变的方向正好相反（即如果一个增加，那么另一个减少）。在对alpha[i]和alpha[j]进行优化之后，给这两个alpha值设置一个常数项b.

最后，在优化过程结束的同时，必须确保在合适的时机结束循环。如果程序执行到for循环的最后一行都不执行continue语句，那么就已经成功地改变了一对alpha，同时可以增加alphaPairsChanged的值。在for循环之外，需要检查alpha值是否做了更新，如果有更新则将iter设为0后继续运行程序。只有在所有数据集上遍历maxIter次，且不再发生任何alpha修改之后，程序才会停止并退出while循环。

利用完整Platt SMO算法加速优化
在几百个点组成的小规模数据集上，简化版SMO算法的运行时没有什么问题的，但是在更大的数据集上的运行速度就会变慢。上面讨论了简化版SMO算法，下面我们就讨论完整的Platt SMO算法。在这两个版本中，实现alpha的更改和代数运算的优化相同。在优化过程中，唯一的不同就是选择alpha的方式。完整的Platt SMO算法应用了一些能够提速的启发方法。
Platt SMO算法是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间进行交替：一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界alpha中实现单遍扫描。而所谓的非边界alpha指的是那些不等于边界0或C的alpha值。对整个数据集的扫描相当容易，而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后在对这个表进行遍历。同时，该步骤会跳过那些已知的不会改变的alpha值。

在选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha值。在优化过程中，会通过最大化步长的方式来获取第二个alpha值。在简化版SMO算法中，我们会在选择j之后计算错误率Ej。但在这里，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说Ei-Ej最大的alpha值。

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)

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 #选择具有最大步长的j
        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

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]

首要的事情就是建立一个数据结构来保存所有的重要值，而这个过程可以通过一个对象来完成。通过数据移到结构中的实现，可以省掉手工输入的麻烦。
eCache的第一列给出的是eCache是否有效的标志位，而第二列给出的是实际的E值。

对于给定的alpha，第一个辅助函数calcEk()能够计算E值并返回。以前，该过程是采用内嵌的方式来完成的，但是由于该过程在这个版本的SMO算法中出现频繁，这里必须要将其单独拎出来。

下一个函数selectJ()用于选择第二个alpha或者说内循环alpha值。回想一下，这里的目标是选择合适的第二个alpha值以保证每次优化中采用最大步长。该函数的误差值与第一个alpha值Ei和下标i有关。首先将输入值Ei在缓存中设置成为有效的。这里的有效意味着它已经计算好了。在eCache中，代码nonzero(oS.eCache[:,0.A])[0]构建出了一个非零表。程序会在所有的值上进行循环并选择其中使得改变最大的那个值。如果是第一次循环，那么久随机选择一个alpha值。

updateEk()用于计算误差值并存入缓存当中。在对alpha值进行优化之后会用这个值。

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

这个函数和smoSimple()中大体是类似的，只是采用了自己的数据结构。
而用SelectJ()替换selectJrand()来选择第二个alpha值。
最后，在alpha值改变时更新Ecache。

完整版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

这里面maxIter变量和函数smoSimple()中的作用有点不同，后者当没有任何alpha发生改变时将整个几何的一次遍历过程记成一次迭代，而这里的一次迭代定义为一次循环过程，而不管该循环具体做了什么事。此时，如果在优化过程中存在波动就会停止，因此这里的做法由于smoSimple()
while循环的内部和smoSimple()中有所不同，一开始for循环在数据集上遍历任何可能的alpha。我们通过调用innerL()来选择第二个alpha，并在可能时对其进行优化处理。如果有任意一对alpha值发生改变，那么会返回1。第二个for循环遍历所有的非边界alpha值，也就是不在边界0或C上的值。

接下来，我们对for循环在非边界循环和完整遍历之间进行切换，并打印出迭代次数。最后程序将会返回常数b和alpha值。

刚刚花了大量时间来计算那些alpha值，但是如何利用它们进行分类呢？
首先，必须基于alpha值得到超平面，它包括了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

有时候我们需要将特征空间转换到另一种特征空间，常常采用的是核函数。
下面为常用的径向基核函数。

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