SVM支持向量机代码实现流程

本文不讲述SVM算法和SMO的推导过程,因为网上已经有很多相关资料,并且我对推导过程还未理解透彻,所以避免误导他人.

本文主要是我通过参考他人的博客和他人实现的代码,总结并实现了SVM算法.其中用到SMO算法来拟合参数.如若有差错误的,期望麻烦给我留言指正.

SMO算法流程:

外循环(优先遍历间隔边界上的支持向量点,即0<αi<C的点,也叫做非边界数据样本):
(1) 查找第一个αi:
      顺序计算一个样本xi
      计算个g(xi)=∑[k从1到N]αk*yk*K(xi,xk)+b,看是否违反KKT条件
      当 αi = 0       <=>  yi*g(xi)>1
           0<αi<C   <=>  yi*g(xi)=1
           αi = C      <=>  yi*g(xi)<1
      选择违反KKT条件最严重的样本点,修改其αi.
      事实上:
      Ei=g(xi)-yi
      yi*Ei=yi*g(xi)-yi^2=yi*g(xi)-1
      所以:
      当KKT条件可改为以下:
      当 αi = 0        <=>   yi*Ei>0
           0<αi<C    <=>    yi*Ei=0
           αi = C       <=>    yi*Ei<0    
内循环:
(2) 查找第二个αj
      计算所有样本的E,E=g(x)-y
      选择最大|Ei-Ej|的样本
