0 写在前面

机器学习强基计划聚焦深度和广度，加深对机器学习模型的理解与应用。“深”在详细推导算法模型背后的数学原理；“广”在分析多个机器学习模型：决策树、支持向量机、贝叶斯与马尔科夫决策、强化学习等。

🚀详情：机器学习强基计划(附几十种经典模型源码合集)

1 为什么需要SMO算法？

详细推导支持向量机SVM原理+Python实现中列出的SVM对偶问题的求解是二次规划问题，可使用二次规划算法进行数值解，但解的复杂度正比于拉格朗日乘子 的维度，造成很大的训练开销。

序列最小优化(Sequential Minimal Optimization)算法是结合SVM算法实际提出的高效优化方法，可以将支持向量机的训练速度提升一个量级。

SMO的优化目标就是上一篇文章介绍的软间隔SVM的优化目标，约束中略去与算法无关的KKT条件重新写在下方

$$\{\begin{array}{l}\underset{\mathbf{\alpha }}{\mathrm{argmax}}{\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^T{\mathbf{x}}_j\\ \begin{array}{r}{\sum }_{i=1}^m{\alpha }_i{y}_i=0\\ 0⩽{\alpha }_i⩽C\end{array}\}\left(i=1,2,\cdots ,m\right),KKT\text{约束}\end{array}$$

注意： 本文理论难度与编程难度较大，请耐心阅读，仔细思考。

2 优化变量的选择

考虑到约束${\sum }_{i=1}^m{\alpha }_i{y}_i=0$，故每次迭代至少更新一对变量$\left({\alpha }_l,{\alpha }_t\right)$，否则会破坏约束条件，设$\left({\alpha }_l,{\alpha }_t\right)$满足：

$$\{\begin{array}{l}{\alpha }_l{y}_l+{\alpha }_t{y}_t=k\\ {\alpha }_{l,u}{y}_l+{\alpha }_{t,u}{y}_t=k\\ {\sum }_{i\ne l,t}^m{\alpha }_i{y}_i=-k\end{array}$$

其中$\left({\alpha }_l,{\alpha }_t\right)$是上一轮迭代的参数值，是常量；$\left({\alpha }_{l,u},{\alpha }_{t,u}\right)$是本轮优化变量，称为未修剪参数，因为其不一定满足KKT约束，需要根据可行域进一步修剪；将$\left({\alpha }_{l,u},{\alpha }_{t,u}\right)$修剪为符合约束的$\left({\alpha }_{l,\ast },{\alpha }_{t,\ast }\right)$即为所求。

3 优化目标的约简

只考虑这一对待优化变量$\left({\alpha }_{l,u},{\alpha }_{t,u}\right)$，其他参数固定为常数，此时优化目标变为

$$\begin{gathered} \Gamma \left({\alpha }_{l,u},{\alpha }_{t,u}\right)=\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^T{\mathbf{x}}_j \\ =\left({\alpha }_{l,u}+{\alpha }_{t,u}\right)+\sum _{i\ne l,t}^m{\alpha }_i-\frac{1}{2}{\alpha }_{l,u}^2{\mathbf{x}}_l^T{\mathbf{x}}_l-\frac{1}{2}{\alpha }_{t,u}^2{\mathbf{x}}_t^T{\mathbf{x}}_t \\ -{\alpha }_{l,u}{\alpha }_{t,u}{y}_l{y}_t{\mathbf{x}}_l^T{\mathbf{x}}_t-{\alpha }_{l,u}{y}_l{\mathbf{x}}_l^T\sum _{i\ne l,t}^m{\alpha }_i{y}_i{\mathbf{x}}_i-{\alpha }_{t,u}{y}_t{\mathbf{x}}_t^T\sum _{i\ne l,t}^m{\alpha }_i{y}_i{\mathbf{x}}_i \end{gathered}$$

记$v_l={\mathbf{x}}_l^T{\sum }_{i\ne l,t}^m{\alpha }_i{y}_i{\mathbf{x}}_i$，$v_t={\mathbf{x}}_t^T{\sum }_{i\ne l,t}^m{\alpha }_i{y}_i{\mathbf{x}}_i$，略去常数项${\sum }_{i\ne l,t}^m{\alpha }_i$，并代入${\alpha }_{l,u}=\left(k-{\alpha }_{t,u}{y}_t\right)y_l$可得

$$\begin{gathered} \Gamma \left({\alpha }_{t,u}\right)=\left(k-{\alpha }_{t,u}{y}_t\right)y_l+{\alpha }_{t,u}-\frac{1}{2}{\left(k-{\alpha }_{t,u}{y}_t\right)}^2{\mathbf{x}}_l^T{\mathbf{x}}_l-\frac{1}{2}{\alpha }_{t,u}^2{\mathbf{x}}_t^T{\mathbf{x}}_t \\ -\left(k-{\alpha }_{t,u}{y}_t\right){\alpha }_{t,u}{y}_t{\mathbf{x}}_l^T{\mathbf{x}}_t-\left(k-{\alpha }_{t,u}{y}_t\right)v_l-{\alpha }_{t,u}{y}_t{v}_t \end{gathered}$$

化为单变量二次规划问题，可直接求导获得极值。SMO算法正是通过将一组参数的组合优化问题分解为一系列变量的单优化问题获得高效性。

$$\begin{gathered} \frac{\partial \Gamma \left({\alpha }_{t,u}\right)}{\partial {\alpha }_{t,u}}=-y_t{y}_l+1+y_t\left(k-{\alpha }_{t,u}{y}_t\right){\mathbf{x}}_l^T{\mathbf{x}}_l-{\alpha }_{t,u}{\mathbf{x}}_t^T{\mathbf{x}}_t-k{y}_t{\mathbf{x}}_l^T{\mathbf{x}}_t+2{\alpha }_{t,u}{\mathbf{x}}_l^T{\mathbf{x}}_t+y_t{v}_l-y_t{v}_t \\ =0 \end{gathered}$$

代入上一轮迭代常数$\left({\alpha }_l,{\alpha }_t\right)$消去$k$

$$\begin{gathered} \left({\mathbf{x}}_l^T{\mathbf{x}}_l+{\mathbf{x}}_t^T{\mathbf{x}}_t-2{\mathbf{x}}_l^T{\mathbf{x}}_t\right){\alpha }_{t,u}=-y_t{y}_l+1+y_t{v}_l-y_t{v}_t+y_t\left({\mathbf{x}}_l^T{\mathbf{x}}_l-{\mathbf{x}}_l^T{\mathbf{x}}_t\right)k \\ =-y_t{y}_l+1+y_t{v}_l-y_t{v}_t+y_t\left({\mathbf{x}}_l^T{\mathbf{x}}_l-{\mathbf{x}}_l^T{\mathbf{x}}_t\right)\left({\alpha }_l{y}_l+{\alpha }_t{y}_t\right) \\ =\left({\mathbf{x}}_l^T{\mathbf{x}}_l+{\mathbf{x}}_t^T{\mathbf{x}}_t-2{\mathbf{x}}_l^T{\mathbf{x}}_t\right){\alpha }_t \\ +y_t\left[{\alpha }_l{y}_l{\mathbf{x}}_l^T{\mathbf{x}}_l+{\alpha }_t{y}_t{\mathbf{x}}_l^T{\mathbf{x}}_t+v_l+b-y_l\right] \\ -y_t\left[{\alpha }_l{y}_l{\mathbf{x}}_l^T{\mathbf{x}}_t+{\alpha }_t{y}_t{\mathbf{x}}_t^T{\mathbf{x}}_t+v_t+b-y_t\right] \\ =\left({\mathbf{x}}_l^T{\mathbf{x}}_l+{\mathbf{x}}_t^T{\mathbf{x}}_t-2{\mathbf{x}}_l^T{\mathbf{x}}_t\right){\alpha }_t \\ +y_t\left[\sum _{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i^T{\mathbf{x}}_l+b-y_l\right]-y_t\left[\sum _{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i^T{\mathbf{x}}_t+b-y_t\right] \\ \stackrel{\mathbf{w}={\sum }_{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i}{=}\left({\mathbf{x}}_l^T{\mathbf{x}}_l+{\mathbf{x}}_t^T{\mathbf{x}}_t-2{\mathbf{x}}_l^T{\mathbf{x}}_t\right){\alpha }_t+y_t\left({\hat{y}}_l-y_l\right)-y_t\left({\hat{y}}_t-y_t\right) \end{gathered}$$

设样本预测类别与真实类别的误差项为

$$\{\begin{array}{l}{E}_l={\mathbf{w}}^T{\mathbf{x}}_l+b-y_l\\ E_t={\mathbf{w}}^T{\mathbf{x}}_t+b-y_t\end{array}$$

即得

$${\alpha }_{t,u}={\alpha }_t+\frac{y_t\left(E_l-E_t\right)}{\left({\mathbf{x}}_l^T{\mathbf{x}}_l+{\mathbf{x}}_t^T{\mathbf{x}}_t-2{\mathbf{x}}_l^T{\mathbf{x}}_t\right)}$$

4 参数可行性修剪

第3节求导计算得到的${\alpha }_{t,u}$有可能不符合约束条件

$${\alpha }_{l,u}{y}_l+{\alpha }_{t,u}{y}_t=k$$

始终需要记住，我们能进行单变量求解，一定是在满足上面的约束的前提下，否则就会破坏KKT条件，解得的结果也不是最优。所以我们需要对不符合约束的变量进行修剪，修剪至可行域

对于二分类情形，无非就两种情况

• $y_l=y_t$，则${\alpha }_{l,u}+{\alpha }_{t,u}=k$
• $y_l\ne y_t$，则${\alpha }_{l,u}-{\alpha }_{t,u}=k$

可视化如图所示

这张图是什么意思呢？由于优化参数的定义域是$[0,C]$，所以它们被限制在以$C$为边长的方形中。以$y_l\ne y_t$的情形为例，实际上就是直线${\alpha }_{t,u}={\alpha }_{l,u}+k$在方形中移动的过程。考虑$k$的正负就有图中红、蓝两种情况，所以${\alpha }_{t,u}$在$y_l\ne y_t$时的取值范围是

$$\{\begin{array}{l}L=\max \left(0,{\alpha }_t-{\alpha }_l\right)\\ H=\min \left(C,C+{\alpha }_t-{\alpha }_l\right)\end{array}$$

另一种情况同理，最终我们可以得到

$${\alpha }_{t,\ast }=\{\begin{array}{l}H,{\alpha }_{t,u}>H\\ {\alpha }_{t,u},L⩽{\alpha }_{t,u}⩽H\\ L,{\alpha }_{t,u}<L\end{array},\text{其中}\begin{array}{c}L=\{\begin{array}{l}\max \left(0,{\alpha }_t-{\alpha }_l\right),y_l\ne y_t\\ \max \left(0,{\alpha }_t+{\alpha }_l-C\right),y_l=y_t\end{array}\\ H=\{\begin{array}{l}\min \left(C,C+{\alpha }_t-{\alpha }_l\right),y_l\ne y_t\\ \min \left(C,{\alpha }_t+{\alpha }_l\right),y_l=y_t\end{array}\end{array}$$

5 权重与偏置更新

接着更新超平面参数$\mathbf{w}={\sum }_{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i^T$与$b$以供下一轮迭代使用。

对于$b$的更新，若$0<{\alpha }_{l,\ast }<C$或$0<{\alpha }_{t,\ast }<C$则样本位于超平面上

$$y_l\left({\mathbf{w}}_{\ast }^T\mathbf{x}+b_{l,\ast }\right)=1\Rightarrow b_{l,\ast }=y_l-{\sum }_{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i^T{\mathbf{x}}_l$$

提出更新过的优化参数可得到迭代式

$$\begin{gathered} b_{l,\ast }=y_l-\sum _{i\ne l,t}{\alpha }_i{y}_i{\mathbf{x}}_i^T{\mathbf{x}}_l-{\alpha }_{l,\ast }{y}_l{\mathbf{x}}_l^T{\mathbf{x}}_l-{\alpha }_{t,\ast }{y}_t{\mathbf{x}}_t^T{\mathbf{x}}_l \\ =y_l-\sum _i{\alpha }_i{y}_i{\mathbf{x}}_i^T{\mathbf{x}}_l-b+b+{\alpha }_l{y}_l{\mathbf{x}}_l^T{\mathbf{x}}_l+{\alpha }_t{y}_t{\mathbf{x}}_t^T{\mathbf{x}}_l-{\alpha }_{l,\ast }{y}_l{\mathbf{x}}_l^T{\mathbf{x}}_l-{\alpha }_{t,\ast }{y}_t{\mathbf{x}}_t^T{\mathbf{x}}_l \\ =b-E_l+\left({\alpha }_l-{\alpha }_{l,\ast }\right)y_l{\mathbf{x}}_l^T{\mathbf{x}}_l+\left({\alpha }_t-{\alpha }_{t,\ast }\right)y_t{\mathbf{x}}_t^T{\mathbf{x}}_l \end{gathered}$$

$b_{t,\ast }$的计算同理，综合为

$$\{\begin{array}{l}{b}_{\ast }=b_{l,\ast },0<{\alpha }_{l,\ast }<C\\ b_{\ast }=b_{t,\ast },0<{\alpha }_{t,\ast }<C\\ b_{\ast }=\left(b_{l,\ast }+b_{t,\ast }\right)/2,\mathrm{otherwise}\end{array}$$

6 收敛性分析

SMO算法的收敛性由Osuna定理保证

若将一个大规模二次规划问题分解为一系列小规模二次规划子问题，且子问题总是加入至少一个违反KKT条件的变量，那么原二次规划问题求解一定收敛。

从算法收敛的角度SMO算法设置了内外两层循环，外循环挑选违反KKT条件的第一个变量${\alpha }_t$，怎么挑选呢？首先我们知道不违反KKT条件的变量满足

$$\{\begin{array}{l}{\alpha }_i=0\Rightarrow y_i\left({\mathbf{w}}^T{\mathbf{x}}_i+b\right)>1\\ 0<{\alpha }_i<C\Rightarrow y_i\left({\mathbf{w}}^T{\mathbf{x}}_i+b\right)=1\\ {\alpha }_i=C\Rightarrow y_i\left({\mathbf{w}}^T{\mathbf{x}}_i+b\right)<1\end{array}$$

在工程上引入精度常数避免浮点误差

$$\{\begin{array}{l}{\alpha }_i=0\Rightarrow y_i\left({\mathbf{w}}^T{\mathbf{x}}_i+b\right)⩾1-\epsilon \\ 0<{\alpha }_i<C\Rightarrow 1-\epsilon ⩽y_i\left({\mathbf{w}}^T{\mathbf{x}}_i+b\right)⩽1+\epsilon \\ {\alpha }_i=C\Rightarrow y_i\left({\mathbf{w}}^T{\mathbf{x}}_i+b\right)⩽1+\epsilon \end{array}\Rightarrow \{\begin{array}{l}{\alpha }_i=0\Rightarrow y_i{E}_i⩾-\epsilon \\ 0<{\alpha }_i<C\Rightarrow -\epsilon ⩽y_i{E}_i⩽\epsilon \\ {\alpha }_i=C\Rightarrow y_i{E}_i⩽\epsilon \end{array}$$

所以不符合KKT约束的变量满足

$$\mathbf{violator}:\left({\alpha }_i<C\mathbf{and}{y}_i{E}_i<-\epsilon \right)\mathbf{or}\left({\alpha }_i>0\mathbf{and}{y}_i{E}_i>\epsilon \right)$$

第二个变量${\alpha }_l$我们可以启发式地选择与${\alpha }_t$相差最远的，这样更新效率最高，即

$$max|E_l-E_t|$$

7 Python实现

7.1 整体算法流程

7.2 挑选优化变量

外层循环选择不满足KKT条件的变量

 for t in range(self.m):
     Et = self.E[t]
     # 判断是否符合KKT条件
     if (self.alpha[t] < self.C and self.y[t] * Et < -self.tol) or \
        (self.alpha[t] > 0 and self.y[t] * Et > self.tol):
         # 启发式选择参数l
         l = self.__selectL(t)
         if l >= 0 and self.updateAlpha(t, l):
             changed = True
     # 符合KKT条件则跳过
     else: continue

内层循环选择满足$max|E_l-E_t|$条件的变量

def __selectL(self, t):
    validIndex = self.alpha.nonzero()[0]
    if validIndex.size > 0:
        deltaE = np.abs(self.E[t] - self.E[validIndex])
        return validIndex[np.argmax(deltaE)]
    else:
        return -1

7.3 裁剪并更新alpha

def updateAlpha(self, t, l):
    # 保留旧值
    alpha_t, alpha_l = self.alpha[t].copy(), self.alpha[l].copy()
    # 计算学习率
    eta = np.dot(self.X[:, l].T, self.X[:, l]) + np.dot(self.X[:, t].T, self.X[:, t]) \
            - 2 * np.dot(self.X[:, l].T, self.X[:, t])
    # 保持其中一个变量的更新方向
    if eta <= 0: 
        return False
    # 计算未裁剪参数
    alpha_tu = self.alpha[t] + self.y[t] * (self.E[l] - self.E[t]) / eta
    # 裁剪取值上下界
    if (self.y[t] != self.y[l]):
        L = max(0, self.alpha[t] - self.alpha[l])
        H = min(self.C, self.C + self.alpha[t] - self.alpha[l])
    else:
        L = max(0, self.alpha[t] + self.alpha[l] - self.C)
        H = min(self.C, self.alpha[t] + self.alpha[l])
    if L==H:
        return False
    # 裁剪
    if alpha_tu < L:
        self.alpha[t] = L
    elif alpha_tu > H:
        self.alpha[t] = H
    else:
        self.alpha[t] = alpha_tu
    self.alpha[l] = self.alpha[l] + self.y[l] * self.y[t] * (alpha_t - self.alpha[t])
    # 计算优化参数变化量
    deltaAlphaT, deltaAlphaL = self.alpha[t] - alpha_t, self.alpha[l] - alpha_l
    if abs(deltaAlphaT) < 0.0001 or abs(deltaAlphaL) < 0.0001:
        return False
    # 更新权重、偏置和误差参数
    self.updateW()
    self.updateB(t, l, deltaAlphaT, deltaAlphaL)
    self.updateE()
    return True

7.4 更新权重与偏置

 def updateW(self):
    self.w = np.sum((self.alpha * self.y).T * self.X, axis=1, keepdims=True)

def updateB(self, t:int, l:int, deltaAlphaT:float, deltaAlphaL:float):
    bL = -self.E[l] - self.y[t] * np.dot(self.X[:, l].T, self.X[:, t]) * (deltaAlphaT) \
         - self.y[l] * np.dot(self.X[:, l].T, self.X[:, l]) * (deltaAlphaL) + self.b
    bT = -self.E[t] - self.y[l] * np.dot(self.X[:, l].T, self.X[:, t]) * (deltaAlphaL) \
         - self.y[t] * np.dot(self.X[:, t].T, self.X[:, t]) * (deltaAlphaT) + self.b
    if 0 < self.alpha[t] < self.C:
        self.b = bT
    elif 0 < self.alpha[l] < self.C:
        self.b = bL
    else:
        self.b = (bL + bT)/2

7.5 可视化

完整代码请联系下方博主名片

🔥 更多精彩专栏：

• 《ROS从入门到精通》
• 《机器人原理与技术》
• 《机器学习强基计划》
• 《计算机视觉教程》
• …