SMO算法求解支持向量机的理论推导与代码实现

阅读提示：1.文章中的斜体符号表示单个元素，加粗斜体表示向量。2.本篇文章公式大多较长，手机端需要左右滑动查看。3.文章中只给出了部分关键代码，如需全部代码可访问我的Github获取：https://github.com/luliu-fighting/Statistical-Learning-Method。代码为之前复现李航的《统计学习方法》的代码，注释没改。

1.前言

  在前一篇文章中我们推导出了非线性可分支持向量机最优化问题的数学表达式：
$$\begin{gathered} {}_{\mathbf{\alpha }}^{min}\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_{j}K({\mathbf{x}}_i,{\mathbf{x}}_j)-\sum _{i=1}^N{\alpha }_i \\ s.t.\sum _{i=1}^N{\alpha }_i{y}_i=0 \\ 0\le {\alpha }_i\le C \end{gathered}$$其中变量${\alpha }_i$的个数与样本点$({\mathbf{x}}_i,y_i)$的个数$N$是相等的，也就是求解的复杂度与训练集的样本个数是成正比的，当训练数据集容量很大时，求解过程会变得十分复杂。为了提高支持向量机的求解效率，Platt在1998年提出了序列最小最优化(sequential minimal optimization, SMO)算法。
  SMO算法的基本思路是：将原来的最优化问题分解为若干个子最优化问题，每个子问题都是两个变量的最优化问题，其他变量是固定的可以看作是常量；子问题的求解是比较简单的，因为子问题可以通过解析方法求得具体的值。当所有的变量的解都满足此最优化问题的KKT条件时，那么此时的解就是要求的最优化问题的解。对于SMO算法为什么高效，我的理解是：SMO算法不是每次都对所有的变量进行优化，而是每次按照某种规则选择两个变量进行优化，当所有的解都满足KKT条件时，问题就求解完成了，无论训练样本个数即变量个数有多少，SMO算法每次都只优化更新两个变量。
  SMO算法可以分为两个部分，一个是每次如何选取子问题的两个变量，另一个是选取的两个变量的子最优化问题的求解。下文将先讲子问题的求解，再讲变量的选取。

2.两个变量的子最优化问题的求解

2.1 ${\alpha }_1^{new},{\alpha }_2^{new}$的求解

  假设此时已经根据变量的选择规则选出了变量${\alpha }_1$和${\alpha }_2$，注意这里的索引1并不是指$\mathbf{\alpha }$的第一个分量，而是指被第一个选择出来要进行优化的分量，也就是说${\alpha }_1$可以是$\mathbf{\alpha }$中的任意分量。此时，其他变量可以视为常量，所以有：
$$\begin{gathered} \sum _{i=1}^N{y}_i{\alpha }_i={\alpha }_1{y}_1+{\alpha }_2{y}_2+\sum _{i=3}^N{y}_i{\alpha }_i=0 \\ {\alpha }_1{y}_1+{\alpha }_2{y}_2=-\sum _{i=3}^N{y}_i{\alpha }_i=k \end{gathered}$$其中$k$为常数，于是我们可以写出子最优化问题的表达式：
$$\begin{gathered} {}_{{\alpha }_1,{\alpha }_2}^{min}W({\alpha }_1,{\alpha }_2)=\frac{1}{2}{K}_{11}{\alpha }_1^2+\frac{1}{2}{K}_{22}{\alpha }_2^2+y_1{y}_2{K}_{12}{\alpha }_1{\alpha }_2-({\alpha }_1+{\alpha }_2) \\ +y_1{\alpha }_1\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i1}+y_2{\alpha }_2\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i2}+Constant \\ s.t.{\alpha }_1{y}_1+{\alpha }_2{y}_2=k \\ 0\le {\alpha }_1\le C \\ 0\le {\alpha }_2\le C \end{gathered}$$在目标函数中，$K_{ij}=K({\mathbf{x}}_i,{\mathbf{x}}_j)$，$Constant$代表一些不含${\alpha }_1,{\alpha }_2$的常数项，不影响最优化问题的求解。
  由约束条件${\alpha }_1{y}_1+{\alpha }_2{y}_2=k$及$y_1^2=1$可得：
