岁虚山

行有不得，反求诸己

0%

SVM 支持向量机

发表于 更新于 分类于

概述

支持向量机（Support Vector Machines，SVM）是一种分类算法。本章内容会对 SVM 的基本概念及原理进行阐述，对其数学原理进行推导。最后将通过 Python 实现一个简单的 SVM 分类器，来帮助理解 SVM 的具体实现过程和原理。

由于对 SVM 的数学推导过程十分复杂，涉及到大量的高等数学和线性代数方面的公式，因此本文尽量使用简单的描述来推导计算过程。

SVM 支持向量机

SVM 实质上是一个分类器，其作用就是对样本进行分类，接下来看一个简单的例子：

简单划分

在平面上有一些橙色和绿色的点，是否能画一条直线将这两类不同的点进行划分，使得不同颜色的点刚好处于直线的两边？

像上图一样，若存在一条直线能够将样本划分开来，我们称之为线性可分。对于线性可分的问题，可以用 SVM 取得很好的分类效果。

若不存在一条直线能够将样本划分完全开来，我们称之为线性不可分。对于线性不可分的问题，我们将在后面引入核函数来将问题转换为线性可分的问题。并且引入软间隔的概念，允许部分样本分在错误的一边，从而达到分类的目的。

本文的原理推导，都是基于样本在特征空间中线性可分的情况来讨论的。直到在软间隔和核函数部分才讨论线性不可分的情况。

SVM 适用于二分类问题，当然也可以经过改造，实现多分类器，但这不在本文的讨论范围内。本文主要讨论二分类问题下 SVM 的推导。

这是一个很简单的问题，我们可以很简单的在图上画出好几条不同的直线来将这些点进行划分，如下：

简单划分

对于上面的 A、B、C、D 四种画法，都能将不同颜色的点划分开来，那么哪一种画法更好一点呢？

直觉上来讲，我们在画这条直线的时候，想让画出的直线尽量处于两类样本的“中间”，这样直线两边都能留有充足的“空间”。这样，当新的样本添加进来的时候，其分类错误的概率比较小。

比如我们添加一些新的样本进去：

02 新加样本

上图中，用三角形代表新加的样本，颜色用来区分各自的类别，可以看到，A、B、C 的画法，会将新的样本进行错误的分类，而 D 中，样本仍然能够正确分类。

当然，上面的举例只是一个简单的二维平面上分类中的例子，我们可以简单的通过视觉和直觉来判断所谓的“中间”的位置在哪。但是当样本数量较大，且其特征空间是四维，或者更高维度时，我们没有办法去画出或者想象样本在高维空间中的分布情况，因此无法通过直觉来画出这条分割线。

同时，计算机并没有人类的视觉能力，我们无法编写一个程序，让它具备这种通过直觉划分样本的能力。

因此，我们需要对“中间”这个概念进行更加严格和清晰的定义，从而能够处理高维特征空间中的分类问题。

那么“中间”这个概念究竟指的是什么呢？我们再次给出一个不算严谨的定义：在特征空间中，存在能将两类样本划分开来的平面，能使得位于分割平面最近的不同分类的样本点距离最大且相等，我们说这个平面处于中间位置。

我们拿上面中的 C 和 D 举例：

03 中间示意

我们把样本到分割线的距离称之为间隔。在 A 中，分割线到两边最近的样本的间隔不一样；图 B 中，虽然分割线到两种分类最近的样本间隔相等，但分割线未能将样本划分成两类；C 中距离分割线最近的样本 b 和 c 的间隔不相等；只有 D 图中，分割线将样本划分为不同的两类，且距离分割线最近的样本 a、b、c 与分割线的间隔是相等的。我们说 D 图中的分割线具处于中间。

从上面的例子来看，“中间” 的位置知识由离分割线最近的几个样本说决定的，也就是说，我们的分割平面所在的位置不是由所有的样本所决定的，而是由部分样本所决定。我们称这些决定分割空间位置的样本向量为：支持向量。

通过支持向量，我们就能很简单得求出分割平面，而无需将所有的样本向量考虑在内。因此，支持向量机并不需要大量的样本进行训练，就能得到很好的分类效果。

原理推导

问题描述

假如我们有 个样本，每个样本的特征空间都是 N 维（每个样本都有 N 个特征值，每个特征代表一个维度）的，这些样本只属于两类，用数学语言来描述就是：

存在训练样本集

其中每一个样本 都是一个 维的向量，即：

我们希望找到一个超平面 能够划分样本，如下图所示：

04 超平面方程

其中， 也是一个 维的向量，即：

推导

我们知道，二维平面上的点 到直线 的距离公式为：

三维空间中，一点 到平面 的距离公式为：

推导到 维空间上，则 维空间中，点 到平面 的距离为：

其中 称之为向量 的二范数（Euclid 范数，用于计算向量长度），即向量元素的平方和再开方。

那么对于样本中的支持向量，则有：

即样本中的支持向量到超平面的距离都是 （距离 肯定是一个正数，上式中出现 的原因是因为未对等号左边的式子取绝对值，这是为了标识 的支持向量处于平面的上方， 的支持向量处于平面的下方）。

换个思路来理解就是将平面 沿着法线的方向上下平移距离 ，得到两个与 平行的两个上下边界平面： ：

此时，标签为 的支持向量位于上边界平面 上，标签 的支持向量位于下边界平面 上，于是对于这些支持向量则有：

支持向量刚好处于上下两个边界平面上，那么 的距离 。

06 上下边界距离示意图

上图中， a、b、c 三个样本就是支持向量，分别位于上下边界上。

那么对于样本中的非支持向量，它们要么处于上边界平面的上方，要么处于边界平面的下方，那么这些非支持向量的样本到平面的距离：

将上面两式合并就得到对于样本集合 中的样本：

将上式进行调整，等号两边同乘以 ，则：

因为距离 是一个正数，因此将等号两边同除以 ：

因为 是一个标量，对于向量而言，向量空间除以一个标量仅仅意味着对向量进行缩放，并不改变向量的方向和相对于其他向量的位置（因为其他向量也被缩放了）。

举个简单的例子，假如我们有一个平面坐标系 ，其中的两个向量 ，如下图：

07 空间缩放示意图

