对SVM的推导和编码实践（二）SMO算法的推导

目标函数和约束条件

minα12∑i,j=1Nαiαjyiyj<xi,xj>−∑i=1Nαis.t.,∑i=1Nαiyi=00≤αi≤C(13)(13)minα12∑i,j=1Nαiαjyiyj<xi,xj>−∑i=1Nαis.t.,∑i=1Nαiyi=00≤αi≤C

SMO简介

SMO算法的目标是求出一系列alpha和b，一旦求出了这些alpha和b，就很容易计算出权重向量w并得到分隔超平面。

SMO算法的工作原理是：每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha，那么就增大其中一个同时减小另一个。这里所谓的“合适”就是指两个alpha必须要符合一定的条件，条件之一就是这两个alpha必须要在间隔边界之外，而其第二个条件则是这两个alpha还没有进行过区间化处理或者不在边界上。

Platt SMO算法中的外循环确定要优化的最佳alpha对。而简化版却会跳过这一部分，首先在数据集上遍历每一个alpha，然后在剩下的alpha集合中随机选择另一个alpha，从而构建alpha对。这里有一点相当重要，就是我们要同时改变两个alpha。之所以这样做是因为我们有一个约束条件：

∑Ni=1αiyi=0(14)(14)∑i=1Nαiyi=0

由于改变一个alpha可能会导致该约束条件失效，因此我们总是同时改变两个alpha。

SMO是一个可以快速解决SVM QP问题而不使用矩阵存储空间和数值优化步的简单算法。SMO使用Qsuna的理论分解QP问题以确保收敛。

SMO在每一步选择尽可能小的优化问题。对标准的SVM QP问题，最小的优化问题涉及到两个拉格朗日乘数，因为拉格朗日乘数必须遵循一个线性等式约束。在每一步SMO选择两个乘数一起优化，寻找最优值，更新SVM以体现这些新的最优值。

SMO的优势体现于解那两个乘数的最优值的时候可以直接计算解析解而不是通过数值优化。此外，SMO不需要额外的空间存储矩阵，因此非常大规模的SVM训练问题也可以装进一台普通的个人电脑的内存里。因为没有涉及到矩阵算法，SMO算法不受数值精度问题的影响。

SMO由两部分组成：

1、解那两个拉格朗日乘数的解析解

2、如何选择那两个拉格朗日乘数进行优化的启发式算法

SMO的基本思路

确保自己理解思路，觉得《统计学习方法》里面这段话很重要，于是抄录如下：

SMO算法是一种启发式算法，其基本思路是：如果所有变量的解都满足此优化问题的KKT条件，那么这个最优化问题的解就得到了。因为KKT条件是该最优化问题的充分必要条件。否则，选择两个变量，固定其他变量，针对这两个变量构建一个二次规划问题。这个二次规划问题关于这两个变量的解应该更接近原始二次规划问题的解，因为这会使得原始二次规划问题的目标函数的值变得更小。重要的是，这时子问题可以通过解析方法求解（注：求导什么的）。这样就可以大大提高整个算法的计算速度。子问题有两个变量，一个是违反kkt条件最严重的那一个，另一个由约束条件自动确定。如此，SMO算法将原问题不断分解为子问题并对子问题求解，进而达到求原问题的目的。

两个变量的二次规划求解方法

得到alpha _1,alpha _2的递推公式

假设已经选择α1,α2α1,α2，其他橙子是固定的，于是式(13)可以写成：

minα1,α212α21<x1,x1>+12α22<x2,x2>+α1α2y1y2<x1,x2>−α1−α2+α1y1∑i=3Nαiyi<xi,x1>+α2y2∑i=3Nαiyi<xi,x2>s.t.,α1y1+α2y2=−∑i=3Nαiyi=ζ0≤αi≤C,i=1,2(14)(14)minα1,α212α12<x1,x1>+12α22<x2,x2>+α1α2y1y2<x1,x2>−α1−α2+α1y1∑i=3Nαiyi<xi,x1>+α2y2∑i=3Nαiyi<xi,x2>s.t.,α1y1+α2y2=−∑i=3Nαiyi=ζ0≤αi≤C,i=1,2

