本文详细介绍了支持向量机（SVM）的基本概念，包括基于最大间隔的分隔数据和寻找最大间隔的优化问题。重点讲解了SVM的SMO优化算法，包括Platt的SMO算法和简化版SMO算法的实现。此外，还讨论了在复杂数据上应用核函数，特别是径向基核函数的使用，以及如何在手写识别问题中应用SVM。最后，总结了SVM的优点和在实际应用中的注意事项。
支持向量机
- 说明：本章节所有代码使用版本为Python3
优点：泛化错误率低，计算开销不大，结果易解释
缺点：对参数调节和核函数的选择敏感，原始分类器不加修改仅适用于处理二类问题
适用数据类型：数值型和标称型数据
（有人说）SVM是最好的现成的分类器，“现成”是指分类器不加修饰即可直接使用。这也就意味着数据上应用基本的SVM分类器就可以得到低错误率的结果。SVM对训练集之外的数据点做出很好的分类决策。
基于最大间隔的分隔数据
两组数据分隔的足够开，很容易用一条直线将两组数据点分开，这组数据被称为线性可分数据，将数据集分隔开来的直线称为分割超平面。因此此时的分割超平面是一条直线。若数据集是三维的，此时用来分隔的是一个平面，更高维度的情况以此类推…n维数据集需要n-1维的对象来进行分割，该对象称为超平面，也就是决策的边界。
支持向量就是离分隔超平面最近的那些点，接下来要试着最大化支持向量到分隔面的距离，需要找到此问题的最优化求解方法。
寻找最大间隔
分割超平面的形式可以写成 w T x + b w^Tx+b wTx+b。计算A点到分隔超平面的距离，必须给出分隔面的法线或垂直的长度，该值为： ∣ w T x + b ∣ ⋅ 1 ∥ w ∥ |w^Tx+b|\cdot \dfrac {1}{\parallel w\parallel} ∣wTx+b∣⋅∥w∥1，这里的常熟b，类似于Logistics回归中的截距 w w w。
1、分类器求解的优化问题
分类器工作原理：输入数据给分类器会输出一个类别的标签，这相当于是一个类似于Sigmoid的函数的作用。下面使用类似海维赛德阶跃函数（单位阶跃函数）最函数
w
T
x
+
b
w^Tx+b
wTx+b作用得到
f
(
w
T
x
+
b
)
f(w^Tx+b)
f(wTx+b)，其中当u<0时
f
(
u
)
f(u)
f(u) 输出-1，反之输出+1。
当计算点到分隔面的距离并确定分隔面的位置时，间隔通过
l
a
b
e
l
∗
(
w
T
x
+
b
)
label*(w^Tx+b)
label∗(wTx+b)来计算，若数据点处于正方向（即+1类）且距离分割超平面很远的位置时，
w
T
x
+
b
w^Tx+b
wTx+b是一个很大的正数，同时
l
a
b
e
l
∗
(
w
T
x
+
b
)
label*(w^Tx+b)
label∗(wTx+b)也是一个很大的正数。若数据点处于负方向（即-1类）且距离分割超平面很远的位置时，分类标签为-1，
l
a
b
e
l
∗
(
w
T
x
+
b
)
label*(w^Tx+b)
label∗(wTx+b)也是一个很大的正数。
现在目标就是找到
w
w
w和
b
b
b。先找到具有最小间隔的数据点，这些数据点也就是前面提到的支持向量。找到具有最小间隔的数据点需要对该间隔最大化，可以写作：
a
r
g
m
a
x
w
,
b
{
m
i
n
n
(
l
a
b
e
l
⋅
(
w
T
x
+
b
)
)
⋅
1
∥
w
∥
}
arg\, \mathop{max}\limits_{w,b} \{ \mathop{min}\limits_{n}(label\cdot (w^Tx+b))\cdot \dfrac {1}{\parallel w\parallel} \}
argw,bmax{nmin(label⋅(wTx+b))⋅∥w∥1}
若令所有支持向量的
l
a
b
e
l
∗
(
w
T
x
+
b
)
label*(w^Tx+b)
label∗(wTx+b)都为1，那么可以通过计算
∥
w
∥
−
1
\parallel w\parallel ^{-1}
∥w∥−1的最大值来得到最终的解。但是并非所有数据点的
l
a
b
e
l
∗
(
w
T
x
+
b
)
label*(w^Tx+b)
label∗(wTx+b)都为1，只有那些离分隔超平面最近的点得到的值才为1。而离超平面越远的数据点
l
a
b
e
l
∗
(
w
T
x
+
b
)
label*(w^Tx+b)
label∗(wTx+b)值越大。
因此该问题是一个带约束条件的优化问题，约束条件是
l
a
b
e
l
∗
(
w
T
x
+
b
)
≥
1.0
label*(w^Tx+b)\geq 1.0
label∗(wTx+b)≥1.0。对于这类问题有一个非常著名的求解方法，即拉格朗日乘子法。通过引入拉格朗日乘子，可以基于约束条件来表达原来的问题。由于约束条件都是基于数据点的，因此可以将超平面写成数据点的形式，最终写成：
m
a
x
a
[
∑
i
=
1
m
a
−
1
2
∑
i
,
j
=
1
m
l
a
b
e
l
(
i
)
⋅
l
a
b
e
l
(
j
)
⋅
a
i
⋅
a
j
⟨
x
(
i
)
,
x
(
j
)
⟩
]
\mathop{max}\limits_{a}[\sum\limits_{i=1}^m a -\dfrac{1}{2}\sum\limits_{i,j=1}^m label^{(i)}\cdot label^{(j)}\cdot a_i\cdot a_j \langle x^{(i)},x^{(j)} \rangle]
amax[i=1∑ma−21i,j=1∑mlabel(i)⋅label(j)⋅ai⋅aj⟨x(i),x(j)⟩]
其约束条件是
a
≥
0
,
∑
i
=
1
m
a
i
⋅
l
a
b
e
l
(
i
)
=
0
a\geq 0 \quad ,\quad \sum\limits_{i=1}^m a_i\cdot label^{(i)}=0
a≥0,i=1∑mai⋅label(i)=0
但是这里有一个假设：数据必须线性可分。松弛变量：允许有些数据点可以处于分隔面的错误一侧，这样优化目标能保持不变，此时约束条件变为：
C
≥
a
≥
0
,
∑
i
=
1
m
a
i
⋅
l
a
b
e
l
(
i
)
=
0
C \geq a\geq 0 \quad ,\quad \sum\limits_{i=1}^m a_i\cdot label^{(i)}=0
C≥a≥0,i=1∑mai⋅label(i)=0
这里的常量C用于控制“最大化间隔”和“保证大部分数据点的函数间隔小于1.0”这两个目标的权重。一旦求出了所有的alpha，那么分隔超平面就可以通过这些alpha来表达。SVM的主要工作就是求解这些alpha。
2、SVM应用的一般框架
SVM的一般流程
- 收集数据：可以使用任何方法
- 准备数据：需要数值型数据
- 分析数据：有助于可视化分隔超平面
- 训练算法：SVM的大部分时间都源于训练，该过程主要实现两个参数的调优
- 测试算法：简单的计算过程就可以实现
- 使用算法：几乎所有分类问题都可以使用SVM，SVM本身是一个二类分类器，对多类问题应用SVM需要对代码做一些修改
SMO高效优化算法
下面对两个式子进行优化：一个是最小化的目标函数，另一个是在优化过程中必须遵守的约束条件。
1、Platt的SMO算法
SMO表示序列最小优化，Platt的SMO算法是将大优化算法分解成小优化问题来求解。这些小优化问题往往很容易求解，并且对他们进行顺序求解的结果与将他们作为整体来求解的结果完全一致。在结果完全一致时，SMO算法求解时间短很多。
SMO算法的目标是求出一系列的alpha和b，一旦求出了这些alpha就很容易求出权重向量w并得到分隔超平面。
SMO算法的工作原理：每次循环中选择两个alpha进行优化，一旦找到一对合适的alpha就增大其中一个，同时减小另一个。
2、应用简化版SMO算法处理小规模数据集
Platt SMO算法的完整实现要大量的代码，所以先实现简化版，过程：首先在数据集上遍历所有的alpha，然后在剩下的alpha集合中随机选择另一个alpha，从而构建alpha对。我们要同时改变两个alpha，因为有一个约束条件：
∑
i
=
1
m
a
i
⋅
l
a
b
e
l
(
i
)
=
0
\sum\limits_{i=1}^m a_i\cdot label^{(i)}=0
i=1∑mai⋅label(i)=0
3、完整示例代码
from operator import mul
import numpy as np
from numpy import *
from numpy.core.defchararray import multiply
# 获取文件中数据
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
# 随机选择J
def selectJrand(i, m):
j = i
while j == i:
j = int(np.random.uniform(0, m))
# 随机选取j，从一个均匀分布[low,high)中随机采样，注意定义域是左闭右开，即包含low，不包含high
return j
# 调整alpha的值
def clipAlpha(aj, H, L):
if aj > H:
aj = H
# alpha的值大于H，则变成H
if L > aj:
aj = L
# alpha的值小于L，则变成L
return aj
# 简化版SMO算法
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
# 5个参数分别是：数据集，类别标签，常数，容错率，最大循环次数
dataMatrix = np.mat(dataMatIn)
# mat()用来创建矩阵dataMatrix
labelMat = np.mat(classLabels).transpose()
# 将classLabels数组转化成二维的矩阵，.transpose()对其进行转置，此时labelMat是一个列向量（不是列表）
b = 0
m, n = shape(dataMatrix)
# 获取dataMatrix的行列数量
alphas = np.mat(zeros((m, 1)))
# 创建一个100(此处m大小)*1的0元素矩阵
iter = 0
# 记录循环次数，没有任何alpha改变的情况下遍历矩阵的次数
while iter < maxIter:
alphaPairsChange = 0
# alphaPairsChange用于记录alpha是否已经优化
for i in range(m):
fXi = float(np.multiply(alphas, labelMat).T *
(dataMatrix*dataMatrix[i, :].T))+b
# np.multiply(alphas,labelMat).T *dataMatrix得到 w 利用拉格朗日乘子法得到
# 第 i 个样本预测的类别，推导见《机器学习》（周志华）公式6.12
Ei = fXi-float(labelMat[i])
# 误差，误差很大，可以对该数据实例所对应的alpha值进行优化
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
# 若alpha可以更改进入优化的过程，不管是正间隔，还是负间隔都会被测试，要保证alpha不能等于0或C
j = selectJrand(i, m)
# 随机选择第二个alpha
fXj = float(np.multiply(alphas, labelMat).T *
(dataMatrix*dataMatrix[j, :].T))+b
# 与fXi一样
Ej = fXj-float(labelMat[j])
# 与Ei相同
alphaIold = alphas[i].copy()
# 保留第一个原先的alpha值，这里需要注意不能直接 alphaIold = alphas[i]
# 否则alphas[i]和alphaIold指向的都是同一内存空间
alphaJold = alphas[j].copy()
# 保留第二个原先的alpha值，方便新旧值之间的比较
if labelMat[i] != labelMat[j]:
# 如果两个alpha对应样本的标签不相同
L = max(0, alphas[j]-alphas[i])
H = min(C, C+alphas[j]-alphas[i])
# 求出相应的上下边界，else同理
else:
L = max(0, alphas[j]+alphas[i]-C)
H = min(C, C+alphas[j]+alphas[i])
# 从上面最近的 if 开始到这几行，是为了保证alpha在 0 到 C 之间
if L == H:
print("L==H")
continue
# 若是 L == H 不做任何改变
eta = 2.0*dataMatrix[i, :]*dataMatrix[j, :].T-dataMatrix[i,\
:]*dataMatrix[i, :].T-dataMatrix[j, :]*dataMatrix[j, :].T
# eta是alphas[j]的最优修改量
if eta >= 0:
# 如果eta>=0,跳出本次循环
print("eta>=0")
continue
# 下面两行代码是计算新的alphas[j]
alphas[j] -= labelMat[j]*(Ei-Ej)/eta
# 更新alpha[j]
alphas[j] = clipAlpha(alphas[j], H, L)
# 修剪alpha[j]，保证alphas[j]在区间[L, H]里面
if abs(alphas[j]-alphaJold) < 0.00001:
print("j not moving enough")
continue
# 如果改变后的alpha[j]值变化不大，跳出本次循环
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold-alphas[j])
# 对alphas[i]进行修改，修改量与 j 相同，但是方向相反
# labelMat[i]与labelMat[j]绝对值均为1，则alphas[i]与alphas[j]改变大小一样
# 保证alpha[i] * labelMal[i] + alpha[j] * labelMal[j] = C
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
# b1与b2都是设置常数项
if 0 < alphas[i] and C > alphas[i]:
b = b1
# 如果0<alpha[i]<C,那么b=b1
elif 0 < alphas[j] and C > alphas[j]:
b = b2
# 如果0<alpha[j]<C,b=b2
else:
b = (b1+b2)/2.0
# 若没有则取均值
alphaPairsChange += 1
# 统计优化次数
print("Iter: %d i: %d, pairs changed %d" %
(iter, i, alphaPairsChange))
if alphaPairsChange == 0:
# 最后判断是否有改变的alpha对，没有就进行下一次迭代
iter += 1
else:
# 否则，迭代次数置0，继续循环
iter = 0
print("iteration number: %d" % iter)
return b, alphas
# 返回最后的b值和alpha向量
dataArr, labelArr = loadDataSet('testSet.txt')
smoSimple(dataArr, labelArr, 0.6, 0.001, 10)
利用完整Platt SMO算法加速优化
在几百个点组成的小规模数据集上，简化版SMO算法的运行是没有问题的，但是在更大的数据集上运行速度就会变慢。完整的Platt SMO算法唯一不同的就是选择alpha的方式。完整版Platt SMO算法应用了一些能够提高启发的方法。
Platt SMO算法是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种之间进行交替：一种是在所有数据集上进行单遍扫描，另一种方式是在非边界alpha中实现单遍扫描，所谓非边界alpha指的是那些不等于边界的0或C的alpha值。对整个数据集的扫描比较容易，而实现非边界alpha值的扫描，首先需要建立这些alpha值的列表，然后在对这个表进行遍历。同时该步骤还会跳过那些已知不会改变的alpha的值。
在选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha值，在优化过程中会通过最大化步长的方式来获取第二个alpha值。
完整示例代码
import numpy as np
from numpy import *
import matplotlib.pyplot as plt
# 创建optStruct对象
class optStruct:
# 在初始化的时候增加dataSetIn,classLabels,C,toler这几个属性分别是：数据集，类别标签，常数，容错率
def __init__(self, dataSetIn, classLabels, C, toler):
self.X = dataSetIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataSetIn)[0]
self.alphas = np.mat(zeros((self.m, 1)))
self.b = 0
self.eCache = np.mat(zeros((self.m, 2)))
# eCache误差缓存m*2的矩阵，在以后第一列存的是标志位，第二列存的是实际E值
# 这里创建对象的目的是为了作为一个数据结构来使用对象
# 获取文件中数据
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
# 完整版Platt SMO算法的支持函数（start）
# 计算E值并返回，实际就是将简化版的这部分代码又进行了一步封装
def calcEK(oS, k):
fXk = float(np.multiply(oS.alphas, oS.labelMat).T*(oS.X*oS.X[k, :].T))+oS.b
# np.multiply(alphas,labelMat).T *dataMatrix得到 w 利用拉格朗日乘子法得到
# 第 i 个样本预测的类别，推导见《机器学习》（周志华）公式6.12
Ek = fXk-float(oS.labelMat[k])
# 误差，误差很大，可以对该数据实例所对应的alpha值进行优化
return Ek
# 随机选择J，用于选择第二个alpha或者是内循环的alpha，具有启发式
def selectJ(i, oS, Ei):
maxK = -1
# 用于保存临时最大索引
maxDeltaE = 0
# 用于保存临时最大差值--->|Ei-Ej|
Ej = 0
# 保存我们需要的Ej误差值
oS.eCache[i] = [1, Ei]
# 将Ei在缓存中设为有效值1
# 重点：这里是把SMO最后一步（根据最新阈值b，来更新Ei)提到第一步来进行了，所以这一步是非常重要的
validECacheList = nonzero(oS.eCache[:, 0].A)[0]
# nonzero(oS.eCache[:,0].A)[0]构建了一个非零表
# nonzero返回了非零E值中对应的alpha值
if len(validECacheList) > 1:
# 如果有效误差缓存长度大于1(因为包括Ei),则正常进行获取j值
for k in validECacheList:
if k == i:
continue
# 相同则不处理
Ek = calcEK(oS, k)
# 开始计算Ek值，进行对比，获取最大差值
deltaE = abs(Ei-Ek)
if deltaE > maxDeltaE:
# 更新Ej及其索引位置
maxK = k
maxDeltaE = deltaE
Ej = Ek
# 以上四行是为了选择具有最大步长的j
return maxK, Ej
else:
# 否则使用selectJradn方法选取一个随机J值
j = selectJrand(i, oS.m)
Ej = calcEK(oS, j)
return j, Ej
# 实现更新Ek操作,因为除了最后我们需要更新Ei之外，我们在内循环中计算alpha1与alpha2时还是需要用到E1与E2，
# 因为每次的E1与E2由于上一次循环中更新了alpha值，所以这一次也是需要更新E1与E2值，所以单独实现一个更新Ek值的方法还是有必要的
def updateEk(oS, k):
Ek = calcEK(oS, k)
oS.eCache[k] = [1, Ek]
# 第一列1，表示为有效标识
# 完整版Platt SMO算法的支持函数（end）
# 随机选择J
def selectJrand(i, m):
j = i
while j == i:
j = int(np.random.uniform(0, m))
# 随机选取j，从一个均匀分布[low,high)中随机采样，注意定义域是左闭右开，即包含low，不包含high
return j
# 调整alpha的值
def clipAlpha(aj, H, L):
if aj > H:
aj = H
# alpha的值大于H，则变成H
if L > aj:
aj = L
# alpha的值小于L，则变成L
return aj
# https://blog.youkuaiyun.com/weixin_43578672/article/details/109988605
# 笔者参考上面的链接的文章写出来的注释，因为innerL函数与smoP函数并不是很理解
# 完整Platt SMO算法优化内循环
# 实现内循环函数，相比于外循环，这里包含了主要的更新操作
def innerL(i, oS):
# 由外循环提供i值（具体选取要违背kkT<这里实现>,使用交替遍历<外循环中实现>）---提供α_1的索引
Ei = calcEK(oS, i)
# 计算E1值，主要是为了下面KKT条件需要使用到
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:
# 注意：对于硬间隔，我们直接和1对比，对于软间隔，我们要和1 +或- ε对比
# 如果下面违背了KKT条件，则正常进行alpha、Ek、b的更新，重点：后面单独说明下面是否满足违反KKT条件
j, Ej = selectJ(i, oS, Ei)
# 开始在内循环中，选取差值最大的alpha2下标索引
alphaIoid = oS.alphas[i].copy()
alphaJoid = oS.alphas[j].copy()
# 分析约束条件（是对所有alpha都适用），一会对我们新的alpha2进行截取纠正，注意：alpha1是由alpha2推出的，所以不需要进行验证了。
if oS.labelMat[i] != oS.labelMat[j]:
# 如果标签1与标签2异号时：
L = max(0, oS.alphas[j]-oS.alphas[i])
H = min(oS.C, oS.C+oS.alphas[j]-oS.alphas[i])
else:
# 标签1与标签2同号时
L = max(0, oS.alphas[j]+oS.alphas[i]-oS.C)
H = min(oS.C, oS.C+oS.alphas[j]+oS.alphas[i])
# 从上面最近的 if 开始到这几行，是为了保证alpha在 0 到 C 之间
if L == H:
print("L==H")
return 0
# 若是 L == H 不做任何改变
eta = 2.0*oS.X[i, :]*oS.X[j, :].T-oS.X[i,
:]*oS.X[i, :].T-oS.X[j, :]*oS.X[j, :].T
# 计算η值=k_11+k_22-2k_12，eta是alphas[j]的最优修改量
if eta >= 0:
# 如果eta>=0,跳出本次循环
print("eta>=0")
return 0
# 当上面所有条件都满足以后，我们开始正式修改alphas2值，并更新对应的Ek值
oS.alphas[j] -= oS.labelMat[j]*(Ei-Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
updateEk(oS, j)
# 更新误差缓存
if abs(oS.alphas[j]-alphaJoid) < 0.00001:
# 查看alpha2是否有足够的变化量，如果没有足够变化量，我们直接返回，不进行下面更新alpha1
# 注意：因为alpha2变化量较小，所以我们没有必要非得把值变回原来的旧值
print("j not moving enough")
return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJoid-oS.alphas[j])
# 对alphas[i]进行修改，修改量与 j 相同，但是方向相反
# labelMat[i]与labelMat[j]绝对值均为1，则alphas[i]与alphas[j]改变大小一样
# 保证alpha[i] * labelMal[i] + alpha[j] * labelMal[j] = C
updateEk(oS, i)
# 更新误差缓存
b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.X[i, :]*oS.X[i,
:].T-oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*oS.X[i, :]*oS.X[j, :].T
b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.X[i, :]*oS.X[j,
:].T-oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*oS.X[j, :]*oS.X[j, :].T
# 以上两行开始更新阈值b,正好使用到了上面更新的Ek值，b1与b2都是设置常数项
if 0 < oS.alphas[i] and oS.C > oS.alphas[i]:
oS.b = b1
# 如果0<alpha[i]<C,那么b=b1
elif 0 < oS.alphas[j] and oS.C > oS.alphas[j]:
oS.b = b2
# 如果0<alpha[j]<C,b=b2
else:
b = (b1+b2)/2.0
# 若没有则取均值
# 注意：这里我们应该根据b_new更新一次Ei，但是我们这里没有写，因为我们将这一步提前到了最开始，即selectJ中
return 1
# 以上全部更新完毕，开始返回标识
else:
return 0
# 完整Platt SMO算法优化外循环
# 交替遍历：
# 交替是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间交替：
# 一种方式是在所有数据集上进行单遍扫描，
# 另一种方式则是在非边界alpha中实现单遍扫描，所谓非边界alpha指的是那些不等于边界0或C的alpha值。
# 对整个数据集的扫描相当容易
# 而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后对这个表进行遍历。
# 同时，该步骤会跳过那些已知不变的alpha值。
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
oS = optStruct(np.mat(dataMatIn), np.mat(
classLabels).transpose(), C, toler)
# 初始化oS对象
iter = 0
entireSet = True
# 标志是否应该遍历整个数据集
alphaPairsChanged = 0
# 标志一次循环中alpha更新的次数
while iter < maxIter and (alphaPairsChanged > 0 or entireSet):
# 开始进行迭代，当iter >= maxIter或者((alphaPairsChanged == 0) and not entireSet)退出循环
# 前半个判断条件很好理解，后面的判断条件中，表示上一次循环中，是在整个数据集中遍历，并且没有α值更新过，则退出
alphaPairsChanged = 0
if entireSet:
# entireSet是true，则在整个数据集上进行遍历
for i in range(oS.m):
alphaPairsChanged += innerL(i, oS)
# 进入内循环，并将返回结果与原来的相加，计算更新次数
print("FullSet, iter: %d i: %d, pairs changes %d" %
(iter, i, alphaPairsChanged))
iter += 1
# 无论是否更新过，我们都计算迭代一次
else:
# 遍历非边界值
nonBoundIs = 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, pair changed %d" %
(iter, i, alphaPairsChanged))
# 非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d
iter += 1
# 无论是否更新过，我们都计算迭代一次
if entireSet:
entireSet = False
elif alphaPairsChanged == 0:
entireSet = True
# 如果是在非边界上，并且alpha更新过。则entireSet还是False,下一次还是在非边界上进行遍历。
# 可以认为这里是倾向于非边界遍历，因为非边界遍历的样本更符合内循环中的违反KKT条件
print("Iteration number: %d" % iter)
# 迭代次数: %d
return oS.b, oS.alphas
def calcWs(alphas, dataArr, classLabels):
# 根据西瓜书6.37求W
X = np.mat(dataArr)
labelMat = np.mat(classLabels).transpose()
m, n = shape(X)
w = zeros((n, 1))
for i in range(m):
w += np.multiply(alphas[i]*labelMat[i], X[i, :].T)
return w
# 画出拟合曲线
def plotFit(dataMat, labelMat, alphas, b):
dataArr = np.array(dataMat)
n = np.shape(dataArr)[0]
xcord1, ycord1 = [], []
xcord2, ycord2 = [], []
for i in range(n):
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i, 0])
ycord1.append(dataArr[i, 1])
else:
xcord2.append(dataArr[i, 0])
ycord2.append(dataArr[i, 1])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
x = np.arange(-1.0, 8.0, 0.1)
weights = np.multiply(alphas, np.mat(
labelMat).transpose()).T * np.mat(dataMat)
weights = np.array(weights)[0]
y = (-np.array(b)[0][0]-weights[0]*x)/weights[1]
ax.plot(x, y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()
dataArr, labelArr = loadDataSet('testSet.txt')
b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
plotFit(dataArr, labelArr, alphas, b)
完
整
版
P
l
a
t
t
S
M
O
算
法
结
果
完整版Platt SMO算法结果
完整版PlattSMO算法结果
在复杂数据上应用核函数
核函数（kernel）目的是：将数据转换成易于分类器理解的形式。
1、利用核函数将数据映射到高维空间
从一个特征空间到另一个特征空间的映射，这种映射会将低维特征空间映射到高维空间。这种从某个特征空间到另一个特征空间的映射是通过核函数来实现的。核函数相当于一个包装器或者接口，它能把数据从某个很难处理的形式转换成一个较容易处理的形式。
SVM优化中一个特别好的地方就是所有的运算都可以写成内积（也称点积）的形式。向量的内积是两个向量相乘，之后得到单个标量或者数值。把内积运算换成核函数的方式称为核技巧或者核变化
2、径向基核函数
径向基核函数是SVM中常用的一个核函数。径向基核函数是一个采用向量作为自变量的函数，能够基于向量距离运算输出一个标量。这个距离可以是从<0，0>向量或者其他向量开始计算的距离。径向基核函数的高斯版本：
k
(
x
,
y
)
=
e
x
p
⟮
∥
x
−
y
∥
2
2
σ
2
⟯
k(x,y)=exp\lgroup\dfrac {{\parallel x-y\parallel}^2}{2\sigma^2} \rgroup
k(x,y)=exp⟮2σ2∥x−y∥2⟯
其中，
σ
\sigma
σ是用户定义的确定达到达率或者说函数值跌落到0的速度参数。
上述高斯核函数将数据从其特征空间映射到更高维的空间，具体说是映射到一个无穷维的空间。
3、测试中使用核函数
构建一个对数据点进行有效分类的分类器该分类器使用了径向基核函数，径向基核函数有一个用户定义的输入 σ \sigma σ，我么需要先确定它的大小，再利用该核函数构建一个分类器。
4、完整示例代码
import numpy as np
from numpy import *
import matplotlib.pyplot as plt
# 核转换函数
def kernelTrans(X, A, kTup):
# kTup元组给出核函数信息，元组第一个是描述核函数类型的字符串
m, n = shape(X)
# 获得矩阵的行列数
K = np.mat(zeros((m, 1)))
# 构建一个m*1的零元素矩阵
if kTup[0] == 'lin':
# 在线性核函数情况下，在所有数据与数据集中的一行两个输入之间展开
K = X*A.T
# 最简单的内积方式
# 根据核函数不同类型的信息，进行相应的操作
# 要是有很多的类型直接在添加if-else实现
elif kTup[0] == 'rbf':
# 在径向基核函数情况下
for j in range(m):
deltaRow = X[j, :]-A
K[j] = deltaRow*deltaRow.T
# 循环对每一个元素计算高斯函数值
K = np.exp(K/(-1*kTup[1]**2))
# 矩阵元素之间的乘法
else:
raise NameError(
'Houston We Have a Problem -- That Kernel is not recognized')
# 前面两者都不是，就会抛出异常
return K
class optStruct:
# 在初始化的时候增加dataSetIn,classLabels,C,toler,kTup这几个属性分别是
# 数据集，类别标签，常数，容错率包含核函数信息的元组
def __init__(self, dataSetIn, classLabels, C, toler, kTup):
self.X = dataSetIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataSetIn)[0]
self.alphas = np.mat(zeros((self.m, 1)))
self.b = 0
self.eCache = np.mat(zeros((self.m, 2)))
# 误差缓存，第一列给出的是eCache是否有效的标志位，第二列给出的是实际的E值。
self.K = np.mat(zeros((self.m, self.m)))
# m行m列的核函数矩阵
for i in range(self.m):
self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)
def innerL(i, oS):
# 由外循环提供i值（具体选取要违背kkT<这里实现>,使用交替遍历<外循环中实现>）---提供α_1的索引
Ei = calcEK(oS, i)
# 计算E1值，主要是为了下面KKT条件需要使用到
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:
# 注意：对于硬间隔，我们直接和1对比，对于软间隔，我们要和1 +或- ε对比
# 如果下面违背了KKT条件，则正常进行alpha、Ek、b的更新，重点：后面单独说明下面是否满足违反KKT条件
j, Ej = selectJ(i, oS, Ei)
# 开始在内循环中，选取差值最大的alpha2下标索引
alphaIoid = oS.alphas[i].copy()
alphaJoid = oS.alphas[j].copy()
# 分析约束条件（是对所有alpha都适用），一会对我们新的alpha2进行截取纠正，注意：alpha1是由alpha2推出的，所以不需要进行验证了。
if oS.labelMat[i] != oS.labelMat[j]:
# 如果标签1与标签2异号时：
L = max(0, oS.alphas[j]-oS.alphas[i])
H = min(oS.C, oS.C+oS.alphas[j]-oS.alphas[i])
else:
# 标签1与标签2同号时
L = max(0, oS.alphas[j]+oS.alphas[i]-oS.C)
H = min(oS.C, oS.C+oS.alphas[j]+oS.alphas[i])
# 从上面最近的 if 开始到这几行，是为了保证alpha在 0 到 C 之间
if L == H:
print("L==H")
return 0
# 若是 L == H 不做任何改变
# ------------------修改的地方（start）①----------------------------
eta = 2.0*oS.K[i, j]-oS.K[i, i]-oS.K[j, j]
# 计算η值=k_11+k_22-2k_12，eta是alphas[j]的最优修改量
# ------------------修改的地方（end）---①---------------------------
if eta >= 0:
# 如果eta>=0,跳出本次循环
print("eta>=0")
return 0
# 当上面所有条件都满足以后，我们开始正式修改alphas2值，并更新对应的Ek值
oS.alphas[j] -= oS.labelMat[j]*(Ei-Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
updateEk(oS, j)
# 更新误差缓存
if abs(oS.alphas[j]-alphaJoid) < 0.00001:
# 查看alpha2是否有足够的变化量，如果没有足够变化量，我们直接返回，不进行下面更新alpha1
# 注意：因为alpha2变化量较小，所以我们没有必要非得把值变回原来的旧值
print("j not moving enough")
return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJoid-oS.alphas[j])
# 对alphas[i]进行修改，修改量与 j 相同，但是方向相反
# labelMat[i]与labelMat[j]绝对值均为1，则alphas[i]与alphas[j]改变大小一样
# 保证alpha[i] * labelMal[i] + alpha[j] * labelMal[j] = C
updateEk(oS, i)
# 更新误差缓存
# -----------------------------------------修改的地方（start）②--------------------------------------------
b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.K[i, i] - \
oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*oS.K[i, j]
b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.K[i, j] - \
oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*oS.K[j, j]
# 以上两行开始更新阈值b,正好使用到了上面更新的Ek值，b1与b2都是设置常数项
# -----------------------------------------修改的地方（end）--②--------------------------------------------
if 0 < oS.alphas[i] and oS.C > oS.alphas[i]:
oS.b = b1
# 如果0<alpha[i]<C,那么b=b1
elif 0 < oS.alphas[j] and oS.C > oS.alphas[j]:
oS.b = b2
# 如果0<alpha[j]<C,b=b2
else:
b = (b1+b2)/2.0
# 若没有则取均值
# 注意：这里我们应该根据b_new更新一次Ei，但是我们这里没有写，因为我们将这一步提前到了最开始，即selectJ中
return 1
# 以上全部更新完毕，开始返回标识
else:
return 0
# 完整版Platt SMO算法的支持函数（start）
# 计算E值并返回，实际就是将简化版的这部分代码又进行了一步封装
def calcEK(oS, k):
# ------------------修改的地方（start）③----------------------------
fXk = float(np.multiply(oS.alphas, oS.labelMat).T*oS.K[:, k]+oS.b)
# np.multiply(alphas,labelMat).T *dataMatrix得到 w 利用拉格朗日乘子法得到
# 第 i 个样本预测的类别，推导见《机器学习》（周志华）公式6.12
# ------------------修改的地方（end）--③----------------------------
Ek = fXk-float(oS.labelMat[k])
# 误差，误差很大，可以对该数据实例所对应的alpha值进行优化
return Ek
# 随机选择J，用于选择第二个alpha或者是内循环的alpha，具有启发式
def selectJ(i, oS, Ei):
maxK = -1
# 用于保存临时最大索引
maxDeltaE = 0
# 用于保存临时最大差值--->|Ei-Ej|
Ej = 0
# 保存我们需要的Ej误差值
oS.eCache[i] = [1, Ei]
# 将Ei在缓存中设为有效值1
# 重点：这里是把SMO最后一步（根据最新阈值b，来更新Ei)提到第一步来进行了，所以这一步是非常重要的
validECacheList = nonzero(oS.eCache[:, 0].A)[0]
# nonzero(oS.eCache[:,0].A)[0]构建了一个非零表
# nonzero返回了非零E值中对应的alpha值
if len(validECacheList) > 1:
# 如果有效误差缓存长度大于1(因为包括Ei),则正常进行获取j值
for k in validECacheList:
if k == i:
continue
# 相同则不处理
Ek = calcEK(oS, k)
# 开始计算Ek值，进行对比，获取最大差值
deltaE = abs(Ei-Ek)
if deltaE > maxDeltaE:
# 更新Ej及其索引位置
maxK = k
maxDeltaE = deltaE
Ej = Ek
# 以上四行是为了选择具有最大步长的j
return maxK, Ej
else:
# 否则使用selectJradn方法选取一个随机J值
j = selectJrand(i, oS.m)
Ej = calcEK(oS, j)
return j, Ej
# 实现更新Ek操作,因为除了最后我们需要更新Ei之外，我们在内循环中计算alpha1与alpha2时还是需要用到E1与E2，
# 因为每次的E1与E2由于上一次循环中更新了alpha值，所以这一次也是需要更新E1与E2值，所以单独实现一个更新Ek值的方法还是有必要的
def updateEk(oS, k):
Ek = calcEK(oS, k)
oS.eCache[k] = [1, Ek]
# 第一列1，表示为有效标识
# 完整版Platt SMO算法的支持函数（end）
# 随机选择J
def selectJrand(i, m):
j = i
while j == i:
j = int(np.random.uniform(0, m))
# 随机选取j，从一个均匀分布[low,high)中随机采样，注意定义域是左闭右开，即包含low，不包含high
return j
# 调整alpha的值
def clipAlpha(aj, H, L):
if aj > H:
aj = H
# alpha的值大于H，则变成H
if L > aj:
aj = L
# alpha的值小于L，则变成L
return aj
# 获取文件中数据
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
# 完整Platt SMO算法优化外循环
# 交替遍历：
# 交替是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间交替：
# 一种方式是在所有数据集上进行单遍扫描，
# 另一种方式则是在非边界alpha中实现单遍扫描，所谓非边界alpha指的是那些不等于边界0或C的alpha值。
# 对整个数据集的扫描相当容易
# 而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后对这个表进行遍历。
# 同时，该步骤会跳过那些已知不变的alpha值。
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
oS = optStruct(np.mat(dataMatIn), np.mat(
classLabels).transpose(), C, toler, kTup)
# 初始化oS对象，与上一个代码不同，这里有kTup参数
iter = 0
entireSet = True
# 标志是否应该遍历整个数据集
alphaPairsChanged = 0
# 标志一次循环中alpha更新的次数
while iter < maxIter and (alphaPairsChanged > 0 or entireSet):
# 开始进行迭代，当iter >= maxIter或者((alphaPairsChanged == 0) and not entireSet)退出循环
# 前半个判断条件很好理解，后面的判断条件中，表示上一次循环中，是在整个数据集中遍历，并且没有α值更新过，则退出
alphaPairsChanged = 0
if entireSet:
# entireSet是true，则在整个数据集上进行遍历
for i in range(oS.m):
alphaPairsChanged += innerL(i, oS)
# 进入内循环，并将返回结果与原来的相加，计算更新次数
print("FullSet, iter: %d i: %d, pairs changes %d" %
(iter, i, alphaPairsChanged))
iter += 1
# 无论是否更新过，我们都计算迭代一次
else:
# 遍历非边界值
nonBoundIs = 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, pair changed %d" %
(iter, i, alphaPairsChanged))
# 非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d
iter += 1
# 无论是否更新过，我们都计算迭代一次
if entireSet:
entireSet = False
elif alphaPairsChanged == 0:
entireSet = True
# 如果是在非边界上，并且alpha更新过。则entireSet还是False,下一次还是在非边界上进行遍历。
# 可以认为这里是倾向于非边界遍历，因为非边界遍历的样本更符合内循环中的违反KKT条件
print("Iteration number: %d" % iter)
# 迭代次数: %d
return oS.b, oS.alphas
def calcWs(alphas, dataArr, classLabels):
# 根据西瓜书6.37求W
X = np.mat(dataArr)
labelMat = np.mat(classLabels).transpose()
m, n = shape(X)
w = zeros((n, 1))
for i in range(m):
w += np.multiply(alphas[i]*labelMat[i], X[i, :].T)
return w
def testRbf(k1=1.3):
dataArr, labelArr = loadDataSet('testSetRBF.txt')
# 获取数据
b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))
# 完整Platt SMO算法优化外循环相比于上一个代码多了一个('rbf',k1)参数
dataMat = np.mat(dataArr)
# 将列表变成矩阵dataMat
labelMat = np.mat(labelArr).transpose()
# 将列表变成矩阵labelMat，并完成转置
svInd = nonzero(alphas.A > 0)[0]
# nonzero作用是：返回数组array中非零元素的位置（数组索引）的函数
# .A作用是矩阵类型转化为数组
# 这里是将大于 1 的元素的位置的横坐标返回给svInd
sVs = dataMat[svInd]
# 构建支持向量矩阵，这里存放的是非零的alpha值
labelSV = labelMat[svInd]
# 找到对应标签的值
print("There are %d Support Vectors" % shape(sVs)[0])
m, n = shape(dataMat)
# 获取dataMat的行列数
errorCount = 0
# 记录错误个数
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], ('rbf', k1))
# 从这一句我们看到传入的第一个参数X就是sVs，往上追溯可知这是参与內积的训练数据点即支持向量
# 第二个参数A，传入的是datMat[i:0]，传入的是一条数据，即待分类的一个数据，原因是每个待分类的数据都需要和支持向量相內积
# kTup[0]第一个元素是选择核函数类型的参数，第二个kTup[1]就更简单了，其实就是径向基的那个参数
predict = kernelEval.T*np.multiply(labelSV, alphas[svInd])+b
# 预测值，和这个svm-simple类似： fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i, :].T)) + b
# 用转化后的数据与前面的alpha以及标签值进行求积，目的就是利用和函数进行分类，multiply数组对应元素位置相乘
if np.sign(predict) != np.sign(labelArr[i]):
# 使用numpy的sign(x)函数可以取数字符号，当predict>0，返回 1
# 当predict=0，返回 0 ；当predict<0，返回 -1
errorCount += 1
# 当预测值与标签不一致时，错误个数加 1
print("The Train error rate is: %f" % (float(errorCount)/m))
dataArr, labelArr = loadDataSet('testSetRBF2.txt')
# 获取数据
errorCount = 0
# 记录错误个数
dataMat = np.mat(dataArr)
# 将列表变成矩阵dataMat
labelMat = np.mat(labelArr).transpose()
# 将列表变成矩阵labelMat，并完成转置
m, n = shape(dataMat)
# 获取dataMat的行列数
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], ('rbf', k1))
# 从这一句我们看到传入的第一个参数X就是sVs，往上追溯可知这是参与內积的训练数据点即支持向量
# 第二个参数A，传入的是datMat[i:0]，传入的是一条数据，即待分类的一个数据，原因是每个待分类的数据都需要和支持向量相內积
# kTup[0]第一个元素是选择核函数类型的参数，第二个kTup[1]就更简单了，其实就是径向基的那个参数
predict = kernelEval.T*np.multiply(labelSV, alphas[svInd])+b
# 预测值，和这个svm-simple类似： fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i, :].T)) + b
# 用转化后的数据与前面的alpha以及标签值进行求积，目的就是利用和函数进行分类，multiply数组对应元素位置相乘
if np.sign(predict) != np.sign(labelArr[i]):
errorCount += 1
# 在测试数据中，当预测值与标签不一致时，错误个数加 1
print("The test error rate is: %f" % (float(errorCount)/m))
testRbf()
示例：手写识别问题回顾
基于SVM的数字识别一般流程
- 收集数据：提供的文本文件
- 准备数据：基于二值图像的构造向量
- 分析数据：对图像向量进行目测
- 训练算法：采用两种不同的核函数，并对径向基核函数采用不同的设置来运行SMO算法
- 测试算法：编写一个函数来测试不同的核函数并计算错误率
- 使用算法：一个图像识别的完整应用还需要一些图像处理的知识
第二章时，介绍了KNN算法来识别数字，但是需要保留样本，采用SVM其需要的样本数量少了很多。
完整示例代码展示
from os import listdir
import numpy as np
from numpy import *
# 将图片加载进入本地
def loadImages(dirName):
hwLabels = []
trainingFileList = listdir(dirName)
# 获取所有文件的名称
m = len(trainingFileList)
# 查看一共有多少个文件
trainingMat = zeros((m, 1024))
# 创建一个和文件数量一样的m*1024的零元素矩阵，因为每一个文件中都是32*32的0,1代码，一共有1024个元素
for i in range(m):
# traingFileList中每一个项为"A_B.txt"，其中A为手写数字，B为样例序号
fileNameStr = trainingFileList[i]
# A_B.txt
fileStr = fileNameStr.split('.')[0]
# A_B
classNumStr = int(fileStr.split('_')[0])
# A
if classNumStr == 9:
hwLabels.append(-1)
# 若是出现数值是9，则打上标签为 -1
else:
hwLabels.append(1)
# 若是出现数值不是9，则打上标签为 1
trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))
# 文件中的每一个文件二进制数据都存放到trainingMat数组中，为以后训练数据使用
return trainingMat, hwLabels
# 返回每一个文件中的原始数据和手写标签
class optStruct:
# 在初始化的时候增加dataSetIn,classLabels,C,toler,kTup这几个属性分别是
# 数据集，类别标签，常数，容错率包含核函数信息的元组
def __init__(self, dataSetIn, classLabels, C, toler, kTup):
self.X = dataSetIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataSetIn)[0]
self.alphas = np.mat(zeros((self.m, 1)))
self.b = 0
self.eCache = np.mat(zeros((self.m, 2)))
# 误差缓存，第一列给出的是eCache是否有效的标志位，第二列给出的是实际的E值。
self.K = np.mat(zeros((self.m, self.m)))
# m行m列的核函数矩阵
for i in range(self.m):
self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)
# 核转换函数
def kernelTrans(X, A, kTup):
# kTup元组给出核函数信息，元组第一个是描述核函数类型的字符串
m, n = shape(X)
# 获得矩阵的行列数
K = np.mat(zeros((m, 1)))
# 构建一个m*1的零元素矩阵
if kTup[0] == 'lin':
# 在线性核函数情况下，在所有数据与数据集中的一行两个输入之间展开
K = X*A.T
# 最简单的内积方式
# 根据核函数不同类型的信息，进行相应的操作
# 要是有很多的类型直接在添加if-else实现
elif kTup[0] == 'rbf':
# 在径向基核函数情况下
for j in range(m):
deltaRow = X[j, :]-A
K[j] = deltaRow*deltaRow.T
# 循环对每一个元素计算高斯函数值
K = np.exp(K/(-1*kTup[1]**2))
# 矩阵元素之间的乘法
else:
raise NameError(
'Houston We Have a Problem -- That Kernel is not recognized')
# 前面两者都不是，就会抛出异常
return K
# 将32*32的二进制矩阵转换成1*1024的向量
def img2vector(filename):
resultVect = np.zeros((1, 1024))
# 创建一个1*1024的矩阵
fr = open(filename)
for i in range(32):
linesStr = fr.readline()
for j in range(32):
resultVect[0, 32*i+j] = int(linesStr[j])
# 将32*32的矩阵变化为1*1024的矩阵，方便后期信息的处理
return resultVect
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
oS = optStruct(np.mat(dataMatIn), np.mat(
classLabels).transpose(), C, toler, kTup)
# 初始化oS对象，与上一个代码不同，这里有kTup参数
iter = 0
entireSet = True
# 标志是否应该遍历整个数据集
alphaPairsChanged = 0
# 标志一次循环中alpha更新的次数
while iter < maxIter and (alphaPairsChanged > 0 or entireSet):
# 开始进行迭代，当iter >= maxIter或者((alphaPairsChanged == 0) and not entireSet)退出循环
# 前半个判断条件很好理解，后面的判断条件中，表示上一次循环中，是在整个数据集中遍历，并且没有α值更新过，则退出
alphaPairsChanged = 0
if entireSet:
# entireSet是true，则在整个数据集上进行遍历
for i in range(oS.m):
alphaPairsChanged += innerL(i, oS)
# 进入内循环，并将返回结果与原来的相加，计算更新次数
print("FullSet, iter: %d i: %d, pairs changes %d" %
(iter, i, alphaPairsChanged))
iter += 1
# 无论是否更新过，我们都计算迭代一次
else:
# 遍历非边界值
nonBoundIs = 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, pair changed %d" %
(iter, i, alphaPairsChanged))
# 非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d
iter += 1
# 无论是否更新过，我们都计算迭代一次
if entireSet:
entireSet = False
elif alphaPairsChanged == 0:
entireSet = True
# 如果是在非边界上，并且alpha更新过。则entireSet还是False,下一次还是在非边界上进行遍历。
# 可以认为这里是倾向于非边界遍历，因为非边界遍历的样本更符合内循环中的违反KKT条件
print("Iteration number: %d" % iter)
# 迭代次数: %d
return oS.b, oS.alphas
def innerL(i, oS):
# 由外循环提供i值（具体选取要违背kkT<这里实现>,使用交替遍历<外循环中实现>）---提供α_1的索引
Ei = calcEK(oS, i)
# 计算E1值，主要是为了下面KKT条件需要使用到
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:
# 注意：对于硬间隔，我们直接和1对比，对于软间隔，我们要和1 +或- ε对比
# 如果下面违背了KKT条件，则正常进行alpha、Ek、b的更新，重点：后面单独说明下面是否满足违反KKT条件
j, Ej = selectJ(i, oS, Ei)
# 开始在内循环中，选取差值最大的alpha2下标索引
alphaIoid = oS.alphas[i].copy()
alphaJoid = oS.alphas[j].copy()
# 分析约束条件（是对所有alpha都适用），一会对我们新的alpha2进行截取纠正，注意：alpha1是由alpha2推出的，所以不需要进行验证了。
if oS.labelMat[i] != oS.labelMat[j]:
# 如果标签1与标签2异号时：
L = max(0, oS.alphas[j]-oS.alphas[i])
H = min(oS.C, oS.C+oS.alphas[j]-oS.alphas[i])
else:
# 标签1与标签2同号时
L = max(0, oS.alphas[j]+oS.alphas[i]-oS.C)
H = min(oS.C, oS.C+oS.alphas[j]+oS.alphas[i])
# 从上面最近的 if 开始到这几行，是为了保证alpha在 0 到 C 之间
if L == H:
print("L==H")
return 0
# 若是 L == H 不做任何改变
# ------------------修改的地方（start）①----------------------------
eta = 2.0*oS.K[i, j]-oS.K[i, i]-oS.K[j, j]
# 计算η值=k_11+k_22-2k_12，eta是alphas[j]的最优修改量
# ------------------修改的地方（end）---①---------------------------
if eta >= 0:
# 如果eta>=0,跳出本次循环
print("eta>=0")
return 0
# 当上面所有条件都满足以后，我们开始正式修改alphas2值，并更新对应的Ek值
oS.alphas[j] -= oS.labelMat[j]*(Ei-Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
updateEk(oS, j)
# 更新误差缓存
if abs(oS.alphas[j]-alphaJoid) < 0.00001:
# 查看alpha2是否有足够的变化量，如果没有足够变化量，我们直接返回，不进行下面更新alpha1
# 注意：因为alpha2变化量较小，所以我们没有必要非得把值变回原来的旧值
print("j not moving enough")
return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJoid-oS.alphas[j])
# 对alphas[i]进行修改，修改量与 j 相同，但是方向相反
# labelMat[i]与labelMat[j]绝对值均为1，则alphas[i]与alphas[j]改变大小一样
# 保证alpha[i] * labelMal[i] + alpha[j] * labelMal[j] = C
updateEk(oS, i)
# 更新误差缓存
# -----------------------------------------修改的地方（start）②--------------------------------------------
b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.K[i, i] - \
oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*oS.K[i, j]
b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIoid)*oS.K[i, j] - \
oS.labelMat[j]*(oS.alphas[j]-alphaJoid)*oS.K[j, j]
# 以上两行开始更新阈值b,正好使用到了上面更新的Ek值，b1与b2都是设置常数项
# -----------------------------------------修改的地方（end）--②--------------------------------------------
if 0 < oS.alphas[i] and oS.C > oS.alphas[i]:
oS.b = b1
# 如果0<alpha[i]<C,那么b=b1
elif 0 < oS.alphas[j] and oS.C > oS.alphas[j]:
oS.b = b2
# 如果0<alpha[j]<C,b=b2
else:
b = (b1+b2)/2.0
# 若没有则取均值
# 注意：这里我们应该根据b_new更新一次Ei，但是我们这里没有写，因为我们将这一步提前到了最开始，即selectJ中
return 1
# 以上全部更新完毕，开始返回标识
else:
return 0
# 完整版Platt SMO算法的支持函数（start）
# 计算E值并返回，实际就是将简化版的这部分代码又进行了一步封装
def calcEK(oS, k):
# ------------------修改的地方（start）③----------------------------
fXk = float(np.multiply(oS.alphas, oS.labelMat).T*oS.K[:, k]+oS.b)
# np.multiply(alphas,labelMat).T *dataMatrix得到 w 利用拉格朗日乘子法得到
# 第 i 个样本预测的类别，推导见《机器学习》（周志华）公式6.12
# ------------------修改的地方（end）--③----------------------------
Ek = fXk-float(oS.labelMat[k])
# 误差，误差很大，可以对该数据实例所对应的alpha值进行优化
return Ek
# 随机选择J，用于选择第二个alpha或者是内循环的alpha，具有启发式
def selectJ(i, oS, Ei):
maxK = -1
# 用于保存临时最大索引
maxDeltaE = 0
# 用于保存临时最大差值--->|Ei-Ej|
Ej = 0
# 保存我们需要的Ej误差值
oS.eCache[i] = [1, Ei]
# 将Ei在缓存中设为有效值1
# 重点：这里是把SMO最后一步（根据最新阈值b，来更新Ei)提到第一步来进行了，所以这一步是非常重要的
validECacheList = nonzero(oS.eCache[:, 0].A)[0]
# nonzero(oS.eCache[:,0].A)[0]构建了一个非零表
# nonzero返回了非零E值中对应的alpha值
if len(validECacheList) > 1:
# 如果有效误差缓存长度大于1(因为包括Ei),则正常进行获取j值
for k in validECacheList:
if k == i:
continue
# 相同则不处理
Ek = calcEK(oS, k)
# 开始计算Ek值，进行对比，获取最大差值
deltaE = abs(Ei-Ek)
if deltaE > maxDeltaE:
# 更新Ej及其索引位置
maxK = k
maxDeltaE = deltaE
Ej = Ek
# 以上四行是为了选择具有最大步长的j
return maxK, Ej
else:
# 否则使用selectJradn方法选取一个随机J值
j = selectJrand(i, oS.m)
Ej = calcEK(oS, j)
return j, Ej
# 实现更新Ek操作,因为除了最后我们需要更新Ei之外，我们在内循环中计算alpha1与alpha2时还是需要用到E1与E2，
# 因为每次的E1与E2由于上一次循环中更新了alpha值，所以这一次也是需要更新E1与E2值，所以单独实现一个更新Ek值的方法还是有必要的
def updateEk(oS, k):
Ek = calcEK(oS, k)
oS.eCache[k] = [1, Ek]
# 第一列1，表示为有效标识
# 完整版Platt SMO算法的支持函数（end）
# 随机选择J
def selectJrand(i, m):
j = i
while j == i:
j = int(np.random.uniform(0, m))
# 随机选取j，从一个均匀分布[low,high)中随机采样，注意定义域是左闭右开，即包含low，不包含high
return j
# 调整alpha的值
def clipAlpha(aj, H, L):
if aj > H:
aj = H
# alpha的值大于H，则变成H
if L > aj:
aj = L
# alpha的值小于L，则变成L
return aj
def testDigits(kTup=('rbf', 10)):
dataArr, labelArr = loadImages('../机器学习02/trainingDigits')
b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
dataMat = np.mat(dataArr)
# 将列表变成矩阵dataMat
labelMat = np.mat(labelArr).transpose()
# 将列表变成矩阵labelMat，并完成转置
svInd = nonzero(alphas.A > 0)[0]
# nonzero作用是：返回数组array中非零元素的位置（数组索引）的函数
# .A作用是矩阵类型转化为数组
# 这里是将大于 1 的元素的位置的横坐标返回给svInd
sVs = dataMat[svInd]
# 构建支持向量矩阵，这里存放的是非零的alpha值
labelSV = labelMat[svInd]
# 找到对应标签的值
print("There are %d Support Vectors" % shape(sVs)[0])
m, n = shape(dataMat)
# 获取dataMat的行列数
errorCount = 0
# 记录错误个数
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], kTup)
# 从这一句我们看到传入的第一个参数X就是sVs，往上追溯可知这是参与內积的训练数据点即支持向量
# 第二个参数A，传入的是datMat[i:0]，传入的是一条数据，即待分类的一个数据，原因是每个待分类的数据都需要和支持向量相內积
# kTup[0]第一个元素是选择核函数类型的参数，第二个kTup[1]就更简单了，其实就是径向基的那个参数
predict = kernelEval.T*np.multiply(labelSV, alphas[svInd])+b
# 预测值，和这个svm-simple类似： fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i, :].T)) + b
# 用转化后的数据与前面的alpha以及标签值进行求积，目的就是利用和函数进行分类，multiply数组对应元素位置相乘
if np.sign(predict) != np.sign(labelArr[i]):
# 使用numpy的sign(x)函数可以取数字符号，当predict>0，返回 1
# 当predict=0，返回 0 ；当predict<0，返回 -1
errorCount += 1
# 当预测值与标签不一致时，错误个数加 1
print("The Train error rate is: %f" % (float(errorCount)/m))
dataArr, labelArr = loadImages('../机器学习02/testDigits')
# 获取数据
errorCount = 0
# 记录错误个数
dataMat = np.mat(dataArr)
# 将列表变成矩阵dataMat
labelMat = np.mat(labelArr).transpose()
# 将列表变成矩阵labelMat，并完成转置
m, n = shape(dataMat)
# 获取dataMat的行列数
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], kTup)
# 从这一句我们看到传入的第一个参数X就是sVs，往上追溯可知这是参与內积的训练数据点即支持向量
# 第二个参数A，传入的是datMat[i:0]，传入的是一条数据，即待分类的一个数据，原因是每个待分类的数据都需要和支持向量相內积
# kTup[0]第一个元素是选择核函数类型的参数，第二个kTup[1]就更简单了，其实就是径向基的那个参数
predict = kernelEval.T*np.multiply(labelSV, alphas[svInd])+b
# 预测值，和这个svm-simple类似： fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i, :].T)) + b
# 用转化后的数据与前面的alpha以及标签值进行求积，目的就是利用和函数进行分类，multiply数组对应元素位置相乘
if np.sign(predict) != np.sign(labelArr[i]):
errorCount += 1
# 在测试数据中，当预测值与标签不一致时，错误个数加 1
print("The test error rate is: %f" % (float(errorCount)/m))
testDigits(('rbf', 20))
# 运行可能需要几分钟
提示：这个代码和上面的代码几乎一样，只是个别地方做了修改
小结
支持向量机是一种分类器。之所以称为 “机” 是因为它会产生一个二值决策结果，即它是一种决策 “机” 。支持向量机的泛化错误率较低，具有良好的学习能力，且学到的结果具有很好的推广性。这些优点使支持向量机十分流行。
支持向量机通过求解一个二次优化问题来最大化分类化间隔。引入SMO算法加快了SVM的训练速度。核方法或者是核技巧将数据（优势是非线性的数据）从一个低维空间映射到高维空间，可以将一个在低维空间中的非线性问题转换成高维空间下的线性问题来求解。
支持向量机是一个二类分类器。当用其解决多类问题时，需要额外的方法对其进行扩展。SVM的效果也对优化参数和所用核函数中的参数敏感。
提醒： 本节的数据较多，不再展示，有需要的话，请自行查找。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
