迈过三重门——详解SVM及其Python实现

来源：互联网发布：java银行项目经验编辑：程序博客网时间：2024/05/13 13:54

本文主要基于李航《统计学习方法》与周志华《机器学习》完成，加入了若干个人推导与注解，文后附Python3源码。跟我推导完，相信你一定会有收获。

初识SVM

支持向量机（Support Vector Machine,SVM）是ML&PR领域具有相当分量的一个算法，经常用在数据的分类与回归分析中，属于监督式、非概率、二分类器。当我们给定一组训练实例，每个实例被贴标签为二类中的一个（如±1），经过训练后，SVM就可以给出分类决策函数f(x)。当我们想要识别一个新样本，就把它扔给SVM得到的决策函数f(x)，从而得到一个判别结果。其可以广泛应用于文本、图像分类、生物信息学等领域。以上是SVM的简介和用途。接下来，进入干货部分。

第一重·线性硬间隔支持向量机

我们先来看看什么是线性可分数据的定义。如果我们给定一个数据集

T = {(x 1, y 1), (x 2, y 2), . . ., (x m, y m)}

其中，

x i ϵ R n, y i ϵ {+ 1, - 1}, i = 1, 2, . . ., m

如果存在一个超平面S

w \cdot x + b = 0

使得数据集的所有正实例点和负实例点被完全正确地划分到超平面两侧，那么称数据集T为线性可分数据集。这里

w=(w1,w2,...,wn)代表超平面的法向量，

b代表截距,为一数值，每一个实例

xi是

n维列向量，其实这里取用行向量还是列向量是无所谓的，矩阵运算自洽即可，周志华《机器学习》中对超平面采用的形式为

w T \cdot x + b = 0

本文采用

w⋅x+b=0的超平面形式。有人可能会问超平面怎么理解，我们以一个1-D数据集为例，如果正负样例线性可分，那么在数轴上一个点便可分开两类实例；对一个2-D的数据集，我们可以在一张平面上画出该数据集，如果线性可分，那么找一条直线便可分开。同理，可拓展至高维。也即，超平面是一比数据维度低一维的存在，是点、线、平面的推广。那么，是不是找到一个超平面就可以了呢？看看下面这幅图。对于这样一组分离的如此任性的样例，我们完全找到无数个超平面来把它们拆分开。

这时候怎么办，随便取一个超平面，确实可以实现对训练数据100%的accuracy，但是这并不能代表是最好的训练模型，如果测试数据在训练数据的基础上有一些偶尔冲动的噪声，那么将出现测试误差。图中最粗的这条线看起来对未来样例出现扰动时候的容忍性最好，即泛化能力最好。支持向量机采用间隔最大化方法来选择这个最优超平面，并且在这时，解是唯一的。我们来看看什么是间隔。试想一下，一个点到分类超平面的远近可以来评估我们预测的值得相信的程度，点距离分类超平面越远，分类结果越值得相信。而(w⋅xi+b)与yi的符号是否一致则代表了分类是否正确。那么，可以定义函数间隔为

γ ˆ i = y i (w \cdot x i + b)

遗憾的是，对于一个由

(w,b)确定的超平面，当我们将两个参数都扩大二倍，超平面没有发生改变，但函数间隔增大了二倍，这显然不是我们想要的结果。那么我们来看一下更直观的几何角度的间隔。努力回忆一下中学学过的解析几何，点到直线距离怎么求？对一个常规表示的直线方程：

A x + B y + C = 0

,点

(x0,y0)到这条直线的距离为：

A x 0 + B y 0 + C A 2 + B 2 - - - - - - - \sqrt

用我们前面定义的超平面（直线，假设数据为2-D）：

w \cdot x + b = 0

那么样例

xi到其距离为：

γ i = y i (w \cdot x i + b ∥ w ∥)

上式前面乘上

yi保证了距离是正数。括号里的内容，分母是二范数，对于一个列向量而言，其二范数

∥x∥=∑x2i−−−−√=xTx−−−√，看看是不是前面点到直线的定义？接下来，定义整个数据集的几何间隔为所有样本点关于超平面的几何间隔最小值

γ = m i n i = 1, 2, . . ., m γ i

函数间隔与几何间隔的区别，就在于后者分母上多了一个

w的二范数，即

