支援向量機|SMO演算法實現
01 起
我們一起回憶一下，支援向量機可以處理三種型別的資料：
線性可分支援向量機——求解策略，硬間隔最大化
線性支援向量機——求解策略，軟間隔最大化
非線性支援向量機——求解策略，核技巧+軟間隔最大化
我們提出一個問題
當資料量很大時，以上提出的演算法求解複雜度呈指數上升，演算法會變得十分低效，該怎麼辦呢？
我們給出瞭解法思路
在對偶問題中，每次只求解最佳化兩個alpha的值，如此遍歷求解的方法與一次求解所有alpha的方法得到的結果是完全一致的。
其實這就是SMO(序列最小最佳化演算法)的原理，今天我們將基於SMO原理，自己寫一個程式碼來實現SMO演算法。
02 SMO演算法原理
在統計學習方法中，給出SMO演算法思路如下
我們可以將SMO演算法過程歸納如下：
初始化所有變數的解ai=0
選取最佳化變數a1,a2，解析求解這兩個變數的二次規劃問題，得到最優解a1',a2'
驗證，所有變數的解ai是否滿足KKT條件，若不滿足，繼續step2，否則step4
得到a=a1',a2'......
SMO演算法求解的是對偶問題的解alpha，在得到alpha值後，可以透過下面的公式求出引數向量w和常數項b，從而得到分離超平面
SMO演算法有兩個關鍵點：
選擇兩個alpha進行最佳化
求解所選alpha的解析解
其中，
選擇第一個alpha(a1)我們稱之為外迴圈，選擇原則是違反KKT條件最嚴重的點(SMO演算法在，選擇足夠違反(>toler)KKT條件的樣本點)，選擇第二個alpha(a2)為內迴圈，選擇原則是與a1變化相反的點(E1-E2最大的樣本點)
-
求解a1,a2解析解的問題，就是兩個變數的二次規劃問題，可以透過一下公式求得
兩個變數的二次規劃問題
a1,a2解析解
下面我們基於python3實現演算法！
03 輔助函式
輔助函式包括，
資料載入函式
根據alphai的i隨機選擇j的函式
根據L H邊界值剪下alpha的函式
根據alpha計算w的函式
繪圖函式
def loadDataSet(filename):
#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,labelMatdef randPickj(i,m):
#i是alphai的i值，整數; m是alpha個數; j不能等於i
j=i while j==i:
j=int(np.random.uniform(0,m)) return jdef clipAlpha(aj,H,L):
if aj>H:
aj=H if aj根據alpha計算w的函式
def weight(data,label,alphas):
dataMatrix=np.mat(data);labelMatrix=np.mat(label).transpose() #這裡labelMatrix形狀為m*1
m,n=dataMatrix.shape
w=np.mat(np.zeros((1,n))) #初始化w，為1行n列的全零矩陣，n為data維度數
"""w=求和(ai*yi*xi),求和物件是支援向量，即，ai>0的樣本點，xi，yi為支援向量對應的label和data"""
for i in range(m): if alphas[i]>0:
w+=labelMatrix[i]*alphas[i]*dataMatrix[i,:] return w.tolist()
繪圖函式
"""
繪製樣本資料以及決策邊界
思路：
1. 將樣本資料根據樣本類別標籤labelMat分別放入不同的座標集中
2. 根據座標集合，分別繪製兩個類別樣本的散點圖
3. 決策邊界即x2=f(x1),由w1*x1+w2*x2+b=0得到x2(即y=(-b-w1x1)/w2)
"""def plotBestFit(weights,b,filename):
dataMat,labelMat=loadDataSet(filename) #載入樣本特徵、樣本類別
dataArr=np.array(dataMat)
n=dataArr.shape[0] #n個樣本
xcord1=[];ycord1=[]
xcord2=[];ycord2=[] #兩個類別的樣本的xy座標值，x對應x1,y對應x2
#將樣本資料根據樣本類別標籤labelMat分別放入不同的座標集中
for i in range(n): if int(labelMat[i])==1: #第i個樣本是1類
xcord1.append(dataArr[i,0]) #第i個樣本的x1值
ycord1.append(dataArr[i,1]) #第i個樣本的x2值
else:
xcord2.append(dataArr[i,0]) #第i個樣本的x1值
ycord2.append(dataArr[i,1]) #第i個樣本的x2值
#繪製兩類樣本的散點圖
fig=plt.figure(figsize=(12,8))
plt.scatter(xcord1,ycord1,c="red",s=50,label="label=1")
plt.scatter(xcord2,ycord2,c="blue",s=50,label="label=-1") #繼續在原圖上作圖
#繪製決策邊界
x=np.arange(-3.0,5.0,0.1)
y=(-b-weights[0][0]*x)/weights[0][1] #由w1*x1+w2*x2+b=0得到x2(即y)=(-b-w1x1)/w2
x.shape=(len(x),1);y.shape=(len(x),1)
plt.plot(x,y,color="darkorange",linewidth=3.0,label="Boarder") #繼續在ax圖上作圖
plt.xlabel("X1",fontsize=16)
plt.ylabel("X2",fontsize=16)
plt.title("SMO BestFit",fontsize=20,fontweight="bold")
plt.legend() #新增圖示註解
plt.show()
04 簡化版SMO演算法
說明：
簡化部分在於省去了選擇alphai的外迴圈過程，省去了更新Ei的過程
寫這個演算法的目的在於深入理解SMO演算法過程
4.1 程式碼
def SMOsimple(data,label,C,toler,maxIter):
"""
data:樣本各屬性值
label：各樣本對應標籤
C：軟間隔最大化的鬆弛係數對應的懲罰因子，也是約束條件中alpha的上界(對於線性可分資料集，C作用不大；對於線性不可分資料集，結果對C敏感)
toler：容錯率，偏離KKT條件的容錯率
maxIter：外層迴圈迭代次數
"""
#初始化alpha=0,b=0,alpha個數為樣本數，一個樣本對應一個alpha
dataMatrix=np.mat(data);labelMatrix=np.mat(label).transpose() #這裡labelMatrix形狀為m*1
b=0;m,n=dataMatrix.shape
alphas=np.mat(np.zeros((m,1)))
iters=0
while iterstoler and alphas[i]>0):
j=randPickj(i,m)
gxj=float(np.multiply(alphas,labelMatrix).transpose()*(dataMatrix*dataMatrix[j,:].transpose()))+b
Ej=gxj-labelMatrix[j]
#儲存alpha初始值，用於後續計算
alphaIold=alphas[i].copy()
alphaJold=alphas[j].copy()
#計算剪下邊界(很簡單的幾何計算，見統計學習方法)
if labelMatrix[i]!=labelMatrix[j]:
L=max(0,alphas[j]-alphas[i]) #這裡alpha[i]仍然等於alphaIold
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 #第一個跳出條件(跳出本次內迴圈，遍歷下一個alpha進行更新)
#計算eta
eta=dataMatrix[i,:]*dataMatrix[i,:].transpose()+dataMatrix[j,:]*dataMatrix[j,:].transpose() -2.0*dataMatrix[i,:]*dataMatrix[j,:].transpose() if eta==0: print ("eta=0") continue #第二個跳出條件(因為eta=0不好處理，且出現情況較少，因此這裡我們不處理，直接跳出)
#根據統計學習方法中的結果公式得到alphaj的解析解
alphas[j]=alphas[j]+labelMatrix[j]*(Ei-Ej)/eta
alphas[j]=clipAlpha(alphas[j],H,L)
#檢驗alphaj與alphaJold是否有足夠大的改變，若改變不夠大，說明與alpha舊值沒有什麼差異，跳出本次內迴圈
if alphas[j]-alphaJold0 and alphas[i]0 and alphas[j]4.2 測試
我們用一個線性可分的資料集進行測試，資料集共100個樣本點，結果如下(線性不可分資料集也進行了測試，效果良好，篇幅限制，這裡只展示線性可分資料集的測試結果)
dataMat1,labelMat1=loadDataSet("SMO_data2.txt") #SMO_data2.txt是線性可分資料集start=time.time()
b1,alphas1=SMOsimple(dataMat1,labelMat1,0.6,0.001,100)print ("n","time used:.{0}s".format(time.time()-start))
w1=weight(dataMat1,labelMat1,alphas1)
plotBestFit(w1,b1,"SMO_data2.txt")
簡化版SMO演算法用時8.67s完成最佳化
得到的分離超平面如下：
05 完整版SMO演算法
在簡化版的SMO演算法測試中，雖然我們只使用了100個樣本點進行測試，但演算法也用了8.67s才完成最佳化，如果是體量更加巨大的資料集，演算法效率將十分低下！
怎麼辦呢？完整版演算法來解決！相比於簡化版，完整版最佳化了：
外迴圈中用啟發方法選擇ai：在全集遍歷迴圈和邊界遍歷迴圈中切換搜尋違反KKT條件超過容錯率toler的點
先遍歷全集，若全集遍歷中未修改任何alpha，說明alpha均滿足KKT條件，迴圈停止；
若全集遍歷中修改了alpha值，則在下一輪迴圈進入邊界遍歷(邊界指0
在邊界遍歷中，修改邊界點的alpha，直到無alpha可修改(此時alphaPairChanged=0)，再次進入全集遍歷；
如此迴圈往復，直到全集遍歷中未修改任何alpha，迴圈停止
內迴圈中透過選擇相較ai具有最大步長(即Ei-Ej)的aj
每次修改ai aj後，緊跟著修改Ei Ej
5.1 程式碼
完整版輔助函式
包括，計算Ei的函式、選擇相較ai具有最大步長(即Ei-Ej)的aj的函式、更新Ei矩陣的函式
class optStruct:
def __init__(self,data,label,C,toler): #全域性變數
self.X=data self.labelMatrix=label self.C=C self.toler=toler self.m=data.shape[0] #m為樣本數
#初始化alpha矩陣、b、Es矩陣
self.alphas=np.mat(np.zeros((self.m,1))) self.Es=np.mat(np.zeros((self.m,2))) #快取誤差，兩列，第一列表示當前Ei是否有效，第二列表示當前的Ei值
self.b=0
def calcEk(oS,k):
gxk=float(np.multiply(oS.alphas,oS.labelMatrix).transpose()*(oS.X*oS.X[k,:].transpose()))+oS.b
Ek=gxk-float(oS.labelMatrix[k]) return Ek
#選擇相較ai具有最大步長(即Ei-Ej)的aj的函式def selectJ(oS,i,Ei):
maxK=-1;maxDeltaE=0;Ej=0 #DeltaE表示Ei-Ej,k表示DeltaE最大的樣本點索引值，最終會將Ek賦值給Ej
oS.Es[i]=[1,Ei] #使Es矩陣第i位有效
validEsList=np.nonzero(oS.Es[:,0].A)[0] #將Es矩陣中有效的Ei對應的索引值選出來，作為挑選j的池子
if len(validEsList)>1: for k in validEsList:
if k==i:
continue
Ek=calcEk(oS,k)
deltaE=abs(Ei-Ek) if deltaE>maxDeltaE:
maxDeltaE=deltaE;maxK=k;Ej=Ek return maxK,Ej else: #若validEsList只有一個Ei有效(初次迴圈)，則隨機選取一個j
j=randPickj(i,oS.m)
Ej=calcEk(oS,j) return j,Ej
def updateEk(oS,k):
Ek=calcEk(oS,k)
oS.Es[k]=[1,Ek]
內迴圈
外迴圈是對ai的迴圈，內迴圈是在ai選定的基礎下對aj的迴圈
程式碼和邏輯與SMO簡化版相似(因為簡化版SMO未對外迴圈做任何最佳化)
def innerL(i,oS):
Ei=calcEk(oS,i)
#判斷Ei是否是違反KKT條件超過toler的點，若是再繼續挑選j
if (oS.labelMatrix[i]*EioS.toler and oS.alphas[i]>0):
j,Ej=selectJ(oS,i,Ei)
alphaIold=oS.alphas[i].copy();alphaJold=oS.alphas[j].copy()
#計算L，H
if oS.labelMatrix[i]!=oS.labelMatrix[j]:
L=max(0,oS.alphas[j]-oS.alphas[i]) #這裡alpha[i]仍然等於alphaIold
H=min(oS.C,oS.C+oS.alphas[j]-oS.alphas[i])
else:
L=max(0,oS.alphas[j]+oS.alphas[i]-oS.C)
H=min(oS.C,oS.alphas[j]+oS.alphas[i]) if L==H: print ("L==H") return 0 #第一個跳出條件(跳出本次內迴圈，遍歷下一個alpha進行更新)
#計算eta
eta=oS.X[i,:]*oS.X[i,:].transpose()+oS.X[j,:]*oS.X[j,:].transpose()-2.0*oS.X[i,:]*oS.X[j,:].transpose() if eta==0: print ("eta=0") return 0 #第二個跳出條件(因為eta=0不好處理，且出現情況較少，因此這裡我們不處理，直接跳出)
#根據統計學習方法中的結果公式得到alphaj的解析解，並更新Ej值
oS.alphas[j]=oS.alphas[j]+oS.labelMatrix[j]*(Ei-Ej)/eta
oS.alphas[j]=clipAlpha(oS.alphas[j],H,L)
updateEk(oS,j) #更新Ej值
#檢驗alphaj與alphaJold是否有足夠大的改變，若改變不夠大，說明與alpha舊值沒有什麼差異，跳出本次內迴圈
if abs(oS.alphas[j]-alphaJold)0 and oS.alphas[i]0 and oS.alphas[j]外迴圈
簡化版未對外迴圈做任何操作，這是完整版SMO與簡化版的差異之一
外迴圈選擇ai的i的邏輯：遍歷全集，若全集無修改alpha對，則說明alpha已符合要求，迴圈停止；否則下一輪進入邊界遍歷修改alpha對，直到邊界遍歷中再無alpha對可修改；則下一輪進入全集遍歷，如此迴圈往復尋找ai，直到全集無修改alpha對時停止
def SMOpro(data,label,C,toler,maxIter,kTup=("lin",0)):
oS=optStruct(np.mat(data),np.mat(label).transpose(),C,toler)
iter=0;entireSet=True;alphaPairsChanged=0
#當迭代次數達到上限(這裡的迭代次數只要完成一次迴圈遍歷就+1，不論該次迴圈遍歷是否修改了alpha對)，或全集再無可修改的alpha對時，迴圈停止，計算完成
while (iter0):
alphaPairsChanged=0
if entireSet: #全集遍歷
for i in range(oS.m):
alphaPairsChanged+=innerL(i,oS) print ("fullset, iter:%d i:%d, pairsChanged: %d" %(iter,i,alphaPairsChanged))
iter+=1 #這裡的迭代次數只要完成一次迴圈遍歷就+1，不論該次迴圈遍歷是否修改了alpha對
else: #邊界遍歷
boundIs=np.nonzero((oS.alphas.A>0)*(oS.alphas.A5.2 測試
使用與簡化版SMO同樣的測試資料集進行測試
data3,label3=loadDataSet("SMO_data2.txt")
start=time.time()
b3,alphas3=SMOpro(data3,label3,0.6,0.001,60)print ("n","time used:.{0}s".format(time.time()-start))
w3=weight(data3,label3,alphas3)
plotBestFit(w3,b3,"SMO_data2.txt")
時間效率相比簡化版提升了40+倍！
得到的分離超平面如下，對比簡化版的分離超平面，分離得好像更好呢：
06 總結
本文介紹了支援向量機SMO的演算法原理，然後基於演算法原理，利用python3實現了SMO簡化版/完整版演算法，旨在深刻理解 SVM原理，SMO過程，希望對你有幫助。
下一期準備嘗試用SVM演算法分類非線性資料集，敬請期待~
作者：鄧莎
連結：
來自 “ ITPUB部落格 ” ，連結：http://blog.itpub.net/4422/viewspace-2815013/，如需轉載，請註明出處，否則將追究法律責任。
相關文章
- 支援向量機（五）SMO演算法演算法
- 支援向量機原理(四)SMO演算法原理演算法
- 《Machine Learning in Action》—— 剖析支援向量機，優化SMOMac優化
- 《機器學習_07_01_svm_硬間隔支援向量機與SMO》機器學習
- SVM 支援向量機演算法-實戰篇演算法
- 支援向量機（SVM）和python實現（二）Python
- 支援向量機python實現（簡易版）Python
- OpenCV筆記（3）實現支援向量機（SVM）OpenCV筆記
- 支援向量機原理(一) 線性支援向量機
- 分類演算法-支援向量機 SVM演算法
- 支援向量機
- SVM 支援向量機演算法-原理篇演算法
- 演算法金 | 再見，支援向量機 SVM！演算法
- 理解支援向量機
- 機器學習實戰6（SMO演算法）機器學習演算法
- 支援向量機（SVM）從原理到python程式碼實現Python
- 機器學習：支援向量機(SVM)機器學習
- 機器學習——支援向量機（SVM）機器學習
- 機器學習演算法筆記之5：支援向量機SVM機器學習演算法筆記
- 機器學習演算法（6）用Python實現用核支援向量機求解非線性問題機器學習演算法Python
- 機器學習之支援向量機原理和sklearn實踐機器學習
- 一文讀懂支援向量機SVM（附實現程式碼、公式）公式
- scikit-learn 支援向量機演算法庫使用小結演算法
- [譯] 支援向量機（SVM）教程
- 機器學習——支援向量機SVM（一）機器學習
- 【機器學習】支援向量機分類機器學習
- svm的smo演算法演算法
- 支援向量機原理(三)線性不可分支援向量機與核函式函式
- 支援向量機原理(二) 線性支援向量機的軟間隔最大化模型模型
- 【Python機器學習實戰】感知機和支援向量機學習筆記（三）之SVM的實現Python機器學習筆記
- 透過現象看本質，圖解支援向量機圖解
- 支援向量機（三）核函式函式
- 【機器學習】支援向量機（個人筆記）機器學習筆記
- 支援向量機原理(五)線性支援迴歸
- 監督學習之支援向量機
- 機器學習之支援向量機的超平面機器學習
- 人工智慧-機器學習-支援向量機SVM人工智慧機器學習
- 支援向量機(SVM)的約束和無約束優化、理論和實現優化