统计学习方法第七章的序列最小最优化算法SMO代码实践

统计学习方法第七章的序列最小最优化算法SMO代码实践

#-*- encoding:utf-8 -*-
from numpy import *


#又调试了半天   写代码需谨慎，另外在内循环的函数中标记了代码所在书中公式的位置，可以参考一下
class SVM(object):
    def __init__(self):
        self.dataMat=[]
        self.label=[]
        self.C=0   #惩罚因子
        self.tol=0  #容错率
        self.b=0
        self.kValue={}  #这个是看函数用那种核函数
        self.maxIter=100
        self.suppVectorIndex=[]   #支持向量的下标
        self.suppVector=[]        #支持向量

    #导入数据到dataMat和label中
    def loadDataSet(self,filename):
        fr=open(filename)
        for line in fr.readlines():
            line=line.strip().split('\t')
            self.dataMat.append([float(line[0]),float(line[1])])
            self.label.append(float(line[-1]))
        self.initparam()    #根据数据集进行初始化部分参数

    def initparam(self):
        self.dataMat=mat(self.dataMat)
        self.labelMat=mat(self.label).T
        m,n=shape(self.dataMat)
        self.m=m
        self.alpha=mat(zeros((self.m,1)))    #初始化拉格朗日乘子向量
        self.eCache=mat(zeros((self.m,2)))   #对误差进行缓存
        self.K=mat(zeros((self.m,self.m)))   #存储核函数计算的向量
        for i in xrange(self.m):
            self.K[:,i]=self.kernels(self.dataMat,self.dataMat[i,:])   #计算核函数中的向量

     #计算核函数向量中的第i列所有值
    def kernels(self,dataMat,A):
        m,n=shape(dataMat)
        K=mat(zeros((m,1)))
        if self.kValue.keys()[0]=='linear':   #如果为线性核，则用求线性核的方法计算向量（这个算内积吗？）
            K=dataMat*A.T   #这个是一个m*n维举证与m*1维向量想乘  得到m个向量的内积组成的向量
        elif self.kValue.keys()[0] == 'Gaussian':
            for j in xrange(m):
                delta=dataMat[j,:]-A
                K[j]=delta*delta.T       #可查看书上P122页对高斯核公式的定义
            K=exp(K/(-self.kValue['Gaussian']**2))
        else:raise NameError('can not identify')
        return K

    def chooseJ(self,i,Ei):
        maxK=-1;maxDelta=0;Ej=0
        self.eCache[i]=[1,Ei]
        validEcacheList=nonzero(self.eCache[:,0].A)[0]
        if len(validEcacheList) >1:
            for k in validEcacheList:   #此循环的作用的是遍历所有误差矩阵每一列第一个数 不为0的列找到与i列类别误差相差最大的哪一行
                if k==i:continue
                Ek=self.calcEk(k)
                deltaE=abs(Ei-Ek)
                if (deltaE>maxDelta):
                    maxK=k;maxDelta=deltaE;Ej=Ek
            return maxK,Ej
        else:       #随机选择一行
            j=self.randJ(i)
            Ej=self.calcEk(j)
        return j,Ej

    def randJ(self,i):
        j=i
        while(j==i):
            j=int(random.rand()*self.m)
        return j

    def calcEk(self,k):
        return float(multiply(self.alpha,self.labelMat).T*self.K[:,k]+self.b)-float(self.labelMat[k])

    def adjustAlpha(self,aj,H,L):  #调节拉格朗日乘子的大小
        if aj>H:aj=H
        if L>aj:aj=L
        return aj

    def train(self):
        step=0
        flag=True
        alphaPairsChanged=0
        while(step<self.maxIter) and ((alphaPairsChanged>0) or (flag)):
            alphaPairsChanged=0
            if flag:
                for i in xrange(self.m):
                    alphaPairsChanged += self.innerLoop(i)
                step += 1
            else:     #提取出非支持向量的所有点
                nonBoundIs=nonzero((self.alpha.A>0) * (self.alpha.A<self.C))[0]
                for i in nonBoundIs:     #对非支持向量的点进行更新
                    alphaPairsChanged += self.innerLoop(i)
                step += 1
            if flag:flag=False    #返回支持向量的索引和向量值
            elif(alphaPairsChanged == 0):flag=True
        self.suppVectorIndex=nonzero(self.alpha.A>0)[0]
        self.suppVector=self.dataMat[self.suppVectorIndex]
        self.suppVectorLabel=self.labelMat[self.suppVectorIndex]

    def innerLoop(self,i):
        Ei=self.calcEk(i)
        #判别公式
        if ((self.labelMat[i]*Ei < -self.tol) and (self.alpha[i]  < self.C)) or ((self.labelMat[i]*Ei > self.tol) and (self.alpha[i]>0)):
            j,Ej=self.chooseJ(i,Ei)
            alphaIold=self.alpha[i].copy();alphaJold=self.alpha[j].copy()
            #下面的判别对应于书上P126的对a2new,unc进行调整的公式
            if(self.labelMat[i] != self.labelMat[j]):
                L=max(0,self.alpha[j]-self.alpha[i]);H=min(self.C,self.C+self.alpha[j]-self.alpha[i])
            else:
                L=max(0,self.alpha[j]+self.alpha[i]-self.C);H=min(self.C,self.alpha[j]+self.alpha[i])
            if H==L:return 0
            W=self.K[i,i]+self.K[j,j]-2.0*self.K[i,j]   #根据书上的公式7.107
            if W<0:return 0
            self.alpha[j]=self.alpha[j]+self.labelMat[j]*(Ei-Ej)/W   #可看书上公式7.106,这里是alpha[j]的未经调节前的值
            self.alpha[j]=self.adjustAlpha(self.alpha[j],H,L)
            self.eCache[j]=[1,self.calcEk(j)] #更新j行对应的误差的缓存
            if (abs(self.alpha[j]-alphaJold)<0.00001):return 0
            self.alpha[i] =self.alpha[i]+self.labelMat[i]*self.labelMat[j]*(alphaJold-self.alpha[j])  #对应于书上公式7.109
            self.eCache[i]=[1,self.calcEk(i)]
            #下面的代码对应于书上的公式7.115,7.116
            bi=-Ei-self.labelMat[i]*self.K[i,i]*(self.alpha[i]-alphaIold)-self.labelMat[j]*self.K[i,j]*(self.alpha[j]-alphaJold)+self.b
            bj=-Ej-self.labelMat[i]*self.K[i,j]*(self.alpha[i]-alphaIold)-self.labelMat[j]*self.K[j,j]*(self.alpha[j]-alphaJold)+self.b
            #对应于书上的判断b的值的说明
            if (self.alpha[i]>0 and self.alpha[i]<self.C):self.b=bi
            elif (self.alpha[j] > 0 and self.alpha[j] < self.C): self.b = bj
            else:self.b=(bi+bj)/2.0
            return 1
        else :return 0

made by zcl at CUMT

I know I can because I have a heart that beats