转载请注明出处: 嚼烂SVM-SMO算法并用python进行实现

本文面向的读者为掌握SVM基础前置知识如读过《统计学习方法》，并希望对SMO细节有更深入了解的人羣。因为笔者想要实现一个简易的SVM,为了搞懂这部分费了不少工夫，所以写下这一篇嚼过的文章，目的是让读者跟着顺序阅读一定能弄懂。自己的行文习惯是简单的地方不能说太复杂，复杂的地方一定会说清楚

一, 推导至SMO要解决的SVM对偶形式

对于含soft margin的支持向量机, 其primitive problem:

minw,b,ξ12∥w∥2+C∑i=1Nξis.t.yi(w⋅xi+b)≥1−ξi,i=1,2,⋯,Nξi≥0,i=1,2,⋯,N(1)(2)(3) (1) min w , b , ξ 1 2 ‖ w ‖ 2 + C ∑ i = 1 N ξ i (2) s . t . y i ( w ⋅ x i + b ) ≥ 1 − ξ i , i = 1 , 2 , ⋯ , N (3) ξ i ≥ 0 , i = 1 , 2 , ⋯ , N

求解此原始问题：

1, 构建Lagrange Function

引入拉格朗日乘子 αi≥0,μi≥0,i=1,2,⋯,N α i ≥ 0 , μ i ≥ 0 , i = 1 , 2 , ⋯ , N 构建Lagrange Function:

L(w,b,ξ,α,μ)=12∥w∥2+C∑i=1Nξi+∑i=1Nαi(−yi(w⋅xi+b)+1−ξi)−∑i=1Nμiξi(4) (4) L ( w , b , ξ , α , μ ) = 1 2 ‖ w ‖ 2 + C ∑ i = 1 N ξ i + ∑ i = 1 N α i ( − y i ( w ⋅ x i + b ) + 1 − ξ i ) − ∑ i = 1 N μ i ξ i

其中，

α=(α1,α2,⋯,αN)T α = ( α 1 , α 2 , ⋯ , α N ) T 以及

μ=(μ1,μ2,⋯,μN)T μ = ( μ 1 , μ 2 , ⋯ , μ N ) T 为
lagrange multiplier , 它们的每个分量都是非负数

2,转化为Lagrange Dual Problem

具体原理与KKT条件推导可看我上一篇博文:SVM之拉格朗日对偶问题与KKT条件推导
现在dual problem:

maxα,μminw,b,ξL(w,b,ξ,α,μ)(5) (5) max α , μ min w , b , ξ L ( w , b , ξ , α , μ )

3,先求内层的min

由于把 α,μ α , μ 都看作常量， 要求最小值就直接求偏导:

∇wL(w,b,ξ,α,μ)=w−∑i=1Nαiyixi=0∇bL(w,b,ξ,α,μ)=−∑i=1Nαiyi=0∇ξiL(w,b,ξ,α,μ)=C−αi−μi=0(6)(7)(8) (6) ∇ w L ( w , b , ξ , α , μ ) = w − ∑ i = 1 N α i y i x i = 0 (7) ∇ b L ( w , b , ξ , α , μ ) = − ∑ i = 1 N α i y i = 0 (8) ∇ ξ i L ( w , b , ξ , α , μ ) = C − α i − μ i = 0

得:

w=∑i=1Nαiyixi∑i=1Nαiyi=0C−αi−μi=0(9)(10)(11) (9) w = ∑ i = 1 N α i y i x i (10) ∑ i = 1 N α i y i = 0 (11) C − α i − μ i = 0

把(9)-(11) 代入(4)可得(这里引入了kernel function):

minw,b,ξL(w,b,ξ,α,μ)=−12∑i=1N∑i=1NαiαjyiyjK(xi,xj)+∑i=1Nαi(12) (12) min w , b , ξ L ( w , b , ξ , α , μ ) = − 1 2 ∑ i = 1 N ∑ i = 1 N α i α j y i y j K ( x i , x j ) + ∑ i = 1 N α i

4, 求解外层的max

对式(12)求解max:

maxα,μ−12s.t.∑i=1N∑i=1NαiαjyiyjK(xi,xj)+∑i=1Nαi∑i=1Nαiyi=0C−αi−μi=0αi≥0μi≥0(13)(14)(15)(16)(17) (13) max α , μ − 1 2 ∑ i = 1 N ∑ i = 1 N α i α j y i y j K ( x i , x j ) + ∑ i = 1 N α i (14) s . t . ∑ i = 1 N α i y i = 0 (15) C − α i − μ i = 0 (16) α i ≥ 0 (17) μ i ≥ 0

