机器学习笔记-SMO序列最小最优化算法

SMO序列最小优化算法

文章目录

• SMO序列最小优化算法
• 前言
• SMO算法介绍
• SMO算法实现
• • 1.变量的选择方法
  • • 1.1 第一个变量的选择
    • 1.2.第二个变量的选择
  • 2.变量的解析方法
  • 3.更新中间变量
• SMO算法的实现

---

前言

  在支持向量机这篇文章中详细叙述了支持向量机的理论推导部分，其中在实现支持向量机的部分，我们说到支持向量机最后就是求解一个二次规划问题，在那里我们使用的是MATLAB中求解二次规划的函数进行求解，这次将会采用一个新的算法来求解支持向量机的二次规划问题。

SMO算法介绍

  SMO算法相比常规求解二次规划问题的算法，是一种非常快速的方法。在支持向量机的对偶问题中，有：

min ⁡ α 1 2 ∑ i = 1 N ∑ j = 1 N α i α j y i y j ( x i ⋅ x j ) − ∑ i = 1 N α i s . t . ∑ i = 1 N α i y i = 0 α i ≥ 0 i = 1 , 2 , ⋯ , N \min\limits_{\alpha}\,\,\,\,\,\frac{1}{2}\sum\limits_{i=1}^N\sum\limits_{j=1}^N{\alpha_i\alpha_jy_iy_j(x_i\cdot x_j)}-\sum\limits_{i=1}^N{\alpha_i}\\ s.t.\,\,\,\,\sum\limits_{i=1}^N{\alpha_iy_i=0}\\ \alpha_i\ge0\,\,\,\,i=1,2,\cdots,N αmin​21​i=1∑N​j=1∑N​αi​αj​yi​yj​(xi​⋅xj​)−i=1∑N​αi​s.t.i=1∑N​αi​yi​=0αi​≥0i=1,2,⋯,N

  在上面的问题中，变量是拉格朗日乘子，虽然只有一种变量，但是变量的个数却有$N$个，且每个变量 α i \alpha_i αi​对应于一个样本点 ( x i , y i ) (x_i,y_i) (xi​,yi​)。SMO是一种启发式算法，它的基本思路就是不断更新迭代使得所有的样本点均满足KKT条件，当所有的点都满足KKT条件时，就找到了此最优化问题的最优解了。具体怎么迭代，是我们要学习的重点，SMO在每次更新时，都会选择两个$\alpha$作为此次迭代的对象，而剩下的$\alpha$均为定值。当确定好需要更新的两个 α 1 \alpha_1 α1​和 α 2 \alpha_2 α2​变量时，SMO算法就会使用一种解析方法来优化两个变量 α 1 \alpha_1 α1​和 α 2 \alpha_2 α2​。所以我们需要掌握的重点就知道了：如何优化 α 1 \alpha_1 α1​, α 2 \alpha_2 α2​两个变量＋如何选取最优的变量。

SMO算法实现

1.变量的选择方法

  在学习SMO算法时，我是对照着《统计学习方法》这本书结合代码来学习的，在这里其实已经是第二遍了，第一遍时我只编出了代码，并没有给出详细的解释，这次打算从0实现一个SMO，顺便介绍我实现代码的思路。
  变量的选择方法是SMO算法最巧妙的地方，在SMO算法中，如果对所有的$N$个变量进行优化，那种速度太慢，SMO采用的方法就是选择两个$\alpha$进行优化，在选择的两个变量中，至少有一个变量是违反了KKT条件的。而这两个变量是一一进行选择的，不是一起把两个的变量都选择，而是一个一个选择，我们把第一个变量的选择称为外层循环，第二个变量的选择称为内层循环。

1.1 第一个变量的选择

  在选择变量之前，我们需要对$N$个$\alpha$赋初值，只有赋完初值后才能选择，初始化所有的$\alpha$均为0，SMO在选择第一个$\alpha$时主要看样本点中违反KKT条件最严重的点，让其作为第一个变量，如何判断样本点违反KKT条件的程度呢？