对于上述优化问题来说常数因子没有意义，故省略了很多不含α1,α2α1,α2的常数项，如α3α4y3y4<x3,x4>α3α4y3y4<x3,x4>。

因为符号太多，想办法做下简化，令

Kij=K(xi,xj)=<xi,xj>(15)(15)Kij=K(xi,xj)=<xi,xj>

如上一章所讲，内积只是xi,xjxi,xj是数据集线性可分的情况下的计算方式，如果引入其他升维核函数将不是内积，这里用KijKij来代换是合适的。

(14)写成：

minα1,α212α21K11+12α22K22+α1α2y1y2K12−α1−α2+α1y1∑i=3NαiyiKi1+α2y2∑i=3NαiyiKi2s.t.,α1y1+α2y2=−∑i=3Nαiyi=ζ0≤αi≤C,i=1,2(16)(16)minα1,α212α12K11+12α22K22+α1α2y1y2K12−α1−α2+α1y1∑i=3NαiyiKi1+α2y2∑i=3NαiyiKi2s.t.,α1y1+α2y2=−∑i=3Nαiyi=ζ0≤αi≤C,i=1,2

另外我们发现∑Ni=1αiyi<xi,x1>∑i=1Nαiyi<xi,x1>不就是wTx1wTx1吗，（16）和式的后两项其实就是对x1,x2预测的计算过程去掉两个样本（用来更新α1,α2α1,α2的那两个样本），有点意思。

我们把决策公式记为：

g(x)=∑i=1NαiyiK(xi,x)+b(17)(17)g(x)=∑i=1NαiyiK(xi,x)+b

引进记号

vj=∑Ni=3αiyiKij=g(xj)−∑i=1,2αiyiKij−b(18)(18)vj=∑i=3NαiyiKij=g(xj)−∑i=1,2αiyiKij−b

要想记住它就这样记忆：vivi是对第i个样本进行预测的向量相乘(wTx−bwTx−b)忽略前两个样本。

则（16）式改写成：

minα1,α212α21K11+12α22K22+α1α2y1y2K12−α1−α2+α1y1v1+α2y2v2s.t.,α1y1+α2y2=−∑i=3Nαiyi=ζ0≤αi≤C,i=1,2(19)(19)minα1,α212α12K11+12α22K22+α1α2y1y2K12−α1−α2+α1y1v1+α2y2v2s.t.,α1y1+α2y2=−∑i=3Nαiyi=ζ0≤αi≤C,i=1,2

为什么要做这些记号，这是因为我们的关注点始终是α1,α2α1,α2这些系数什么的只要和α1,α2α1,α2无关，那么在求解的过程中只是被带来带去而已，简写之后看起来简洁，计算不易出错，什么时候要恢复成原样，把记号替换了就行。马上就可以看到疗效。

符号撸完之后，现在开始求解：

1、把α1α1用α2α2表示：α1=y1(ζ−α2y2)α1=y1(ζ−α2y2)并带入19，就只剩下变量α2α2了：

minα1,α212K11(ζ−y2α2)2+12K22α22+y2K12(ζ−y2α2)α2−y1(ζ−y2α2)−α2+v1(ζ−y2α2)+y2v2α2(20)(20)minα1,α212K11(ζ−y2α2)2+12K22α22+y2K12(ζ−y2α2)α2−y1(ζ−y2α2)−α2+v1(ζ−y2α2)+y2v2α2

注意，y21=1y12=1所以被消掉，为了好看我们把α2α2的系数全部写在了前面

2、（20）式min后面的函数对α2α2求导并令为0，得

−K11ζy2+K11α2+K22α2+K12y2ζ−2K12α2+y1y2−1−y2v1+y2v2=0⇒K11α2+K22α2−2K12α2−K11ζy2+K12ζy2+y1y2−1−v1y2+y2v2=0⇒(K11+K22−2K12)α2=y2(K11ζ−K12ζ+y2−y1+v1−v2)−K11ζy2+K11α2+K22α2+K12y2ζ−2K12α2+y1y2−1−y2v1+y2v2=0⇒K11α2+K22α2−2K12α2−K11ζy2+K12ζy2+y1y2−1−v1y2+y2v2=0⇒(K11+K22−2K12)α2=y2(K11ζ−K12ζ+y2−y1+v1−v2)

