支持向量机的SMO算法原理推导

用SMO和Python实现支持向量机的分类算法

在前面的课程中，我们曾经介绍了一种暴力算法如何在一个平面（二维）空间中找到一些样本点的分割平面。那个算法是很直接的，能够方便大家更好地理解支持向量机的概念和一些参数的意义，但是显然不是一个成熟的、具有扩展性、具有实际使用价值的方法。
九十年代末的时候，如我们前面提到的，John Platt发表了一篇重要的论文，如何用一个叫做SMO的算法实现SVM问题的快速求解。当然，如我们前面所介绍用QP、Convex optimization也可以求解SVM的问题的，但是SMO有诸多优点。我们不讨论各种算法的优缺点，今天我们一起来看看用Python如何实现一个SMO的算法。其实在John的论文里已经阐述的比较清楚，而且他还给出了一个Pseudo Code，但是不得不说的是SMO理解起来仍然是非常困难，好在网上能找到很多资料可以作为理解这一算法的很好的补充。

先让我们简单回顾一下上次的讲座我们得到的结论：

$$min_\alpha\;\psi(\alpha) = \frac{1}{2}\sum_{i,j=1}^m y^{(i)}y^{(j)}\alpha_i\alpha_jK(x_{i}, x_{j})-\sum_{i=1}^m \alpha_i$$
$$约束条件：0 \leq \alpha_i \leq C, \; i = 1,...,m$$
$$\sum_{i=1}^m\alpha_iy^{(i)} = 0$$
$$(J11)$$

其中：
$m$ 是样本数量；
$x^{(i)}$ 是第i个样本的属性（feature）向量；
$\langle x^{(i)}, x^{(j)} \rangle$ 表示第i个样本feature向量和第j个样本feature向量的内积；
$y^{(i)}$ 是第i个样本所属于的class，1或-1；
$\alpha_i$ 是附着于第i个样本的拉格朗日乘子；
$C$ 是正则化超参.

当存在属性mapping函数时，用核函数$\phi(x)$的时候，内积可以表示成$\langle \phi(x^{(i)}), \phi(x^{(j)}) \rangle$。

当使用dual form 时，判别函数将表示成：
$$f(x) = \sum^m_{i=1} \alpha_i y^{(i)} \langle x_{i}, x \rangle + b\space\space\space\space\space(J10)$$
其中b是标量，称为bias。

我们的目的就是找到”一堆“ $\alpha_i = { \alpha_1,\alpha_2,\alpha_3,...,\alpha_n \ }$ 使得目标函数最优（最小或最大）。SMO的思路是：每次选取两个$\alpha$, $\alpha_1$和$\alpha_2$，固定其它$\alpha $乘子，使得目标函数只是关于$\alpha_1$和$\alpha_2$的函数，求最优解。然后在这个基础上迭代，直到找到所有的$\alpha $，得到目标函数的最优解。$\alpha_1$和$\alpha_2$的选取是heuristically，优化是analytically。

我们先来看看有了$\alpha_1$和$\alpha_2$，其它$\alpha$固定的情况下，我们如何通过找到最优的$\alpha_1$和$\alpha_2$进而找到目标函数的最优解，这个也是SMO算法的核心。为了表达方便，先定义如下几个简化形式：

$$K(x_{i},x_{j})= K_{ij} \\
\eta=2K_{12}-K_{11}-K_{22}\space\space(J15)$$

$$f(x) = \sum\limits_{i = 1}^m {{\alpha _i}{y^{(i)}}K({x^{(i)}},x)}  + b  \space\space\space(1)$$           

$${E_i} = f({x_{i}}) - {y^{(i)}} = (\sum\limits_{i = 1}^m \alpha_iy^{(i)}K_{ij}+b) - y^{(i)} \space \space      (2)$$    

$${v_i} = \sum\limits_{j = 3}^m {{\alpha _j}{y^{(j)}}{K_{ij}}}  = f({x_{i}}) - \sum\limits_{j = 1}^2 {{\alpha _j}{y^{(j)}}{K_{ij}}}  - b   \space\space\space(3)$$      

对于$\alpha_1$,$\alpha_2$，固定其它$\alpha$后，目标函数为表示如下：