式(15)-(17) 可简化为:

0≤αi≤C(18) (18) 0 ≤ α i ≤ C

5,最终形式:一个凸二次规划的对偶问题

minαs.t.12∑i=1N∑j=1NαiαjyiyjK(xi,xj)−∑i=1Nαi∑i=1Nαiyi=00≤αi≤C,i=1,2,⋯,N(19)(20)(21) (19) min α 1 2 ∑ i = 1 N ∑ j = 1 N α i α j y i y j K ( x i , x j ) − ∑ i = 1 N α i (20) s . t . ∑ i = 1 N α i y i = 0 (21) 0 ≤ α i ≤ C , i = 1 , 2 , ⋯ , N

收敛的值需满足的
KKT condition(求解SMO时有用):

∇wL(w∗,b∗,ξ∗,α∗,μ∗)∇bL(w∗,b∗,ξ∗,α∗,μ∗)∇ξiL(w∗,b∗,ξ∗,α∗,μ∗)α∗i1−ξ∗i+yi(α∗i(1−ξ∗i+yiμ∗iξ∗iμ∗iξ∗i=w∗−∑i=1Nα∗iyixi=0=−∑i=1Nα∗iyi=0=C−α∗i−μ∗i=0≥0w∗⋅x∗+b∗)≤0(w∗⋅x∗+b∗))=0≥0≤0=0(22)(23)(24)(25)(26)(27)(28)(29)(30) (22) ∇ w L ( w ∗ , b ∗ , ξ ∗ , α ∗ , μ ∗ ) = w ∗ − ∑ i = 1 N α i ∗ y i x i = 0 (23) ∇ b L ( w ∗ , b ∗ , ξ ∗ , α ∗ , μ ∗ ) = − ∑ i = 1 N α i ∗ y i = 0 (24) ∇ ξ i L ( w ∗ , b ∗ , ξ ∗ , α ∗ , μ ∗ ) = C − α i ∗ − μ i ∗ = 0 (25) α i ∗ ≥ 0 (26) 1 − ξ i ∗ + y i ( w ∗ ⋅ x ∗ + b ∗ ) ≤ 0 (27) α i ∗ ( 1 − ξ i ∗ + y i ( w ∗ ⋅ x ∗ + b ∗ ) ) = 0 (28) μ i ∗ ≥ 0 (29) ξ i ∗ ≤ 0 (30) μ i ∗ ξ i ∗ = 0

一共3个偏导为0 + 2*3个乘子与不等式约束间满足的条件 = 9个KKT 条件

二,切入正题:SMO算法完整推导

SMO算法，首先选择 α α 的两个分量来求解式(19)的子问题，思路是不断迭代求解原问题的子问题来逼近原始问题的解, 具体怎么选取这两个分量待会儿再说

1, 选取两个 α α 的分量,求解式(19)子问题，目的是获取对这两个分量进行更新的方法

minα1,α2W(α1,α2)=12α21K11+12α22K22+α1α2y1y2K12+α1y1∑i=3NαiyiK1i+α2y2∑i=3NαiyiK2i−α1−α2s.t.α1y1+α2y2=−∑i=3Nyiαi=ς0≤αi≤C,i=1,2(31)(32)(33) (31) min α 1 , α 2 W ( α 1 , α 2 ) = 1 2 α 1 2 K 11 + 1 2 α 2 2 K 22 + α 1 α 2 y 1 y 2 K 12 + α 1 y 1 ∑ i = 3 N α i y i K 1 i + α 2 y 2 ∑ i = 3 N α i y i K 2 i − α 1 − α 2 (32) s . t . α 1 y 1 + α 2 y 2 = − ∑ i = 3 N y i α i = ς (33) 0 ≤ α i ≤ C , i = 1 , 2

1.1 换元

有两个变量，那先通过式(32)的变形式 α1=(ς−α2y2)y1 α 1 = ( ς − α 2 y 2 ) y 1 代入式(31)换为只包含 α2 α 2 的式子:

minα2W(α2)=12(ς−α2y2)2K11+12α22K22+(ς−α2y2)α2y2K12+(ς−α2y2)∑i=3NαiyiK1i+α2y2∑i=3NαiyiK2i−(ς−α2y2)y1−α2(34) (34) min α 2 W ( α 2 ) = 1 2 ( ς − α 2 y 2 ) 2 K 11 + 1 2 α 2 2 K 22 + ( ς − α 2 y 2 ) α 2 y 2 K 12 + ( ς − α 2 y 2 ) ∑ i = 3 N α i y i K 1 i + α 2 y 2 ∑ i = 3 N α i y i K 2 i − ( ς − α 2 y 2 ) y 1 − α 2

