机器学习系列 10：支持向量机 02 - SMO（序列最小化）

---

针对 “支持向量机” 将分为以下几个部分进行介绍：

1. 支持向量机 01 - 线性可分支持向量机和线性支持向量机
2. 支持向量机 02 - SMO（序列最小化）（本篇）
3. 支持向量机 03 - 非线性支持向量机
4. 支持向量机 04 - SVR（支持向量回归）

---

  本内容将介绍 SMO（序列最小化）算法，包含详细公式推导以及 Python 代码实现。

  在学习本内容前，需要对支持向量机基本理论知识有一定了解，如果您还不了解，可以参阅 支持向量机 01 - 线性可分支持向量机和线性支持向量机。

  我们知道，支持向量机的学习问题可以形式化为求解凸二次规划问题。这样的凸二次规划问题具有全局最优解，并且有许多最优化算法可以用于这一问题的求解。但是当训练样本容量很大时，这些算法往往变得非常低效，以致无法使用。

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

一、SMO 算法描述

  整个 SMO 算法包括两个部分：求解两个变量二次规划的解析方法和选择变量的启发式方法。

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

1.1.1 问题转化

  在 支持向量机 01 - 线性可分支持向量机和线性支持向量机 中讲解了 SVM 的基本原理，了解到 SMO 算法要解如下凸二次规划的对偶问题

$$\begin{array}{rl}& \underset{\alpha }{\max }\left(\sum _{i=1}^m{\alpha }_i-\frac{1}{2}\sum _{i=1}^m\sum _{j=1}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}\right)\\ & \\ & \begin{array}{rl}s.t.& \sum _{i=1}^m{\alpha }_i{y}^{(i)}=0\\ & \\ & {\alpha }_i\ge 0,i=1,2,\cdots ,m\end{array}\end{array}\text{\hspace{1em}(1)}$$

  添加一个负号，将最大值问题转换成最小值问题

$$\begin{array}{rl}& \underset{\alpha }{\min }\left(\frac{1}{2}\sum _{i=1}^m\sum _{j=1}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}-\sum _{i=1}^m{\alpha }_i\right)\\ & \\ & \begin{array}{rl}s.t.& \sum _{i=1}^m{\alpha }_i{y}^{(i)}=0\\ & \\ & {\alpha }_i\ge 0,i=1,2,\cdots ,m\end{array}\end{array}\text{\hspace{1em}(2)}$$

1.1.2 转换为二元函数

  假设选择的两个变量为 ${\alpha }_1$ 和 ${\alpha }_2$，其他变量 ${\alpha }_i(i=3,4,\cdots ,m)$ 是固定的，式（2）可写为

$$\begin{array}{rl}W({\alpha }_1,{\alpha }_2)& =\frac{1}{2}\sum _{i=1}^m\sum _{j=1}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}-\sum _{i=1}^m{\alpha }_i\\ & \\ & =\frac{1}{2}\sum _{i=1}^m\left(\sum _{j=1}^2{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}+\sum _{j=3}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}\right)-{\alpha }_i-{\alpha }_2-\sum _{i=3}^m{\alpha }_i\\ & \\ & =\frac{1}{2}\sum _{i=1}^2\left(\sum _{j=1}^2{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}+\sum _{j=3}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}\right)\\ & +\frac{1}{2}\sum _{i=3}^m\left(\sum _{j=1}^2{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}+\sum _{j=3}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}\right)-{\alpha }_i-{\alpha }_2-\sum _{i=3}^m{\alpha }_i\\ & \\ & =\frac{1}{2}\sum _{i=1}^2\sum _{j=1}^2{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}+\sum _{i=1}^2\sum _{j=3}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}\\ & +\frac{1}{2}\sum _{i=3}^m\sum _{j=3}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}-{\alpha }_1-{\alpha }_2-\sum _{i=3}^m{\alpha }_i\\ & \\ & =\frac{1}{2}{\alpha }_1^2{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}+\frac{1}{2}{\alpha }_2^2{\mathbf{x}}^{(2)}{\mathbf{x}}^{(2)}+y^{(1)}{y}^{(2)}{\alpha }_1{\alpha }_2{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)}+y^{(1)}{\alpha }_1\sum _{j=3}^m{\alpha }_j{y}^{(j)}{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(j)}\\ & +y^{(2)}{\alpha }_2\sum _{j=3}^m{\alpha }_j{y}^{(j)}{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(j)}-{\alpha }_1-{\alpha }_2+\frac{1}{2}\sum _{i=3}^m\sum _{j=3}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}-\sum _{i=3}^m{\alpha }_i\end{array}\text{\hspace{1em}(3)}$$

  为了描述方便，我们定义如下符号：