把v1，v2的原式带入：

(K11+K22−2K12)α2=y2[K11ζ−K12ζ+y2−y1+(g(x1)−∑i=1,2αiyiKi1−b)−(g(x2)−∑i=1,2αiyiKi2−b)]=y2[(g(x1)−y1)−(g(x2)−y2)+K11ζ−K12ζ−∑i=1,2αiyiKi1+∑i=1,2αiyiKi2](K11+K22−2K12)α2=y2[K11ζ−K12ζ+y2−y1+(g(x1)−∑i=1,2αiyiKi1−b)−(g(x2)−∑i=1,2αiyiKi2−b)]=y2[(g(x1)−y1)−(g(x2)−y2)+K11ζ−K12ζ−∑i=1,2αiyiKi1+∑i=1,2αiyiKi2]

接着做两个动作

a.记Ei=g(xi)−yiEi=g(xi)−yi表示预测分类和真实分类的误差

b.把ζ=∑i=1,2αoldiyiζ=∑i=1,2αioldyi带进去

(K11+K22−2K12)αnew,ucp2=y2(E1−E2)+y2(K11∑i=1,2αiyi−K12∑i=1,2αiyi−∑i=1,2αiyiKi1+∑i=1,2αiyiKi2)=y2(E1−E2)+y2(K11α2y2−K12α2y2−α2y2K21+α2y2K22)=y2(E1−E2)+(K11−2K12+K22)αold2(K11+K22−2K12)α2new,ucp=y2(E1−E2)+y2(K11∑i=1,2αiyi−K12∑i=1,2αiyi−∑i=1,2αiyiKi1+∑i=1,2αiyiKi2)=y2(E1−E2)+y2(K11α2y2−K12α2y2−α2y2K21+α2y2K22)=y2(E1−E2)+(K11−2K12+K22)α2old

漂亮，终于松一口气，接着

令η=K11−2K12+K22η=K11−2K12+K22，则有：

αnew,ucp2=y2(E1−E2)η+αold2(21)(21)α2new,ucp=y2(E1−E2)η+α2old

你可能对ucp这个标记不解，我们随后会解释.

这一步是求解α1α2α1α2的里程碑，这里意味着我们有了一个递推公式，就像梯度下降一样，对于旧值我们知道了旧值→新值的步长。如果你对之前的推导和计算毫无兴趣，那么这个递推公式可以直接用于编程了

如何计算α1α1，因为无论新值旧值都要受到约束，因此

αnew1y1+αnew2y2=αold1y1+αold2y2⇒α1newy1+α2newy2=α1oldy1+α2oldy2⇒

αnew1=αold1y1+αold2y2−αnew2y2y1=αold1−y1y2(αnew2−αold2)(22)(22)α1new=α1oldy1+α2oldy2−α2newy2y1=α1old−y1y2(α2new−α2old)

可以看到a1和a2变化的绝对值是一样的，只是方向相反。

alpha _1，alpha _2的修剪问题

这里有一个很讨厌的问题，就是递推式(21)计算出来的αnew,ucp2α2new,ucp可能太大或太小，换句话说它是有范围的。

直到现在我们只使用了约束条件中的求和等式（14），但是不等式还没有用。

我们首先应该意识到0≤αnew,ucp2≤C0≤α2new,ucp≤C，另外在迭代之前由变量旧值决定的ζ=∑i=1,2αoldiyiζ=∑i=1,2αioldyi仍然约束着新值。所以αnew,ucp2α2new,ucp的边界应综合考虑这两条。

SMO作者论文中给出了这样一张图：