(3) 更新优化αi和αj
      计算约束值,边界值:
      如果yi != yj  则:L=max(0,αj-αi),H=min(C,C+αj-αi)
      如果yi = yj   则:L=max(0,αi+αj-C),H=min(C,αi+αj)
      αj=αj-yi(Ei-Ej)/η   η=2K<xi,xj>-K<xi,xi>-K<xj,xj>
      αj = H   if αj>H
          = αj   if L<=αj<=H
          = L    if αj<L
      αi=αi+yiyj(αj'old-αj)
(4) 计算b值
      b1=b-Ei-yi*(αi-αi'old)K<xi,xi>-yj*(αj-αj'old)K<xi,xj>
      b2=b=Ej-yi*(αi-αi'old)K<xi,xj>-yj*(αj-αj'old)K<xj,xj>
     如果0<αi<C和0<αi<C,那么b1=b2
     如果αi=0或者αi=C,同时αj=0或者αj=C,则取中间值b=(b1+b2)/2
     总之:
      b = b1          if 0<αi<C
         = b2          if 0<αj<C
         = (b1+b2)/2   otherwise

循环直到所有样本都满足KKT条件,则结束循环,但一般设置个容忍值,所有样本在容忍值范围内则可以结束

SVM代码(Python实现):

# -*- coding: utf-8 -*-from numpy import *import timeimport matplotlib.pyplot as plt import operator#核函数(点积函数)def calcKernelValue(x1,x2):    value=0    for i in range(len(x1)):        value=value+x1[i]*x2[i]    return float(value)#计算核函数矩阵,存放K<x1,x2>的值def calcKernelMatrix(dataSetIn):    kernelMatrix=zeros((len(dataSetIn),len(dataSetIn)))    for i in range(len(dataSetIn)):        for j in range(len(dataSetIn)):            kernelMatrix[i][j]=calcKernelValue(dataSetIn[i],dataSetIn[j])    return kernelMatrixclass SVM:    def __init__(self,dataSetIn,labelsIn,maxiter,C,error_admit):        self.X=dataSetIn        self.Y=labelsIn        self.m=len(labelsIn)#样本数        self.alphas = [0]*self.m#阿尔法参数数组        self.b=0        self.maxiter=maxiter#最大迭代次数        self.C=C#惩罚因子        self.error_admit=error_admit#容错值        self.kernelValue = calcKernelMatrix(dataSetIn)#核函数矩阵        self.fx = [0]*self.m        self.E = [0]*self.m#计算g(xi)=∑[k从1到N]αk*yk*K(xi,xk)+b    def calcfx(svm):    for i in range(svm.m):        svm.fx[i]=sum(map(operator.mul,map(operator.mul,svm.alphas,svm.Y),svm.kernelValue[i]))+svm.b#计算E(xi)def calcE(svm):    for k in range(svm.m):        svm.E[k]=svm.fx[k]-svm.Y[k]#通过比较|Ei-Ek|,获取第二个阿尔法参数下标jdef findSecondAlphas(svm,i):    maxEdiff=-99999    j=-1    for k in range(svm.m):        if i!=k and abs(svm.E[i]-svm.E[k])>maxEdiff:            maxEdiff=abs(svm.E[i]-svm.E[k])            j=k    j=int(random.uniform(0, svm.m))

    #这里如果用比较|Ei-Ek|的方法获取第二个拉格朗日乘子的话,会迭代进入死循环,所以我改为用随机获取第二个乘子,则结果好很多.最终的支持向量也分布在分割面附近    while i==j:        j=int(random.uniform(0, svm.m))    return j#构建SVM分类器def structure_svm(dataSet,labels,maxiter=1000,C=0.6,error_admit=0):    svm=SVM(dataSet,labels,maxiter,C,error_admit)    continue_iter=True    iter_num=1    while iter_num<svm.maxiter and continue_iter:        continue_iter=False        iter_num+=1        #更新fx(x)和E(x)        calcfx(svm)        calcE(svm)        #KKT条件:        #当 αi=0 <=>  yi*Ei>0        #0<αi<C  <=>  yi*Ei=0        #αi=C  <=>  yi*Ei<0         #所以          #如果yi*Ei < 0并且αi < C 则违反了KKT条件         #因为原本yi*Ei < 0 应该对应的是αi = C         #同理，yi*Ei > 0并且αi > 0则违反了KKT条件         #因为原本yi*Ei  > 0对应的应该是αi =0         i=-1        for k in range(svm.m):            yi_Ei=svm.Y[k]*svm.E[k]            #先遍历间隔边界上的支持向量点,即0<αi<C的点,也叫做非边界数据样本            if svm.alphas[k]>0 and svm.alphas[k]<C and yi_Ei!=0 :                i=k                continue_iter=True                break        #在遍历整个样本集,查找违反KKT条件的样本        if not continue_iter:            for k in range(svm.m):                yi_Ei=svm.Y[k]*svm.E[k]                if (yi_Ei<-error_admit and svm.alphas[k]<C) or (yi_Ei>error_admit and svm.alphas[k]>0) or (yi_Ei==0 and (svm.alphas[k]==0 or svm.alphas[k]==C)):                    i=k                    continue_iter=True                    break                        if continue_iter and i>=0:            j=findSecondAlphas(svm,i)            #确定边界            #如果yi!=yj  则:L=max(0,αj-αi),H=min(C,C+αj-αi)            #如果yi=yj   则:L=max(0,αi+αj-C),H=min(C,αi+αj)            #αj=αj-yi(Ei-Ej)/η   η=2K<xi,xj>-K<xi,xi>-K<xj,xj>            #αj= H   if αj>H            #  = αj  if L<=αj<=H            #  = L   if αj<L            #αi=αi+yiyj(αj'old-αj)            L,H=0,0            if svm.Y[i]!=svm.Y[j]:                L=max(0,svm.alphas[j]-svm.alphas[i])                H=min(C,C+svm.alphas[j]-svm.alphas[i])            else:                L=max(0,svm.alphas[j]+svm.alphas[i]-C)                H=min(C,svm.alphas[j]+svm.alphas[i])                            eta=2*svm.kernelValue[i,j]-svm.kernelValue[i,i]-svm.kernelValue[j,j]            #更新alpha_j            #print(svm.E[i],svm.E[j],svm.E[i]-svm.E[j],eta)            alphas_j_new=svm.alphas[j]-svm.Y[j]*(svm.E[i]-svm.E[j])/eta            if alphas_j_new>H:                alphas_j_new=H            elif alphas_j_new<L:                alphas_j_new=L            else:                alphas_j_new=alphas_j_new            #更新alpha_i            #print(svm.alphas[j],alphas_j_new,svm.alphas[j]-alphas_j_new)            alphas_i_new=svm.alphas[i]+svm.Y[i]*svm.Y[j]*(svm.alphas[j]-alphas_j_new)            #更新b值            #b1=b-Ei-yi*(αi-αi'old)K<xi,xi>-yj*(αj-αj'old)K<xi,xj>            #b2=b=Ej-yi*(αi-αi'old)K<xi,xj>-yj*(αj-αj'old)K<xj,xj>            #如果0<αi<C和0<αi<C,那么b1=b2            #如果αi=0或者αi=C,同时αj=0或者αj=C,则取中间值b=(b1+b2)/2            #总之:            #b= b1          if 0<αi<C            # = b2          if 0<αj<C            # = (b1+b2)/2   otherwise            b1=svm.b-svm.E[i]-svm.Y[i]*(alphas_i_new-svm.alphas[i])*svm.kernelValue[i,i]-svm.Y[j]*(alphas_j_new-svm.alphas[j])*svm.kernelValue[i,j]            b2=svm.b-svm.E[j]-svm.Y[i]*(alphas_i_new-svm.alphas[i])*svm.kernelValue[i,j]-svm.Y[j]*(alphas_j_new-svm.alphas[j])*svm.kernelValue[j,j]            if alphas_i_new>0 and alphas_i_new<C:                b=b1            elif alphas_j_new>0 and alphas_j_new<C:                b=b2            else:                b=(b1+b2)/2            alphas_i_old=svm.alphas[i]            alphas_j_old=svm.alphas[j]            b_old=svm.b            svm.alphas[i]=alphas_i_new            svm.alphas[j]=alphas_j_new            svm.b=b    print svm.alphas    return svm #画图def showSVM(svm):    #画出所有样本点    for i in xrange(svm.m):        if svm.Y[i] == -1:            plt.plot(svm.X[i][0], svm.X[i][1], 'or' ,c='red')        elif svm.Y[i] == 1:            plt.plot(svm.X[i][0], svm.X[i][1], 'or' ,c='green')    for i in xrange(svm.m):        if svm.alphas[i] != 0:            plt.plot(svm.X[i][0], svm.X[i][1], 'or' ,c='blue')     w1,w2=0,0    min_x,max_x=99999,-99999    #计算w,和最小,最大的x1    for i in xrange(svm.m):        w1 += svm.alphas[i]*svm.Y[i]*svm.X[i][0]        w2 += svm.alphas[i]*svm.Y[i]*svm.X[i][1]        if min_x>svm.X[i][0]:            min_x=svm.X[i][0]        if max_x<svm.X[i][0]:            max_x=svm.X[i][0]    #计算x2            y_min_x = float(-svm.b - w1 * min_x) / w2    y_max_x = float(-svm.b - w1 * max_x) / w2    #画出分界线    plt.plot([min_x, max_x], [y_min_x, y_max_x], '-g')    plt.show()    dataSet = []labels = []fileIn = open('testSet.txt')for line in fileIn.readlines():    lineArr = line.strip().split(' ')    dataSet.append([float(lineArr[0]), float(lineArr[1])])    labels.append(float(lineArr[2]))svm=structure_svm(dataSet,labels)showSVM(svm)

训练样本示例(整个样本已上传附件):

3.542485 1.977398 -1
3.018896 2.556416 -1
7.551510 -1.580030 1
2.114999 -0.004466 -1