在KKT条件中有下面的式子：
α i = 0 ⇔ y i g ( x i ) ≥ 1 0 < α i < C ⇔ y i g ( x i ) = 1 α i = C ⇔ y i g ( x i ) ≤ 1 \alpha_i=0\Leftrightarrow y_ig(x_i)\ge1\\ 0<\alpha_iαi​=0⇔yi​g(xi​)≥10<αi​<C⇔yi​g(xi​)=1αi​=C⇔yi​g(xi​)≤1
其中， g ( x i ) = ∑ j = 1 N α j y j K ( x i , x j ) + b g(x_i)=\sum\limits_{j=1}^N\alpha_jy_jK(x_i,x_j)+b g(xi​)=j=1∑N​αj​yj​K(xi​,xj​)+b
  对于一组给定的$\alpha$，我们就可以计算任意样本点的 g ( x i ) g(x_i) g(xi​)。上述的KKT违反条件其实是分为三种情况，例如，当 α i = 0 \alpha_i=0 αi​=0时，如果 y i g ( x i ) < 1 y_ig(x_i)<1 yi​g(xi​)<1那么此时的 α i \alpha_i αi​就是违反KKT条件的点，而我们把 1 − y i g ( x i ) 1-y_ig(x_i) 1−yi​g(xi​)称为违反KKT条件的程度，遍历所有样本点，找出违反KKT条件最大的样本点的 α i \alpha_i αi​作为第一个变量。
  还有一个需要注意的点就是判断某点是否违反KKT条件是在精度$\epsilon$ 下进行的，这是因为计算机在计算优化 α i \alpha_i αi​时不能做到绝对的优化，只能保证优化过程在无限接近最优解，所以我们需要取一个精度范围： α i = 0 ⇔ y i g ( x i ) ≥ 1 − ε 0 < α i < C ⇔ y i g ( x i ) ∈ [ 1 − ε , 1 + ε ] α i = C ⇔ y i g ( x i ) ≤ 1 + ε \alpha_i=0\Leftrightarrow y_ig(x_i)\ge1-\varepsilon\\ 0<\alpha_iαi​=0⇔yi​g(xi​)≥1−ε0<αi​<C⇔yi​g(xi​)∈[1−ε,1+ε]αi​=C⇔yi​g(xi​)≤1+ε
  而精度$\epsilon$是我们开始赋值的，我取了0.01。在遍历所有样本时，也不是采用按顺序遍历，这里首先遍历所有满足条件 0 < α i < C 0<\alpha_i0<αi​<C的样本点，因为这些样本点肯定在间隔边界上对支持向量机的学习起决定性作用，首先检验它们是否满足KKT条件。如果支持向量都满足KKT条件，那么再遍历剩下的所有样本点，检验它们是否满足KKT条件，如果不满足，则取出作为第一个变量 α 1 \alpha_1 α1​。

1.2.第二个变量的选择

  当第一个变量选取之后，我们就可以选取第二个变量了，如果说第一个变量选取的标准是违反KKT条件最大的点，那么第二个变量选取的标准就是能够使 α 2 \alpha_2 α2​有很大的变化，所谓$\alpha$变化最大，就是指在选择 α 2 \alpha_2 α2​之后进行优化，我们新的 α 2 n e w \alpha_2^{new} α2new​改变幅度要大。至于如何判断$\alpha$的变化幅度有多大，需要了解新的 α 2 n e w \alpha^{new}_2 α2new​是如何产生的，由于变量的解析方法在下一节介绍，这里只简单引入一下，那就是通过比较 ∣ E 1 − E 2 ∣ |E_1-E_2| ∣E1​−E2​∣的值来判断改变的幅度大小，而$E$又是什么呢？$E$的计算公式为：
E i = g ( x i ) − y i E_i=g(x_i)-y_i Ei​=g(xi​)−yi​
  所以可以看 E i E_i Ei​的值应该可以说是已知的，是可以计算出来的。而第一个变量已经确定下来，所以 E 1 E_1 E1​的值固定，于是只需要比较所有样本点的 E 2 E_2 E2​的大小即可。
  到这里，我们就已经找到第二个变量的选取标准了：找到最大的$\Delta E$。但是问题又来了，$E$的值有正有负，当 E 1 E_1 E1​的值为正时，我们取$E$中最小的值对应的$\alpha$作为第二个变量，当 E 1 E_1 E1​的值为负时，我们取$E$中最大的值对应的$\alpha$作为第二个变量。另外，如果所有的样本点$\alpha$已知，那么其实$E$的值也是已知的，于是我们可以在每一次迭代之后都计算一遍$E$的值，然后将$E$的值存储在一个列表中，随时调用即可。
  和第一个变量选择一样，第二个变量选择的方法也不是绝对的，因为如果我们通过两层遍历选择了一组变量 α 1 、 α 2 \alpha_1、\alpha_2 α1​、α2​，但是此时的目标函数并没有足够的下降，这时我们就舍去这组变量，重新选择 α 2 \alpha_2 α2​，这个时候选择 α 2 \alpha_2 α2​顺序和选择 α 1 \alpha_1 α1​是一样的，首先遍历间隔边界上的点，然后再遍历其它点，直到使得目标函数有足够的下降。若已经遍历完所有的样本点仍然找不到合适的变量，那就放弃第一个变量，重新再选取第一个变量，且此时要记住选取过的第一个变量不能再选取了。