$$f({\mathbf{x}}^{(i)})=\sum _{j=1}^m{\alpha }_j{y}^{(j)}{\mathbf{x}}^{(j)}\cdot {\mathbf{x}}^{(i)}+b=\sum _{j=1}^m{\alpha }_j{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}+b\text{\hspace{1em}(4)}$$

$$v_i=\sum _{j=3}^m{\alpha }_j{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}=f({\mathbf{x}}^{(i)})-\sum _{j=1}^2{\alpha }_j{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}-b\text{\hspace{1em}(5)}$$

  根据式（4）和（5），目标函数式（3）可写为

$$\begin{array}{rl}W({\alpha }_1,{\alpha }_2)& =\frac{1}{2}{\alpha }_1^2{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}+\frac{1}{2}{\alpha }_2^2{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(2)}+y^{(1)}{y}^{(2)}{\alpha }_1{\alpha }_2{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)}\\ & +y^{(1)}{\alpha }_1{v}_1+y^{(2)}{\alpha }_2{v}_2-{\alpha }_1-{\alpha }_2+constant\end{array}\text{\hspace{1em}(6)}$$

其中

$$constant=\frac{1}{2}\sum _{i=3}^m\sum _{j=3}^m{\alpha }_i{\alpha }_j{y}^{(i)}{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}-\sum _{i=3}^m{\alpha }_i\text{\hspace{1em}(7)}$$

  因为 $constant$ 部分对 ${\alpha }_1$ 和 ${\alpha }_2$ 来说，属于常数项，在求导的时候，直接变为 0，所以我们不需要关心这部分的内容。

1.1.3 转换为一元函数

  因为有 ${\sum }_{i=1}^m{\alpha }_i{y}^{(i)}=0$，可得

$${\alpha }_1{y}^{(1)}+{\alpha }_2{y}^{(2)}=-\sum _{i=3}^m{\alpha }_i{y}^{(i)}\text{\hspace{1em}(8)}$$

  对变量 ${\alpha }_1$ 和 ${\alpha }_2$ 来说，${\alpha }_i$ 和 $y^{(i)}$（$i=3,4,\cdots ,m$ ）可看作常数项，因此有

$${\alpha }_1{y}^{(1)}+{\alpha }_2{y}^{(2)}=B\text{\hspace{1em}(9)}$$

其中，$B$ 为一个常数。将等式（9）两边同时乘以 $y^{(1)}$，得

$${\alpha }_1=B{y}^{(1)}-{\alpha }_2{y}^{(1)}{y}^{(2)}=\left(B-{\alpha }_2{y}^{(2)}\right)y^{(1)}\text{\hspace{1em}(10)}$$

将其带入式（6），得到只含变量 ${\alpha }_2$ 的目标函数

$$\begin{array}{rl}W({\alpha }_2)& =\frac{1}{2}(B-{\alpha }_2{y}^{(2)}{)}^2{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}+\frac{1}{2}{\alpha }_2^2{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(2)}+y^{(2)}(B-{\alpha }_2{y}^{(2)}){\alpha }_2{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)}\\ & +(B-{\alpha }_2{y}^{(2)})v_1+y^{(2)}{\alpha }_2{v}_2-(B-{\alpha }_2{y}^{(2)})y^{(1)}-{\alpha }_2+constant\end{array}\text{\hspace{1em}(11)}$$

1.1.4 求解 ${\alpha }_2^{new,unclipped}$

  将式（11）对 ${\alpha }_2$ 求导，得

$$\begin{array}{rl}\frac{\partial W({\alpha }_2)}{\partial {\alpha }_2}& ={\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}{\alpha }_2+{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(2)}{\alpha }_2-2{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)}{\alpha }_2\\ & -B{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}{y}^{(2)}+B{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)}{y}^{(2)}+y^{(1)}{y}^{(2)}-v_1{y}^{(2)}+v_2{y}^{(2)}-1\end{array}\text{\hspace{1em}(12)}$$

