跟我一起学scikit-learn19：支持向量机算法

支持向量机（SVM，Support Vector Machine）算法是一种常见的分类算法，在工业界和学术界都有广泛的应用。特别是针对数据集较小的情况下，往往其分类效果比神经网络好。

1.SVM算法原理

SVM的原理就是使用分隔超平面来划分数据集，并使得支持向量（数据集中离分隔超平面最近的点）到该分隔超平面的距离最大。其最大特点是能构造出最大间距的决策边界，从而提高分类算法的鲁棒性。

1.大间距分类算法

假设要对一个数据集进行分类，如下图所示，可以构造一个分隔线把圆形的点和方形的点分开。这个分隔线称为分隔超平面（Separating Hyperplane）。

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
class1 = np.array([[1, 1], [1, 3], [2, 1], [1, 2], [2, 2]])
class2 = np.array([[4, 4], [5, 5], [5, 4], [5, 3], [4, 5], [6, 4]])
plt.figure(figsize=(8, 6), dpi=144)

plt.title('Decision Boundary')

plt.xlim(0, 8)
plt.ylim(0, 6)
ax = plt.gca()                                  # gca 代表当前坐标轴，即 'get current axis'
ax.spines['right'].set_color('none')            # 隐藏坐标轴
ax.spines['top'].set_color('none')