γi=γˆi∥w∥。正是这个分母的存在，使得当我们将

(w,b)两个参数扩大二倍，函数间隔也随之扩大了二倍，而几何间隔并不会。理所应当地，几何间隔成为了我们希望优化的目标。约束条件为：所有样本到超平面的距离不小于该几何间隔。具体可写成下式：

m a x w, b γ

s . t . y i (w x i + b ∥ w ∥) \geq γ, i = 1, 2, . . ., m

考虑到几何间隔与函数间隔的关系，写为下式：

m a x w, b γ ˆ ∥ w ∥

s . t . y i (w x i + b) \geq γ ˆ, i = 1, 2, . . ., m

又，前面说过，函数间隔的取值对优化问题没有影响，那么我们取

γˆ=1，将优化目标等价修改，可得到下式：

m i n w, b 1 2 ∥ w ∥ 2

s . t . y i (w x i + b - 1) \geq 0, i = 1, 2, . . ., m

OK，说到这里，我们终于完成了支持向量机的最基础思想。我们定义上式中使得等号成立的点为支持向量。关于这一点，我们还可以通过取

γˆ=1推出几何间隔为

γ=1∥w∥来得到。支持向量机中的支持向量一定存在，或者说，支持向量不是存不存在的问题，正是这些支持向量决定了支持向量机。当我们修改其他实例点，甚至去掉其他点，模型不会发生改变。因此，现在，对于支持向量机这五个字，读者是否有了更加清晰的认识？对于上式这种最优化问题，常常利用拉格朗日对偶性将其原始问题转化为对偶问题来进行求解。引入m个拉格朗日乘子

αi,i=1,2,...,m，定义拉格朗日函数：

L (w, b, α) = 1 2 ∥ w ∥ 2 - \sum α i y i (w \cdot x i + b) + \sum a i

注意到，我们要求的是拉格朗日函数对

w、b的极小，将

L函数分别对

w、b求偏导，并令其为0，可得

w = \sum α i y i x i

0 = \sum α i y i

将这两个式子带入到拉格朗日函数中，可得

L (w, b, α) = 1 2 ∥ w ∥ 2 - \sum α i y i (w \cdot x i + b) + \sum a i

= 1 2 \sum \sum α i α j y i y j x T i x j - \sum \sum α i α j y i y j x T i x j + \sum a i

= \sum a i - 1 2 \sum \sum α i α j y i y j x T i x j

根据拉格朗日对偶性，求拉格朗日函数对

w、b的极小,等同于求该函数对

α的极大。于是，我们终于推出了原始问题的对偶问题：

m a x α \sum a i - 1 2 \sum \sum α i α j y i y j x T i x j

s . t . \sum α i y i = 0

α i \geq 0, i = 1, 2, . . ., m

很好，我们对原始问题求出了它的对偶问题。由于我们的优化是强对偶的（主问题为凸优化问题，且可行域中至少有一点使得不等式约束严格成立），那么如果我们求出了对偶问题的解，就可以得到原始问题的解。根据拉格朗日对偶性，这句话还有一个充分必要条件。充分必要条件是拉格朗日函数的几个参数满足KKT条件。这里我觉得很有意思，我们先来看KKT条件在我们写出的拉格朗日函数下的具体表示：

▽ w L = 0 ▽ b L = 0 ▽ α L = 0 α i \geq 0 α i (y i (w x i + b) - 1) = 0 y i (w x i + b) - 1 \geq 0

当然，

i还是我们熟知的样本序号，

i=1,2,...,m。以上六个要同时成立。前三个，很显然了，前两个是

L对

w、b的偏导求取结果，第三个是对

α，代表了原始问题和对偶问题的优化目标。来看第四个到第六个式子。

αi我们定义过，大于等于0。那么对任一样例，其对应的

αi要么为0，要么大于0。当

αi大于0，那么

yi(wxi+b)−1)=0，这是SVM很多重要性质的来源，比如支持向量。接下来看

αi=0的情况，可以看到，其对

yi(wxi+b)−1)的结果没有任何要求，只要求样本被正确归类即可，而

αi=0代表了该i值对应的样例不会出现在优化目标中，也就是说，这个样本是一个没有存在感的样本，取多大，有没有，对SVM模型没有影响。