$${\alpha }_1=y_1(k-{\alpha }_2{y}_2)$$将目标函数中的${\alpha }_1$全部用${\alpha }_2$代替可得：
$$\begin{gathered} {}_{{\alpha }_2}^{min}W({\alpha }_2)=\frac{1}{2}{K}_{11}(k-{\alpha }_2{y}_2{)}^2+\frac{1}{2}{K}_{22}{\alpha }_2^2+y_2{K}_{12}(k-{\alpha }_2{y}_2){\alpha }_2-y_1(k-{\alpha }_2{y}_2)-{\alpha }_2 \\ +(k-{\alpha }_2{y}_2)\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i1}+y_2{\alpha }_2\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i2}+Constant \end{gathered}$$设$g(\mathbf{x})$为对样本$(\mathbf{x},y)$的预测标签值，即分离决策函数
$f(\mathbf{x})$去掉了符号函数输出:
$$g(\mathbf{x})=\sum _{i=1}^N{\alpha }_i{y}_{i}K({\mathbf{x}}_i,\mathbf{x})+b$$由此可以得出：
$$\begin{gathered} \sum _{i=3}^N{y}_i{\alpha }_i{K}_{i1}=g({\mathbf{x}}_1)-y_1{\alpha }_1{K}_{11}-y_2{\alpha }_2{K}_{22}-b \\ \sum _{i=3}^N{y}_i{\alpha }_i{K}_{i2}=g({\mathbf{x}}_2)-y_1{\alpha }_1{K}_{11}-y_2{\alpha }_2{K}_{22}-b \end{gathered}$$对${\alpha }_2$求导数并令其等于0：
$$\begin{gathered} \frac{dW}{d{\alpha }_2}=K_{11}{\alpha }_2+K_{22}{\alpha }_2-2K_{12}{\alpha }_2-K_{11}k{y}_2+K_{12}k{y}_2+y_1{y}_2-1 \\ -y_2\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i1}+y_2\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i2}=0 \\ (K_{11}+K_{22}-2K_{12}){\alpha }_2=y_2(y_2-y_1+k{K}_{11}-k{K}_{12}+\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i1}-\sum _{i=3}^N{y}_i{\alpha }_i{K}_{i2}) \end{gathered}$$将
$$\begin{gathered} k={\alpha }_1^{old}{y}_1+{\alpha }_2^{old}{y}_2 \\ \sum _{i=3}^N{y}_i{\alpha }_i{K}_{i1}=g({\mathbf{x}}_1)-y_1{\alpha }_1^{old}{K}_{11}-y_2{\alpha }_2^{old}{K}_{22}-b \\ \sum _{i=3}^N{y}_i{\alpha }_i{K}_{i2}=g({\mathbf{x}}_2)-y_1{\alpha }_1^{old}{K}_{11}-y_2{\alpha }_2^{old}{K}_{22}-b \end{gathered}$$代入上式：
$$(K_{11}+K_{22}-2K_{12}){\alpha }_2^{new,unc}=(K_{11}+K_{22}-2K_{12}){\alpha }_2^{old}+y_2(g({\mathbf{x}}_1)-y_1-(g({\mathbf{x}}_2)-y_2))$$令：
$$\begin{gathered} \eta =K_{11}+K_{22}-2K_{12} \\ {E}_i=g({\mathbf{x}}_i)-y_i \end{gathered}$$则：
$${\alpha }_2^{new,unc}={\alpha }_2^{old}+\frac{y_2(E_1-E_2)}{\eta }$$此时求得的${\alpha }_2^{new,unc}$仅仅是根据求导并令导数等于0得到的，还并未考虑到最优化问题的约束条件，其有可能不在约束条件的取值范围内，所以要得到最后的${\alpha }_2^{new}$，还应该知道${\alpha }_2$根据约束条件的取值范围。
  根据约束条件${\alpha }_1{y}_1+{\alpha }_2{y}_2=k$和$0\le {\alpha }_2\le C$可知，${\alpha }_2$的取值范围不仅与其本身有关，而且与${\alpha }_1,y_1,y_2$有关。
