嚼烂SVM-SMO算法并用python进行实现

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

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

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

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

(1)minw,b,ξ12w2+Ci=1Nξi(2)s.t.yi(wxi+b)1ξi,i=1,2,,N(3)ξi0,i=1,2,,N

求解此原始问题:

1, 构建Lagrange Function

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

(4)L(w,b,ξ,α,μ)=12w2+Ci=1Nξi+i=1Nαi(yi(wxi+b)+1ξi)i=1Nμiξi

其中,α=(α1,α2,,αN)T以及μ=(μ1,μ2,,μN)Tlagrange multiplier , 它们的每个分量都是非负数

2,转化为Lagrange Dual Problem

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

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

3,先求内层的min

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

(6)wL(w,b,ξ,α,μ)=wi=1Nαiyixi=0(7)bL(w,b,ξ,α,μ)=i=1Nαiyi=0(8)ξiL(w,b,ξ,α,μ)=Cαiμi=0

得:
(9)w=i=1Nαiyixi(10)i=1Nαiyi=0(11)Cαiμi=0

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

(12)minw,b,ξL(w,b,ξ,α,μ)=12i=1Ni=1NαiαjyiyjK(xi,xj)+i=1Nαi

4, 求解外层的max

对式(12)求解max:

(13)maxα,μ12i=1Ni=1NαiαjyiyjK(xi,xj)+i=1Nαi(14)s.t.i=1Nαiyi=0(15)Cαiμi=0(16)αi0(17)μi0

式(15)-(17) 可简化为:
(18)0αiC

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

(19)minα12i=1Nj=1NαiαjyiyjK(xi,xj)i=1Nαi(20)s.t.i=1Nαiyi=0(21)0αiC,i=1,2,,N

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

(22)wL(w,b,ξ,α,μ)=wi=1Nαiyixi=0(23)bL(w,b,ξ,α,μ)=i=1Nαiyi=0(24)ξiL(w,b,ξ,α,μ)=Cαiμi=0(25)αi0(26)1ξi+yi(wx+b)0(27)αi(1ξi+yi(wx+b))=0(28)μi0(29)ξi0(30)μiξi=0

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

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

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

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

(31)minα1,α2W(α1,α2)=12α12K11+12α22K22+α1α2y1y2K12+α1y1i=3NαiyiK1i+α2y2i=3NαiyiK2iα1α2(32)s.t.α1y1+α2y2=i=3Nyiαi=ς(33)0αiC,i=1,2

1.1 换元

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

(34)minα2W(α2)=12(ςα2y2)2K11+12α22K22+(ςα2y2)α2y2K12+(ςα2y2)i=3NαiyiK1i+α2y2i=3NαiyiK2i(ςα2y2)y1α2

1.2 求极值点

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

(35)Wα2=ςy2K11+K11α2+K22α2+ςy2K122K12α2y2i=3NαiyiK1i+y2i=3NαiyiK2i+y1y21(36)=(K11+K222K12)α2+y2(ςK11+ςK12i=3NαiyiK1i+i=3NαiyiK2i+y1y2)

这里要设置一些方便的符号:
1, 模型对x的预测值
(37)g(x)=i=1NαiyiK(xi,x)+b

2, 预测值减去真实值
(38)Ei=g(xi)yi=(j=1NαjyjK(xj,xi)+b)yi,i=1,2

3,式(36)中比较难处理的那一坨
(39)vi=j=3NαjyjKij=g(xi)j=12αjyjKijb,i=1,2(40)=Ei+yij=12αjyjKijb,i=1,2

还要注意α1y1+α2y2=i=12αiyi=ς
令式(36)为0, 并代入上面的符号:
(41)(K11+K222K12)α2new,unc=y2(ςK11ςK12+i=3NαiyiK1ii=3NαiyiK2iy1+y2)(42)=y2(i=12αiyiK11i=12αiyiK12+v1v2y1+y2)(43)=y2(i=12αiyiK11i=12αiyiK12+E1+y1i=12αiyiK1ibE2y2+i=12αiyiK2i+by1+y2)(44)=y2(E1E2+α2y2K112α2y2K12+α2y2K22)(45)=(K112K12+K22)α2+y2(E1E2)

