分类问题 Python 支持向量机 手写识别系统（学习笔记）

手写识别系统可以使用kNN实现且效果不错，但是使用kNN占用的内存太大，而且必须在保持性能不变的同时使用较少的内存。对于SVM而言，只需保留较少的支持向量就可以实现手写识别系统。
SVM的泛化错误率较低，具有良好的学习能力，而且几乎所有的分类问题都可以使用SVM解决。

'''SVM实现手写识别系统'''
from numpy import *
from os import listdir
import math

#计算核
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,:]-A
            K[j]=deltaRow*deltaRow.T
        K=exp(K/(-1*kTup[1]**2))
    else:#如果遇到无法识别的远足，程序抛出异常
        raise NameError("无法识别")
    return K    

#保存所有重要值，构建一个数据结构
class optStruct:
    def __init__(self,dataMatIn,classLabels,C,toler,kTup):
        self.X=dataMatIn
        self.labelMat=classLabels
        self.C=C
        self.tol=toler
        self.m=shape(dataMatIn)[0]
        self.alphas=mat(zeros((self.m,1)))#mat()转化为矩阵
        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)

#选择第二个alpha的下标
def selectJrand(i,m):#i表示alpha的下标，m表示alpha的总数
    j=i
    while(j==i):
        j=int(random.uniform(0,m))
    return 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

#更新误差值
def updateEk(oS,k):
    Ek=calcEk(oS,k)
    oS.eCache[k]=[1,Ek]

#选择第二个alpha
def selectJ(i,oS,Ei):
    maxK=-1
    maxDeltaE=0
    Ej=0
    oS.eCache[i]=[1,Ei]
    
    validEcacheList=nonzero(oS.eCache[:,0].A)[0]#.A转化为数组，nonzero()非零元素的位置
    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

#调整第二个alpha的范围
def clipAlpha(aj,H,L):
    if aj>H:
        aj=H
    if L>aj:
        aj=L
    return aj

#计算参数,如果任意一对alpha有变化，返回1
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)#第二个alpha
        
        alphaIold=oS.alphas[i].copy()
        alphaJold=oS.alphas[j].copy()
        
        #调整L H
        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
        
        #计算alpha
        eta=2.0*oS.K[i,j]-oS.K[i,i]-oS.K[j,j]
        if eta>=0:
            print("eta>=0")
            return 0 
        oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/eta#计算第二个alpha
        oS.alphas[j]=clipAlpha(oS.alphas[j],H,L)#调整alpha
        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])#计算第一个alpha
        updateEk(oS,i)#更新误差值
        
        #计算b
        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]    
        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

#Platt SMO算法 计算b与alpha
def smoP(dataMatIn,classLabels,C,toler,maxIter,kTup=('lin',0)):
    oS=optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler,kTup)
    
    iter=0
    entireSet=True
    alphaPairsChanged=0
    while(iter<maxIter) and ((alphaPairsChanged>0) or (entireSet)):#迭代次数超过指定最大值或遍历整个集合都没有对任意alpha进行修改
        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("迭代次数：%d"%iter)
    return oS.b,oS.alphas

测试SMO算法

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

dataArr,labelArr=loadDataSet("testSetRBF2.txt")
b,alphas=smoP(dataArr,labelArr,0.6,0.001,40)

输出结果：
……

#把二值化图像转化为向量
def img2vector(filename):
    returnVect=zeros((1,1024))
    fr=open(filename)
    for i in range(32):
        lineStr=fr.readline()
        for j in range(32):
            returnVect[0,32*i+j]=int(lineStr[j])
    return returnVect

#导入数据
def loadImages(dirName):#路径
    hwLabels=[]#分类标签
    trainingFileList=listdir(dirName)#返回指定路径下的文件和文件夹列表
    m=len(trainingFileList)
    trainingMat=zeros((m,1024))
    for i in range(m):
        fileNameStr=trainingFileList[i]
        filrStr=fileNameStr.split('.')[0]
        classNumStr=int(filrStr.split('_')[0])
        if classNumStr==9:
            hwLabels.append(-1)
        else:
            hwLabels.append(1)
        trainingMat[i,]=img2vector("%s/%s"%(dirName,fileNameStr))#路径+文件名
    return trainingMat,hwLabels

#实现基于SVM的手写数字识别
def testDigits(kTup=('rbf',10)):
    #训练集
    dataArr,labelArr=loadImages('trainingDigits')
    b,alphas=smoP(dataArr,labelArr,200,0.0001,10000,kTup)
    datMat=mat(dataArr)
    labelMat=mat(labelArr).transpose()
    
    svInd=nonzero(alphas.A>0)[0]#支持向量的下标
    sVs=datMat[svInd]
    labelSV=labelMat[svInd]
    print("有%d支持向量"%shape(sVs)[0])
    
    m,n=shape(datMat)    
    errorCount=0
    for i in range(m):
        kernelEval=kernelTrans(sVs,datMat[i,:],kTup)#支持向量、数据集、核函数类型
        predict=kernelEval.T*multiply(labelSV,alphas[svInd])+b #预测
        if sign(predict)!=sign(labelArr[i]):
            errorCount+=1
    print("训练数据错误率是: %f"%(float(errorCount)/m))
    
    #测试集
    dataArr,labelArr=loadImages('testDigits')
    datMat=mat(dataArr)
    labelMat=mat(labelArr).transpose()
    m,n=shape(datMat)
    errorCount=0
    for i in range(m):
        kernelEval=kernelTrans(sVs,datMat[i,:],kTup)#支持向量、数据集、核函数类型
        predict=kernelEval.T*multiply(labelSV,alphas[svInd])+b #预测
        if sign(predict)!=sign(labelArr[i]):
            errorCount+=1
        #print("分类结果为：",sign(predict))
    print("测试数据错误率是: %f"%(float(errorCount)/m))

testDigits(('rbf',20))

输出结果：
……

……