《機器學習實戰》學習總結（七）——支持向量機SVM（2）

摘要

1. 線性支持向量機

2. 非線性支持向量機和核函數

3. SMO演算法

4. 支持向量回歸

5. 代碼實現與注釋

上一節我們一起了解了支持向量機的基本概念、拉格朗日乘子法和對偶問題。這一節我們一起深入到支持向量機的具體問題及其解法。

線性支持向量機

1. 線性可分支持向量機

線性可分支持向量機是SVM中最簡單最基本的形式，如下圖（圖片來自網路，侵刪）：

即通過一個超平面就可以完全將兩個類別分開。

這種線性可分支持向量機的基本型：

首先我們利用拉格朗日乘子法引入拉格朗日函數：

上一節講過，原始問題的對偶問題就是極大極小問題：

我們先求極小，再求極大：

（1）求min L(w,b,α)

將(1)式帶入拉格朗日函數，並考慮(2)式的約束，得到：

即：

(2) 求L(w,b,α)對α的極大，即得：

即可得到下列優化問題：

將目標函數由求最大值轉換為求最小值，得到：

原始問題就轉化為了上述對α的優化問題。

這時，我們假設求得了α的最優解：

根據KKT條件，得：

所以得：

其中，至少有一個αj>0，(因為若所有的α都等於零，則w*等於零，則超平面不存在，矛盾)，對此，對j有：

帶入w*，得：

兩邊同乘yj，並考慮到(yj)2=1，得：

所以，我們可以由下式求得原始最優問題的解w*和b*：

那分離超平面就可以寫成：

這就是說，分離超平面只依賴於輸入x和訓練樣本的內積。

我們進一步可以發現，w*和b*只依賴於α>0的樣本點，而其他樣本點對w*和b*沒有影響。我們又通過KKT約束條件：

可知，對於α>0的樣本i，其必定滿足：

也就是說，樣本點i位於間隔邊界上，我們稱它們為支持向量。這就得到了，分離超平面只與支持向量有關。

通過上述討論，我們知道了線性可分數據的求解方式，那對於線性不可分數據我們應該怎麼求解呢？

2.線性不可分支持向量機

其實在實際情況中，數據並不都是線性可分的，如下圖(圖片來自網路，侵刪)：

對於這樣的數據集，我們剛剛的不等式約束不能對所有的樣本點都成立。這時候，我們需要引入軟間隔最大化的概念，將SVM推廣到線性不可分數據問題。

線性不可分意味著有某些樣本點(xi,yi)不能滿足函數間隔大於等於1的約束條件，為了解決這個問題，我們為每個樣本點引入一個鬆弛變數ξi>=0，使得函數間隔加上鬆弛變數大於等於1，得：

同時對每一個鬆弛變數支付一個代價，得到新的目標函數：

這裡C>0，稱為懲罰係數，C值大時表示對於錯誤分類的懲罰增大。

這時，線性不可分的問題就變成了如下優化問題：

上述優化問題的拉格朗日函數為：

和線性可分支持向量機一樣，對偶問題是極大極小問題。首先求L對w,b,ξ的極小，得：

得到：

帶入，得極小問題為：

再對其求極大，得到對偶問題：

消去μ，得：

上式就是線性不可分問題的對偶問題。

同線性可分問題的求解方式，我們得到原始問題的解：

這時候可以得到分離超平面：

以上就是求解線性分類問題的線性支持向量機的方法。

非線性支持向量機和核函數

在之前的討論中，我們用的數據都是線性可分的，即存在一個劃分超平面可以將訓練樣本正確分類。但是有時候，並不存在這樣的超平面能正確劃分兩類樣本，如下圖左圖所示(圖片來自於網路，侵刪)：

對於這樣的問題，我們可以將輸入空間映射到一個更高維的特徵空間，使得樣本在特徵空間內線性可分。如上圖的變換過程。

同時，我們可以證明，如果輸入空間是有限維的，即屬性有限，那麼一定存在一個高維空間使得屬性可分。

同樣，類似上面的線性可分數據集，我們得到如下優化問題：

但由於特徵空間的維數可能會很高，所以計算內積會比較困難。為了解決這個問題，我們定義了這樣一個函數：

即xi和xj在特徵空間的內積等於它們在輸入空間中通過函數k(；)計算的結果。函數k(；)就稱為核函數。有了這個函數，我們就不用去計算在高維空間的內積。