现在我们将平面坐标系进行 0.5 倍的缩放，得到新的平面坐标系 ，则相对的，在新的坐标系中，原向量 变成了 ，这两个向量的方向和相对位置并未改变，只是距离变成了原本的 。

那么，现在我们对 维空间缩放到原本的 ，则在新的空间中：

所有的样本向量 都变成了新的向量 ，平面 也变为了：

那么样本中的向量到平面的距离：

将 提取出来：

即：

我们对特征空间进行缩放，于是得到了：

因为我们只关心样本向量中的支持向量，换句话说，平面 的确定只和支持向量相关。我们的目标是让样本中的支持向量的距离最大，而支持向量恰好位于特征空间中的上下边界平面上，回过头来，看一下最开始上下边界平面的定义：

因为我们的特征空间进行了缩放，因此在缩放后的特征空间中，上下边界平面为：

08 空间缩放后上下边界距离示意图

那么位于上下边界平面上的支持向量则：

那么这些向量到平面 的距离（为了不产生混淆，在缩放后的特征空间中，我们用 来表示距离，并且在以后不再对缩放后的特征空间进行特殊的说明，在后续的推导中所提到的特征空间都是缩放后的特征空间）：

又因为支持向量都处于上下两个边界平面上，因此联立上面两式，得到：

即分处于上下两个边界平面上的特征向量到平面 的距离为 ，同时，这也是上下两个边界平面之间的距离。

我们的目标就是让 最大化，即求

因为 ，作为分母不好进行运算，因此我们求 的等价问题：

因为求 的最大值就是求 的最小值，而 ，求 的最小值就是求 的最小值。为了简化后续的运算，这里不直接求 的最大值，而是求 的最大值。

而我们之前提过，对于样本向量集合中的样本，支持向量位于上下边界平面之上，非支持位于上下边界平面以外，即：

将上面的式子等号两边同乘以 并进行合并得到：

我们对问题重新梳理一遍。

首先将特征空间进行缩放后，存在一个平面 ，和上下两个边界平面 ，而这些样本向量则处于上下边界平面之上和上下边界平面之外，即：，我们的目标就是想让上下两个边界的距离（也是不同类别的支持向量到平面 的距离） 最大。

即在满足 的条件下，求得平面 ，使 最大，等同于求 。

转为数学语言来描述就是：

是满足于的意思，即：Subject to

求解

目标函数：

上述目标函数的求解，实际上是求解约束条件下的最优化问题，需要用到三个知识点：拉格朗日乘子、 KKT条件和对偶性，本文不对这三个知识点做过多的阐述，只使用其结论来帮助我们求解，感兴趣的可以自行查阅相关资料。

在这里给出不算严谨的结论：

对于不等式约束的凸优化问题，最优解满足 KKT 条件，具备强对偶性，可以通过求解其对偶问题来求解原问题。

目标函数的约束条件是不等式约束，因为样本集合 D 中的每一个样本向量都要满足 的条件，因此我们的约束条件的个数与样本的的数量是相同的（每一个样本都满足约束条件），共有 M 个约束条件。

因此我们为每个约束条件引入一个变量 作为拉格朗日乘子（），将原函数转换为拉格朗日函数：

注：拉格朗日函数中，，即 是样本向量，而我们的样本是已知的，因此每个 都是常数，不能算作变量。同样的， 是每个样本对应的标签，也是已知的，也不能算作变量。因此转化后的拉格朗日函数是只关于未知量 的函数。

为了简化书写，我们将约束条件 转化为函数表达：

此时，拉格朗日函数为：

好了，我们原本的目标是求 两个变量，现在引入了一个 ，成功地将问题变成了求三个变量 ，难度又增加了。

但是先别急，现在我们通过引入 ，将原本的约束条件添加进作为目标函数的一部分，形成了新的目标函数，这是一个含有不等式约束的凸二次规划问题，可以通过求解其对偶问题来解。

我们先来分析一下 的情况：

在样本向量集合中，对于支持向量，肯定满足 ，因此 ；

对于非支持向量，肯定满足 ，因此 。

因此，可以得出：

上式是什么意思呢？按照常理来说，通过 我们应该得到 ：，才对啊，为什么直接认定了等于 0 呢？

这里要注意的是， 是变量， 这个式子的含义是寻找使得 取最大值的 ，那么最简单的，只要让 全为 0， 自然等于 0。当然，这样取值没有意义，我们可以分析一下这种情况：

当然，在这里我们认为所有的样本都是线性可分的，即对于任意样本向量，都有 ，对于不可以线性可分的样本集，我们可以使用核函数或软间隔来划分，这部分在后面会陆续讲到。

假如数据集不是线性可分的，那么存在其中的部分不满足约束条件的越界样本，使得

对于这部分样本，，则导致 ，那么：

此时，我们得到了新的目标函数：

这个目标函数是怎么转化得到的呢？

我们首先来理解一下目标函数的含义：

这个式子的含义是先固定住 ，将其当作常数（如记作 ），那么目标函数就变成了关于 的函数，然后对该函数取最大值，记作 ：

这里得到的 实际上是一个关于 的函数，然后在求出最小的 ，即：

可能比较绕，举个例子来理解上面的过程。

从同一个年级的每个班里面，选取每个班里个子最高的学生，然后再从这些选出来的学生中，选出个子最小的。

如果用伪代码的形式来理解的话就是：

1
2
3
4
5
6
7
8
9
target
for w_i, b_i in w, b:
 max_a
 for a_i in a:
     if max_a < L(w_i, b_i, a_i):
         max_a < L(w_i, b_i, a_i)
 if target > max_a:
     target = max_a
return target

对于线性可分的样本集合，恒有 ，因此：

则：

这与原问题是一样的。

对于线性不可分的样本集合，则有：

故：

对于我们构造的拉格朗日函数，要想使用强对偶性来求解，需要满足两个条件：

1、优化问题是凸优化问题；

2、满足 KKT 条件。

在凸优化问题中，满足 KKT 条件的点是全局最小值点，在该点上满足强对偶性。而我们的目标函数 的 KKT 条件为：

满足上述 KKT 条件的驻点就是函数的全局最小值的解，在该驻点上满足强对偶性，因此我们得到原问题的对偶问题：

根据 KKT 条件可知：