$$min_\alpha\;\psi(\alpha) = \frac{1}{2}\sum_{i,j=1}^m y^{(i)}y^{(j)}\alpha_i\alpha_jK(x_{i}, x_{j}) - \sum_{i=1}^m \alpha_i $$
$$\psi(\alpha_1,\alpha_2) =-[\alpha_1 + \alpha_2 -\frac{1}{2}K_{11}\alpha_1^{2} -\frac{1}{2}K_{22}\alpha_2^{2} -\alpha_1 \alpha_2y^{(1)}y^{(2)}K_{12}-y^{(1)}\alpha_1\sum\limits_{j=3}^{m}y^{(j)}\alpha_jK_{(1j)}-y^{(2)}\alpha_2\sum\limits_{j=3}^{m}y^{(j)}\alpha_jK_{2j}+\sum\limits_{i=3}^{m}\alpha_i-\frac{1}{2}\sum\limits_{i,j=3}^{m}\alpha_i\alpha_j y^{(i)}y^{(j)}K_{ij}]$$ 
$$\psi(\alpha_1,\alpha_2) =-[ \alpha_1 + \alpha_2 -\frac{1}{2}K_{11}\alpha_1^{2} -\frac{1}{2}K_{22}\alpha_2^{2} -\alpha_1 \alpha_2y^{(1)}y^{(2)}K_{12}-y^{(1)}\alpha_1v_1-y^{(2)}\alpha_2v_2+constant] \space\space\space (3_1)$$ 

根据等式约束，

$$\sum\limits_{i = 1}^m {{\alpha _i^{(old)}}{y^{(i)}}} = \sum\limits_{i = 1}^m {{\alpha _i}{y^{(i)}}} = 0$$

我们有：

$$\alpha_1 + s\alpha_2 = constant = \alpha_1^{(old)} + s\alpha_2^{(old)} = \gamma \space\space\space其中 s = y^{(1)}y^{(2)}\space\space\space(4)$$

$$从而有，\alpha_1 = \gamma - s\alpha_2,\space\space\space\space\space\space 后面我们用这个公式在知道\alpha_2^{(new)}的情况下计算\alpha_1^{(new)}。$$

把上式带入式（$3_1$）使之成为只针对$\alpha_2$的优化问题，如下：

$$\psi(\alpha_2) = -[\gamma - s\alpha_2+\alpha_2 - \frac{1}{2}K_{11}(\gamma - s\alpha_2)^{2}-\frac{1}{2}K_{22}\alpha_2^{2}-sK_{12}(\gamma - s\alpha_2)\alpha_2 - y^{(1)}(\gamma - s\alpha_2)v_1 - y^{(2)}\alpha_2v_2+constant]$$

$$= -[-\frac{1}{2}K_{11}(\gamma^{2}-2\gamma s\alpha_2+\alpha_2^{2})-\frac{1}{2}K_{22}\alpha_2^{2}+s^{2}K_{12}\alpha_2^{2}+(1-s-sK_{12}\gamma)\alpha_2 + \alpha_2y^{(2)}(v_1-v_2) + constant]$$

$$=-[\frac{1}{2}(2K_{12} - K_{11} - K_{22})\alpha_2^{2} + (1 - s+K_{11}s\gamma - K_{12}s\gamma + y^{(2)}v_1 - y^{(2)}v_2)\alpha_2 -\frac{1}{2}K_{11}\gamma^2+ constant]\space\space\space\space\space (5)$$

这个二维函数的驻点（stationary points）满足：

$$\frac{\mathrm dW(\alpha_2)}{\mathrm d\alpha_2} = (K_{11} + K_{22}-2K_{12})\alpha_2 - (1 - s+K_{11}s\gamma - K_{12}s\gamma + y^{(2)}v_1 - y^{(2)}v_2) = 0$$

$$从而有，\alpha_2^{(new,unc)}(K_{11}+K_{22}-2K_{12}) = 1 - s + K_{11}s\gamma - K_{12}s\gamma + y^{(2)}v_1 - y^{(2)}v_2 \\ = 1-s + (K_{11}-K_{12})s\gamma + y^{(2)}(v_1-v_2)\space\space\space\space（6）$$

式（6）两边同乘$y^{(2)}$，则有：

$$\alpha_2^{(new,unc)}\eta y^{(2)} = y^{(2)} - y^{(1)} + (K_{11}-K_{12})y^{(1)}\gamma+v_1-v_2 \\=y^{(2)} - y^{(1)} + (K_{11}-K_{12})y^{(1)}\gamma + (f(x_1)-\sum\limits_{j=1}^{2}y^{(j)}\alpha_jK_{1j}-b) - (f(x_2)-\sum\limits_{j=1}^{2}y^{(j)}\alpha_jK_{2j}-b)\space\space\space\space(7)$$