那麼優化函數就可以重寫為：

通過之前的討論，在不知道特徵映射的形式時，我們也就不知道什麼樣的核函數是合適的。於是，核函數的選取也就至關重要，如果核函數選取不佳，那麼分類器的分類效果也就不會很好。這時，經驗就比較重要。給大家一些常用的核函數有：線性核、高斯核、多項式核、拉普拉斯核、Sigmoid核等。

以上就是非線性支持向量機和核函數的基本概念。

SMO演算法

大家發現，我們前面在解優化問題中的w*和b*時，總是先假設了參數α*已經求得。那參數α*怎麼求呢？這就用到了SMO演算法。

《機器學習實戰》這本書中用代碼實現了SMO演算法，但是對於演算法原理沒有過多解釋，所以讀者可能會對代碼中的公式很疑惑，在這裡我給大家推一遍SMO演算法，大家敲代碼時可以對照著這部分的公式推導看。

序列最小最優化演算法(SMO演算法)是快速求解支持向量機的一種方法。它的實現手段是重複下面兩步直到收斂：

1. 選取一對需要更新的變數αi和αj；

2. 固定αi和αj之外的參數，求解獲得更新後的αi和αj。

我們要解如下優化問題：

我們首先選取兩個變數，不失一般性，選取變數α1和α2，其他變數αi固定，於是上述優化問題的子問題可以寫成：

我們首先分析一下約束條件，因為有兩個變數，約束可以用二維平面中的圖形表示,如下：

不等式約束使得兩個變數在[0,C]*[0,C]內。因此要求目標函數在一條平行於對角線的線段上取最優值。這就使得兩個變數的最優化問題實質上變成了單變數的最優化問題。不妨考慮為α2的最優化問題。

假設變數的初始解為α1old和α2old，最優解為α1new，α2new，所以：

由於α2new需要滿足不等式約束，所以其取值範圍需要滿足：

其中，L和H是上圖對角端點的界，由上圖可知：

首先先考慮沿著約束方向未經剪輯即未考慮L和H約束時α2的最優解α2{new,unc}；然後再求剪輯後α2的最優解α2{new}.

令：

這時候目標函數可以寫為：

將α1帶入上式，得到只有變數α2的目標函數：

對α2求導得：

令其為零，得到：

因為α2要滿足約束條件L<=α2<=H,得到α2{new}如下表達式：

之後再得到α1{new}：

這時候就得到了兩個參數的最優化解α1{new}和α2{new}。

變數選擇

1. 第一個變數的選擇

對於軟間隔支持向量機，其KKT條件為：

我們需要尋找的第一個變數就是違反KKT條件最嚴重的樣本點。我們先遍歷所有滿足條件0<αi<C的樣本點，即在間隔邊界上的支持向量點，看他們是否滿足KKT條件。如果這些樣本點都滿足，那麼遍歷整個數據集，找出不滿足KKT條件的點。

2. 第二個變數的選擇

當我們已經選好第一個變數α1，這時第二個變數選擇的標準就是希望能使α2有足夠大的變化。

由之前得到：

α2依賴於|E1-E2|，所以一種簡單的做法就是選擇能使|E1-E2|最大的α2。

在特殊情況時，若上述方法選擇的α2不能使目標函數有足夠的下降，那麼採用如下啟發式的規則繼續選擇α2。首先遍歷在間隔邊界上的支持向量點，直到找到一個α2能使目標函數有足夠的下降。若找不到合適的α2，那麼遍歷整個數據集；若仍找不到合適的α2，則放棄第一個α1，重新尋找一個α1.

3. 計算b和差值Ei

當我們完成了兩個變數的優化後，需要對b的值進行更新。

當我們更新完兩個變數後，Ei的值同樣需要更新：

重複上述步驟，直到α收斂，得到最優的α取值。

以上就是SMO演算法的主要部分。

支持向量回歸（SVR）

支持向量機不僅僅可以用於分類，其同樣可以用於回歸問題。

其和傳統回歸問題不同的地方在於，傳統的回歸模型是基於模型輸出f(x)和真實輸出y之間的差來計算損失，當且僅當f(x)和y完全相同時，損失才為零。

支持向量回歸假設我們能容忍f(x)和y之間有ε的偏差，當且僅當f(x)和y之間的差值的絕對值大於ε時，才計算損失。