将其代入拉格朗日函数 :

于是，在该驻点上，问题转换为了：

到达这一步转换，我们发现，本来是要求超平面 的未知数 的问题，现在变成了求 的问题了。

求出 后，再通过 KKT 条件中的 判断哪些 不为 0，其对应的样本向量就是支持向量，有了支持向量，就可以很方便得求出 和 。

则根据 KKT 条件中的条件可以计算出：

由此可知在 中，至少存在一个 （若 全为 0，则 ，即不存在能将线性可分的样本分割的平面，与已知矛盾），与之对应的支持向量为 ，对于支持向量 ，有：

然后就可以得到平面 了。

SMO 算法

SMO 算法与启发式算法想法相似，在面对多个变量的求解问题时选择固定一个变量，其他的视为常量，再一个个推导其他的变量。

根据式 的 KKT 条件，我们有约束条件：

如果我们根据上面的约束条件来直接固定 为变量， 为常量时，则可以直接推导出：，无法解决问题。因此我们稍微进行一下变化，一次固定两个变量 ，以此计算 的推导式。

根据式 ：

固定 为变量，其余的 为常量，则可以得到一个关于 两个变量的的函数：

因为 或 ，所以 ，我们令 ，将上式化简得到：

又因为 为常量，因此 为常数，我们令其为 ，则：

再根据我们式 的约束条件 ，可知：

等号两边同乘以 ，并令 得：

将式 代入式 来消掉 ：

化简得：

使用导数求 的极值：

化简得：

令导数 ：

接下来我们使用迭代的方法来求解 。

根据式 的约束条件：

可知 恒等于 ，所以在第 次和第 次迭代的时候，有：

而 恒为 或 ，在迭代的过程中不变，故：

这里是什么意思呢，也就是说在第 次迭代的时候，我们根据此时的 的值，来求第 次计算的最优解的 即可。

因此我们可以看到，在下面的计算中，我们会通过将第 次的 的值代入来计算第 次迭代的的 值。

将 中的 代入 ：

解释一下这样做的原因，在第 次计算时，我们可以认为此时的 的值是已知的，那么我们就通过 式中的偏导为 0 的情况，来求出在 次的条件下的 最大值点 （最优点），将其作为第 次迭代的新 ，以此来完成迭代。

换句话说， 在最开始的时候，我们给定 一个初始值，然后根据这一系列的初始值，将 视作变量，然后通过求偏导的方式，找到使得 的最大值点，然后将其视为下一次迭代的初始值。

现在想办法化简一下上面的式子。

在式 中，有：

那么令函数 ， 则有：

因为 是一个标量，且 ，所以 。

因为 的值是由 和 来决定的，而在每次迭代过程中， 和 的值都会更新，因此 的值也是不断迭代的。

因此在 次迭代中，有：

将式 代入 式 ：

根据式 ，可得：

即：

现在，我们得到了 的迭代关系：

而：

故可知最终的迭代公式：

在每次对 进行更新后，还需要同时对 进行更新。

根据对支持向量的定义，可以知道当 ，对应的 是支持向量，则有：

若在第 次迭代时求得 ，则：

而由 的定义可知：

因为在迭代的过程中 是保持不变的，因此 是个常数，也在迭代过程中保持不变，因此：

故化简为：

同理求得当 时：

故此可知，若 都大于 0 时，那么 与 均是有效的，取其和的一半即可；

若 同时为 0，即 则对应的 不是支持向量，不需要更新 的值。

于是，得到关于 的迭代更新公式：

软间隔

在实际的样本中，样本向量往往并不是那么恰巧可分的，或者说如果我们过于追求一定将样本分成两类，那么我们最后求的的超平面的间隔可能会比较小，发生过拟合，这样会导致我们求出的超平面的分类效果变差。

··· 待补充

因此我们在实际操作过程中，引入“软间隔”的概念。

软间隔这个概念是相对于硬间隔而言的。

所谓的硬间隔是指我们的超平面的上下边界能够严格将样本向量进行区分，而软间隔是指我们划分的超平面允许一部分样本在上下边界内，甚至允许一部分样本划分在错误的边界内，如下图：

···待补充

对于硬间隔，我们原本的约束条件是：

而对于软间隔，我们引入松弛变量 ，每一个样本都有一个对应的松弛变量，用来表示该样本不满足约束条件的程度。

对于 的含义该如何理解呢？

我们知道，支持向量其实一共决定了三个平面：

如下图：

08 空间缩放后上下边界距离示意图

如果我们想让上下边界往外移一点，让间隔变大，如下图：

09 软间隔

那么可能会导致部分样本在上下边界之内，也就是说，现在的间隔没能将全部样本划分到上下边界之外（如样本 a、b、d），暂且将这些在边界之内的样本称之为跑偏了（将样本分在错误的一边（如样本点 d）也是跑偏了的一种）。

那么 就是用来描述这些跑偏了的样本和它原来该呆在的边界的偏离程度，也就是说这些跑偏的样本离正确的上边界或下边界越远，这些样本对应的 越大。

我们的上下边界向外移动了，那么上下边界的距离还是 2 吗？

答案是：是的，上下边界的距离还是 2。

则是因为我们可以在移动上下边界之后，仍然可以使用缩放特征空间的方法，仍然可以使得上下边界的距离等于 2。

因此形成新的约束条件：

此时，我们的目标函数也就变成了：

其中， 是常数，是用来控制惩罚项的系数， 值越大，对错误分类的惩罚越大。

我们来分析一下引入的惩罚系数 会对最终结果产生的影响。

我们前面讲到了，将上下边界平面往外移动，会造成部分样本被包含在上下边界之内，而 衡量的是每个样本跨过其对应的上下边界的距离，而

也就是说当一个样本在其正确的上下边界之外时，那么它所对应的 ；当一个样本不在其正确的上下边界之外时，它对应的 ，且其大小与跑偏的距离成正比。

那么 就是所有跑偏的样本跑偏的程度的和。

当我们的上下边界间隔越小，跑偏的样本越少， 越小；上下边界间隔越大，跑偏的样本越多， 越大。

我们可以将 理解成损失，上下边界间隔越小，损失越小，否则损失越大。