令其为 $0$，得

$${\alpha }_2^{new}=\frac{y^{(2)}\left(y^{(2)}-y^{(1)}+B({\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}-{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)})+v_1-v_2\right)}{{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}+{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(2)}-2{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)}}\text{\hspace{1em}(13)}$$

  定义如下符号

$$E_i=f({\mathbf{x}}^{(i)})-y^{(i)}\text{\hspace{1em}(14)}$$

$$\eta ={\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}+{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(2)}-2{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)}\text{\hspace{1em}(15)}$$

其中 $E_i$ 为误差项，$\eta$ 为学习速率。

  由于已知

$$B={\alpha }_1^{old}{y}^{(1)}+{\alpha }_2^{old}{y}^{(2)}\text{\hspace{1em}(16)}$$

$$v_i=f\left({\mathbf{x}}^{(i)}\right)-\sum _{j=1}^2{\alpha }_j{y}^{(j)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(j)}-b\text{\hspace{1em}(17)}$$

将式（16）和（17）带入式（13）， 将 ${\alpha }_2^{new,unclipped}$ 化简后得

$${\alpha }_2^{new,unclipped}={\alpha }_2^{old}+\frac{y^{(2)}(E_1-E_2)}{\eta }\text{\hspace{1em}(18)}$$

1.1.5 求解 ${\alpha }_2^{new}$（对原始 ${\alpha }_2^{new,unclipped}$ 解进行修剪）

  上面求出的 ${\alpha }_2^{new,unclipped}$ 解是没有经过修剪的。我们知道 ${\alpha }_i$ 存在以下约束条件

