經過前幾篇文章的學習，SVM的優化目標，SMO演算法的基本實現步驟，模型對應引數的選擇，我們已經都有了一定的理解，結合《機器學習實戰》，動手實踐一個基本的SVM支援向量機，來完成一個簡單的二分類任務。

建立模型之前，首先看一下我們的資料，然後再用支援向量機實現分類：

這裡只擷取部分，我們的資料是一個txt檔案，包含100個數據點，前兩行為x1，x2，負一行為標籤值，取值{-1,+1}.

讀取檔案資料

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,labelMat#資料集，標籤值
 

LoadDataSet函式的作用是讀取資料檔案，return得到兩個列表，一個包含二維資料點，一個包含點對應的標籤值

隨機選擇j

import random
def selectJrand(i,m):
    j=i 
    while (j==i):
        j = int(random.uniform(0,m))
    return j

selectJrand函式通過while函式判斷 j==i，從而達到隨機生成不等於 i 的 j ，因為這裡是簡化的SVM演算法，如果在複雜情況下，j的選擇會和這裡j所的選擇方法有所不同

調整alpha

def clipAlpha(aj,H,L):
    if aj > H: 
        aj = H
    if L > aj:
        aj = L
    return aj
 

clipAlpha函式根據（4）的alpha的界限，調整迭代得到的alpha

簡化版SMO演算法

import numpy as np

def smoSimple(dataMatIn,classLabels,C,toler,maxIter):
	#引數分別為資料集，類別標籤，常數c，容錯率和最大迴圈次數
	#外迴圈
	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
	iter = 0#初始化迭代數
	while (iter < maxIter):
		#內迴圈
		alphaParisChanged = 0
		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)):#判斷alpha是否滿足KKT條件
				j = selectJrand(i,m)#隨機選取j
				fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b #計算預測類別
				Ej = fXj - float(labelMat[j])#計算誤差e
				alphaIold = alphas[i].copy()#為alphaIold和alphaJold分配新的記憶體
				alphaJold = alphas[j].copy()
				if (labelMat[i] != labelMat[j]):#alpha的上限與下限
					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
				eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
				#計算smo演算法中的η
				if eta > 0 : print("eta>=0");continue

				#對i進行修改，修改量與j相同，但方向相反
				alphas[j] -= labelMat[j]*(Ei-Ej)/eta
				#通過L H的值優化
				alphas[j] = clipAlpha(alphas[j],H,L)
                                #判斷alphaJ的優化程度，更新alphaI，並根據新的alpha得到新的b
				if (abs(alphas[j]-alphaJold) < 0.0001):
					print("j not moving enough");continue
				alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[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[i]):
					b = b1
				elif (0 < alphas[j]) and (C > alphas[j]):
					b = b2
				else:
					b = (b1 + b2) / 2.0
				alphaParisChanged += 1
				print("iter : %d i:%d,pairs changed %d" %(iter,i,alphaParisChanged))
		if (alphaParisChanged == 0):#判斷是否對alpha進行了更新
			iter += 1
		else:
			iter = 0
		print("iterration number :%d" % iter)
	return b,alphas#返回超平面需要的引數alpha，b
 

SmoSimple函式是smo演算法的簡單實現版本，與原始演算法思想基本一致，需要計算的是誤差 ei = f(xi)-yi，二次項η (eta)，以及alpha的上界H下界L，需要判斷的是開頭的KKT條件是否滿足，判斷yi與yj的類別是否一致從而決定上下界的形式，判斷alpha的更新程度以及值的大小從而決定b，需要更新的是alpha的值，b的值，還有演算法迭代次數，最後返回我們超平面構建需要的引數alpha，b . Tips： KKT條件的判斷與西瓜書，《統計學習方法》有所不同，在（3）中解釋了程式碼中的KKT條件是如何判斷的.

繪製樣本資料點，超平面與支援向量

import matplotlib.pyplot as plt

def plot_point(dataMat,labelMat,alphas,b):

	for i in range(np.shape(dataMat)[0]):#繪製shape(dataMat)[0]個數據點
		plt.scatter(dataMat[i][0],dataMat[i][1])

	alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)#轉為陣列求w
	sum = 0#初始化w
	for i in range(np.shape(dataMat)[0]):
		sum += alphas[i]*labelMat[i]*dataMat[i].T#通過求導後得到的w公式，求和得到w
	# print(sum)

	for i, alpha in enumerate(alphas):#根據kkt條件，選取alpha不為0的點作為支援向量
	    if abs(alpha) > 0:
	        x, y = dataMat[i]
	        plt.scatter(x, y, s=100, c = '', alpha=0.5, linewidth=1.5, edgecolor='red')

	x = np.arange(0,10,0.01)
	y = (sum[0]*x+b)/(-1*sum[1])#通過原始公式推匯出超平面

	plt.scatter(x,y,s=5,marker = 'h')
	
	plt.show()

plot_point函式是負責繪製超平面的，書上只給了示意圖，所以就只能自己動手豐衣足食了，首先繪製了所有資料點，由於資料點比較標準(無明顯越界)，所以為了好看就沒有根據類別設定顏色，其次根據對偶問題得到w的計算公式

通過SmoSimple返回的alpha計算w，再加上b，就得到了最優分割超平面的公式，由於是二維資料點，所以稍加整理就可以得到超平面的二維表示式，然後我們設定適當的x範圍，繪製超平面就ok了；還有一點就是支援向量的標註，因為alpha≠0時，(wi.T).x+b=0，對應的點在分割邊界上，所以對應的點即為支援向量Support Vector，最後用c=' '空心標註出支援向量，就能得到書中的圖了.

主函式

if __name__ == '__main__':
	dataMat, labelMat = loadDataSet('testSet.txt')
	b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
	plot_point(dataMat,labelMat,alphas,b)

可以自己加print語句，看一下w，b：

b =  [[-3.86378686]]
w =  [ 0.79293619 -0.23660327]
[Finished in 12.9s]

然後看一下得到的分隔超平面：

總結：

第一次看書上的程式碼會有些懵，不過看一遍SMO演算法推導，程式碼部分就基本無壓力了，這裡運用的是簡化版演算法，用時較長，不過得到的結果比較滿意，之後還會用完整的SMO演算法和Sklearn實現SVM，執行速度會大大提升，同時結果也比較滿意。