# -*-coding:utf-8 -*- import matplotlib.pyplot as plt import numpy as np import random """ Author: Jack Cui Blog: http://blog.csdn.net/c406495762 Zhihu: https://www.zhihu.com/people/Jack--Cui/ Modify: 2017-10-03 """ class optStruct: """ 数据结构，维护所有需要操作的值 Parameters： dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 kTup - 包含核函数信息的元组,第一个参数存放核函数类别，第二个参数存放必要的核函数需要用到的参数 """ def __init__(self, dataMatIn, classLabels, C, toler, kTup): self.X = dataMatIn #数据矩阵 self.labelMat = classLabels #数据标签 self.C = C #松弛变量 self.tol = toler #容错率 self.m = np.shape(dataMatIn)[0] #数据矩阵行数 self.alphas = np.mat(np.zeros((self.m,1))) #根据矩阵行数初始化alpha参数为0 self.b = 0 #初始化b参数为0 self.eCache = np.mat(np.zeros((self.m,2))) #根据矩阵行数初始化虎误差缓存，第一列为是否有效的标志位，第二列为实际的误差E的值。 self.K = np.mat(np.zeros((self.m,self.m))) #初始化核K for i in range(self.m): #计算所有数据的核K self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup) def kernelTrans(X, A, kTup): """ 通过核函数将数据转换更高维的空间 Parameters： X - 数据矩阵 A - 单个数据的向量 kTup - 包含核函数信息的元组 Returns: K - 计算的核K """ m,n = np.shape(X) K = np.mat(np.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 = np.exp(K/(-1*kTup[1]**2)) #计算高斯核K else: raise NameError('核函数无法识别') return K #返回计算的核K def loadDataSet(fileName): """ 读取数据 Parameters: fileName - 文件名 Returns: dataMat - 数据矩阵 labelMat - 数据标签 """ 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 def calcEk(oS, k): """ 计算误差 Parameters： oS - 数据结构 k - 标号为k的数据 Returns: Ek - 标号为k的数据误差 """ fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b) Ek = fXk - float(oS.labelMat[k]) return Ek def selectJrand(i, m): """ 函数说明:随机选择alpha_j的索引值 Parameters: i - alpha_i的索引值 m - alpha参数个数 Returns: j - alpha_j的索引值 """ j = i #选择一个不等于i的j while (j == i): j = int(random.uniform(0, m)) return j def selectJ(i, oS, Ei): """ 内循环启发方式2 Parameters： i - 标号为i的数据的索引值 oS - 数据结构 Ei - 标号为i的数据误差 Returns: j, maxK - 标号为j或maxK的数据的索引值 Ej - 标号为j的数据误差 """ maxK = -1; maxDeltaE = 0; Ej = 0 #初始化 oS.eCache[i] = [1,Ei] #根据Ei更新误差缓存 validEcacheList = np.nonzero(oS.eCache[:,0].A)[0] #返回误差不为0的数据的索引值 if (len(validEcacheList)) > 1: #有不为0的误差 for k in validEcacheList: #遍历,找到最大的Ek if k == i: continue #不计算i,浪费时间 Ek = calcEk(oS, k) #计算Ek deltaE = abs(Ei - Ek) #计算|Ei-Ek| if (deltaE > maxDeltaE): #找到maxDeltaE maxK = k; maxDeltaE = deltaE; Ej = Ek return maxK, Ej #返回maxK,Ej else: #没有不为0的误差 j = selectJrand(i, oS.m) #随机选择alpha_j的索引值 Ej = calcEk(oS, j) #计算Ej return j, Ej #j,Ej def updateEk(oS, k): """ 计算Ek,并更新误差缓存 Parameters： oS - 数据结构 k - 标号为k的数据的索引值 Returns: 无 """ Ek = calcEk(oS, k) #计算Ek oS.eCache[k] = [1,Ek] #更新误差缓存 def clipAlpha(aj,H,L): """ 修剪alpha_j Parameters: aj - alpha_j的值 H - alpha上限 L - alpha下限 Returns: aj - 修剪后的alpah_j的值 """ if aj > H: aj = H if L > aj: aj = L return aj def innerL(i, oS): """ 优化的SMO算法 Parameters： i - 标号为i的数据的索引值 oS - 数据结构 Returns: 1 - 有任意一对alpha值发生变化 0 - 没有任意一对alpha值发生变化或变化太小 """ #步骤1：计算误差Ei Ei = calcEk(oS, i) #优化alpha,设定一定的容错率。 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)): #使用内循环启发方式2选择alpha_j,并计算Ej j,Ej = selectJ(i, oS, Ei) #保存更新前的aplpha值，使用深拷贝 alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy(); #步骤2：计算上下界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 #步骤3：计算eta eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] if eta >= 0: print("eta>=0") return 0 #步骤4：更新alpha_j oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta #步骤5：修剪alpha_j oS.alphas[j] = clipAlpha(oS.alphas[j],H,L) #更新Ej至误差缓存 updateEk(oS, j) if (abs(oS.alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小") return 0 #步骤6：更新alpha_i oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j]) #更新Ei至误差缓存 updateEk(oS, i) #步骤7：更新b_1和b_2 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] #步骤8：根据b_1和b_2更新b 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)): """ 完整的线性SMO算法 Parameters： dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 maxIter - 最大迭代次数 kTup - 包含核函数信息的元组 Returns: oS.b - SMO算法计算的b oS.alphas - SMO算法计算的alphas """ oS = optStruct(np.mat(dataMatIn), np.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) #使用优化的SMO算法 print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged)) iter += 1 else: #遍历非边界值 nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] #遍历不在边界0和C的alpha for i in nonBoundIs: alphaPairsChanged += innerL(i,oS) print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged)) iter += 1 if entireSet: #遍历一次后改为非边界遍历 entireSet = False elif (alphaPairsChanged == 0): #如果alpha没有更新,计算全样本遍历 entireSet = True print("迭代次数: %d" % iter) return oS.b,oS.alphas #返回SMO算法计算的b和alphas def img2vector(filename): """ 将32x32的二进制图像转换为1x1024向量。 Parameters: filename - 文件名 Returns: returnVect - 返回的二进制图像的1x1024向量 """ returnVect = np.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): """ 加载图片 Parameters: dirName - 文件夹的名字 Returns: trainingMat - 数据矩阵 hwLabels - 数据标签 """ from os import listdir hwLabels = [] trainingFileList = listdir(dirName) m = len(trainingFileList) trainingMat = np.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, hwLabels def testDigits(kTup=('rbf', 10)): """ 测试函数 Parameters: kTup - 包含核函数信息的元组 Returns: 无 """ dataArr,labelArr = loadImages('trainingDigits') b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10, kTup) datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose() svInd = np.nonzero(alphas.A>0)[0] sVs=datMat[svInd] labelSV = labelMat[svInd]; print("支持向量个数:%d" % np.shape(sVs)[0]) m,n = np.shape(datMat) errorCount = 0 for i in range(m): kernelEval = kernelTrans(sVs,datMat[i,:],kTup) predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1 print("训练集错误率: %.2f%%" % (float(errorCount)/m)) dataArr,labelArr = loadImages('testDigits') errorCount = 0 datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose() m,n = np.shape(datMat) for i in range(m): kernelEval = kernelTrans(sVs,datMat[i,:],kTup) predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1 print("测试集错误率: %.2f%%" % (float(errorCount)/m)) if __name__ == '__main__': testDigits()