1.3 获取α2new,unc的迭代方式

由式(45)可得:

(46)α2new,unc=α2old+y2(E1E2)K112K12+K22

1.4 根据α2的定义域裁剪得到α2new

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

  • y1=y2
    根据式(32),设α1+α2=k
    根据式(33),可得:
    {0α2C0kα2C{0α2CkCα2k{0α2Cα1old+α2oldCα2α1old+α2old

    α2上界为H,下界为L, 则有:
    (47)L=max(0,α1old+α2oldC)(48)H=min(C,α1old+α2old)
  • y1y2
    根据式(32),设α1α2=k
    根据式(33),可得:
    {0α2C0α2+kC{0α2Ckα2Ck{0α2Cα2oldα1oldα2C+α2oldα1old

    可得新的上下界:
    (49)L=max(0,α2oldα1old)(50)H=min(C,C+α2oldα1old)

裁剪方法:

α2new={H,α2new,unc>Hα2new,unc,Lα2new,uncCL,α2new,unc<L

1.5 根据约束条件得到α1new

根据式(32)可得:

(51)α1oldy1+α2oldy2=α1newy1+α2newy2

根据上式,有:
(52)α1new=(α1oldy1+α2oldy2α2newy2)y1(53)=α1old+(α2oldα2new)y1y2

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

2.1 计算阈值b

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

(54)b=yij=1NαjyjKij

如果α1满足此条件,则
(55)b1new=y1i=3NαiyiK1iα1newy1K11α2newy2K12

上面的公式有一部分可以用E1进行替换:
(56)E1=g(x1)y1=i=3NαiyiK1i+α1oldy1K11+α2oldy2K21+boldy1

结合式(55)与式(56),将E1引入式(55)可得:
(57)b1new=E1+y1K11(α1oldα1new)+y2K12(α2oldα2new)+bold

每次计算的时候,存下Ei可以极大的方便计算
同理,如果0<α2new<C, 则:
(58)b2new=E2+y1K12(α1oldα1new)+y2K22(α2oldα2new)+bold

下面讨论bnew的最终取值:

  • 0<α1<C, 0<α2<C:
    bnew=b1new=b2new (此时x1x2都在间隔边界上)

  • 若只有一个0<αi<C,i{1,2}
    bnew=binew

  • α1,α2{0,C}
    bnew=12(b1new+b2new), (若αi=0说明xi不是支持向量,yig(xi)1, xi在正确分类的间隔一侧, αi=C 说明 yig(xi)1, 这些都可以从式(22)-式(30)的KKT条件推出,下面还会推导)

2.2 更新Ei,方便下一次的 b 计算

(59)Ei=sαjyjKijyi

其中s是所有>0的αj, 即所有支持向量

3 α选取策略

3.1 通过满足KKT条件与否选择α1

因为收敛后的最优解是满足KKT条件的,所以第一次选择最不满足KKT条件的α:
从式(22)-(30)可得:

  • αi=0
    (1), 根据Cαiui=0 可得 ui=C>0
    (2), 根据uiξi=0 可得 ξi=0
    (3), 根据yi(wxi+b)1ξi 可得 yi(wxi+b)1
    (4), 综上,(60)αi=0yig(xi)1

  • 0<αi<C
    (1), 根据Cαiui=0 可得 ui>0
    (2), 根据uiξi=0 可得 ξi=0
    (3), 根据αi(yi(wxi+b)1+ξ)=0 及上面一条可得 yi(wxi+b)1=0
    (4), 综上, (61)0<αi<Cyig(xi)=1

  • αi=C
    (1). 根据Cαiui=0 可得 ui=0
    (2). 根据ui=0uiξi=0, ξi0可得ξi0
    (3). 根据αi(yi(wxi+b)1+ξi)=0αi=C>0可得