1.2 求极值点

明显地，对 α2 α 2 求偏导并令其为0:

∂Wα2=−ςy2K11+K11α2+K22α2+ςy2K12−2K12α2−y2∑i=3NαiyiK1i+y2∑i=3NαiyiK2i+y1y2−1=(K11+K22−2K12)α2+y2(−ςK11+ςK12−∑i=3NαiyiK1i+∑i=3NαiyiK2i+y1−y2)(35)(36) (35) ∂ W α 2 = − ς y 2 K 11 + K 11 α 2 + K 22 α 2 + ς y 2 K 12 − 2 K 12 α 2 − y 2 ∑ i = 3 N α i y i K 1 i + y 2 ∑ i = 3 N α i y i K 2 i + y 1 y 2 − 1 (36) = ( K 11 + K 22 − 2 K 12 ) α 2 + y 2 ( − ς K 11 + ς K 12 − ∑ i = 3 N α i y i K 1 i + ∑ i = 3 N α i y i K 2 i + y 1 − y 2 )

这里要设置一些方便的符号:

1, 模型对

x x 的预测值

g(x)=∑i=1NαiyiK(xi,x)+b(37) (37) g ( x ) = ∑ i = 1 N α i y i K ( x i , x ) + b

2, 预测值减去真实值

Ei=g(xi)−yi=(∑j=1NαjyjK(xj,xi)+b)−yi,i=1,2(38) (38) E i = g ( x i ) − y i = ( ∑ j = 1 N α j y j K ( x j , x i ) + b ) − y i , i = 1 , 2

3,式(36)中比较难处理的那一坨

vi=∑j=3NαjyjKij=g(xi)−∑j=12αjyjKij−b,i=1,2=Ei+yi−∑j=12αjyjKij−b,i=1,2(39)(40) (39) v i = ∑ j = 3 N α j y j K i j = g ( x i ) − ∑ j = 1 2 α j y j K i j − b , i = 1 , 2 (40) = E i + y i − ∑ j = 1 2 α j y j K i j − b , i = 1 , 2

还要注意

α1y1+α2y2=∑2i=1αiyi=ς α 1 y 1 + α 2 y 2 = ∑ i = 1 2 α i y i = ς

令式(36)为0, 并代入上面的符号:

(K11+K22−2K12)αnew,unc2=y2(ςK11−ςK12+∑i=3NαiyiK1i−∑i=3NαiyiK2i−y1+y2)=y2(∑i=12αiyiK11−∑i=12αiyiK12+v1−v2−y1+y2)=y2(∑i=12αiyiK11−∑i=12αiyiK12+E1+y1−∑i=12αiyiK1i−b−E2−y2+∑i=12αiyiK2i+b−y1+y2)=y2(E1−E2+α2y2K11−2α2y2K12+α2y2K22)=(K11−2K12+K22)α2+y2(E1−E2)(41)(42)(43)(44)(45) (41) ( K 11 + K 22 − 2 K 12 ) α 2 n e w , u n c = y 2 ( ς K 11 − ς K 12 + ∑ i = 3 N α i y i K 1 i − ∑ i = 3 N α i y i K 2 i − y 1 + y 2 ) (42) = y 2 ( ∑ i = 1 2 α i y i K 11 − ∑ i = 1 2 α i y i K 12 + v 1 − v 2 − y 1 + y 2 ) (43) = y 2 ( ∑ i = 1 2 α i y i K 11 − ∑ i = 1 2 α i y i K 12 + E 1 + y 1 − ∑ i = 1 2 α i y i K 1 i − b − E 2 − y 2 + ∑ i = 1 2 α i y i K 2 i + b − y 1 + y 2 ) (44) = y 2 ( E 1 − E 2 + α 2 y 2 K 11 − 2 α 2 y 2 K 12 + α 2 y 2 K 22 ) (45) = ( K 11 − 2 K 12 + K 22 ) α 2 + y 2 ( E 1 − E 2 )

1.3 获取 αnew,unc2 α 2 n e w , u n c 的迭代方式

由式(45)可得:

αnew,unc2=αold2+y2(E1−E2)K11−2K12+K22(46) (46) α 2 n e w , u n c = α 2 o l d + y 2 ( E 1 − E 2 ) K 11 − 2 K 12 + K 22

1.4 根据 α2 α 2 的定义域裁剪得到 αnew2 α 2 n e w

上式左边有unc, 意味着没有进行cut,现在我们讨论下 α2 α 2 的取值范围:
关于 α1,α2 α 1 , α 2 的约束条件一共就式(32),(33)两个式子
那么我们根据 y1,y2 y 1 , y 2 进行分类讨论:

