1. 程式人生 > >【機器學習】演算法原理詳細推導與實現(五):支援向量機(下)

【機器學習】演算法原理詳細推導與實現(五):支援向量機(下)

【機器學習】演算法原理詳細推導與實現(五):支援向量機(下)

上一章節介紹了支援向量機的生成和求解方式,能夠根據訓練集依次得出\(\omega\)、\(b\)的計算方式,但是如何求解需要用到核函式,將在這一章詳細推導實現。

核函式

在講核函式之前,要對上一章節得到的結果列舉出來。之前需要優化的凸函式為:

\[ min_{\gamma,\omega,b}->\frac{1}{2}||\omega||^2 \]

\[ y^{(i)}(\omega^Tx^{(i)}+b) \geq 1 ,i=1,2,...,m \]

這裡假設資料是線性可分隔的,對於這個優化專案,給定一個訓練集合,這個問題的演算法會找到一個數據集合的最優間隔分類器,可以使訓練樣本的幾何間隔最大化。

在上一章節【機器學習】演算法原理詳細推導與實現(四):支援向量機(上)中,我們推出了這個問題的對偶問題,也就是要使這個式子最大化:

\[ max_{\alpha}\Gamma(\omega,b,\alpha)=\sum_{i=1}^m\alpha_i-\frac{1}{2}\sum^m_{i=1,j=1}\alpha_i\alpha_jy^{(i)}y^{(j)}<x^{(i)},x^{(j)}> \]

\[ \alpha_i \geq 0,i=1,2,...,m \]

\[ \sum_{i=1}^m\alpha_iy^{(i)}=0 \]

上面是我們的原始問題,且根據拉格朗日對偶步驟計算得到引數\(\omega\):

\[ \omega=\sum^m_{i=1}\alpha_iy^{(i)}x^{(i)} \]

\[ b=\frac{max_{i:y^{(i)}=-1}\omega^Tx^{(i)}+min_{i:y^{(i)}=1}\omega^Tx^{(i)}}{2} \]

當需要做分類預測時,需要對新來的輸入值\(x\)進行計算,計算其假設的值是否大於零,也就是做一次線性運算來判斷是正樣本還是負樣本,有如下計算函式:

\[ \begin{split} h_{\omega,b}(x)&=g(\omega^Tx+b) \\ &=g(\sum^m_{i=1}\alpha_iy^{(i)}<x^{(i)},x>+b) \end{split} \]

核函式概念

接下來要介紹“核”的概念,這個概念具有這樣的性質:

演算法對於x的依賴僅僅侷限於這些內積的計算,甚至在整個演算法中,都不會直接使用到向量x的值,而是隻需要用到訓練樣本與輸入特徵向量的內積

而“核”的概念是這樣的,考慮到最初在【機器學習】演算法原理詳細推導與實現(一):線性迴歸中提出的問題,比如有一個輸入\(x\in R\)是房屋的面積,\(y\)是房子的價格。假設我們從樣本點的分佈中看到\(x\)和\(y\)符合3次曲線,那麼我們會希望使用\(x\)的三次多項式來逼近這些樣本點。首先將特徵\(x\)擴充套件到三維\((x,x^2,x^3)\),這裡將這種特徵變換稱作特徵對映,對映函式為\(\varphi(x)\):

\[ \varphi(x)=\begin{bmatrix} x \\ x^2 \\ x^3 \end{bmatrix} \]

用\(\varphi(x)\)代表原來的特徵\(x\)對映成的,這裡希望得到對映後的特徵應用於svm分類,而不是最初的一維特徵,只需要將前面\(\omega^Tx+b\)公式中的內積從\(<x^{(i)},x^{(j)}>\)對映到\(<\varphi(x)^{(i)},\varphi(x)^{(j)}>\)。至於為什麼需要對映後的特徵而不是最初的特徵來參與計算,上面提到的一個原因:為了更好的擬合,另外一個原因是樣本可能存線上性不可分的情況,而特徵對映到高維過後往往就可分了。