我们来看KKT的充分性（拉格朗日函数的几个参数满足KKT条件→可以从原始问题的解求对偶问题的解）：取一个样例，其对应的αi>0，那么根据KKT条件，yi(wxi+b)−1)=0，左右两边同乘yi，可得yj((wxj+b)−yj)=0，进而，b=yj−wxj。而根据之前的求偏导，

w = \sum α i y i x i

代入可得

b = y j - \sum α i y i (x i \cdot x j)

好，到这里，我们讲完了线性硬间隔支持向量机的原理及其求解。来梳理一下，我们通过拉格朗日对偶性，将最原始的最小化几何间隔：

m i n w, b 1 2 ∥ w ∥ 2

s . t . y i (w x i + b - 1 \geq 0), i = 1, 2, . . ., m

转化为

m a x α \sum a i - 1 2 \sum \sum α i α j y i y j x T i x j

s . t . \sum α i y i = 0

α i \geq 0, i = 1, 2, . . ., m

关于

αi的求解，我最后会讲到。随后，通过

w = \sum α i y i x i

b = y j - \sum α i y i (x i \cdot x j)

来得到超平面的

w、b值，进而得到SVM训练模型。除此之外，我们通过KKT条件，深入了解了支持向量的内涵及SVM的本质，这对我们理解SVM有很大帮助。接下来，我们将在线性硬间隔支持向量机的基础上，继续深入一小步。

第二重·线性软间隔支持向量机

在前面的讨论中，我们始终假设样例在特征空间中是线性可分的。在有些情况下，训练样例中大部分是线性可分的，偶有一些特异点，那么我们引入软间隔支持向量机的意义，就是允许我们的模型在这些特异点上出错。相比于硬间隔支持向量机所要求的所有样本都要划分正确，软间隔支持向量机对样本的判别结果给予一些宽容，那么在某种程度上便可以实现对非线性样例的判别。所谓线性不可分，即是某些样例不能满足我们设定的函数间隔大于等于1的约束条件。为了soften这个间隔，我们给它加一个松弛变量，使得间隔变为弹性的，那么我们之前的约束条件就变为：

y i (w \cdot x i + b) \geq 1 - ξ i

由于松弛变量的存在，我们的优化目标函数要修改一下。我们仍然希望间隔尽量大，但是由于约束条件允许了一些样例犯错，那么我们在优化目标里，要为每个松弛变量都支付一个代价，从而在优化的过程中保证犯错的样例个数尽量少：

1 2 ∥ w ∥ 2 + C \sum ξ i

这里的

C>0为惩罚参数，C越大，对误分类的惩罚也就越大，一般由具体的应用问题来决定，或者看你的心情（和经验）。这种形式其实很类似于正则化了，一项为间隔大小，一项代表训练集上的误差。笔者之前接触过一些LASSO和ElasticNet等正则化方法，后面有时间会说一说这方面的内容。基于上面的思路，我们的线性软间隔支持向量机变为如下优化问题：

m i n w, b 1 2 ∥ w ∥ 2 + C \sum ξ i s . t . y i (w x i + b \geq 1 - ξ i) ξ i \geq 0 i = 1, 2, . . ., m

祭出神器拉格朗日函数的时候到了。将上式写成：

L (w, b, ξ, α, μ) = 1 2 ∥ w ∥ 2 + C \sum ξ i + \sum α i (1 - ξ i - y i (w x i + b)) - \sum μ i ξ i

其中，

αi、μi是拉格朗日乘子。令

L对

w,b,ξ的偏导为0，可得：

w = \sum α i y i x i 0 = \sum α i y i C = α i + μ i

将上式代入拉格朗日函数中，可得

L (w, b, ξ, α, μ) = 1 2 ∥ w ∥ 2 + C \sum ξ i + \sum α i (1 - ξ i - y i (w x i + b)) - \sum μ i ξ i = 1 2 \sum \sum α i α j y i y j x i x j + (α i + α j) \sum ξ i + \sum α i - \sum α i ξ i - \sum μ i ξ i - \sum \sum α i α j y i y j x i x j = \sum α i - 1 2 \sum \sum α i α j y i y j x i x j

根据拉格朗日对偶性，求拉格朗日函数对

ξ、b的极小等同于求其对

α的极大。于是可以将线性软间隔支持向量机的原始问题写为：