plt.scatter(class1[:, 0], class1[:, 1], marker='o')
plt.scatter(class2[:, 0], class2[:, 1], marker='s')
plt.plot([1, 5], [5, 1], '-r')
plt.arrow(4, 4, -1, -1, shape='full', color='r')
plt.plot([3, 3], [0.5, 6], '--b')
plt.arrow(4, 4, -1, 0, shape='full', color='b', linestyle='--')
plt.annotate(r'margin 1',
             xy=(3.5, 4), xycoords='data',
             xytext=(3.1, 4.5), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'margin 2',
             xy=(3.5, 3.5), xycoords='data',
             xytext=(4, 3.5), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'support vector',
             xy=(4, 4), xycoords='data',
             xytext=(5, 4.5), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'support vector',
             xy=(2, 2), xycoords='data',
             xytext=(0.5, 1.5), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
figure8_1.png

从上图可以明显地看出，实线的分隔线比虚线的分隔线更好，因为使用实线的分隔线进行分类时，离分隔线最近的点到分隔线的距离最大，即margin2>margin1。这段距离的2倍称为间距（margin）。那些离分隔超平面最近的点，称为支持向量（Support Vector）。为了达到最好的分类效果，SVM算法的原理就是要找到一个分隔超平面，它能把数据集正确地分类，并且间距最大。

首先，我们来看怎么计算间距。在二维空间里，可以使用方程w_1x_1+w_2x_2+b=0来表示分隔超平面。针对高维度空间，可以写成一般的向量化形式，即w^Tx+b=0。我们画出与分隔超平面平行的两条直线，分别穿过两个类别的支持向量（离分隔超平面距离最近的点）。这两条直线的方程分别为w^Tx+b=-1和w^Tx+b=1。如下图所示。

plt.figure(figsize=(8, 6), dpi=144)

plt.title('Support Vector Machine')

plt.xlim(0, 8)
plt.ylim(0, 6)
ax = plt.gca()                                  # gca 代表当前坐标轴，即 'get current axis'
ax.spines['right'].set_color('none')            # 隐藏坐标轴
ax.spines['top'].set_color('none')

plt.scatter(class1[:, 0], class1[:, 1], marker='o')
plt.scatter(class2[:, 0], class2[:, 1], marker='s')
plt.plot([1, 5], [5, 1], '-r')
plt.plot([0, 4], [4, 0], '--b', [2, 6], [6, 2], '--b')
plt.arrow(4, 4, -1, -1, shape='full', color='b')
plt.annotate(r'$w^T x + b = 0$',
             xy=(5, 1), xycoords='data',
             xytext=(6, 1), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'$w^T x + b = 1$',
             xy=(6, 2), xycoords='data',
             xytext=(7, 2), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'$w^T x + b = -1$',
             xy=(3.5, 0.5), xycoords='data',
             xytext=(4.5, 0.2), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'd',
             xy=(3.5, 3.5), xycoords='data',
             xytext=(2, 4.5), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'A',
             xy=(4, 4), xycoords='data',
             xytext=(5, 4.5), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
figure8_2.png

根据点到直线的距离公式，可以容易地算出支持向量A到分隔超平面的距离为：
d=\frac{|w^TA+b|}{||w||}

由于点A在直线w^Tx+b=1上，因此w^TA+b=1，代入即可算出支持向量A到分隔超平面的距离为d=\frac{1}{||w||}。为了使得距离最大，我们只需要找到合适的参数w和b，使得\frac{1}{||w||}最大即可。||w||是向量w的L2范数，其计算公式为：
||w||=\sqrt{\sum_{i=1}^{n}w_i^2}
由此可得，求\frac{1}{||w||}的最大值，等价于求||w||^2的最小值：
||w||^2=\sum_{i=1}^{n}w_i^2

其中n为向量w的维度。除了间距最大外，我们选出来的分隔超平面还要能正确地把数据集分类。问题来了，怎样在数学上表达出“正确地把数据集分类”这个描述呢？

根据上图可以看出，针对方形的点，必定满足w^Tx+b\geq1的约束条件。类别是离散的值，我们使用-1来表示圆形的类别，用1来表示方形的类别，即y \in \left\{-1,1\right\}。针对数据集中的所有样本x^{(i)}，y^{(i)}，只要它们都满足以下的约束条件，则由w和b定义的分隔超平面即可正确地把数据集分类：
y^{(i)}(w^Tx^{(i)}+b)\geq1

等等，怎么得出这个数学表达式的呢？
其技巧在于使用1和-1来定义类别标签。针对y^{(i)}=1的情况，由于其满足w^Tx^{(i)}+b\geq1的约束，两边都乘以y^{(i)}后，大于等于号保持不变。针对y^{(i)}=-1的情况，由于其满足w^Tx^{(i)}+b\leq-1的约束，两边都乘以y^{(i)}后，负负得正，小于等于号变成了大于等于号。这样，我们就可以使用一个公式来表达针对两个不同类别的约束函数了。

在逻辑回归算法里，使用0和1作为类别标签，而这里我们使用-1和1作为类别标签。其目的都是为了让数学表达尽量简洁。

总结：求解SVM算法，就是在满足约束条件y^{(i)}(w^Tx^{(i)}+b)\geq1的前提下，求解||w||^2的最小值。

2.松弛系数

针对线性不可分的数据集，上面介绍的方法就失灵了。因为无法找到最大间距的分隔超平面，如下图左所示。

from sklearn.datasets import make_blobs

plt.figure(figsize=(13, 6), dpi=144)

# sub plot 1
plt.subplot(1, 2, 1)

X, y = make_blobs(n_samples=100, 
                  n_features=2, 
                  centers=[(1, 1), (2, 2)], 
                  random_state=4, 
                  shuffle=False,
                  cluster_std=0.4)

plt.title('Non-linear Separatable')

plt.xlim(0, 3)
plt.ylim(0, 3)
ax = plt.gca()                                  # gca 代表当前坐标轴，即 'get current axis'
ax.spines['right'].set_color('none')            # 隐藏坐标轴
ax.spines['top'].set_color('none')

plt.scatter(X[y==0][:, 0], X[y==0][:, 1], marker='o')
plt.scatter(X[y==1][:, 0], X[y==1][:, 1], marker='s')
plt.plot([0.5, 2.5], [2.5, 0.5], '-r')

# sub plot 2
plt.subplot(1, 2, 2)

class1 = np.array([[1, 1], [1, 3], [2, 1], [1, 2], [2, 2], [1.5, 1.5], [1.2, 1.7]])
class2 = np.array([[4, 4], [5, 5], [5, 4], [5, 3], [4, 5], [6, 4], [5.5, 3.5], [4.5, 4.5], [2, 1.5]])

plt.title('Slack Variable')

plt.xlim(0, 7)
plt.ylim(0, 7)
ax = plt.gca()                                  # gca 代表当前坐标轴，即 'get current axis'
ax.spines['right'].set_color('none')            # 隐藏坐标轴
ax.spines['top'].set_color('none')

plt.scatter(class1[:, 0], class1[:, 1], marker='o')
plt.scatter(class2[:, 0], class2[:, 1], marker='s')
plt.plot([1, 5], [5, 1], '-r')
plt.plot([0, 4], [4, 0], '--b', [2, 6], [6, 2], '--b')
plt.arrow(2, 1.5, 2.25, 2.25, shape='full', color='b')
plt.annotate(r'violate margin rule.',
             xy=(2, 1.5), xycoords='data',
             xytext=(0.2, 0.5), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'normal sample. $\epsilon = 0$',
             xy=(4, 5), xycoords='data',
             xytext=(4.5, 5.5), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'$\epsilon > 0$',
             xy=(3, 2.5), xycoords='data',
             xytext=(3, 1.5), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
figure8_3.png

解决这个问题的办法是引入一个参数\varepsilon，称为松弛系数。然后把优化的目标函数变为：
argmin||w||^2+R\sum_{i=1}^{m}\varepsilon_i

其中，m为数据集的个数，R为算法参数。其约束条件相应地变为：
y^{(i)}(w^Tx^{(i)}+b)\geq1-\varepsilon_i

怎么理解松弛系数呢？我们可以把\varepsilon_i理解为数据样本x^{(i)}违反最大间距规则的程度，如上图右所示。针对大部分正常的样本，即满足约束条件的样本\varepsilon_i=0。而对部分违反最大间距规则的样本\varepsilon_i>0。而参数R则表示对违反最大间距规则的样本的“惩罚”力度。当R选择一个很大的值时，我们的目标函数对违反最大间距规则的点的惩罚力度将变得很大。当R选择一个较小的值时，针对那些违反最大间距规则的样本的惩罚力度较小，我们的模型就会倾向于允许部分点违反最大间距规则。我们可以把y^{(i)}(w^Tx^{(i)}+b)作为横坐标，把样本由于违反约束条件所付出的代价J_i作为纵坐标，画出二者的关系如下图所示。

plt.figure(figsize=(8, 4), dpi=144)

plt.title('Cost')

plt.xlim(0, 4)
plt.ylim(0, 2)
plt.xlabel('$y^{(i)} (w^T x^{(i)} + b)$')
plt.ylabel('Cost')
ax = plt.gca()                                  # gca 代表当前坐标轴，即 'get current axis'
ax.spines['right'].set_color('none')            # 隐藏坐标轴
ax.spines['top'].set_color('none')

plt.plot([0, 1], [1.5, 0], '-r')
plt.plot([1, 3], [0.015, 0.015], '-r')
plt.annotate(r'$J_i = R \epsilon_i$ for $y^{(i)} (w^T x^{(i)} + b) \geq 1 - \epsilon_i$',
             xy=(0.7, 0.5), xycoords='data',
             xytext=(1, 1), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate(r'$J_i = 0$ for $y^{(i)} (w^T x^{(i)} + b) \geq 1$',
             xy=(1.5, 0), xycoords='data',
             xytext=(1.8, 0.2), fontsize=10,
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
figure8_4.png

从图中可以清楚地看出，针对那些没有违反约束条件y^{(i)}(w^Tx^{(i)}+b)\geq1的样本，其成本为0。而针对那些违反了约束条件y^{(i)}(w^Tx^{(i)}+b)\geq1-\varepsilon_i的样本，其成本与\varepsilon_i成正比。

从这里的描述可知，引入松弛系数类似于逻辑回归算法里的成本函数引入正则项，目的都是为了纠正过拟合问题，让支持向量机对噪声数据有更强的适应性。当出现一些违反最大间距规则的噪声样本时，仍然希望我们的分隔超平面是原来的样子，这就是松弛系数的作用。

2.核函数

什么是核函数？核函数是特征转换函数。

1.最简单的核函数

根据上面的介绍，我们的任务是找出合适的参数w和b，使得由它们决定的分隔超平面间距最大，且能正确地对数据集进行分类。间距最大是我们的优化目标，正确地对数据集进行分类是约束条件。用数学来表达，在满足约束条件y^{(i)}(w^Tx^{(i)}+b)\geq1，即y^{(i)}(w^Tx^{(i)}+b)-1\geq0的前提下，求\frac{1}{2}||w||^2的最小值。

拉格朗日乘子法是解决约束条件下，求函数极值的理想方法。其方法是引入非负系数\alpha作为约束条件的权重：
L=\frac{1}{2}||w||^2-\sum_{i=1}^{m}\alpha_i(y^{(i)}(w^Tx^{(i)}+b)-1)

公式中，针对数据集中的每个样本x^{(i)}，y^{(i)}，都有一个系数\alpha_i与之对应。根据微积分的知识，极值处的偏导数为0。我们先求L对w的偏导数，并令其偏导数为0：
\frac{\partial L}{\partial w}=w-\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)}=0

从而可以得到w和\alpha的关系：
w=\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)}

至此，我们就知道了，为什么要把求\frac{2}{||w||}的最大值转换为求\frac{1}{2}||w||^2的最小值，其目的就是为了使得w的数学表达尽量简洁优美。接着我们继续求L对b的偏导数，并令其偏导数为0：
\frac{\partial L}{\partial b}=\sum_{i=1}^{m}y^{(i)}\alpha_i=0

把w=\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)}和\sum_{i=1}^{m}y^{(i)}\alpha_i=0代入L，通过代数运算可得：
L=\sum_{i=1}^{m}\alpha_i-\frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{m}\alpha_i\alpha_jy^{(i)}y^{(j)}x^{(i)T}x^{(j)}

