手撕SVM（公式推导代码实现）（三）
前言
前面我们进行了很多的理论性研究，下面我们开始用代码进行实现。编程求解线性SVM可视化数据集import matplotlib.pyplot as plt import numpy as np #读取数据 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 #数据可视化 def showDataSet(dataMat,labelMat): data_plus = [] data_minus = [] for i in range(len(dataMat)): if labelMat[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]) #正样本 plt.scatter(np.transpose(data_minus_np)[0],np.transpose(data_minus_np)[1]) #负样本 plt.show() dataMat,labelMat = loadDataSet(＂testSet.txt＂) showDataSet(dataMat,labelMat)
这个数据集显然线性可分。 应用简化版SMO算法处理小规模数据集
import random def selectJrand(i,m): j=i while(j==i): #选择一个不等于i的j j = int(random.uniform(0,m)) return j def clipAlpha(aj,H,L): if aj > H: aj = H if L > aj: aj = L return aj dataArr,labelArr =loadDataSet(＂testSet.txt＂) labelArr
[-1.0,
-1.0,
1.0,
-1.0,
1.0,
1.0,
1.0,
-1.0,
-1.0,
-1.0,
-1.0,
-1.0,
-1.0,
1.0,
-1.0,
1.0,
1.0,
-1.0,
1.0,
-1.0,
-1.0,
-1.0,
1.0,
-1.0,
-1.0,
1.0,
1.0,
-1.0,
-1.0,
-1.0,
-1.0,
1.0,
1.0,
1.0,
1.0,
-1.0,
1.0,
-1.0,
-1.0,
1.0,
-1.0,
-1.0,
-1.0,
-1.0,
1.0,
1.0,
1.0,
1.0,
1.0,
-1.0,
1.0,
1.0,
-1.0,
-1.0,
1.0,
1.0,
-1.0,
1.0,
-1.0,
-1.0,
-1.0,
-1.0,
1.0,
-1.0,
1.0,
-1.0,
-1.0,
1.0,
1.0,
1.0,
-1.0,
1.0,
1.0,
-1.0,
-1.0,
1.0,
-1.0,
1.0,
1.0,
1.0,
1.0,
1.0,
1.0,
1.0,
-1.0,
-1.0,
-1.0,
-1.0,
1.0,
-1.0,
1.0,
1.0,
1.0,
-1.0,
-1.0,
-1.0,
-1.0,
-1.0,
-1.0,
-1.0] data = np.mat(dataArr) data[2,:]matrix([[ 7.55151, -1.58003]])
可以看出来，这里使用的类别标签是-1和1
SMO算法的伪代码： 创建一个alpha向量并将其初始化为 0 向量 当迭代次数小于最大迭代次数时 (外循环) 对数据集中的每个数据向量 (内循环): 如果该数据向量可以被优化: 随机选择另外一个数据向量 同时优化这两个向量 如果两个向量都不能被优化, 退出内循环 如果所有向量都没被优化, 增加迭代数目, 继续下一次循环#简化版SMO算法 def smoSimple(dataMatIn,classLabels,C,toler,maxIter): dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose() b = 0; m,n = np.shape(dataMatrix) alphas = np.mat(np.zeros((m,1)))#初始化alpha参数，设置为0 iterSmo = 0 #初始化迭代次数 while(iterSmo < maxIter): alphaPairsChanged = 0#用于记录alpha是否已经进行优化 #步骤1. 计算误差Ei for i in range(m): fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b Ei = fXi - float(labelMat[i]) if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > -toler) and (alphas[i] > 0)): j = selectJrand(i,m) #步骤1. 计算误差Ej fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b Ej = fXj - float(labelMat[j]) #保存更新前的alpha值，使用浅拷贝 alphaIold = alphas[i].copy() alphaJold = alphas[j].copy() #步骤2：计算上界H和下界L 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 #步骤5：修剪alpha_j alphas[j] = clipAlpha(alphas[j],H,L) if (abs(alphas[j] - alphaJold) < 0.00001) : print(＂j not moving enough＂) ; continue #步骤6：更新alpha_i alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) #步骤7：更新b_1和b_2 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 #步骤8：根据b_1和b_2更新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 #统计优化次数 #如果程序执行到for循环的最后一行都不执行continue语句，那么就已经成功地改变了一对alpha，同时可以增加alphaPairsChanged的值 alphaPairsChanged += 1 #打印统计信息 print(＂第%d次迭代 样本:%d, alpha优化次数:%d＂ % (iterSmo,i,alphaPairsChanged)) #更新迭代次数 #在for循环之外，需要检查alpha值是否做了更新，如果有更新则将iterSmo设为0后继续运行程序。只有在所有数据集上遍历maxIter次，且不再发生任何alpha修改之后，程序才会停止并退出while循环 if (alphaPairsChanged == 0): iterSmo += 1 else: iterSmo = 0 print(＂迭代次数：%d＂ % iterSmo) return b,alphas b,alphas = smoSimple(dataArr,labelArr,0.6,0.001,500)
L==H
第0次迭代 样本:1, alpha优化次数:1
第0次迭代 样本:3, alpha优化次数:2
第0次迭代 样本:5, alpha优化次数:3
L==H
第0次迭代 样本:8, alpha优化次数:4
L==H
j not moving enough
j not moving enough
L==H
L==H
j not moving enough
L==H
第0次迭代 样本:30, alpha优化次数:5
第0次迭代 样本:31, alpha优化次数:6
L==H
L==H
第0次迭代 样本:54, alpha优化次数:7
L==H
L==H
第0次迭代 样本:71, alpha优化次数:8
L==H
L==H
L==H
第0次迭代 样本:79, alpha优化次数:9
L==H
第0次迭代 样本:92, alpha优化次数:10
j not moving enough
L==H
迭代次数：0
第0次迭代 样本:1, alpha优化次数:1
j not moving enough
j not moving enough
j not moving enough
j not moving enough
j not moving enough
L==H
L==H
j not moving enough
j not moving enough
j not moving enough
j not moving enough
L==H
j not moving enough
L==H
L==H
j not moving enough
j not moving enough
第0次迭代 样本:37, alpha优化次数:2
第0次迭代 样本:39, alpha优化次数:3
第0次迭代 样本:52, alpha优化次数:4
j not moving enough
j not moving enough
j not moving enough
j not moving enough
j not moving enough
第0次迭代 样本:71, alpha优化次数:5
j not moving enough
j not moving enough
j not moving enough
j not moving enough
j not moving enough
迭代次数：0
j not moving enough
j not moving enough
j not moving enough
第0次迭代 样本:8, alpha优化次数:1
L==H
j not moving enough
第0次迭代 样本:23, alpha优化次数:2
L==H
j not moving enough
j not moving enough
L==H
j not moving enough
j not moving enough
j not moving enough
第0次迭代 样本:39, alpha优化次数:3
L==H
j not moving enough
第0次迭代 样本:52, alpha优化次数:4
j not moving enough
第0次迭代 样本:55, alpha优化次数:5
L==H
L==H
L==H
L==H
L==H
j not moving enough
第0次迭代 样本:79, alpha优化次数:6
第0次迭代 样本:92, alpha优化次数:7
迭代次数：0
j not moving enough
L==H
j not moving enough
j not moving enough
L==H
j not moving enough
第0次迭代 样本:23, alpha优化次数:1
j not moving enough
j not moving enough
j not moving enough
j not moving enough
j not moving enough
j not moving enough
j not moving enough
L==H
L==H
第0次迭代 样本:51, alpha优化次数:2
j not moving enough
j not moving enough
j not moving enough
j not moving enough
L==H
第0次迭代 样本:69, alpha优化次数:3
L==H
j not moving enough
第0次迭代 样本:94, alpha优化次数:4
j not moving enough
j not moving enough
迭代次数：0
j not moving enough
j not moving enough
j not moving enough
j not moving enough
j not moving enough
...
迭代次数：497
j not moving enough
j not moving enough
j not moving enough
迭代次数：498
j not moving enough
j not moving enough
j not moving enough
迭代次数：499
j not moving enough
j not moving enough
j not moving enough
迭代次数：500 b matrix([[-3.83785102]]) alphas[alphas>0] matrix([[0.1273855 , 0.24131542, 0.36872064]]) np.shape(alphas[alphas>0]) (1, 3) for i in range(100): if alphas[i] > 0.0: print(dataArr[i],labelArr[i]) [4.658191, 3.507396] -1.0 [3.457096, -0.082216] -1.0 [6.080573, 0.418886] 1.0 #分类结果可视化 def showClassifer(dataMat, w, b): #绘制样本点 data_plus = [] #正样本 data_minus = [] #负样本 for i in range(len(dataMat)): if labelMat[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() #计算w def get_w(dataMat, labelMat, alphas): alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat) w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas) return w.tolist() w = get_w(dataMat,labelMat,alphas) showClassifer(dataMat,w,b)