如果原始特徵的內積為\(<x,z>\),對映後為\(<\varphi(x),\varphi(z)>\),那麼一般核函式定義為:

\[ K(x,z)=\varphi(x)^T\varphi(z) \]

為什麼會那麼定義核函式?有些時候\(\varphi(x)\)的維度將會非常的高,可能會包含非常高維的多項式特徵,甚至會到無限維。當\(\varphi(x)\)的維度非常高時,可能無法高效的計算內積,甚至無法計算。如果要求解前面所提到的凸函式,只需要先計算\(\varphi(x)\),然後再計算\(\varphi(x)^T\varphi(z)\)即可,但是這種常規方法是很低效的,比如最開始的特徵是\(n\)維,並將其對映到\(n^2\)維度,這時候計算需要\(O(n^2)\)的時間複雜度。這裡假設\(x\)和\(z\)都是\(n\)維的:

\[ K(x,z)=(x^Tz)^2 \]

展開後得到:

\[ \begin{split} K(x,z)&=(x^Tz)^2 \\ &=(\sum^n_{i=1}x_iz_i)(\sum^n_{j=1}x_jz_j) \\ &=\sum^n_{i=1}\sum^n_{j=1}x_ix_jz_iz_j \\ &=\sum^n_{i=1}\sum^n_{j=1}(x_ix_j)(z_iz_j) \\ &=\varphi(x)^T\varphi(z) \end{split} \]

也就是說,如果開始的特徵是\(n\)維,並將其對映到\(n^2\)維度後,其對映後的計算量為\(O(n^2)\)。而如果只是計算原始特徵\(x\)和\(z\)的內積平方,時間複雜度還是\(O(n)\),其結果等價於對映後的特徵內積。

回到之前的假設,當\(n=3\)時,這個核\(K(x,z)\)對應的特徵對映\(\varphi(x)\)為:

\[ \varphi(x)=\begin{bmatrix} x_1x_1 \\ x_1x_2 \\ x_1x_3 \\ x_2x_1 \\ x_2x_2 \\ x_2x_3 \\ x_3x_1 \\ x_3x_2 \\ x_3x_3 \end{bmatrix} \]

這是時間複雜度為\(O(n^2)\)計算方式,而如果不計算\(\varphi(x)\),直接計算\(<x,z>\)從而得到<\(\varphi(x)\),\(\varphi(z)\)>的內積,時間複雜度將縮小\(O(n)\)。

同理將核函式定義為:

\[ \begin{split} K(x,z)&=(x^Tz+c) \\ &=\sum^n_{i,j=1}(x_ix_j)(z_iz_j)+\sum^n_{i=1}(\sqrt{2cx_i})(\sqrt{2cx_j})+c^2 \end{split} \]

當\(n=3\)時,這個核\(K(x,z)\)對應的特徵對映\(\varphi(x)\)為:

\[ \varphi(x)=\begin{bmatrix} x_1x_1 \\ x_1x_2 \\ x_1x_3 \\ x_2x_1 \\ x_2x_2 \\ x_2x_3 \\ x_3x_1 \\ x_3x_2 \\ x_3x_3 \\ \sqrt{2c}x_1 \\ \sqrt{2c}x_2 \\ \sqrt{2c}x_3 \\ c \end{bmatrix} \]

總結來說,核的一種一般化形式可以表示為:

\[ K(x,z)=(x^Tz+c)^d \]

對應著\(\begin{bmatrix} n+d \\ d \end{bmatrix}\) 個特徵單項式,即特徵維度。

假如給定一組特徵\(x\),將其轉化為一個特徵向量\(\varphi(x)\);給定一組特徵\(z\),將其轉化為一個特徵向量\(\varphi(z)\),所以核計算就是兩個向量的內積\(<\varphi(x),\varphi(z)>\)。如果\(\varphi(x)\)和\(\varphi(z)\)向量夾角越小,即兩個向量越相似(餘弦定理),那麼\(\varphi(x)\)和\(\varphi(z)\)將指向相同的方向,因此內積會比較大;相反的如果\(\varphi(x)\)和\(\varphi(z)\)向量夾角越大,即兩個向量相似度很低,那麼\(\varphi(x)\)和\(\varphi(z)\)將指向不同的方向,因此內即將會比較小。