这个公式看起来很复杂。我们解释一下公式里各个变量的含义。其中，m是数据集的个数，\alpha是拉格朗日乘子法引入的一个系数，针对数据集中的每个样本x^{(i)}，都有对应的\alpha_i。x^{(i)}是数据集中第i个样本的输入，它是一个向量，y^{(i)}是数据集第i个样本的输出标签，其值为y^{(i)} \in \left \{-1,1 \right \}。

怎么求这个公式的最小值，是数学分析要解决的问题。这是一个典型的二次规划问题。目前广泛应用的是一个称为SMO（序列最小优化）的算法。这里不再展开。

最后求解出来的\alpha有个明显的特点，即大部分\alpha_i=0，这个结论背后的原因很直观，因为只有那些支持向量所对应的样本，直接决定了间隙的大小，其他离分隔超平面太远的样本，对间隙大小根本没有影响。

实际上，推导这个公式的目的就是为了引入支持向量机的另外一个核心概念：核函数。

我们注意到L里的x^{(i)T}x^{(j)}部分，其中x^{(i)}是一个特征向量，所以x^{(i)T}x^{(j)}是一个数值，它是两个输入特征向量的内积。另外，我们的预测函数为：
\hat y = w^Tx+b=\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)T}x+b

当\hat y > 0时预测类别为1，当\hat y < 0时预测类别为-1。注意到预测函数里也包含了式子x^{(i)T}x。我们把K(x^{(i)},x^{(j)})=x^{(i)T}x^{(j)}称为核函数。x^{(i)T}x^{(j)}是两个向量的内积，它的物理含义是衡量两个向量的相似性。典型的，当这两个向量相互垂直时，即完全线性无关，此时x^{(i)T}x^{(j)}=0。引入核函数之后，预测函数就变成：
\hat y=\sum_{i=1}^{m}\alpha_iy^{(i)}K(x^{(i)},x)+b

