1. 程式人生 > >【機器學習】SVM基礎知識+程式碼實現

【機器學習】SVM基礎知識+程式碼實現

 

1. 基本知識

二分類:通過分離超平面對資料點進行分類,訓練分離超平面。

原理:最大化支援向量到分離超平面的距離。支援向量:離分離超平面最近的點。

2. 完全線性可分(硬間隔)

2.1 SVM基本型

分離超平面:w^Tx + b。(訓練中更新w和b,或alpha,使得分離超平面分類效果最佳)

某點到分離超平面的函式距離:y_i(w^T\mathbf{x}_i + b)

某點到分離超平面的幾何距離:\frac{y_i(w^T\mathbf{x}_i + b)}{||w||}, ||w||為w的L2範數。

點集到分離超平面的幾何距離 => 距離超平面最近的點到其的距離:\min_{n}{\frac{y_i(w^T\mathbf{x}_i + b)}{||w||}}

SVM的目標:

arg\max_{w,b}\min_{n}{\frac{y_i(w^T\mathbf{x}_i + b)}{||w||}}

引入約束條件,假設所有點到超平面的距離都大於等於1,其中裡分離超平面距離為1的點稱為‘支援向量’(即等號成立時的點)。則問題轉化為:

arg\max_{w,b}{\frac{1}{||w||}}, \ y(w^T\mathbf{x}+b) \geq 1

arg\min_{w,b}{\frac{1}{2}}||w||^2, \ y(w^T\mathbf{x}+b) \geq 1

2.2 對偶問題

求解:arg\min_{w,b}{\frac{1}{2}}||w||^2, \ y(w^T\mathbf{x}+b) \geq 1,以求得分離超平面,為凸二次規劃問題。

引入拉格朗日乘子\alpha_i \geq 0,則拉格朗日函式為:

L(w,b,\alpha) = \frac{1}{2}||w||^2 + \sum_{i=1}^{n}\alpha_i(1-y_i(w^T\mathbf{x}_i + b) ) , \ \alpha \geq 0

L(w,b,\alpha) = \frac{1}{2}||w||^2 - \sum_{i=1}^{n}y_i(w^T\mathbf{x}_i + b) + \sum_{i=1}^n\alpha_i, \ \alpha \geq 0

分別對w,b求偏導,使導數為零:

\\w = \sum_{i=1}^n\alpha_i y_i \mathbf{x}_i \\0 = \sum_{i=1}^n \alpha_iy_i

將w和b由alpha表示,並新增約束條件,最後使得問題轉化為:

\\ {\color{DarkBlue} \max_{\alpha} -\frac{1}{2}\sum_{i=1}^n\sum_{j=1}^n\alpha_i\alpha_jy_iy_j(\mathbf{x}_i ^T\mathbf{x}_j) + \sum_{i=1}^n\alpha_i } \\ {\color{DarkBlue} s.t. \ \ \alpha\geq 0, \ \sum_{i=1}^n\alpha_iy_i = 0}

最後得到:

f(\mathbf{x}) = \mathbf{w}^T\mathbf{x}+b = \sum_{i=1}^n\alpha_iy_i\mathbf{x}_i \mathbf{x} + b

上述過程滿足KKT條件:

\\KKT: \\ {\color{DarkBlue} \alpha_i \geq 0} \\ {\color{DarkBlue} y_if(\mathbf{x}_i) \geq 1} \\ {\color{DarkBlue} \alpha_i(y_if(\mathbf{x}_i)-1) = 0 }

上述KKT條件說明,當alpha為0時,對應的資料點不參與w的計算,而當alpha不為0時,y_if(\mathbf{x}_i)為1,說明這些點到分離超平面距離為1,為支援向量。這說明:訓練完成後,只需要保留支援向量的樣本即可。

3. SMO演算法

SMO:sequential minimal optimization 序列最小優化。優化求解alpha。

基本思路:固定\alpha_i之外的所有引數,然後求\alpha_i上的極值。由於約束條件:\sum_{i=1}^n\alpha_iy_i = 0,則固定\alpha_i之外的其他變數後,\alpha_i可以直接求解。

