<学习笔记> 支持向量机

2018-01-22

书中说理解支持向量机需要掌握一定的理论知识，而且读懂的它的代码难度很大，瞬间我的兴趣就来了。
支持向量机（Support Vector Machines）被有些人称为最好的现成的分类器，不加修改，就可以得到低错误率的结果。
SVM有很多种实现，书中只介绍了最流行的一种-序列最小优化算法，那么我就先把这种算法做一个笔记哈。
1.基础概念
考虑我们分隔数据的方法，将我们的数据点放入一个二维的坐标图中，如果我们的数据可以经由一条直线很容易的分开，
这组数据就可以被称为线性可分，这条直线称为分隔超平面。将我们的这个平面图放入三维的图中，很明显我们的这条直线其实
是一个平面的截面，更高维以此类推，顾名思义。
如果数据集是一个1024维的，那么我们就需要一个1023维的某对象进行分隔。这个某对象就是我们的超平面，也就是分类的决策边界。
我们希望我们基于一下思想来构造分类器：
如果数据点离我们的超平面距离越远，那么预测结果也就越可信。
我们的方案是找到离超平面最近的点，确保它们离分隔面的距离尽可能远。这些点到分隔面的距离我们称为间隔。
支持向量就是离超平面最近的那些点，接着就要最大化支持向量到分隔面的距离，我们就要找这种优化方法。
分隔超平面的形式我们可以很轻易的写出：

1	y=(w^T)x+b

那么根据点到直线的距离公式，可写出任意一点a距离：

1	L=\|(w^T)a+b\|/\|w\|

常数b类似与logistic回归中的截距，向量w和截距b一同描述了超平面。
2.分类器工作原理
在logistic回归中，我们使用的输出标签是0和1，这里我们将会使用-1和+1。换用这两个是因为它们只相差一个符号，方便数学上的处理。
我们可以通过一个同=统一的公式来表示间隔或者数据点到超平面之间的距离，不必担心一些无关紧要的问题。
当计算分隔面的距离并确定分隔面的放置位置的时候，间隔通过

1	label*[(w^T)x+b]

来计算，这时就能看出来我们标签选择-1和+1的好处了。
如果数据点处于正方向或者负方向并且远离分隔面的时候，我们得到的结果都会是很大的正数。
我们的目标就是要找出能够描述边界的向量w和截距b。
接下来，我们就可以将问题转化为寻找具有最小间隔的数据点，也就是我们的支持向量。然后我们就对该间隔进行最大化处理：

1	arg max{min(label(w^Tx+b))(1/\|w\|)}

我们采用固定其中一个因子然后最大化其他因子的方式进行优化。令我们的所有的支持向量label[(w^T)x+b]均为1，然后
通过求我们的|w|^(-1)的最大值来求解。但是并不似每一个数据点的label[(w^T)x+b]都为1，只有支持向量才能是1，
那么我们就给定一些约束条件，

1	label*[(w^T)x+b]>=1

在此约束条件下，我们采用一种著名的求解方法，拉格朗日乘子法。通过引入拉格朗日乘子，就可以进行优化求解了。
那么我们的优化函数就变了另外一种形式：

1	max[∑α-(1/2)label(i)label(j)ai*aj<x(i),x(j)>]

其中最后尖括号的两项表示xi与xj的内积。
其约束条件为：

1	α>=0 and ∑αi*label(i)=0

鉴于我们要的目标，数据必须完全可分。那么我们此时得到的数据用书里的话说就是不那么干净，此时通过引入松弛变量来
放宽松我们的标准，使得我们可以允许一些小的错误的分类。
那么我们的约束条件要进行修改：

1	C>α>=0 and ∑αi*label(i)=0

这里的常熟C用来控制“最大化间隔”和“保证大部分点的函数间隔小于1.0”的权重。
我们的问题接下来就又转化为求参数α，我们的支持向量机的工作就是如此。
3.SVM应用的基本框架
SMO高效优化算法
工作原理：
每次循环选择两个α的值进行优化，一旦找到合适的一对，那么就增大其中一个减小另外一个。
所谓合适，就是它们要符合两个条件：
(1)两个alpha的值要在间隔边界之外；
(2)两个alpha未进行过区间化处理或者不在边界上。
优化最佳alpha对：在数据集上遍历每一个alpha，然后在剩下的alpha集合中随机选择另一个，构建alpha对。
由于改变一个alpha会使约束条件失效，因此我们同时改变两个。
为此，我么将会构造一个辅助函数，用于在某个区间内选择一个随机的整数。
同时还有另一个辅助函数，用于在数值太大时进行调整。
下面是这两个辅助函数以及数据载入函数的python实现。