這一部分由於不是今天的重點，只是想提醒大家，支持向量機同樣可以用於回歸問題。想要詳細了解的讀者可以參閱周志華老師《機器學習》P133.

以上就是支持向量機全部的原理部分。

代碼實現與注釋

這部分代碼內容涉及公式較多，建議配合SMO演算法推導的公式一起看會很容易理解。

1.SMO演算法中輔助函數

程序清單

from numpy import *# 讀取數據def loadDataSet(filename): dataMat=[];labelMat=[] fr=open(filename) for line in fr.readlines(): lineArr=line.strip().split( ) dataMat.append([float(lineArr[0]),float(lineArr[1])]) labelMat.append(float(lineArr[2])) return dataMat,labelMat# 選擇參數i之外的一個參數並返回def selectJrand(i,m): j=i while(j==i): j=int(random.uniform(0,m)) return j# 調整aj的值def clipAlpha(aj,H,L): if(aj>H): aj=H if(L>aj): aj=L return aj

2.簡化版SMO演算法

程序清單：

# 簡化版SMO演算法def smoSimple(dataMatIn,classLabels,C,toler,maxIter): # 得到矩陣，將labelMat矩陣轉置 dataMatrix=mat(dataMatIn);labelMat=mat(classLabels).transpose() # 得到dataMatrix的大小 b=0;m,n=shape(dataMatrix) # 初始化m*1的零矩陣 alphas=mat(zeros((m,1))) iter=0 while(iter<maxIter): # 初始化 alphaPairsChanged=0 # 對整個集合遍歷 for i in range(m): # 為第i個樣本的預測類別。這一行代碼看不懂的給大家一個鏈接，講的很清楚 # http://www.bubuko.com/infodetail-694615.html fxi=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:]))+b # 計算預測出來的和真實的誤差 Ei=fxi-float(labelMat[i]) # 如果誤差比較大且alpha在0到C之間 if(((labelMat[i]*Ei<-toler)and(alphas[i]<C))or((labelMat[i]*Ei>toler) and(alphas[i]>0))): # 隨機選取一個j j=selectJrand(i,m) # 同理，計算誤差 fxj=float(multiply(alphas,labelMat).T* (dataMatrix*dataMatrix[j,:]))+b Ej = fxj - float(labelMat[j]) alphaIold=alphas[i].copy() alphaJold=alphas[j].copy() # 保證alpha在0到C之間 if(labelMat[i]!=labelMat[j]): L=max(0,alphas[i]-alphas[j]) 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，下面這些公式看不懂沒關係，會在正文部分給大家推導，大家對照著看 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]*(Ei-Ej)/eta alphas[j]=clipAlpha(alphas[j],H,L) if(abs(alphas[i]-alphaJold)<0.00001): print(j is 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 alphaPairsChanged+=1 print("iter: %d i %d,pairs changed %d"%(iter,i,alphaPairsChanged)) if(alphaPairsChanged==0):iter+=1 else:iter=0 print(itertion number: %d% iter) return b,alphas

3.完整版SMO支持函數

程序清單：

# 完整版SMO支持函數class optStruct: def __init__(self,dataMatIn,classLabels,C,toler,Ktup): self.X=dataMatIn self.labelMat=classLabels self.C=C self.m=shape(dataMatIn)[0] self.alphas=mat(zeros((self.m,1))) self.b=0 # eCache中第一列為標誌位，第二列是誤差E值 self.eCache=mat(zeros((self.m,2))) self.tol=toler self.K=mat(zeros(self.m,self.m)) for i in range(self.m): self.K[:,i]=kernelTrans(self.X,self.X[i,:],kup)def calcEk(oS,k): # 計算預測值 fXk=float(multiply(oS.alphas,oS.labelMat).T* (oS.X*oS.X[k,:].T))+oS.b # 計算誤差 Ek=fXk-float(oS.labelMat[k]) return Ek# 選擇第二個alphadef selectJ(i,oS,Ei): maxK=-1;maxDeltaE=0;Ej=0 # 將輸入值在eCache中設置為有效 oS.eCache[i]=[1,Ei] # 返回非零E值的位置 validEcacheList=nonzero(oS.eCache[:,0].A)[0] if((len(validEcacheList))>1): # 選擇誤差E最大的 for k in validEcacheList: if k==i:continue Ek=calcEk(oS,k) deltaE=abs(Ei-Ek) if(deltaE>maxDeltaE): maxK=k maxDeltaE=deltaE Ej=Ek return maxK,Ej # 如果是第一次循環，就隨機挑選一個j else: j=selectJrand(i,oS.m) Ej=calcEk(oS,j) return j,Ej# 更新eCache中的誤差值def updateEk(oS,k): Ek=calcEk(oS,k) oS.eCache[k]=[1,Ek]