m a x α \sum α i - 1 2 \sum \sum α i α j y i y j x i x j s . t . \sum α i y i = 0 C - α i - μ i = 0 α i \geq 0 μ i \geq 0, i = 1, 2, . . ., m

后三个约束条件可以写成

0 \leq α i \leq C

综上，原始问题转化为对偶问题：

m a x α \sum α i - 1 2 \sum \sum α i α j y i y j x i x j s . t . \sum α i y i = 0 0 \leq α i \leq C ， i = 1, 2, . . ., m

可以证明，

w的解是唯一的，但

b的解可能不唯一。相比于前面讲的线性硬间隔支持向量机，

w的求解没有发生变化，而

b发生了变化，其求解公式我们会在后面的SMO算法中讲到。

第三重·非线性支持向量机

常年混迹于机器学习圈的大佬一定听过核函数这个（paper灌水）神器，直白理解，核函数采用一种非线性变换，将原空间的线性不可分的数据映射到新空间，然后在新空间中采用线性方法完成模型的训练。我们定义ϕ(x)为输入空间到特征空间的映射，那么如果存在一个函数K(x,z)满足

K (x, z) = ϕ (x) \cdot ϕ (z)

那么我们称

K(x,z)为核函数，

ϕ(x)⋅ϕ(z)为二者内积。这里藏了一个trick，我们通常不显式地定义映射函数

ϕ，而是直接定义核函数

K(x,z)，为什么这么搞呢？设想一下，我们如果先将两个样例分别映射到高维（甚至无穷维）的特征空间中，再去计算高维空间的内积，计算量可想而知。通常，我们定义的核函数有线性核、高斯核等，篇幅所限，读者可自行wiki相关关键字。在核函数的使用上，有一些基本经验，比如，对文本数据通常可采用线性核，情况不明时可先尝试高斯核等等。有了核函数，我们可以将线性不可分的数据映射到高维空间，在高维空间内应用线性支持向量机来完成模型的训练。注意到，在前面讲过的两种线性支持向量机中，我们的目标函数都只涉及到输入实例与实例的内积，那么，我们可以直接将其目标函数中的内积

xi⋅xj用核函数

K(xi,xj)=ϕ(xi)⋅ϕ(xj)来代替。那么显而易见的，我们的非线性支持向量机要解决的最优化问题如下：