import random
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
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 L>aj:
        aj=L
    return aj

以下代码是我们展示的一个smo算法的示例：

def smosimple(datamatin,classlabels,c,toler,maxiter):
    datamatrix=mat(datamatin)
    labelmatrix=mat(classlabels).transpose()
    b=0
    m,n=shape(datamatrix)
    alphas=mat(zeros(m,1))
    iter=0
    while(iter<maxiter):
        alpahpairschanged=0
        for i in range(m):
            fxi=float(multiply(alphas,labelmatrix).T*(datamatrix*datamatrix[i,:].T))+b
            ei=fxi-float(labelmatrix[i])
            if((labelmatrix[i]*ei<-toler)and(alphas[i]<c))\
                or((labelmatrix[i]*ei>toler)and(alphas[i]>0)):
                j=selectjrand(i,m)
                fxj=float(multiply(alphas,labelmatrix.T*(datamatrix*datamatrix[j,:]).T))+b
                ej=fxj-float(labelmatrix[j])
                alphaiold=alphas[i].copy()
                alphajold=alphas[j].copy()
                if (labelmatrix[i]!=labelmatrix[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]-=labelmatrix[j]*(ei-ej)/eta
                alphas[j]=clipalpha(alphas[j],H,L)
                if (abs(alphas[j]-alphajold)<0.00001):
                    print("j not move enough")
                    continue
                alphas[i]+=labelmatrix[j]*labelmatrix[i]*(alphajold-alphas[j])
                b1=b-ei-labelmatrix[i]*(alphas[i]-alphaiold)*\
                    datamatrix[i,:]*datamatrix[i,:].T-labelmatrix[j]*(alphas[j]-alphajold)*\
                    datamatrix[i,:]*datamatrix[j,:].T
                b2=b-ei-labelmatrix[i]*(alphas[i]-alphaiold)*\
                    datamatrix[i,:]*datamatrix[j,:].T-labelmatrix[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
                alpahpairschanged+=1
                print("iter:%d i:%d,pairs changed %d"%(iter,i,alpahpairschanged))
        if (alpahpairschanged==0):
            iter+=1
        else:
            iter=0
        print("iteration number is %d"%iter)
    return b,alphas

该函数共有5个输入参数：数据集，类别标签，常数c，容错率，最大循环次数。
上述的函数将我们输入的列表转化成numpy矩阵来进行计算。我们将列表进行了装置，因此相当于将行向量变成了列向量。
每次循环当中，将我们的alphapairschanged先设置为0，然后对整个集合进行遍历，这个变量的作用是记录是否对alpha
进行了优化。fxi是我们要预测的类别，ei是我们的误差，误差大，就要对对应的alpha进行优化。在过程当中，我们要对
alpha的值进行检查，保证不能等于0或者c。在后续的处理当中，我们将不在0到c之间的值归为0或者c，那么如果我们发现
它们等于0或者等于c时，它们就不能够进行减小或者增大，不值得再进行优化了。程序当中的L和H便是用来将我们的值控制
在0和c之间。
eta是alpha[j]的最优修改量，如果eta为0，那么我们需要退出当前的for迭代过程。
我们在过程党总检测alpha[j]是否发生了微小的改变，如果有，那么退出for循环。alpha[i]和alpha[j]的改变的大小相同，
但方向是相反的。对它们优化之后，设置一个常数项b。

4.利用完整的platt smo算法进行加速优化
上述我们学习了简化的smo算法，然而当面对非常大的数据量的时候，简化版就会变得非常慢，此时我们就需要完整版了。
二者在优化过程当中唯一的不同在于对alpha的选择方式。
plattsmo算法是通过一个外循环来选择第一个alpha，并且在选择过程当中在两种方式之间交替：
一种方式在所有的数据集上进行扫描，另一种是在非边界的alpha值中实现单遍扫描。
非边界就是不是0也不是c的值。
在实现非边界alpha值的扫描时，首先要建立这些alpha值的列表，然后再对这个表进行遍历，这个过程当中还要跳过那些
不会变化的alpha值。
选择完第一个alpha值之后，通过一个内循环来选择第二个alpha，简化版中，我们在选择j之后计算错误率ej，在加强版中，
我们通过一个全局缓存来保存误差值，然后依据步长最大的alpha值。
下面这段代码包含了一个用于清理代码的数据结构，和三个对e进行缓存的辅助函数。

class optStruct:    #存储重要数据的结构；
    def __init__(self,datamatin,classlebels,c,toler，kTup):
        self.x=datamatin    
        self.labelmat=classlebels
        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))
        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.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=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 updateek(os,k):
    ek=calcek(os,k)
    os.ecache[k]=[1,ek]