那么 就是我们对损失的容忍程度， 越大，代表我们越不能容忍偏移，自然最后求出的上下边界间隔越小； 越小，代表我们越嫩容忍更多的偏移，自然最后求出的上下边界间隔越大。

因此，我们的拉格朗日函数变成了：

同样的，我们的目标函数变为了:

此时， 最小值点需要满足的 条件为：

在改点，利用强对偶性求对偶问题：

将 KKT 条件中的信息代入 的表达式得：

可以看到，这里化简出来的式子和硬间隔是一样的，只是多了个 约束条件。

因为 ，因此，而最终的化简式中并不包含 ，因此可以将 的约束条件化简为： 。

因此：

这里，只需要继续使用 SMO 算法来求解 即可。但是要注意的是，因为多了个 约束条件，我们在更新的时候需要注意不能打破约束边界。

也就是说当我们得到 的迭代公式（为了便于查看，这里写的是式 ，最终的迭代公式请见式 ）：

此时， 只是初步值，还要根据约束条件进行裁剪。

根据约束条件 ：

同时，我们根据 的约束条件，得到：

由约束条件可知， 只能落在 的区域内，又因为 ，其中 的值为 -1 或 1，因此：

即 在一条斜率为 1 或 -1 的直线上，因此 只能在 范围内部的线段上。

异号：

10 y 异号

同号：

11 y 同号

因此，根据迭代公式算出来的 还需要将其定位在上面描述的线段上，即对于 ，有上下边界：

其中， 是线段的两个端点，且：

因此，裁剪后的 的值应该是：

有了 ，即可求得 ：

同样的，求出 后也需要对 进行更新。

但是

由于引入了软间隔的概念，因此我们可以得到：

由硬间隔的 的迭代公式可知，当 时：

当 时：

当 同时满足 时，；当 都为 0 或 时，那么所有在 之间内的值都满足 KKT 条件，此时取两者的均值作为新的 。

于是得到引入软间隔后 的迭代公式：

核函数

前面提到，当面对线性不可分的问题时，解决的方法一种是引入软间隔，允许部分样本在上下边界之内或在错误的一边，另一种就是使用核函数。

对于在当前特征空间中不可线性划分的样本，我们可以通过将其映射到高维空间中来将其划分。

这里举个简单的例子：

一维空间中有四个点：

12 一维分割

如果说要在直线上找一个点，能将不同颜色的点分开，这显然是不可能的。

现在，我们给这个一维空间 再添加一个维度 ，形成了一个二维空间。那么原本的点在另一个维度上的值该是多少呢？

我们定义一个映射函数 ，并且点在 这一维度上： ，于是，在二维空间中，样本点的坐标就变成了：。

在新的二维空间中：

13 二维分割

现在就可以用一条直线将其划分开了。

值的注意的是，将原本的样本升维后，划分样本的分割平面也是在高维空间中的。

根据我们推导出的目标公式：

可以看到，公式里面关于样本 的计算只有 这是一个内积计算，结果是一个标量。因此面对线性不可分的问题，我们可以将样本向量映射到线性可分的高维空间中去，然后再计算内积。

如，假设我们有高维映射函数 能够将低维的向量 映射到高维空间中去（即： ），那么我们完全可以将其带入到目标函数中，得到新的目标函数：

为了简化书写，我们暂且将 记作 ，因此我们的映射的含义就是：

将原本需要计算 的问题变成了在更高维度上计算 的问题。

那么这里的映射是什么意思呢？

举个简单的例子，假如说将二维空间中的向量，映射到三维空间中：

这里只是一个很简单的例子，三维空间中的新增的一维是通过对二维空间中的不同维度相加得到的。

那么我们看一下高维映射后的内积的计算。

在没有进行映射之前，二维向量 的内积为：

而将其映射到三维空间后计算：

现在看来，对于线性不可分的问题，我们似乎已经通过将样本映射到更高维的空间中来变成线性可分的问题了。

但是事情真的解决了吗？如果原始数据的特征空间是二维的，我们可以将其映射到五维空间来进行计算，如果是三维的，可以将其映射到十九维，，但如果唯独比较大，几百维、几千维，那么映射到超级高的维度中进行计算，计算量太大了。而且如果遇到无穷维的情况，根本无法计算。

因此，核函数便应运而生。

假设 是一个从低维（ 维）的空间 到高维（ 维）空间 的映射，那么如果存在函数 ，对于任意的 都有：

那么对于函数 ，我们称之为核函数，或者核技巧。而对于 的计算，是在低维度下进行计算的，无需真的将向量映射到高维空间中去，因此大大的减少了计算量。

下面将结合高斯核函数来帮助理解核函数的作用。

高斯核函数：

其中， 是一个常数，用来缩放 的值而已。

为什么高斯核函数能够完成上面我们所说的事情呢？

因为 ，由于 为常数，令其为 ，则有：

我们都知道，对指数函数 可以进行泰勒展开：，因此将高斯核函数进行泰勒展开得到：

从这里可以看出，高斯核函数化简后，相当于是一个关于 的 阶多项式和，相当于将向量映射到了 维空间中。

使用核函数，就巧妙的将特征空间映射到更高维度的空间中，且内积计算仍然保留在低维空间中计算。

那么，核函数该在哪些地方使用呢？

我们可以将核函数简单理解成在高维空间中对两个向量相似度的衡量，其结果应该是一个标量。我们上面推导的公式里面最主要的就是对 的迭代的计算。而其推导中，涉及到两向量的运算的只有 ，因此答案就很明显了，在计算 的时候使用核函数就行了。

因此，我们可以将 由原本简单的内积计算替换成核函数计算即可，即： 。

代码实现

简单 SVM 实现

了解了 SVM 的原理和推导之后，接下来就能开始代码实现部分的内容了。

样本图像

首先，我们要创建要用到的样本：

1
2
3
4
5
6
7
8
9
10
11
import numpy as np
from sklearn.datasets import make_blobs

def CreateSamples(M):
    samples, labels = make_blobs(n_samples = M, centers = 2, cluster_std = 1.8, random_state = 6)
    for i in range(len(labels)):
        if labels[i] == 0:
            labels[i] = -1
    return samples, labels

samples, labels = CreateSamples(40)