（1）当$y_1\ne y_2$时，因为有$0\le {\alpha }_1\le C$，所以有：
$${\alpha }_2^{new}={\alpha }_1^{old}-k\in [0-k,C-k]=[{\alpha }_2^{old}-{\alpha }_1^{old},C-{\alpha }_1^{old}+{\alpha }_2^{old}]$$即
$${\alpha }_2^{old}-{\alpha }_1^{old}\le {\alpha }_2^{new}\le C-{\alpha }_1^{old}+{\alpha }_2^{old}$$再结合${\alpha }_2^{new}$自己的约束$0\le {\alpha }_2^{new}\le C$取交集可得：
$$L\le {\alpha }_2^{new}\le H$$其中：
$$\begin{gathered} L=max(0,{\alpha }_2^{old}-{\alpha }_1^{old}) \\ H=min(C,C-{\alpha }_1^{old}+{\alpha }_2^{old}) \end{gathered}$$
（2）当$y_1=y_2$时，有：
$${\alpha }_2^{new}=k-{\alpha }_1^{old}\in [k-C,k]=[{\alpha }_2^{old}+{\alpha }_1^{old}-C,{\alpha }_1^{old}+{\alpha }_2^{old}]$$即
$${\alpha }_2^{old}-{\alpha }_1^{old}\le {\alpha }_2^{new}\le C-{\alpha }_1^{old}+{\alpha }_2^{old}$$再结合${\alpha }_2^{new}$自己的约束$0\le {\alpha }_2^{new}\le C$取交集可得：
$$L\le {\alpha }_2^{new}\le H$$其中：
$$\begin{gathered} L=max(0,{\alpha }_2^{old}+{\alpha }_1^{old}-C) \\ H=min(C,{\alpha }_1^{old}+{\alpha }_2^{old}) \end{gathered}$$将根据约束条件求得的${\alpha }_2^{new}$的取值范围与前面求导得到的${\alpha }_2^{new}$结合，就可以得到${\alpha }^{new}$的更新公式：
$${\alpha }_2^{new}=\{\begin{array}{ll}H,& \text{if}{\alpha }^{new,unc}>H\\ {\alpha }^{new,unc},& \text{if}L\le {\alpha }^{new,unc}\le H\\ L,& \text{if}{\alpha }^{new,unc}<L\end{array}$$根据公式：
$$y_1{\alpha }_1^{new}+y_2{\alpha }_2^{new}=k=y_1{\alpha }_1^{old}+y_2{\alpha }_2^{old}$$可得${\alpha }_1^{new}$的更新公式为：
$${\alpha }_1^{new}={\alpha }_1^{old}+y_1{y}_2({\alpha }_2^{old}-{\alpha }_2^{new})$$至此，子最优化问题的解$({\alpha }_1^{new},{\alpha }_2^{new})$就得到了。
$({\alpha }_1^{new},{\alpha }_2^{new})$的更新过程代码如下：

            alpha2_new = self.alpha[i2] + self.trainlabel[i2]*(E1 - E2)/eta
            if alpha2_new < L:
                alpha2_new = L
            elif alpha2_new > H:
                alpha2_new = H
            #式（7.109）
            alpha1_new = self.alpha[i1] + self.trainlabel[i1]*self.trainlabel[i2]*(self.alpha[i2] - alpha2_new)

2.2 常数b的求解

  当$0<{\alpha }_1^{new}<C$时，根据KKT条件可知：
$$\sum _{i=1}^N{\alpha }_i{y}_i{K}_{i1}+b=y_1$$所以：
$$b_1^{new}=y_1-\sum _{i=3}^N{\alpha }_i{y}_i{K}_{i1}-{\alpha }_1^{new}{y}_1{K}_{11}-{\alpha }_2^{new}{y}_2{K}_{21}$$由
$$E_1=\sum _{i=3}^N{\alpha }_i{y}_i{K}_{i1}+{\alpha }_1^{old}{y}_1{K}_{11}+{\alpha }_2^{old}{y}_2{K}_{21}+b^{old}-y_1$$可得：
$$y_1-\sum _{i=3}^N{\alpha }_i{y}_i{K}_{i1}=-E_1+{\alpha }_1^{old}{y}_1{K}_{11}+{\alpha }_2^{old}{y}_2{K}_{21}+b^{old}$$于是可得：
$$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}$$
同样，当$ 0 < \alpha^{new}_2 < C$时可以得到：
$$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}$$当$0<{\alpha }_1^{new}<C$并且$0<{\alpha }_2^{new}<C$时，$b_1^{new}$和$b_2^{new}$是相同的，所以$b^{new}=b_1^{new}=b_2^{new}$，但是若${\alpha }_1^{new},{\alpha }_2^{new}$是0或$C$时，二者就不再相等，这时$b^{new}$应取二者的平均值，其实在前一种情况中也是取二者的平均值，因为二者是相同的，所以$b$的更新公式就是：
$$b^{new}=\frac{b_1^{new}+b_2^{new}}{2}$$参数b的更新过程代码如下：

            b1_new = -E1 - self.trainlabel[i1] * self.k[1, 1] * (alpha1_new - self.alpha[i1]) - \
                     self.trainlabel[i2] * self.k[2, 1] * (alpha2_new - self.alpha[i2]) + self.b
            #式（7.116）
            b2_new = -E2 - self.trainlabel[i1] * self.k[1, 2] * (alpha1_new - self.alpha[i1]) - \
                     self.trainlabel[i2] * self.k[2, 2] * (alpha2_new - self.alpha[i2]) + self.b

            if 0 < alpha1_new < self.C:
                b_new = b1_new
            elif 0 < alpha2_new < self.C:
                b_new = b2_new
            else:
                # 选择中点，P148倒数第二段
                b_new = (b1_new + b2_new) / 2

