【機器學習實戰-python3】支援向量機（Support Vecrtor Machines SVM）

工具：PythonCharm 書中的程式碼是python2的，而我用的python3，結合實踐過程，這裡會標註實踐時遇到的問題和針對python3的修改。
實踐程式碼和訓練測試資料可以參考這裡
https://github.com/stonycat/ML-in-Action
（原書作者也提供了原始碼，但是存在一些問題，且在python3中有部分修改）

有人認為SVM是最好的現成的分類器，“現成”指的是分類器不加修改即可直接使用，意味著直接應用SVM可以取得較低的錯誤率，對訓練集之外的資料點做出很好的分類決策。
SVM有許多實現，這裡介紹其中一種最流行的實現，即序列最小優化（SMO）演算法，然後新增kernel函式將SVM拓展到更多資料集。
SVM是基於最大間隔分隔資料，若所給資料是二維的，則分隔線為一條直線，若資料為三維的，則分割線為一個平面，依次類推，隨著資料維數的增加，分隔面就變成了超平面。而最大間隔，是讓離分隔面最近的點，確保他們離分隔面儘可能遠。SVM本身是一個二類分類器，若要解決多類問題，需要修改SVM。

一、尋找最大間隔
SVM中為了找到劃分資料集的最佳分隔，確保最近的點到分隔的垂線最短，因此轉化為了一個優化問題。這裡分隔面可以定義為：wT+b，將輸入資料給分類器，會輸出分隔標籤，這裡採用sigmoid函式，即f（wT+b），可輸出1/-1（與1或0無異）。間隔則通過label∗f（wT+b）來計算。找到最小間隔的資料點，最大化其與間隔的距離。

argmaxw,b{minn(label⋅(wT+b))⋅1||w||} 這裡
{minn(label⋅(wT+b))⋅1||w||}是點到間隔的幾何間隔，與上面所說的label∗f（wT+b）的函式間隔意義相同。這裡固定label∗f（wT+b）

>=1.0 去優化1||w||，採用拉格朗日乘子法。關於SVM的詳細推到公式，可參考其他相關材料：如PRML原書pdf版本

二、SMO高效優化演算法
John Platt所提出的SMO演算法用於訓練SVM，將最大優化問題分解為多個小優化問題，求解的結果是完全一致的，且SMO演算法求解的時間短很多。
為了得到分隔的超平面，就需要得到最優的alpha值，上面所提到的最優化式子可化為下式：

∑αi⋅labeli=0SMO演算法的目標就是求出一系列alpha和b，一旦求出alpha就很容易計算出權重向量從而得到分個超平面。SMO演算法的工作原理是，每次迴圈對兩個alpha進行優化處理，得到一對合適的alpha後，就增大其中一個同時減小另一個。這裡的“合適”指的是這兩個alpha必須在間隔邊界之外，並且兩個alpha還沒有經過區間化處理或者不在邊界上。下面是匯入資料的程式碼。

def 
 loadDataSet(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,labelMat

def selectJrand(i,m): #i表示alpha的下標，m表示alpha的總數
    j=i
    while(j==i):
        j=int(random.uniform(0,m)) #簡化版SMO，alpha隨機選擇
    return j

def clipAlpha(aj,H,L): #輔助函式，用於調整alpha範圍
    if aj>H:
        aj=H
    if L>aj:
        aj=L
    return aj

可以看出，類別標籤為1與-1

這裡寫圖片描述

1、下面進入簡化版的SMO
虛擬碼大致如下：（不同層次的迴圈用顏色匹配）
建立一個alpha向量並將其初始化為0向量
當迭代次數小於最大迭代次數(外迴圈)
對資料集中每個資料向量(內迴圈)：
如果該資料向量可以被優化：
隨機選擇另外一個數據向量，同時優化這個向量
如果兩個向量都不能被優化，退出內迴圈
如果所有向量都未被優化，增加迭代數，進入下次迴圈。
對SMO原理感興趣的可參考：http://blog.csdn.net/on2way/article/details/47730367

原始碼：

#SMO simple algorithm
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):#toler表示容錯率 常數C
    dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
    b = 0; m,n = shape(dataMatrix)
    alphas = mat(zeros((m,1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0   #  標記alpha是否被優化
        for i in range(m):
            # fXi是預測的類別
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            # Ei表示誤差
            Ei = fXi - float(labelMat[i])# 預測結果和真實結果比對，計算誤差
            # 對alpha進行優化，同時檢查alpha的值滿足兩個條件：if
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i,m)# 隨機選擇第二個alpha
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy(); # 分配記憶體 稍後比較誤差
                if (labelMat[i] != labelMat[j]): # 計算L H用於將alpha[j]調整到0—C之間
                    L = max(0, alphas[j] - alphas[i])
                    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為alpha[j]的最優修改量
                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)
                # 檢查alpha[j]是否有輕微改變
                if (abs(alphas[j] - alphaJold) < 0.00001): print ("j not moving enough"); continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
                                                                        #the update is in the oppostie direction
                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 ("iteration number: %d" % iter)
    return b,alphas

test result：
這裡寫圖片描述

2、下面利用完整的Platt SMO演算法加速優化
簡化版SMO對於小規模資料集可用，但對於大規模資料，執行速度會變慢。完整的SMO演算法通過一個外迴圈來選擇alpha值，選擇過程在兩種方式之間交替：
1-在所有資料集上進行單遍掃描； 2-在非邊界（不等於邊界0/C的值）alpha中實現單遍掃描。實現非邊界值的掃描時，需要建立alpha列表，然後對錶進行遍歷，跳過那些已知不會改變的alpha值。
選擇第一個alpha後，通過內迴圈來選擇第二個alpha，優化過程中選擇最大步長的方式獲取第二個alpha。
在寫改進程式碼前，需要對之前的簡化版SMO的程式碼進行清理：
下面包含用於清理程式碼的資料結構和3個用對E快取的輔助函式。

class optStruct:  #儲存所有重要值，實現對成員變數的填充
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #誤差快取
        #self.K = mat(zeros((self.m,self.m)))
        #for i in range(self.m):
         #   self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

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

def selectJ(i,oS,Ei): #選擇第二個alpha值以保證每次優化的最大步長（內迴圈）change English~
    maxK=-1; maxDeltaE=0; Ej=0
    oS.eCache[i]=[1,Ei] #input Ei
    validEcacheList=nonzero(oS.eCache[:,0].A)[0] #set Ei valid in cache，nonzero is a list
    #choose the num which change most, if the first choice: use random
    if(len(validEcacheList))>1:
        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
    else:
        j=selectJrand(i,oS.m)
        Ej=calcEk(oS,j)
    return j,Ej

def updateEk(oS, k):#after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]

以上程式碼本身作用不大，當和優化過程及外迴圈組合在一起時，能組成SMO演算法。
3、優化過程：（選擇第二個alpha）與之前簡化版smoSimple差別不大，不過添加了自己的資料結構，在oS中傳遞，使用selectJ代替selectJrand來選擇則第二個alpha，alpha改變時更新Ecache值。

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) #this has been changed from selectJrand
        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.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
        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) #added this for the Ecache
        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])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        b1 = 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]
        b2 = 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]
        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