$$根据（4）有：y^{(1)}\gamma = y^{(1)}(\alpha_1+s\alpha_2)=\alpha_1y^{(1)}+\alpha_2y^{(2)}\\ 又因为：\sum\limits_{j=1}^{2}y^{(j)}\alpha_jK_{2j}-\sum\limits_{j=1}^{2}y^{(j)}\alpha_jK_{1j} = y^{(1)}\alpha_1K_{21}+y^{(2)}\alpha_2K_{22} - y^{(1)}\alpha_1K_{11}-y^{(2)}\alpha_2K_{12}\space\space 两式带入（7）有：$$

$$\alpha_2^{(new,unc)}\eta y^{(2)} = y^{(2)}-y^{(1)}+(K_{11}-K_{12})(\alpha_1 y^{(1)}+\alpha_2 y^{(2)})+y^{(1)}\alpha_1K_{21}+y^{(2)}\alpha_2K_{22}-y^{(1)}\alpha_1K_{11}-y^{(2)}\alpha_2K_{12}+f(x_1)-f(x_2) \\
=y^{(2)}-y^{(1)}+f(x_1)-f(x_2)+y^{(2)}\alpha_2K_{11}-y^{(2)}\alpha_2K_{12}+y^{(2)}\alpha_2K_{22}-y^{(2)}\alpha_2K_{12} \\ = y^{(2)}\alpha_2\eta + (f(x_1)-y^{(1)})-(f(x_2)-y^{(2)})$$

从而我们就得出一个非常关键的计算$\alpha_2^{(new)}$公式：

$$\alpha_2^{(new,unc)} = \alpha_2^{(old)} + \frac{y^{(2)}\{E_2^{(old)}-E_1^{(old)}\}}{\eta} \space\space\space(J16)$$ 

同时根据式子（4）我们可以得出计算$\alpha_1^{(new)}$的公式如下：

$$\alpha_1^{(new)} = \alpha_1^{(old)}+y^{(1)}y^{(2)}(\alpha_2^{(old)}-\alpha_2^{(new)})  \space\space\space(J18)$$

上面的（8）式并没有考虑$\alpha_2$的取值范围，我们需要把它修剪(“夹”可能更准确)（clipped）到它的取值范围以内。那么什么是它的取值范围呢？让我们看下两图。

根据不等约束，$\alpha$在0和C之间（inclusive），有根据相等约束有$\alpha_1 y^{(1)}+\alpha_2 y^{(2)}=\zeta$。$\alpha_2$的取值范围要么是Line1上的【$\alpha_2 -\alpha_1$, C】，要么是line2上的【0，$C+\alpha_2-\alpha_1$】，如Mr.Platt的论文，我们需要将$\alpha_2$的取值范围表示成$L\leq\alpha_2\leq H$的形式。怎么做呢？令
$$L=max(0,\alpha_2-\alpha_1)，\\H=min(C,C+\alpha_2-\alpha_1)\space\space (J13)$$，当L=0时，说明$\alpha_2-\alpha_1$<0，所以C+$\alpha_2-\alpha_1$<C，所以H=C+$\alpha_2-\alpha_1$，则$L\leq\alpha_2\leq H$表示$\alpha_2$位于Line2上，它的取值范围是$0\leq\alpha_2\leq C+\alpha_2-\alpha_1$。当L=$\alpha_2-\alpha_1$时，说明$\alpha_2-\alpha_1$>0，所以C+$\alpha_2-\alpha_1$>C，所以H=C，则$L\leq\alpha_2\leq H$表示$\alpha_2$位于Line1上，它的取值范围是$\alpha_2-\alpha_1\leq\alpha_2\leq C$。
所以当$y^{(1)}$和$y^{(2)}$异号时，L=max(0,$\alpha_2-\alpha_1$)，H=min(C,C+$\alpha_2-\alpha_1$)，用$L\leq\alpha_2\leq H$就能够很方便地表示$\alpha_2$的取值范围。$y^{(1)}$和y^{(2)}同号时，情况类似，不再赘述，但给出结论如下：
$$L=max(0,\alpha_2+\alpha_1-C), H=min(C,\alpha_2,\alpha_1)\space\space(J14)$$
从而裁剪后，$\alpha_2^{new}$表示如下：