思考：在前面我们把方形类别的约束定义为w^Tx+b\geq1，把圆形类别的约束定义为w^Tx+b\leq-1。而这里的预测函数，又以0为分界点，即针对输入特征向量x，当w^Tx+b＞0时，预测为方形类别，这是为什么呢？

2.相似性函数

为什么我们要引入核函数？假设我们有一个数据集，只有一个输入特征，要对这个数据集进行分类。由于只有一个输入特征，这些训练样本分布在一条直线上，此时我们很难找出一个分隔超平面来分隔这个数据集。如左图所示。

plt.figure(figsize=(13, 6), dpi=144)

class1 = np.array([[1, 1], [1, 2], [1, 3], [2, 1], [2, 2], [3, 2], [4, 1], [5, 1]])
class2 = np.array([[2.2, 4], [1.5, 5], [1.8, 4.6], [2.4, 5], [3.2, 5], [3.7, 4], [4.5, 4.5], [5.4, 3]])

# sub plot 1
plt.subplot(1, 2, 1)

plt.title('Non-linear Separatable in Low Dimension')

plt.xlim(0, 6)
plt.ylim(0, 6)
plt.yticks(())
plt.xlabel('X1')
ax = plt.gca()                                  # gca 代表当前坐标轴，即 'get current axis'
ax.spines['right'].set_color('none')            # 隐藏坐标轴
ax.spines['top'].set_color('none')
ax.spines['left'].set_color('none')

plt.scatter(class1[:, 0], np.zeros(class1[:, 0].shape[0]) + 0.05, marker='o')
plt.scatter(class2[:, 0], np.zeros(class2[:, 0].shape[0]) + 0.05, marker='s')

# sub plot 2
plt.subplot(1, 2, 2)

plt.title('Linear Separatable in High Dimension')

plt.xlim(0, 6)
plt.ylim(0, 6)
plt.xlabel('X1')
plt.ylabel('X2')
ax = plt.gca()                                  # gca 代表当前坐标轴，即 'get current axis'
ax.spines['right'].set_color('none')            # 隐藏坐标轴
ax.spines['top'].set_color('none')

plt.scatter(class1[:, 0], class1[:, 1], marker='o')
plt.scatter(class2[:, 0], class2[:, 1], marker='s')
plt.plot([1, 5], [3.8, 2], '-r')
figure8_5.png

为了解决这个问题，我们可以想个办法，用一定的规则把这些无法进行线性分隔的样本，映射到更高维度的空间里，然后在高维度空间里找出分隔超平面。针对这个例子，把一维空间里的样本映射到二维空间，这样很容易就能找出一个分隔超平面把这些样本分离开。如上面右图所示。

SVM的核函数就是为了实现这种相似性映射。从前面的介绍我们知道，最简单的核函数是K(x^{(i)},x^{(j)})=x^{(i)T}x^{(j)}，它衡量的是两个输入特征向量的相似性。可以通过定义核函数K(x^{(i)},x^{(j)})来重新定义相似性，从而得到想要的映射。例如，在基因测试领域，我们需要根据DNA分子的特征来定义相似性函数，即核函数。在文本处理领域，也可以自己定义核函数来衡量两个词之间的相似性。

怎样把低维度的空间映射到高维度的空间呢？我们是否还记得之前介绍过的一个解决欠拟合的方法，就是使用多项式来增加特征数，这个本质上就是从低维度到高维度。针对上图左中的例子，我们的输入特征是一维的，即只有[x_1]变量，如果我们要变成二维的，一个方法是把输入特征变为[x_1,x_1^2]，此时的输入特征就变成了一个二维向量。定义这种特征映射的函数为\Phi(x)，称之为相似性函数。针对一个输入特征向量x，经过\Phi(x)作用之后，会变成一个新的、更高维度的输入特征向量。这样在原来低维度计算相似性的运算x^{(i)T}x^{(j)}，就可以转换为高维度空间里进行相似性运算\Phi(x^{(i)})^T\Phi(x^{(j)})。

思考：核函数K(x^{(i)},x^{(j)})和相似性函数\Phi(x)有什么关系？
相似性函数是特征映射函数，比如针对二维的特征向量[x_1,x_2]，我们可以定义相似性函数为\Phi(x)=[x_1,x_2,x_1x_2,x_1^2,x_2^2]，经过相似性函数转换后，二维的特征向量就变成了五维的特征向量。而核函数定义为特征向量的内积，经过相似性函数\Phi(x)转换后，核函数即变为两个五维特征向量的内积，即K(x^{(i)},x^{(j)})=\Phi(x^{(i)})^T\Phi(x^{(j)})。

这里我们介绍相似性函数\Phi(x)的目的，是为了帮助理解核函数的生成过程有其背后的思想。在实际计算过程中，我们不会计算相似性函数及其映射值，因为这样做的计算效率很低。例如，我们把二维空间映射到n维空间，如果n非常大，要在n维空间里计算两个向量的内积，需要n*n次运算才可以完成，这个计算成本是非常高的。

