前言

前面我們進(jìn)行了很多的理論性研究，下面我們開(kāi)始用代碼進(jìn)行實(shí)現(xiàn)。

編程求解線(xiàn)性SVM

• 可視化數(shù)據(jù)集

import matplotlib.pyplot as pltimport numpy as np#讀取數(shù)據(jù)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])]) #添加數(shù)據(jù) labelMat.append(float(lineArr[2])) #添加標(biāo)簽 return dataMat,labelMat#數(shù)據(jù)可視化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) #轉(zhuǎn)換為numpy矩陣 data_minus_np = np.array(data_minus)#轉(zhuǎn)換為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]) #負(fù)樣本 plt.show()dataMat,labelMat = loadDataSet(‘testSet.txt’)showDataSet(dataMat,labelMat)

這個(gè)數(shù)據(jù)集顯然線(xiàn)性可分。

應(yīng)用簡(jiǎn)化版SMO算法處理小規(guī)模數(shù)據(jù)集

import randomdef selectJrand(i,m): j=i while(j==i): #選擇一個(gè)不等于i的j j = int(random.uniform(0,m)) return jdef clipAlpha(aj,H,L): if aj > H: aj = H if L > aj: aj = L return ajdataArr,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]])

可以看出來(lái)，這里使用的類(lèi)別標(biāo)簽是-1和1

SMO算法的偽代碼：

創(chuàng)建一個(gè)alpha向量并將其初始化為 0 向量 當(dāng)?shù)?span id="zq7hqw5" class="wpcom_tag_link">次數(shù)小于最大迭代次數(shù)時(shí) (外循環(huán)) 對(duì)數(shù)據(jù)集中的每個(gè)數(shù)據(jù)向量 (內(nèi)循環(huán)): 如果該數(shù)據(jù)向量可以被優(yōu)化: 隨機(jī)選擇另外一個(gè)數(shù)據(jù)向量 同時(shí)優(yōu)化這兩個(gè)向量 如果兩個(gè)向量都不能被優(yōu)化, 退出內(nèi)循環(huán)如果所有向量都沒(méi)被優(yōu)化, 增加迭代數(shù)目, 繼續(xù)下一次循環(huán)#簡(jiǎn)化版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參數(shù)，設(shè)置為0 iterSmo = 0 #初始化迭代次數(shù) while(iterSmo < maxIter): alphaPairsChanged = 0#用于記錄alpha是否已經(jīng)進(jìn)行優(yōu)化 #步驟1. 計(jì)算誤差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 0)): j = selectJrand(i,m) #步驟1. 計(jì)算誤差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：計(jì)算上界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：計(jì)算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：根據(jù)b_1和b_2更新b if (0 alphas[i]) : b = b1 elif (0 alphas[j]): b = b2 else: b = (b1 + b2)/2.0 #統(tǒng)計(jì)優(yōu)化次數(shù) #如果程序執(zhí)行到for循環(huán)的最后一行都不執(zhí)行continue語(yǔ)句，那么就已經(jīng)成功地改變了一對(duì)alpha，同時(shí)可以增加alphaPairsChanged的值 alphaPairsChanged += 1 #打印統(tǒng)計(jì)信息 print("第%d次迭代 樣本:%d, alpha優(yōu)化次數(shù):%d" % (iterSmo,i,alphaPairsChanged)) #更新迭代次數(shù) #在for循環(huán)之外，需要檢查alpha值是否做了更新，如果有更新則將iterSmo設(shè)為0后繼續(xù)運(yùn)行程序。只有在所有數(shù)據(jù)集上遍歷maxIter次，且不再發(fā)生任何alpha修改之后，程序才會(huì)停止并退出while循環(huán) if (alphaPairsChanged == 0): iterSmo += 1 else: iterSmo = 0 print("迭代次數(shù)：%d" % iterSmo) return b,alphasb,alphas = smoSimple(dataArr,labelArr,0.6,0.001,500)

