《机器学习实战》之支持向量机--基于Python3

Posted 2021-04-06 华少的知识宝典

书上使用SMO算法实现的SVM真的是看的头痛，代码简直又臭又长。看了半天后，我放弃了挣扎，还是直接用sklearn来实现吧。不过书上的代码我也照着敲了一遍，也写了一些注释，所以也写一下留作参考吧！（虽然有些地方我也没太搞懂）

支  持  向  量  机

1、SMO的高效优化算法  

1996年， John Platt发布了一个称为SMO的强大算法，用于训练SVM。SMO表示序列最小优化（Sequential Minimal Optimization）。

• SMO目标：SMO算法的目标是求出一系列alpha和b，一旦求出了这些alpha，就很容易计算出权重向量w并得到分隔超平面

• SMO思想：Platt的SMO算法是将大优化问题分解为多个小优化问题来求解的。这些小优化问题往往很容易求解，并且对它们进行顺序求解的结果与将它们作为整体来求解的结果是完全一致的。在结果完全相同的同时， SMO算法的求解时间短很多。

• SMO原理：每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha，那么就增大其中一个同时减小另一个.

• 这里指的合适必须要符合一定的条件

  1. 这两个 alpha 必须要在间隔边界之外
  

  2. 这两个 alpha 还没有进行过区间化处理或者不在边界上。

• 之所以要同时改变2个 alpha；原因是我们有一个约束条件：

  

          如果只是修改一个 alpha，很可能导致约束条件失效。

2、简化版SMO      

这里使用书上提供的 testSet.txt 文件里的数据来实现支持向量机。

2.1 获得数据

跟以往一样的方法导入数据，当然，还有导入模块。

from numpy import * from time import sleep
"""函数说明: 加载数据集
Parameters: fileName - 文件名Returns: dataMat - 数据矩阵，list类型 labelMat - 数据标签，list类型"""def loadDataSet(fileName): dataMat = []; labelMat = [] fr = open(fileName) for line in fr.readlines(): lineArr = line.strip().split('	') dataMat.append([float(lineArr[0]), float(lineArr[1])]) labelMat.append(float(lineArr[2])) return dataMat,labelMat

2.2 辅助函数

这里定义两个辅助函数，用于后面的SMO算法。