3.常用的核函数

核函数一般和应用场景相关，例如在基因测试领域和文本处理领域，它们的核函数可能是不一样的，有专门针对特定应用领域进行核函数开发和建模的科研人员在从事这方面的研究。虽然核函数和应用场景相关，但实际上还是有一些通用的、“万金油”式的核函数。常用的核函数有两种，一种是多项式核函数，顾名思义，是对输入特征向量增加多项式的一种相似性映射函数，其数学表达为：
K(x^{(i)},x^{(j)})=(\gamma x^{(i)T}x^{(j)}+c)^n

其中\gamma为正数，c为非负数。我们介绍过的线性核函数K(x^{(i)},x^{(j)})=x^{(i)T}x^{(j)}是多项式核函数在n=1, \gamma=1, c=0处的一种特例。在二维空间里，K(x^{(i)},x^{(j)})=x^{(i)T}x^{(j)}只能表达直线的分隔超平面，而多项式核函数K(x^{(i)},x^{(j)})=(\gamma x^{(i)T}x^{(j)}+c)^n在n>1时，可以表达更复杂的、非直线的分隔超平面。

另一个常用的核函数是高斯核函数，其数学表达式为：
K(x^{(i)},x^{(j)})=exp\big (-\frac{(x^{(i)}-x^{(j)})^2}{2\sigma^2}\big)

如果我们的输入特征是一维的标量，那么高斯核函数对应的形状就是一个反钟形的曲线，其参数\sigma控制反钟形的宽度，如下所示。

def gaussian_kernel(x, mean, sigma):
    return np.exp(- (x - mean)**2 / (2 * sigma**2))

x = np.linspace(0, 6, 500)
mean = 1
sigma1 = 0.1
sigma2 = 0.3

plt.figure(figsize=(10, 3), dpi=144)

# sub plot 1
plt.subplot(1, 2, 1)
plt.title('Gaussian for $\sigma={0}$'.format(sigma1))

plt.xlim(0, 2)
plt.ylim(0, 1.1)
ax = plt.gca()                                  # gca 代表当前坐标轴，即 'get current axis'
ax.spines['right'].set_color('none')            # 隐藏坐标轴
ax.spines['top'].set_color('none')

plt.plot(x, gaussian_kernel(x, mean, sigma1), 'r-')

# sub plot 2
plt.subplot(1, 2, 2)
plt.title('Gaussian for $\sigma={0}$'.format(sigma2))

plt.xlim(0, 2)
plt.ylim(0, 1.1)
ax = plt.gca()                                  # gca 代表当前坐标轴，即 'get current axis'
ax.spines['right'].set_color('none')            # 隐藏坐标轴
ax.spines['top'].set_color('none')

plt.plot(x, gaussian_kernel(x, mean, sigma2), 'r-')
figure8_6.png

由于K(x^{(i)},x^{(j)})=\Phi(x^{(i)})^T\Phi(x^{(j)})，经过合适的数学变换，可得高斯核函数对应的特征转换函数为：
\Phi(x)=\sum_{i=0}^{\infty}exp(-x^2)\sqrt{\frac{2^i}{i!}}x^i

注意前面无限多项的累加器\sum_{i=0}^{\infty}，其物理意义就是把特征向量转换到无限多维向量空间里，即高斯核函数可以把输入特征向量扩展到无限维空间里。

接下来看一下高斯核函数对应的预测函数：
\hat y=\sum_{i=1}^{m}\alpha_iy^{(i)}K(x^{(i)},x)+b

其中，K(x^{(i)},x)是高斯核函数，而\alpha_i只在支持向量对应的样本处不为0，其他的样本为0。由此得知，预测函数是中心点在支持向量处的高斯函数的线性组合，其线性组合的系数为\alpha_iy^{(i)}。因此，高斯核函数也称为径向基核函数（RBF，Radial Basic Function），即反钟形函数的线性组合。

4.核函数的对比

我们对几个常用的核函数进行对比，看看它们各有哪些优缺点。

1.线性核函数
K(x^{(i)},x^{(j)})=x^{(i)T}x^{(j)}

这是最简单的核函数，它直接计算两个输入特征向量的内积。它的优点是简单，运算效率高，因为不涉及复杂的变换；结果容易解释，因为总是能生成一个最简洁的线性分隔超平面。它的缺点也很明显，即对线性不可分的数据集没有很好的办法。

2.多项式核函数
K(x^{(i)},x^{(j)})=(\gamma x^{(i)T}x^{(j)}+c)^n

多项式核函数通过多项式来作为特征映射函数，它的优点是可以拟合出复杂的分隔超平面。它的缺点是可选的参数太多，有\gamma，c，n这三个参数要选择，在实践过程中，选择一组合适的参数会变得比较困难；另外一个缺点是，多项式的阶数n不宜太高，否则会给模型的求解带来一些困难。典型地，当x^{(i)T}x^{(j)}<1时，经过n次方运算后会接近于0，而当x^{(i)T}x^{(j)}>1时，经过n次方运算后，又会变得非常大，这样核函数就会变得不够稳定。