4、完整版SMO外迴圈程式碼（包含了kernel函式後面介紹）

def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
    m,n = shape(X)
    K = mat(zeros((m,1)))
    if kTup[0]=='lin': K = X * A.T   #linear kernel
    elif kTup[0]=='rbf':
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab
    else: raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')
    return K

class optStruct:  #儲存所有重要值，實現對成員變數的填充
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #誤差快取
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

def calcEk(oS, k):#計算誤差
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

def selectJ(i,oS,Ei): #選擇第二個alpha值以保證每次優化的最大步長（內迴圈）change English~
    maxK=-1; maxDeltaE=0; Ej=0
    oS.eCache[i]=[1,Ei] #input Ei
    validEcacheList=nonzero(oS.eCache[:,0].A)[0] #set Ei valid in cache，nonzero is a list
    #choose the num which change most, if the first choice: use random
    if(len(validEcacheList))>1:
        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
    else:
        j=selectJrand(i,oS.m)
        Ej=calcEk(oS,j)
    return j,Ej

def updateEk(oS, k):#after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]

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) #this has been changed from selectJrand
        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.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
        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) #added this for the Ecache
        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])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        b1 = 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]
        b2 = 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]
        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

#full Platt SMO
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    # 當迭代超過maxIter或者整個資料集都未對任意alpha進行修改時，退出迴圈
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   #go over all
            for i in range(oS.m):
                alphaPairsChanged += innerL(i,oS)#use innerL to choose the second alpha
                print ("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1
        else:#go over non-bound (railed) alphas
            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 #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True
        print ("iteration number: %d" % iter)
    return oS.b,oS.alphas

#由alpha得到權值w

def calcWs(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

迭代結果和顯示權值w：
這裡寫圖片描述

5、利用kernel函式將資料對映
低維空間無法用線性函式分類的資料，需要使用kernel trick/kernel substation將資料對映到高維空間，然後可採用線性的方法，等價於在低維空間處理非線性資料。SVM所有運算可用內積表示，易轉化為kernel函式。核函式不僅僅應用於SVM，後續機器學習演算法也有涉及。
這裡選擇了一種比較流行的核函式，徑向基核函式

k(x,y)=exp∥x−y∥22

【機器學習實戰-python3】支援向量機（Support Vecrtor Machines SVM）

【機器學習實戰-python3】支援向量機（Support Vecrtor Machines SVM）

【機器學習實戰】第6章支援向量機（Support Vector Machine / SVM）

【機器學習筆記17】支援向量機

【機器學習演算法推導】支援向量機

【機器學習實戰-python3】線性迴歸

機器學習實戰筆記5—支援向量機

【深度學習基礎-06】支援向量機SVM（下）-線性不可分

【深度學習基礎-05】支援向量機SVM（上）-線性可分

機器學習筆記09：支援向量機（二）-核函式(Kernels)

機器學習之支援向量機： Support Vector Machines (SVM)

支援向量機（Support Vector machines）

【機器學習】支援向量機（4）——非線性支援向量機（核函式）

【機器學習實戰之一】：C++實現K-近鄰演算法KNN

【機器學習實戰系列】讀書筆記之AdaBoost演算法公式推導和例子講解（一）

[完]機器學習實戰第六章支援向量機（Support Vector Machine）

【機器學習實戰系列】讀書筆記之KNN演算法（三）

Andrew Ng 機器學習筆記 11 ：支援向量機(Support Vector Machine)

Stanford機器學習第八講支援向量機SVM

機器學習之旅：支援向量機通俗導論（理解SVM的三層境界）

自學機器學習之sklearn實現支援向量機

【機器學習實戰-python3】支援向量機（Support Vecrtor Machines SVM）

相關推薦