如果有一個核函式如下:

\[ K(x,z)=exp^{(-\frac{||x-z||^2}{2\sigma^2})} \]

如果\(x\)和\(z\)很相近(\(||x-z||\approx0\)),那麼核函式的值為1;如果\(x\)和\(z\)相差很大(\(||x-z||>>0\)),那麼核函式的值約等於0。這個核函式類似於高斯分佈,所以稱為高斯核函式,能夠把原始特徵對映到無窮維。

在前面說了:為什麼需要對映後的特徵而不是最初的特徵來參與計算?

上面提到了兩個原因:

  1. 為了更好的擬合
  2. 樣本可能存線上性不可分的情況,而特徵對映到高維過後往往就可分了

第二種情況如下所示:

左邊使用線性的時候,使用svm學習出\(\omega\)和\(b\)後,新來樣本\(x\)就可以代入到\(\omega^Tx+b\)中進行判斷,但是像圖中所示是無法判斷的;如果使用了核函式過後,\(\omega^Tx+b\)變成了\(\omega^T\varphi(x)+b\),直接可以用下面的方式計算:

\[ \begin{split} \omega^Tx+b&=(\sum^m_{i=1}\alpha_iy^{(i)}x^{(i)})^Tx+b \\ &=\sum^m_{i=1}\alpha_iy^{(i)}<x^{(i)},x>+b \\ &=\sum^m_{i=1}\alpha_iy^{(i)}K(x^{(i)})+b \end{split} \]

只需要將\(<x^{(i)},x>\)替換成\(K(x^{(i)})\)就能將低維特徵轉化為高維特徵,將線性不可分轉化成高維可分。

規則化和不可分情況處理

我們之前討論的情況都是建立在樣例線性可分的假設上,當樣例線性不可分時,我們可以嘗試使用核函式來將特徵對映到高維,這樣很可能就可分了。然而,對映後我們也不能100%保證可分。那怎麼辦呢,我們需要將模型進行調整,以保證在不可分的情況下,也能夠儘可能地找出分隔超平面。

看下面的圖可以解釋:

在右邊的圖可以可以看到上面一個離群點(可能是噪聲),會造成超平面的移動改變,使集合間隔的間隔距離縮小,可見以前的模型對噪聲非常敏感。再有甚者,如果離群點在另外一個類中,那麼這時候就是線性不可分了。 這時候我們應該允許一些點遊離在模型中違背限制條件(函式間隔大於 1)。我們設計得到新的模型如下(也稱軟間隔):

\[ min_{\gamma,\omega,b}->\frac{1}{2}||\omega||^2+C\sum^m_{i=1}\xi_i \]

\[ y^{(i)}(\omega^Tx^{(i)}+b) \geq 1-\xi_i ,i=1,2,...,m \]

\[ \xi_i \geq 0,i=1,...,m \]

引入非負引數\(\xi_i\)(鬆弛變數)過後,也就意味著允許某些樣本的函式間隔小於1,甚至是負數,負數就代表樣本點在對方區域中,如上方右邊圖的虛線作為超平面,一個空心圓點的函式間隔為負數。

增加新的條件後,需要重新調整目標函式,增加對離群點進行處罰,也就是在求最小值的目標函式後面加上\(C\sum^m_{i=1}\xi_i\),因為定義\(\xi_i \geq 0\),所以離群點越多,那麼目標函式的值越大,就等於違背求最小值的初衷。而\(C\)是離群點的權重,\(C\)越大表明離群點對於目標函式的影響越大,也就是越不希望看到離群點。

修改目標函式後,原式子變成:

\[ \Gamma(\omega,b,\xi,\alpha,r)=\frac{1}{2}\omega^T\omega+C\sum^m_{i=1}\xi_i-\sum^m_{i=1}\alpha_i[y^{(i)}(x^T\omega+b)-1+\xi_i]-\sum^m_{i=1}r_i\xi_i \]