"""函数说明: 随机选择alpha
Parameters: i - 第一个alpha的下标 m - 所有alpha的数目Returns: j - 不等于i的j"""def selectJrand(i,m): j=i #we want to select any J not equal to i while (j==i): j = int(random.uniform(0,m)) return j    """函数说明: 调整大于H或小于L的alpha值
Parameters: aj - alpha值 H - alpha上界 L - alpha下界Returns: aj - alpha值"""def clipAlpha(aj,H,L): if aj > H:  aj = H if L > aj: aj = L return aj

2.3 简化版的SMO

"""函数说明: 简化版的SMO算法
Parameters: dataMatIn - 数据矩阵 classLabels - 数据标签 C - 惩罚参数 toler - 松弛变量 maxIter - 最大迭代次数Returns: b - 模型的常量值 alphas - 拉格朗日乘子"""def smoSimple(dataMatIn, classLabels, C, toler, maxIter): # 转换为numpy的matrix存储 dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()  # 初始化参数 b，统计 dataMatirx 的 维度 b = 0; m,n = shape(dataMatrix) alphas = mat(zeros((m,1))) # 初始化 alpha 参数为 0 iter = 0 #初始化迭代次数  # 限制循环迭代次数，也就是在数据集上遍历 maxIter 次，且不再发生任何 alpha 修改，则循环停止 while (iter < maxIter): # 记录alpha是否已经进行优化，每次循环时设为0，然后再对整个集合顺序遍历 alphaPairsChanged = 0 for i in range(m): # 第 i 个样本的预测类别，其实就是公式的实现, 不过求和是以向量相乘的方式实现的 fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b  # 预测结果与真实结果比对，计算误差 Ei Ei = fXi - float(labelMat[i])            #如果误差很大就对该数据对应的alpha进行优化，正、负间隔都会被测试，同时检查alpha值 if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): j = selectJrand(i,m) # 随机选择不等于i的 0-m 的 第二个 alpha值  # 步骤1：计算误差Ej, 和上面的步骤一样的 fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b Ej = fXj - float(labelMat[j]) # 保存更新前的alpha值, 以备后用(深拷贝) alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();  # 步骤2: 计算上界L和下界H # 这里是对SMO最优化问题的子问题的约束条件的分析 # L和H用于将alphas[j]调整到0-C之间。如果L==H，就不做任何改变，直接执行continue语句 # 如果labelMat[i] != labelMat[j] 表示异侧，就相减，否则是同侧，就相加。 if (labelMat[i] != labelMat[j]): L = max(0, alphas[j] - alphas[i]) H = min(C, C + alphas[j] - alphas[i]) else: L = max(0, alphas[j] + alphas[i] - C) H = min(C, alphas[j] + alphas[i]) if L == H:  print("L==H") continue  # 步骤3:计算 eta eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T if eta >= 0:  print("eta>=0") continue  # 步骤4: 更新alpha[j] alphas[j] -= labelMat[j]*(Ei - Ej)/eta # 并使用辅助函数, 对L和H进行调整 alphas[j] = clipAlpha(alphas[j],H,L) # 检查alpha[j]是否只是轻微的改变，如果是的话，就退出for循环 if (abs(alphas[j] - alphaJold) < 0.00001):  print("j not moving enough") continue  # 步骤5: 更新alpha[i] alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) # update i by the same amount as j # the update is in the oppostie direction # 步骤6: 更新b1和b2 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  # 步骤7:根据b1和b2更新b if (0 < alphas[i]) and (C > alphas[i]): b = b1 elif (0 < alphas[j]) and (C > alphas[j]): b = b2 else: b = (b1 + b2)/2.0  # 程序能运行到这里（没有通过continue提前结束循环），说明已经成功改变了一对alpha了 # 优化次数加1 alphaPairsChanged += 1 print("iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))  # 检查 alpha 是否做了更新 if (alphaPairsChanged == 0): iter += 1 # alphaPairsChanged仍然等于0，说明 alpha 没有 更新，则直接 +1 ，继续往后遍历数据集 else: iter = 0 # 说明有更新，则将 iter 设为 0 后继续运行此程序---重新遍历数据集 print("iteration number: %d" % iter) return b,alphas
if __name__ == '__main__': dataMatIn, classLabels = loadDataSet('testSet.txt') b, alphas = smoSimple(dataMatIn, classLabels, 0.6, 0.001, 40)    print(b)  # 输出算出的 b    print(alphas[alphas > 0])  # 输出算出的 大于0 的 alpha

运行结果如下。这里只显示了一小部分，由于里面的迭代是很长的，所以最后结果一般会很长。

3、完整版SMO         

使用简化版SMO在大的数据集上运行速度会很慢，所以对算法进行优化，使用一些能够提速的启动方式。

3.1 定义一个类

首先，定义一个对象来保存所有的重要值。这里的这个类只是作为一个数据结构来使用。

"""类说明: 数据结构, 维护所有需要操作的值
Parameters: dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率"""class optStruct: def __init__(self,dataMatIn, classLabels, C, toler): # Initialize the structure with the parameters  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))) # 误差缓存

3.2 辅助的支持函数

同样，这里也定义三个辅助的支持函数，用于后面的SMO算法。

"""函数说明: 计算误差
Parameters: oS - 数据结构 k - 标号为k的数据Returns: Ek - 标号为k的数据的误差""" def calcEk(oS, k): fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b Ek = fXk - float(oS.labelMat[k]) return Ek """函数说明: 内循环启发方式2Parameters: i - 标号为i的数据的索引值 oS - 数据结构 Ei - 标号为i的数据误差Returns: j, maxK - 标号为j或maxK的数据的索引值 Ej - 标号为j的数据误差
"""def selectJ(i, oS, Ei): #this is the second choice -heurstic, and calcs Ej # 初始化 maxK = -1; maxDeltaE = 0; Ej = 0 # 根据Ei更新误差缓存 oS.eCache[i] = [1,Ei] #set valid #choose the alpha that gives the maximum delta E # 返回误差不为0的数据的索引值 validEcacheList = nonzero(oS.eCache[:,0].A)[0] if (len(validEcacheList)) > 1: for k in validEcacheList: #loop through valid Ecache values and find the one that maximizes delta E if k == i: continue #don't calc for i, waste of time Ek = calcEk(oS, k) deltaE = abs(Ei - Ek) if (deltaE > maxDeltaE): maxK = k; maxDeltaE = deltaE; Ej = Ek return maxK, Ej else: #in this case (first time around) we don't have any valid eCache values j = selectJrand(i, oS.m) Ej = calcEk(oS, j) return j, Ej

"""函数说明: 计算Ek, 并更新误差缓存
Parameters: oS - 数据结构 k - 标号为k的数据的索引值Returns: 无"""def updateEk(oS, k):#after any alpha has changed update the new value in the cache Ek = calcEk(oS, k) oS.eCache[k] = [1,Ek]

3.3 优化的SMO算法

"""函数说明: 优化的SMO算法
Parameters: i - 标号为i的数据的索引值 oS - 数据结构Returns: 1 - 有任意一对alpha值发生变化 0 - 没有任意一对alpha值发生变化或变化太小"""def innerL(i, oS): # 步骤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) #this has been changed from selectJrand # 保存更新前的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.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].T 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) #added this for the Ecache if (abs(oS.alphas[j] - alphaJold) < 0.00001):  print("j not moving enough") return 0  # 步骤6: 更新alpha_i oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j # 更新Ei至误差缓存 updateEk(oS, i) #added this for the Ecache #the update is in the oppostie direction  # 步骤7: 更新b_1和b_2 b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T  # 步骤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

  3.4 完整的线性SMO算法

"""函数说明: 完整的线性SMO算法
Parameters: dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 maxIter - 最大迭代次数Returns: oS.b - SMO算法计算的b oS.alphas - SMO算法计算的alphas"""def smoP(dataMatIn, classLabels, C, toler, maxIter): # 初始化数据结构 oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler) iter = 0 # 初始化当前迭代次数 entireSet = True; alphaPairsChanged = 0  # 若遍历整个数据集都alpha也没有更新或者超过最大迭代次数,则退出循环 while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): 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] 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): entireSet = True # 如果alpha没有更新,计算全样本遍历 print("迭代次数: %d" % iter) return oS.b,oS.alphas

3.5 计算w

到上面，我们已经得到了b 和 alpha 的值，要得到超平面，还需要知道w的值。

"""函数说明：计算w
Parameters: dataArr - 数据矩阵 classLabels - 数据标签 alphas - alphas值Returns: w - 计算得到的w"""def 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) return w

3.6 使用算法

到这里，算法基本完成了，接下来使用书中提供的数据来测试一下算法，并可视化分类结果。

"""函数说明: 分类结果可视化
Parameters: dataMat - 数据矩阵 w - 直线法向量 b - 直线解决Returns: 无"""def showClassifer(dataMat, classLabels, w, b): import matplotlib.pyplot as plt
 #绘制样本点 data_plus = [] #正样本 data_minus = [] #负样本 for i in range(len(dataMat)): if classLabels[i] > 0: data_plus.append(dataMat[i]) else: data_minus.append(dataMat[i]) data_plus_np = np.array(data_plus) #转换为numpy矩阵 data_minus_np = np.array(data_minus) #转换为numpy矩阵 plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7) #正样本散点图 plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) #负样本散点图 #绘制直线 x1 = max(dataMat)[0] x2 = min(dataMat)[0] a1, a2 = w b = float(b) a1 = float(a1[0]) a2 = float(a2[0]) y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2 plt.plot([x1, x2], [y1, y2]) #找出支持向量点 for i, alpha in enumerate(alphas): if abs(alpha) > 0: x, y = dataMat[i] plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red') plt.show()
if __name__ == '__main__': dataArr, classLabels = loadDataSet('testSet.txt') b, alphas = smoP(dataArr, classLabels, 0.6, 0.001, 40) w = calcWs(alphas, dataArr, classLabels) showClassifer(dataArr, classLabels, w, b)

得到的分类效果图，其中红圈圈部分为支持向量。

4、4

4、使用核函数       

上面用到的数据是一个线性可分的，但对于分线性可分的，就要使用核函数将数据从低维映射到高维，然后再找到超平面。

要使用核函数，需要对上面的优化的SMO算法进行一些改动，包括其中的类和辅助支持函数，其 主要是加入了一个包含核函数信息的元组。而且，还需要有一个核转换函数以告诉程序使用哪种核函数方法。代码这里直接就都放在这里了，不再分开描述。

'''函数说明：核转换函数
Parameters: X - 原始数据 A - 参数 kTup - 包含核函数信息的元组Returns:    K - 转换后的数据矩阵'''def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space m,n = shape(X) K = mat(zeros((m,1))) if kTup[0]=='lin': K = X * A.T #linear kernel 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)) #divide in NumPy is element-wise not matrix like Matlab else: raise NameError('Houston We Have a Problem: That Kernel is not recognized') return K
"""类说明: 数据结构, 维护所有需要操作的值
Parameters: dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 kTup - 包含核函数信息的元组"""class optStructK: def __init__(self,dataMatIn, classLabels, C, toler, kTup): # Initialize the structure with the parameters  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)  """函数说明: 计算误差
Parameters: oS - 数据结构 k - 标号为k的数据Returns: Ek - 标号为k的数据的误差""" def calcEkK(oS, k): fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b) Ek = fXk - float(oS.labelMat[k]) return Ek
"""函数说明: 内循环启发方式2Parameters: i - 标号为i的数据的索引值 oS - 数据结构 Ei - 标号为i的数据误差Returns: j, maxK - 标号为j或maxK的数据的索引值 Ej - 标号为j的数据误差
"""def selectJK(i, oS, Ei): #this is the second choice -heurstic, and calcs Ej # 初始化 maxK = -1; maxDeltaE = 0; Ej = 0 # 根据Ei更新误差缓存 oS.eCache[i] = [1,Ei] #set valid #choose the alpha that gives the maximum delta E # 返回误差不为0的数据的索引值 validEcacheList = nonzero(oS.eCache[:,0].A)[0] if (len(validEcacheList)) > 1: for k in validEcacheList: #loop through valid Ecache values and find the one that maximizes delta E if k == i: continue #don't calc for i, waste of time Ek = calcEkK(oS, k) deltaE = abs(Ei - Ek) if (deltaE > maxDeltaE): maxK = k; maxDeltaE = deltaE; Ej = Ek return maxK, Ej else: #in this case (first time around) we don't have any valid eCache values j = selectJrand(i, oS.m) Ej = calcEkK(oS, j) return j, Ej
"""函数说明: 计算Ek, 并更新误差缓存
Parameters: oS - 数据结构 k - 标号为k的数据的索引值Returns: 无"""def updateEkK(oS, k):#after any alpha has changed update the new value in the cache Ek = calcEkK(oS, k) oS.eCache[k] = [1,Ek]    """函数说明: 优化的SMO算法
Parameters: i - 标号为i的数据的索引值 oS - 数据结构Returns: 1 - 有任意一对alpha值发生变化 0 - 没有任意一对alpha值发生变化或变化太小"""def innerLK(i, oS): # 步骤1: 计算误差Ei Ei = calcEkK(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 = selectJK(i, oS, Ei) #this has been changed from selectJrand # 保存更新前的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] #changed for kernel 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至误差缓存 updateEkK(oS, j) if (abs(oS.alphas[j] - alphaJold) < 0.00001):  print("j not moving enough") return 0  # 步骤6: 更新alpha_i oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j # 更新Ei至误差缓存 updateEkK(oS, i) #added this for the Ecache #the update is in the oppostie direction  # 步骤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    """函数说明: 完整的线性SMO算法
Parameters: dataMatIn - 数据矩阵 classLabels - 数据标签 C - 松弛变量 toler - 容错率 maxIter - 最大迭代次数 kTup - Returns: oS.b - SMO算法计算的b oS.alphas - SMO算法计算的alphas"""def smoPK(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #full Platt SMO # 初始化数据结构 oS = optStructK(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup) iter = 0 # 初始化当前迭代次数 entireSet = True; alphaPairsChanged = 0  # 若遍历整个数据集都alpha也没有更新或者超过最大迭代次数,则退出循环 while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): alphaPairsChanged = 0 if entireSet: # 遍历整个数据集 for i in range(oS.m):  alphaPairsChanged += innerLK(i,oS) #使用优化的SMO算法 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 # 如果alpha没有更新,计算全样本遍历  print("iteration number: %d" % iter) return oS.b,oS.alphas '''函数说明：测试使用测函数
Parameters: k1 - Returns: 无'''def testRbf(k1=1.3): dataArr,labelArr = loadDataSet('testSetRBF.txt') b,alphas = smoPK(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important datMat=mat(dataArr); labelMat = mat(labelArr).transpose() svInd=nonzero(alphas.A>0)[0] sVs=datMat[svInd] #get matrix of only support vectors 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,:],('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') errorCount = 0 datMat=mat(dataArr); labelMat = mat(labelArr).transpose() m,n = shape(datMat) 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 test error rate is: %f" % (float(errorCount)/m))  if __name__ == "__main__": testRbf()

    

5、手写数字识别       

这里使用之前的手写数字识别的数据来使用一下算法，由于支持向量机本质是一个二类分类器，所以这里只是用其中的 1 和 9 两个数字。其中的将图像转换成向量的函数和之前的一样。

def loadImages(dirName): from os import listdir hwLabels = [] trainingFileList = listdir(dirName) #load the training set m = len(trainingFileList) trainingMat = zeros((m,1024)) for i in range(m): fileNameStr = trainingFileList[i] fileStr = fileNameStr.split('.')[0] #take off .txt 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 
'''函数说明:将32x32的二进制图像转换称1x1024向量
Parameters: filename - 文件名Returns: returnVect - 返回的二进制图像的1x1024向量'''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 testDigits(kTup=('rbf', 10)): dataArr,labelArr = loadImages('trainingDigits') b,alphas = smoPK(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)) 
if __name__ == "__main__": testDigits()

最终测试的结果：

总结放到下一篇文章里吧，这里使用sklearn才是最舒服的，也是我最想写的。