實驗步驟

step1:選取一對需要更新的\alpha_i\alpha_j

step2:固定\alpha_i\alpha_j以外的引數,求解上述藍色部分的對偶問題,獲得更新後的\alpha_i\alpha_j

\alpha_i的選取規則:第一個變數選取違背KKT條件程度最大的變數(KKT越違背,則更新alpha後目標函式的增值越大)。

\alpha_j的選取規則:第二個變數選取使目標數值增長最快的變數。

SMO中啟發式方法:選取的兩個變數\alpha_i\alpha_j對應的兩個樣本(\mathbf{x}_i,y_i),\ (\mathbf{x}_j,y_j)之間的間隔最大,因為對差別大的兩個變數進行更新,會帶給目標函式值更大的變化。

當計算得到\alpha後,因為對於所有支援向量(\mathbf{x}_i,y_i)y_if(\sum_{j=1}^S\alpha_jy_j\mathbf{x}_jx_i + b) = 1, 其中S為所有支援向量的集合。因此可以將一個支援向量帶入求得b值,或者求所有支援向量的b值然後取平均。

總結:

1. 初始化\alpha為全1向量。

2. 第一輪兩個\alpha只能隨機選擇,因為都為1.

3. 以後幾輪:更新\alpha時,外迴圈遍歷所有資料集或非邊界\alpha(不等於0或C的\alpha),選取第一個\alpha_i。然後根據最大化誤差步長來選取第二個\alpha_j, 根據目標函式及其約束條件,更新這兩個\alpha

4. 重複上述步驟知道最大迭代步數,或所有\alpha都不違背KKT條件。

4. 核函式 

當資料線性不可分時,利用核函式將低維度下線性不可分的資料對映到高維度下線性可分的資料。即將資料從舊特徵空間對映到新的特徵空間。因為SVM中的運算可寫成內積的形式(內積後得數為一個標量or數值),所以可以直接用核函式來代替原內積。

若資料完全線性不可分,則引入核函式,即對輸入空間進行轉換到特徵空間:

K(\mathbf{x}_i, \mathbf{x}_j) = f(\mathbf{x}_i)f(\mathbf{x}_j),將原式中的x的內積轉換成核函式。

線性核函式:

k(\mathbf{x}_i, \mathbf{x}_j) = \mathbf{x}_i^T\mathbf{x}_j

多項式核函式:

k(\mathbf{x}_i, \mathbf{x}_j) = (\mathbf{x}_i^T\mathbf{x}_j)^d,d 為多項式次數

高斯核(徑向基核函式):

k(\mathbf{x}_i, \mathbf{x}_j) = exp(-\frac{||\mathbf{x}_i - \mathbf{x}_j||^2}{2\sigma^2}), sigma > 0, 為高斯核的頻寬

拉普拉斯核:

k(\mathbf{x}_i, \mathbf{x}_j) = exp(-\frac{||\mathbf{x}_i - \mathbf{x}_j||}{\sigma}),sigma > 0

sigmoid 核:

k(\mathbf{x}_i, \mathbf{x}_j) = tanh(\beta\mathbf{x}_i^T\mathbf{x}_j + \theta),\ \beta > 0, \theta < 0

5. 非完全線性可分(軟間隔)

若資料非完全線性可分,則引入鬆弛變數,,允許少部分資料點處於分割面錯誤一側。

優化目標:

\\ arg\min_{w,b}{\frac{1}{2}}||w||^2 + C\sum_{i=1}^n l_{0/1}(y_i(\mathbf{w}^T\mathbf{x}_i + b)-1), \\ s.t. \ y(w^T\mathbf{x}+b) \geq 1

其中l_{0/1}為“0/1損失函式”:即當某點到分離超平面距離小於1時,考慮常數C,否則還按照線性可分的硬間隔來做。