$$\alpha_2^{new} = \left\{ {\begin{array}{*{20}{c}}{H,\quad \alpha _2^{new,unc} > H}\\{\alpha_2^{new,unc},\quad L\leq \alpha _2^{new,unc}\leq H}\\{L,\quad \alpha_2^{new,unc}<L}\end{array}} \right.\\(J17)$$

得到了新的$\alpha_2^{new}和\alpha_1^{new}$后，我们还要计算并更新$b^{new}和E_i^{new}$，根据最上面的（1），（2），（3）我们有：

$$b_1^{new}=y^{(1)}-\sum\limits_{i=3}^{m}\alpha_iy^{(i)}K_{i1}-\alpha_1^{new}y^{(1)}K_{11}-\alpha_2^{new}y^{(2)}K_{21}\space\space\space\space(10)$$

$$E_1=\sum\limits_{i=3}^m\alpha_iy^{(i)}K_{i1}+\alpha_1^{old}y^{(1)}K_{11}+\alpha_2^{old}y^{(2)}K_{21}+b^{old}-y^{(1)}\space\space\space\space(11)$$

合并以上（10）和（11），有：

$$b_1^{new}=-E_1-y^{(1)}K_{11}(\alpha_1^{new}-\alpha_1^{old})-y^{(2)}K_{21}(\alpha_2^{new}-\alpha_2^{old})+b^{old}\\(J20)$$

从而有：$$b_2^{new}=-E_2-y^{(1)}K_{12}(\alpha_1^{new}-\alpha_1^{old})-y^{(2)}K_{22}(\alpha_2^{new}-\alpha_2^{old})+b^{old}\\(J21)$$

更新b的规则:如果$0<\alpha_1^{new}<C$，则$b^{new} = b_1^{new}$；如果$0<\alpha_2^{new}<C$，则$b^{new} = b_2^{new}$；如果都不满足则$b^{new} = (b_1^{new} + b_2^{new})/2$，然后更新$E_i$:

$$E_i^{new} = \sum\limits_S y^{(j)}\alpha_jK(x_i,x_j)+b^{new}-y^{(i)}\space\space其中S是所有支持向量X_j的集合\\(12)$$

在特殊情况下，$\eta$的值可能为负，即如果核函数K不遵从（not obey）Mercer定理的时候，负的$\eta值$就会出现，导致目标函数非正定（indefinite）。当一个以上的训练样本有相同的向量x时，$\eta值还可能为0。SMO算法在这些特殊情况下目标函数将人为地收敛到端点的值，即H或L。这样将$\alpha_2=L或\alpha_2=H$分别代入目标函数中即可求得极值。具体计算公式如下：

$$f_1=y^{(1)}(E_1+b)-\alpha_1K(x_1,x_1)-s\alpha_2K(x_1,x_2),\\
f_2=y^{(2)}(E_2+b)-s\alpha_1K(x_1,x_2)-\alpha_2K(x_2,x_2),\\
L_1 = \alpha_1+s(\alpha_2-L),\\
H_1 = \alpha_1+s(\alpha_2-H),\\
\psi_L = L_1f_1+Lf_2+\frac{1}{2}L_1^{2}K(x_1,x_1)+\frac{1}{2}L^2K(x_2,x_2)+sLL_1K(x_1,x_2),\\
\psi_H = H_1f_1+Hf_2+\frac{1}{2}H_1^{2}K(x_1,x_1)+\frac{1}{2}H^2K(x_2,x_2)+sHH_1K(x_1,x_2).\\
(J19)$$

最后我们看一下算法中如何选择$\alpha_2和\alpha_1$，看一下细化的KKT条件：

$$\alpha_i=0\leftrightarrow y^{(i)}f(x_i)\geq1\\0<\alpha_i<C\leftrightarrow y^{(i)}f(x_i)=1\\\alpha_i =C\leftrightarrow y^{(i)}f(x_i)\leq 1\\(J12)$$

这个KKT条件说明，在两条间隔线外面的点所对应的拉格朗日乘子$\alpha_i$为0，在两条间隔线里面的对应乘子$\alpha_i$为C，只有在两条间隔线上对应乘子的值不为0，在0和C之间。
首先选择违反$0<\alpha_i<C\leftrightarrow y^{(i)}f(x_i)=1$这个条件最严重的点。如果$0<\alpha_i<C$对应的点都满足KKT条件，再选择违反$\alpha_i=0\leftrightarrow y^{(i)}f(x_i)\geq1$ 和 $\alpha_i =C\leftrightarrow y^{(i)}f(x_i)\leq 1$的点。第二个变量的选择标准是让$|E_1−E_2|$有足够大的变化。由于$\alpha_1$定了的时候,$E_1$也确定了，所以要想$|E_1−E_2|$最大，只需要在$E_1$为正时，选择最小的$E_i$作为$E_2$， 在$E_1$为负时，选择最大的$E_i$作为$E_2$，可以将所有的$E_i$保存下来加快迭代。