【Python機器學習實戰】感知機和支援向量機學習筆記（三）之SVM的實現

Uniqe發表於2021-08-11
Python機器學習筆記

前面已經對感知機和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演算法進行半監督學習的演算法，後面看到這一塊會進一步完善這一塊內容。 

 

 

 

 

相關文章