l_{0/1}(z) = \left\{\begin{matrix} 1 & if \ z < 0 \\ 0& otherwise \end{matrix}\right.

一般使用一些替代損失來代替上述“非凸,非連續的0/1損失”,一般選擇有:

hinge損失:l_{hinge}(z) = max(0, 1-z)

指數損失:l_{exp} (z) = exp(-z)

對數損失:l_{log} (z) = log(1 + exp(-z))

若引入鬆弛變數,目標優化函式為:

\\ arg\min_{w,b}{\frac{1}{2}}||w||^2 + C\sum_{i=1}^n \xi _i \\ s.t. \ y(w^T\mathbf{x}+b) \geq 1-\xi_i,\ \xi_i > 0

因此對應的拉格朗日函式:

\\ L(\mathbf{w}, b, \mathbf{\xi}, \mathbf{\alpha}, \mathbf{\mu}) = \frac{1}{2} ||w||^2 + \sum_{i=1}^n \alpha_i(1-\xi_i-y_i(\mathbf{w}^T\mathbf{x}_i + b)) - \sum_{i=1}^n\mu_i\xi_i \\ s.t. \ \alpha \geq 0, \ \mu \geq 0, \xi \geq 0

w,\ b, \ \xi求偏導為0,求得:

\\ \mathbf{w} = \sum_{i=1}^n \alpha_iy_i\mathbf{x}_i; \\ 0 = \sum_{i=1}^n \alpha_iy_i; \\ C = \alpha_i + \mu_i

因此原目標函式的對偶問題為:

\\ {\color{DarkBlue} \max_{\alpha} -\frac{1}{2}\sum_{i=1}^n\sum_{j=1}^n\alpha_i\alpha_jy_iy_j(\mathbf{x}_i ^T\mathbf{x}_j) + \sum_{i=1}^n\alpha_i }\\ {\color{DarkBlue} s.t. \ \ C\geq \alpha\geq 0, \ \sum_{i=1}^n\alpha_iy_i = 0}

軟間隔支援向量機的KKT條件為:

\\ \alpha_i \geq 0, \ \mu_i \geq 0 \\ y_i(\mathbf{w}\mathbf{x}_i + b) - 1 - \xi_i \geq 0 \\ \alpha_i(y_i(\mathbf{w}\mathbf{x}_i + b) - 1 - \xi_i) = 0 \\ \xi_i \geq 0, \ \mu_i\xi_i = 0

6. 合頁損失函式形式的目標函式

原目標函式:

\\ \min_{w,\xi,b} \sum_{i=1}^n \frac{1}{2}||w||^2 + C\sum_{i=1}^n \xi_i \\ s.t. \ y_i(\mathbf{w}^T\mathbf{x}_i + b) \geq 1 - \xi_i,\ \xi \geq 0

新目標函式:

\min ([1 - y_i(\mathbf{w}^T\mathbf{x}_i + b)]_{+} + \lambda||w||^2)

其中:

\lambda = \frac{1}{2C}

合頁損失函式是指:

[z ]_+= \left\{\begin{matrix} z & z > 0\\ 0 & otherwise \end{matrix}\right.

 

7. 程式碼實現

參考:《機器學習實戰》

原始碼地址以及資料:https://github.com/JieruZhang/MachineLearninginAction_src

簡化SMO:

from numpy import *
import random
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

#在簡單版SMO中,alpha是隨機選擇的
def selectJrand(i,m):
    j=i
    while j==i:
        j = int(random.uniform(0,m))
    return j

#輔助函式,調整過大過小值 
def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    elif aj < L:
        aj = L
    return aj

#簡單版的SMO演算法,隨機選擇兩個alpha,若滿足優化條件,則進行優化
def smoSimple(dataMatIn, classLabels, C, tol, maxIter):
    '''
    dataMatIn: 輸入資料集
    classLabels: 類別標籤
    C: 常數
    tol: 容錯率
    maxIter:最大迴圈數
    '''
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    b = 0
    m,n = shape(dataMatrix)
    alphas = mat(zeros((m,1)))
    iters = 0
    #外迴圈
    while iters < maxIter:
        alphaPairChanged = 0
        #內迴圈,對於每一個例項,對應的alpha為第一個alpha
        for i in range(m):
            #預測類別
            fxi = float(multiply(alphas,labelMat).T * (dataMatrix*dataMatrix[i,:].T)) + b
            #分類誤差
            Ei = fxi- float(labelMat[i])
            #當誤差和alpha的值滿足優化條件時,進行優化,正負間隔均已考慮
            if ((labelMat[i]*Ei < -tol) and (alphas[i] < C) or (labelMat[i]*Ei > tol) and (alphas[i]) > 0):
                #因為alphas的和需要為零,所以要再選一個alpha進行優化
                #隨機選擇第二個alpha
                j = selectJrand(i,m)
                fxj = float(multiply(alphas,labelMat).T * (dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fxj- float(labelMat[j])
                #儲存舊的alpha的值,因為python中傳的是引用,所以使用copy函式分配新的記憶體。
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                #設定L和H的值,用於確保alpha在0到C之間
                if labelMat[i] != labelMat[j]:
                    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])
                #如果L和H相等,則不做任何改變,跳出迴圈
                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
                #如果eta不滿足條件,則退出當前迴圈
                if eta >= 0:
                    print('eta>=0')
                    continue
                #更新alpha[j]
                alphas[j] -= labelMat[j]*(Ei -Ej)/eta
                #限制alpha[j]在0到C之間
                alphas[j] = clipAlpha(alphas[j], H, L)
                #如果alpha[j]該變數較小,則不做改變,直接退出迴圈,重新隨機選擇j
                if abs(alphas[j]-alphaJold) < 0.00001:
                    print('j is not moving enough')
                    continue
                #更新alpha[i], 與j的更新方向相反
                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 alphas[i] > 0 and alphas[i] < C:
                    b = b1
                elif alphas[j] > 0 and alphas[j] < C:
                    b = b2
                else:
                    b = (b1+b2)/2.0
                alphaPairChanged += 1
                print('iter: %d i:%d, pairs changed %d'%(iters, i, alphaPairChanged))
        #如果alpha未更新過,則進行下一次迭代,知道alpha更新或達到最大迭代次數
        if alphaPairChanged == 0:
            iters += 1
        else:
            iters = 0
        print('iteration number: %d'%iters)
    return b, alphas

#測試上述函式
dataArr, labelArr = loadDataSet('testSet6.txt')
b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)

 

原版SMO:

#完整版SMO,啟發式方法選取待更新的兩個alpha

#核轉換函式
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
        #用於快取誤差,幫助尋找最佳的第二個引數alphas[j]
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

#已知alpha,計算誤差E,
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

#啟發式選取使得誤差最大的第二個變數alpha[j]
def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    #找出合法的誤差E值,即非零值
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]
    #如果誤差快取中非零值存在,則找到使得誤差最大的第二個變數的下標j
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            if k == i: continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    #若是第一次執行,則誤差快取中為空,只能先隨機選取alpha[j]
    else:   #in this case (first time around) we don't have any valid eCache values
        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]

