python机器学习案例系列教程——支持向量机SVM、核函数

全栈工程师开发手册（作者：栾鹏）

python数据挖掘系列教程

线性函数、线性回归、线性分类

参考：http://blog.csdn.net/luanpeng825485697/article/details/78933084

线性回归：

线性回归就是计算回归系数，通过回归系数线性组合属性预测结果值。

f(x)=wTx+bf(x)=w^Tx+bf(x)=wTx+b

二分类的线性分类：

只不过在线性回归的基础上条件了大于和小于0的判断。

y={1,f(x)≥0−1,f(x)<0y= \begin{cases} 1, & \text {f(x)$\geq$0} \\ -1, & \text{f(x)<0} \end{cases} y={1,−1,f(x)≥0f(x)<0

使用线性分类时存在线性可分和线性不可分。

比如：下图为线性可分的

当然，也存在着许多线性不可分的情况，例如下图所示

即使是线性可分的，分类边界也不一定是确定的。比如下面的几个分类边界都是可以实现分类的，哪个更好呢。

这就需要用到支持向量机SVM算法了。

支持向量机/SVM

SVM（Support Vector Machines）是分类算法中应用广泛、效果不错的一类。由简至繁SVM可分类为三类：线性可分（linear SVM in linearly separable case）的线性SVM、线性不可分的线性SVM、非线性（nonlinear）SVM。

线性可分

对于二类分类问题，训练集T=(x1,y1),(x2,y2),⋯,(xN,yN)T={(x_1,y_1),(x_2,y_2),⋯,(x_N,y_N)}T=(x1,y1),(x2,y2),⋯,(xN,yN)，其中xix_ixi为输入样本向量，类别yi∈−1,1y_i∈{-1,1}yi∈−1,1，线性SVM通过学习得到分离超平面（hyperplane）:

w⋅x+b=0w⋅x+b=0w⋅x+b=0

以及相应的分类决策函数：

f(x)=sign(w⋅x+b)f(x)=sign(w⋅x+b)f(x)=sign(w⋅x+b)

其中signsignsign为符号函数，将正数映射为1，负数映射为-1。

有如下图所示的分离超平面，哪一个超平面的分类效果更好呢？

直观上，超平面B1的分类效果更好一些。将距离分离超平面最近的两个不同类别的样本点称为支持向量（support vector）的，构成了两条平行于分离超平面的长带，二者之间的距离称之为margin。显然，margin更大，则分类正确的确信度更高（与超平面的距离表示分类的确信度，距离越远则分类正确的确信度越高）。通过计算得到（计算过程参考《统计学习方法》115页）

margin=2∥w∥margin=\frac{2}{∥w∥}margin=∥w∥2

从上图中可观察到：margin以外的样本点对于确定分离超平面没有贡献，换句话说，SVM是有很重要的训练样本（支持向量）所确定的。至此，SVM分类问题可描述为在全部分类正确的情况下，最大化2∥w∥\frac{2}{∥w∥}∥w∥2（等价于最小化0.5∥w∥20.5∥w∥^20.5∥w∥2）；线性分类的约束最优化问题转化为如下的凸优化问题。

min⁡w,b12∣∣w∣∣2\min_{w,b}\frac{1}{2}||w||^2w,bmin21∣∣w∣∣2 s.t.yi(w⋅xi+b)−1≥0s.t.\quad \quad y_i(w\cdot x_i+b)-1 \geq 0s.t.yi(w⋅xi+b)−1≥0

在线性可分情况下，训练数据集的样本点中与分离超平面距离最近的样本点的实例称为支持向量。

支持向量满足yi(w⋅xi+b)−1=0y_i(w\cdot x_i+b)-1 =0yi(w⋅xi+b)−1=0，即正样本点满足w⋅xi+b=1w\cdot x_i+b=1w⋅xi+b=1，负样本点满足w⋅xi+b=−1w\cdot x_i+b=-1w⋅xi+b=−1

线性支持向量机

线性可分是理想情形，大多数情况下，由于噪声或特异点等各种原因，训练样本是线性不可分的。因此，需要更一般化的学习算法。

这时我们就可以通过引入所谓的松弛变量(slack variable)，来允许有些数据点可以处于超平面的错误的一侧。这样我们的优化目标就能保持仍然不变，但是此时我们的约束条件有所改变。

凸优化问题

（a）无约束优化问题，可以写为：

minf(x)min\quad f(x)minf(x)

（b）有等式约束的优化问题，可以写为：

minf(x)s.t.hi(x)=0,i=0,1,2...nmin\quad f(x) \\ s.t. \quad h_i(x)=0 ,i=0,1,2...nminf(x)s.t.hi(x)=0,i=0,1,2...n

（c）有不等式约束的优化问题，可以写为：

minf(x)s.t.gi(x)≤0,i=0,1,2...nhj(x)=0,j=0,1,2...mmin\quad f(x) \\ s.t. \quad g_i(x) \leq 0 ,i=0,1,2...n\\ h_j(x)=0 ,j=0,1,2...mminf(x)s.t.gi(x)≤0,i=0,1,2...nhj(x)=0,j=0,1,2...m

对于第(a)类的优化问题，尝尝使用的方法就是费马大定理(Fermat)，即使用求取函数f(x)的导数，然后令其为零，可以求得候选最优值，再在这些候选值中验证；如果是凸函数，可以保证是最优解。这也就是我们高中经常使用的求函数的极值的方法。

对于第(b)类的优化问题，常常使用的方法就是拉格朗日乘子法（Lagrange Multiplier) ，即把等式约束h_i(x)用一个系数与f(x)写为一个式子，称为拉格朗日函数，而系数称为拉格朗日乘子。通过拉格朗日函数对各个变量求导，令其为零，可以求得候选值集合，然后验证求得最优值。

对于第©类的优化问题，常常使用的方法就是KKT条件。同样地，我们把所有的等式、不等式约束与f(x)写为一个式子，也叫拉格朗日函数，系数也称拉格朗日乘子，通过一些条件，可以求出最优值的必要条件，这个条件称为KKT条件。

必要条件和充要条件如果不理解，可以看下面这句话：

事件M的必要条件就是M可以推出的结论
事件M的充分条件就是可以推出M的前提

了解到这些，现在让我们再看一下我们的最优化问题：

min⁡w12∣∣w∣∣2s.t.yi(wTxi+b)≥1,i=0,1,2...n\min_w \quad \frac{1}{2} ||w||^2 \\ s.t. \quad y_i(w^Tx_i+b) \geq 1 ,i=0,1,2...n wmin21∣∣w∣∣2s.t.yi(wTxi+b)≥1,i=0,1,2...n

现在，我们的这个对优化问题属于哪一类？很显然，它属于第©类问题。因为，在学习求解最优化问题之前，我们还要学习两个东西：拉格朗日函数和KKT条件。

拉格朗日函数

首先，我们先要从宏观的视野上了解一下拉格朗日对偶问题出现的原因和背景。

我们知道我们要求解的是最小化问题，所以一个直观的想法是如果我能够构造一个函数，使得该函数在可行解区域内与原目标函数完全一致，而在可行解区域外的数值非常大，甚至是无穷大，那么这个没有约束条件的新目标函数的优化问题就与原来有约束条件的原始目标函数的优化问题是等价的问题。这就是使用拉格朗日方程的目的，它将约束条件放到目标函数中，从而将有约束优化问题转换为无约束优化问题。

随后，人们又发现，使用拉格朗日获得的函数，使用求导的方法求解依然困难。进而，需要对问题再进行一次转换，即使用一个数学技巧：拉格朗日对偶。
所以，显而易见的是，我们在拉格朗日优化我们的问题这个道路上，需要进行下面二个步骤：

将有约束的原始目标函数转换为无约束的新构造的拉格朗日目标函数
使用拉格朗日对偶性，将不易求解的优化问题转化为易求解的优化

下面是拉格朗日对偶问题的公式，不想推导的记住就行了。

想看推导的参考：https://www.zhihu.com/question/58584814

我们考虑优化问题如下，记作问题（P）。

z∗=min⁡xf(x)s.t.gi(x)≤0,∀i=1,…,m,x∈Xz^* = \min_x f(x) \\ {s.t. } g_i(x)\leq 0 ,~\forall~ i=1,\ldots,m,x\in Xz∗=xminf(x)s.t.gi(x)≤0, ∀ i=1,…,m,x∈X

问题（P）的拉格朗日对偶问题（D）写作

v∗=max⁡α≥0min⁡x∈Xf(x)+αTg(x)⏟L(x,α)v^* = \max_{α\geq 0} \min_{x\in X} \underbrace{f(x)+α^T g(x)}_{L(x,α)}v∗=α≥0maxx∈XminL(x,α)f(x)+αTg(x)

其中的函数L(x,α)就是我们熟知的拉格朗日函数。

下面，先将有约束的原始目标函数转换为无约束的新构造的拉格朗日目标函数

公式变形如下：

L(w,b,α)=12∣∣w∣∣2−∑i=1nαi(yi(wTxi+b)−1)L(w,b,α) = \frac{1}{2}||w||^2-\sum_{i=1}^nα_i(y_i(w^Tx_i+b)-1)L(w,b,α)=21∣∣w∣∣2−i=1∑nαi(yi(wTxi+b)−1)

其中αiα_iαi是拉格朗日乘子，αiα_iαi大于等于0，是我们构造新目标函数时引入的系数变量。

其中αiαiαi是拉格朗日乘子，αiαiαi大于等于0，是我们构造新目标函数时引入的系数变量(我们自己设置)。现在我们令：

θ(w)=max⁡αi≥0L(w,b,α)\theta(w) = \max_{α_i≥0} L(w,b,α) θ(w)=αi≥0maxL(w,b,α)

当样本点不满足约束条件时，即在可行解区域外：

yi(wTxi+b)<1y_i(w^Tx_i+b)<1yi(wTxi+b)<1

此时，我们将αiαiαi设置为正无穷，此时θ(w)θ(w)θ(w)显然也是正无穷。

当样本点满足约束条件时，即在可行解区域内：

yi(wTxi+b)≥1y_i(w^Tx_i+b)≥1yi(wTxi+b)≥1

此时，我们设在yi(wTxi+b)>1y_i(w^Tx_i+b)>1yi(wTxi+b)>1时的αiαiαi为0，在yi(wTxi+b)=1y_i(w^Tx_i+b)=1yi(wTxi+b)=1时(也就是支持向量点)的αiαiαi为大于0的数，那么显然θ(w)θ(w)θ(w)为原目标函数本身。我们将上述两种情况结合一下，就得到了新的目标函数：

θ(w)={12∣∣w∣∣2,x在可行区域+∞x在非可行区域θ(w) = \begin{cases} \frac{1}{2}||w||^2, & \text {x在可行区域} \\ +\infty& \text{x在非可行区域} \end{cases} θ(w)={21∣∣w∣∣2,+∞x在可行区域x在非可行区域

此时，再看我们的初衷，就是为了建立一个在可行解区域内与原目标函数相同，在可行解区域外函数值趋近于无穷大的新函数，现在我们做到了。

现在，我们的问题变成了求新目标函数的最小值，即：

p∗=min⁡w,bθ(w)=min⁡w,bmax⁡αi≥0L(w,b,α)p^* = \min_{w,b}θ(w) = \min_{w,b} \max_{αi≥0}L(w,b,α)p∗=w,bminθ(w)=w,bminαi≥0maxL(w,b,α)

这里用p*表示这个问题的最优值，且和最初的问题是等价的。

接下来，我们进行第二步：将不易求解的优化问题转化为易求解的优化

我们看一下我们的新目标函数，先求最大值，再求最小值。这样的话，我们首先就要面对带有需要求解的参数w和b的方程，而αiαiαi又是不等式约束，这个求解过程不好做。所以，我们需要使用拉格朗日函数对偶性，将最小和最大的位置交换一下，这样就变成了：

d∗=max⁡α≥0min⁡w,bL(w,b,α)d^* = \max_{α\geq 0} \min_{w,b} L(w,b,α)d∗=α≥0maxw,bminL(w,b,α)

交换以后的新问题是原始问题的对偶问题，这个新问题的最优值用d∗d^*d∗来表示。而且d∗≤p∗d^*≤p^*d∗≤p∗。我们关心的是d=pd=pd=p的时候，这才是我们要的解。需要什么条件才能让d=pd=pd=p呢？

首先必须满足这个优化问题是凸优化问题。
其次，需要满足KKT条件。

凸优化问题的定义是：求取最小值的目标函数为凸函数的一类优化问题。目标函数是凸函数我们已经知道，这个优化问题又是求最小值。所以我们的最优化问题就是凸优化问题。

接下里，就是探讨是否满足KKT条件了。

KKT条件

我们已经使用拉格朗日函数对我们的目标函数进行了处理，生成了一个新的目标函数。通过一些条件，可以求出最优值的必要条件，这个条件就是接下来要说的KKT条件。一个最优化模型能够表示成下列标准形式：

min⁡f(x)s.t.hj(x)=0,j=1,2,3...p,s.t.gk(x)≤0,k=1,2,....,q,x∈X⊂Rn\min f(x) \\ s.t. \,\,\, h_j(x)=0, \,\,\,\,\,\, j=1,2,3...p, \\ s.t. \,\,\,\, g_k(x)\leq 0 ,k=1,2,....,q, \\ x \in X \subset R^n minf(x)s.t.hj(x)=0,j=1,2,3...p,s.t.gk(x)≤0,k=1,2,....,q,x∈X⊂Rn

优化模型的拉格朗日函数为

L(X,λ,μ)=f(X)+∑j=1pλjhj(X)+∑k=1qμkgk(X)L(X,λ,μ)=f(X)+\sum_{j=1}^p λ_jh_j(X)+\sum_{k=1}^qμ_kg_k(X)L(X,λ,μ)=f(X)+j=1∑pλjhj(X)+k=1∑qμkgk(X)

其中f(x)f(x)f(x)是原目标函数，hj(x)h_j(x)hj(x)是第jjj个等式约束条件，λjλ_jλj是对应的约束系数，gkg_kgk是不等式约束，uku_kuk是对应的约束系数。

KKT条件的全称是Karush-Kuhn-Tucker条件，KKT条件是说最优值条件必须满足以下条件：

其中X∗X^*X∗表示样本点。这些求解条件就是KKT条件。(1)是对拉格朗日函数取极值时候带来的一个必要条件，(2)是拉格朗日系数约束（同等式情况），(3)是不等式约束情况，(4)是互补松弛条件，(5)、(6)是原约束条件。

对于一般的任意问题而言，KKT条件是使一组解成为最优解的必要条件，当原问题是凸问题的时候，KKT条件也是充分条件。

对于我们的优化问题：

min⁡w12∣∣w∣∣2s.t.yi(wTxi+b)≥1,i=0,1,2...n\min_w \quad \frac{1}{2} ||w||^2 \\ s.t. \quad y_i(w^Tx_i+b) \geq 1 ,i=0,1,2...n wmin21∣∣w∣∣2s.t.yi(wTxi+b)≥1,i=0,1,2...n

显然，条件二已经满足了。另外两个条件为啥也满足呢？

这里原谅我省略一系列证明步骤，感兴趣的可以移步这里：点我查看

这里已经给出了很好的解释。现在，凸优化问题和KKT都满足了，问题转换成了对偶问题。而求解这个对偶学习问题，可以分为三个步骤：首先要让L(w,b,α)关于w和b最小化，然后求对α的极大，最后利用SMO算法求解对偶问题中的拉格朗日乘子。现在，我们继续推导。

对偶问题求解

第一步：根据上述推导已知：

L(w,b,α)=12∣∣w∣∣2−∑i=1nαi(yi(wTxi+b)−1)L(w,b,α) = \frac{1}{2}||w||^2-\sum_{i=1}^nα_i(y_i(w^Tx_i+b)-1)L(w,b,α)=21∣∣w∣∣2−i=1∑nαi(yi(wTxi+b)−1)

首先固定ααα，要让L(w,b,α)L(w,b,α)L(w,b,α)关于www和bbb最小化，我们分别对www和bbb偏导数，令其等于0，即：

dLdw=0⟹w=∑i=1nαiyixi\frac{d L}{d w} = 0 \implies w = \sum_{i=1}^n α_iy_ix_idwdL=0⟹w=i=1∑nαiyixi
dLdb=0⟹∑i=1nαiyi=0\frac{d L}{d b} = 0 \implies \sum_{i=1}^n α_iy_i=0dbdL=0⟹i=1∑nαiyi=0

将上述结果带回L(w,b,α)L(w,b,α)L(w,b,α)得到

L(w,b,α)=∑i=1nαi−0.5∗∑i,j=1nαiαjyiyjxiTxjL(w,b,α) = \sum_{i=1}^n α_i-0.5*\sum_{i,j=1}^n α_iα_jy_iy_jx_i^Tx_jL(w,b,α)=i=1∑nαi−0.5∗i,j=1∑nαiαjyiyjxiTxj

从上面的最后一个式子，我们可以看出，此时的L(w,b,α)L(w,b,α)L(w,b,α)函数只含有一个变量，即αiα_iαi。

第二步：

现在内侧的最小值求解完成，我们求解外侧的最大值，从上面的式子得到

max⁡α∑i=1nαi−12∑i,j=1nαiαjyiyjxiTxjs.t.αi≥0,i=1,2...n,∑i=1nαiyi=0\max_α\sum_{i=1}^n α_i-\frac{1}{2}\sum_{i,j=1}^n α_iα_jy_iy_jx_i^Tx_j \\ s.t. \,\,\,\, α_i\geq0,i=1,2...n, \\ \sum_{i=1}^n α_iy_i=0αmaxi=1∑nαi−21i,j=1∑nαiαjyiyjxiTxjs.t.αi≥0,i=1,2...n,i=1∑nαiyi=0

现在我们的优化问题变成了如上的形式。对于这个问题，我们有更高效的优化算法，即序列最小优化（SMO）算法。我们通过这个优化算法能得到α，再根据α，我们就可以求解出w和b，进而求得我们最初的目的：找到超平面，即"决策平面"。

而上式可以等价于
min⁡α12∑i,j=1nαiαjyiyjxiTxj−∑i=1nαis.t.αi≥0,i=1,2...n,∑i=1nαiyi=0\min_α \frac{1}{2}\sum_{i,j=1}^n α_iα_jy_iy_jx_i^Tx_j -\sum_{i=1}^n α_i\\ s.t. \,\,\,\, α_i\geq0,i=1,2...n, \\ \sum_{i=1}^n α_iy_i=0αmin21i,j=1∑nαiαjyiyjxiTxj−i=1∑nαis.t.αi≥0,i=1,2...n,i=1∑nαiyi=0

总结一句话：我们为啥使出吃奶的劲儿进行推导？因为我们要将最初的原始问题，转换到可以使用SMO算法求解的问题，这是一种最流行的求解方法。

SMO算法

SMO算法的目标是求出一系列αiα_iαi和bbb，一旦求出了这些αiα_iαi，就很容易计算出权重向量w并得到分隔超平面。

SMO算法的工作原理是：每次循环中选择两个αiα_iαi进行优化处理。一旦找到了一对合适的αiα_iαi，那么就增大其中一个同时减小另一个。这里所谓的"合适"就是指两个αiα_iαi必须符合以下两个条件，条件之一就是两个αiα_iαi必须要在间隔边界之外，而且第二个条件则是这两个αiα_iαi还没有进行过区间化处理或者不在边界上。

SMO算法的解法

先来定义特征到结果的输出函数为：

u=wTx+bu=w^Tx+bu=wTx+b
接着，我们回忆一下原始优化问题，如下：

min⁡w12∣∣w∣∣2s.t.yi(wTxi+b)≥1,i=0,1,2...n\min_w \quad \frac{1}{2} ||w||^2 \\ s.t. \quad y_i(w^Tx_i+b) \geq 1 ,i=0,1,2...n wmin21∣∣w∣∣2s.t.yi(wTxi+b)≥1,i=0,1,2...n

求导得：

w=∑i=1nαiyixiw = \sum_{i=1}^n α_iy_ix_iw=i=1∑nαiyixi
将上述公式带入输出函数中：

u=∑i=1nαiyixiTx+bu = \sum_{i=1}^n α_iy_ix_i^Tx+bu=i=1∑nαiyixiTx+b

与此同时，拉格朗日对偶后得到最终的目标化函数：

min⁡α12∑i,j=1nαiαjyiyjxiTxj−∑i=1nαis.t.αi≥0,i=1,2...n,∑i=1nαiyi=0\min_α \frac{1}{2}\sum_{i,j=1}^n α_iα_jy_iy_jx_i^Tx_j -\sum_{i=1}^n α_i\\ s.t. \,\,\,\, α_i\geq0,i=1,2...n, \\ \sum_{i=1}^n α_iy_i=0αmin21i,j=1∑nαiαjyiyjxiTxj−i=1∑nαis.t.αi≥0,i=1,2...n,i=1∑nαiyi=0

线性支持下的SMO

实际上，对于上述目标函数，是存在一个假设的，即数据100%线性可分。但是，目前为止，我们知道几乎所有数据都不那么"干净"。这时我们就可以通过引入所谓的松弛变量(slack variable)，来允许有些数据点可以处于超平面的错误的一侧。这样我们的优化目标就能保持仍然不变，但是此时我们的约束条件有所改变：

s.t.C≥αi≥0,i=1,2...n,∑i=1nαiyi=0s.t.\quad C\geqα_i\geq0,i=1,2...n, \\ \sum_{i=1}^n α_iy_i=0s.t.C≥αi≥0,i=1,2...n,i=1∑nαiyi=0

根据KKT条件可以得出其中αiαiαi取值的意义为：

αi=0⟺yiui≥1α_i = 0 \iff y_iu_i\geq1αi=0⟺yiui≥1
0<αi<C⟺yiui=10<α_i <C \iff y_iu_i=10<αi<C⟺yiui=1
αi=C⟺yiui≤1α_i = C \iff y_iu_i\leq1αi=C⟺yiui≤1

对于第1种情况，表明αiαiαi是正常分类，在边界内部；
对于第2种情况，表明αiαiαi是支持向量，在边界上；
对于第3种情况，表明αiαiαi是在两条边界之间。

而最优解需要满足KKT条件，即上述3个条件都得满足，以下几种情况出现将会不满足

yiui≤1αi<Cy_iu_i\leq1 \quad α_i <C yiui≤1αi<C
yiui≥1αi>0y_iu_i\geq1 \quad α_i > 0yiui≥1αi>0
yiui=1αi=0或者αi=Cy_iu_i=1 \quad α_i = 0 或者 α_i =C yiui=1αi=0或者αi=C

也就是说，如果存在不能满足KKT条件的αiαiαi，那么需要更新这些αiαiαi，这是第一个约束条件。此外，更新的同时还要受到第二个约束条件的限制，即：

∑i=1nαiyi=0\sum_{i=1}^n α_iy_i=0i=1∑nαiyi=0

因为这个条件，我们同时更新两个ααα值，因为只有成对更新，才能保证更新之后的值仍然满足和为0的约束，假设我们选择的两个乘子为α1α_1α1和α2α_2α2：
α1newy1+α2newy2=α1oldy1+α2oldy2=ζα_1^{new}y_1+α_2^{new}y_2 = α_1^{old}y_1+α_2^{old}y_2 = \zetaα1newy1+α2newy2=α1oldy1+α2oldy2=ζ

其中， $ \zeta为常数。因为两个因子不好同时求解，所以可以先求第二个乘子为常数。因为两个因子不好同时求解，所以可以先求第二个乘子为常数。因为两个因子不好同时求解，所以可以先求第二个乘子α_2的解的解的解(α_2^{{new})$，得到$α_2$的解$(α_2}{new})之后，再用之后，再用之后，再用α_2的解的解的解(α_2^{{new})$表示$α_1$的解$(α_1}{new})。为了求解。为了求解。为了求解(α_2^{new})$ ，得先确定(α2new)(α_2^{new})(α2new)的取值范围。假设它的上下边界分别为H和L，那么有：

L≤α2new≤HL\leqα_2^{new}\leq HL≤α2new≤H

接下来，综合下面两个条件：

C≥αi≥0,i=1,2...n,α1newy1+α2newy2=α1oldy1+α2oldy2=ζC\geqα_i\geq0,i=1,2...n, \\ α_1^{new}y_1+α_2^{new}y_2 = α_1^{old}y_1+α_2^{old}y_2 = \zetaC≥αi≥0,i=1,2...n,α1newy1+α2newy2=α1oldy1+α2oldy2=ζ

当y1不等于y2时，即一个为正1，一个为负1的时候，可以得到：

α1old−α2old=ζα_1^{old}-α_2^{old}= \zetaα1old−α2old=ζ

所以有：

L=max(0,ζ),H=min(C,C−ζ)L=max(0, \zeta),H=min(C,C-\zeta)L=max(0,ζ),H=min(C,C−ζ)

此时，取值范围如下图所示：

当y1等于y2时，即两个都为正1或者都为负1，可以得到：

α1old+α2old=ζα_1^{old}+α_2^{old}= \zetaα1old+α2old=ζ

所以有：

L=max(0,ζ−C),H=min(C,ζ)L=max(0, \zeta-C),H=min(C,\zeta)L=max(0,ζ−C),H=min(C,ζ)

此时，取值范围如下图所示：

如此，根据y1和y2异号或同号，可以得出α2newα_2^{new}α2new的上下界分别为：
L=max(0,α2old−α1old),H=min(C,C+α2old−α1old)ify1≠y2L=max(0, α_2^{old}-α_1^{old}),H=min(C,C+α_2^{old}-α_1^{old}) \qquad if \quad y_1 \neq y_2L=max(0,α2old−α1old),H=min(C,C+α2old−α1old)ify1=y2
L=max(0,α1old+α2old−C),H=min(C,α1old+α2old)ify1=y2L=max(0, α_1^{old}+α_2^{old}-C),H=min(C,α_1^{old}+α_2^{old}) \qquad if \quad y_1=y_2L=max(0,α1old+α2old−C),H=min(C,α1old+α2old)ify1=y2

这个界限就是编程的时候需要用到的。已经确定了边界，接下来，就是推导迭代式，用于更新 ααα值。

我们已经知道，更新ααα的边界，接下来就是讨论如何更新ααα值。我们依然假设选择的两个乘子为α1α_1α1和α2α_2α2。(推导略)

α2new,clipped={H,α2new>Hα2new,L≤α2new≤HL,α2new<Lα_2^{new,clipped}= \begin{cases} H, & α_2^{new}>H \\ α_2^{new}, & L\leq α_2^{new}\leq H \\ L, & α_2^{new}<L \end{cases} α2new,clipped=⎩⎪⎨⎪⎧H,α2new,L,α2new>HL≤α2new≤Hα2new<L

α1new=α1old+y1y2(α2old−α2new,clipped)α_1^{new}= α_1^{old}+y_1y_2(α_2^{old}-α_2^{new,clipped}) α1new=α1old+y1y2(α2old−α2new,clipped)

这样，我们就知道了怎样计算α1α_1α1和α2α_2α2了，也就是如何对选择的α进行更新。

我们要根据$ α$ 的取值范围，去更正b的值，使间隔最大化。当α1newα_1^{new}α1new在0和C之间的时候，根据KKT条件可知，这个点是支持向量上的点。因此，满足下列公式：

y1(wTx1+b)=1y_1(w^Tx_1+b)= 1 y1(wTx1+b)=1

公式两边同时乘以y1y_1y1得(y1×y1=1y_1\times y_1=1y1×y1=1)：

∑i=1nαiyixix1+b=y1\sum_{i=1}^nα_iy_ix_ix_1+b=y_1i=1∑nαiyixix1+b=y1

因为我们是根据α1α_1α1和α2α_2α2的值去更新b，所以单独提出i=1和i=2的时候，整理可得：

b1new=y1−∑i=3nαiyixiTx1−α1newy1x1Tx1−α2newy2x2Tx1b_1^{new} = y_1-\sum_{i=3}^nα_iy_ix_i^Tx_1-α_1^{new}y_1x_1^Tx_1-α_2^{new}y_2x_2^Tx_1b1new=y1−i=3∑nαiyixiTx1−α1newy1x1Tx1−α2newy2x2Tx1

其中前两项为：

y1−∑i=3nαiyixiTx1=−E1+α1oldy1x1Tx1+α2oldy2x2Tx1+boldy_1-\sum_{i=3}^nα_iy_ix_i^Tx_1=-E_1+α_1^{old}y_1x_1^Tx_1+α_2^{old}y_2x_2^Tx_1+b^{old}y1−i=3∑nαiyixiTx1=−E1+α1oldy1x1Tx1+α2oldy2x2Tx1+bold

将上述两个公式，整理得：

b1new=bold−E1−y1(α1new−α1old)x1Tx1−y2(α2new−α2old)x2Tx1b_1^{new} =b^{old}-E_1-y_1(α_1^{new}-α_1^{old})x_1^Tx_1-y_2(α_2^{new}-α_2^{old})x_2^Tx_1b1new=bold−E1−y1(α1new−α1old)x1Tx1−y2(α2new−α2old)x2Tx1

同理可得b2newb_2^{new}b2new为：

b2new=bold−E2−y1(α1new−α1old)x1Tx2−y2(α2new−α2old)x2Tx2b_2^{new} =b^{old}-E_2-y_1(α_1^{new}-α_1^{old})x_1^Tx_2-y_2(α_2^{new}-α_2^{old})x_2^Tx_2b2new=bold−E2−y1(α1new−α1old)x1Tx2−y2(α2new−α2old)x2Tx2

当我们更新了α1α_1α1和α2α_2α2之后，需要重新计算阈值b，因为b关系到了我们f(x)f(x)f(x)的计算，也就关系到了误差Ei的计算。

b={b1,0<α1new<Cb2,0<α2new<C≤H(b1+b2)/2,otherwiseb= \begin{cases} b_1, & 0<α_1^{new}<C \\ b_2, & 0<α_2^{new}<C \leq H \\ (b_1+b_2)/2, & otherwise \end{cases} b=⎩⎪⎨⎪⎧b1,b2,(b1+b2)/2,0<α1new<C0<α2new<C≤Hotherwise

现在，让我们梳理下SMO算法的步骤：

步骤1：计算误差：

Ei=f(xi)−yi=∑j=1nαjyjxiTxj+b−yiE_i=f(x_i)-y_i=\sum_{j=1}^nα_jy_jx_i^Tx_j+b-y_iEi=f(xi)−yi=j=1∑nαjyjxiTxj+b−yi

步骤2：计算上下界L和H：

L=max(0,αjold−αiold),H=min(C,C+αjold−αiold)ifyi≠yjL=max(0, α_j^{old}-α_i^{old}),H=min(C,C+α_j^{old}-α_i^{old}) \qquad if \quad y_i \neq y_jL=max(0,αjold−αiold),H=min(C,C+αjold−αiold)ifyi=yj
L=max(0,αjold+αiold−C),H=min(C,αjold+αiold)ifyj=yiL=max(0, α_j^{old}+α_i^{old}-C),H=min(C,α_j^{old}+α_i^{old}) \qquad if \quad y_j=y_iL=max(0,αjold+αiold−C),H=min(C,αjold+αiold)ifyj=yi

步骤3：计算ηηη：

η=xiTxi+xjTxj−2xiTxjη = x_i^Tx_i+x_j^Tx_j-2x_i^Tx_jη=xiTxi+xjTxj−2xiTxj

步骤4：更新αjα_jαj：

αjnew=αjold+yj(Ei−Ej)ηα_j^{new}=α_j^{old}+\frac{y_j(E_i-E_j)}{η }αjnew=αjold+ηyj(Ei−Ej)

步骤5：根据取值范围修剪αjα_jαj

αjnew,clipped={H,αjnew>Hαjnew,L≤αjnew≤HL,αjnew<Lα_j^{new,clipped}= \begin{cases} H, & α_j^{new}>H \\ α_j^{new}, & L\leq α_j^{new}\leq H \\ L, & α_j^{new}<L \end{cases} αjnew,clipped=⎩⎪⎨⎪⎧H,αjnew,L,αjnew>HL≤αjnew≤Hαjnew<L

步骤6：更新αiα_iαi：

αinew=αiold+yiyj(αjold−αjnew,clipped)α_i^{new}= α_i^{old}+y_iy_j(α_j^{old}-α_j^{new,clipped}) αinew=αiold+yiyj(αjold−αjnew,clipped)

步骤7：更新b1和b2：

b1new=bold−Ei−yi(αinew−αiold)xiTxi−yj(αjnew−αjold)xjTxib_1^{new} =b^{old}-E_i-y_i(α_i^{new}-α_i^{old})x_i^Tx_i-y_j(α_j^{new}-α_j^{old})x_j^Tx_ib1new=bold−Ei−yi(αinew−αiold)xiTxi−yj(αjnew−αjold)xjTxi

b2new=bold−Ej−yi(αinew−αiold)xiTxj−yj(αjnew−αjold)xjTxjb_2^{new} =b^{old}-E_j-y_i(α_i^{new}-α_i^{old})x_i^Tx_j-y_j(α_j^{new}-α_j^{old})x_j^Tx_jb2new=bold−Ej−yi(αinew−αiold)xiTxj−yj(αjnew−αjold)xjTxj

步骤8：根据b1和b2更新b：

python实现

样本集下载：https://github.com/626626cdllp/data-mining/blob/master/SVM/testSet.txt

# -*- coding:UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import random# 简化版smo# 函数说明:读取数据
def loadDataSet(fileName):alldata = np.loadtxt(fileName)dataMat = alldata[:,0:2]   #添加数据labelMat = alldata[:,2]   #.astype(int).reshape(-1,1)  #添加标签return dataMat,labelMat"""
函数说明:随机选择alphaParameters:i - alpha_i的索引值m - alpha参数个数
Returns:j - alpha_j的索引值"""
def selectJrand(i, m):j = i                                 #选择一个不等于i的jwhile (j == i):j = int(random.uniform(0, m))return j"""
函数说明:修剪alphaParameters:aj - alpha_j值H - alpha上限L - alpha下限
Returns:aj - alpah值"""
def clipAlpha(aj,H,L):if aj > H:aj = Hif L > aj:aj = Lreturn aj# 函数说明:数据可视化
def showDataSet(dataMat, labelMat):place_plus = np.where(labelMat==1)[0]   # 正样本的位置place_minus = np.where(labelMat==-1)[0]  # 负样本的位置data_plus = dataMat[place_plus]    #正样本data_minus = dataMat[place_minus]  #负样本plt.scatter(np.transpose(data_plus)[0], np.transpose(data_plus)[1])   #正样本散点图plt.scatter(np.transpose(data_minus)[0], np.transpose(data_minus)[1]) #负样本散点图plt.show()"""
函数说明:简化版SMO算法Parameters:dataMatIn - 数据矩阵classLabels - 数据标签C - 松弛变量toler - 容错率maxIter - 最大迭代次数
"""
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):#转换为numpy的mat存储dataMatrix = np.mat(dataMatIn)labelMat = np.mat(classLabels).transpose()#初始化b参数，统计dataMatrix的维度b = 0; m,n = np.shape(dataMatrix)#初始化alpha参数，设为0alphas = np.mat(np.zeros((m,1)))#初始化迭代次数iter_num = 0#最多迭代matIter次while (iter_num < maxIter):alphaPairsChanged = 0for i in range(m):#步骤1：计算误差EifXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + bEi = fXi - float(labelMat[i])#优化alpha，设定一定的容错率。if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):#随机选择另一个与alpha_i成对优化的alpha_jj = selectJrand(i,m)#步骤1：计算误差EjfXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + bEj = fXj - float(labelMat[j])#保存更新前的aplpha值，使用深拷贝alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();#步骤2：计算上下界L和Hif (labelMat[i] != labelMat[j]):L = max(0, alphas[j] - alphas[i])H = min(C, C + alphas[j] - alphas[i])else:L = max(0, alphas[j] + alphas[i] - C)H = min(C, alphas[j] + alphas[i])if L==H: print("L==H"); continue#步骤3：计算etaeta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].Tif eta >= 0: print("eta>=0"); continue#步骤4：更新alpha_jalphas[j] -= labelMat[j]*(Ei - Ej)/eta#步骤5：修剪alpha_jalphas[j] = clipAlpha(alphas[j],H,L)if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小"); continue#步骤6：更新alpha_ialphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#步骤7：更新b_1和b_2b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].Tb2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T#步骤8：根据b_1和b_2更新bif (0 < alphas[i]) and (C > alphas[i]): b = b1elif (0 < alphas[j]) and (C > alphas[j]): b = b2else: b = (b1 + b2)/2.0#统计优化次数alphaPairsChanged += 1#打印统计信息print("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num,i,alphaPairsChanged))#更新迭代次数if (alphaPairsChanged == 0): iter_num += 1else: iter_num = 0print("迭代次数: %d" % iter_num)return b,alphas"""
函数说明:分类结果可视化Parameters:dataMat - 数据矩阵w - 直线法向量b - 直线解决
"""
def showClassifer(dataMat, w, b):# 绘制样本点place_plus = np.where(labelMat==1)[0]   # 正样本的位置place_minus = np.where(labelMat==-1)[0]  # 负样本的位置data_plus = dataMat[place_plus]    #正样本data_minus = dataMat[place_minus]  #负样本plt.scatter(np.transpose(data_plus)[0], np.transpose(data_plus)[1],s=30, alpha=0.7)   #正样本散点图plt.scatter(np.transpose(data_minus)[0], np.transpose(data_minus)[1], s=30, alpha=0.7) #负样本散点图#绘制直线x1 = max(dataMat[:,0])  # 第一个属性的最大值x2 = min(dataMat[:,0])  # 第一个属性的最小值a1, a2 = wb = float(b)a1 = float(a1[0])a2 = float(a2[0])y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2plt.plot([x1, x2], [y1, y2])#找出支持向量点for i, alpha in enumerate(alphas):if abs(alpha) > 0:x, y = dataMat[i]plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')plt.show()"""
函数说明:计算wParameters:dataMat - 数据矩阵labelMat - 数据标签alphas - alphas值
"""
def get_w(dataMat, labelMat, alphas):alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)return w.tolist()if __name__ == '__main__':dataMat, labelMat = loadDataSet('testSet.txt')showDataSet(dataMat,labelMat)b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)w = get_w(dataMat, labelMat, alphas)showClassifer(dataMat, w, b)

3、非线性向量机

1 核技巧

我们已经了解到，SVM如何处理线性可分的情况，而对于非线性的情况，SVM的处理方式就是选择一个核函数。简而言之：在线性不可分的情况下，SVM通过某种事先选择的非线性映射（核函数）将输入变量映到一个高维特征空间，将其变成在高维空间线性可分，在这个高维空间中构造最优分类超平面。

线性可分的情况下，可知最终的超平面方程为：
f(x)=∑i=1nαiyixiTx+bf(x)=\sum_{i=1}^nα_iy_ix_i^Tx+bf(x)=i=1∑nαiyixiTx+b

将上述公式用内积来表示：
f(x)=∑i=1nαiyi<xi,x>+bf(x)=\sum_{i=1}^nα_iy_i<x_i,x>+bf(x)=i=1∑nαiyi<xi,x>+b

对于线性不可分，我们使用一个非线性映射，将数据映射到特征空间，在特征空间中使用线性学习器，分类函数变形如下：
f(x)=∑i=1nαiyi<ϕ(xi),ϕ(x)>+bf(x)=\sum_{i=1}^nα_iy_i<ϕ(x_i),ϕ(x)>+bf(x)=i=1∑nαiyi<ϕ(xi),ϕ(x)>+b

其中ϕϕϕ从输入空间(X)到某个特征空间(F)的映射，这意味着建立非线性学习器分为两步：

首先使用一个非线性映射将数据变换到一个特征空间F；
然后在特征空间使用线性学习器分类。

如果有一种方法可以在特征空间中直接计算内积<ϕ(xi),ϕ(x)><ϕ(x_i),ϕ(x)><ϕ(xi),ϕ(x)>，就像在原始输入点的函数中一样，就有可能将两个步骤融合到一起建立一个分线性的学习器，这样直接计算的方法称为核函数方法。

这里直接给出一个定义：核是一个函数k，对所有x,z∈Xx,z∈Xx,z∈X，满足k(x,z)=<ϕ(xi),ϕ(x)>k(x,z)=<ϕ(x_i),ϕ(x)>k(x,z)=<ϕ(xi),ϕ(x)>，这里ϕ(⋅)ϕ(·)ϕ(⋅)是从原始输入空间X到内积空间F的映射。

简而言之：如果不是用核技术，就会先计算线性映ϕ(x1)和ϕ(x2)ϕ(x_1)和ϕ(x_2)ϕ(x1)和ϕ(x2)，然后计算这它们的内积，使用了核技术之后，先把ϕ(x1)ϕ(x_1)ϕ(x1)和ϕ(x2)ϕ(x_2)ϕ(x2)的一般表达式<ϕ(x1),ϕ(x2)>=k(<ϕ(x1),ϕ(x2)>)<ϕ(x_1),ϕ(x_2)>=k(<ϕ(x_1),ϕ(x_2) >)<ϕ(x1),ϕ(x2)>=k(<ϕ(x1),ϕ(x2)>)计算出来，这里的<⋅，⋅><·，·><⋅，⋅>表示内积，k(⋅，⋅)k(·，·)k(⋅，⋅)就是对应的核函数，这个表达式往往非常简单，所以计算非常方便。

这种将内积替换成核函数的方式被称为核技巧(kernel trick)。

核函数的数学要求

核函数有严格的数学要求，所以设计一个核函数是很困难的。K(x,z)是正定核的充要条件是：K（x,z）对应的Gram矩阵实半正定矩阵。

Gram矩阵：矩阵对应点的内积。KTK, KKT

半正定矩阵：设A是实对称矩阵。如果对任意的实非零列矩阵X有XTAX≥0，就称A为半正定矩阵。

当检验一个K是否为正定核函数，要对任意有限输入集{xi…}验证K对应的Gram矩阵实是否为半正定矩阵。

LIBSVM中提供的核函数

线性核函数：
K(x,z)=x⋅zK(x,z)=x\cdot zK(x,z)=x⋅z
线性核，主要用于线性可分的情况，我们可以看到特征空间到输入空间的维度是一样的，其参数少速度快，对于线性可分数据，其分类效果很理想，因此我们通常首先尝试用线性核函数来做分类，看看效果如何，如果不行再换别的

多项式核函数：
K(x,z)=(x⋅z+1)pK(x,z)=(x\cdot z+1)^pK(x,z)=(x⋅z+1)p
对应的支持向量机为p次多项式分类器。
多项式核函数可以实现将低维的输入空间映射到高纬的特征空间，但是多项式核函数的参数多，当多项式的阶数比较高的时候，核矩阵的元素值将趋于无穷大或者无穷小，计算复杂度会大到无法计算。

RBF核函数（高斯核函数）：
K(x,z)=exp(−∣∣x−z∣∣22σ2)K(x,z)=exp(-\frac{||x-z||^2}{2\sigma^2})K(x,z)=exp(−2σ2∣∣x−z∣∣2)
对应的支持向量机为高斯径向基函数分类器。
高斯径向基函数是一种局部性强的核函数，其可以将一个样本映射到一个更高维的空间内，该核函数是应用最广的一个，无论大样本还是小样本都有比较好的性能，而且其相对于多项式核函数参数要少，因此大多数情况下在不知道用什么核函数的时候，优先使用高斯核函数。

sigmoid核函数：
κ(x,z)=tanh(η<x,z>+θ)\kappa(x,z) = tanh(\eta<x,z> + \theta)κ(x,z)=tanh(η<x,z>+θ)
采用sigmoid核函数，支持向量机实现的就是一种多层神经网络。

SVM的核函数如何选取

最常用的是Linear核与RBF核。需要注意的是需要对数据归一化处理。

1、Linear核：主要用于线性可分的情形。参数少，速度快，对于一般数据，分类效果已经很理想了。
2、RBF核：主要用于线性不可分的情形。参数多，分类结果非常依赖于参数。

有很多人是通过训练数据的交叉验证来寻找合适的参数，不过这个过程比较耗时。我个人的体会是：使用libsvm，默认参数，RBF核比Linear核效果稍差。通过进行大量参数的尝试，一般能找到比linear核更好的效果。如果特征的提取的好，包含的信息量足够大，很多问题都是线性可分的。当然，如果有足够的时间去寻找RBF核参数，应该能达到更好的效果。

这些函数中应用最广的应该就是RBF核了，无论是小样本还是大样本，高维还是低维等情况，RBF核函数均适用，它相比其他的函数有一下优点：1）RBF核函数可以将一个样本映射到一个更高维的空间，而且线性核函数是RBF的一个特例，也就是说如果考虑使用RBF，那么就没有必要考虑线性核函数了。2）与多项式核函数相比，RBF需要确定的参数要少，核函数参数的多少直接影响函数的复杂程度。另外，当多项式的阶数比较高时，核矩阵的元素值将趋于无穷大或无穷小，而RBF则在上，会减少数值的计算困难。3）对于某些参数，RBF和sigmoid具有相似的性能。

因此，在选用核函数的时候，如果我们对我们的数据有一定的先验知识，就利用先验来选择符合数据分布的核函数；如果不知道的话，通常使用交叉验证的方法，来试用不同的核函数，误差最下的即为效果最好的核函数，或者也可以将多个核函数结合起来，形成混合核函数。在吴恩达的课上，也曾经给出过一系列的选择核函数的方法：

下面是吴恩达的见解：
1、如果Feature的数量很大，跟样本数量差不多，这时候选用LR或者是Linear Kernel的SVM
2、如果Feature的数量比较小，样本数量一般，不算大也不算小，选用SVM+Gaussian Kernel
3、如果Feature的数量比较小，而样本数量很多，需要手工添加一些feature变成第一种情况

总结

SVM的优缺点

优点

可用于线性/非线性分类，也可以用于回归，泛化错误率低，也就是说具有良好的学习能力，且学到的结果具有很好的推广性。
可以解决小样本情况下的机器学习问题，可以解决高维问题，可以避免神经网络结构选择和局部极小点问题。
SVM是最好的现成的分类器，现成是指不加修改可直接使用。并且能够得到较低的错误率，SVM可以对训练集之外的数据点做很好的分类决策。

缺点

对参数调节和和函数的选择敏感。