这里我们使用一个对象来存储所有的重要值。这里使用对象的目的并不是要进行面向对象编程，而只是
将其作为一个数据结构来使用。在将值传给函数的时候，我们就可以将所有的数据整合为一个对象结构
当中去实现，从而省去了手动输入的麻烦。
ecache矩阵的第一列给出的是ecache是否有效的标志位，第二列给出的是实际的e值。
对于给定的alpha值，辅助函数calcek()能够计算e的值并返回。我们其实是将之前嵌入到程序当中的一部分
提了出来，因为我们要在多个地方使用这部分，将其函数独立化会更方便。
辅助函数selectj()用于内循环的alpha值的选择。核心思想是保证步长最大。
nonzero()语句返回的是非零e值对应的alpha值，程序将所有的值遍历选中其中最大的那个。
辅助函数updateek()计算误差值并且存入缓存当中。上面这段代码自身的作用并不是特别大，但和我们接下里
介绍的优化过程以及外循环组合起来就是强大的smo算法。
用于寻找决策边界的优化例程

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[i]) and (os.c>os.alphas[j]):
            os.b=b2
        else:
            os.b=(b1+b2)/2.0
        return 1
    else:
        return 0

我们可以看到这段代码和我们上面的简化版smo大同小异，但是这里我们已经将重要的数据
整合到了结构当中，这个结构利用对象os进行传递。
除此之外，还有一个重要的修改就是我们使用上面的selectJ()函数来选择alpha的值。
最后在alpha值改变的时候更新缓存。

def smop(datamatin,classlabels,c,toler,maxiter,ktup=('lin',0)):
    os=optStruct(mat(datamatin),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))
            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算法的基本框架，也就是选择第一个alpha值的外循环。
函数一开始构建一个数据结构os来将重要的数据整合起来，退出主循环的条件有两个，一个是当迭代次数超过指定的
最大值，另一个是遍历整个集合未对alpha对进行修改。
这里的maxiter的作用和简化版中的稍有不同，后者当没有alpha值发生变化时会将整个集合的一次遍历过程记作一次迭代，
而在这里的一次迭代则为一次循环过程，而不管它到底做了什么事。
如果再优化过程当中发生波动就会停止，这里的做法优于简化版。
c值给出的是不同优化问题的权重，常数c一方面要保障所有的样例之间的间隔不小于1.0，另一方面又要使得分类间隔尽可能
的大，并且要在这两个方面平衡。
得到了alpha的值以后，我们就要进行真正的分类了，首先必须基于alpha得到超平面，
下面的函数可以实现我们所需要的，

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

上面的函数其实就相当于多个数的乘积，其中alpha表中的值其实大部分都是0，而非零的对应的alpha就是支持向量，
可想而知我们可以通过上面的方式舍弃其中对我们分类没什么作用的0值，留下有用的支持向量。
有时候我们的数据集并不是线性可分的，那么我们就要借助核函数了。
核函数的计算方法：

def kerneltrans(X,A,kTup):
    m,n=shape(X)
    K=mat(zeros(m,1))
    if kTup[0]=='lin':
        K=X*A.T
    elif kTup[0]=='rbf':
        for j in range(m):
            deltaRow=X[j,:]
            K[j]=deltaRow*deltaRow.T
        K=exp(K/(-1*kTup[1]**2))
    else:
        raise NameError('Houston we have a problem\
                        that kernel is not recognized')
    return K

并且在optstruct结构中加入K；
那么当我们需要核函数的时候，调用就好了。
当使用核函数的时候，需要对innerL()和calcEk()两个函数进行修改：

def innerL():
    ...
    eta=2.0*os.K[i,j]-os.K[i,i]-os.K[j,j]
    ...
    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]
    ...
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