這裡的\(\alpha\)和\(r\)都是拉格朗日運算元,根據上一章節拉格朗日的求解步驟:

  1. 構造出拉格朗日函式後,將其看作是變數\(\omega\)和\(b\)的函式
  2. 分別對其求偏導,得到\(\omega\)和\(b\)的表示式
  3. 然後帶入上述拉格朗日式子中,求帶入後式子的極大值

最後化簡得到的結果是:

\[ max_{\alpha}W(\alpha)=\sum_{i=1}^m\alpha_i-\frac{1}{2}\sum^m_{i=1,j=1}\alpha_i\alpha_jy^{(i)}y^{(j)}<x^{(i)},x^{(j)}> \]

\[ C \geq \alpha_i \geq 0,i=1,2,...,m \]

\[ \sum_{i=1}^m\alpha_iy^{(i)}=0 \]

這裡唯一不同的地方是限制條件多了一個離群點的權重\(C\)。

SMO優化演算法

SMO 是一個求解對偶問題的優化演算法,目前還剩下最後的對偶問題還未解決:

\[ max_{\alpha}W(\alpha)=\sum_{i=1}^m\alpha_i-\frac{1}{2}\sum^m_{i=1,j=1}\alpha_i\alpha_jy^{(i)}y^{(j)}<x^{(i)},x^{(j)}> \]

\[ C \geq \alpha_i \geq 0,i=1,2,...,m \]

\[ \sum_{i=1}^m\alpha_iy^{(i)}=0 \]

我們需要根據上述問題設計出一個能夠高效解決的演算法,步驟如下:

  1. 首先選擇兩個要改變的\(\alpha\)值\(\alpha_i\)、\(\alpha_j\)
  2. 其次保持除了\(\alpha_i\)、\(\alpha_j\)之外的所有引數固定
  3. 最後同時相對於這兩個引數使\(\omega\)取最優,且同時滿足所有約束條件

怎樣在滿足所有約束條件的情況下,相對於選出來的兩個引數\(\alpha_i\)、\(\alpha_j\)使\(\omega\)取最優值?SMO優化演算法能夠高效完成這個工作。SMO演算法非常的高效,只需要更多次數的迭代以達到收斂,而且每次迭代所需要的代價都非常小。

為了推出這個步驟,我們需要相對於\(\alpha_i\)、\(\alpha_j\)進行更新,假設取值是\(\alpha_1\)、\(\alpha_2\),即假設\(\alpha_1\)、\(\alpha_2\)不再是變數(可以由其他值推出),可以根據約束條件推導得到:

\[ \begin{split} \sum_{i=1}^m\alpha_iy^{(i)}&=0 \\ \alpha_1y_1+\alpha_2y_2&=-\sum_{i=3}^m\alpha_iy^{(i)} \end{split} \]

由於\(\alpha_3\)、\(\alpha_4\)、...、\(\alpha_m\)都是已知固定值,因此為了方便將等式右邊,可將等式右邊標記成\(\zeta\):

\[ \alpha_1y^{(1)}+\alpha_2^{(2)}=\zeta \]

還有一個約束條件:

\[ C \geq \alpha_i \geq 0,i=1,2,...,m \]

這個約束條件被稱作為“方形約束”,如果將\(\alpha_1\)、\(\alpha_2\)畫出來:

那麼\(\alpha_1\)、\(\alpha_2\)表示的值應該都在\([0,C]\)之間,也就是在方框裡面,這意味著:

\[ \alpha=\frac{\zeta-\alpha_2y^{(2)}}{y^{(1)}} \]

然後帶入到需要求解的式子中:

\[ W(\alpha_1,\alpha_2,...,\alpha_m)=W(\frac{\zeta-\alpha_2y^{(2)}}{y^{(1)}},\alpha_2,...,\alpha_m) \]