L==H第0次迭代 樣本:1, alpha優(yōu)化次數(shù):1第0次迭代 樣本:3, alpha優(yōu)化次數(shù):2第0次迭代 樣本:5, alpha優(yōu)化次數(shù):3L==H第0次迭代 樣本:8, alpha優(yōu)化次數(shù):4L==Hj not moving enoughj not moving enoughL==HL==Hj not moving enoughL==H第0次迭代 樣本:30, alpha優(yōu)化次數(shù):5第0次迭代 樣本:31, alpha優(yōu)化次數(shù):6L==HL==H第0次迭代 樣本:54, alpha優(yōu)化次數(shù):7L==HL==H第0次迭代 樣本:71, alpha優(yōu)化次數(shù):8L==HL==HL==H第0次迭代 樣本:79, alpha優(yōu)化次數(shù):9L==H第0次迭代 樣本:92, alpha優(yōu)化次數(shù):10j not moving enoughL==H迭代次數(shù)：0第0次迭代 樣本:1, alpha優(yōu)化次數(shù):1j not moving enoughj not moving enoughj not moving enoughj not moving enoughj not moving enoughL==HL==Hj not moving enoughj not moving enoughj not moving enoughj not moving enoughL==Hj not moving enoughL==HL==Hj not moving enoughj not moving enough第0次迭代 樣本:37, alpha優(yōu)化次數(shù):2第0次迭代 樣本:39, alpha優(yōu)化次數(shù):3第0次迭代 樣本:52, alpha優(yōu)化次數(shù):4j not moving enoughj not moving enoughj not moving enoughj not moving enoughj not moving enough第0次迭代 樣本:71, alpha優(yōu)化次數(shù):5j not moving enoughj not moving enoughj not moving enoughj not moving enoughj not moving enough迭代次數(shù)：0j not moving enoughj not moving enoughj not moving enough第0次迭代 樣本:8, alpha優(yōu)化次數(shù):1L==Hj not moving enough第0次迭代 樣本:23, alpha優(yōu)化次數(shù):2L==Hj not moving enoughj not moving enoughL==Hj not moving enoughj not moving enoughj not moving enough第0次迭代 樣本:39, alpha優(yōu)化次數(shù):3L==Hj not moving enough第0次迭代 樣本:52, alpha優(yōu)化次數(shù):4j not moving enough第0次迭代 樣本:55, alpha優(yōu)化次數(shù):5L==HL==HL==HL==HL==Hj not moving enough第0次迭代 樣本:79, alpha優(yōu)化次數(shù):6第0次迭代 樣本:92, alpha優(yōu)化次數(shù):7迭代次數(shù)：0j not moving enoughL==Hj not moving enoughj not moving enoughL==Hj not moving enough第0次迭代 樣本:23, alpha優(yōu)化次數(shù):1j not moving enoughj not moving enoughj not moving enoughj not moving enoughj not moving enoughj not moving enoughj not moving enoughL==HL==H第0次迭代 樣本:51, alpha優(yōu)化次數(shù):2j not moving enoughj not moving enoughj not moving enoughj not moving enoughL==H第0次迭代 樣本:69, alpha優(yōu)化次數(shù):3L==Hj not moving enough第0次迭代 樣本:94, alpha優(yōu)化次數(shù):4j not moving enoughj not moving enough迭代次數(shù)：0j not moving enoughj not moving enoughj not moving enoughj not moving enoughj not moving enough

…

迭代次數(shù)：497j not moving enoughj not moving enoughj not moving enough迭代次數(shù)：498j not moving enoughj not moving enoughj not moving enough迭代次數(shù)：499j not moving enoughj not moving enoughj not moving enough迭代次數(shù)：500

bmatrix([[-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#分類(lèi)結(jié)果可視化def showClassifer(dataMat, w, b): #繪制樣本點(diǎn) data_plus = [] #正樣本 data_minus = [] #負(fù)樣本 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) #轉(zhuǎn)換為numpy矩陣 data_minus_np = np.array(data_minus) #轉(zhuǎn)換為numpy矩陣 plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7) #正樣本散點(diǎn)圖 plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) #負(fù)樣本散點(diǎn)圖 #繪制直線(xiàn) 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]) #找出支持向量點(diǎn) 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()#計(jì)算wdef 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)