学习支持向量机时参考了很多大神的博客，和经典著作，从公式推导到代码实现，亲历亲为。

遇到很多疑惑，也是各种百度+Google，虽然最后也差不多都解决了，但终归是因为数学基础不扎实。

在这里做一个总结，以便以后复习。

有些地方加入了我自己的理解，仅作参考。

当然，如有疑问，请各位不吝赐教。

线性不可分的情况和核函数的部分本文没有涉及。（主要是作者本人也不是很懂。。。

参考文献及完整代码在文末。

码字不易，不如点赞支持一下这个推公式推到头秃的博主？

---

1. SVM解决的问题

SVM解决的问题是经典的二元分类问题。给出一个分类标准使得样本集可以被最好地分类。在样本特征是二维的情况下，可以用下图表示

上图是我用Python实现的SVM分类器的最终结果图。

中间的实线是我们最终需要的分割线，在三维及以上的情况下叫做划分超平面。

画圈的样本是距离分割线最近的样本，叫做支持向量，这也是支持向量机名字的由来。两条黄色的虚线之间的距离叫做间隔，显然，这个间隔越大，也就代表两边的样本离分割线越远，我们得到的分割线就越鲁棒。

SVM的最终目的就是求出分割线的所有参数，在这个过程中还要确定样本中的支持向量是哪些。

2. 目标函数

我们的目的是找到使间隔最大的支持向量和分割平面，那么下面就要找到间隔的数学表达，也就是目标函数。

首先，样本集表示为：
D = { ( x 1 , y i ) , ( x 2 , y 2 ) . . . ( x m , y m ) } D=\{(\boldsymbol x_1,y_i), (\boldsymbol x_2, y_2)...(\boldsymbol x_m,y_m)\} D={(x1​,yi​),(x2​,y2​)...(xm​,ym​)}
其中，${\mathbf{x}}_i\in {\mathrm{R}}^n$, y i ∈ { + 1 , − 1 } , i = 1 , 2 , . . . , m y_i\in\{+1, -1\}, i=1,2,...,m yi​∈{+1,−1},i=1,2,...,m.
${\mathbf{x}}_i$ 是n维向量，为了便于可视化，我们以二维为例。
划分超平面表示为：
$$f(x)={\mathbf{w}}^T\mathbf{x}+b$$
我们将这个超平面标记为 $(\mathbf{w},b)$
样本点到超平面的距离表示为：
$$\gamma =\frac{|{\mathbf{w}}^T\mathbf{x}+b|}{||\mathbf{w}||}$$
假设超平面 $(\mathbf{w},b)$ 可以将样本点正确分类，那么每个样本集中的点到超平面的距离 ${\gamma }_i$ 应该都大于等于支持向量（假设表示为 ${\mathbf{x}}_0$）到超平面的距离 $\hat{\gamma }$:
$$\frac{|{\mathbf{w}}^T{\mathbf{x}}_i+b|}{||\mathbf{w}||}\ge \frac{|{\mathbf{w}}^T{\mathbf{x}}_0+b|}{||\mathbf{w}||},i=1,2,...m$$
那么支持向量与超平面的间隔就是不等式右边的部分。但是它过于复杂，所以这里对$(\mathbf{w},b)$ 进行缩放变换，使得 $|{\mathbf{w}}^T{\mathbf{x}}_0+b|=1$。这样处理对结果不会有影响，原因如下：

假设 ${\mathbf{w}}^T{\mathbf{x}}_0+b=c,c\in R$，因为支持向量不在分割超平面上，所以 $c\ne 0$。因此只要两边同除 $c$ 就可以了。

变换之后的不等式为：
$$\frac{|{\mathbf{w}}^T{\mathbf{x}}_i+b|}{||\mathbf{w}||}\ge \frac{1}{||\mathbf{w}||},i=1,2,...m$$