这里我们定义了一个样本创建函数，该函数使用 sklearn.datasets 库 make_blobs 函数按照指定的样本数 M，创建了两类样本和标签，并将原标签 [0, 1] 处理成 [-1, 1] 的形式。值的注意的是，使用 make_blobs 函数创建出的每一类样本都有 M 个，换句话说，sample 的大小为 2M。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 从 0-N 中随机选择
def SelectIndexJRand(index_i, M):
    index_j = index_i
    while index_j == index_i:
        index_j = int(np.random.uniform(0, M))
    return index_j

# 裁剪 alpha
def ClipAlpha(alpha, H, L):
    if alpha > H:
        alpha = H
    if L > alpha:
        alpha = L
    return alpha

接下来，我们定义两个函数，SelectIndexJRand 函数用于根据样本数量 M，随机选取不等于 i 的下标。因为我们的 SMO 算法需要同时操作两个 ，这里就是使用 SelectIndexJRand 来选取另一个 的下标。

ClipAlpha 用于通过 H、L 来修建 的值，保证 。

接下来，就是 SVM 的实现代码了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def SmoSVM(samples, labels, C, toler, MaxIter):
    samples = np.mat(samples)
    labels = np.mat(labels).transpose()

    M, N = np.shape(samples)
    b = 0
    alphas = np.mat(np.zeros((M, 1)))
    
    iter = 0
    while iter < MaxIter:
        alphas_pairs_changed = 0
        for i in range(M):
            # h(x_i)
            h_x_i = float(np.multiply(alphas, labels).T * (samples * samples[i, :].T)) + b
            error_i = h_x_i - float(labels[i])	# h(x_i) - y_i 
            if (((labels[i] * error_i < -toler) and (alphas[i] < C)) or ((labels[i] * error_i > toler) and (alphas[i] > 0))):
                j = SelectIndexJRand(i, M)	# 随机选择另一个 alpha_j
                # h(x_j)
                h_x_j = float(np.multiply(alphas, labels).T * (samples * samples[j, :].T)) + b
                error_j = h_x_j - float(labels[j])	#  h(x_j) - y_j

                alpha_i_old = alphas[i].copy()
                alpha_j_old = alphas[j].copy()
                
                # 计算上下边界
                if labels[i] != labels[j]:
                    L = max(0, float(alphas[j] - alphas[i]))
                    H = min(C, C + float(alphas[j] - alphas[i]))
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, float(alphas[j]) + alphas[i])
                if L == H:
                    continue

                eta = -(samples[i, :] - samples[j, :]) * (samples[i, :] - samples[j, :]).T
                if eta >= 0:
                    continue

                # 更新 alpha_i、 alpha_j
                alphas[j] -= labels[j] * (error_i - error_j) / eta
                alphas[j] = ClipAlpha(alphas[j], H, L)
                if abs(alphas[j] - alpha_j_old) < 0.00001:
                    continue

                alphas[i] += labels[j] * labels[i] *(alpha_j_old - alphas[j])

                # 更新 b
                b_i = b - error_i- labels[i] * (alphas[i] - alpha_i_old) * samples[i, :] * samples[i, :].T - labels[j] * (alphas[j] - alpha_j_old) * samples[i, :] * samples[j, :].T
                b_j = b - error_j- labels[i] * (alphas[i] - alpha_i_old) * samples[i, :] * samples[j, :].T - labels[j] * (alphas[j] - alpha_j_old) * samples[j, :] * samples[j, :].T
                if 0 < alphas[i] and C > alphas[i]:
                    b = b_i
                elif 0 < alphas[j] and C > alphas[j]:
                    b = b_j
                else:
                    b = (b_i + b_j) / 2.0

                alphas_pairs_changed += 1

        if alphas_pairs_changed == 0:
            iter += 1
        else:
            iter = 0
    return b, alphas

首先，函数接收 5 个参数：

samples：样本

labels：样本对应的标签

C：惩罚系数

toler：容错率

MaxIter：最大迭代次数

代码的第 2-3 行，将输入的 samples、labels 转换为矩阵（numpy.mat），并且将 labels 进行转置，从而将 的矩阵（行向量）转换为 的矩阵（列向量）（M 为样本的个数，在这里，代表 samples 和原始的 labels 的行数）。此时，类别标签中的每行元素都和数据矩阵 samples 中的行一一对应

代码的第 5 行，定义并初始化 。

代码的第 5 行，获取 samples 矩阵的形状（ 为矩阵的行数，代表样本的个数； 为矩阵的列数，代表着样本特征空间的维度；在该示例中，）。

代码的第 6-7 行，定义并初始化了 $b = 0$ 和 矩阵，其大小为 ，初始值全为 0。

代码的第 9-10 行定义了记录迭代次数的变量 iter 并开始进行迭代。

代码的第 11 行，定义了 alphas_pairs_changed 变量，用于记录该次迭代是否对 进行了优化更新，若该次迭代对 进行了更新，则认为更新后的值为初始值，重新开始迭代。

从代码的第 12 行开始，遍历每一个 进行更新迭代（每次循环都固定住 并挑选出另一个 来进行计算，就是前面 SMO 算法中的 的概念）。通过遍历的方式，将 全部进行更新迭代。

在代码的第 14 行，定义了变量 hx_i 用于存储 的值，$$h(x_i) = w^Tx_i + b = \sum{j = 1}^{m} \alpha_j y_j x_j^T x_i + b$$ 。

我们知道，分隔超平面也可以叫做决策超平面。当分隔超平面建立完成后，则会有：

因为我们已经对 进行了初始化，并在 while 循环中不断优化 的值，而 ，因此相当于在不断优化分隔超平面。

上式的意思就是根据当前的分隔超平面，来预测样本 的类别。

在代码的第 15 行定义了变量 error_i，用于式 中对 的迭代更新使用。

变量 ，前面提到过 的结果相当于是对样本 所属类别的预测，而 是样本 的实际类别。

因此 实际上就是对样本 的预测值与实际值的误差。

在代码的第 16 行，通过判断 是否在误差范围内来决定是否继续进行计算并更新 。

这里解释一下这个判断的依据和原理是什么。

我们知道，SMO 算法中，主要是通过迭代来更新计算 的值，最后再通过大于 0 的 来分辨出支持向量，最后再计算出分隔超平面 中的 和 的值。

