SVM详解(二)

一、问题引入

接上一篇关于SVM的文章SVM详解(一),上一节我们讨论到加入了松弛变量的SVM问题最终演化为求下列带约束条件的问题:

现在我们来看看最初的约束条件是什么样子的:

现在有多少个约束条件就会有多少个αi,根据KKT条件,有:

由于αi≥0 ,而后面那个小于等于0,那么很明显,要使他们的乘积为0,则他们两个中间至少有一个为0。
假设我们现在正在训练一个SVM分类器,在某次迭代后,它对给定样本的分类结果如下图:

其中中间的蓝色的线为当前的分界线,其左右为上下边界,紫色的为正确的分界线。我们通过观察很容易发现图中存在着三类点:

  1. 在上下边界外的点,即yi(Wxi+b)>1,它们是分对了的点。
  2. 在上下边界上的点,即yi(Wxi+b)=1,它们是支持向量。
  3. 在上下边界间的点,即yi(Wxi+b)<1,它们应当要满足松弛变量的条件。

我们应当调整以下两类点对应的αi以使得整个推导过程仍满足KKT条件:

  • yi(Wxi+b)>1但是αi>0,这个很容易理解,因为根据前面的描述,这两个量中至少有一个为0,既然yi(Wxi+b)>1,那么αi肯定需要等于0。

  • yi(Wxi+b)<1但是αi<C,这个就比较难理解了。建议先回去翻看一下关于SVM的第一篇文章的最后一页纸的推导,然后我给出下面的解释:

那么α应当怎么调整呢?比如说某一次迭代完后发现有10个点不满足,也就是有10个α需要调整,那么10个α在一起,你怎么知道去增加或者减少哪一个或者哪几个α呢?这就是下面要介绍的SMO算法需要解决的问题。

二、SMO算法思路

SMO(Sequential Minimal Optimization)即序列最小优化算法的工作原理是:每次循环中选择两个α进行优化处理。假设选取的两个α分别为α1,α2。看看前面有一个条件α1y1+α2y2+...+αn*yn=0,所以在调整前后满足:

至于ϵ是多少,管它多少,又不用它。也就是说我们只要知道了α1,α2中一个调整以后的值,根据条件另一个值不用算就知道。那么怎么求呢?假设我们来求α2吧。观察上面那个等式,y1,y2是标签,要么1要么-1。而两个α>=0。所以新的α是有范围的。假设现在y1=y2=1或−1,先=1吧,那么αnew1+αnew2=αold1+αold2=ϵ 因为αnew1是在0-C之间,所以假设αnew1=0,那么αnew2可以取到最大值为ϵ,也就是αold1+αold2。而αnew2又不能大于C,所以其最大取值为min(C,αold1+αold2)。同理如果αnew1=C,那么αnew2可以取到最小值为ϵ−C,也就是αold1+αold2−C,而αnew2最小不能小于0,那玩意这个值比0小怎么办?所以αnew2的下限值就为max(0,αold1+αold2−C)。这是相等取1,那么相等取-1呢?同理吧,不过是 −αnew1−αnew2=−αold1−αold2=ϵ两边乘以-1,就是αnew1+αnew2=αold1+αold2=−ϵ,因为你也不知道−ϵ是多少,不过一个字母而已,那么用ϵ代替−ϵ吧,接下来的讨论过程一模一样。从而属于同类别的时候αnew2上下限就有了。同理去计算下不同类(1,-1)的时候的上下限。最终也能得到一个类似的结果,这里就省了,给出最后的结果:

其中L为αnew2的下限,H为αnew2的上限。到这只是给出了αnew2的范围,那么它怎么解呢?这个解的过程实在是太复杂了,估计要写好多页纸,建议参考大佬的证明,这里我们直接给出更新公式:

上程序之前先规整一下步骤吧:
(1)根据αold1,αold2和上面的公式找到αnew2的上下限。
(2)计算一个二阶导数,也就是上面的
η=2K(x1,x2)−K(x1,x1)−K(x2,x2),K表示里面两个向量的内积。
(3)根据公式更新αnew2,矫正它的范围,再更新αnew1
(4)按照公式计算两个b,再更新b。

三、SMO算法代码

以下代码来源于Peter Harrington所著的《机器学习实战》一书,为其中关于SMO最简单的一种实现。
两个辅助函数:

#为了使markdown代码高亮显示
import os

#随机选择两个不同的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

#用于调整alpha的值使其处于正确的取值范围内
def clipAlpha(aj,H,L):
    if aj > H: 
        aj = H
    if L > aj:
        aj = L
    return aj

首先给出SMO函数的伪代码:

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

Python实现:

from numpy import *

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]):
                    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

该函数有五个参数,分别是:数据集、类别标签、常数C、容错率和退出前最大的循环次数。该函数基本上就是对前面数学推导结果的复现,看懂应该不是很困难。