針對 「支持向量機」 將分爲如下幾個部分進行介紹:html
- 支持向量機 01 - 線性可分支持向量機和線性支持向量機
- 支持向量機 02 - SMO(序列最小化)(本篇)
- 支持向量機 03 - 非線性支持向量機
- 支持向量機 04 - SVR(支持向量迴歸)
本內容將介紹 SMO(序列最小化)算法,包含詳細公式推導以及 Python 代碼實現。
python
在學習本內容前,須要對支持向量機基本理論知識有必定了解,若是您還不瞭解,能夠參閱 支持向量機 01 - 線性可分支持向量機和線性支持向量機。
web
咱們知道,支持向量機的學習問題能夠形式化爲求解凸二次規劃問題。這樣的凸二次規劃問題具備全局最優解,而且有許多最優化算法能夠用於這一問題的求解。可是當訓練樣本容量很大時,這些算法每每變得很是低效,以至沒法使用。算法
SMO 算法是一種啓發式算法,其基本思路是:若是全部變量的解都知足此最優化問題的 KKT 條件,那麼這個最優化問題的解就獲得了。由於 KKT 條件是該最優化問題的充要條件。不然,選擇兩個變量,固定其餘變量,針對這兩個變量構建一個二次規劃問題。這個二次規劃問題關於這兩個變量的解應該更接近原始二次規劃問題的解,由於這會使得原始二次規劃問題的目標函數值變得更小。重要的是,這時子問題能夠經過解析方法求解,這樣就能夠大大提升整個算法的計算速度。子問題有兩個變量,一個是違反 KKT 條件最嚴重的那一個,另外一個由約束條件自動肯定。如此,SMO 算法將原問題不斷分解爲子問題並對子問題求解,進而達到求解原問題的目的。
緩存
1、SMO 算法描述
整個 SMO 算法包括兩個部分:求解兩個變量二次規劃的解析方法和選擇變量的啓發式方法。數據結構
1.1 兩個變量二次規劃的求解方法
1.1.1 問題轉化
在 支持向量機 01 - 線性可分支持向量機和線性支持向量機 中講解了 SVM 的基本原理,瞭解到 SMO 算法要解以下凸二次規劃的對偶問題app
αmax(i=1∑mαi−21i=1∑mj=1∑mαiαjy(i)y(j)x(i)⋅x(j))s.t.i=1∑mαiy(i)=0αi≥0,i=1,2,⋯,m(1)dom
添加一個負號,將最大值問題轉換成最小值問題機器學習
αmin(21i=1∑mj=1∑mαiαjy(i)y(j)x(i)⋅x(j)−i=1∑mαi)s.t.i=1∑mαiy(i)=0αi≥0,i=1,2,⋯,m(2)ide
1.1.2 轉換爲二元函數
假設選擇的兩個變量爲
α1 和
α2,其餘變量
αi(i=3,4,⋯,m) 是固定的,式(2)可寫爲
W(α1,α2)=21i=1∑mj=1∑mαiαjy(i)y(j)x(i)⋅x(j)−i=1∑mαi=21i=1∑m(j=1∑2αiαjy(i)y(j)x(i)⋅x(j)+j=3∑mαiαjy(i)y(j)x(i)⋅x(j))−αi−α2−i=3∑mαi=21i=1∑2(j=1∑2αiαjy(i)y(j)x(i)⋅x(j)+j=3∑mαiαjy(i)y(j)x(i)⋅x(j))+21i=3∑m(j=1∑2αiαjy(i)y(j)x(i)⋅x(j)+j=3∑mαiαjy(i)y(j)x(i)⋅x(j))−αi−α2−i=3∑mαi=21i=1∑2j=1∑2αiαjy(i)y(j)x(i)⋅x(j)+i=1∑2j=3∑mαiαjy(i)y(j)x(i)⋅x(j)+21i=3∑mj=3∑mαiαjy(i)y(j)x(i)⋅x(j)−α1−α2−i=3∑mαi=21α12x(1)⋅x(1)+21α22x(2)x(2)+y(1)y(2)α1α2x(1)⋅x(2)+y(1)α1j=3∑mαjy(j)x(1)⋅x(j)+y(2)α2j=3∑mαjy(j)x(2)⋅x(j)−α1−α2+21i=3∑mj=3∑mαiαjy(i)y(j)x(i)⋅x(j)−i=3∑mαi(3)
爲了描述方便,咱們定義以下符號:
f(x(i))=j=1∑mαjy(j)x(j)⋅x(i)+b=j=1∑mαjy(j)x(i)⋅x(j)+b(4)
vi=j=3∑mαjy(j)x(i)⋅x(j)=f(x(i))−j=1∑2αjy(j)x(i)⋅x(j)−b(5)
根據式(4)和(5),目標函數式(3)可寫爲
W(α1,α2)=21α12x(1)⋅x(1)+21α22x(2)⋅x(2)+y(1)y(2)α1α2x(1)⋅x(2)+y(1)α1v1+y(2)α2v2−α1−α2+constant(6)
其中
constant=21i=3∑mj=3∑mαiαjy(i)y(j)x(i)⋅x(j)−i=3∑mαi(7)
由於
constant 部分對
α1 和
α2 來講,屬於常數項,在求導的時候,直接變爲 0,因此咱們不須要關心這部分的內容。
1.1.3 轉換爲一元函數
由於有
∑i=1mαiy(i)=0,可得
α1y(1)+α2y(2)=−i=3∑mαiy(i)(8)
對變量
α1 和
α2 來講,
αi 和
y(i)(
i=3,4,⋯,m )可看做常數項,所以有
α1y(1)+α2y(2)=B(9)
其中,
B 爲一個常數。將等式(9)兩邊同時乘以
y(1),得
α1=By(1)−α2y(1)y(2)=(B−α2y(2))y(1)(10)
將其帶入式(6),獲得只含變量
α2 的目標函數
W(α2)=21(B−α2y(2))2x(1)⋅x(1)+21α22x(2)⋅x(2)+y(2)(B−α2y(2))α2x(1)⋅x(2)+(B−α2y(2))v1+y(2)α2v2−(B−α2y(2))y(1)−α2+constant(11)
1.1.4 求解
α2new,unclipped
將式(11)對
α2 求導,得
∂α2∂W(α2)=x(1)⋅x(1)α2+x(2)⋅x(2)α2−2x(1)⋅x(2)α2−Bx(1)⋅x(1)y(2)+Bx(1)⋅x(2)y(2)+y(1)y(2)−v1y(2)+v2y(2)−1(12)
令其爲
0,得
α2new=x(1)⋅x(1)+x(2)⋅x(2)−2x(1)⋅x(2)y(2)(y(2)−y(1)+B(x(1)⋅x(1)−x(1)⋅x(2))+v1−v2)(13)
定義以下符號
Ei=f(x(i))−y(i)(14)
η=x(1)⋅x(1)+x(2)⋅x(2)−2x(1)⋅x(2)(15)
其中
Ei 爲偏差項,
η 爲學習速率。
因爲已知
B=α1oldy(1)+α2oldy(2)(16)
vi=f(x(i))−j=1∑2αjy(j)x(i)⋅x(j)−b(17)
將式(16)和(17)帶入式(13), 將
α2new,unclipped 化簡後得
α2new,unclipped=α2old+ηy(2)(E1−E2)(18)
1.1.5 求解
α2new(對原始
α2new,unclipped 解進行修剪)
上面求出的
α2new,unclipped 解是沒有通過修剪的。咱們知道
αi 存在如下約束條件
⎩⎪⎨⎪⎧0<αi<Cα1y(1)+α2y(2)=B(19)
可知
α2new 的取值須要知足必定條件,假設爲
L≤α2new≤H(20)
圖-1 在二維平面上直觀地表達了式(19)中的約束條件,從而可知
- 當
y(1)̸=y(2) 時,存在
α2old−α1old=k,得
L=max(0,α2old−α1old),H=min(C,C+α2old−α1old)(21)
- 當
y(1)=y(2) 時,存在
α2old+α1old=k,得
L=max(0,α2old+α1old−C),H=min(C,α2old+α1old)(22)
圖-1
則通過修剪後,
α2new 的解爲
α2new=⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧H,α2new,unclipped,L,α2new,unclipped>HL≤α2new,unclipped≤Hα2new,unclipped<L(23)
1.1.6 求解
α1new
由於存在如下等式條件
α1oldy(1)+α2oldy(2)=α1newy(1)+α2newy(2)(24)
從而得出
α1new=α1old+y(1)y(2)(α2old−α2new)(25)
1.1.7 求解閾值
bnew
在求得
α1new 和
α2new 後,須要根據它們更新
b。
-
若是
0<α1new<C
由 KKT 條件可知(在 支持向量機 01 - 線性可分支持向量機和線性支持向量機 中 2.2 節有進行介紹),當
0<αi<C 時,這個樣本點是支持向量。所以,知足
y(1)(wTx(1)+b)=1(26)
上面公式兩邊同時乘以
y(1) 得
wTx(1)+b=y(1)(27)
將
w=∑i=1mαiy(i)x(i) 代入,得
i=1∑mαiy(i)x(i)⋅x(1)+b=y(1)(28)
由於咱們是根據
α1new 和
α2new 的值更新
b,整理可得:
b1new=y(1)−i=3∑mαiy(i)x(i)⋅x(1)−α1newy(1)x(1)⋅x(1)−α2newy(2)x(2)⋅x(1)(29)
其中有
y(1)−i=3∑mαiy(i)x(i)⋅x(1)=y(1)−i=1∑mαiy(i)x(i)⋅x(1)−bold+α1oldy(1)x(1)⋅x(1)+α2oldy(2)x(2)⋅x(1)+bold=y(1)−f(x(1))+α1oldy(1)x(1)⋅x(1)+α2oldy(2)x(2)⋅x(1)+bold=−E1+α1oldy(1)x(1)⋅x(1)+α2oldy(2)x(2)⋅x(1)+bold(30)
將式(30)代入式(29)得
b1new=bold−E1−y(1)(α1new−α1old)x(1)⋅x(1)−y(2)(α2new−α2old)x(2)x(1)(31)
-
若是
0<α2new<C
按照上面的步驟,一樣可求得
b2new
b2new=bold−E2−y(1)(α1new−α1old)x(1)⋅x(2)−y(2)(α2new−α2old)x(2)x(2)(32)
若是
α1new 和
α2new 同時知足條件
0<αinew<C,則
bnew=b1new=b2new;若是
α1new 和
α2new 是
0 或者
C,那麼
b1new 和
b2new 以及它們之間的數都是符合 KKT 條件的閾值,這是選擇它們的重點做爲
bnew。則
bnew 爲
bnew=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧b1new,b2new,2(b1new+b2new),0<α1new<C0<α2new<Cotherwise(33)
1.1.8 變量更新算法的步驟
根據上面講解的內容,咱們來整理一下變量更新算法的步驟:
E1=f(x(1))−y(1)=i=1∑mαiy(i)x(i)⋅x(1)+b−y(1)
E2=f(x(2))−y(2)=i=1∑mαiy(i)x(i)⋅x(2)+b−y(2)
- 步驟 2:計算
α2new 取值範圍
⎩⎪⎪⎨⎪⎪⎧if(y(1)̸=y(2))L=max(0,α2old−α1old),H=min(C,C+α2old−α1old)if(y(1)=y(2))L=max(0,α2old+α1old−C),H=min(C,α2old+α1old)
η=x(1)⋅x(1)+x(2)⋅x(2)−2x(1)x(2)
- 步驟 4:求解
α2new,unclipped
α2new,unclipped=α2old+ηy(2)(E1−E2)
- 步驟 5:求解
α2new
α2new=⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧H,α2new,unclipped,L,α2new,unclipped>HL≤α2new,unclipped≤Hα2new,unclipped<L
- 步驟 6:求解
α1new
α1new=α1old+y(1)y(2)(α2old−α2new)
- 步驟 7:求解
b1new 和
b2new
b1new=bold−E1−y(1)(α1new−α1old)x(1)⋅x(1)−y(2)(α2new−α2old)x(2)⋅x(1)
b2new=bold−E2−y(1)(α1new−α1old)x(1)⋅x(2)−y(2)(α2new−α2old)x(2)⋅x(2)
- 步驟 8:求解
bnew
bnew=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧b1new,b2new,2(b1new+b2new),0<α1new<C0<α2new<Cotherwise
前面進行介紹時,咱們有假設選擇兩個變量
α1 和
α2,可是在實際實現算法時,咱們應該如何選擇這兩個變量呢?下面就來介紹一下。
1.2 變量
α1 和
α2 的選擇方法
1.2.1 第 1 個變量的選擇
SMO 稱選擇第 1 個變量的過程爲外層循環。外層循環在訓練樣本中選取違反 KKT 條件最嚴重的樣本點,並將其對應的變量做爲第 1 個變量。具體地,檢驗訓練樣本點是否知足 KKT 條件,即(在 支持向量機 01 - 線性可分支持向量機和線性支持向量機 中 2.2 節有進行介紹)
αi=0⇔y(i)f(x(i))≥10<αi<C⇔y(i)f(x(i))=1αi=C⇔y(i)f(x(i))≤1(34)
從而可得知,如下幾種狀況將會不知足 KKT 條件
y(i)f(x(i))≥1andαi>0y(i)f(x(i))=1and(αi=0orαi=C)y(i)f(x(i))≤1andαi<C(35)
在檢驗過程當中,外層循環首先遍歷全部知足條件
0<αi<C 的樣本點,即在間隔邊界上的支持向量點,檢驗它們是否知足 KKT 條件。若是這些樣本點都知足 KKT 條件,那麼遍歷整個訓練集,檢驗它們是否知足 KKT 條件。
1.2.2 第 2 個變量的選擇
SMO 稱選擇第 2 個變量的過程爲內循環。假設在外層循環中已經找到第 1 個變量
α1,如今要在內層循環中找第 2 個變量
α2。第 2 個變量選擇的標準是但願能使
α2 有足夠大的變化。
由式(18)可知,
α2new 是依賴於
∣E1−E2∣ 的,爲了加快計算速度,一種簡單的作法是選擇
α2,使其對應的
∣E1−E2∣ 最大。爲了節省計算時間,能夠將全部
Ei 值保存在一個列表中。
在特殊狀況下,若是內層循環經過以上方法選擇的
α2 不能使目標函數有足夠的降低(等價於
∣α2new−α2∣ 很小),那麼能夠採用如下啓發式規則繼續選擇
α2。遍歷在間隔邊界上的支持向量點,依次將其對應的變量做爲
α2 試用,知道目標函數有足夠的降低。若仍然找不到合適的
α2,那麼遍歷訓練數據集;若仍是找不到合適的
α2,則放棄第 1 個
α2,再經過外層循環尋求另外的
α1。
由 Osuna 定理可知,只需選取的
α1 和
α2 中有一個不知足 KKT 條件,目標函數就會在迭代後減少。
1.3 SMO 算法的步驟
咱們總結一下 SMO 算法的步驟:
- 步驟 1:初始化
α 和
b,好比初始化
α 爲零向量,
b 爲 0。(注意這裏的
α 是一個列向量)
- 步驟 2:選取優化變量
α1 和
α2,而後更新
α1、
α2 和
b。
- 步驟 3:若是知足中止條件(即前面提到的「全部變量的解都知足此最優化問題的 KKT 條件」)則結束;不然,跳到步驟 2。
2、Python 代碼實現
如下是 Python 3 的代碼實現:
import numpy as np
import random
import matplotlib.pyplot as plt
class OptStruct:
""" 數據結構,存儲須要操做的數據 """
def __init__(self, input_mat, label_mat, c, toler):
""" :param input_mat: 樣本特徵值 :param label_mat: 樣本標籤值 :param c: 懲罰因子 :param toler: 容錯率 """
self.x = input_mat
self.label_mat = label_mat
self.c = c
self.toler = toler
self.m = np.shape(input_mat)[0]
self.alphas = np.mat(np.zeros((self.m, 1)))
self.b = 0
self.e_cache = np.mat(np.zeros((self.m, 2)))
class SVM:
def __init__(self):
self.alphas = None
self.b = None
self.w = None
pass
@staticmethod
def calc_ek(os, k):
""" 計算偏差 E 值 :param os: :param k: :return: """
fxk = float(np.multiply(os.alphas, os.label_mat).T *\
(os.x * os.x[k, :].T)) + os.b
ek = fxk - float(os.label_mat[k])
return ek
@staticmethod
def select_j_rand(i, m):
""" 隨機選擇第二個 alpha :param i: 第一個 alpha 對應的索引值 :param m: alpha 參數的個數(即全部訓練樣本的總數) :return: 第二個 alpha 對應的索引值 """
j = i
while j == i:
j = int(random.uniform(0, m))
return j
@staticmethod
def select_j(i, os, ei):
""" 選擇第二個 alpha :param i: 第一個 alpha 的索引值 :param os: 數據結構 :param ei: 第一個 alpha 對應樣本點的偏差 :return: 返回第二個 alpha 的索引值 和 對應的偏差值 """
max_k = -1
max_delta_e = 0
ej = 0
os.e_cache[i] = [1, ei]
valid_e_cache_list = np.nonzero(os.e_cache[:, 0].A)[0]
if (len(valid_e_cache_list)) > 1:
for k in valid_e_cache_list:
if k == i:
continue
ek = SVM.calc_ek(os, k)
delta_e = abs(ei - ek)
if delta_e > max_delta_e:
max_k = k
max_delta_e = delta_e
ej = ek
return max_k, ej
else:
j = SVM.select_j_rand(i, os.m)
ej = SVM.calc_ek(os, j)
return j, ej
@staticmethod
def update_ek(os, k):
""" 更新偏差緩存 :param os: 數據結構 :param k: 須要更新偏差項對應的索引值 """
ek = SVM.calc_ek(os, k)
os.e_cache[k] = [1, ek]
@staticmethod
def clip_alpha(alpha_j, high, low):
""" 修剪 alpha_j :param alpha_j: alpha_j 的未修剪的值 :param high: alpha_j 的上限 :param low: alpha_j 的下限 :return: alpha_j 修剪後的值 """
if alpha_j > high:
alpha_j = high
elif alpha_j < low:
alpha_j = low
return alpha_j
def inner_l(self, i, os):
""" 選擇第二個 alpha,並更新 alpha 和 b :param i: 第一個 alpha 對應樣本點的索引值 :param os: 數據結構 :return: 是否有更新 alpha 對 """
ei = SVM.calc_ek(os, i)
if ((os.label_mat[i] * ei < -os.toler) and (os.alphas[i] < os.c)) or \
((os.label_mat[i] * ei > os.toler) and (os.alphas[i] > 0)):
j, ej = SVM.select_j(i, os, ei)
alpha_iold = os.alphas[i].copy()
alpha_jold = os.alphas[j].copy()
if os.label_mat[i] != os.label_mat[j]:
low = max(0, os.alphas[j] - os.alphas[i])
high = min(os.c, os.c + os.alphas[j] - os.alphas[i])
else:
low = max(0, os.alphas[j] + os.alphas[i] - os.c)
high = min(os.c, os.alphas[j] + os.alphas[i])
if low == high:
print('low == high')
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.label_mat[j] * (ei - ej) / eta
os.alphas[j] = SVM.clip_alpha(os.alphas[j], high, low)
SVM.update_ek(os, j)
if abs(os.alphas[j] - alpha_jold) < 0.00001:
print('j not moving enough')
return 0
os.alphas[i] += os.label_mat[j] * os.label_mat[i] * \
(alpha_jold - os.alphas[j])
SVM.update_ek(os, i)
b1 = os.b - ei - os.label_mat[i] * (os.alphas[i] - alpha_iold) * \
os.x[i, :] * os.x[i, :].T