似乎我们只是在最后一步中，在计算出 的值之后才开始计算分隔超平面。

但其实不是这样的，SMO 算法中，看似只是通过初始的 的值来不断更新迭代 ，但是要注意的是，在第一次，也就是算法最开始的时候，我们不仅给定了 的初始值为 0，还给定了初始值 。而我们在计算过程中始终满足于式 (5) 的 KKT 条件，这些条件中有一条就是 ，也就是 。

我们给定了 的初始值，其实也就是给定了 的初始值和 的初始值，也就是给定了分隔超平面 。

只不过最开始的超平面只是 是特征空间的原点而已。我们在后续的不断更新 和 ，其实也是在不断调整分隔超平面 。

理解了这个，就能明白这个判断条件和真正的 SMO 算法做的事情了。

当我们进入这次迭代的时候，我们暂且称这一次迭代为第 次迭代，在这一次迭代中，我们要做的事情是根据将 次迭代后的 作为初始值（暂且记为 ，那么第 次的超平面记作 ）通过迭代计算出 （也计算出了 ）。

那么在第 次迭代的时候，我们有了 ，我们知道，对于分割平面还有与之平行的上下边界，对于所有的样本，应该有:

也就是说，所有的样本都应该位于上下边界之上或者之外。我们暂且将 称之为样本 到分隔超平面的距离（这个距离并不是严格意义上的样本向量到平面的距离，而是为了方便理解，将其称作距离而已），那么上式的含义就可与理解成：对于所有的样本，到分隔超平面的距离大于等于 1。

对于距离 我们记作 。

而由于我们引入了软间隔这个概念，因此有了：

接下来，我们结合代码来看一下这个条件判断：

1
if (((labels[i] * error_i < -toler) and (alphas[i] < C)) or ((labels[i] * error_i > toler) and (alphas[i] > 0))):

这行代码的条件中其实包含了四个条件的判断。

我们现看 or 条件两边的条件判断：

1
(labels[i] * error_i < -toler) and (alphas[i] < C)

前面讲过 labels[i] * error_i 其实就是:

那么上式就可以表示为：

我们先看 ，这是引入软间隔后的 KKT 条件之一，我们初始化的时候，将所有的 都初始化为 0，在进行迭代更新的时候，也会对 进行裁剪，因此 是始终满足的。

对于 ，我们知道，通过 KKT 条件可知 ，而根据 的条件判断，当样本不满足该条 KKT 条件时， 成立。也就是说，对于样本 ，不满足当前的分隔超平面 的 KKT 条件，因此需要更新迭代超平面的位置。

接下来在看另一半条件：

1
(labels[i] * error_i > toler) and (alphas[i] > 0)

和上面一样的理解，上式可以表示为：

现看 ，我们知道， 代表着其对应的样本 为支持向量，那么支持向量应该处于软间隔之内，即 ，而 成立的时候，显然说明支持向量在软间隔之外，即 不满足作为支持向量应该具备的条件，因此需要更新 。

这两个条件合在一起就是相当于判断该样本点是否满足当前的分隔超平面 的 KKT 条件，若不满足，就更新 。

当然，对于这个条件判断，我们也可以有另一种理解：

是样本 的预测类别与实际类别的误差，在 if 判断语句中，不管是正间隔还是负间隔都会被测试。同时也要检查 的值，以保证其不能等于 0 或 （因为每次更新 时都会对其进行裁剪，小于 0 或 大于 时将被裁剪为 0 或 ）。若 等于 0 或 ，表明 已经在 的边界上了，不能对其再见效或增大，因此也就不值得再对它们进行优化了。

代码的第 17 行，通过随机的方式选择了与 不同的下标 。

在 19-20 行，同样的，分别计算 和 。

在代码的 22-23 行，分别将 保存在变量 alpha_i_old 和 alpha_j_old 中。

在代码的 26-33 行，通过 来计算 （详见式 ）。

代码第 35 行，计算 :

详见式 中的迭代公式中的分母。

代码的第 40 行：

1
alphas[j] -= labels[j] * (error_i - error_j) / eta

label[j] =

error_i =

error_j =

error_i - error_j =

eta =

即是通过式 来计算 的值。

接下来，代码的第 41 行，对新计算出的 进行裁剪。

代码的第 43-44 行，通过判断该次更新的幅度，来决定这次更新是不是一次有效的更新，从而决定是否要更新剩下的 。

代码的第 45 行，利用式 来更新 ：

代码的 47-55 行，就是利用式 来更新 的值。

代码的 57 行用于记录该次 for 循环是否对 这一对和 进行了改变。

代码的 59-62 行，会判断在该次 while 循环中，是否对 的值有过更新，如果有更新则将 iter 迭代计数器设为 0 之后继续执行程序。只有在所有数据集上遍历 maxiter 次，且不再发生任何 的更新之后，程序才会停止并退出 while 循环。

将上面的函数转化为伪代码的话大致如下：

1
2
3
4
5
6
7
初始化 alpha、b
while iter < max_iter:
    for x_i in samples:
        if x_i 不满足当前的分隔超平面的 KKT 条件:
            更新分隔超平面
    if 分隔超平面被更新了:
    	重置迭代次数 iter

接下来，就是使用 SmoSVM 函数计算出的 来计算 ：

1
2
3
4
5
6
7
8
def CalculateW(alphas, samples, labels):
    samples = np.mat(samples) 
    labels = np.mat(labels).transpose()
    M, N = np.shape(samples)
    w = np.zeros((N,1))
    for i in range(M):
        w += np.multiply(alphas[i] * labels[i], samples[i,:].T)
    return w

下面，我们对代码进行测试：

1
2
3
4
5
6
7
8
b, alphas = SmoSVM(samples, labels, 0.01, 0.001, 40)
w = CalculateW(alphas, samples, labels)
for i in range(len(alphas)): 
    if alphas[i] > 0:
        print(i, end = ',')
print()
print("b = ",b)
print("w = ",w)

image-20230221112251082

我们在将样本和分隔超平面画出来如下图所示：

分割样本图像

上图中，形状为五角星的即为支持向量，蓝色分割线将样本分为两部分。

完整代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
import numpy as np
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt

# M ：每一类的样本数
def CreateSamples(M):
    samples, labels = make_blobs(n_samples = M, centers = 2, cluster_std = 1.8, random_state = 6)
    for i in range(len(labels)):
        if labels[i] == 0:
            labels[i] = -1
    return samples, labels



def DrawSplitedSamples(samples, labels, w, b):
    # 画出支持向量
    support_vector = []
    support_vector_label = []

    normal_vector = []
    normal_vector_label = []
    for i in range(len(alphas)):
        if alphas[i] > 0:
            support_vector.append(samples[i])
            support_vector_label.append(labels[i])
        else:
            normal_vector.append(samples[i])
            normal_vector_label.append(labels[i])

    support_vector = np.array(support_vector)
    normal_vector = np.array(normal_vector)
    plt.scatter(normal_vector[:,0], normal_vector[:,1], c = normal_vector_label, cmap='summer')
    plt.scatter(support_vector[:,0], support_vector[:,1], c = support_vector_label, marker='*', cmap='summer')

    # 画出分割线
    ax = plt.gca()
    xlim = ax.get_xlim()    # 获取数据点x坐标的最大值和最小值
    # 根据坐标轴生成等差数列
    xx = np.linspace(xlim[0], xlim[1],30)
    yy = - (b + w[0] * xx) / w[1]
    yy = yy.T
    plt.plot(xx, yy)

    plt.savefig("分割样本图像.svg",format='svg',dpi=15, transparent=True,)
    plt.show()

# 从 0-N 中随机选择
def SelectIndexJRand(index_i, M):
    index_j = index_i
    while index_j == index_i:
        index_j = int(np.random.uniform(0, M))
    return index_j

# 裁剪 alpha
def ClipAlpha(alpha, H, L):
    if alpha > H:
        alpha = H
    if L > alpha:
        alpha = L
    return alpha

def SmoSVM(samples, labels, C, toler, MaxIter):
    samples = np.mat(samples)
    labels = np.mat(labels).transpose()

    b = 0

    M, N = np.shape(samples)
    alphas = np.mat(np.zeros((M, 1)))

    iter = 0

    while iter < MaxIter:
        alphas_pairs_changed = 0
        for i in range(M):
            # h(x_i)
            h_x_i = float(np.multiply(alphas, labels).T * (samples * samples[i, :].T)) + b
            # h(x_i) - y_i 
            error_i = h_x_i - float(labels[i])

            if (((labels[i] * error_i < -toler) and (alphas[i] < C)) or ((labels[i] * error_i > toler) and (alphas[i] > 0))):

                # 随机选择另一个 alpha_j
                j = SelectIndexJRand(i, M)
                # h(x_j)
                h_x_j = float(np.multiply(alphas, labels).T * (samples * samples[j, :].T)) + b
                
                # y_j - h(x_j)
                error_j = h_x_j - float(labels[j])

                alpha_i_old = alphas[i].copy()
                alpha_j_old = alphas[j].copy()
                
                # 计算上下边界
                if labels[i] != labels[j]:
                    L = max(0, float(alphas[j] - alphas[i]))
                    H = min(C, C + float(alphas[j] - alphas[i]))
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, float(alphas[j]) + alphas[i])
                if L == H:
                    continue

                eta = -(samples[i, :] - samples[j, :]) * (samples[i, :] - samples[j, :]).T
                if eta >= 0:
                    continue

                # 更新 alpha_i、 alpha_j
                alphas[j] -= labels[j] * (error_i - error_j) / eta
                alphas[j] = ClipAlpha(alphas[j], H, L)
                if abs(alphas[j] - alpha_j_old) < 0.00001:
                    continue

                alphas[i] += labels[j] * labels[i] *(alpha_j_old - alphas[j])

                # 更新 b
                b_i = b - error_i- labels[i] * (alphas[i] - alpha_i_old) * samples[i, :] * samples[i, :].T - labels[j] * (alphas[j] - alpha_j_old) * samples[i, :] * samples[j, :].T
                b_j = b - error_j- labels[i] * (alphas[i] - alpha_i_old) * samples[i, :] * samples[j, :].T - labels[j] * (alphas[j] - alpha_j_old) * samples[j, :] * samples[j, :].T

                if 0 < alphas[i] and C > alphas[i]:
                    b = b_i
                elif 0 < alphas[j] and C > alphas[j]:
                    b = b_j
                else:
                    b = (b_i + b_j) / 2.0

                alphas_pairs_changed += 1

        if alphas_pairs_changed == 0:
            iter += 1
        else:
            iter = 0
    return b, alphas

# 计算 w
def CalculateW(alphas, samples, labels):
    samples = np.mat(samples) 
    labels = np.mat(labels).transpose()
    M, N = np.shape(samples)
    w = np.zeros((N,1))
    for i in range(M):
        w += np.multiply(alphas[i] * labels[i], samples[i,:].T)
    return w

samples, labels = CreateSamples(40)
b, alphas = SmoSVM(samples, labels, 0.01, 0.8, 40)
w = CalculateW(alphas, samples, labels)
DrawSplitedSamples(samples, labels, w, b)

使用核函数的 SVM 实现

在前面我们提到了，对于线性不可分的问题，可以使用核函数将其转换为高维度特征空间中的线性可分问题。

我们已经实现了简单的 SVM 代码，现在将核函数应用到我们的代码之中。

在应用核函数之前，我们首先要知道，该在哪里添加使用核函数来进行计算。

前面讲过，将推导式中的 替换成 核函数来计算即可。

那么我们的代码中，出现 的地方有哪些呢？

答案是：

1
2
3
4
5
eta = -(samples[i, :] - samples[j, :]) * (samples[i, :] - samples[j, :]).T
···
b_i = b - error_i- labels[i] * (alphas[i] - alpha_i_old) * samples[i, :] * samples[i, :].T - labels[j] * (alphas[j] - alpha_j_old) * samples[i, :] * samples[j, :].T
···
b_j = b - error_j- labels[i] * (alphas[i] - alpha_i_old) * samples[i, :] * samples[j, :].T - labels[j] * (alphas[j] - alpha_j_old) * samples[j, :] * samples[j, :].T

我们将 eat 的计算展开，得到：

1
eta = 2 * samples[i, :] * samples[j, :] - samples[i, :] * samples[i, :].T - samples[j, :] * samples[j, :].T