3.高斯核函数
K(x^{(i)},x^{(j)})=exp\big (-\frac{(x^{(i)}-x^{(j)})^2}{2\sigma^2}\big)
高斯核函数可以把输入特征映射到无限多维，所以它会比线性核函数功能上要强大很多，并且没有多项式核函数的数值计算那么困难，因为它的核函数计算出来的值永远都在[0,1]之间。高斯核函数还有一个优点是参数容易选择，因为它只有一个参数\sigma。它的缺点是不容易解释，因为映射到无限多维向量空间这个事情显得太不直观；计算速度比较慢；容易造成过拟合，原因是映射到无限维向量空间，这是个非常复杂的模型，它会试图去拟合所有的样本，从而造成过拟合。

在实践中怎么选择核函数呢？更进一步，逻辑回归算法也可以用来解决分类问题，到底是用逻辑回归算法还是用SVM算法呢？假设n是特征个数，m是训练数据集的样本个数，一般情况下可以按照下面的规则来选择算法。

如果n相对m来说比较大，例如n=10000，m=10 ~ 1000，如文本处理问题，这个时候使用逻辑回归算法或者线性核函数的SVM算法都可以；如果n比较小，m中等大小，例如n=1 ~ 1000，m=10 ~ 10000，那么可以使用高斯核函数的SVM算法；如果n比较小，m比较大，例如n=1 ~ 1000，m=50000+，那么一般需要增加特征，此时需要使用多项式核函数或者高斯核函数的SVM算法。

更一般的算法选择原则是，针对数据量很大的问题，我们可以选择复杂一点的模型。虽然复杂模型容易造成过拟合，但由于数据量很大，可以有效地弥补过拟合问题。如果数据量比较小，一般需要选择简单一点的模型，否则很容易造成过拟合，此时要特别注意模型是否欠拟合，如果出现了欠拟合，可以使用增加多项式特征的方法纠正欠拟合问题。

3.scikit-learn里的SVM

在scikit-learn里对SVM的算法实现都在包sklearn.svm下面，其中SVC类是用来进行分类任务的，SVR类是用来进行数值回归任务的。我们可能会有疑问，SVM不是用来进行分类的算法吗，为什么可以用来进行数值回归？实际上这只是数学上的一些扩展而已，在计算机里，可以用离散的数值计算来代替连续的数值回归。我们在K-近邻算法中已经看到过这种扩展实现。

我们以SVC为例。首先需要选择SVM的核函数，由参数kernel来指定，其中值linear表示线性核函数，它只能产生直线形状的分隔超平面；值poly表示多项式核函数，用它可以构建出复杂形状的分隔超平面；值rbf表示径向基核函数，即高斯核函数。

不同的核函数需要指定不同的参数。针对线性核函数，只需要指定参数C，它表示对不符合最大间距规则的样本的惩罚力度，即前面介绍的系数R。针对多项式核函数，除了参数C之外，还需要指定degree，它表示多项式的阶数。针对高斯核函数，除了参数C之外，还需要指定gamma值，这个值对应的是高斯核函数公式中的\frac{1}{2 \sigma^2}的值。

下面先来看一个最简单的例子。我们生成一个有两个特征、包含两种类别的数据，然后用线性核函数的SVM算法进行分类：

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
def plot_hyperplane(clf, X, y, 
                    h=0.02, 
                    draw_sv=True, 
                    title='hyperplan'):
    # create a mesh to plot in
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
                         np.arange(y_min, y_max, h))

    plt.title(title)
    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())
    plt.xticks(())
    plt.yticks(())
    
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    plt.contourf(xx, yy, Z, cmap='hot', alpha=0.5)

    markers = ['o', 's', '^']
    colors = ['b', 'r', 'c']
    labels = np.unique(y)
    for label in labels:
        plt.scatter(X[y==label][:, 0], 
                    X[y==label][:, 1], 
                    c=colors[label], 
                    marker=markers[label])
    if draw_sv:
        sv = clf.support_vectors_
        plt.scatter(sv[:, 0], sv[:, 1], c='y', marker='x')

from sklearn import svm
from sklearn.datasets import make_blobs

X, y = make_blobs(n_samples=100, centers=2, 
                  random_state=0, cluster_std=0.3)
clf = svm.SVC(C=1.0, kernel='linear')
clf.fit(X, y)

plt.figure(figsize=(12, 4), dpi=144)
plot_hyperplane(clf, X, y, h=0.01, 
                title='Maximum Margin Hyperplan')

输出的图形如下所示，其中带有x标记的点即为支持向量，它保存在模型的support_vectors里。


figure8_7.png

此处需要注意的是plot_hyperplane()函数，其主要功能是画出样本点，同时画出分类区间。它的主要原理是使用numpy.meshgrid()生成一个坐标矩阵，最后用contourf()函数为坐标矩阵中不同类别的点填充不同的颜色。其中，contourf()函数是画等高线并填充颜色的函数。

接着来看另外一个例子。我们生成一个有两个特征、包含三种类别的数据集，然后分别构造出4个SVM算法来拟合数据集，分别是线性核函数、三阶多项式核函数、gamma=0.5的高斯核函数，以及gamma=0.1的高斯核函数。最后把这4个SVM算法拟合出来的分隔超平面画出来。

