本文共 12331 字,大约阅读时间需要 41 分钟。
接下来的py代码里面实现了简易的纯暴力smo算法和基于启发式的smo算法,并测试了rbf核函数,并对手写识别问题的分类进行了再次的考证。
# -*- coding: utf-8 -*-"""照葫芦画瓢完成于2017.5.11 19:55算法名称 : 支持向量机作者: zzt941006"""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,labelMatdef selectJrand(i,m): j = i while(j == i): j = int(random.uniform(0,m)) return jdef clipAlpha(aj,H,L): if aj > H: aj = H if L > aj: aj = L return ajdef smoSimple(dataMatIn, classLabels, C, toler, maxIter): dataMatrix = mat(dataMatIn) # 100 * 2 labelMat = mat(classLabels).transpose() # 1 * 100 b = 0 m,n = shape(dataMatrix) alphas = mat(zeros((m,1))) # 100 * 1 iter = 0 while (iter < maxIter): alphaPairsChanged = 0#标志位用于判定这一对alpha是否已经优化 for i in range(m): fXi = float(multiply(alphas,labelMat).T * (dataMatrix * dataMatrix[i,:].T)) + b #dataMatrix[i,:]取得是矩阵中第i行的所有元素 #multiply做的是元素间的点乘,需要两个矩阵维度一模一样,在此alphas,labelMat点乘后得到100 * 1 维度,转置后得到1 * 100 #而dataMatrix * dataMatrixp[i,:] 是整个矩阵乘以该矩阵的某一行的转置,所以最终维度是 100 * 1 #点乘矩阵转置维度变为1 * 100 在跟后者进行矩阵相乘,得到一个标量(1 * 1) 维度 #在此通过这一些列的计算获取fXi,用来代表我们预测的类别 Ei = fXi - float(labelMat[i])#与真实标签的误差计算,如果误差很大,就可以基于此对alpha值进行优化 if((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * i > toler) and (alphas[i] > 0)):#调优检验 #不在边界上的点,且误差很大时,可以对值进行进一步的优化 j = selectJrand(i,m) fXj = float(multiply(alphas,labelMat).T * (dataMatrix * dataMatrix[j,:].T)) + b #随机选取第二个alpha并计算 Ej = fXj - float(labelMat[j]) alphaIold = alphas[i].copy() alphaJold = alphas[j].copy()#要用copy方法来赋值,这样可以让其有新的内存,否则无法发现改变 if(labelMat[i] != labelMat[j]):#边界调整,防止alpha[j] > C 或者 < 0 L = max(0, alphas[j] - alphas[i]) H = min(C, C + alphas[j] - alphas[i]) else: L = max(0, alphas[i] + alphas[j] - 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 #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) if(abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j]) #更新alpha[i],修改量与j相同??? 方向相反显然 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 (0 < alphas[i]) and (C > alphas[j]) : b = b1 elif (0 < alphas[j]) and (C > alphas[j]): 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,alphasdef 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('Houston We Have a Problem -- That Kernel is not recognized') return Kclass 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))) 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) # 调用该函数,对矩阵K进行填充def calcEk(oS,k): fXk = float(multiply(oS.alphas,oS.labelMat).T * oS.K[:,k] + oS.b) Ek = fXk - float(oS.labelMat[k]) return Ekdef 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 = deltaE Ej = Ek return maxK,Ej else : j = selectJrand(i,oS.m) Ej = calcEk(oS,j) return j,Ejdef updateEk(oS,k): Ek = calcEk(oS,k) oS.eCache[k] = [1,Ek]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.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 oS.alphas[j] = clipAlpha(oS.alphas[j],H,L) updateEk(oS,j) if(abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enouth"; 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.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 0def 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)): 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 #toggle entire set loop elif (alphaPairsChanged == 0): entireSet = True print "iteration number: %d" % iter return oS.b,oS.alphasdef 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) #大部分alpha为零向量,因此非零向量即为支持向量 return wdef testRbf(k1 = 1.30): dataArr,labelArr = loadDataSet('testSetRBF.txt') b,alphas = smoP(dataArr,labelArr,200,0.0001,10000,('rbf',k1)) datMat = mat(dataArr) labelMat =mat(labelArr).transpose() svInd = nonzero(alphas.A > 0)[0] sVs = datMat[svInd] labelSV = labelMat[svInd] print "there ar %d Support Vectors" % shape(sVs)[0] m,n = shape(datMat) errorCount = 0 for i in range(m): kernelEval = kernelTrans(sVs,datMat[i,:],('rbf',k1)) predict = kernelEval.T * multiply(labelSV,alphas[svInd]) + b if sign(predict) != sign(labelArr[i]): errorCount += 1 print "the training error rate is: %f" % (float(errorCount) / m) dataArr,labelArr = loadDataSet('testSetRBF2.txt') b,alphas = smoP(dataArr,labelArr,200,0.0001,10000,('rbf',k1)) datMat = mat(dataArr) labelMat =mat(labelArr).transpose() svInd = nonzero(alphas.A > 0)[0] sVs = datMat[svInd] labelSV = labelMat[svInd] print "there ar %d Support Vectors" % shape(sVs)[0] m,n = shape(datMat) errorCount = 0 for i in range(m): kernelEval = kernelTrans(sVs,datMat[i,:],('rbf',k1)) predict = kernelEval.T * multiply(labelSV,alphas[svInd]) + b if sign(predict) != sign(labelArr[i]): errorCount += 1 print "the training error rate is: %f" % (float(errorCount) / m)def img2vector(filename): returnVect = zeros((1,1024));#将32*32的二进制图像转化为1*1024维度的矩阵 fr = open(filename) for i in range(32): lineStr = fr.readline() # print lineStr for j in range(32): returnVect[0,32*i+j] = int(lineStr[j])#第一行为0-31 第二行加在后面为32-63....直到第32行的32个数字映射为 992-1023 return returnVectdef loadImapge(dirName): from os import listdir hwLabels = [] trainingFileList = listdir(dirName) m = len(trainingFileList) trainingMat = zeros((m,1024)) for i in range(m): fileNameStr = trainingFileList[i] fileStr = fileNameStr.split('_')[0] classNumStr = int(fileStr.split('_')[0]) if classNumStr == 9:hwLabels.append(-1) else: hwLabels.append(1) trainingMat[i,:] = img2vector('%s/%s' %(dirName,fileNameStr)) return trainingMat,hwLabelsdef 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 "there are %d Support Vectors" % 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 "the training error rate is: %f" % (float(errorCount)/m) dataArr,labelArr = loadImages('testDigits') errorCount = 0 datMat=mat(dataArr); labelMat = mat(labelArr).transpose() m,n = shape(datMat) 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 "the test error rate is: %f" % (float(errorCount)/m)
转载地址:http://zupws.baihongyu.com/