#內迴圈,選擇第二個alpha[j]並更新所選的兩個變數
#大體上和simpleSMO內迴圈類似,只是選擇alpha[j]時使用啟發式方法,而不是隨機選擇。此外,更新alpha後,會將其對應的誤差存入快取,用作以後使用。
#這裡使用了核函式
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)):
        #啟發式選擇alpha[j]
        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)
        #更新alpha後,將對應的誤差存入快取
        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
        #更新alpha後,將對應的誤差存入快取
        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

#外迴圈
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
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        #遍歷整個資料集對應的alpha[i]
        if entireSet:   #go over all
            for i in range(oS.m):        
                alphaPairsChanged += innerL(i,oS)
                print ("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1
        #僅遍歷不在邊界上的alpha[i],即大於0小於C
        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))
    #雖然遍歷了整個資料集,但是隻有支援向量發揮了作用,只有支援向量的alpha不為0
    for i in range(m):
        w += multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w

#測試一下資料
dataArr, labelArr = loadDataSet('testSet6.txt')
b, alphas = smoP(dataArr,labelArr,0.6,0.001,40)
dataMat = mat(dataArr)
ws = calcWs(alphas,dataArr,labelArr)
print('The predicted label is:  ',dataMat[0]*mat(ws)+b)
print('The actual label is: ', labelArr[0])