由于两个自变量的线性组合固定为常数，所以ζ=∑i=1,2αoldiyiζ=∑i=1,2αioldyi决定着一条直线，加上α本身的范围决定着一个线段。新的α取值仍然要保证在这条线段上，由于αnew1α1new由αnew2α2new计算出，所以只需考虑αnew2α2new的范围控制在线段的断点即可，具体来说，

对于y1≠y2:y1≠y2:

下界L=max(0,α2−α1)L=max(0,α2−α1)

上界H=min(C,C+α2−α1)H=min(C,C+α2−α1)

对于y1=y2:y1=y2:

下界L=max(0,α2+α1−C)L=max(0,α2+α1−C)

上界H=min(C,α2+α1)H=min(C,α2+α1)

这里面的α全部是旧值。

基于上下界限定，我们需要对αnew,ucp2α2new,ucp进行修剪：

αnew2=⎧⎩⎨⎪⎪Lαnew,ucp2Hif:αnew,ucp2<Lif:L≤αnew,ucp2≥Hif:αnew,ucp2>H(23)(23)α2new={Lif:α2new,ucp<Lα2new,ucpif:L≤α2new,ucp≥HHif:α2new,ucp>H

b的更新（递推公式）

α有n个，但是b只有1个，怎么确定b的递推公式呢？

考虑（17）式，当0<αnew1<C0<α1new<C有

∵g(x1)new=∑i=1NαnewiyiK(xi,x1)+bnew1=y1∴bnew1=y1−αnew1y1K11−αnew2y2K21−∑3NαiyiKi1∵E1=αold1y1K11+αold2y2K21+∑3NαiyiKi1−y1+bold∴bnew1=−E1+αold1y1K11+αold2y2K21−αnew1y1K11−αnew2y2K21+bold∴bnew1=−E1+y1K11(αold1−αnew1)+y2K21(αold2−αnew2)+bold∵g(x1)new=∑i=1NαinewyiK(xi,x1)+b1new=y1∴b1new=y1−α1newy1K11−α2newy2K21−∑3NαiyiKi1∵E1=α1oldy1K11+α2oldy2K21+∑3NαiyiKi1−y1+bold∴b1new=−E1+α1oldy1K11+α2oldy2K21−α1newy1K11−α2newy2K21+bold∴b1new=−E1+y1K11(α1old−α1new)+y2K21(α2old−α2new)+bold

也就是说x1在边界上（x1是支持向量）时能求出这个递推式，同样如果有0<αnew2<C0<α2new<C则：

bnew2=−E2+y1K12(αold1−αnew1)+y2K22(αold2−αnew2)+boldb2new=−E2+y1K12(α1old−α1new)+y2K22(α2old−α2new)+bold

如果0<αnew1<C0<α1new<C，0<αnew2<C0<α2new<C，有bnew1=bnew2b1new=b2new。

if:α1=0或α1=C或α2=0或α2=Cif:α1=0或α1=C或α2=0或α2=C那么bnew1,bnew2b1new,b2new以及它们之间的任意一个数都是符合KKT条件的，这时选择它们的中点作为bnewbnew。

综上：

bnew=⎧⎩⎨⎪⎪−E1+y1K11(αold1−αnew1)+y2K21(αold2−αnew2)+boldbnew1+bnew22if:x1,x2 are both support vectorothers(24)(24)bnew={−E1+y1K11(α1old−α1new)+y2K21(α2old−α2new)+boldif:x1,x2 are both support vectorb1new+b2new2others

小结

至此所有的迭代公式都已经写出。简单总结如下：

SMO的基本思路是选出两个橙子，其余视为常量，重写目标函数是只保留α1,α2α1,α2相关项，重写目标函数限定∑i=1,2αiyi∑i=1,2αiyi

各种数学处理加替换，形成只有α2α2作为自变量的目标函数，然后求导志玲得到α2α2的迭代或称递推公式，自然得到α1α1的递推公式。

α2α2的递推公式是第一个里程碑。

下一个里程碑是修剪超出范围的α1,α2α1,α2。

最后一个里程碑是b的迭代，这需要考察KKT条件。

我们所有的计算都是基于KKT条件及其演变的。

接下来，我们就可以写代码来实现这一部分了。