在前面我們認為\(\alpha_3\)、\(\alpha_4\)、...、\(\alpha_m\)都是已知固定值,只有\(\alpha_1\)、\(\alpha_2\)是未知需要求解的。那麼把\(W(\frac{\zeta-\alpha_2y^{(2)}}{y^{(1)}},\alpha_2,...,\alpha_m)\)展開後可以表示成\(a\alpha_2^2+b\alpha_2+c\)的形式,其中\(a\)、\(b\)、\(c\)是由\(\alpha_3\)、\(\alpha_4\)、...、\(\alpha_m\)表示出來,即\(W\)是一個二次函式。而其實對於所有的\(\alpha\),如果保持其他引數都固定的話,都可以表示成\(W\)關於某個\(\alpha\)的一元二次函式:

\[ \begin{split} W(\alpha_1,\alpha_2,...,\alpha_m)&=W(\frac{\zeta-\alpha_2y^{(2)}}{y^{(1)}},\alpha_2,...,\alpha_m) \\ &=a\alpha_2^2+b\alpha_2+c \end{split} \]

由於上面式子是一個標準的一元二次函式,所以很容易求解出最優值,從而可以得到\(\alpha_2\)的最優值,而這個最優值一定會在上圖中\(\alpha_1-\alpha_2=\zeta\)這條線上,且在“方形約束”中。按照這種方式解除\(\alpha_2\)後,之後根據\(\alpha_1\)和\(\alpha_2\)的關係求解出\(\alpha_1\),這樣子就求解出了相對於\(\alpha_1\)和\(\alpha_2\)關於\(W\),且滿足所有約束條件的最優值,該演算法的關鍵是對一個一元二次函式求最優解,這個求解非常簡單,這就使得SMO演算法的內嵌計算非常高效。

如何求解\(\alpha_2\)的值呢?只需要對式子進行求導\(a\alpha_2^2+b\alpha_2+c\),即對\(W\)進行求導,然而要保證\(\alpha_2\)即在方形約束內,也在\(\alpha_1-\alpha_2=\zeta\)這條線上,那麼就要保證\(H \geq \alpha_2 \geq L\),這裡使用\(\alpha_2^{new,unclipped}\)來表示求匯出來的\(\alpha_2\),然後最後\(\alpha_2^new\)的迭代更新方式如下所示:

\[ \alpha_2^new=\begin{cases} H, & \text {if $\alpha_2^{new,unclipped}>H$} \\ \alpha_2^{new,unclipped}, & \text{if $H \geq \alpha_2^{new,unclipped} \geq L$} \\ L, & \text{if $\alpha_2^{new,unclipped} < L$} \end{cases} \]

得到\(\alpha_2\)後,由此可以返回求解\(\alpha_1\)得到新值\(\alpha_1\),這裡就是SMO優化演算法的核心思想。根據SMO優化演算法的核心思想:

  1. 首先選擇兩個要改變的\(\alpha\)值\(\alpha_i\)、\(\alpha_j\)
  2. 其次保持除了\(\alpha_i\)、\(\alpha_j\)之外的所有引數固定
  3. 最後同時相對於這兩個引數使\(\omega\)取最優,且同時滿足所有約束條件

可以求解出所有的\(\alpha\),使得\(W\)取得最大值,即原問題將得到解決:

\[ max_{\alpha}W(\alpha)=\sum_{i=1}^m\alpha_i-\frac{1}{2}\sum^m_{i=1,j=1}\alpha_i\alpha_jy^{(i)}y^{(j)}<x^{(i)},x^{(j)}> \]

\[ C \geq \alpha_i \geq 0,i=1,2,...,m \]

\[ \sum_{i=1}^m\alpha_iy^{(i)}=0 \]

總結