两个异类支持向量之间的间隔表示为：$2/||\mathbf{w}||$，就是目标函数。对上式两边同时乘以 $||\mathbf{w}||$ 可得，对样本集 $D$ 的所有点满足：
$$\{\begin{array}{l}{\mathbf{w}}^T{\mathbf{x}}_i+b\ge +1,y=+1\\ {\mathbf{w}}^T{\mathbf{x}}_i+b\le -1,y=-1\end{array}$$
合并一下就是：$y_{i}f(x_i)=y_i({\mathbf{w}}^T{\mathbf{x}}_i+b)\ge 1$

《机器学习》书中对这部分的描述如下：

目标函数的数学表达：
$$\underset{\mathbf{w},b}{\max }\frac{2}{||\mathbf{w}||}$$

为了方便计算，将上式取倒数，并对 $||\mathbf{w}||$ 取平方，得到最终需要优化的目标函数：
$$\begin{array}{rl}& \underset{\mathbf{w},b}{\min }\frac{1}{2}||\mathbf{w}|{|}^2\\ & \text{s.t. }{y}_i({\mathbf{w}}^T{\mathbf{x}}_i+b)\ge 1,i=1,2,...m\end{array}$$

3. 目标函数的优化

如果没有约束条件，这就是一个简单的多项式求极值的问题，只需令一阶导数等于零，二阶导数小于零求出$\mathbf{w}$的所有元素即可。
但是有了约束条件就要使用些特殊方法，比如拉格朗日乘子法。

为什么要用拉格朗日乘子法？
在特征维度低的线性问题中是可以不使用的，具体参考：https://www.zhihu.com/question/36694952

3.1 拉格朗日乘子法（Lagrange Multiplier）

这里简单介绍一下拉格朗日乘子法及KKT条件的原理和使用。
对有一个等式约束的凸函数优化问题：
$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})\\ & \text{s.t. }h(\mathbf{x})=0\end{array}$$
在三维的情况下，$h(\mathbf{x})=0$ 是一个曲面，这个曲面在坐标轴 x-y 平面上的投影是一条曲线，也就是，看起来和 x-y 平面 “垂直”。
$f(\mathbf{x})$是一个凸函数，其与$h(\mathbf{x})$相交。

上图中红色代表$f(\mathbf{x})$的等值线，蓝色线代表$h(\mathbf{x})$在x-y平面上的投影。

而两个曲面相交的部分是一条曲线，这条曲线在$f(\mathbf{x})$曲面上的最低点（就是图中黄色的点）就是$f(\mathbf{x})$取最小值的点。

这个点一定与$f(\mathbf{x})$的某一条等值线相切。且在切点处两函数的梯度共线（可以同向或反向）。即存在$\lambda \ne 0$，使得在最小值点处：
$$\nabla f({\mathbf{x}}^{\ast })+\lambda \nabla h({\mathbf{x}}^{\ast })=0$$
其中，${\mathbf{x}}^{\ast }$是最小值点。此时，令拉格朗日函数为：
$$L(\mathbf{x},\lambda )=f(\mathbf{x})+\lambda h(\mathbf{x})$$
对该函数求导得：
$$\nabla L(\mathbf{x},\lambda )=\nabla f(\mathbf{x})+\lambda \nabla h(x)$$
显然，当$\mathbf{x}={\mathbf{x}}^{\ast }$时，$\nabla L(\mathbf{x},\lambda )=0$。也就是说，$L(\mathbf{x},\lambda )$与$f(\mathbf{x})$同时取得极小值。

也就是说，带有等式约束的优化问题转化为了不带约束条件的优化问题。

令 $L(\mathbf{x},\lambda )$ 对 $\mathbf{x},\lambda$ 的导数都为0就可以得到最小值。

3.2 KKT条件

以上方法适用于在约束为等式的情况。当约束为不等式时，即:
$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})\\ & \text{s.t. }g(\mathbf{x})\le 0\end{array}$$

极小值点的位置存在以下两种情况（以一个不等式约束为例）：

• 极小值点在 $g(\mathbf{x})<0$ 区域内，此时可以忽略约束条件直接对$f(\mathbf{x})$求极小值即可，相当于把上节中的 $\lambda$ 置零；
• 极小值点在边界 $g(\mathbf{x})=0$ 上，此时与3.1节所述情况一样。唯一不同的是，在上图所示的极值点处，$f(\mathbf{x})$ 与 $g(\mathbf{x})$ 的梯度方向一定相反。（如果它们梯度方向相同，就可以继续向两者减小的方向优化）

