支持向量機—SVM原理代碼實現
本文系作者原創，轉載請注明出處:https://www.cnblogs.com/further-further-further/p/9596898.html
1. 解決什么問題？
最基本的應用是數據分類，特別是對於非線性不可分數據集。支持向量機不僅能對非線性可分數據集進行分類，對於非線性不可分數據集的也可以分類
（我認為這才是支持向量機的真正魅力所在，因為現實場景中，樣本數據往往是非線性不可分的）。
現實場景一 ：樣本數據大部分是線性可分的，但是只是在樣本中含有少量噪聲或特異點，去掉這些噪聲或特異點后線性可分 => 用支持向量機的軟間隔方法進行分類；
現實場景二：樣本數據完全線性不可分 => 引入核函數，將低維不可分的非線性數據集轉化為高維可分的數據集，用支持向量機的軟間隔方法進行分類；
一定要明白一點：對分類器准確性有影響只是樣本中的支持向量，因此，其他樣本在計算分類器的權重矩陣時可以直接過濾掉，大大節省運行時間。
目標：獲取離邊界最近樣本點到超平面最遠距離
2. SVM的思想演化？
博客 https://www.cnblogs.com/zhizhan/p/4430253.html 關於SVM的演化說得很透徹，也很形象，這里借用一下。
2.1 硬間隔支持向量機
SVM中最關鍵的思想之一就是引入和定義了“間隔”這個概念。這個概念本身很簡單，以二維空間為例，就是點到分類直線之間的距離。假設直線為y=wx+b，那么只要使所有正分類點到該直線的距離與所有負分類點到該直線的距離的總和達到最大，這條直線就是最優分類直線。這樣，原問題就轉化為一個約束優化問題，可以直接求解。這叫做硬間隔最大化，得到的SVM模型稱作硬間隔支持向量機。
2.2 軟間隔支持向量機
但是新問題出現了，在實際應用中，我們得到的數據並不總是完美的線性可分的，其中可能會有個別噪聲點，他們錯誤的被分類到了其他類中。如果將這些特異的噪點去除后，可以很容易的線性可分。但是，我們對於數據集中哪些是噪聲點卻是不知道的，如果以之前的方法進行求解，會無法進行線性分開。是不是就沒辦法了呢？假設在y=x+1直線上下分為兩類，若兩類中各有對方的幾個噪點，在人的眼中，仍然是可以將兩類分開的。這是因為在人腦中是可以容忍一定的誤差的，仍然使用y=x+1直線分類，可以在最小誤差的情況下進行最優的分類。同樣的道理，我們在SVM中引入誤差的概念，將其稱作“松弛變量”。通過加入松弛變量，在原距離函數中需要加入新的松弛變量帶來的誤差，這樣，最終的優化目標函數變成了兩個部分組成：距離函數和松弛變量誤差。這兩個部分的重要程度並不是相等的，而是需要依據具體問題而定的，因此，我們加入權重參數C，將其與目標函數中的松弛變量誤差相乘，這樣，就可以通過調整C來對二者的系數進行調和。如果我們能夠容忍噪聲，那就把C調小，讓他的權重降下來，從而變得不重要；反之，我們需要很嚴格的噪聲小的模型，則將C調大一點，權重提升上去，變得更加重要。通過對參數C的調整，可以對模型進行控制。這叫做軟間隔最大化，得到的SVM稱作軟間隔支持向量機。
2.3 非線性支持向量機
之前的硬間隔支持向量機和軟間隔支持向量機都是解決線性可分數據集或近似線性可分數據集的問題的。但是如果噪點很多，甚至會造成數據變成了線性不可分的，那該怎么辦？最常見的例子是在二維平面笛卡爾坐標系下，以原點(0,0)為圓心，以1為半徑畫圓，則圓內的點和圓外的點在二維空間中是肯定無法線性分開的。但是，學過初中幾何就知道，對於圓圈內（含圓圈）的點：x^2+y^2≤1，圓圈外的則x^2+y^2＞1。我們假設第三個維度：z=x^2+y^2，那么在第三維空間中，可以通過z是否大於1來判斷該點是否在圓內還是圓外。這樣，在二維空間中線性不可分的數據在第三維空間很容易的線性可分了。這就是非線性支持向量機。
實際中，對某個實際問題函數來尋找一個合適的空間進行映射是非常困難的，幸運的是，在計算中發現，我們需要的只是兩個向量在新的映射空間中的內積結果，而映射函數到底是怎么樣的其實並不需要知道。這一點不太好理解，有人會問，既然不知道映射函數，那怎么能知道映射后在新空間中的內積結果呢？答案其實是可以的。這就需要引入了核函數的概念。核函數是這樣的一種函數：仍然以二維空間為例，假設對於變量x和y，將其映射到新空間的映射函數為φ，則在新空間中，二者分別對應φ(x)和φ(y)，他們的內積則為<φ(x),φ(y)>。我們令函數Kernel(x,y)=<φ(x),φ(y)>=k(x,y)，可以看出，函數Kernel(x,y)是一個關於x和y的函數！而與φ無關！這是一個多么好的性質！我們再也不用管φ具體是什么映射關系了，只需要最后計算Kernel(x,y)就可以得到他們在高維空間中的內積，這樣就可以直接帶入之前的支持向量機中計算！
核函數不是很好找到，一般是由數學家反向推導出來或拼湊出來的。現在知道的線性核函數有多項式核函數、高斯核函數等。其中，高斯核函數對應的支持向量機是高斯徑向基函數（RBF），是最常用的核函數。
RBF核函數可以將維度擴展到無窮維的空間，因此，理論上講可以滿足一切映射的需求。為什么會是無窮維呢？我以前都不太明白這一點。后來老師講到，RBF對應的是泰勒級數展開，在泰勒級數中，一個函數可以分解為無窮多個項的加和，其中，每一個項可以看做是對應的一個維度，這樣，原函數就可以看做是映射到了無窮維的空間中。這樣，在實際應用中，RBF是相對最好的一個選擇。當然，如果有研究的話，還可以選用其他核函數，可能會在某些問題上表現更好。但是，RBF是在對問題不了解的情況下，對最廣泛問題效果都很不錯的核函數。因此，使用范圍也最廣。
3. SVM原理
3.1 目標函數
求使幾何間隔最大的分離超平面
以及相應的分類決策函數
也就是求出w，b
3.2 幾何間隔定義
3.3 學習的對偶算法
為了求解線性可分支持向量機的最優化問題，將它作為原始最優化問題，應用拉格朗日對偶性，通過求解對偶問題得到原始問題最優解。
3.4 最大軟間隔分類器（線性可分支持向量學習的最優化問題）
這就是線性可分支持向量機的對偶算法，這樣做的優點，一是對偶問題往往更容易求解；二是自然引入核函數，進而推廣到非線性分類問題。
的求解就轉化為 
的求解。上述公式就是對不同變量求偏導，具體推導過程見《統計學習方法》第7章 或者相關博客。
3.5 KKT條件
上圖的總結有點小問題，
分類不能合並，正確的分類描述應該是：
在超平面上或誤分的樣本點是不能正確分類的。
4. SMO算法（sequential minimal optimization 序列最小最優化）
高效實現支持向量機的算法是SMO算法，其基本思路是：如果所有變量的解都滿足此最優化的KKT條件，那么這個最優化問題的解就得到了。因為KKT條件是該優化問題的充分必要條件。否則，
選擇兩個變量，固定其他變量，針對這兩個變量構建一個二次規划問題，這個二次規划問題關於這兩個變量的解應該更接近原始二次規划問題的解，因為這會使原始二次規划問題的目標函數值變得更小。
重要的是，這時子問題可以通過解析方法求解，這樣可以大大提高整個算法的計算速度。
子問題有兩個變量，一個是違反KKT條件最嚴重的那一個，另一個由約束條件自動確定。如此，SMO算法將原問題不斷分解為子問題並對子問題求解，進而達到求解原問題的目的。
由
並且 
推出
4.1 范圍
由
假定
不固定，其他拉格朗日因子固定（也就是常量），得到
常量 ，
為變量，
為因變量，L，H的范圍就是 的 范圍
4.2 計算拉格朗日乘子 
，
第一個 公式如下，推導過程見《統計學習方法》第7章或者相關博客。
選擇一個違背KKT的數據項，根據下面結論得到優化后 ，
s表示
，從而可以得到另一個優化后的 ，同樣也需要進行約束。
4.3 計算b1，b2
在每次完成兩個變量的優化后，都要重新計算閾值b。具體推導見《統計學習方法》
5. 代碼實現
5.1 輸入數據
這里使用兩個原始數據文件 trainingData.txt，testData.txt。
trainingData.txt
1 -0.214824 0.662756 -1.000000 2 -0.061569 -0.091875 1.000000 3 0.406933 0.648055 -1.000000 4 0.223650 0.130142 1.000000 5 0.231317 0.766906 -1.000000 6 -0.748800 -0.531637 -1.000000 7 -0.557789 0.375797 -1.000000 8 0.207123 -0.019463 1.000000 9 0.286462 0.719470 -1.000000 10 0.195300 -0.179039 1.000000 11 -0.152696 -0.153030 1.000000 12 0.384471 0.653336 -1.000000 13 -0.117280 -0.153217 1.000000 14 -0.238076 0.000583 1.000000 15 -0.413576 0.145681 1.000000 16 0.490767 -0.680029 -1.000000 17 0.199894 -0.199381 1.000000 18 -0.356048 0.537960 -1.000000 19 -0.392868 -0.125261 1.000000 20 0.353588 -0.070617 1.000000 21 0.020984 0.925720 -1.000000 22 -0.475167 -0.346247 -1.000000 23 0.074952 0.042783 1.000000 24 0.394164 -0.058217 1.000000 25 0.663418 0.436525 -1.000000 26 0.402158 0.577744 -1.000000 27 -0.449349 -0.038074 1.000000 28 0.619080 -0.088188 -1.000000 29 0.268066 -0.071621 1.000000 30 -0.015165 0.359326 1.000000 31 0.539368 -0.374972 -1.000000 32 -0.319153 0.629673 -1.000000 33 0.694424 0.641180 -1.000000 34 0.079522 0.193198 1.000000 35 0.253289 -0.285861 1.000000 36 -0.035558 -0.010086 1.000000 37 -0.403483 0.474466 -1.000000 38 -0.034312 0.995685 -1.000000 39 -0.590657 0.438051 -1.000000 40 -0.098871 -0.023953 1.000000 41 -0.250001 0.141621 1.000000 42 -0.012998 0.525985 -1.000000 43 0.153738 0.491531 -1.000000 44 0.388215 -0.656567 -1.000000 45 0.049008 0.013499 1.000000 46 0.068286 0.392741 1.000000 47 0.747800 -0.066630 -1.000000 48 0.004621 -0.042932 1.000000 49 -0.701600 0.190983 -1.000000 50 0.055413 -0.024380 1.000000 51 0.035398 -0.333682 1.000000 52 0.211795 0.024689 1.000000 53 -0.045677 0.172907 1.000000 54 0.595222 0.209570 -1.000000 55 0.229465 0.250409 1.000000 56 -0.089293 0.068198 1.000000 57 0.384300 -0.176570 1.000000 58 0.834912 -0.110321 -1.000000 59 -0.307768 0.503038 -1.000000 60 -0.777063 -0.348066 -1.000000 61 0.017390 0.152441 1.000000 62 -0.293382 -0.139778 1.000000 63 -0.203272 0.286855 1.000000 64 0.957812 -0.152444 -1.000000 65 0.004609 -0.070617 1.000000 66 -0.755431 0.096711 -1.000000 67 -0.526487 0.547282 -1.000000 68 -0.246873 0.833713 -1.000000 69 0.185639 -0.066162 1.000000 70 0.851934 0.456603 -1.000000 71 -0.827912 0.117122 -1.000000 72 0.233512 -0.106274 1.000000 73 0.583671 -0.709033 -1.000000 74 -0.487023 0.625140 -1.000000 75 -0.448939 0.176725 1.000000 76 0.155907 -0.166371 1.000000 77 0.334204 0.381237 -1.000000 78 0.081536 -0.106212 1.000000 79 0.227222 0.527437 -1.000000 80 0.759290 0.330720 -1.000000 81 0.204177 -0.023516 1.000000 82 0.577939 0.403784 -1.000000 83 -0.568534 0.442948 -1.000000 84 -0.011520 0.021165 1.000000 85 0.875720 0.422476 -1.000000 86 0.297885 -0.632874 -1.000000 87 -0.015821 0.031226 1.000000 88 0.541359 -0.205969 -1.000000 89 -0.689946 -0.508674 -1.000000 90 -0.343049 0.841653 -1.000000 91 0.523902 -0.436156 -1.000000 92 0.249281 -0.711840 -1.000000 93 0.193449 0.574598 -1.000000 94 -0.257542 -0.753885 -1.000000 95 -0.021605 0.158080 1.000000 96 0.601559 -0.727041 -1.000000 97 -0.791603 0.095651 -1.000000 98 -0.908298 -0.053376 -1.000000 99 0.122020 0.850966 -1.000000 100 -0.725568 -0.292022 -1.000000
testData.txt
1 3.542485 1.977398 -1 2 3.018896 2.556416 -1 3 7.551510 -1.580030 1 4 2.114999 -0.004466 -1 5 8.127113 1.274372 1 6 7.108772 -0.986906 1 7 8.610639 2.046708 1 8 2.326297 0.265213 -1 9 3.634009 1.730537 -1 10 0.341367 -0.894998 -1 11 3.125951 0.293251 -1 12 2.123252 -0.783563 -1 13 0.887835 -2.797792 -1 14 7.139979 -2.329896 1 15 1.696414 -1.212496 -1 16 8.117032 0.623493 1 17 8.497162 -0.266649 1 18 4.658191 3.507396 -1 19 8.197181 1.545132 1 20 1.208047 0.213100 -1 21 1.928486 -0.321870 -1 22 2.175808 -0.014527 -1 23 7.886608 0.461755 1 24 3.223038 -0.552392 -1 25 3.628502 2.190585 -1 26 7.407860 -0.121961 1 27 7.286357 0.251077 1 28 2.301095 -0.533988 -1 29 -0.232542 -0.547690 -1 30 3.457096 -0.082216 -1 31 3.023938 -0.057392 -1 32 8.015003 0.885325 1 33 8.991748 0.923154 1 34 7.916831 -1.781735 1 35 7.616862 -0.217958 1 36 2.450939 0.744967 -1 37 7.270337 -2.507834 1 38 1.749721 -0.961902 -1 39 1.803111 -0.176349 -1 40 8.804461 3.044301 1 41 1.231257 -0.568573 -1 42 2.074915 1.410550 -1 43 -0.743036 -1.736103 -1 44 3.536555 3.964960 -1 45 8.410143 0.025606 1 46 7.382988 -0.478764 1 47 6.960661 -0.245353 1 48 8.234460 0.701868 1 49 8.168618 -0.903835 1 50 1.534187 -0.622492 -1 51 9.229518 2.066088 1 52 7.886242 0.191813 1 53 2.893743 -1.643468 -1 54 1.870457 -1.040420 -1 55 5.286862 -2.358286 1 56 6.080573 0.418886 1 57 2.544314 1.714165 -1 58 6.016004 -3.753712 1 59 0.926310 -0.564359 -1 60 0.870296 -0.109952 -1 61 2.369345 1.375695 -1 62 1.363782 -0.254082 -1 63 7.279460 -0.189572 1 64 1.896005 0.515080 -1 65 8.102154 -0.603875 1 66 2.529893 0.662657 -1 67 1.963874 -0.365233 -1 68 8.132048 0.785914 1 69 8.245938 0.372366 1 70 6.543888 0.433164 1 71 -0.236713 -5.766721 -1 72 8.112593 0.295839 1 73 9.803425 1.495167 1 74 1.497407 -0.552916 -1 75 1.336267 -1.632889 -1 76 9.205805 -0.586480 1 77 1.966279 -1.840439 -1 78 8.398012 1.584918 1 79 7.239953 -1.764292 1 80 7.556201 0.241185 1 81 9.015509 0.345019 1 82 8.266085 -0.230977 1 83 8.545620 2.788799 1 84 9.295969 1.346332 1 85 2.404234 0.570278 -1 86 2.037772 0.021919 -1 87 1.727631 -0.453143 -1 88 1.979395 -0.050773 -1 89 8.092288 -1.372433 1 90 1.667645 0.239204 -1 91 9.854303 1.365116 1 92 7.921057 -1.327587 1 93 8.500757 1.492372 1 94 1.339746 -0.291183 -1 95 3.107511 0.758367 -1 96 2.609525 0.902979 -1 97 3.263585 1.367898 -1 98 2.912122 -0.202359 -1 99 1.731786 0.589096 -1 100 2.387003 1.573131 -1
5.2 SMO算法實現
定義數據結構體optStruct，用於緩存，提高運行速度。SMO算法具體實現如下（mySVMMLiA.py）
每個方法的作用，以及每行代碼的作用，我都做了詳細的注解，希望對大家的理解有幫助。
1 # -*- coding: utf-8 -*- 2 """ 3 Created on Tue Sep 4 16:58:16 2018 4 支持向量機代碼實現 5 SMO(Sequential Minimal Optimization)最小序列優化 6 @author: weixw 7 """ 8 import numpy as np 9 #核轉換函數（一個特征空間映射到另一個特征空間，低維空間映射到高維空間） 10 #高維空間解決線性問題，低維空間解決非線性問題 11 #線性內核 = 原始數據矩陣（100*2）與原始數據第一行矩陣轉秩乘積（2*1） =>（100*1） 12 #非線性內核公式：k(x,y) = exp(-||x - y||**2/2*(e**2)) 13 #1.原始數據每一行與原始數據第一行作差， 14 #2.平方 15 def kernelTrans(dataMat, rowDataMat, kTup): 16 m,n=np.shape(dataMat) 17 #初始化核矩陣 m*1 18 K = np.mat(np.zeros((m,1))) 19 if kTup[0] == 'lin': #線性核 20 K = dataMat*rowDataMat.T 21 elif kTup[0] == 'rbf':#非線性核 22 for j in range(m): 23 #xi - xj 24 deltaRow = dataMat[j,:] - rowDataMat 25 K[j] = deltaRow*deltaRow.T 26 #1*m m*1 => 1*1 27 K = np.exp(K/(-2*kTup[1]**2)) 28 else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized') 29 return K 30 31 #定義數據結構體，用於緩存，提高運行速度 32 class optStruct: 33 def __init__(self, dataSet, labelSet, C, toler, kTup): 34 self.dataMat = np.mat(dataSet) #原始數據，轉換成m*n矩陣 35 self.labelMat = np.mat(labelSet).T #標簽數據 m*1矩陣 36 self.C = C #懲罰參數，C越大，容忍噪聲度小，需要優化；反之，容忍噪聲度高，不需要優化； 37 #所有的拉格朗日乘子都被限制在了以C為邊長的矩形里 38 self.toler = toler #容忍度 39 self.m = np.shape(self.dataMat)[0] #原始數據行長度 40 self.alphas = np.mat(np.zeros((self.m,1))) # alpha系數，m*1矩陣 41 self.b = 0 #偏置 42 self.eCache = np.mat(np.zeros((self.m,2))) # 保存原始數據每行的預測值 43 self.K = np.mat(np.zeros((self.m,self.m))) # 核轉換矩陣 m*m 44 for i in range(self.m): 45 self.K[:,i] = kernelTrans(self.dataMat, self.dataMat[i,:], kTup) 46 47 #計算原始數據第k項對應的預測誤差 1*m m*1 =>1*1 48 #oS：結構數據 49 #k： 原始數據行索引 50 def calEk(oS, k): 51 #f(x) = w*x + b 52 fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b) 53 Ek = fXk - float(oS.labelMat[k]) 54 return Ek 55 56 #在alpha有改變都要更新緩存 57 def updateEk(oS, k): 58 Ek = calEk(oS, k) 59 oS.eCache[k] = [1, Ek] 60 61 62 #第一次通過selectJrand()隨機選取j,之后選取與i對應預測誤差最大的j（步長最大） 63 def selectJ(i, oS, Ei): 64 #初始化 65 maxK = -1 #誤差最大時對應索引 66 maxDeltaE = 0 #最大誤差 67 Ej = 0 # j索引對應預測誤差 68 #保存每一行的預測誤差值 1相對於初始化為0的更改 69 oS.eCache[i] = [1,Ei] 70 #獲取數據緩存結構中非0的索引列表(先將矩陣第0列轉化為數組) 71 validEcacheList = np.nonzero(oS.eCache[:,0].A)[0] 72 #遍歷索引列表，尋找最大誤差對應索引 73 if len(validEcacheList) > 1: 74 for k in validEcacheList: 75 if k == i: 76 continue 77 Ek = calEk(oS, k) 78 deltaE = abs(Ei - Ek) 79 if(deltaE > maxDeltaE): 80 maxK = k 81 maxDeltaE = deltaE 82 Ej = Ek 83 return maxK, Ej 84 else: 85 #隨機選取一個不等於i的j 86 j = selectJrand(i, oS.m) 87 Ej = calEk(oS, j) 88 return j,Ej 89 90 #隨機選取一個不等於i的索引 91 def selectJrand(i, m): 92 j = i 93 while (j == i): 94 j = int(np.random.uniform(0, m)) 95 return j 96 97 #alpha范圍剪輯 98 def clipAlpha(aj, L, H): 99 if aj > H: 100 aj = H 101 if aj < L: 102 aj = L 103 return aj 104 105 #從文件獲取特征數據，標簽數據 106 def loadDataSet(fileName): 107 dataSet = []; labelSet = [] 108 fr = open(fileName) 109 for line in fr.readlines(): 110 #分割 111 lineArr = line.strip().split('\t') 112 dataSet.append([float(lineArr[0]), float(lineArr[1])]) 113 labelSet.append(float(lineArr[2])) 114 return dataSet, labelSet 115 116 #計算 w 權重系數 117 def calWs(alphas, dataSet, labelSet): 118 dataMat = np.mat(dataSet) 119 #1*100 => 100*1 120 labelMat = np.mat(labelSet).T 121 m, n = np.shape(dataMat) 122 w = np.zeros((n, 1)) 123 for i in range(m): 124 w += np.multiply(alphas[i]*labelMat[i], dataMat[i,:].T) 125 return w 126 #計算原始數據每一行alpha,b，保存到數據結構中，有變化及時更新 127 def innerL(i, oS): 128 #計算預測誤差 129 Ei = calEk(oS, i) 130 #選擇第一個alpha，違背KKT條件2 131 #正間隔，負間隔 132 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)): 133 #第一次隨機選取不等於i的數據項，其后根據誤差最大選取數據項 134 j, Ej = selectJ(i, oS, Ei) 135 #初始化，開辟新的內存 136 alphaIold = oS.alphas[i].copy() 137 alphaJold = oS.alphas[j].copy() 138 #通過 a1y1 + a2y2 = 常量 139 # 0 <= a1,a2 <= C 求出L,H 140 if oS.labelMat[i] != oS.labelMat[j]: 141 L = max(0, oS.alphas[j] - oS.alphas[i]) 142 H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i]) 143 else: 144 L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C) 145 H = min(oS.C, oS.alphas[j] + oS.alphas[i]) 146 if L == H : 147 print ("L == H") 148 return 0 149 #內核分母 K11 + K22 - 2K12 150 eta = oS.K[i, i] + oS.K[j, j] - 2.0*oS.K[i, j] 151 if eta <= 0: 152 print ("eta <= 0") 153 return 0 154 #計算第一個alpha j 155 oS.alphas[j] += oS.labelMat[j]*(Ei - Ej)/eta 156 #修正alpha j的范圍 157 oS.alphas[j] = clipAlpha(oS.alphas[j], L, H) 158 #alpha有改變，就需要更新緩存數據 159 updateEk(oS, j) 160 #如果優化后的alpha 與之前的alpha變化很小，則舍棄，並重新選擇數據項的alpha 161 if (abs(oS.alphas[j] - alphaJold) < 0.00001): 162 print ("j not moving enough, abandon it.") 163 return 0 164 #計算alpha對的另一個alpha i 165 # ai_new*yi + aj_new*yj = 常量 166 # ai_old*yi + ai_old*yj = 常量 167 # 作差=> ai = ai_old + yi*yj*(aj_old - aj_new) 168 oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j]) 169 #alpha有改變，就需要更新緩存數據 170 updateEk(oS, i) 171 #計算b1,b2 172 # y(x) = w*x + b => b = y(x) - w*x 173 # w = aiyixi(i= 1->N求和) 174 #b1_new = y1_new - (a1_new*y1*k11 + a2_new*y2*k21 + ai*yi*ki1(i = 3 ->N求和 常量)) 175 #b1_old = y1_old - (a1_old*y1*k11 + a2_old*y2*k21 + ai*yi*ki1(i = 3 ->N求和 常量)) 176 #作差=> b1_new = b1_old + (y1_new - y1_old) - y1*k11*(a1_new - a1_old) - y2*k21*(a2_new - a2_old) 177 # => b1_new = b1_old + Ei - yi*(ai_new - ai_old)*kii - yj*(aj_new - aj_old)*kij 178 #同樣可推得 b2_new = b2_old + Ej - yi*(ai_new - ai_old)*kij - yj*(aj_new - aj_old)*kjj 179 bi = oS.b - Ei - oS.labelMat[i]*(oS.alphas[i] - alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j] - alphaJold)*oS.K[i,j] 180 bj = oS.b - Ej - oS.labelMat[i]*(oS.alphas[i] - alphaIold)*oS.K[i,j] - oS.labelMat[j]*(oS.alphas[j] - alphaJold)*oS.K[j,j] 181 #首選alpha i，相對alpha j 更准確 182 if (0 < oS.alphas[i]) and (oS.alphas[i] < oS.C): 183 oS.b = bi 184 elif (0 < oS.alphas[j]) and (oS.alphas[j] < oS.C): 185 oS.b = bj 186 else: 187 oS.b = (bi + bj)/2.0 188 return 1 189 else: 190 return 0 191 192 #完整SMO核心算法，包含線性核核非線性核，返回alpha,b 193 #dataSet 原始特征數據 194 #labelSet 標簽數據 195 #C 凸二次規划參數 196 #toler 容忍度 197 #maxInter 循環次數 198 #kTup 指定核方式 199 #程序邏輯： 200 #第一次全部遍歷，遍歷后根據alpha對是否有修改判斷， 201 #如果alpha對沒有修改，外循環終止；如果alpha對有修改，則繼續遍歷屬於支持向量的數據。 202 #直至外循環次數達到maxIter 203 #相比簡單SMO算法，運行速度更快，原因是： 204 #1.不是每一次都全量遍歷原始數據，第一次遍歷原始數據， 205 #如果alpha有優化，就遍歷支持向量數據，直至alpha沒有優化，然后再轉全量遍歷，這是如果alpha沒有優化，循環結束； 206 #2.外循環不需要達到maxInter次數就終止； 207 def smoP(dataSet, labelSet, C, toler, maxInter, kTup = ('lin', 0)): 208 #初始化結構體類，獲取實例 209 oS = optStruct(dataSet, labelSet, C, toler, kTup) 210 iter = 0 211 #全量遍歷標志 212 entireSet = True 213 #alpha對是否優化標志 214 alphaPairsChanged = 0 215 #外循環 終止條件：1.達到最大次數 或者 2.alpha對沒有優化 216 while (iter < maxInter) and ((alphaPairsChanged > 0) or (entireSet)): 217 alphaPairsChanged = 0 218 #全量遍歷 ，遍歷每一行數據 alpha對有修改，alphaPairsChanged累加 219 if entireSet: 220 for i in range(oS.m): 221 alphaPairsChanged += innerL(i, oS) 222 print ("fullSet, iter: %d i:%d, pairs changed %d" %(iter, i, alphaPairsChanged)) 223 iter += 1 224 else: 225 #獲取(0，C)范圍內數據索引列表，也就是只遍歷屬於支持向量的數據 226 nonBounds = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] 227 for i in nonBounds: 228 alphaPairsChanged += innerL(i, oS) 229 print ("non-bound, iter: %d i:%d, pairs changed %d" %(iter, i, alphaPairsChanged)) 230 iter += 1 231 #全量遍歷->支持向量遍歷 232 if entireSet: 233 entireSet = False 234 #支持向量遍歷->全量遍歷 235 elif alphaPairsChanged == 0: 236 entireSet = True 237 print ("iteation number: %d"% iter) 238 print ("entireSet :%s"% entireSet) 239 print ("alphaPairsChanged :%d"% alphaPairsChanged) 240 return oS.b,oS.alphas 241 242 #繪制支持向量 243 def drawDataMap(dataArr,labelArr,b,alphas): 244 import matplotlib.pyplot as plt 245 #alphas.A>0 獲取大於0的索引列表，只有>0的alpha才對分類起作用 246 svInd=np.nonzero(alphas.A>0)[0] 247 #分類數據點 248 classified_pts = {'+1':[],'-1':[]} 249 for point,label in zip(dataArr,labelArr): 250 if label == 1.0: 251 classified_pts['+1'].append(point) 252 else: 253 classified_pts['-1'].append(point) 254 fig = plt.figure() 255 ax = fig.add_subplot(111) 256 #繪制數據點 257 for label,pts in classified_pts.items(): 258 pts = np.array(pts) 259 ax.scatter(pts[:, 0], pts[:, 1], label = label) 260 #繪制分割線 261 w = calWs(alphas, dataArr, labelArr) 262 #函數形式：max( x ,key=lambda a : b ) # x可以是任何數值，可以有多個x值 263 #先把x值帶入lambda函數轉換成b值，然后再將b值進行比較 264 x1, _=max(dataArr, key=lambda x:x[0]) 265 x2, _=min(dataArr, key=lambda x:x[0]) 266 a1, a2 = w 267 y1, y2 = (-b - a1*x1)/a2, (-b - a1*x2)/a2 268 #矩陣轉化為數組.A 269 ax.plot([x1, x2],[y1.A[0][0], y2.A[0][0]]) 270 271 #繪制支持向量 272 for i in svInd: 273 x, y= dataArr[i] 274 ax.scatter([x], [y], s=150, c ='none', alpha=0.7, linewidth=1.5, edgecolor = '#AB3319') 275 plt.show() 276 277 #alpha>0對應的數據才是支持向量，過濾不是支持向量的數據 278 sVs= np.mat(dataArr)[svInd] #get matrix of only support vectors 279 print ("there are %d Support Vectors.\n" % np.shape(sVs)[0]) 280 281 #訓練結果 282 def getTrainingDataResult(dataSet, labelSet, b, alphas, k1=1.3): 283 datMat = np.mat(dataSet) 284 #100*1 285 labelMat = np.mat(labelSet).T 286 #alphas.A>0 獲取大於0的索引列表，只有>0的alpha才對分類起作用 287 svInd=np.nonzero(alphas.A>0)[0] 288 sVs=datMat[svInd] 289 labelSV = labelMat[svInd]; 290 m,n = np.shape(datMat) 291 errorCount = 0 292 for i in range(m): 293 kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1)) 294 # y(x) = w*x + b => b = y(x) - w*x 295 # w = aiyixi(i= 1->N求和) 296 predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b 297 if np.sign(predict)!=np.sign(labelSet[i]): errorCount += 1 298 print ("the training error rate is: %f" % (float(errorCount)/m)) 299 300 def getTestDataResult(dataSet, labelSet, b, alphas, k1=1.3): 301 datMat = np.mat(dataSet) 302 #100*1 303 labelMat = np.mat(labelSet).T 304 #alphas.A>0 獲取大於0的索引列表，只有>0的alpha才對分類起作用 305 svInd=np.nonzero(alphas.A>0)[0] 306 sVs=datMat[svInd] 307 labelSV = labelMat[svInd]; 308 m,n = np.shape(datMat) 309 errorCount = 0 310 for i in range(m): 311 kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1)) 312 # y(x) = w*x + b => b = y(x) - w*x 313 # w = aiyixi(i= 1->N求和) 314 predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b 315 if np.sign(predict)!=np.sign(labelSet[i]): errorCount += 1 316 print ("the test error rate is: %f" % (float(errorCount)/m)) 317 318
5.3 測試代碼（testMySVMMLiA.py）
1 # -*- coding: utf-8 -*- 2 """ 3 Created on Wed Sep 5 15:22:26 2018 4 5 @author: weixw 6 """ 7 8 9 import mySVMMLiA as sm 10 11 #通過訓練數據計算 b, alphas 12 dataArr,labelArr = sm.loadDataSet('trainingData.txt') 13 b, alphas = sm.smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', 0.10)) 14 sm.drawDataMap(dataArr,labelArr,b,alphas) 15 sm.getTrainingDataResult(dataArr, labelArr, b, alphas, 0.10) 16 dataArr1,labelArr1 = sm.loadDataSet('testData.txt') 17 #測試結果 18 sm.getTestDataResult(dataArr1, labelArr1, b, alphas, 0.10)
通過訓練數據計算出 b， 權重矩陣，從而分類超平面和決策分類函數就明確了，然后測試數據以決策分類函數進行預測。
這里采用高斯核RBF。
5.4 運行結果
“j not moving enough, abandon it” 表示數據項對應的 
和 
非常接近，不需要優化；
“fullSet” 表示全量數據遍歷；
“non-bound” 表示非邊界遍歷，也就是只遍歷屬於支持向量的數據項。
另外我將支持向量的數據項繪制出來了，這樣更直觀。
可以看出，有77個支持向量，訓練差錯率是0，測試差錯率6%。
6. 參考文獻
《統計學習方法》第 七 章
《機器學習實戰》第 六 章
《從零推導支持向量機》
博客：https://www.cnblogs.com/bentuwuying/p/6444516.html
博客：http://www.cnblogs.com/jerrylead/archive/2011/03/18/1988419.html
博客：https://www.cnblogs.com/zhizhan/p/4430253.html
博客：https://www.cnblogs.com/xxrxxr/p/7538430.html
不要讓懶惰占據你的大腦，不要讓妥協拖垮了你的人生。青春就是一張票，能不能趕上時代的快車，你的步伐就掌握在你的腳下。