svm的步驟總結如下:

  1. 先確定間隔器,這裡svm一般預設是幾何間隔
  2. 由間隔器確定間隔函式
  3. 從間隔函式檢視是否包含不等式約束形式
  4. 根據拉格朗日對偶步驟計算得到引數w、b
  5. 規則化不可分的引數,即在原對偶式子中加入離群點權重\(C\),問題轉換為\(max_{\alpha}W(\alpha)\)
  6. 利用SMO優化演算法求解\(W(\alpha)\)最優值,首先選擇兩個要改變的\(\alpha\)值\(\alpha_i\)、\(\alpha_j\)
  7. 其次保持除了\(\alpha_i\)、\(\alpha_j\)之外的所有引數固定
  8. 最後同時相對於這兩個引數使\(\omega\)取最優,且同時滿足所有約束條件,最後確定選取的這兩個\(\alpha_i\)、\(\alpha_j\)的值
  9. 重複步驟6-9直到所有引數\(\alpha\)求解完成

svm在神經網路出來之前一直是最優的演算法。相比於之前的演算法推導複雜一些,但是邏輯並不難,它不想邏輯迴歸那樣去擬合樣本點,而是根據幾何空間去尋找最優的分割超平面,為了判斷哪個超平面最好,引入幾個平面間隔最大化目標,從而求解出結果。

例項

有一份資料svm_data1,載入讀取:

import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
from sklearn import svm

# 載入data1
raw_data = loadmat('./svm_data1.mat')
# print(raw_data)

# 讀取data1的資料
data = pd.DataFrame(raw_data['X'], columns=['X1', 'X2'])
data['y'] = raw_data['y']

positive = data[data['y'].isin([1])]
negative = data[data['y'].isin([0])]
print(positive.shape)
print(negative.shape)

# 檢視data1的資料分佈
fig, ax = plt.subplots(figsize=(12, 8))
ax.scatter(positive['X1'], positive['X2'], s=50, marker='x', label='Positive')
ax.scatter(negative['X1'], negative['X2'], s=50, marker='o', label='Negative')
ax.legend()
plt.show()

資料分佈如下所示:

可以看到資料分在兩邊很好區分,用一般的分類器例如邏輯迴歸、樸素貝葉斯即可區分,這裡就用svm的線性核進行分類,設定離群點的權重\(C=1\),即不區分離群點:

svc = svm.LinearSVC(C=1, loss='hinge', max_iter=1000)

svc.fit(data[['X1', 'X2']], data['y'])
data1_score_1 = svc.score(data[['X1', 'X2']], data['y'])
print(data1_score_1)

得到的準確率為0.980392156863,分類的圖如下:

可以看到左上角有一個點原來是正樣本,但是被分類為藍色(負樣本),所以正樣本21個,負樣本30個,被誤分的概率剛好是\(\frac{1}{51}=‭0.01960784313‬\),所以準確率是\(1-‭0.01960784313‬=0.980392156863\),剛好對的上。現在這裡設定離群點的權重\(C=100\)用以區分離群點,得到的準確率為1.0,分類影象為:

再看第二份資料分佈圖如下:

這次就不能用線性核分類,需要用到RBF核分類:

# 做svm分類,使用RBF核
svc = svm.SVC(C=100, gamma=10, probability=True)
svc.fit(data[['X1', 'X2']], data['y'])
data['Probability'] = svc.predict_proba(data[['X1', 'X2']])[:, 0]

分類的結果圖如下所示:

結果得到的準確率只有0.769228287521,因此設定了網格調參:

# 簡單的網格調參
C_values = [0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100]
gamma_values = [0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100]

best_score = 0
best_params = {'C': None, 'gamma': None}

# 網格調參開始
for C in C_values:
    for gamma in gamma_values:

        # 做svm分類,使用RBF核
        svc = svm.SVC(C=C, gamma=gamma, probability=True)
        svc.fit(data[['X1', 'X2']], data['y'])

        # 交叉驗證
        data2_score = cross_validation.cross_val_score(svc, data[['X1', 'X2']], data['y'], scoring='accuracy', cv=3)
        print(data2_score.mean())

最後準確率提高到0.858437379017,調整到的最優引數為{'C': 10, 'gamma': 100}

資料和程式碼下載請關注公眾號【 機器學習和大資料探勘 】,後臺回覆【 機器學習 】即可獲