3.变量的选择方法

  SMO算法在每个子问题中选择两个变量进行优化，其中至少一个变量时违反KKT条件的(具体可参考李航《统计学习方法》P128，需要pdf版的可以私聊获取)。

3.1 第一个变量的选择

  选择第一个变量的过程为外层循环，外层循环在训练样本中选取违反KKT条件最严重的样本点，并将其对应的变量作为第一个变量。即检验样本$({\mathbf{x}}_i,y_i)$是否满足KKT条件(由李航《统计学习方法》P128推出来的)：
$$\{\begin{array}{l}{\alpha }_i=0,y_{i}g({\mathbf{x}}_i)\ge 1\\ 0<{\alpha }_i<C,y_{i}g({\mathbf{x}}_i)=1\\ {\alpha }_i=C,y_{i}g({\mathbf{x}}_i)\le 1\end{array}$$其中$g({\mathbf{x}}_i)={\sum }_{j=1}^N{\alpha }_j{y}_{j}K({\mathbf{x}}_i,{\mathbf{x}}_j)+b$。
  在检验过程中，外层循环首先遍历在间隔边界上的样本点，即$0<{\alpha }_i<C$的样本点,因为这些点比较容易违反KKT条件；如果这些点都满足KKT条件，再遍历剩下的点。判断是否符合KKT条件的代码如下：

 #判断样本是否满足KKT条件，用于SMO中第一个变量的选择
    def isSatisfyKKT(self, i):
        y_g = self.calcGxi(i) * self.trainlabel[i]
        #式（7.111）-式（7.113）
        if (self.alpha[i] == 0) and (y_g >= 1):
            return True
        elif (0 < self.alpha[i] < self.C) and (y_g == 1):
            return True
        elif (self.alpha[i] == self.C) and (y_g <= 1):
            return True
        else:
            return False

3.2 第二个变量的选择

  选择第二个变量的过程为内层循环，假设外层循环已经找到第1个变量${\alpha }_1$，内层循环则是要找到第二个变量${\alpha }_2$。第二个变量的选择标准时希望能使${\alpha }_2$有足够大的变化。在前面的推导中，我们可以看到${\alpha }_2^{new}$是依赖于$|E_1-E_2|$的，当$|E_1-E_2|$较大时，${\alpha }_2$的更新速度也会较快，这样就加快了计算速度。所以在选择第二个样本点时，可以使$|E_1-E_2|$的值最大的${\alpha }_2$。因为已经确定了${\alpha }_1$，这样$E_1$的值就确定了。如果$E_1$是正的，那么选择最小的$E_i$作为$E_2$，如果$E_1$是负的，那么选择最大的$E_i$作为$E_2$，这样就能最大化$|E_1-E_2|$了，因为在每次迭代中都需要判断$|E_1-E_2|$的值，所以可以将所有$E_i$值保存在一个列表中，当然，在每次迭代后，因为${\alpha }_1$和${\alpha }_2$都变了，因此也需要对$E_i$列表进行更新。计算$E_i$的代码如下：

    #计算g(xi)
    def calcGxi(self, i):
        #式（7.104）
        gxi = self.b
        for j in range(self.m):
            #在“7.2.3支持向量”开头第一句话有说到“对应于α > 0的样本点
            # (xi, yi)的实例xi称为支持向量”。也就是说只有支持向量的α是大于0的，在求和式内的
            # 对应的αi*yi*K(xi, xj)不为0，非支持向量的αi*yi*K(xi, xj)必为0，也就不需要参与
            # 到计算中。也就是说，在g(xi)内部求和式的运算中，只需要计算α>0的部分，其余部分可
            # 忽略。因为支持向量的数量是比较少的，这样可以再很大程度上节约时间
            # 从另一角度看，抛掉支持向量的概念，如果α为0，αi*yi*K(xi, xj)本身也必为0，从数学
            # 角度上将也可以扔掉不算
            if self.alpha[j] == 0:
                continue
            else:
                gxi += self.alpha[j] * self.trainlabel[j] * self.k[j, i]

        return gxi

    #计算序列最小最优化算法的（SMO）中的E(i)
    def calcEi(self, i):
        #式（7.105）
        return self.calcGxi(i) - self.trainlabel[i]

两个变量的选择过程代码如下：

#SMO中两个变量的选择
    def getAlpha(self):
        #根据P147的中间那段话
        #先找出满足条件0<alpha[i]<C的样本点，检验是否满足KKT条件
        index_list = np.array([i for i in range(self.m) if 0 < self.alpha[i] < self.C])
        #如果index_list中的样本均满足KKT条件，则遍历剩下的样本点，检验是否满足KKT条件
        index_other = np.array([i for i in range(self.m) if i not in index_list])
        index = np.hstack([index_list, index_other])
        for i in index:
            #在Debug过程中发现hstack指令会将int型变为float型，所以要进行调整
            i = int(i)
            if self.isSatisfyKKT(i):
                continue
            #选择样本中违反KKT条件最严重的样本点作为第一个变量
            E1 = self.E[i]
            #根据P147倒数第二段选择第二个变量
            if E1 >= 0:
                j = self.E.index(min(self.E))
            else:
                j = self.E.index(max(self.E))
            return i, j

4.参考资料

1.李航《统计学习方法》
2.https://www.jianshu.com/p/eef51f939ace
3.https://zhuanlan.zhihu.com/p/248862271
4.https://github.com/Dod-o/Statistical-Learning-Method_Code
5.https://github.com/fengdu78/lihang-code