5-支持向量机
🍧 支持向量机 SVM
本章首先讲解 支持向量机(Support Vector Machines, SVM) 的基本概念，书中会引入一些关键术语。SVM 有很多实现，但是本章只关注其中最流行的一种实现，即序列最小化 Sequential Minimal Optimization，SMO 算法。在此之后，将介绍如何使用一种称为核函数 kernel 的方式将 SVM 扩展到更多数据集上。最后，会回顾我们在 kNN 算法中手写识别的例子，利用 SVM 算法提高识别的效果。
1. 线性可分支持向量机与硬间隔最大化
① 基于最大间隔分割数据
在介绍 SVM 之前，先解释几个概念。
-
线性可分 linearly separable
下图是线性可分的数据集：
如下图便是线性不可分的数据集：
-
分割超平面 separating hyperplane
上述将数据集分割开来的直线称为分割超平面，也就是分类的决策边界。由于上述例子中的数据点都在二维平面上，所以此时分割超平面就只是一条直线。如果所给的数据集是三维的，那么此时用来分割数据的就是一个平面。依次类推，如果数据集是 n 维的，那么分割超平面就是 n-1 维的。
-
间隔 margin
如果数据点离决策边界越远，那么其最后的预测结果也就越可信。最佳方案就是我们确保离分割超平面最近的数据点，离分割面的距离尽可能的远。
这里点到分割面的距离就是间隔。我们希望间隔尽可能的大。
💡 硬间隔 hard margin 就是所有样本必须划分正确，软间隔就是允许错误分类。
-
支持向量 support vector：
支持向量就是离分割超平面最近的那些点。⭐ 我们的目标就是最大化支持向量到分割面的距离。
② SVM 优化目标与约束条件
如何求解线性可分数据集的最佳分割直线？
分割超平面的形式可以写成 $wx + b$。要计算点 A 到分割超平面的距离，就必须给出点到分割面的法线或垂线的长度，该值为 $\frac{|wA +b|}{||w||}$
对分割超平面 $wx + b$ 作用得到函数 $f(x) = sign(wx + b)$ ，该函数实现分类：当$wx + b < 0$ 时函数 f 输出 -1，当$wx + b > 0$ 时函数 f 输出 +1。
函数 $sign(wx + b)$ 就是分类器 / 分类决策函数。
💡 类别标签用 -1、1，是因为仅差一个符号，方便数学上的处理。我们可以用一个统一的公式来表示间隔或者数据点到分割超平面的距离。
针对所有的训练数据，SVM 都希望：
这时就能体现 -1 和 +1 类的好处了，因为无论间隔是多少，都可以依靠伸缩 w 和 b 约为 1。
🚩 在这个图中，决策边界两边的线之间的距离（最大间隔）为：
它是这样求出来的，数据中的支持向量在影响着最大间隔，那么假设两个支持向量 x1 和 x2 分别为正负，最大间隔就应该是 x2-x1 在法向量上的投影：
⭐ 所以求取最大间隔 d 的过程为：
这就是 SVM 的优化目标，它想要找到 max(d) 最大间隔，然后大家可能发现了，这个目标里面没有 b，只和 w 有关，❓ 那么是不是任意的 b 都可以呢？
显然不是的，这个优化有一个约束条件，因为推导的过程就有假设条件是两个支持向量要求在两侧，所以这个 ⭐ 约束条件可以写成：（ i = 1…….n，n为样本的个数）
💡 在 $wx_i + b = ± 1$ 边界的点就是支持向量
🚀 显然，为了最大化间隔，仅需最大化 $||w||^{-1}$，最大化 $||w||^{-1}$ 等价于最小化 $||w||^2$，以下便是 SVM 的最优化问题：
③ 学习的对偶算法
Ⅰ 对偶问题
在上述优化问题中，我们给定了一些约束条件然后求最优值，因此该问题就是一个带约束条件的优化问题。这里的约束条件就是 $y_i * (wx_i + b) ≥ 1.0$，目标就是使得 $\frac{1}{2} ||w||^2$ 尽可能的小。对于这类优化问题，有一个著名的求解方法：拉格朗日乘子法。
根据拉格朗日乘子法：就是求函数 f(x1,x2,…) 在 g(x1,x2,…)=0 的约束条件下的极值的方法。其主要思想是将约束条件函数与原函数联系到一起，使能配成与变量数量相等的等式方程，从而求出得到原函数极值的各个变量的解。即可以列出下式：
其中 $α = (a_1,a_2,...,a_n)^T$ 就是拉格朗日乘子法引入的新的参数，也就是拉格朗日乘子。
那么问题就变成了：
所谓的对偶问题就是：
做这种转换是为了后面的求解方便。✍ 所以，为了得到对偶问题的解，需要先对 $L(w,b,\alpha)$ 的 w，b 求极小，再对 $\alpha$ 求极大
1）求 $min_{w,b}L(w,b,\alpha)$
将拉格朗日函数对分别对 w，b 求偏导 =0 并令其等于 0 ：
💡 偏导符号也可写成 $▽_w$ ，表示对 w 偏导
在这里求出了两个结果，带入到 L(w,b,α) 中消除 w 和 b，将其转换成 α 单变量问题：
🔈 $(x_i,x_j)$ 表示两个向量的内积，即 $x_i^Tx_j$
所以问题被转化成为：
2）求 $min_{w,b}L(w,b,\alpha)$ 对 α 的极大，即是对偶问题：
⭐ 将上式目标函数由极大转换成极小，就得到下面与之等价的对偶最优化问题：
Ⅱ 线性可分支持向量机学习算法
假设上述对偶最优化问题对 α 的解为 $\alpha^* = (a_1*,a_2...a_N*)T$，可以由 $α^$ 求得对 w, b 的解 $w*，b*$，有下面的定理：
🏃 关于此定理的证明可以参考《统计学习方法》
由此定理可知，分离超平面可以写成：
分类决策函数可以写成：
也就是说，分类决策函数只依赖于输入 x 和训练样本输入的内积。上式分类决策函数称为线性可分支持向量机的对偶形式。
🚩 上述过程就是线性可分 SVM 学习算法：
2. 线性支持向量机与软间隔最大化
① 引入松弛变量 slack variable
从前面内容来看，SVM理论可以完美的找到正负样本间的最大的分类间隔，这意味着不仅仅可以实现对训练数据的分类，还保证了决策平面是最理想的。那么SVM为什么还要引入松弛因子呢？我们可以看下面这个例子：
实线的决策平面分开了所有的正负样本，但是造成的结果是最大分类间隔非常小，而虚线虽然没有完全区分开所有的正负样本，但是保证了分类间隔更大，这个例子可以说明，不一定实线完全分类的决策平面就是最好的。
同时，还有一种情况，如果左上角的圆再向右侧移动的话，那么样本的数据本身就是线性不可分的。
综合上面两点，我们可以通过引入一个松弛变量 $ε_i ≥ 0$ 来允许数据点可以处于分隔面错误的一侧。这样就能保证我们的优化目标仍能保持不变，此时的约束条件变为：
很显然，由于松弛变量是一个正数，那么新的约束条件一定没有原来的条件“严格”，也就是松弛了。
同时，对每个松弛变量 $ε_i$ ，支付一个 代价 C，这样，目标函数由原来的 $1/2 * ||w||^2$ 变成：
除了一阶软间隔分类器，SVM的目标函数还有二阶软间隔分类器的形式：
而无论是上述哪种情况，实际上都是为目标函数引入一个损失，而上面的参数 C > 0 被称之为惩罚因子。
惩罚因子的大小决定了对离群点的容忍程度（松弛的程度），比如，如果 C 是一个非常大的值，那么很小的松弛因子 $ε_i$ 都会让目标函数的值变得很大，由于目标是在求一个最小值，这就意味着哪怕是个很小的松弛都不愿意容忍。
可以拿一个很直观的例子说明惩罚因子C的影响，C 越大意味着对训练数据而言能得到很好的分类结果，但是同时最大分类间隔就会变小，毕竟我们做模型不是为了在训练数据上有个多么优异的结果。相反的，如果 C 比较小，那么间隔就会变大，模型也就有了相对而言更好的泛化能力。
💡 例如: 正类有10000个样本，而负类只给了100个（C 越大表示100个负样本的影响越大，就会出现过度拟合，所以 C 决定了负样本对模型拟合程度的影响，C 就是一个非常关键的优化点！）
② 优化目标与约束条件
⭐ 最后新的目标函数（采用一阶范式）：
🚩 最小化目标 函数包含两层含义：使 $1/2 * ||w||^2$ 尽量小即间隔尽量大，同时使误分类点的个数尽量小， C 是调和二者的系数。
由上式求出的分类决策函数模型称为训练样本线性不可分时的线性支持向量机，简称为线性支持向量机。显然，线性支持向量机包含线性可分支持向量机。
③ 学习的对偶算法
Ⅰ 对偶问题
和前面一样，有了目标函数之后，后面的工作就是引入拉格朗日乘子，对偶问题转化后求导，解出拉格朗日乘子后反带回去求w和b，这个讨论对于引得目标函数来说是一样的，而且推导的过程在前面四部分中解释的已经非常清楚，对于新的问题，只是加了两个参数：
-
目标函数本身多出来的一个松弛因子 $ε$
-
拉格朗日乘子法为新的约束条件引入的参数 C
需要注意的是，参数 C 是一个已知量，是人为设定的！
OK，轻车熟路，我们的拉格朗日函数如下：
其中，α ≥ 0 和 u ≥ 0 是拉格朗日乘子，对偶后就变成了：
1）首先求 L 对 w，b，$ε$ 的极小：
将三式带入到 L 中求得：
2）再对 minL 求 α 的极大，即得对偶问题：
此时的约束条件为：
这个约束条件中的后三个是可以转化的，因为 u 是一个大于等于0的数（拉格朗日性质决定的），同时由于 $C-a=u$，即 $C-a>=0$。同时由于 $a>=0$，后三个条件最后就变成：
⭐ 那么最后整理得到的对偶问题就是：
可以转变符号改为 min 问题：
这个结果和线性可分 SVM 的对偶问题的最后的结果很像，只是多出了些约束。
Ⅱ 线性支持向量机学习算法
假设上述对偶最优化问题对 α 的解为 $\alpha^* = (a_1*,a_2...a_N*)T$，可以由 $α^$ 求得对 w, b 的解 $w*，b*$，有下面的定理：
🏃 关于此定理的证明可以参考《统计学习方法》
上述过程就是线性 SVM 学习算法：
3. 序列最小最优化算法 SMO
对于我们上述讨论的问题，SVM 有很多种实现，最流行的一种实现是: 序列最小优化(Sequential Minimal Optimization, SMO) 算法。
SMO 算法要解决的对偶问题是：
此处我们将坐标点上升到高维空间，利用核技巧，用核函数 $K(x_i,x_j)$ 代替 $(x_i,x_j)$ （核函数见下文详解），对于整个推导结果没有影响。即：
在这个问题中，变量是拉格朗日乘子，⭐ 一个变量 $\alpha_i$ 对应于一个样本点 $(x_i,y_i)$ ，变量的总数等于训练样本容量 N。
🚩 SMO 算法：
-
创建作者: John Platt
-
创建时间: 1996年
-
SMO用途: 用于训练 SVM
-
SMO 目标: 求出一系列 α 和 b, 一旦求出 α，就很容易计算出权重向量 w 并得到分隔超平面。
-
SMO 思想: 将大优化问题分解为多个小优化问题。
-
⭐ SMO 原理: 每次循环选择两个 α 进行优化处理，一旦找出一对合适的 α，那么就增大一个同时减少一个。这里指的合适必须要符合一定的条件：
-
这两个 α 必须要在间隔边界之外；
-
这两个 α 还没有进行过区间化处理或者不在边界上。
之所以要同时改变 2 个 α，原因是我们有一个约束条件: $\sum_{i=1}^{m} a_i·y_i=0$：
如果只是修改一个 α1 ，然后在 α1 上求极值。即固定 α1 以外的所有参数，由上面这个约束条件可以知道，α1 将不再是变量（可以由其他值推出），因为我们可以计算出来：
因此，我们需要一次选取两个参数做优化，比如 αi 和 αj，此时 αi 可以由 αj 和其他参数表示出来。
这样回代入W 中，W 就只是关于 αj 的函数了，这时候就可以只对 αj 进行优化了。在这里就是对 αj 进行求导，令导数为 0 就可以解出这个时候最优的 αj 了。然后也可以得到 αi。这就是一次的迭代过程，一次迭代只调整两个拉格朗日乘子 αi 和 αj。
-
⭐ 所以，SMO的优化过程为，重复以下三步，直至收敛：
- 选择两个拉格朗日乘子 αi 和 αj ；
- 固定其他拉格朗日乘子 αk (k 不等于 i 和 j )，只对 αi 和 αj 优化，将 αi 表示成有关 αj 和 αk (k 不等于 i 和 j )的式子，则最终 w(α) 只是有关 αj 的式子，从而对 αj 求导，得到优化后的 αj，再解出 αi；
- 根据优化后的 αi 和 αj，更新阈值 b；
OK，下面详细解析这三步 👇：
① 优化两个拉格朗日乘子
假设选择的两个拉格朗日乘子是 $α_1，α_2$，其他变量 $α_i(i=3,4...N)$ 是固定的。于是 SMO 的最优化问题
可以写成：
其中，$K(i,j)$ 表示 $K(x_i,x_j)$
由于只有两个变量 $α_1，α_2$，因此约束可以用二维空间中的图形表示：
🏃 关于此定理的证明可以参考《统计学习方法》
② 选择两个拉格朗日乘子
Ⅰ 第 1 个变量的选择
⭐ SMO 称选择第一个变量的过程为外层循环。外层循环在训练样本中选取 违反 KKT 条件最严重的样本点，并将其对应的变量作为第一个变量。具体地，检验训练样本点 $(x_i,y_i)$ 是否满足 KKT 条件，即：
其中：
该检验是在松弛变量 ε 范围内进行的，外层循环首先变量所有满足条件 $0 < α_i < C$ 的样本点，记载间隔边界上的支持向量点，检验它们是否满足 KKT 条件。如果这些样本点都满足 KKT 条件，那么遍历整个训练集，检验它们是否满足 KKT 条件。
Ⅱ 第 2 个变量的选择
SMO 称选择第 2 个变量的过程为 内层循环。假设外层循环中已经找到第 1 个变量 $α_1$，现在要在内层循环中找第 2 个变量 $α_2$。第 2 个变量选择的标准是希望使 $α_2$ 有足够大的变化。
③ 计算阈值 b 和差值 Ei
每次完成两个变量的优化之后，都要重新计算阈值 b 和 $E_i$：
④ SMO 算法
⑤ 应用简化版 SMO 算法处理小规模数据集
Platt SMO 算法的完整实现需要大量代码，在接下来的一个例子中，我们将会对算法进行简化处理，以便了解算法的基本工作思路，之后再基于简化版给出完整版。
Ⅰ 辅助函数
Platt SMO 算法中的外循环确定要优化的最佳 α 对，而简化版会跳过这一部分，首先在数据集上遍历每一个 α，然后在剩下的 α 集合中随机挑选另一个 α，从而构建 α 对。
为此，我们将构建一个辅助函数，用于在某个区间范围内随机选择一个整数。同时，我们还需要另一个辅助函数，用于在数值太大时对其进行调整。
import numpy as np
def loadDataSet(fileName):
dataMat = []
labelMat = []
fr = open(fileName)
for line in fr.readlines():
lineArr = line.strip().split('\t')
dataMat.append([float(lineArr[0]), float(lineArr[1])])
labelMat.append(float(lineArr[2]))
return dataMat, labelMat
# 辅助函数，用于在某个区间范围内随机选择一个整数
def selectJrand(i, m):
"""
Desc:
用于从0~m的区间范围内选择除 i 外的一个随机整数
Args:
i: 已选择的 α 的下标
m: 所有 α 的数目
"""
j = i
while(j == i): # 希望选择的随机数j不等于 i
j = int(np.random.uniform(0,m))
return j
# 辅助函数，用于在α数值太大时对其进行调整。
def clipAlpha(aj, H, L):
"""
Desc:
用于调整大于 H 或者小于 L 的 α值
"""
if aj > H:
aj = H
if aj < L:
aj = L
return aj
我们的数据集如下，保存了上图 6-3 的数据。
dataArr, labelArr = loadDataSet('testSet.txt')
可以看的出来，这里采用的类别标签是 -1 和 +1，而不是 0 和 1 了。
OK，上述工作完成 后，就可以使用 SMO 算法的第一个版本了。
Ⅱ 简化版 SMO 函数代码
简化版 SMO 函数的伪代码大致如下:
创建一个 alpha 向量并将其初始化为0向量
当迭代次数小于最大迭代次数时(外循环)
对数据集中的每个数据向量(内循环):
如果该数据向量可以被优化
随机选择另外一个数据向量
同时优化这两个向量
如果两个向量都不能被优化，退出内循环
如果所有向量都没被优化，增加迭代数目，继续下一次循环
具体 Python 实现：
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
"""smoSimple
Args:
dataMatIn 特征集合
classLabels 类别标签
C 松弛变量(常量值)，允许有些数据点可以处于分隔面的错误一侧。
控制最大化间隔和保证大部分的函数间隔小于1.0这两个目标的权重。
可以通过调节该参数达到不同的结果。
toler 容错率（是指在某个体系中能减小一些因素或选择对某个系统产生不稳定的概率。）
maxIter 退出前最大的循环次数
Returns:
b 模型的常量值
alphas 拉格朗日乘子
"""
dataMatrix = np.mat(dataMatIn)
# 矩阵转置
labelMat = np.mat(classLabels).transpose()
m, n = np.shape(dataMatrix) # 100 x 3
# 初始化 b 和 α
b = 0
alphas = np.mat(np.zeros((m,1)))
iter = 0 # 没有任何 α 改变的情况下 遍历数据的次数
while(iter < maxIter):
# 记录 α 是否已经进行优化，每次 循环时设为 0，然后再对整个集合顺序遍历
alphaPairsChanged = 0
for i in range(m):
# 我们预测的类别 y[i] = w^T * x[i] + b
# 其中 w = Σ(1~n) α[n]*y[n]*x[n]
# fXi 即我们预测的类别
fXi = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
# 预测结果与真实结果比对，计算误差Ei
Ei = fXi - float(labelMat[i])
"""
如果误差很大，那么就对该数据实例所对应的 α进行优化
labelMat[i]*Ei 表示发生错误的概率: 如果超出了容错率 toler， 才需要优化。至于正负号，我们考虑绝对值就对了
"""
if((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
# 如果满足优化的条件，我们就随机选取非i的一个点作为第二个α，进行优化比较
j = selectJrand(i, m)
# 预测 j 的结果
fXj = float(np.multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[j, :].T)) + b
Ej = fXj - float(labelMat[j])
# 两个 α 如下
alphaIold = alphas[i].copy()
alphaJold = alphas[j].copy()
# L和H用于将alphas[j]调整到0-C之间。如果L==H，就不做任何改变，直接执行continue语句
# labelMat[i] != labelMat[j] 表示异侧，就相减，否则是同侧，就相加。
if (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
# eta是alphas[j]的最优修改量，如果eta==0，需要退出for循环的当前迭代过程
eta = 2.0 * dataMatrix[i, :]*dataMatrix[j, :].T - dataMatrix[i, :]*dataMatrix[i, :].T - dataMatrix[j, :]*dataMatrix[j, :].T
if eta >= 0:
print("eta>=0")
continue
# 计算出一个新的alphas[j]值
alphas[j] -= labelMat[j]*(Ei - Ej)/eta
# 并使用辅助函数，以及L和H对其进行调整
alphas[j] = clipAlpha(alphas[j], H, L)
# 检查alpha[j]是否只是轻微的改变，如果是的话，就退出for循环。
if (abs(alphas[j] - alphaJold) < 0.00001):
print("j not moving enough")
continue
# 然后alphas[i]和alphas[j]同样进行改变，虽然改变的大小一样，但是改变的方向正好相反
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
# 在对alpha[i], alpha[j] 进行优化之后，给这两个alpha值设置一个常数b。
# w= Σ[1~n] ai*yi*xi => b = yj- Σ[1~n] ai*yi(xi*xj)
b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i, :]*dataMatrix[i, :].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i, :]*dataMatrix[j, :].T
b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i, :]*dataMatrix[j, :].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j, :]*dataMatrix[j, :].T
if (0 < alphas[i]) and (C > alphas[i]):
b = b1
elif (0 < alphas[j]) and (C > alphas[j]):
b = b2
else:
b = (b1 + b2)/2.0
alphaPairsChanged += 1 # α优化次数+1
print("iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
# 在for循环外，检查alpha值是否做了更新，如果更新则将iter设为0后继续运行程序
# 直到更新完毕后，iter次循环无变化，才退出循环。
if (alphaPairsChanged == 0):
iter += 1
else:
iter = 0
print("iteration number: %d" % iter)
return b, alphas
我们可以直接观察 α 矩阵本身，由于其中 0 元素太多，我们只观察其中大于 0 的元素：
输出支持向量：
在原始数据集上对这些支持向量画圈后结果如下：
Ⅲ 将简化版 SMO 算法应用于分类器
OK，现在，我们已经获得了 α 和 b 的值，我们的分类超平面是 $w*x + b$ ，而函数 $f(x) = sign(w^Tx + b)$ 便是分类器：当 $w^Tx + b < 0$ 时函数 f 输出 -1 即表示该数据属于类别 -1，当$w^Tx + b > 0$ 时函数 f 输出 +1 即表示该数据属于类别 +1。
首先，我们需要基于 alpha 值得到超平面：$w*x + b$，上面已经推导过：
， α 和 b 已经知晓，求出 w 即可：
# 基于 alpha 值得到超平面中的参数 w 的值
def calcWs(alphas,dataArr,classLables):
X = np.mat(dataArr)
labelMat = np.mat(classLables).transpose()
m,n = np.shape(X)
w = np.zeros((n,1))
for i in range(m):
w += np.multiply(alphas[i]*labelMat[i],X[i,:].T)
return w
ws = calcWs(alphas,dataArr,labelArr)
接下来利用分类超平面进行分类，如果分类超平面 wx + b 的值小于 0，说明该数据属于类别 -1，如果分类超平面 wx + b 的值大于 0，说明该数据属于类别 +1。
我们利用第一个数据进行测试，首先看看我们数据集中的第一个数据：
显然，它属于类别 -1，接下来，看看我们分类器的效果：
dataMat = np.mat(dataArr)
dataMat[0] * np.mat(ws) + b # 分类超平面 wx +b
结果 <0，即属于类别 -1。正确 🎉
⑥ 利用完整 Platt SMO 算法加速优化
在小型数据集上，简化版算法的运行是没啥问题的，但是在更大的数据集上，运行速度就会变慢。下面我们讨论完整版的 Platt SMO 算法。两者之间唯一的不同就是选择 alpha 的方式。
Platt SMO 算法是通过一个外循环来选择第一个 alpha 值的，并且在其选择过程中会在两种方式之间进行交替：一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界 alpha 中实现单遍扫描（所谓非边界 alpha 就是那些不等于边界 0 或 C 的 alpha 值）。
选择完第一个 alpha 值后，算法会通过一个内循环来选择第二个 alpha 值。在优化过程中，通过最大化步长的方式来获得第二个 alpha 值。在简化版 SMO 算法中，我们会在选择 j 之后计算错误率 Ej。但在这里，我们会建立了一个全局的缓存用于保存误差值，并从中选择使得步长或者说 Ei-Ej 最大的 alpha 值。
完整版的代码实在太长了，不想看，日后再来 👻
class optStruct:
def __init__(self,dataMatIn, classLabels, C, toler, kTup): # Initialize the structure with the parameters
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = np.shape(dataMatIn)[0]
self.alphas = np.mat(np.zeros((self.m,1)))
self.b = 0
self.eCache = np.mat(np.zeros((self.m,2))) #first column is valid flag
self.K = np.mat(np.zeros((self.m,self.m)))
for i in range(self.m):
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
def calcEk(oS, k):
fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek
def selectJ(i, oS, Ei): #this is the second choice -heurstic, and calcs Ej
maxK = -1; maxDeltaE = 0; Ej = 0
oS.eCache[i] = [1,Ei] #set valid #choose the alpha that gives the maximum delta E
validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]
if (len(validEcacheList)) > 1:
for k in validEcacheList: #loop through valid Ecache values and find the one that maximizes delta E
if k == i: continue #don't calc for i, waste of time
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if (deltaE > maxDeltaE):
maxK = k; maxDeltaE = deltaE; Ej = Ek
return maxK, Ej
else: #in this case (first time around) we don't have any valid eCache values
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej
def updateEk(oS, k):#after any alpha has changed update the new value in the cache
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
def innerL(i, oS):
Ei = calcEk(oS, i)
if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
if (oS.labelMat[i] != oS.labelMat[j]):
L = max(0, oS.alphas[j] - oS.alphas[i])
H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
else:
L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
H = min(oS.C, oS.alphas[j] + oS.alphas[i])
if L==H: print("L==H"); return 0
eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
if eta >= 0: print("eta>=0"); return 0
oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
updateEk(oS, j) #added this for the Ecache
if (abs(oS.alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
updateEk(oS, i) #added this for the Ecache #the update is in the oppostie direction
b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]
b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]
if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
else: oS.b = (b1 + b2)/2.0
return 1
else: return 0
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #full Platt SMO
oS = optStruct(np.mat(dataMatIn),np.mat(classLabels).transpose(),C,toler, kTup)
iter = 0
entireSet = True; alphaPairsChanged = 0
while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
alphaPairsChanged = 0
if entireSet: #go over all
for i in range(oS.m):
alphaPairsChanged += innerL(i,oS)
print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
iter += 1
else:#go over non-bound (railed) alphas
nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerL(i,oS)
print("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
iter += 1
if entireSet: entireSet = False #toggle entire set loop
elif (alphaPairsChanged == 0): entireSet = True
print("iteration number: %d" % iter)
return oS.b,oS.alphas
6. 非线性支持向量机与核函数
考虑下图给出的数据，显而易见，在该数据中存在某种可以识别的模式。我们能否像线性可分的情况一样，利用强大的工具来捕捉数据种的这种模式？
显然，答案是肯定的，接下来，我们使用核函数 kernel 将 数据转换成易于分类器理解的形式。
本节首先解释核函数的概念，并介绍它们在 SVM 中的使用 方法，然后，介绍一种称为 径向基函数 radial bias function 的最流行的核函数。最后，将核函数应用到我们的分类器中。
① 利用核函数将数据映射到高维空间
在上图中，数据点可以用一个圆来分割（非线性模型），人类的大脑可以意识到这一点。然后，对于分类器而言，它只能识别分类器的结果是大于 0 还是小于 0。如果只在 x 和 y 轴构成的坐标系中插入直线进行分类的话，我们并不会得到理想的结果。我们可以对圆中的数据进行某种形式的转换，从而得到某些新的变量来表示数据。
在这个例子中，我们将数据从一个特征空间转换到另一个特征空间。这个过程称为映射。通常情况下，一般将低纬特征空间映射到高维特征空间。
通过核函数来实现映射。可以把核函数想象成一个包装器(wrapper)或者是接口(interface)，它能将数据从某个很难处理的形式转换成为另一个较容易处理的形式。
📜 核函数的定义如下：
👍 经过空间转换后: 低维需要解决的非线性问题，就变成了高维需要解决的线性问题。
② 核技巧 — 将内积替换成核函数
SVM 优化特别好的地方，在于所有的运算都可以写成内积 (inner product: 是指2个向量相乘，得到单个标量 或者 数值)；将内积替换成核函数的方式被称为核技巧(kernel trick)或者核"变电"(kernel substation)
比如，我们可以将线性可分 SVM 中的目标函数
利用核技巧转换为：
同样，分类决策函数中的内积也能用核函数代替：
其实也就等价于在新的特征空间中训练模型。
③ 常用的核函数
对于高斯核函数来说，其中 σ 是用户定义的用于确定到达率 reach 或者或函数值跌落到 0 的速度参数。（关于 σ 下文会有更直观的解释）
高斯核函数对应的支持向量机是**高斯径向基函数分类器 Radial Basis Function，RBF **。分类决策函数为：
# 核转换函数
def kernelTrans(X,A,kTup):
"""
核转换函数
Args:
X dataMatIn数据集
A dataMatIn数据集的第i行的数据
kTup 核函数的信息
Returns:
"""
m,n = np.shape(X)
K = np.mat(np.zeros((m,1)))
# 元组kTup给出的是核函数的信息
# 元组的第一个参数是描述所用核函数类型的一个字符串
# 其他两个参数则是核函数可能需要的可选参数
if kTup[0] == 'lin': # 线性核函数
K = X * A.T
elif kTup[0] == 'rbf': # 径向基核函数
for j in range(m):
deltaRow = X[j, :] - A # 对应公式中的 x-y
K[j] = deltaRow * deltaRow.T # 对应公式中的 ||x-y||^2
# 径向基函数的高斯版本
K = np.exp(K / (-1 * kTup[1] ** 2))
else: # 遇到无法识别的元组，程序就会抛出异常
raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
return K
④ 非线性支持向量机学习算法
利用核技巧将线性分类的学习方法应用到非线性分类问题中，将线性 SVM 扩展到 非线性 SVM，只需将 线性 SVM 对偶形式中的内积转换成核函数。
⑤ 在测试中使用核函数
接下来我们将构建一个对图 6-6 中的数据点进行有效分类的分类器，该分类器使用了高斯核函数。前面 提到的高斯核函数有一个用户定义的输入 σ，首先我们要确定它的大小，然后利用该核函数构建出一个分类器。
我们的数据集如下：
# 利用核函数进行分类的径向基测试函数
def testRbf(k1=1.3):
dataArr,labelArr = loadDataSet('testSetRBF.txt')
b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
datMat=np.mat(dataArr)
labelMat = np.mat(labelArr).transpose()
svInd = np.nonzero(alphas.A>0)[0]
sVs=datMat[svInd] # get matrix of only support vectors
labelSV = labelMat[svInd]
print("there are %d Support Vectors" % np.shape(sVs)[0])
m,n = np.shape(datMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
if np.sign(predict)!=np.sign(labelArr[i]): errorCount += 1
print("the training error rate is: %f" % (float(errorCount)/m))
dataArr,labelArr = loadDataSet('testSetRBF2.txt')
errorCount = 0
datMat=np.mat(dataArr); labelMat = np.mat(labelArr).transpose()
m,n = np.shape(datMat)
for i in range(m):
kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
if np.sign(predict)!=np.sign(labelArr[i]): errorCount += 1
print("the test error rate is: %f" % (float(errorCount)/m))
可以尝试更换不同的 k1 参数（即 σ）以观察测试错误率、训练错误率、支持向量个数的变化。
下图给出了 k1 非常小（0.1）时的支持向量的个数：
显然，错误率上升了。
当我们增大σ 的时候：
⭐ 由此可得：RBF 中参数 σ 越大，支持向量的个数越少，参数 σ 越小，支持向量的个数越多。
σ 存在一个最优值，降低它，训练错误率就会降低，减少它，测试错误率就会上升。
⭐ 即支持向量的数目存在一个最优值（或者说 σ 存在一个最优值，因为 σ 直接决定支持向量的数目），如果支持向量太少，就会得到一个很差的决策边界。如果支持向量太多，就相当于每次都用整个数据集进行分类，那这个时候的分类方法其实就是 kNN 了。
📚 References
-
《Machine Learning in Action》
-
《统计学习方法 - 李航》
-
《机器学习 - 周志华》