4.完整的SMO優化常式

程序清單：

# 完整的SMO的優化常式# 其和上面簡化版的基本一致，公式部分看不懂看正文def innerL(i,oS): Ei=calcEk(oS,i) if((oS.labelMat[i]*Ei<-oS.tol)and(oS.alphas[i]<oS.C)or (oS.labelMat[i] * Ei>oS.tol) and (oS.alphas[i] >0)): j,Ej=selectJ(i,oS,Ei) alphaIold=oS.alphas[i].copy() alphaJold=oS.alphas[j].copy() if(oS.labelMat[i]!=oS.labelMat[j]): L=max(0,oS.alphas[j]-oS.alphas[i]) 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 eta=(2.0*oS.X[i,:]*oS.X[j,:].T-oS.X[i,:]*oS.X[i,:].T- oS.X[j,:]*oS.X[j,:].T) if(eta>=0):print("eta>=0");return 0 oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/eta oS.alphas[j]=clipAlpha(oS.alphas[j],H,L) updateEk(oS,j) if((abs(oS.alphas[j]-alphaJold)<0.00001)): print(j not moving enough) return 0 oS.alphas[i]+=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j]) updateEk(oS,i) b1=(oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIold)* oS.X[i,:]*oS.X[i,:].T-oS.labelMat[j]* (oS.alphas[j]-alphaJold)*oS.X[i,:]*oS[j,:].T) b2 = (oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.X[j, :] * oS[j, :].T) if((0<oS.alphas[i])and(oS.C>oS.alphas[i])): oS.b=b1 elif((0<oS.alphas[j])and(oS.C>oS.alphas[j])):oS.b=b2 else:oS.b=(b1+b2)/2.0 return 1 else: return 0# 完整版SMO的外循環代碼def smoP(dataMatIn,classLabels,C,toler,maxIter,kTup=(lin,0)): oS.optStruct(mat(dataMatin),mat(classLabels).transpose(),C,toler) iter=0 entireSet=True;alphaPairsChanged=0 while((iter<maxIter)and((alphaPairsChanged>0)or(entireSet))): alphaPairsChanged=0 if(entireSet): for i in range(oS.m): alphaPairsChanged+=innerL(i,oS) print(fullSet iter : %d i:%d,pairs changed %d% (iter,i,alphaPairsChanged)) iter+=1 else: nonBoundIs=nonzero((oS.alphas.A>0)*(oS.alphas.A<C))[0] for i in nonBoundIs: alphaPairsChanged+=innerL(i,oS) print("non-bound,iter:%d i:%d,pairs changed %d"% (iter,i,alphaPairsChanged)) iter+=1 if(entireSet): entireSet=False elif(alphaPairsChanged==0): entireSet=True print(iteration number: %d%iter) return oS.b,oS.alphasdef calWs(alphas,dataArr,classLabels): X=mat(dataArr);labelMat=mat(classLabels).transpose() m.n=shape(X) w=zeros((n,1)) for i in range(m): w+=multiply(alphas[i]*labelMat[i],X[i,:].T) return w

以上就是支持向量機（SVM）全部的內容。這個演算法推導部分公式較多，但其實自己推導一遍並不是很難，同時也會加深自己對於演算法的理解。推導過程中一些注意要點我也都寫了出來，希望能給大家提供一點幫助。下一節一起梳理AdaBoost元演算法。

聲明

最後，所有資料均本人自學整理所得，如有錯誤，歡迎指正，有什麼建議也歡迎交流，讓我們共同進步！轉載請註明作者與出處。

以上原理部分主要來自於《機器學習》—周志華，《統計學習方法》—李航，《機器學習實戰》—Peter Harrington。代碼部分主要來自於《機器學習實戰》，代碼用Python3實現，這是機器學習主流語言，本人也會儘力對代碼做出較為詳盡的注釋。