将以上两种情况结合可得，存在$\mu \ge 0$ 使得：
$$\mu g(\mathbf{x})=0$$
由此引出$f(\mathbf{x})$ 取得极小值的充分必要条件，即KKT(Karush-Kuhn-Tucker)条件：
$$\{\begin{array}{l}g(\mathbf{x})\le 0\\ \mu \ge 0\\ \mu g(\mathbf{x})=0\end{array}$$
原函数 $f(\mathbf{x})$ 转化为与上节形式相同的拉格朗日函数。

将这几种情况结合起来，推广到具有多个不等式或等式约束的优化问题上：

对有有限个约束的最小化问题：
$$\begin{array}{rl}& \underset{\mathbf{x}}{\min }f(\mathbf{x})\\ \text{s.t. }& h_i(\mathbf{x})=0,i=1,,2,...m\\ & g_j(\mathbf{x})\le 0,j=1,,2,...n\end{array}$$
其拉格朗日函数为：
$$L(\mathbf{x},{\lambda }_i)=f(\mathbf{x})+\sum _{i=1}^m{\lambda }_i{h}_i(\mathbf{x})+\sum _{j=1}^n{\mu }_j{g}_j(\mathbf{x})$$

KKT条件为：
$$\{\begin{array}{l}{h}_i(\mathbf{x})=0\\ g_j(\mathbf{x})\le 0\\ {\lambda }_i,{\mu }_j\ge 0\\ {\mu }_j{g}_j(\mathbf{x})=0\end{array}$$

3.3 对偶问题

现在我们回到第2节末尾，再看这个优化问题：
$$\begin{array}{rl}& \underset{\mathbf{w},b}{\min }\frac{1}{2}||\mathbf{w}|{|}^2\\ & \text{s.t. }{y}_i({\mathbf{w}}^T{\mathbf{x}}_i+b)\ge 1,i=1,2,...m\end{array}$$
这就是个适用于不等式约束的优化。拉格朗日函数如下：
$$L(\mathbf{w},b,\alpha )=\frac{1}{2}||\mathbf{w}|{|}^2+\sum _{i=1}^m{\alpha }_i(1-y_i({\mathbf{w}}^T{\mathbf{x}}_i+b))，{\alpha }_i\ge 0$$
若要求$L(\mathbf{w},b,\alpha )$ 的最小值则首先把 $\mathbf{w},b$ 看作常量，对$L(\mathbf{w},b,\alpha )$求关于$\alpha$的最大值。因为 $1-y_i({\mathbf{w}}^T{\mathbf{x}}_i+b)$ 是小于等于0的，在${\alpha }_i\ge 0$的前提下，其最大值就是：$\frac{1}{2}||\mathbf{w}|{|}^2$。

因此，原问题可写为：
$$\underset{\mathbf{w},b}{\min }\underset{{\alpha }_i}{\max }L(\mathbf{w},b,\alpha )$$
直接对该问题求解比较困难，因此这里引入对偶问题。

使用对偶问题求解的原因不止是为了计算方便，还是为线性不可分时核函数的引入作铺垫（本文不涉及这部分）。

先上结论：
$$\underset{\mathbf{w},b}{\min }\underset{{\alpha }_i}{\max }L(\mathbf{w},b,\alpha )=\underset{{\alpha }_i}{\max }\underset{\mathbf{w},b}{\min }L(\mathbf{w},b,\alpha )$$

也就是求取最大最小值的顺序发生了变化。从先对$\alpha$求最大值变为先对 $\mathbf{w},b$ 求最小值。

等式两边的问题互为对偶问题。

当然两个对偶问题同时取得最优解是有条件的，具体证明可以参考：https://www.cnblogs.com/90zeng/p/Lagrange_duality.html

转化为对偶问题后，令 $\mathbf{w},b$ 的导数为0：
$$\{\begin{array}{l}\frac{\partial L}{\partial \mathbf{w}}=\mathbf{w}-{\sum }_{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i=0\\ \frac{\partial L}{\partial b}=-{\sum }_{i=1}^m{\alpha }_i{y}_i=0\end{array}$$
将上式代入 $L(\mathbf{w},b,\alpha )$ 并化简可得：
$$L(\mathbf{w},b,\alpha )=\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$$
化简过程如下：（公式太多还是手写吧emmm…）

化简后就只用求 $\alpha$ 即可。这样问题的复杂度就从与样本的维度有关变为只与样本数量有关了。

周志华的《机器学习》一书对这个最优解的几何意义有一段描述可以参考：

3.4 SMO(Sequential Minimal Optimazation)方法

上节最后得到的关于 $\alpha$ 的式子是一个二次规划问题，使用通用算法开销比较大。SMO方法是1998年提出的，用于求取SVM中的拉格朗日乘子，比通用算法更加高效。

其主要思想为，选取两个变量${\alpha }_i,{\alpha }_j$，固定其他参数，对这两个参数进行优化，然后重复这个过程。

Note：选取${\alpha }_i,{\alpha }_j$的标准：

不过博主有点懒，所以直接随机选取了。。。。需要深入了解可以参考《机器学习》（文末有下载方法~）

每次选取两个的原因是这些参数之间有约束 ${\sum }_{i=1}^m{\alpha }_i{y}_i=0$，只改变一个约束会被打破。

为了表示方便，记选取的两个参数为${\alpha }_1,{\alpha }_2$，固定除${\alpha }_1,{\alpha }_2$以外的参数：
$$\begin{gathered} {\alpha }_1{y}_1+{\alpha }_2{y}_2=-\sum _{i\ne 1,2}{\alpha }_i{y}_i=\xi \\ {\alpha }_1=y_1(\xi -{\alpha }_2{y}_2) \end{gathered}$$

带入$L(\mathbf{w},b,\alpha )$得：

为简便表示，下式中，$K_{ij}={\mathbf{x}}_i^T{\mathbf{x}}_j$

$$\begin{array}{rl}& L(\mathbf{w},b,\alpha )\\ =& \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{K}_{ij}\\ =& \sum _{i\ge 3}^m{\alpha }_i+y_1\xi -{\alpha }_2{y}_1{y}_2+{\alpha }_2-\frac{1}{2}\sum _{i\ge 3}^m\sum _{j\ge 3}^m{\alpha }_i{\alpha }_j{y}_i{y}_j{K}_{ij}-\frac{1}{2}(y_1\xi -{\alpha }_2{y}_1{y}_2{)}^2{K}_{11}-y_2{\alpha }_2(\xi -{\alpha }_2{y}_2)K_{12}-\frac{1}{2}{\alpha }_2^2{K}_{22}\\ & -\sum _{i\ge 3}^m(\xi -{\alpha }_2{y}_2){\alpha }_i{y}_i{K}_{1i}-\sum _{i\ge 3}^m{\alpha }_2{y}_2{\alpha }_i{y}_i{K}_{2i}\\ =& -{\alpha }_2{y}_1{y}_2+{\alpha }_2-\frac{1}{2}(y_1\xi -{\alpha }_2{y}_1{y}_2{)}^2{K}_{11}-y_2{\alpha }_2(-{\alpha }_2{y}_2)K_{12}-\frac{1}{2}{\alpha }_2^2{K}_{22}\\ & -\sum _{i\ge 3}^m(\xi -{\alpha }_2{y}_2){\alpha }_i{y}_i{K}_{1i}-\sum _{i\ge 3}^m{\alpha }_2{y}_2{\alpha }_i{y}_i{K}_{2i}+constant\\ =& -{\alpha }_2{y}_1{y}_2+{\alpha }_2-\frac{1}{2}{\alpha }_2^2{K}_{11}+{\alpha }_2{y}_2\xi K_{11}-{\alpha }_2{y}_2\xi K_{12}+{\alpha }_2^2{K}_{12}-\frac{1}{2}{\alpha }_2^2{K}_{22}\\ & +\sum _{i\ge 3}^m{\alpha }_2{y}_2{\alpha }_i{y}_i{K}_{1i}-\sum _{i\ge 3}^m{\alpha }_2{y}_2{\alpha }_i{y}_i{K}_{2i}+constant\end{array}$$
因为我们的最终目的是让 $L(\mathbf{w},b,\alpha )$ 对 ${\alpha }_2$ 求导，而其他 $\alpha$ 已被固定，视为常数。且 $x,y$ 也都是已知常数，所以将上式中不含 ${\alpha }_2$ 的常数项合并为 $constant$。

下一步，令其 ${\alpha }_2$ 的导数为0：

$$\frac{\partial L(\mathbf{w},b,\alpha )}{\partial {\alpha }_2}=-y_1{y}_2+1-{\alpha }_2{K}_{11}+y_2\xi K_{11}-y_2\xi K_{12}+2{\alpha }_2{K}_{12}-{\alpha }_2{K}_{22}+\sum _{i\ge 3}^m{y}_2{\alpha }_i{y}_i{K}_{1i}-\sum _{i\ge 3}^m{y}_2{\alpha }_i{y}_i{K}_{2i}=0$$

由上节可知：
$$\begin{gathered} \mathbf{w}=\sum _{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i \\ f({\mathbf{x}}_1)={\mathbf{w}}^T{x}_1+b=\sum _{i=1}^m{\alpha }_i{y}_i{K}_{1i}+b \\ f({\mathbf{x}}_2)={\mathbf{w}}^T{x}_2+b=\sum _{i=1}^m{\alpha }_i{y}_i{K}_{2i}+b \end{gathered}$$

注意到这个形式和上面的最后两项很相似，所以可以带进去计算。

但是需要注意的是这里的$f(\mathbf{x})$是常量（否则带入一个变量并不利于求导结果的计算），所以这里的$f(\mathbf{x})$中包含的 ${\alpha }_i(i=1,2,...,m)$ 是初始化时赋予的值，将其记作 ${\alpha }_2^{old}$；把其他的 ${\alpha }_2$ 看作待更新的变量，记作 ${\alpha }_2^{new}$。

$$\begin{array}{rl}\frac{\partial L(\mathbf{w},b,\alpha )}{\partial {\alpha }_2}=& -y_1{y}_2+1-{\alpha }_2^{new}{K}_{11}+y_2\xi K_{11}-y_2\xi K_{12}+2{\alpha }_2^{new}{K}_{12}-{\alpha }_2^{new}{K}_{22}\\ & +y_2[f(x_1)-(\xi -{\alpha }_2^{old}{y}_2)K_{11}-{\alpha }_2^{old}{y}_2{K}_{12}-b]-y_2[f(x_2)-(\xi -{\alpha }_2^{old}{y}_2)K_{12}-{\alpha }_2^{old}{y}_2{K}_{22}-b]\\ =& (-K_{11}-K_{22}+2K_{12}){\alpha }_2^{new}+(K_{11}+K_{22}-2K_{12}){\alpha }_2^{old}+1+y_1{y}_2+y_2[f(x_1)-f(x_2)]\\ =& 0\end{array}$$

令 $\eta =K_{11}+K_{22}-2K_{12}$，得：
$$\begin{gathered} \frac{\partial L(\mathbf{w},b,\alpha )}{\partial {\alpha }_2}=-\eta {\alpha }_2^{new}+\eta {\alpha }_2^{old}+1+y_1{y}_2+y_2[f({\mathbf{x}}_1)-f({\mathbf{x}}_2)]=0 \\ {\alpha }_2^{new}={\alpha }_2^{old}+\frac{1+y_1{y}_2+y_2[f({\mathbf{x}}_1)-f({\mathbf{x}}_2)]}{\eta } \end{gathered}$$
同时，因为需要求最大值，所以还要求对 ${\alpha }_2$ 的二阶导小于零：

$$\begin{gathered} \frac{{\partial }^{2}L(\mathbf{w},b,\alpha )}{{\partial }^2{\alpha }_2}=-K_{11}-K_{22}+2K_{12}\le 0 \\ \eta =K_{11}+K_{22}-2K_{12}\ge 0 \end{gathered}$$

到这里为止，${\alpha }_2$ 的更新已经完成，${\alpha }_1$ 可以通过两者之间的约束求得。
$${\alpha }_1^{new}{y}_1+{\alpha }_2^{new}{y}_2={\alpha }_1^{old}{y}_1+{\alpha }_2^{old}{y}_2=\xi$$
最后一步，就是求解参数 $b$。

本文只采用了前一种方法，感兴趣的童鞋可以用（6.18）的方法试试，欢迎在评论区讨论~

具体操作如下：

假设更新过后的 ${\alpha }_2^{new}>0$，即该 ${\alpha }_2^{new}$ 对应的样本为支持向量，则满足 $y_2{f}^{new}({\mathbf{x}}_2)=1$，结合更新之前的 $f({\mathbf{x}}_2)$：

$$\{\begin{array}{l}f({\mathbf{x}}_2)={\sum }_{i\ge 3}^m{\alpha }_i{y}_i{K}_{2i}+{\alpha }_1^{old}{y}_1{K}_{11}+{\alpha }_2^{old}{y}_2{K}_{12}+b^{old}\\ y_2({\sum }_{i\ge 3}^m{\alpha }_i{y}_i{K}_{2i}+{\alpha }_1^{new}{y}_1{K}_{11}+{\alpha }_2^{new}{y}_2{K}_{12}+b^{new})=1\end{array}$$
解得：
$$b^{new}=y_2-f({\mathbf{x}}_2)-({\alpha }_1^{new}-{\alpha }_1^{old})y_1{K}_{11}+({\alpha }_2^{new}-{\alpha }_2^{old})y_2{K}_{12}$$

4 软间隔

由于不是所有的数据集都可以完美的用一条线分割开，在两个类别中间可能会有个别样本点超过分割平面，使得数据集不能被完美分开。即便可以被分开，也有可能是过拟合造成的。

这时我们的算法需要对这些异常点进行容错，使其不影响数据集正常划分。

异常点的表示就是其不满足约束条件：

$$y_i({\mathbf{w}}^T{\mathbf{x}}_i+b)\ge 1$$

不过我们希望不满足约束的样本越少越好，因此将原问题改为：

$$\underset{\mathbf{w},b}{\min }L(\mathbf{w},b)=\underset{\mathbf{w},b}{\min }\frac{1}{2}||\mathbf{w}|{|}^2+C\sum _{i=1}^m{l}_{0/1}(y_i({\mathbf{w}}^T{\mathbf{x}}_i+b)-1)$$

$C>0$ 是一个常数。当其取较大的值时，就迫使更多的样本趋于满足约束条件，取无穷大时就和原问题等价。

$l_{0/1}$叫做“0/1损失函数”，数学表示为：

$$l_{0/1}(z)=\{\begin{array}{ll}1,& \text{if z<0 }\\ 0,& \text{else}\end{array}$$

但是它数学性质不好，不连续，所以我们用合页损失函数（Hinge loss）代替：

$$l_{hinge}(z)=\max (0,1-z)$$

如下图（图片来自《机器学习》）：

此时，优化目标变为：

$$\underset{\mathbf{w},b}{\min }\frac{1}{2}||\mathbf{w}|{|}^2+C\sum _{i=1}^m\max (0,1-y_i({\mathbf{w}}^T{\mathbf{x}}_i+b))$$

这里有一个问题，按照合页损失函数的定义，把前面式子中的求和项带入，得到的结果应该是： $$1-z=1-(y_i({\mathbf{w}}^T{\mathbf{x}}_i+b)-1)=2-y_i({\mathbf{w}}^T{\mathbf{x}}_i+b)$$
个人理解是，根据我们要解决的问题对hinge损失函数做了改动：
$$l_{hinge}(z)=\max (0,-z)$$

优化目标里含有max()，显然不便于计算，所以这里引入“松弛变量(slack variables)” ${\xi }_i\ge 0$，表示每个样本不满足约束的程度。若样本 $x_i$ 满足约束，则对应的${\xi }_i=0$；若不满足，则 $1-y_i({\mathbf{w}}^T{\mathbf{x}}_i+b)\le {\xi }_i$。优化目标可以重新写为：
$$\begin{array}{rl}& \underset{\mathbf{w},b}{\min }\frac{1}{2}||\mathbf{w}|{|}^2+C\sum _{i=1}^m{\xi }_i\\ \text{s.t. }& 1-y_i({\mathbf{w}}^T{\mathbf{x}}_i+b)\le {\xi }_i\\ & {\xi }_i\ge 0\end{array}$$

和前面一样，重新使用拉格朗日乘子法。拉格朗日函数：

$$\begin{gathered} L(\mathbf{w},b,\alpha ,\mu )=\frac{1}{2}||\mathbf{w}|{|}^2+C\sum _{i=1}^m{\xi }_i+\sum _{i=1}^m{\alpha }_i(1-y_i({\mathbf{w}}^T{\mathbf{x}}_i+b)-{\xi }_i)-\sum _{i=1}^m{\mu }_i{\xi }_i， \\ {\alpha }_i,{\mu }_i\ge 0 \end{gathered}$$

求解过程也和前面一样，贴个图：

和前面得到的结果对比一下，唯一不同的地方在于 ${\alpha }_i$ 的取值范围有了上界。这个上界来源于图片上的式（6.39），意义就是表示了算法对于样本点的容错程度。

5 代码

下面，我们将把以上的算法用Python程序实现。
首先我们对上面的算法步骤作一个总结，帮助我们理清思路：

整个SMO算法分为以下几个阶段：

• 选择 ${\alpha }_1,{\alpha }_2$
• 判断 $\eta =K_{11}+K_{22}-2K_{12}\ge 0$
• 计算 ${\alpha }_2^{new}={\alpha }_2^{old}+\frac{1+y_1{y}_2+y_2[f(x_1)-f(x_2)]}{\eta }$
• 对 ${\alpha }_1,{\alpha }_2$ 做上下界限定，范围为 $[0,C]$
• 判断 ${\alpha }_2$ 的改变量是否足够大，太小视为不变
• 计算 ${\alpha }_1^{new}={\alpha }_1^{old}+{\alpha }_2^{old}{y}_1{y}_2-{\alpha }_2^{new}{y}_1{y}_2$
• 计算 $b$
• 重复以上步骤

5.1 数据集生成及一些辅助函数

写了一个简单的函数用来随机生成数据集。

指定一条曲线 $f(x)$，本例中是 $f(x)=x$，随机均匀生成100对数据，取值为[0, 10]，标为两个类别并画图表示。（为了突出两数据集的间隔，这里舍弃与分割线垂直距离小于1的点）

import matplotlib.pyplot as plt
import numpy as np

def generate_dataset():
    # np.random.seed(3)
    class1 = []; class2 = []
    label1 = []; label2 =[]
    # define decision line
    def f(x):
        return 1*x
    for _ in range(100):
        x = np.random.rand() * 10
        y = np.random.rand() * 10
        if y-f(x) > 1:
            class1.append([x, y])
            label1.append(1)
        elif y-f(x) < -1:
            class2.append([x, y])
            label2.append(-1)
    x1, y1 = zip(*[data for data in class1])
    x2, y2 = zip(*[data for data in class2])
    plt.plot(x1, y1, 'ro')
    plt.plot(x2, y2, 'bo')
    plt.axis([0,10,0,10])
    plt.show()
    return class1+class2, label1+label2

随机选择第2个参数的下标：

def select_j(i, m):
    j = i
    while(j == i):
        j = np.random.randint(0, m)
    return j

限制 $\alpha$ 的上下界：

def clip_alpha(aj, H: float, L:float):
    return min(max(L, aj), H)

已知 $\alpha$ 求 $\mathbf{w}$：

def get_W(alphas, dataset, label):
    a, x, y = map(np.array, [alphas, dataset, label])
    W = np.dot(a * y, x)
    # transform W form np.array to Python List
    return W.tolist()

5.2 主函数

简化版的SMO算法：

def smo_simple(dataset, labels, C, max_iter):
    data = np.array(dataset, dtype=np.float)
    label = np.array(labels, dtype=np.float)
    b = 0
    m, n = np.shape(data)
    alphas = np.zeros(m, dtype=np.float)
    iter = 0
    while iter < max_iter:
        alpha_pair_changed = 0
        for i in range(m):
            x_i, y_i = data[i], label[i]
            fx_i = np.dot(alphas*label, np.dot(data, x_i)) + b
            e_i = fx_i - y_i
            j = select_j(iter, m)
            x_j, y_j = data[j], label[j]
            fx_j = np.dot(alphas*label, np.dot(data, x_j)) + b
            e_j = fx_j - y_j

            # calculate a_j_new
            eta = np.dot(x_i, x_i) + np.dot(x_j, x_j) - 2*np.dot(x_i, x_j)
            if eta <= 0:
                print("eta <= 0, continue")
                continue
            a_i, a_j = alphas[i], alphas[j]
            a_j_new = a_j + y_j * (e_i - e_j) / eta

            # limit a_j_new to [0, C]
            if y_i == y_j:
                L = max(0., a_i + a_j - C)
                H = min(a_i + a_j, C)
            else:
                L = max(0., a_j - a_i)
                H = min(C - a_i + a_j, C)
            if L < H:
                a_j_new = clip_alpha(a_j_new, H, L)
            else:
                print("L >= H. (L, H) =", (L, H))
                continue

            # judge if a_j moves an enough distance
            if abs(a_j_new - a_j) < 0.00001:
                print("a_j has not moved enough, a_j_new - a_j = %f" % (a_j_new - a_j))
                continue

            # calculate a_i_new and update a_i, a_j
            a_i_new = (a_j - a_j_new)*y_i*y_j + a_i
            alphas[i], alphas[j] = a_i_new, a_j_new
            alpha_pair_changed += 1

            #calculate b
            b_i = -e_i + (a_i - a_i_new)*y_i*np.dot(x_i, x_i) + (a_j - a_j_new)*y_j*np.dot(x_i, x_j) + b
            b_j = -e_j + (a_i - a_i_new)*y_i*np.dot(x_i, x_j) + (a_j - a_j_new)*y_j*np.dot(x_j, x_j) + b
            # b = b_i if b_i == b_j else (b_i + b_j)/2
            if 0 < a_i_new < C:
                b = b_i
            elif 0 < a_j_new < C:
                b = b_j
            else:
                b = (b_i + b_j)/2

            print("(a_i, a_j) moved from (%f, %f) to (%f, %f)." % (a_i, a_j, a_i_new, a_j_new))


        if alpha_pair_changed == 0:
            print("Iteration %d of max_iter %d" % (iter+1, max_iter))
            iter += 1
        else:
            iter = 0

    return alphas, b

main函数：

dataset, label = generate_dataset()
# print(label)
alphas, b = smo_simple(dataset, label, C=6, max_iter=40)
print(alphas)
W = get_W(alphas, dataset, label)
print(W, b)

*5.3 使用matplotlib绘图

如果想要可视化计算结果，可以在main函数后添加以下代码：

class1, class2 = [], []
for data in zip(dataset, label):
    if data[1] == 1.0:
        class1.append(data[0])
    elif data[1] == -1.0:
        class2.append(data[0])
x11, x12 = zip(*class1)
x21, x22 = zip(*class2)
plt.plot(x11, x12, 'ro')
plt.plot(x21, x22, 'bo')

x = np.linspace(0, 10, 50)
y = -(W[0]*x + b)/W[1]
plt.plot(x, y)
for i in range(len(dataset)):
    if alphas[i] > 1e-3:
        xi_1, xi_2 = dataset[i][0], dataset[i][1]
        plt.scatter(xi_1, xi_2, s=150, c='none', linewidths=1.5, edgecolors='#1f77b4')
        x = np.linspace(0, 10, 50)
        y = -W[0]/W[1] * x + (xi_2 + W[0]/W[1] * xi_1)
        plt.plot(x, y, 'y--')
plt.show()

6 github 地址

https://github.com/chenyr0021/machine_learning_exercises

7 参考文献

1. 《机器学习》，周志华
2. 《机器学习实战》，Peter Harrington
3. 简易解说拉格朗日对偶（Lagrange duality）
4. 支持向量机系列（5）——SMO算法解对偶问题
5. SVM为什么能够求解对偶问题，求解对偶问题为什么和原问题一样？为什么要求解对偶问题？svm的公式是什么？如果线性不可分怎么办？
6. 为什么支持向量机要用拉格朗日对偶算法来解最大化间隔问题？
7. 机器学习算法实践-SVM中的SMO算法
8. 支持向量机通俗导论（理解SVM的三层境界）

都看到这里了，不点赞支持一下头秃博主吗

---

知乎：@陈小白233
公众号：一本正经的搬砖日常

扫码关注公众号即可获赠《机器学习》、《机器学习实战》电子书~