为了加速计算，避免在循环过程中造成重复计算，我们可以用一个 大小的矩阵 kernel_mat 提前保存核函数的结果，即：

核函数计算函数：

1
2
3
4
5
6
7
8
9
def GetKernelMat(samples, gamma):
    M, N = np.shape(samples)
    kernel_mat = np.mat(np.zeros((M, M)))
    for i in range(M):
        for j in range(M):
            delta = (samples[i, :] - samples[j, :])
            k = np.exp(-gamma / 2 * delta * delta.T)[0, 0]
            kernel_mat[i, j] = k
    return kernel_mat

带有核函数的 SVM 代码实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
'''
Descripttion: 
version: 
Author: xudapeng
Date: 2023-02-20 10:50:17
'''
import numpy as np
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt

# M ：每一类的样本数
def CreateSamples(M):
    samples, labels = make_blobs(n_samples = M, centers = 2, cluster_std = 1.8, random_state = 6)
    for i in range(len(labels)):
        if labels[i] == 0:
            labels[i] = -1

    return samples, labels

def DrawSplitedSamples(samples, labels, w, b):
    # 画出支持向量
    support_vector = []
    support_vector_label = []

    normal_vector = []
    normal_vector_label = []
    for i in range(len(alphas)):
        if alphas[i] > 0:
            support_vector.append(samples[i])
            support_vector_label.append(labels[i])
        else:
            normal_vector.append(samples[i])
            normal_vector_label.append(labels[i])

    support_vector = np.array(support_vector)
    normal_vector = np.array(normal_vector)
    plt.scatter(normal_vector[:,0], normal_vector[:,1], c = normal_vector_label, cmap='summer')
    plt.scatter(support_vector[:,0], support_vector[:,1], c = support_vector_label, marker='*', cmap='summer')

    # 画出分割线
    ax = plt.gca()
    xlim = ax.get_xlim()    # 获取数据点x坐标的最大值和最小值
    # 根据坐标轴生成等差数列
    xx = np.linspace(xlim[0], xlim[1],30)
    yy = - (b + w[0] * xx) / w[1]
    yy = yy.T
    plt.plot(xx, yy)
    plt.show()

# 从 0-N 中随机选择
def SelectIndexJRand(index_i, M):
    index_j = index_i
    while index_j == index_i:
        index_j = int(np.random.uniform(0, M))
    return index_j

# 裁剪 alpha
def ClipAlpha(alpha, H, L):
    if alpha > H:
        alpha = H
    if L > alpha:
        alpha = L
    return alpha

def GetKernelMat(samples, gamma):
    M, N = np.shape(samples)
    kernel_mat = np.mat(np.zeros((M, M)))

    for i in range(M):
        for j in range(M):
            delta = (samples[i, :] - samples[j, :])
            k = np.exp(-gamma / 2 * delta * delta.T)[0, 0]
            kernel_mat[i, j] = k

    return kernel_mat

def KernelSmoSVM(samples, labels, C, toler, MaxIter):
    samples = np.mat(samples)
    labels = np.mat(labels).transpose()

    b = 0

    M, N = np.shape(samples)
    alphas = np.mat(np.zeros((M, 1)))
    kernel_mat = GetKernelMat(samples, 1.3)

    iter = 0

    while iter < MaxIter:
        alphas_pairs_changed = 0
        for i in range(M):
            # h(x_i)
            h_x_i = float(np.multiply(alphas, labels).T * (samples * samples[i, :].T)) + b
            # h(x_i) - y_i 
            error_i = h_x_i - float(labels[i])

            if (((labels[i] * error_i < -toler) and (alphas[i] < C)) or ((labels[i] * error_i > toler) and (alphas[i] > 0))):

                # 随机选择另一个 alpha_j
                j = SelectIndexJRand(i, M)
                # h(x_j)
                h_x_j = float(np.multiply(alphas, labels).T * (samples * samples[j, :].T)) + b
                
                # y_j - h(x_j)
                error_j = h_x_j - float(labels[j])

                alpha_i_old = alphas[i].copy()
                alpha_j_old = alphas[j].copy()
                
                # 计算上下边界
                if labels[i] != labels[j]:
                    L = max(0, float(alphas[j] - alphas[i]))
                    H = min(C, C + float(alphas[j] - alphas[i]))
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, float(alphas[j]) + alphas[i])

                if L == H:
                    continue

                eta = 2 * kernel_mat[i, j] - kernel_mat[i, i] - kernel_mat[j, j]
                if eta >= 0:
                    continue

                # 更新 alpha_i、 alpha_j
                alphas[j] -= labels[j] * (error_i - error_j) / eta
                alphas[j] = ClipAlpha(alphas[j], H, L)
                if abs(alphas[j] - alpha_j_old) < 0.00001:
                    continue

                alphas[i] += labels[j] * labels[i] *(alpha_j_old - alphas[j])

                # 更新 b
                b_i = b - error_i- labels[i] * (alphas[i] - alpha_i_old) * kernel_mat[i, i] - labels[j] * (alphas[j] - alpha_j_old) * kernel_mat[i, j]
                b_j = b - error_j- labels[i] * (alphas[i] - alpha_i_old) * kernel_mat[i, j] - labels[j] * (alphas[j] - alpha_j_old) * kernel_mat[j, j]

                if 0 < alphas[i] and C > alphas[i]:
                    b = b_i
                elif 0 < alphas[j] and C > alphas[j]:
                    b = b_j
                else:
                    b = (b_i + b_j) / 2.0

                alphas_pairs_changed += 1

        if alphas_pairs_changed == 0:
            iter += 1
        else:
            iter = 0
    return b, alphas

# 计算 w
def CalculateW(alphas, samples, labels):
    samples = np.mat(samples) 
    labels = np.mat(labels).transpose()
    M, N = np.shape(samples)
    w = np.zeros((N,1))
    for i in range(M):
        w += np.multiply(alphas[i] * labels[i], samples[i,:].T)
    return w

samples, labels = CreateSamples(40)

b, alphas = KernelSmoSVM(samples, labels, 0.01, 0.8, 40)
w = CalculateW(alphas, samples, labels)
DrawSplitedSamples(samples, labels, w, b)
-------------本文结束 感谢阅读-------------
打赏一包辣条