m a x α \sum α i - 1 2 \sum \sum α i α j y i y j K (x i, x j ） s . t . \sum α i y i = 0 0 \leq α i \leq C ， i = 1, 2, . . ., m

将目标函数取反，求最小值。如下：

m i n α 1 2 \sum \sum α i α j y i y j K (x i, x j) - \sum α i s . t . \sum α i y i = 0 0 \leq α i \leq C ， i = 1, 2, . . ., m

迈门利器·序列最小最优化算法

OK了，介绍完三种支持向量机模型，我们发现所有的问题都剑指α的求解。这里我们来介绍一种高效求解凸二次规划问题的算法——Sequential Minimal Optimization(SMO)算法，1998年由Platt提出。观察我们的求解目标：

m i n α 1 2 \sum \sum α i α j y i y j K (x i, x j) - \sum α i s . t . \sum α i y i = 0 0 \leq α i \leq C ， i = 1, 2, . . ., m

问题的变量是拉格朗日乘子，也就是要求出样例个数个

α值。SMO算法的基本思想是固定

αi之外的所有参数，通过对

αi求偏导取极值，来更新

αi、b值，最后计算

w。由于约束

∑αiyi=0的存在，如果固定

α1、α2之外的所有变量，那么可推出：

α 1 = - y 1 \sum i = 2 m α i y i

如果

α2确定，那么

α1也随之确定，所以SMO同时更新两个变量。当固定其他变量

αi,i=3,4,...,m后，略去不含

α1、α2的常数项，SMO的最优化问题转化为

m i n α 1, α 2 W (α 1, α 2) = 1 2 K 11 α 21 + 1 2 K 22 α 22 + y 1 y 2 K 12 α 1 α 2 + y 1 α 1 \sum i = 3 m y i α i K i 1 + y 2 α 2 \sum i = 3 m y i α i K i 2 - (α 1 + α 2) s . t . α 1 y 1 + α 2 y 2 = - \sum i = 3 m y i α i = ζ 0 \leq α i \leq C i = 1, 2

上式中，

K11=K(xi,xj)。由约束条件

α1y1+α2y2=ζ，可得

α 1 = (ζ - y 2 α 2) y 1

将

α1代入到目标函数

W中，对

α2求导，令其为0，可得

α n e w, u n c u t 2 = α o l d 2 + y 2 ( E 1 - E 2 ) η

其中，我们定义

ζ=αold1y1+αold2y2，

η=K11+K22−2K12，

Ei=∑mj=1αjyjK(xj,xi)+b−yi i=1,2，即训练时输入样本的预测误差。要使

αnew,uncut2满足条件约束，需将其限制在某一区域

[L,H]内。从而得到剪辑之后的

α2值。由于只有两个变量

α1、α2，约束可以用下面的框框表示。

由于约束条件α1y1+α2y2=ζ的存在，使得最优值出现在一条平行于对角线的线段上。当y1≠y2时，线段斜率为正，否则为负。为方便理解，我们以上图为例，假设αold2>αold1，初始可行解为 (αold1,αold2)，令αold2−αold1=K,可以看到，αnew2的取值范围必将在[K,C]之间，考虑所有情况，可以得到如下范围：
如果y1≠y2：

L = m a x (0, α o l d 2 - α o l d 1), H = m i n (C, C + α o l d 2 - α o l d 1)

如果

y1=y2：

L = m a x (0, α o l d 2 + α o l d 1 - C), H = m i n (C, α o l d 2 + α o l d 1)

有了这样的范围，我们可以得到经剪辑之后的解

αnew2:

α n e w 2 = ⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ H α n e w, u n c u t 2 > H α n e w, u n c u t 2 L \leq α n e w, u n c u t 2 \leq H L α n e w, u n c u t 2 < L

另外，根据：

α n e w 1 y 1 + α n e w 2 y 2 = α o l d 1 y 1 + α o l d 2 y 2

可得

α n e w 1 = α o l d 1 + y 1 y 2 (α o l d 2 - α n e w 2)

那么到这里，我们完成了SMO最优化问题的解

α1、α2的求取。但是变量这么多，如何来决定每次选择哪两个变量来优化呢？SMO将其分为两个循环：通过外循环找到第一个变量，再通过内循环寻找第二个变量。

外循环

外层循环在训练样本中选择违反KKT条件最严重的样本点，将改点对应的变量作为第1个变量。具体来看，就是检验某个样例(xi,yi)是否满足KKT条件，即

α i = 0 \Leftrightarrow y i g (x i) \geq 1 0 < α i < C \Leftrightarrow y i g (x i) = 1 α i = C \Leftrightarrow y i g (x i) \leq = 1

其中，

g(xi)=∑mj=1αjyjK(xi,xj)+b。上述KKT条件的推导可参看线性支持向量机中KKT条件的讨论。具体操作时，可以选一个

ε为精度。一般直接遍历所有满足条件

0<αi<C的样本点，即处于间隔边界上的支持向量，检验其是否满足KKT条件，如果都满足，那么遍历整个训练集，检验其是否满足KKT条件。

内循环

假设已经在外循环中找到了第一个变量α1，那么对第二个变量的要求是希望其能有足够大的变化。由αnew,uncut2=αold2+y2(E1−E2)η可知，αnew2是依赖于|E1−E2|的，那么我们直接选择对应|E1−E2|最大的α2即可。如果α2变化极小，那么放弃α1，重新选择。

更新b，Ei

在完成一组两变量的优化后，需要重新计算b，Ei值。如果0<αnew1<C，那么该样例处于间隔边界上。由KKT条件，有y1g(x1)=1，左右同乘y1，有

\sum α i y i K i 1 + b = y 1

进而有

b n e w 1 = y 1 - \sum i = 3 m α i y i K i 1 - α n e w 1 y 1 K 11 - α n e w 2 y 2 K 21

根据

Ei的定义

E 1 = \sum i = 3 m α i y i K i 1 + α o l d 1 y 1 K 11 + α o l d 2 y 2 K 21 + b o l d - y 1

代入至上上式中，可得

b n e w 1 = - E 1 - y 1 K 11 (α n e w 1 - α o l d 1) - y 2 K 21 (α n e w 2 - α o l d 2) + b o l d

同理，若

0<αnew2<C，则有

b n e w 2 = - E 2 - y 1 K 12 (α n e w 1 - α o l d 1) - y 2 K 22 (α n e w 2 - α o l d 2) + b o l d

如果

α1、α2均满足

0<αnewi<C，那么两点均在间隔边界上，求得的截距

b相同，即

bnew1=bnew2，否则取二者中点作为

bnew。随后，更新

Ei：

E i = \sum y j α j K (x i, x j) + b (n e w) - y i

至此，完成了所有

α值的优化。若在精度范围内满足前述KKT停机条件，则算法结束。

迈门演示·Python源码

注：所用数据为UCI葡萄酒数据集，采用Class1、Class2来验证程序。数据集见葡萄酒数据集下载链接。

# !/usr/bin/env python3# coding=utf-8"""Support Vector Machine,SVMAuthor  :Chai ZhengBlog    :http://blog.csdn.net/chai_zheng/Github  :https://github.com/Chai-Zheng/Machine-LearningEmail   :zchaizju@gmail.comDate    :2017.10.3"""import timeimport numpy as npimport matplotlib.pyplot as pltfrom sklearn import preprocessingclass SVM():    def __init__(self,dataset,labels,C,toler,kernelOption):        self.train_x = dataset        self.train_y = labels        self.C = C        self.toler = toler        self.numSamples = dataset.shape[0]        self.alphas = np.zeros((self.numSamples,1))        self.b = 0        self.errorCache = np.zeros((self.numSamples,2))        self.kernelOpt = kernelOption        self.kernelMat = calcKernelMatrix(self.train_x,self.kernelOpt)def calcKernelValue(matrix_x,sample_x,kernelOption):    kernelType = kernelOption[0]    numSamples = matrix_x.shape[0]    kernelValue = np.zeros((1,numSamples))    if kernelType == 'linear':            kernelValue = np.dot(matrix_x,sample_x.T)    elif kernelType == 'rbf':            sigma = kernelOption[1]            if sigma == 0:                sigma =1            for i in range(numSamples):                diff = matrix_x[i,:] - sample_x                kernelValue[0,i] = np.exp((np.dot(diff,diff.T)/(-2.0*sigma**2)))    else:        raise NameError('Not support kernel type! You should use linear or rbf!')    return kernelValuedef calcKernelMatrix(train_x,kernelOption):    numSamples = train_x.shape[0]    kernelMatrix = np.zeros((numSamples,numSamples))    for i in range(numSamples):        kernelMatrix[i,:] = calcKernelValue(train_x,train_x[i,:],kernelOption)    return kernelMatrixdef calcError(svm,alpha_k):     #计算第k个样本的误差，k∈[1，m]    output_k = float(np.dot((svm.alphas*svm.train_y).T,svm.kernelMat[:,alpha_k])+svm.b)    error_k = output_k-float(svm.train_y[alpha_k])    return error_kdef updateError(svm,alpha_k):   #更新误差    error = calcError(svm,alpha_k)    svm.errorCache[alpha_k] = [1,error]def innerLoop(svm,alpha_1,error_1,train_x): #内循环，根据alpha1确定alpha2    svm.errorCache[alpha_1] =[1,error_1]    candidateAlphaList = np.nonzero(svm.errorCache[:,0])[0]    maxStep = 0    alpha_2 = 0    error_2 = 0    numSample = train_x.shape[0]    if len(candidateAlphaList)>1:        #找出|E2-E1|最大的alpha2        for alpha_k in candidateAlphaList:            if alpha_k == alpha_1:                continue            error_k = calcError(svm,alpha_k)            if abs(error_1-error_k)>maxStep:                maxStep = abs(error_1-error_k)                alpha_2 = alpha_k                error_2 = error_k    else:   #第一次进入，随机选择alpha2        while alpha_2 == alpha_1:   #alpha_2不能等于alpha_1            alpha_2 = np.random.randint(svm.numSamples)        error_2 = calcError(svm,alpha_2)        #采用下述方式来初始化alpha_2位置，可稳定结果。采用上述方法会稍微不稳定，但这是正常的。        # if alpha_1 == numSample:        #     alpha_2 = numSample - 1        # else:        #     alpha_2 = alpha_1 + 1        # error_2 = calcError(svm,alpha_2)    return alpha_2,error_2def outsideLoop(svm,alpha_1,train_x):    error_1 = calcError(svm,alpha_1)    #检查alpha_1是否违背KKT条件    if ((svm.alphas[alpha_1]<svm.C) and (svm.alphas[alpha_1]>0) and ((svm.train_y[alpha_1]*error_1>svm.toler) or (svm.train_y[alpha_1]*error_1< -svm.toler)))\                or((svm.train_y[alpha_1]*error_1<-svm.toler) and (svm.alphas[alpha_1]<svm.C))\            or((svm.train_y[alpha_1]*error_1>svm.toler) and (svm.alphas[alpha_1]>0)):        #固定alpha1，求alpha2        alpha_2,error_2 = innerLoop(svm,alpha_1,error_1,train_x)        alpha_1_old = svm.alphas[alpha_1].copy()    #拷贝，分配新的内存        alpha_2_old = svm.alphas[alpha_2].copy()        #alpha2的取值范围，其中L=<alpha2<=H，参见李航《统计学习方法》P126        if svm.train_y[alpha_1] != svm.train_y[alpha_2]:            L = max(0,alpha_2_old-alpha_1_old)            H = min(svm.C,svm.C+alpha_2_old-alpha_1_old)        else:            L = max(0,alpha_2_old+alpha_1_old-svm.C)            H = min(svm.C,alpha_2_old+alpha_1_old)        eta = svm.kernelMat[alpha_1,alpha_1]+svm.kernelMat[alpha_2,alpha_2]-2.0*svm.kernelMat[alpha_1,alpha_2]        svm.alphas[alpha_2] += svm.train_y[alpha_2]*(error_1-error_2)/eta  #计算alpha2_new        #对alpha2进行剪辑        if svm.alphas[alpha_2]>H:            svm.alphas[alpha_2] = H        elif svm.alphas[alpha_2]<L:            svm.alphas[alpha_2] = L        #如果alpha2无变化，返回，重选alpha1        if abs(svm.alphas[alpha_2]-alpha_2_old)<0.00001:            updateError(svm,alpha_2)            return 0        #更新alpha1        svm.alphas[alpha_1] += svm.train_y[alpha_1]*svm.train_y[alpha_2]*(alpha_2_old-svm.alphas[alpha_2])        #更新b        b1 = svm.b-error_1-svm.train_y[alpha_1]*svm.kernelMat[alpha_1,alpha_1]*(svm.alphas[alpha_1]-alpha_1_old)\                          -svm.train_y[alpha_2]*svm.kernelMat[alpha_2,alpha_1]*(svm.alphas[alpha_2]-alpha_2_old)        b2 = svm.b-error_2-svm.train_y[alpha_1]*svm.kernelMat[alpha_1,alpha_2]*(svm.alphas[alpha_1]-alpha_1_old)\                          -svm.train_y[alpha_2]*svm.kernelMat[alpha_2,alpha_2]*(svm.alphas[alpha_2]-alpha_2_old)        #alpha2经剪辑，始终在(0,C)内。若1也满足，那么b1=b2；若1不满足，取均值        if (svm.alphas[alpha_1]>0) and (svm.alphas[alpha_1]<svm.C):            svm.b = b1        else:            svm.b = (b1+b2)/2.0        updateError(svm,alpha_1)        updateError(svm,alpha_2)        return 1    else:        return 0def SVMtrain(train_x,train_y,C,toler,maxIter,kernelOption=('rbf',1.0)):    startTime = time.time()    svm = SVM(train_x,train_y,C,toler,kernelOption)    alphaPairsChanged = 1    iterCount = 0    # 迭代终止条件：    #     1.到达最大迭代次数    #     2.迭代完所有样本后alpha不再变化，也就是所有alpha均满足KTT条件    while(iterCount<maxIter) and (alphaPairsChanged>0):        alphaPairsChanged = 0   #标记在该次循环中，alpha有无被优化        SupportAlphaList = np.nonzero((svm.alphas>0)*(svm.alphas<svm.C))[0] #支撑向量序号列表        for i in SupportAlphaList:          #遍历支持向量            alphaPairsChanged += outsideLoop(svm,i,train_x)        for i in range(svm.numSamples):     #遍历所有样本            alphaPairsChanged += outsideLoop(svm,i,train_x)        iterCount += 1    print('---Training Completed.Took %f s.Using %s kernel.'%((time.time()-startTime),kernelOption[0]))    return svmdef SVMtest(svm,test_x,test_y):    numTestSamples = test_x.shape[0]    matchCount = 0    for i in range(numTestSamples):        kernelValue = calcKernelValue(svm.train_x,test_x[i,:],svm.kernelOpt)        predict = np.dot(kernelValue,svm.train_y*svm.alphas)+svm.b        if np.sign(predict) == np.sign(test_y[i]):            matchCount += 1    accuracy = float(matchCount/numTestSamples)    return accuracydef SVMvisible(svm):    #仅针对二变量样本可视化，即被注释掉的训练数据，非葡萄酒数据    w = np.zeros((2,1))    for i in range(svm.numSamples):        if svm.train_y[i] == -1:            plt.plot(svm.train_x[i,0],svm.train_x[i,1],'or')        elif svm.train_y[i] ==1:            plt.plot(svm.train_x[i,0],svm.train_x[i,1],'ob')        w += (svm.alphas[i]*svm.train_y[i]*svm.train_x[i,:].T).reshape(2,1)    supportVectorIndex = np.nonzero(svm.alphas>0)[0]    for i in supportVectorIndex:        plt.plot(svm.train_x[i,0],svm.train_x[i,1],'oy')    min_x = min(svm.train_x[:,0])    max_x = max(svm.train_x[:,0])    min_y = float((-svm.b-w[0]*min_x)/w[1])    max_y = float((-svm.b-w[0]*max_x)/w[1])    plt.plot([min_x,max_x],[min_y,max_y],'-g')    plt.show()if __name__ =='__main__':    print('Step 1.Loading data...')    #构建10个训练样本，6个测试样本，线性可分,若采用被注释的数据，可将本程序的最后一行取消注释，从而可视化结果    # train_data = np.array([[2.95,6.63,1],[2.53,7.79,1],[3.57,5.65,1],[2.16,6.22,-1],[3.27,3.52,-1],[3,7,1],[3,8,1],[3,2,-1],[2,9,1],[2,4,-1]])    # test_data = np.array([[3.16,5.47,1],[2.58,4.46,-1],[2,2,-1],[3,4,-1],[5,100,1],[6,1000,1]])    #    # train_x = train_data[:,0:2]    # train_y = train_data[:,2].reshape(10,1)    # test_x = test_data[:,0:2]    # test_y = test_data[:,2].reshape(6,1)    #数据集下载http://download.csdn.net/download/chai_zheng/10009314    train_data = np.loadtxt("Wine_Train.txt",delimiter=',') #载入葡萄酒数据集    test_data = np.loadtxt("Wine_Test.txt",delimiter=',')    train_x = train_data[:,1:14]    scaler = preprocessing.StandardScaler().fit(train_x)    train_x = scaler.transform(train_x) #数据标准化    train_y = train_data[:,0].reshape(65,1)    for i in range(len(train_y)):        if train_y[i] == 1:     #修改标签为±1            train_y[i] = -1        if train_y[i] == 2:            train_y[i] = 1    test_x = test_data[:,1:14]    test_x = scaler.transform(test_x)   #数据标准化    test_y = test_data[:,0].reshape(65,1)    for i in range(len(test_y)):        if test_y[i] == 1:            test_y[i] = -1        if test_y[i] == 2:            test_y[i] = 1    print('---Loading completed.')    print('Step 2.Training...')    C = 0.6    toler = 0.001    maxIter = 100    svmClassifier = SVMtrain(train_x,train_y,C,toler,maxIter,kernelOption=('rbf',2))    print('Step 3.Testing...')    accuracy = SVMtest(svmClassifier,test_x,test_y)    print('---Testing completed.Accuracy: %.3f%%'%(accuracy*100))    # SVMvisible(svmClassifier)

运行效果如下：使用rbf核，标准差为2：

使用线性核：

参考文献

[1]李航. 统计学习方法[M]. 清华大学出版社, 2012.
[2]周志华. 机器学习[M]. 清华大学出版社, 2016.
[3]http://blog.csdn.net/xianlingmao/article/details/7919597
[4]http://blog.csdn.net/zouxy09/article/details/17291543

阅读全文

2 0