$$\{\begin{array}{c}\begin{array}{rl}& 0<{\alpha }_i<C\\ & \\ & {\alpha }_1{y}^{(1)}+{\alpha }_2{y}^{(2)}=B\end{array}\end{array}\text{\hspace{1em}(19)}$$

可知 ${\alpha }_2^{new}$ 的取值需要满足一定条件，假设为

$$L\le {\alpha }_2^{new}\le H\text{\hspace{1em}(20)}$$

  图-1 在二维平面上直观地表达了式（19）中的约束条件，从而可知

• 当 $y^{(1)}̸=y^{(2)}$ 时，存在 ${\alpha }_2^{old}-{\alpha }_1^{old}=k$，得

$$L=max(0,{\alpha }_2^{old}-{\alpha }_1^{old}),H=min(C,C+{\alpha }_2^{old}-{\alpha }_1^{old})\text{\hspace{1em}(21)}$$

• 当 $y^{(1)}=y^{(2)}$ 时，存在 ${\alpha }_2^{old}+{\alpha }_1^{old}=k$，得

$$L=max(0,{\alpha }_2^{old}+{\alpha }_1^{old}-C),H=min(C,{\alpha }_2^{old}+{\alpha }_1^{old})\text{\hspace{1em}(22)}$$

图-1

  则经过修剪后，${\alpha }_2^{new}$ 的解为

$${\alpha }_2^{new}=\{\begin{array}{c}\begin{array}{rlr}& H,& {\alpha }_2^{new,unclipped}>H\\ & & \\ & {\alpha }_2^{new,unclipped},& L\le {\alpha }_2^{new,unclipped}\le H\\ & & \\ & L,& {\alpha }_2^{new,unclipped}<L\end{array}\end{array}\text{\hspace{1em}(23)}$$

1.1.6 求解 ${\alpha }_1^{new}$

  因为存在以下等式条件

$${\alpha }_1^{old}{y}^{(1)}+{\alpha }_2^{old}{y}^{(2)}={\alpha }_1^{new}{y}^{(1)}+{\alpha }_2^{new}{y}^{(2)}\text{\hspace{1em}(24)}$$

从而得出

$${\alpha }_1^{new}={\alpha }_1^{old}+y^{(1)}{y}^{(2)}({\alpha }_2^{old}-{\alpha }_2^{new})\text{\hspace{1em}(25)}$$

1.1.7 求解阈值 $b^{new}$

  在求得 ${\alpha }_1^{new}$ 和 ${\alpha }_2^{new}$ 后，需要根据它们更新 $b$。

1. 如果 $0<{\alpha }_1^{new}<C$
     由 KKT 条件可知（在 支持向量机 01 - 线性可分支持向量机和线性支持向量机 中 2.2 节有进行介绍），当 $0<{\alpha }_i<C$ 时，这个样本点是支持向量。因此，满足
   $$y^{(1)}({\mathbf{w}}^T{\mathbf{x}}^{(1)}+b)=1\text{\hspace{1em}(26)}$$

   上面公式两边同时乘以 $y^{(1)}$ 得
   $${\mathbf{w}}^T{\mathbf{x}}^{(1)}+b=y^{(1)}\text{\hspace{1em}(27)}$$

   将 $\mathbf{w}={\sum }_{i=1}^m{\alpha }_i{y}^{(i)}{\mathbf{x}}^{(i)}$ 代入，得

   $$\sum _{i=1}^m{\alpha }_i{y}^{(i)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(1)}+b=y^{(1)}\text{\hspace{1em}(28)}$$

     因为我们是根据 ${\alpha }_1^{new}$ 和 ${\alpha }_2^{new}$ 的值更新 $b$，整理可得：

   $$b_1^{new}=y^{(1)}-\sum _{i=3}^m{\alpha }_i{y}^{(i)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(1)}-{\alpha }_1^{new}{y}^{(1)}{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}-{\alpha }_2^{new}{y}^{(2)}{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(1)}\text{\hspace{1em}(29)}$$

   其中有

   $$\begin{array}{rl}{y}^{(1)}-\sum _{i=3}^m{\alpha }_i{y}^{(i)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(1)}& =y^{(1)}-\sum _{i=1}^m{\alpha }_i{y}^{(i)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(1)}-b^{old}\\ & +{\alpha }_1^{old}{y}^{(1)}{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}+{\alpha }_2^{old}{y}^{(2)}{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(1)}+b^{old}\\ & \\ & =y^{(1)}-f\left({\mathbf{x}}^{(1)}\right)+{\alpha }_1^{old}{y}^{(1)}{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}+{\alpha }_2^{old}{y}^{(2)}{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(1)}+b^{old}\\ & \\ & =-E_1+{\alpha }_1^{old}{y}^{(1)}{\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}+{\alpha }_2^{old}{y}^{(2)}{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(1)}+b^{old}\end{array}\text{\hspace{1em}(30)}$$

   将式（30）代入式（29）得

   $$b_1^{new}=b^{old}-E_1-y^{(1)}({\alpha }_1^{new}-{\alpha }_1^{old}){\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}-y^{(2)}({\alpha }_2^{new}-{\alpha }_2^{old}){\mathbf{x}}^{(2)}{\mathbf{x}}^{(1)}\text{\hspace{1em}(31)}$$

2. 如果 $0<{\alpha }_2^{new}<C$
     按照上面的步骤，同样可求得 $b_2^{new}$
   $$b_2^{new}=b^{old}-E_2-y^{(1)}({\alpha }_1^{new}-{\alpha }_1^{old}){\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)}-y^{(2)}({\alpha }_2^{new}-{\alpha }_2^{old}){\mathbf{x}}^{(2)}{\mathbf{x}}^{(2)}\text{\hspace{1em}(32)}$$

     如果 ${\alpha }_1^{new}$ 和 ${\alpha }_2^{new}$ 同时满足条件 $0<{\alpha }_i^{new}<C$，则 $b^{new}=b_1^{new}=b_2^{new}$；如果 ${\alpha }_1^{new}$ 和 ${\alpha }_2^{new}$ 是 $0$ 或者 $C$，那么 $b_1^{new}$ 和 $b_2^{new}$ 以及它们之间的数都是符合 KKT 条件的阈值，这是选择它们的重点作为 $b^{new}$。则 $b^{new}$ 为

$$b^{new}=\{\begin{array}{ccc}& b_1^{new},& 0<{\alpha }_1^{new}<C\\ & & \\ & b_2^{new},& 0<{\alpha }_2^{new}<C\\ & & \\ & \frac{(b_1^{new}+b_2^{new})}{2},& otherwise\end{array}\text{\hspace{1em}(33)}$$

1.1.8 变量更新算法的步骤

  根据上面讲解的内容，我们来整理一下变量更新算法的步骤：

• 步骤 1：计算误差

$$E_1=f({\mathbf{x}}^{(1)})-y^{(1)}=\sum _{i=1}^m{\alpha }_i{y}^{(i)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(1)}+b-y^{(1)}$$

$$E_2=f({\mathbf{x}}^{(2)})-y^{(2)}=\sum _{i=1}^m{\alpha }_i{y}^{(i)}{\mathbf{x}}^{(i)}\cdot {\mathbf{x}}^{(2)}+b-y^{(2)}$$

• 步骤 2：计算 ${\alpha }_2^{new}$ 取值范围

$$\{\begin{array}{c}\begin{array}{rl}& if(y^{(1)}̸=y^{(2)})L=max(0,{\alpha }_2^{old}-{\alpha }_1^{old}),H=min(C,C+{\alpha }_2^{old}-{\alpha }_1^{old})\\ & \\ & if(y^{(1)}=y^{(2)})L=max(0,{\alpha }_2^{old}+{\alpha }_1^{old}-C),H=min(C,{\alpha }_2^{old}+{\alpha }_1^{old})\end{array}\end{array}$$

• 步骤 3：计算 $\eta$

$$\eta ={\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}+{\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(2)}-2{\mathbf{x}}^{(1)}{\mathbf{x}}^{(2)}$$

• 步骤 4：求解 ${\alpha }_2^{new,unclipped}$

$${\alpha }_2^{new,unclipped}={\alpha }_2^{old}+\frac{y^{(2)}(E_1-E_2)}{\eta }$$

• 步骤 5：求解 ${\alpha }_2^{new}$

$${\alpha }_2^{new}=\{\begin{array}{c}\begin{array}{rlr}& H,& {\alpha }_2^{new,unclipped}>H\\ & & \\ & {\alpha }_2^{new,unclipped},& L\le {\alpha }_2^{new,unclipped}\le H\\ & & \\ & L,& {\alpha }_2^{new,unclipped}<L\end{array}\end{array}$$

• 步骤 6：求解 ${\alpha }_1^{new}$

$${\alpha }_1^{new}={\alpha }_1^{old}+y^{(1)}{y}^{(2)}({\alpha }_2^{old}-{\alpha }_2^{new})$$

• 步骤 7：求解 $b_1^{new}$ 和 $b_2^{new}$

$$b_1^{new}=b^{old}-E_1-y^{(1)}({\alpha }_1^{new}-{\alpha }_1^{old}){\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(1)}-y^{(2)}({\alpha }_2^{new}-{\alpha }_2^{old}){\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(1)}$$

$$b_2^{new}=b^{old}-E_2-y^{(1)}({\alpha }_1^{new}-{\alpha }_1^{old}){\mathbf{x}}^{(1)}\cdot {\mathbf{x}}^{(2)}-y^{(2)}({\alpha }_2^{new}-{\alpha }_2^{old}){\mathbf{x}}^{(2)}\cdot {\mathbf{x}}^{(2)}$$

• 步骤 8：求解 $b^{new}$

$$b^{new}=\{\begin{array}{ccc}& b_1^{new},& 0<{\alpha }_1^{new}<C\\ & & \\ & b_2^{new},& 0<{\alpha }_2^{new}<C\\ & & \\ & \frac{(b_1^{new}+b_2^{new})}{2},& otherwise\end{array}$$

  前面进行介绍时，我们有假设选择两个变量 ${\alpha }_1$ 和 ${\alpha }_2$，但是在实际实现算法时，我们应该如何选择这两个变量呢？下面就来介绍一下。

1.2 变量 ${\alpha }_1$ 和 ${\alpha }_2$ 的选择方法

1.2.1 第 1 个变量的选择

  SMO 称选择第 1 个变量的过程为外层循环。外层循环在训练样本中选取违反 KKT 条件最严重的样本点，并将其对应的变量作为第 1 个变量。具体地，检验训练样本点是否满足 KKT 条件，即（在 支持向量机 01 - 线性可分支持向量机和线性支持向量机 中 2.2 节有进行介绍）

$$\begin{array}{r}{\alpha }_i=0\iff y^{(i)}f({\mathbf{x}}^{(i)})\ge 1\\ \\ 0<{\alpha }_i<C\iff y^{(i)}f({\mathbf{x}}^{(i)})=1\\ \\ {\alpha }_i=C\iff y^{(i)}f({\mathbf{x}}^{(i)})\le 1\end{array}\text{\hspace{1em}(34)}$$

  从而可得知，以下几种情况将会不满足 KKT 条件

$$\begin{array}{rl}& y^{(i)}f\left({\mathbf{x}}^{(i)}\right)\ge 1and{\alpha }_i>0\\ & \\ & y^{(i)}f\left({\mathbf{x}}^{(i)}\right)=1and({\alpha }_i=0or{\alpha }_i=C)\\ & \\ & y^{(i)}f\left({\mathbf{x}}^{(i)}\right)\le 1and{\alpha }_i<C\end{array}\text{\hspace{1em}(35)}$$

  在检验过程中，外层循环首先遍历所有满足条件 $0<{\alpha }_i<C$ 的样本点，即在间隔边界上的支持向量点，检验它们是否满足 KKT 条件。如果这些样本点都满足 KKT 条件，那么遍历整个训练集，检验它们是否满足 KKT 条件。

1.2.2 第 2 个变量的选择

  SMO 称选择第 2 个变量的过程为内循环。假设在外层循环中已经找到第 1 个变量 ${\alpha }_1$，现在要在内层循环中找第 2 个变量 ${\alpha }_2$。第 2 个变量选择的标准是希望能使 ${\alpha }_2$ 有足够大的变化。

  由式（18）可知，${\alpha }_2^{new}$ 是依赖于 $|E_1-E_2|$ 的，为了加快计算速度，一种简单的做法是选择 ${\alpha }_2$，使其对应的 $|E_1-E_2|$ 最大。为了节省计算时间，可以将所有 $E_i$ 值保存在一个列表中。

  在特殊情况下，如果内层循环通过以上方法选择的 ${\alpha }_2$ 不能使目标函数有足够的下降（等价于 $|{\alpha }_2^{new}-{\alpha }_2|$ 很小），那么可以采用以下启发式规则继续选择 ${\alpha }_2$。遍历在间隔边界上的支持向量点，依次将其对应的变量作为 ${\alpha }_2$ 试用，知道目标函数有足够的下降。若仍然找不到合适的 ${\alpha }_2$，那么遍历训练数据集；若还是找不到合适的 ${\alpha }_2$，则放弃第 1 个 ${\alpha }_2$，再通过外层循环寻求另外的 ${\alpha }_1$。

  由 Osuna 定理可知，只需选取的 ${\alpha }_1$ 和 ${\alpha }_2$ 中有一个不满足 KKT 条件，目标函数就会在迭代后减小。

1.3 SMO 算法的步骤

我们总结一下 SMO 算法的步骤：

• 步骤 1：初始化 $\alpha$ 和 $b$，比如初始化 $\alpha$ 为零向量，$b$ 为 0。（注意这里的 $\alpha$ 是一个列向量）
• 步骤 2：选取优化变量 ${\alpha }_1$ 和 ${\alpha }_2$，然后更新 ${\alpha }_1$、${\alpha }_2$ 和 $b$。
• 步骤 3：如果满足停止条件（即前面提到的“所有变量的解都满足此最优化问题的 KKT 条件”）则结束；否则，跳到步骤 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]
        # 初始化 alphas 为零向量，初始化 b 为 0
        self.alphas = np.mat(np.zeros((self.m, 1)))
        self.b = 0
        # e_cache 用于缓存误差值
        # 第一列给出 e_cache 是否有效的标志位，第二列给出实际的误差值
        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
        # 保存最大的 |ei - ej|
        max_delta_e = 0
        ej = 0
        # 将 ei 更新到缓存 e_cache 中
        os.e_cache[i] = [1, ei]
        # 返回 os.e_cache 中 e 为非零的索引值列表
        valid_e_cache_list = np.nonzero(os.e_cache[:, 0].A)[0]
        # 第一次执行这个函数时，会执行 else 语句
        if (len(valid_e_cache_list)) > 1:
            # 通过遍历找到使 |ei - ej| 最大的第二个 alpha
            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:
            # 随机选择第二个 alpha
            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 对
        """
        # 步骤 1：计算第一个 alpha 对应的样本点的误差
        ei = SVM.calc_ek(os, i)
        # 如果 (yi * f(xi) < 1 && alpha_i < c) 或者 (yi * f(xi) > 1 && alpha_i > 0)
        # 说明不满足 KKT 条件，可以将其作为第一个 alpha
        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)):
            # 步骤 1：选择第二个 alpha，并计算误差
            j, ej = SVM.select_j(i, os, ei)
            alpha_iold = os.alphas[i].copy()
            alpha_jold = os.alphas[j].copy()
            # 步骤 2：求取 alphas_j 取值的上下限
            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
            # 步骤 3：计算 2*x1*x2 - x1*x2 - x2*x2
            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
            # 步骤 4：求解未修剪的 alpha_2
            os.alphas[j] -= os.label_mat[j] * (ei - ej) / eta
            # 步骤 5：求解修剪后的 alpha_2
            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
            # 步骤 6：求解 alpha_1
            os.alphas[i] += os.label_mat[j] * os.label_mat[i] * \
                            (alpha_jold - os.alphas[j])
            SVM.update_ek(os, i)
            # 步骤 7：求解 b1 和 b2
            b1 = os.b - ei - os.label_mat[i] * (os.alphas[i] - alpha_iold) * \
                os.x[i, :] * os.x[i, :].T - os.label_mat[j] * \
                 (os.alphas[j] - alpha_jold) * os.x[i, :] * os.x[j, :].T
            b2 = os.b - ej - os.label_mat[i] * (os.alphas[i] - alpha_iold) * \
                os.x[i, :] * os.x[j, :].T - os.label_mat[j] * \
                 (os.alphas[j] - alpha_jold) * os.x[j, :] * os.x[j, :].T
            # 步骤 8：求解 b
            if (os.alphas[i] > 0) and (os.alphas[i] < os.c):
                os.b = b1
            elif (os.alphas[j] > 0) and (os.alphas[j] < os.c):
                os.b = b2
            else:
                os.b = (b1 + b2) / 2.0
            return 1
        else:
            return 0

    def smo_p(self, input_data, label_data, c, toler, max_iter):
        """
        SMO 算法实现
        :param input_data: 样本特征值
        :param label_data: 样本标签值
        :param c: 惩罚因子
        :param toler: 容错率
        :param max_iter: 最大迭代次数
        :return:
        """
        self.b = None
        self.w = None
        # 将需要操作的数据存放到数据结构 OptStruct 中
        os = OptStruct(np.mat(input_data), \
                       np.mat(label_data).transpose(), c, toler)
        iter = 0
        entire_set = True
        # 用于标记是否有更新 alpha 对
        alpha_pairs_changed = 0
        # 当迭代次数超过指定的最大值，
        # 或者遍历整个训练集都未对任意 alpha 对进行更新时，退出循环
        while (iter < max_iter) and ((alpha_pairs_changed > 0) or entire_set):
            alpha_pairs_changed = 0
            if entire_set:
                # 遍历整个训练集
                # 这里的 i 为第一个 alpha 对应的索引值
                for i in range(os.m):
                    alpha_pairs_changed += self.inner_l(i, os)
                print('fullSet, iter: {} i: {},' + \
                      'paris changed {}'.format(iter, i, alpha_pairs_changed))
                iter += 1
            else:
                # 遍历所有 0 < alpha < c 的样本点
                non_bound_i_s = \
                    np.nonzero((os.alphas.A > 0) * (os.alphas.A < os.c))[0]
                # 这里的 i 为第一个 alpha 对应的索引值
                for i in non_bound_i_s:
                    alpha_pairs_changed += self.inner_l(i, os)
                    print('non-bound, iter: {} i: {},' +\
                          'pairs changed {}'.format(iter, i, alpha_pairs_changed))
                iter += 1
            if entire_set:
                entire_set = False
            elif alpha_pairs_changed == 0:
                entire_set = True
            print('iteration number: {}'.format(iter))
        self.b = os.b
        self.alphas = os.alphas
        return os.b, os.alphas

    def calc_ws(self, input_data, label_data):
        """
        计算权值 w
        :param input_data: 样本特征值
        :param label_data: 样本标签值
        :return: 权值w
        """
        x = np.mat(input_data)
        label_mat = np.mat(label_data).transpose()
        m, n = np.shape(x)
        w = np.zeros((n, 1))
        for i in range(m):
            w += np.multiply(self.alphas[i] * label_mat[i], x[i, :].T)
        self.w = w
        # self.alphas = None
        return w

    def predict(self, input_feature):
        """
        进行预测
        :param input_feature: 需要预测的特征值
        :return: 预测标签
        """
        fx = np.mat(input_feature) * self.w + self.b
        output_label = -1
        if fx > 0:
            output_label = 1
        return output_label


def load_data_set(file_name):
    """
    从文件中获取数据集
    :param file_name: 文件名
    :return: 返回从文件中获取的数据集
             input_data 存储特征值，label_data 存储标签值
    """
    input_data, label_data = [], []
    fr = open(file_name)
    for line in fr.readlines():
        cur_line = line.strip().split('\t')
        input_data.append([float(cur_line[0]), float(cur_line[1])])
        label_data.append(float(cur_line[2]))
    return input_data, label_data


def show_plot(input_data, label_data, w, b, alphas):
    # 绘制训练集中的样本点
    data01_x1, data01_x2, data02_x1, data02_x2 = [], [], [], []
    for i in range(len(input_data)):
        if label_data[i] > 0:
            data01_x1.append(input_data[i][0])
            data01_x2.append(input_data[i][1])
        else:
            data02_x1.append(input_data[i][0])
            data02_x2.append(input_data[i][1])
    plt.scatter(data01_x1, data01_x2, s=30)
    plt.scatter(data02_x1, data02_x2, s=30)

    # 绘制分割线
    x_max = max(input_data)[0]
    x_min = min(input_data)[0]
    x = np.arange(x_min, x_max, 0.1)
    y = (-float(b) - float(w[0])*x) / float(w[1])
    plt.plot(x, y)

    # 标记支持向量
    for i, alpha in enumerate(alphas):
        if abs(alpha) > 0:
            x1, x2 = input_data[i]
            plt.scatter([x1], [x2], s=100,\
                        c='none', linewidths=1.5, edgecolors='red')
    plt.show()


def test_svm():
    input_data, label_data = load_data_set('testSet.txt')
    svm = SVM()
    svm.smo_p(input_data, label_data, 0.6, 0.001, 40)
    for i in range(len(svm.alphas)):
        if svm.alphas[i] > 0.0:
            print('support vector: {}, {}; alpha = {}'.\
                  format(input_data[i], label_data[i], svm.alphas[i]))
    print('b = {}'.format(svm.b))
    svm.calc_ws(input_data, label_data)
    print('w = {}'.format(svm.w))
    show_plot(input_data, label_data, svm.w, svm.b, svm.alphas)


if __name__ == "__main__":
    test_svm()

运行以上代码，将打印如下信息（只列出部分打印信息）及绘制如下图形：

support vector: [3.542485, 1.977398], -1.0; alpha = [[0.06961952]]
support vector: [2.114999, -0.004466], -1.0; alpha = [[0.0169055]]
support vector: [8.127113, 1.274372], 1.0; alpha = [[0.0169055]]
support vector: [4.658191, 3.507396], -1.0; alpha = [[0.0272699]]
support vector: [8.197181, 1.545132], 1.0; alpha = [[0.04522972]]
support vector: [7.40786, -0.121961], 1.0; alpha = [[0.0272699]]
support vector: [6.960661, -0.245353], 1.0; alpha = [[0.0243898]]
support vector: [6.080573, 0.418886], 1.0; alpha = [[0.06140181]]
support vector: [3.107511, 0.758367], -1.0; alpha = [[0.06140181]]
b = [[-2.89901748]]
w = [[ 0.65307162]
 [-0.17196128]]

图-2

参考：
[1] 周志华《机器学习》
[2] 李航《统计学习方法》
[3] 《机器学习实战》
[4] https://blog.youkuaiyun.com/luoshixian099/article/details/51227754
[5] https://zhuanlan.zhihu.com/p/29604517
[6] https://zhuanlan.zhihu.com/p/29872905