from sklearn import svm
from sklearn.datasets import make_blobs

X, y = make_blobs(n_samples=100, centers=3, 
                  random_state=0, cluster_std=0.8)
clf_linear = svm.SVC(C=1.0, kernel='linear')
clf_poly = svm.SVC(C=1.0, kernel='poly', degree=3)
clf_rbf = svm.SVC(C=1.0, kernel='rbf', gamma=0.5)
clf_rbf2 = svm.SVC(C=1.0, kernel='rbf', gamma=0.1)

plt.figure(figsize=(10, 10), dpi=144)

clfs = [clf_linear, clf_poly, clf_rbf, clf_rbf2]
titles = ['Linear Kernel', 
          'Polynomial Kernel with Degree=3', 
          'Gaussian Kernel with $\gamma=0.5$', 
          'Gaussian Kernel with $\gamma=0.1$']
for clf, i in zip(clfs, range(len(clfs))):
    clf.fit(X, y)
    plt.subplot(2, 2, i+1)
    plot_hyperplane(clf, X, y, title=titles[i])

输出的图形如下所示，其中带有x标记的点即为支持向量。


figure8_8.png

左上角是线性核函数，它只能拟合出直线分隔超平面。右上角是三阶多项式核函数，它能拟合出复杂曲线分隔超平面。左下角是gamma=0.5的高斯核函数，右下角是gamma=0.1的高斯核函数，通过调整参数gamma的值，可以调整分隔超平面的形状。典型地，gamma值太大，越容易造成过拟合，gamma值太小，高斯核函数会退化成线性核函数。我们把代码中的gamma值改为100和0.01后看一下输出图形是什么样的。

思考：左下角gamma=0.5的高斯核函数的图片，带有x标记的点是支持向量。我们之前介绍过，离分隔超平面最近的点是支持向量，为什么很多离分隔超平面很远的点，也是支持向量呢？

原因是高斯核函数把输入特征向量映射到了无限维的向量空间里，在映射后的高维向量空间里，这些点其实是离分隔超平面最近的点。当回到二维向量空间中时，这些点“看起来”就不像是距离分隔超平面最近的点了，但实际上它们就是支持向量。

4.示例：乳腺癌检测

之前我们使用逻辑回归算法进行过乳腺癌检测模型的学习和训练。这里我们再使用支持向量机算法来解决这个问题。首先我们载入数据：

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from sklearn.model_selection import train_test_split
# 载入数据
from sklearn.datasets import load_breast_cancer
cancer = load_breast_cancer()
X = cancer.data
y = cancer.target
print('data shape: {0}; no. positive: {1}; no. negative: {2}'.format(
    X.shape, y[y==1].shape[0], y[y==0].shape[0]))
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

输出如下：

data shape: (569, 30); no. positive: 357; no. negative: 212

可以看出，我们的数据集很小。高斯核函数太复杂，容易造成过拟合，模型效果应该不会太好。我们先用高斯核函数来试一下，看与我们的猜测是否一致。

from sklearn.svm import SVC
clf = SVC(C=1.0, kernel='rbf', gamma=0.1)
clf.fit(X_train, y_train)
train_score = clf.score(X_train, y_train)
test_score = clf.score(X_test, y_test)
print('train score: {0}; test score: {1}'.format(train_score, test_score))

输出如下：

train score: 1.0; test score: 0.5877192982456141

训练数据集分数为1.0，交叉验证数据集分数只有0.52，这是典型的过拟合现象。这里gamma=0.1，这个值相对已经比较小了。我们可以把gamma改的更小如0.0001看看什么结果。

当然，我们完全可以使用前面介绍过的GridSearchCV来自动选择最优参数。我们看看如果使用高斯模型，最优的gamma参数值是多少，其对应的模型交叉验证评分是多少。

from sklearn.model_selection import learning_curve
import numpy as np
def plot_learning_curve(plt, estimator, title, X, y, ylim=None, cv=None, n_jobs=1, train_sizes=np.linspace(.1, 1.0, 5)):
    plt.title(title)
    if ylim is not None:
        plt.ylim(*ylim)
    plt.xlabel("Training examples")
    plt.ylabel("Score")
    train_sizes, train_scores, test_scores = learning_curve(estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes)
    train_scores_mean = np.mean(train_scores, axis=1)
    train_scores_std = np.std(train_scores, axis=1)
    test_scores_mean = np.mean(test_scores, axis=1)
    test_scores_std = np.std(test_scores, axis=1)
    plt.grid()
    plt.fill_between(train_sizes, train_scores_mean - train_scores_std, train_scores_mean + train_scores_std, alpha=0.1, color="r")
    plt.fill_between(train_sizes, test_scores_mean - test_scores_std, test_scores_mean + test_scores_std, alpha=0.1, color="g")
    plt.plot(train_sizes, train_scores_mean, 'o--', color="r", label="Training score")
    plt.plot(train_sizes, test_scores_mean, 'o-', color="g", label="Cross-validation score")
    plt.legend(loc="best")
    return plt

