支持向量机

SVM

分隔超平面（separating hyperplane）：将数据集分隔开来的直线，简称分割面。
间隔：点到分隔面的距离。
支持向量（support vector）：离分隔超平面最近的那些点。
显然，构建分类器的时候，数据点离决策边界越远，那么其最后的预测结果也就越可信。即让间隔尽可能大，提高分类器健壮性。

分类器求解的优化问题

使用类似海维赛德阶跃函数（即单位阶跃函数）的函数对wTx+b作用得到f(wTx+b)，其中当u<0时f(u)输出1，反之则输出+1。
即：

如果数据点处于离分隔超平面很远的位置时，不管是+1类还是-1类，label * (wTx+b)仍然是一个很大的正数。

现在的目标就是找出分类器定义中的w和b。为此，我们必须找到具有最小间隔的数据点，而这些数据点也就是前面提到的支持向量。一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化。
这就可以写作：
                                    
进行优化（固定其中一个因子而最大化其他因子，得到约束条件，拉格朗日乘子法）：
                                      

SVM 应用的一般框架

**

SMO 高效优化算法

**

简化版SMO（运行速度慢）

Eta是alpha[j]的最优修改量，该过程对真实SMO算法进行了简化处理。如果eta为0，那么计算新的alpha[j]就比较麻烦了，因此忽略这一部分。

from numpy import*
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

#随机选择另一alpha
def selectJrand(i,m):
    j=i
    while(j==i):
        j=int(random.uniform(0,m)) #random.uniform(),随机生成（包括小数）
    return j

#调整大于H小于L的alpha
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):#数据集、类别标签、C、容错率、迭代次数
    dataMatrix=mat(dataMatIn)
    labelMat=mat(classLabels).transpose()  #标签列向量
    b=0
    m,n=shape(dataMatrix)
    alphas=mat(zeros((m,1)))  #alpha列矩阵
    iter=0          #迭代计数
    while(iter<maxIter): 
        alphaPairsChanged=0  #al... 记录alpha是否优化，初始为0
        for i in range(m):
            fXi=float(multiply(alphas,labelMat).T* \
                      (dataMatrix*dataMatrix[i,:].T))+b  #alpha*labelMat,对应位置alpha与分类乘积，dataMatrix*dataMatrix[i,:].T，每一行与第i行乘积之和组成新的行
            Ei=fXi-float(labelMat[i])   #误差
            if((labelMat[i]*Ei<-toler) and (alphas[i]<C)) or \
               ((labelMat[i]*Ei>toler)and \
                (alphas[i]>0)):
                j=selectJrand(i,m)               #随机另一个alpha
                fXj=float(multiply(alphas,labelMat).T* \
                      (dataMatrix*dataMatrix[j,:].T))+b
                Ej=fXj-float(labelMat[j])        
                alphaIold=alphas[i].copy()
                alphaJold=alphas[j].copy()
                #计算L、H,用于将alpha[j]调整到0到C之间
                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  #continue,结束本次循环，直接下一次for循环
                eta=2.0*dataMatrix[i,:]*dataMatrix[j,:].T- \
                     dataMatrix[i,:]*dataMatrix[i,:].T- \
                     dataMatrix[j,:]*dataMatrix[j,:].T   #eta是alpha[j]的最优修改量
                if eta>=0: print("eta>=0");continue  
                alphas[j]-=labelMat[j]*(Ei-Ej)/eta    #计算新alpha[j]
                alphas[j]=clipAlpha(alphas[j],H,L)    #调整alpha[j]
                if(abs(alphas[j]-alphaJold)<0.00001):
                    print("j not moving enough");continue
                alphas[i]+=labelMat[j]*labelMat[i]* \
                            (alphaJold-alphas[j])       #修正alpha[i],对他进行同样改变,修正大小相等，方向相反
                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

利用完整 Platt SMO 算法加速优化

class optStruct:
#误差缓存？？？？？？？？？？
    def __init__(self,dataMatIn,classLabels,C,toler):
        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)))   #eCache的第一列给出的是eCache是否有效的标志位，而第二列给出的是实际的E值。

#对于给定的alpha,返回E值
def clacEk(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

#选择第二个alpha,与Ei下标i有关
def selectJ(i,oS,Ei):
    maxK=-1
    maxDeltaE=0
    Ej=0
    oS.eCache[i]=[1,Ei]
    validEcacheList=nonzero(oS.eCache[:,0].A)[0] #返回非零E值所对应的alpha值
    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=daltaE;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=clacEk(oS,i)
    if((oS.labelMat[i]*Ei<-oS.tol)and(oS.alphas[i]<oS.C)) or \
        ((oS.labelMat[i]*Ei>oS.til)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 movong 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.S[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(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
    iter=0
    entireSet=True
    alphaPairsChanged=0
    while(iter<maxIter)and((alphaPairsChanged>0)or(entireSet)):
        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

由于进入了考试周，比较忙，后续内容会慢慢补上。。。。