• y1=y2 y 1 = y 2
  根据式(32),设 α1+α2=k α 1 + α 2 = k
  根据式(33),可得:
  {00≤α2≤C≤k−α2≤C⇒{0≤α2≤Ck−C≤α2≤k⇒{0≤α2≤Cαold1+αold2−C≤α2≤αold1+αold2 { 0 ≤ α 2 ≤ C 0 ≤ k − α 2 ≤ C ⇒ { 0 ≤ α 2 ≤ C k − C ≤ α 2 ≤ k ⇒ { 0 ≤ α 2 ≤ C α 1 o l d + α 2 o l d − C ≤ α 2 ≤ α 1 o l d + α 2 o l d
  设 α2 α 2 上界为 H H ,下界为 L L , 则有:
  L=max(0,αold1+αold2−C)H=min(C,αold1+αold2)(47)(48) (47) L = max ( 0 , α 1 o l d + α 2 o l d − C ) (48) H = min ( C , α 1 o l d + α 2 o l d )
• y1≠y2 y 1 ≠ y 2
  根据式(32),设 α1−α2=k α 1 − α 2 = k
  根据式(33),可得:
  {00≤α2≤C≤α2+k≤C⇒{0−k≤α2≤C≤α2≤C−k⇒{0αold2−αold1≤α2≤C≤α2≤C+αold2−αold1 { 0 ≤ α 2 ≤ C 0 ≤ α 2 + k ≤ C ⇒ { 0 ≤ α 2 ≤ C − k ≤ α 2 ≤ C − k ⇒ { 0 ≤ α 2 ≤ C α 2 o l d − α 1 o l d ≤ α 2 ≤ C + α 2 o l d − α 1 o l d
  可得新的上下界:
  L=max(0,αold2−αold1)H=min(C,C+αold2−αold1)(49)(50) (49) L = max ( 0 , α 2 o l d − α 1 o l d ) (50) H = min ( C , C + α 2 o l d − α 1 o l d )

裁剪方法:

αnew2=⎧⎩⎨⎪⎪H,αnew,unc2,L,αnew,unc2>HL≤αnew,unc2≤Cαnew,unc2<L α 2 n e w = { H , α 2 n e w , u n c > H α 2 n e w , u n c , L ≤ α 2 n e w , u n c ≤ C L , α 2 n e w , u n c < L

1.5 根据约束条件得到 αnew1 α 1 n e w

根据式(32)可得:

αold1y1+αold2y2=αnew1y1+αnew2y2(51) (51) α 1 o l d y 1 + α 2 o l d y 2 = α 1 n e w y 1 + α 2 n e w y 2

根据上式，有:

αnew1==(αold1y1+αold2y2−αnew2y2)y1αold1+(αold2−αnew2)y1y2(52)(53) (52) α 1 n e w = ( α 1 o l d y 1 + α 2 o l d y 2 − α 2 n e w y 2 ) y 1 (53) = α 1 o l d + ( α 2 o l d − α 2 n e w ) y 1 y 2

2, 通过对这两个 α α 分量的更新,获取其它变量的更新

2.1 计算阈值b

原理,通过支持向量,即正好在间隔边界的点来进行计算(此时 0<αi<C,yig(xi)=1 0 < α i < C , y i g ( x i ) = 1 ):

b=yi−∑j=1NαjyjKij(54) (54) b = y i − ∑ j = 1 N α j y j K i j

如果

α1 α 1 满足此条件,则

bnew1=y1−∑i=3NαiyiK1i−αnew1y1K11−αnew2y2K12(55) (55) b 1 n e w = y 1 − ∑ i = 3 N α i y i K 1 i − α 1 n e w y 1 K 11 − α 2 n e w y 2 K 12

上面的公式有一部分可以用

E1 E 1 进行替换:

E1=g(x1)−y1=∑i=3NαiyiK1i+αold1y1K11+αold2y2K21+bold−y1(56) (56) E 1 = g ( x 1 ) − y 1 = ∑ i = 3 N α i y i K 1 i + α 1 o l d y 1 K 11 + α 2 o l d y 2 K 21 + b o l d − y 1

结合式(55)与式(56),将

E1 E 1 引入式(55)可得:

bnew1=−E1+y1K11(αold1−αnew1)+y2K12(αold2−αnew2)+bold(57) (57) b 1 n e w = − E 1 + y 1 K 11 ( α 1 o l d − α 1 n e w ) + y 2 K 12 ( α 2 o l d − α 2 n e w ) + b o l d