def plot_param_curve(plt, train_sizes, cv_results, xlabel):
    train_scores_mean = cv_results['mean_train_score']
    train_scores_std = cv_results['std_train_score']
    test_scores_mean = cv_results['mean_test_score']
    test_scores_std = cv_results['std_test_score']
    plt.title('parameters turning')
    plt.grid()
    plt.xlabel(xlabel)
    plt.ylabel('score')
    plt.fill_between(train_sizes, train_scores_mean - train_scores_std, train_scores_mean + train_scores_std,alpha=0.1, color="r")
    plt.fill_between(train_sizes, test_scores_mean - test_scores_std,test_scores_mean + test_scores_std, alpha=0.1, color="g")
    plt.plot(train_sizes, train_scores_mean, '.--', color="r", label="Training score")
    plt.plot(train_sizes, test_scores_mean, '.-', color="g", label="Cross-validation score")
    plt.legend(loc="best")
    return plt

from sklearn.model_selection import GridSearchCV
gammas = np.linspace(0, 0.0003, 30)
param_grid = {'gamma': gammas}
clf = GridSearchCV(SVC(), param_grid, cv=5)
clf.fit(X, y)
print("best param: {0}\nbest score: {1}".format(clf.best_params_,clf.best_score_))
plt.figure(figsize=(10, 4), dpi=144)
plot_param_curve(plt, gammas, clf.cv_results_, xlabel='gamma');

输出如下：

best param: {'gamma': 0.00011379310344827585}
best score: 0.9367311072056239
figure8_8_2.png

由此可见，即使在最好的gamma参数下，其平均最优得分也只是0.9367311072056239。我们选择在gamma为0.01时，画出学习曲线，更直观地观察模型拟合情况。

import time
from sklearn.model_selection import ShuffleSplit
cv = ShuffleSplit(n_splits=10, test_size=0.2, random_state=0)
title = 'Learning Curves for Gaussian Kernel'
start = time.clock()
plt.figure(figsize=(10, 4), dpi=144)
plot_learning_curve(plt, SVC(C=1.0, kernel='rbf', gamma=0.01), title, X, y, ylim=(0.5, 1.01), cv=cv)
print('elaspe: {0:.6f}'.format(time.clock()-start))

输出如下：

elaspe: 0.680074

画出来的图形如下所示：


figure8_9.png

这是明显的过拟合现象，交叉验证数据集的评分非常低，且离训练数据集评分非常远。

接下来换一个模型，使用二阶多项式核函数的SVM来拟合模型，看看结果如何。

from sklearn.svm import SVC
clf = SVC(C=1.0, kernel='poly', degree=2)
clf.fit(X_train, y_train)
train_score = clf.score(X_train, y_train)
test_score = clf.score(X_test, y_test)
print('train score: {0}; test score: {1}'.format(train_score, test_score))

输出如下：

train score: 0.9758241758241758; test score: 0.9736842105263158

看起来结果好多了。作为对比，我们画出一阶多项式核函数的SVM和二阶多项式核函数的SVM的学习曲线，观察模型的拟合情况。

import time
from sklearn.model_selection import ShuffleSplit
cv = ShuffleSplit(n_splits=5, test_size=0.2, random_state=0)
title = 'Learning Curves with degree={0}'
degrees = [1, 2]
start = time.clock()
plt.figure(figsize=(12, 4), dpi=144)
for i in range(len(degrees)):
    plt.subplot(1, len(degrees), i + 1)
    plot_learning_curve(plt, SVC(C=1.0, kernel='poly', degree=degrees[i]),
                        title.format(degrees[i]), X, y, ylim=(0.8, 1.01), cv=cv, n_jobs=4)
print('elaspe: {0:.6f}'.format(time.clock()-start))

输出如下：

elaspe: 90.780260

输出的图形如下所示：


figure8_10.png

从图中可以看出，二阶多项式核函数SVM的拟合效果更好。平均交叉验证数据集评分可以达到0.950，最高时达到0.975（绿色实线为平均值，绿色阴影上边界为最大值，下边界为最小值）。从输出的运行时间可以看出，二阶多项式核函数计算代价很高。

前面我们使用逻辑回归算法来处理乳腺癌检测问题时，使用二阶多项式增加特征，同时使用L1范数作为正则项，其拟合效果比这里的支持向量机效果好。更重要的是，逻辑回归算法的运算效率远远高于二阶多项式核函数的支持向量机算法。当然，这里的支持向量机算法的效果还是比使用L2范数作为正则项的逻辑回归算法好的。由此可见，模型选择和模型参数调优，在工程实践中有着非常重要的作用的。

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 199,340评论 5赞 467
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 83,762评论 2赞 376
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 146,329评论 0赞 329
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 53,678评论 1赞 270
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 62,583评论 5赞 359
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 47,995评论 1赞 275
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,493评论 3赞 390
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,145评论 0赞 254
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,293评论 1赞 294
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,250评论 2赞 317
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,267评论 1赞 328
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 32,973评论 3赞 316
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,556评论 3赞 303
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,648评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,873评论 1赞 255
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,257评论 2赞 345
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 41,809评论 2赞 339

推荐阅读更多精彩内容