整个过程中外层遍历和内层遍历选择的点是不重复的。

2.变量的解析方法

  当两个变量都选取之后，我们就可以开始更新这两个变量 α 1 , α 2 \alpha_1,\alpha_2 α1​,α2​了，具体更新的方法也叫做解析方法。首先在更新前我们需要注意，对于选择的两个变量 α 1 , α 2 \alpha_1,\alpha_2 α1​,α2​如果改变必须两个都要变，因为这两组变量之间是相互制约的，如下：
α 1 = − y 1 ∑ i = 2 N α i y i \alpha_1=-y_1\sum\limits_{i=2}^N{\alpha_iy_i} α1​=−y1​i=2∑N​αi​yi​
如果说 α 1 \alpha_1 α1​在变化， y 1 y_1 y1​是不会改变的，所以 α 2 \alpha_2 α2​必定改变，如果 α 2 \alpha_2 α2​不变，那么就不满足最优化问题的约束条件，因为整个过程必须满足：
∑ i = 1 N α i y i = 0 \sum\limits_{i=1}^N{\alpha_iy_i=0} i=1∑N​αi​yi​=0
所以 α 1 \alpha_1 α1​与 α 2 \alpha_2 α2​是一起改变的。
  假设选取的两个变量分别为 α 1 , α 2 \alpha_1,\alpha_2 α1​,α2​，其他的变量 α i ( i = 3 , 4 , ⋯ , N ) \alpha_i(i=3,4,\cdots,N) αi​(i=3,4,⋯,N)是固定不变的。于是我们就可以开始推导：
