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演算法分類非線性資料集，敬請期待~


作者：鄧莎
連結：

支援向量機|SMO演算法實現

01 起

02 SMO演算法原理

03 輔助函式

04 簡化版SMO演算法

4.1 程式碼

4.2 測試

05 完整版SMO演算法

5.1 程式碼

5.2 測試

06 總結

相關文章