前面已經對感知機和SVM進行了簡要的概述，本節是SVM演算法的實現過程用於輔助理解SVM演算法的具體內容，然後藉助sklearn對SVM工具包進行實現。
SVM演算法的核心是SMO演算法的實現，首先對SMO演算法過程進行實現，先對一些輔助函式進行定義：
1 # 先定義一些輔助函式
2 # 選取第二變數函式
3 def select_J_rand(i, m):
4 j=i
5 while(j==i):
6 j = int(random.uniform(0, m))
7 return j
8
9 # 定義對α進行裁剪的函式
10 def clip_alpha(aj, H, L):
11 if aj > H:
12 aj=H
13 if L > aj:
14 aj = L
15 return aj
然後實現一個簡化版的SMO演算法：
"""
Input：dataX, dataY, C(懲罰因子), toler（容忍度）, iter_num
Output: alpha、b
"""
def smo_simple(dataX, dataY, C, toler, iter_num):
dataMatrix = mat(dataX); labelMat = dataY.transpose()
# 初始化引數
b = 0; m, n = np.shape(dataMatrix)
alphas = mat(np.zeros((m, 1)))
iter = 0
# 當超過迭代次數停止
while iter < iter_num:
# 記錄α1和α2變化次數
alphaPairsChanged = 0
for i in range(m):
# 計算f(xi),預測類別
fXi = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
# 計算誤差
Ei = fXi - float(labelMat[i])
# 當不滿足條件時，選取變數j，這裡要判斷α是否在[0,C]內，若超出範圍則不再進行優化
if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and alphas[i] > 0):
j = select_J_rand(i, m)
# 計算x2的預測值y2
fXj = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
Ej = fXj - float(labelMat[j])
alpha_I_old = alphas[i].copy()
alpha_J_old = alphas[j].copy()
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[i] + alphas[j] - C)
H = min(C, alphas[i] + alphas[j])
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
if eta > 0:
print("eta > 0")
continue
alphas[j] -= labelMat[j] * (Ej - Ei)/eta
alphas[j] = clip_alpha(alphas[j], H, L)
# 當α2不再變化
if (abs(alphas[j]-alpha_J_old) < 0.00001):
print("j not moving enough")
continue
# 更新α1
alphas[i] += labelMat[i] * labelMat[j] * (alpha_J_old - alphas[j])
# 計算b1和b2
b1 = b - Ei - labelMat[i] * (alphas[i] - alpha_I_old) * dataMatrix[i, :] * dataMatrix[i, :].T - labelMat[j] * (alphas[j] - alpha_J_old) * dataMatrix[i, :] * dataMatrix[j, :].T
b2 = b - Ej - labelMat[i] * (alphas[i] - alpha_I_old) * dataMatrix[i, :] * dataMatrix[j, :].T - labelMat[j] * (alphas[j] - alpha_J_old) * dataMatrix[j, :] * dataMatrix[j, :].T
if (alphas[i] > 0) and (alphas[i] < C):
b = b1
elif (alphas[j] > 0) and (alphas[j] < C):
b = b2
else:
b = (b1 + b2)/2
alphaPairsChanged += 1
if alphaPairsChanged == 0:
iter += 1
else:
iter = 0
print("iteration number: %d"%iter)
return b, alphas
SMO演算法具有一定的隨機性，因此每次執行的結果不一定相同。上面就是一個簡單的SMO演算法的實現部分，對於小批量資料可以滿足需求，但當資料量過於龐大時，上面的演算法的效率將會很慢，這是因為在α的選擇問題上，下面提供一種改進的SMO演算法，改進的SMO演算法會通過一個外迴圈選擇第一個α的值，選擇方法是在遍歷所有樣本（資料集）和非邊界α中進行掃描，所謂非邊界α是指那些不等於0或者C的α值，建立這些α值的列表，在列表中進行遍歷，並在掃描時跳過不再改變的α進行遍歷。下面是具體實現過程
首先定義輔助函式用於儲存和更新E，並建立一個數據結構儲存變數
1 # 首先建立3個輔助函式用於對E進行快取，以及1種用於儲存資料的資料結構
2 # 儲存變數的資料結構
3 class optStruct:
4 def __init__(self, dataX, dataY, C, toler):
5 self.X = dataX
6 self.Y = dataY
7 self.C = C
8 self.toler = toler
9 self.m = np.shape(dataX)[0]
10 self.alphas = np.mat(zeros((self.m, 1)))
11 self.b = 0
12 # cache第一列為有效性標誌位，第二列為E值
13 self.eCache = np.mat(np.zeros((self.m, 2)))
14
15 # 計算E值並返回，由於頻繁使用單獨寫成一個函式
16 def calcEk(oS, k):
17 fXk = float(np.multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k, :].T)) + oS.b
18 Ek = fXk - float(oS.labelMat[k])
19 return Ek
20
21 # 用於選擇第二個α的值，保證每次優化採用最大的步長
22 def select_J(i, oS, Ei):
23 maxK = -1; maxDeltaE = 0; Ej = 0
24 oS.eCache[i] = [1, Ei]
25 validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
26 if len(validEcacheList) > 1:
27 for k in validEcacheList:
28 if k == i:
29 continue
30 Ek = calcEk(oS, k)
31 deltaE = abs(Ei - Ek)
32 # 選擇變化最大的那個
33 if deltaE > maxDeltaE:
34 maxK = k
35 maxDeltaE = deltaE
36 Ej = Ek
37 return maxK, Ej
38 else:
39 j = select_J_rand(i, oS.m)
40 Ej = calcEk(oS, j)
41 return j, Ej
42
43
44 def updateEk(oS, k):
45 Ek = calcEk(oS, k)
46 oS.eCache[k] = [1, Ek]
接下來就是SMO演算法的改進版本
1 # 與simpleSMO一致，更新的alpha存入cache中
2 def innerL(i, oS):
3 Ei = calcEk(oS, i)
4 if ((oS.labelMat[i] * Ei < -oS.toler) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i] * Ei > oS.toler) and (oS.alphas[i] > 0)):
5 j, Ej = select_J(i, oS, Ei)
6 alpha_I_old = oS.alphas[i].copy()
7 alpha_J_old = oS.alphas[j].copy()
8 if oS.labelMat[i] != oS.labelMat[j]:
9 L = max(0, oS.alphas[j] - oS.alphas[i])
10 H = min(oS.C, oS.C+ oS.alphas[j] - oS.alphas[i])
11 else:
12 L = min(0, oS.alphas[i] + oS.alphas[j] - oS.C)
13 H = max(oS.C, oS.alphas[i] + oS.alphas[j])
14 if H == L:
15 return 0
16 eta = 2.0 * oS.X[i, :] * oS.X[j, :].T - oS.X[i, :] * oS.X[i, :].T - oS.X[j, :] * oS.X[j, :].T
17 if eta >= 0:
18 return 0
19 oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta
20 oS.alphas[j] = clip_alpha(oS.alphas[j], H, L)
21 updateEk(oS, j)
22 if abs(oS.alphas[j] - alpha_J_old) < 0.00001:
23 return 0
24 oS.alphas[i] -= oS.labelMat[i] * oS.labelMat[j] * (oS.alphas[j] - alpha_J_old)
25 updateEk(oS, i)
26 b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alpha_I_old) * oS.X[i, :] * oS.X[i, :].T - oS.labelMat[j] * (oS.alphas[j] - alpha_J_old) * oS.X[i, :] * oS.X[j, :].T
27 b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alpha_I_old) * oS.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * (oS.alphas[j] - alpha_J_old) * oS.X[j, :] * oS.X[j, :].T
28 if oS.alphas[i] > 0 and oS.alphas[i] < oS.C:
29 oS.b = b1
30 elif oS.alphas[j] > 0 and oS.alphas[j] < oS.C:
31 oS.b = b2
32 else:
33 os.b = (b1 + b2)/2
34 return 1
35 else:
36 return 0
37
38
39 def smoP(dataX, labelMat, C, toler, maxIter, kTup=('lin', 0)):
40 oS = optStruct(mat(dataX), mat(labelMat).transpose(), C, toler)
41 iter = 0
42 entireSet = True
43 alphaPairsChanged = 0
44 while (iter < maxIter) and alphaPairsChanged > 0 or entireSet:
45 alphaPairsChanged = 0
# 搜尋第一個變數的值，採用兩個方法交替進行的方式，利用entireSet變數控制
46 # 第一種遍歷全體樣本
47 if entireSet:
48 for i in range(oS.m):
49 alphaPairsChanged += innerL(i, oS)
50 iter += 1
51 # 第二種遍歷非邊界樣本
52 else:
53 nonBoundIs = nonzero((oS.alphas.A > 0) * oS.alphas.A < C)[0]
54 for i in nonBoundIs:
55 alphaPairsChanged += innerL(i, oS)
56 iter += 1
57 if entireSet:
58 entireSet = False
59 elif alphaPairsChanged == 0:
60 entireSet = True
61 return oS.alphas, oS.b
獲取到α的值後，則可以進一步求出模型的w和b，具體實現過程為:
1 def calW(alphas, dataArr, classLabels):
2 X = mat(dataArr)
3 labelMat = mat(classLabels).transpose()
4 m, n = shape(X)
5 w = zeros((m, 1))
6 for i in range(m):
7 w += multiply(alphas[i] * labelMat[i], X[i, :])
8 return w
上述即為SVM的實現過程，是對線性可分的資料的分類的實現過程，當對於非線性資料，需要運用到核函式，在實現過程中也較為簡單，只需只需將（xi·xj）替換為核函式即可，具體實現過程如下
1 # 首先原先的資料結構中要計算核矩陣，將下述程式碼加入資料結構程式碼部分即可
2 self.K = mat(zeros(self.m, self.m))
3 for i in range(self.m):
4 self.K[:, i] = kernelTrans(self.X, self.X[i, :], self.kTup)
接下來實現核轉換函式
1 def kernelTrans(X, A, kTup):
2 m, n = shape(X)
3 K = mat(zeros((m, 1)))
4 if kTup['0'] == 'lin':
5 K = X * A.T
6 elif kTup[0] == 'rbf':
7 for j in range(m):
8 deltaRow = X[j, :] - A
9 K[j] = deltaRow * deltaRow.T
10 K = exp(K/(-1 * kTup[1] ** 2))
11 else:
12 raise NameError('沒有定義核函式')
13 return K
然後需要在引數計算的函式中將對應的xi*xj替換為核函式即可：
1 # 首先是innerL中
2 eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
3 b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alpha_I_old) * oS.K[i, i] - oS.labelMat[j] * (oS.alphas[j] - alpha_J_old) * oS.K[i, j]
4 b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alpha_I_old) * oS.K[i, j] - oS.labelMat[j] * (oS.alphas[j] - alpha_J_old) * oS.K[j, j]
5
6 # 然後是calcEk
7 fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b)
至此，SVM演算法的實現過程基本已經完成了，接下來我們利用MNIST的資料集，對手寫數字進行分類和辨識，MNIST手寫辨識資料在網上就可以找到。
首先需要寫一個讀取資料的函式以及將函式影象轉化為矩陣的函式：
1 def img2Vector(filename):
2 returnVec = zeros((1, 1024))
3 fr = open(filename)
4 for i in range(32):
5 lineStr = fr.readline()
6 for j in range(32):
7 returnVec[0, 32*i + j] = int(lineStr[j])
8 return returnVec
1 def loadImages(dir):
2 hwLabels = []
3 trainingFileList = os.listdir(dir)
4 m = len(trainingFileList)
5 for i in range(m):
6 fileNameStr = trainingFileList[i]
7 fileStr = fileNameStr.split('.')[0]
8 classNumStr = int(fileStr.split('_')[0])
# 只考慮二分類問題
9 if classNumStr == 9:
10 hwLabels.append(-1)
11 else:
12 hwLabels.append(1)
13 trainingMat[i, :] = img2Vector('%s/%s' % (dir, fileNameStr))
14 return trainingMat, hwLabels
然後編寫主程式，用於分類和測試。
1 def testDigits(kTup=('rbf', 10)):
# 獲取資料
2 dataArr, labelArr = loadImages('trainingDigits')
# 利用SMO演算法求解出α和b
3 alphas, b = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
4 dataMat = mat(dataArr)
5 labelMat = mat(labelArr).transpose()
# 獲取支援向量的索引
6 svInd = nonzero(alphas.A > 0)[0]
# 獲取支援向量
7 svs = dataMat[svInd]
8 labelSv = labelMat[svInd]
9 m, n = shape(dataMat)
10 errorCount = 0
11 for i in range(m):
12 kernelEval = kernelTrans(svs, dataMat[i, :], kTup)
13 predict = kernelEval.T * multiply(labelSv, alphas[svInd]) + b
14 if sign(predict) != sign(labelArr[i]):
15 errorCount += 1
16 print('there are %d Support Vectors'%shape(svs)[0])
17 print('the error rate is %f' % (errorCount / (len(labelMat))))
18 test_dataArr, test_labelArr = loadImages('testDigits')
19 test_dataMat = mat(test_dataArr)
20 test_labelMat = mat(test_labelArr).transpose()
21 m1, n1 = shape(test_dataMat)
22 test_errorCount = 0
23 for i in range(m1):
24 kernelEval = kernelTrans(svs, test_dataMat[i, :], kTup)
25 predict = kernelEval.T * multiply(labelSv, alphas[svInd]) + b
26 if sign(predict) != sign(test_labelArr[i]):
27 errorCount += 1
28 print('the error rate is %f' % (test_errorCount / (len(test_labelMat))))
輸出結果如下：
至此，上述是一個完整的SVM演算法用於二分類的實現過程，核函式可以進行修改和替換，同時，對於多分類情況相當於建立多個分類器進行分類，過程這裡不再贅述，接下來，使用python中sklearn包對Mnist的資料集分類實現一遍。
首先是從sklearn匯入svm包，並讀取資料
1 from sklearn import svm
2 train_dataArr, train_labelArr = loadImages(dir)
然後是模型的初始化，這裡模型中有一些引數，其具體說明如下
1 model = svm.SVC(C=200, kernel='rbf', tol=1e-4, max_iter=-1, degree=3, gamma='auto_deprecated', coef0=0, shrinking=True,
2 probability=False, cache_size=200, verbose=False, class_weight=None, decision_function_shape='ovr')
3 """
4 C: 懲罰因子，預設為1.0，C越大，相當於懲罰鬆弛變數，希望鬆弛變數接近0，即對誤分類的懲罰增大，趨向於對訓練集全分對的情況，這樣對訓練集測試時準確率很高，但泛化能力弱。C值小，對誤分類的懲罰減小，允許容錯，將他們當成噪聲點，泛化能力較強。
5 tol: 容忍度閾值
6 max_iter: 迭代次數
7 kernel：核函式，包括:
8 linear(線性核)：u*v
9 poly(多項式核):(gamma * u * v + coef0)^degree
10 rbf(高斯核): exp(-gamma|u-v|^2)
11 sigmoid核: tanh(gamma*u*v + coef0)
12 degree: 多項式核中的維度
13 gamma： 核函式中的引數，預設為“auto”，選擇1/n_features
14 coef： 多項式核和simoid核中的常數項，僅對這兩個核函式有效
15 probability: 是否採用概率估計，預設為False
16 shrinking： 是否採用shrinking heuristic方法，預設為true
17 cache_size: 核函式的快取大小，預設為200
18 verbose: 是否允許冗餘輸出
19 class_weight: 類別權重
20 decision_function_shape: 可以取'ovo'一對一和'ovr'一對其他
21 """
然後就是資料進行訓練，並檢視訓練準確率
1 model.fit(mat(train_dataArr), mat(train_labelArr).transpose())
2 train_score = model.score(mat(train_dataArr), mat(train_labelArr).transpose())
3 print('訓練集上的準確率為%s'%train_score)
4
5 test_dataArr, test_labelArr = loadImages('E:\資料\PythonML_Code\Charpter 5\data\\testDigits'.format(dir))
6 test_score = model.score(mat(test_dataArr), mat(test_labelArr).transpose())
7 print('測試集上的準確率為%f' % test_score)
最終結果為：
至此，SVM的基本內容已基本完畢，此外還有利用SVM進行線性迴歸的演算法以及採用SVM演算法進行半監督學習的演算法，後面看到這一塊會進一步完善這一塊內容。