将对偶问题拆解为含 α 1 、 α 2 \alpha_1、\alpha_2 α1​、α2​的子问题：
min ⁡ α 1 , α 2 W ( α 1 , α 2 ) = 1 2 K 11 α 1 2 + 1 2 K 22 α 2 2 + y 1 y 2 K 12 α 1 α 2 − ( α 1 + α 2 ) + y 1 α 1 ∑ i = 3 N y i α i K i 1 + y 2 α 2 ∑ i = 3 N y i α i K i 2 s . t . α 1 y 1 + α 2 y 2 = − ∑ i = 3 N y i α i = k 0 ≤ α i ≤ C , i = 1 , 2 \min\limits_{\alpha_1,\alpha_2}\,\,\,\,\,\,\,W(\alpha_1,\alpha_2)=\frac{1}{2}K_{11}\alpha_1^2+\frac{1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1\alpha_2-\\(\alpha_1+\alpha_2)+y_1\alpha_1\sum\limits_{i=3}^N{y_i\alpha_iK_{i1}}+y_2\alpha_2\sum\limits_{i=3}^Ny_i\alpha_iK_{i2}\\ s.t.\,\,\,\,\,\alpha_1y_1+\alpha_2y_2=-\sum\limits_{i=3}^N{y_i\alpha_i}=k\\ 0\le\alpha_i\le C,\,\,\,\,\,i=1,2 α1​,α2​min​W(α1​,α2​)=21​K11​α12​+21​K22​α22​+y1​y2​K12​α1​α2​−(α1​+α2​)+y1​α1​i=3∑N​yi​αi​Ki1​+y2​α2​i=3∑N​yi​αi​Ki2​s.t.α1​y1​+α2​y2​=−i=3∑N​yi​αi​=k0≤αi​≤C,i=1,2
  在子问题中我们可以发现，此时的约束变量只有 α 1 \alpha_1 α1​和 α 2 \alpha_2 α2​这两个。其中$K$是核函数。如果想求解上述问题的最优 α 1 , α 2 \alpha_1,\alpha_2 α1​,α2​，需要针对约束条件进行分析，我们可以发现约束条件非常有特点，两个$\alpha$的最大值最小值均相等，且 α 1 \alpha_1 α1​和 α 2 \alpha_2 α2​满足线性关系，这是因为 α 1 y 1 + α 2 y 2 = k \alpha_1y_1+\alpha_2y_2=k α1​y1​+α2​y2​=k等于定值，而 y 1 y_1 y1​和 y 2 y_2 y2​的取值是1或-1，这时的$\alpha$值肯定时存在某种关系的。
y 1 y_1 y1​与 y 2 y_2 y2​的组合有两种情况，要么相等，要么不相等，此时可以将 α 1 \alpha_1 α1​与 α 2 \alpha_2 α2​的值使用二维图来表现出来：
  至于为什么 α 1 \alpha_1 α1​和 α 2 \alpha_2 α2​的值可以表示为上述关系，这是因为 y 1 , y 2 y_1,y_2 y1​,y2​的取值决定的，当 y 1 = y 2 y_1=y_2 y1​=y2​时我们可以对 α 1 y 1 + α 2 y 2 = − ∑ i = 3 N y i α i = k \alpha_1y_1+\alpha_2y_2=-\sum\limits_{i=3}^N{y_i\alpha_i}=k α1​y1​+α2​y2​=−i=3∑N​yi​αi​=k约去$y$参数。结果就变成了 α 1 + α 2 = k \alpha_1+\alpha_2=k α1​+α2​=k这种形式，所以就可以在二维坐标系中画出了，对于 y 1 ≠ y 2 y_1\ne y_2 y1​​=y2​的情况也是如此，我们想更新 α 1 \alpha_1 α1​和 α 2 \alpha_2 α2​的值，就可以先更新 α 2 \alpha_2 α2​然后再更新 α 1 \alpha_1 α1​的值。
下面先讨论 y 1 ≠ y 2 y_1\ne y_2 y1​​=y2​的情况：
当 y 1 ≠ y 2 y_1\ne y_2 y1​​=y2​时，约束条件变成： α 1 − α 2 = k \alpha_1-\alpha_2=k α1​−α2​=k，由于这里没有给出 α 1 \alpha_1 α1​与 α 2 \alpha_2 α2​的大小并没有给出，所以$k$的值有正有负，我们在$k$取正取负时均画在图中：

可以发现在 y 1 ≠ y 2 y_1\ne y_2 y1​​=y2​时， α 2 \alpha_2 α2​的最大值和最小值是不一定的，因为$k$的取值不同。
但是我们可以设新的 α 2 \alpha_2 α2​为 α 2 n e w \alpha_2^{new} α2new​，此时变得到最优值 α 2 n e w \alpha_2^{new} α2new​的取值范围为：
L ≤ α 2 n e w ≤ H L\le\alpha_2^{new}\le H L≤α2new​≤H
其中$L$与$H$为 α 2 n e w \alpha_2^{new} α2new​所在的对角线段端点的界，且一定有：
L = max ⁡ ( 0 , α 2 o l d − α 1 o l d ) , H = min ⁡ ( C , C + α 2 o l d − α 1 o l d ) L=\max(0,\alpha_2^{old}-\alpha_1^{old}),\,\,\,\,\,H=\min(C,C+\alpha_2^{old}-\alpha_1^{old}) L=max(0,α2old​−α1old​),H=min(C,C+α2old​−α1old​)
当 y 1 = y 2 y_1 =y_2 y1​=y2​时，则：
L = max ⁡ ( 0 , α 2 o l d + α 1 o l d − C ) , H = min ⁡ ( C , α 2 o l d + α 1 o l d ) L=\max(0,\alpha_2^{old}+\alpha_1^{old}-C),\,\,\,\,\,H=\min(C,\alpha_2^{old}+\alpha_1^{old}) L=max(0,α2old​+α1old​−C),H=min(C,α2old​+α1old​)
$L$和$H$只是我们对 α 2 n e w \alpha_2^{new} α2new​的值进行修正作用的，就是当新的 α 2 n e w \alpha_2^{new} α2new​超出范围之后，我们要对其进行修正，但是更新 α 2 n e w \alpha_2^{new} α2new​的方法却不是借助这个。更新 α 2 \alpha_2 α2​的方法如下：
记：
g ( x i ) = ∑ j = 1 N α j y j K ( x i , x j ) + b g(x_i)=\sum\limits_{j=1}^N\alpha_jy_jK(x_i,x_j)+b g(xi​)=j=1∑N​αj​yj​K(xi​,xj​)+b
令

E i = g ( x i ) − y i = ( ∑ j = 1 N α j y j K ( x i , x j ) + b ) − y i , i = 1 , 2 E_i=g(x_i)-y_i=(\sum\limits_{j=1}^N\alpha_jy_jK(x_i,x_j)+b)-y_i,\,\,\,\,i=1,2 Ei​=g(xi​)−yi​=(j=1∑N​αj​yj​K(xi​,xj​)+b)−yi​,i=1,2
当$i=1,2$时， E i E_i Ei​为函数$g(x)$对输入 x i x_i xi​的预测值与真实值输出 y i y_i yi​之差。
下面直接给出更新 α 2 \alpha_2 α2​的定理
α 2 n e w , u n c = α 2 o l d + y 2 ( E 1 − E 2 ) η \alpha_2^{new,unc}=\alpha_2^{old}+\frac{y_2(E_1-E_2)}{\eta} α2new,unc​=α2old​+ηy2​(E1​−E2​)​
其中：
η = K 11 + K 22 − 2 K 12 = ∣ ∣ ϕ ( x 1 ) − ϕ ( x 2 ) ∣ ∣ 2 \eta=K_{11}+K_{22}-2K_{12}=||\phi(x_1)-\phi(x_2)||^2 η=K11​+K22​−2K12​=∣∣ϕ(x1​)−ϕ(x2​)∣∣2
其中$\varphi (x)$表示映射关系，但是在核函数中不需要计算映射，直接计算$K(x,x)$即可。
使用上述公式求解得到的 α 2 n e w , u n c \alpha_2^{new,unc} α2new,unc​是没有经过修正的解，后续我们要对其进行修正：
α 2 n e w = { H , α 2 n e w , u n c > H α 2 n e w , u n c L ≤ α 2 n e w , u n c ≤ H L α 2 n e w , u n c < L \alpha_2^{new}=\left\{ \begin{array}{l} H,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\alpha_2^{new,unc}>H\\ \alpha_2^{new,unc}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,L\le\alpha_2^{new,unc}\le H\\ L\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\alpha_2^{new,unc}α2new​=⎩⎨⎧​H,α2new,unc​>Hα2new,unc​L≤α2new,unc​≤HLα2new,unc​<L​
上述式是用来修正 α 2 n e w , u n c \alpha_2^{new,unc} α2new,unc​为 α 2 n e w \alpha_2^{new} α2new​的，在求解得到 α 2 n e w \alpha_2^{new} α2new​之后，就可以利用 α 1 o l d y 1 + α 2 o l d y 2 = α 1 n e w y 1 + α 2 n e w y 2 \alpha_1^{old}y_1+\alpha_2^{old}y_2=\alpha_1^{new}y_1+\alpha_2^{new}y_2 α1old​y1​+α2old​y2​=α1new​y1​+α2new​y2​式子来求解 α 1 n e w \alpha_1^{new} α1new​:
α 1 n e w = α 1 o l d + y 1 y 2 ( α 2 o l d − α 2 n e w ) \alpha_1^{new}=\alpha_1^{old}+y_1y_2(\alpha_2^{old}-\alpha_2^{new}) α1new​=α1old​+y1​y2​(α2old​−α2new​)
到这里我们就可以把 α 1 n e w \alpha_1^{new} α1new​和 α 2 n e w \alpha_2^{new} α2new​全部求出了。

3.更新中间变量

  其实SMO的主体部分已经学习完了，包括新变量的选择和变量的解析方法。在完成一次求解子问题后，我们还需要对一些中间变量进行更新，例如误差$E$，阈值$b$等。
首先来看，在满足KKT条件的情况下有 y i g ( x i ) = 1 y_ig(x_i)=1 yi​g(xi​)=1，两边同除以 y i y_i yi​得到：
∑ i = 1 N α i y i K i , 1 + b = y 1 \sum\limits_{i=1}^N{\alpha_iy_iK_{i,1}+b=y_1} i=1∑N​αi​yi​Ki,1​+b=y1​
于是有：
b 1 n e w = y 1 − ∑ i = 3 N α i y i K i 1 − α 1 n e w y 1 K 11 − α 2 n e w y 2 K 21 b_1^{new}=y_1-\sum\limits_{i=3}^N{\alpha_iy_iK_{i1}-\alpha_1^{new}y_1K_{11}-\alpha_2^{new}y_2K_{21}} b1new​=y1​−i=3∑N​αi​yi​Ki1​−α1new​y1​K11​−α2new​y2​K21​
由 E 1 E_1 E1​的定义式得到：
E 1 = ∑ i = 3 N α i y i K i 1 + α 1 o l d y 1 K 11 + α 2 o l d y 2 K 21 + b o l d − y 1 E_1=\sum\limits_{i=3}^N\alpha_iy_iK_{i1}+\alpha_1^{old}y_1K_{11}+\alpha_2^{old}y_2K_{21}+b^{old}-y_1 E1​=i=3∑N​αi​yi​Ki1​+α1old​y1​K11​+α2old​y2​K21​+bold−y1​
将上式进行移项得到：
y 1 − ∑ i = 3 N α i y i K i 1 = − E 1 + α 1 o l d y 1 K 11 + α 2 o l d y 2 K 21 + b o l d y_1-\sum\limits_{i=3}^N{\alpha_iy_iK_{i1}}=-E_1+\alpha_1^{old}y_1K_{11}+\alpha_2^{old}y_2K_{21}+b^{old} y1​−i=3∑N​αi​yi​Ki1​=−E1​+α1old​y1​K11​+α2old​y2​K21​+bold
带入 b 1 n e w b^{new}_1 b1new​的求解式中，可以得到 b 1 n e w b^{new}_1 b1new​的更新式：
b 1 n e w = − E 1 − y 1 K 11 ( α 1 n e w − α 1 o l d ) − y 2 K 21 ( α 2 n e w − α 2 o l d ) + b o l d b^{new}_1=-E_1-y_1K_{11}(\alpha_1^{new}-\alpha_1^{old})-y_2K_{21}(\alpha_2^{new}-\alpha_2^{old})+b^{old} b1new​=−E1​−y1​K11​(α1new​−α1old​)−y2​K21​(α2new​−α2old​)+bold
如果 0 < α 2 n e w < C 0<\alpha^{new}_20<α2new​<C，那么：
b 2 n e w = − E 2 − y 1 K 12 ( α 1 n e w − α 1 o l d ) − y 2 K 22 ( α 2 n e w − α 2 o l d ) + b o l d b^{new}_2=-E_2-y_1K_{12}(\alpha_1^{new}-\alpha_1^{old})-y_2K_{22}(\alpha_2^{new}-\alpha_2^{old})+b^{old} b2new​=−E2​−y1​K12​(α1new​−α1old​)−y2​K22​(α2new​−α2old​)+bold
如果 α 1 n e w , α 2 n e w \alpha_1^{new},\alpha_2^{new} α1new​,α2new​同时满足条件 0 < α i n e w < C 0<\alpha_i^{new}0<αinew​<C的条件，那么 b 1 n e w = b 2 n e w b_1^{new}=b_2^{new} b1new​=b2new​，如果 α 1 n e w \alpha_1^{new} α1new​和 α 2 n e w \alpha_2^{new} α2new​是0或者$C$那么 b 1 n e w b_1^{new} b1new​和 b 2 n e w b_2^{new} b2new​以及它们之间的数都是符合KKT条件的，所以这时候选择中点作为 b n e w b^{new} bnew。
除了$b$的值需要更新外，另一个需要更新的值就是$E$，我们每求解一次子问题时，都需要对$E$的值进行更新。更新方式如下：
E i n e w = ∑ S y j α j K ( x i , x j ) + b n e w − y i E^{new}_i=\sum\limits_{S}y_j\alpha_jK(x_i,x_j)+b^{new}-y_i Einew​=S∑​yj​αj​K(xi​,xj​)+bnew−yi​
其中$S$是所有支持向量的集合。

SMO算法的实现

  SMO的理论部分终于结束，我相信绝大多数人都不一定会喜欢这么多的数学公式的，但是没办法，想要实现算法必须得了解每一道公式得含义，既然已经学习完理论部分，我们就可以上手编程，使用的软件为MATLAB。整理编程主要就三大版块：选择变量，更新变量，更新参数
代码如下：

clc;
clear;
tic
%% 数据提取
data = [27    53    -149    37    -156    39    -128    60    -168    75     157    69     164    62     177    68     170    54     156    63     125    41    -166    34     155    79     177    31    -146    66     130    23    -121    45    -168    42    -143    43    -156    59     179    68     160    34    -149    32    -180    79     177    46     126    66     129    29    -177    34     120    71    -149    25    -158    65     133    57    -131    79     120    78     177    37    -173    34    -160    26    -177    66     171    75     135    36    -149    61     126    37    -142    73     136    50    -166    73     171    43     133    62     143    41    -142    29    -158    20    -1];
Alpha = [];
% save('data.mat') % 保存数据，不需要运行 
X = data(:,1:2); % 提取特征向量
y = data(:,3); % 提取标签% data_1 =[];data_2 = [];
% for i  = 1:length(data)
%     if data(i,3)==1
%         data_1 = [data_1;data(i,1:2)];
%     else
%         data_2 = [data_2;data(i,1:2)];
%     end
% end
% scatter(data_1(:,1),data_1(:,2),'*r')
% hold on
% scatter(data_2(:,1),data_2(:,2),'*b')
%% SMO序列最小优化算法
b = 0; % 初始化b
[m,~] = size(X); % 获得特征向量的大小
alphas = zeros(m,1); % 构造初始解
iter = 0; % 初始化迭代次数
err = 0.01; % 这是误差允许范围
Iter_max = 300; % 最大迭代次数
E = zeros(m,1); % 初始化存储矩阵
i_isok =zeros(m,1); % 记录无法选择的变量
j_isok = zeros(m,1); % 记录无法选择的变量
MSE = zeros(Iter_max,1);
C = 0.01; % 惩罚参数
while (iter <= Iter_max)disp(['当前迭代的次数为：',num2str(iter)])% 选择第一个变量violatedMax = 0; % 定义中间变量I = 0; % 初始化第一个变量的下标for i =1:m % 遍历循环，选择第一个变量G_X = (X(i,:)*X')*(alphas .* y) + b; % 先计算G_X的值E(i) = G_X - y(i); % 计算E的值% 判断KKT条件% 如果违反第一个条件if (alphas(i) == 0 && y(i) * G_X < 1 - err)err_i = abs(1 - err - y(i) * G_X); % 计算误差的绝对值大小% 如果误差比最大误差还大并且在i的记录中没有当前i时i，选择此i为第一个优化变量if (err_i > violatedMax && i_isok(i) == 0)violatedMax = err_i; % 更新最大误差I = i; % 更新当前i为第一个变量endend% 如果违反第二个条件if (alphas(i) > 0 && alphas(i) < C && (y(i) * G_X < 1 - err) || (y(i) * G_X > 1 + err) )err_i = max(abs(1 - err - y(i) * G_X),abs(1 + err - y(i) * G_X)); % 计算此时的误差绝对值大小% 如果误差比最大误差还大并且在i的记录中没有当前i时i，选择此i为第一个优化变量if (err_i > violatedMax && i_isok(i) == 0)violatedMax = err_i; % 更新最大误差I = i; % 更新当前i为第一个变量endend% 如果违反第三个条件if (alphas(i) == C && y(i) * G_X > 1 + err)err_i = abs(1 + err - y(i) * G_X); % 计算误差的绝对值大小% 如果误差比最大误差还大并且在i的记录中没有当前i时i，选择此i为第一个优化变量if (err_i > violatedMax && i_isok(i) == 0)violatedMax = err_i; % 更新最大误差I = i; % 更新当前i为第一个变量endend% 第一个变量选择的结果是要么I = i,要么I = -1;end% 找到第一个变量后，下面找第二个变量J = 0; % 初始化第二个变量的标记if (I~=0) % 如果第一个变量已经被选% 第一种情况，详细参考《统计学习方法》147页if E(I) > 0E2_value = max(E); % 初始化E2的中间值为E的最大，后续逐渐找到最小值for j = 1:m% 如果当前的j满足条件，则选择当前的jif (E(j) < E2_value && j ~= I && j_isok(j)==0)E2_value = E(j); % 将更小的值赋给E2_valueJ = j; % 并将选择次j为第二个变量endendend% 第二种情况if E(I) < 0E2_value = min(E); % 初始化E2_value的值for j = 1:m% 如果当前的j满足条件，则选择当前的jif (E(j) > E2_value && j ~= I &&  j_isok(j)==0)E2_value = E(j); % 将更大的值赋给E2_valueJ = j; % 并将选择次j为第二个变量endendend% 如果第一个变量没有选，第二个变量也没有选，则说明循环结束if I==0 && J==0 break;end% 在此I的情况下未找到合适的J，需要重新寻找Iif I ~= 0 && J == 0j_isok = zeros(m,1); % 清空禁搜表i_isok(I) = 1; % 说明这个I的情况下没有合适的j值，需要重新挑选continueend% 下面计算变量的更新方法alphaI_old = alphas(I);alphaJ_old = alphas(J);% 计算H和L，具体参考《统计学习方法》的145页if (y(I) ~= y(J))L = max(0,alphaJ_old-alphaI_old);H = min(C,C + alphaJ_old-alphaI_old);elseL = max(0,alphaJ_old + alphaI_old - C);H = min(C,alphaJ_old + alphaI_old);endif L == Hj_isok(J) = 1; % 并将当前的第二个变量放进禁搜表中continueend% 如果当前的第二个变量已经选中，则进行后续的计算eta = X(I,:) * X(I,:)' + X(J,:) * X(J,:)' -  2 * X(I,:) * X(J,:)';alphasJ_unCut = alphaJ_old + y(J) * (E(I) - E(J)) / eta; % 计算未剪辑的解% 这里是通过将剪辑后的解与剪辑前的解做对比，如果变化小，则重新选择第二个变量if (abs(clipAlpha(alphasJ_unCut,H,L) - alphaJ_old) < 0.0001)j_isok(J) = 1;continue;end% 更新I，Jalphas(J) = clipAlpha(alphasJ_unCut,H,L);alphas(I) = alphaI_old + y(I) * y(J) * (alphaJ_old - alphas(J));% 更新b1b1 = b - E(I) - y(I) * X(I,:) * X(I,:)' * (alphas(I) - alphaI_old) - y(J) * X(J,:) * X(I,:)'* (alphas(J) - alphaJ_old);% 更新b2b2 = b - E(J) - y(I) * X(I,:) * X(J,:)' * (alphas(I) - alphaI_old) - y(J) * X(J,:) * X(J,:)' * (alphas(J) - alphaJ_old);% 更新总的bif (0 < alphas(I) && C > alphas(I))b = b1;elseif (0 < alphas(J) && C > alphas(J))b = b2;elseb = (b1 + b2) / 2;endend% 如果所有的I都是满足KKT条件的，就结束循环iter = iter + 1; % 迭代次数加一i_isok = zeros(m,1); % 清空禁选表j_isok = zeros(m,1); % 清空禁选表
%     draw(alphas,X,y)
%     pause(0.1)
%     hold off
MSE(iter) = Distance(alphas,X,y); % 计算迭代误差
Alpha = [Alpha,alphas];
end

修正函数：

function [alpha] = clipAlpha(alpha,H,L)
%CLIPALPHA
% 对alpha的值进行剪枝
if alpha > Halpha = H;
end
if L > alphaalpha = L;
end
end

绘图函数：

function  draw(alphas,X,y)
x = alphas;  % 赋值
% 利用求解得到x求解系数w
w = [0,0]; % 初始化系数，方便绘图
[a,~] = find(x ~= 0); % 找到可以求解b的值
temp = 0;
% 计算偏置系数
for i = 1:length(X)w = w + x(i) * y(i) * X(i,:);temp = temp + x(i) * y(i) * X(i,:) * X(a(1),:)';
end
% 计算偏置系数
b = y(a(1)) - temp;
% 数据可视化
k = - w(1)/w(2); % 构造截距式
b_ = -b/w(2); % 构造系数b
m = 10:2:90; % 生成一些点
n = k*m+b_; % 计算点的坐标
plot(m,n,'--') % 绘图
n_2 = k*m+b_+1/w(2);
hold on
plot(m,n_2,'--')
n_3 = k*m+b_-1/w(2);
plot(m,n_3,'--');
legend('类别1','类别二','分界线')
% 分离出两类样本点
X_1 = [];
X_2 = [];
for i =1:length(y)if y(i)==-1X_1 = [X_1;X(i,:)];elseX_2 = [X_2;X(i,:)];end
end
% 绘图
scatter(X_1(:,1),X_1(:,2),'filled','r')
hold on
scatter(X_2(:,1),X_2(:,2),'filled','b')
toc
end

  为了便于以后理解，上面大部分代码都已经加好了注释，并且就是按照文章实现的思路进行的。
其中部分环节没有和SMO算法叙述的完全相同，因为在第二个变量选择时，我是直接遍历所有的样本然后找到 ∣ E 1 − E 2 ∣ |E_1-E_2| ∣E1​−E2​∣最大的值作为第二个变量的，没有考虑目标函数下降的问题，实验表明这样也是可以的。

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！