Python3：《机器学习实战》之支持向量机（4）核函数及其实现

Python3：《机器学习实战》之支持向量机（4）核函数及其实现

转载请注明作者和出处：http://blog.csdn.net/u011475210

代码地址：https://github.com/WordZzzz/ML/tree/master/Ch06

操作系统：WINDOWS 10

软件版本：python-3.6.2-amd64

编  者：WordZzzz

Python3机器学习实战之支持向量机4核函数及其实现
前言

核函数简介

径向基函数

核函数实现

手写识别问题回顾

总结

前言：

  前面三篇讲的都是SVM线性分类器，如果数据集非线性可分（如下图所示），那就需要做些修改了。



  显而易见，在该数据中存在某种可以识别的模式。接下来，我们就需要使用一种称为核函数（kernel）的工具将数据转换成易于分类器理解的形式。在本篇博文中，我们先对核函数进行简单的了解，然后重点研究径向基函数（radial basis function，最流行的核函数）及其实现，最后在我们之前的手写数字识别问题上进行实践。

核函数简介

  再看本文章开头的图片，我们似乎可以用一个圆来吧数据划分开来，但是对于线性分类器来说，这好像很难实现。我们或许可以对数据进行某种形式的转换，从而得到某些新的变量来表示数据。在这个例子中，我们将数据从一个特征空间转换到另一个特征空间，在新空间下，我们可以很容易利用已有的工具对数据进行处理。这个过程，学过数学的都知道，即从一个特征空间到另一个特征空间的映射。通常情况下，核函数实现的是低维到高维的映射，以后其他算法涉及到的PCA等，则是高维到低微的映射。经过空间转换之后，我们可以在高维空间解决线性问题，这也就等价于在低维空间中解决非线性问题。

  SVM优化中一个特别好的地方就是，所有的运算都可以写成内积（inner product，也叫点积）的形式。向量的内积指的是两个向量相乘，之后得到单个标量或者数值。我们可以把内积运算替换成核函数，而不必做简化处理。将内积替换成核函数的方式被称为核技巧（kernel trick）或者核“变电”（kernel substation）。

  当然，核函数并不仅仅应用于SVM中，很多其他的机器学习算法也都用到核函数。接下来，我们就来介绍一个流行的核函数，那就是径向基函数。

径向基函数

  径向基函数是一个采用向量作为自变量的函数，能够基于向量距离运算输出一个标量。这个距离可以使从<0,0>向量或者其他向量开始计算的距离。我们用到的径向基函数的高斯版本公式为：

k(x,y)=exp(−∥x−y∥22σ2)

  其中，σ是用户定义的用于确定到达率（reach）或者说是函数值跌落到0的速度参数。

  上述高斯核函数将数据从其特征空间映射到跟高维的空间，具体来说这里是映射到一个无穷维的空间。在该数据集上，使用高斯核函数得到一个很好的结果，当然，该函数也可以用于许多其他的数据集，并且也能够得到低错误率的结果。

核函数实现

  如果在svmMLiA.py文件中添加一个函数并稍作修改，那么我们就能在已有代码中使用核函数了（所有与核函数实现相关的函数，函数名末尾都是K）。其中主要区分代码在innerLK()和calcEk()中，我已经重点标记出。

  新添加的核转换函数，主要用于填充结构体和后续的计算：

def kernelTrans(X, A, kTup):
"""
Function：   核转换函数

Input：      X：数据集
A：某一行数据
kTup：核函数信息

Output： K：计算出的核向量
"""
#获取数据集行列数
m, n = shape(X)
#初始化列向量
K = mat(zeros((m, 1)))
#根据键值选择相应核函数
#lin表示的是线性核函数
if kTup[0] == 'lin': K = X * A.T
#rbf表示径向基核函数
elif kTup[0] == 'rbf':
for j in range(m):
deltaRow = X[j,:] - A
K[j] = deltaRow * deltaRow.T
#对矩阵元素展开计算，而不像在MATLAB中一样计算矩阵的逆
K =  exp(K/(-1*kTup[1]**2))
#如果无法识别，就报错
else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
#返回计算出的核向量
return K

  其他函数：

class optStructK:
"""
Function：   存放运算中重要的值

Input：      dataMatIn：数据集
classLabels：类别标签
C：常数C
toler：容错率
kTup：速度参数

Output： X：数据集
labelMat：类别标签
C：常数C
tol：容错率
m：数据集行数
b：常数项
alphas：alphas矩阵
eCache：误差缓存
K：核函数矩阵
"""
def __init__(self, dataMatIn, classLabels, C, toler, kTup):
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataMatIn)[0]
self.alphas = mat(zeros((self.m, 1)))
self.b = 0
self.eCache = mat(zeros((self.m, 2)))

""" 主要区分 """
self.K = mat(zeros((self.m, self.m)))
for i in range(self.m):
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
""" 主要区分 """

def calcEkK(oS, k):
"""
Function：   计算误差值E

Input：      oS：数据结构
k：下标

Output： Ek：计算的E值
"""

""" 主要区分 """
#计算fXk，整个对应输出公式f(x)=w`x + b
#fXk = float(multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k,:].T)) + oS.b
fXk = float(multiply(oS.alphas, oS.labelMat).T*oS.K[:, k] + oS.b)
""" 主要区分 """

#计算E值
Ek = fXk - float(oS.labelMat[k])
#返回计算的误差值E
return Ek

def selectJK(i, oS, Ei):
"""
Function：   选择第二个alpha的值

Input：      i：第一个alpha的下标
oS：数据结构
Ei：计算出的第一个alpha的误差值

Output： j：第二个alpha的下标
Ej：计算出的第二个alpha的误差值
"""
#初始化参数值
maxK = -1; maxDeltaE = 0; Ej = 0
#构建误差缓存
oS.eCache[i] = [1, Ei]
#构建一个非零列表，返回值是第一个非零E所对应的alpha值，而不是E本身
validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
#如果列表长度大于1，说明不是第一次循环
if (len(validEcacheList)) > 1:
#遍历列表中所有元素
for k in validEcacheList:
#如果是第一个alpha的下标，就跳出本次循环
if k == i: continue
#计算k下标对应的误差值
Ek = calcEkK(oS, k)
#取两个alpha误差值的差值的绝对值
deltaE = abs(Ei - Ek)
#最大值更新
if (deltaE > maxDeltaE):
maxK = k; maxDeltaE = deltaE; Ej = Ek
#返回最大差值的下标maxK和误差值Ej
return maxK, Ej
#如果是第一次循环，则随机选择alpha，然后计算误差
else:
j = selectJrand(i, oS.m)
Ej = calcEkK(oS, j)
#返回下标j和其对应的误差Ej
return j, Ej

def updateEkK(oS, k):
"""
Function：   更新误差缓存

Input：      oS：数据结构
j：alpha的下标

Output： 无
"""
#计算下表为k的参数的误差
Ek = calcEkK(oS, k)
#将误差放入缓存
oS.eCache[k] = [1, Ek]

def innerLK(i, oS):
"""
Function：   完整SMO算法中的优化例程

Input：      oS：数据结构
i：alpha的下标

Output： 无
"""
#计算误差
Ei = calcEkK(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)):
#启发式选择第二个alpha值
j, Ej = selectJK(i, oS, Ei)
#利用copy存储刚才的计算值，便于后期比较
alphaIold = oS.alphas[i].copy(); alpahJold = oS.alphas[j].copy();
#保证alpha在0和C之间
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.X[i, :]*oS.X[j, :].T - oS.X[i, :]*oS.X[i, :].T - oS.X[j, :]*oS.X[j, :].T
eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
""" 主要区分 """

#如果最优修改量大于0，则不做处理直接跳出本次循环，这里对真实SMO做了简化处理
if eta >= 0: print("eta>=0"); return 0
#计算新的alphas[j]的值
oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
#对新的alphas[j]进行阈值处理
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
#更新误差缓存
updateEkK(oS, j)
#如果新旧值差很小，则不做处理跳出本次循环
if (abs(oS.alphas[j] - alpahJold) < 0.00001): print("j not moving enough"); return 0
#对i进行修改，修改量相同，但是方向相反
oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alpahJold - oS.alphas[j])
#更新误差缓存
updateEkK(oS, i)

""" 主要区分 """
#更新常数项
#b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :]*oS.X[i, :].T - oS.labelMat[j] * (oS.alphas[j] - alpahJold) * oS.X[i, :]*oS.X[j, :].T
#b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :]*oS.X[j, :].T - oS.labelMat[j] * (oS.alphas[j] - alpahJold) * oS.X[j, :]*oS.X[j, :].T
b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, i] - oS.labelMat[j] * (oS.alphas[j] - alpahJold) * 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] - alpahJold) * oS.K[j, j]
""" 主要区分 """

#谁在0到C之间，就听谁的，否则就取平均值
if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[i]): oS.b = b2
else: oS.b = (b1 + b2) / 2.0
#成功返回1
return 1
#失败返回0
else: return 0

def smoPK(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin', 0)):
"""
Function：   完整SMO算法

Input：      dataMatIn：数据集
classLabels：类别标签
C：常数C
toler：容错率
maxIter：最大的循环次数
kTup：速度参数

Output： b：常数项
alphas：数据向量
"""
#新建数据结构对象
oS = optStructK(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
#初始化迭代次数
iter = 0
#初始化标志位
entireSet = True; alphaPairsChanged = 0
#终止条件：迭代次数超限、遍历整个集合都未对alpha进行修改
while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
alphaPairsChanged = 0
#根据标志位选择不同的遍历方式
if entireSet:
#遍历任意可能的alpha值
for i in range(oS.m):
#选择第二个alpha值，并在可能时对其进行优化处理
alphaPairsChanged += innerLK(i, oS)
print("fullSet, iter: %d i: %d, pairs changed %d" % (iter, i, alphaPairsChanged))
#迭代次数累加
iter += 1
else:
#得出所有的非边界alpha值
nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
#遍历所有的非边界alpha值
for i in nonBoundIs:
#选择第二个alpha值，并在可能时对其进行优化处理
alphaPairsChanged += innerLK(i, oS)
print("non-bound, iter: %d i: %d, pairs changed %d" % (iter, i, alphaPairsChanged))
#迭代次数累加
iter += 1
#在非边界循环和完整遍历之间进行切换
if entireSet: entireSet = False
elif (alphaPairsChanged == 0): entireSet =True
print("iteration number: %d" % iter)
#返回常数项和数据向量
return oS.b, oS.alphas

  接下来我们写测试函数。整个代码中最重要的是for循环开始的那两行，他们给出了如何利用核函数进行分类。首先利用结构初始化方法中使用过的kernelTrans()函数，得到转换后的数据。然后，再用其与前面的alpha及类别标签值求积。特别需要注意观察的是，我们是如何做到只需要支持向量数据就可以进行分类的。

def testRbf(k1 = 1.3):
"""
Function：   利用核函数进行分类的径向基测试函数

Input：      k1：径向基函数的速度参数

Output： 输出打印信息
"""
#导入数据集
dataArr, labelArr = loadDataSet('testSetRBF.txt')
#调用Platt SMO算法
b, alphas = smoPK(dataArr, labelArr, 200, 0.00001, 10000, ('rbf', k1))
#初始化数据矩阵和标签向量
datMat = mat(dataArr); labelMat = mat(labelArr).transpose()
#记录支持向量序号
svInd = nonzero(alphas.A > 0)[0]
#读取支持向量
sVs = datMat[svInd]
#读取支持向量对应标签
labelSV = labelMat[svInd]
#输出打印信息
print("there are %d Support Vectors" % shape(sVs)[0])
#获取数据集行列值
m, n = shape(datMat)
#初始化误差计数
errorCount = 0
#遍历每一行，利用核函数对训练集进行分类
for i in range(m):
#利用核函数转换数据
kernelEval = kernelTrans(sVs, datMat[i,:], ('rbf', k1))
#仅用支持向量预测分类
predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
#预测分类结果与标签不符则错误计数加一
if sign(predict) != sign(labelArr[i]): errorCount += 1
#打印输出分类错误率
print("the training error rate is: %f" % (float(errorCount)/m))
#导入测试数据集
dataArr, labelArr = loadDataSet('testSetRBF2.txt')
#初始化误差计数
errorCount = 0
#初始化数据矩阵和标签向量
datMat = mat(dataArr); labelMat = mat(labelArr).transpose()
#获取数据集行列值
m, n = shape(datMat)
#遍历每一行，利用核函数对测试集进行分类
for i in range(m):
#利用核函数转换数据
kernelEval = kernelTrans(sVs, datMat[i,:], ('rbf', k1))
#仅用支持向量预测分类
predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
#预测分类结果与标签不符则错误计数加一
if sign(predict) != sign(labelArr[i]): errorCount += 1
#打印输出分类错误率
print("the test error rate is: %f" % (float(errorCount)/m))

  上述代码分别在训练集和测试集上进行性能测试，打印输出如下：

>>> reload(svmMLiA)
<module 'svmMLiA' from 'E:\\机器学习实战\\mycode\\Ch06\\svmMLiA.py'>
>>> svmMLiA.testRbf()
L==H
fullSet, iter: 0 i: 0, pairs changed 0
fullSet, iter: 0 i: 1, pairs changed 1
fullSet, iter: 0 i: 2, pairs changed 2
fullSet, iter: 0 i: 3, pairs changed 3
···
fullSet, iter: 6 i: 96, pairs changed 0
fullSet, iter: 6 i: 97, pairs changed 0
fullSet, iter: 6 i: 98, pairs changed 0
fullSet, iter: 6 i: 99, pairs changed 0
iteration number: 7
there are 27 Support Vectors
the training error rate is: 0.030000
the test error rate is: 0.040000

  当然，大家也可以尝试更换不同的k1参数以观察测试错误率、训练错误率、支持向量个数随k1的变化情况。下面个两张图是书上的举例。





  我们会发现，支持向量的数目存在一个最优值。SVM的优点在于它能对数据进行高效分类。如果支持向量太少，就可能会得到一个很差的决策边界；如果支持向量太多，也就是相当于每次都利用整个数据集进行分类，这种分类方法成为kNN（多么熟悉）。

手写识别问题回顾

  SVM对kNN的优点在于，SVM只需要保留支持向量就可以获得可比的效果，占用内存大大减小。下面，我们就用SVM来对手写数字进行分类识别（不加修改的SVM只能用于二分类问题，在这里，我们规定，如果是数字9，则为-1，否则为+1）。

def img2vector(filename):
"""
Function：   32*32图像转换为1*1024向量

Input：      filename：文件名称字符串

Output： returnVect：转换之后的1*1024向量
"""
#初始化要返回的1*1024向量
returnVect = zeros((1, 1024))
#打开文件
fr = open(filename)
#读取文件信息
for i in range(32):
#循环读取文件的前32行
lineStr = fr.readline()
for j in range(32):
#将每行的头32个字符存储到要返回的向量中
returnVect[0, 32*i+j] = int(lineStr[j])
#返回要输出的1*1024向量
return returnVect

def loadImages(dirName):
"""
Function：   加载图片

Input：      dirName：文件路径

Output： trainingMat：训练数据集
hwLabels：数据标签
"""
from os import listdir
#初始化数据标签
hwLabels = []
#读取文件列表
trainingFileList = listdir(dirName)
#读取文件个数
m = len(trainingFileList)
#初始化训练数据集
trainingMat = zeros((m,1024))
#填充数据集
for i in range(m):
#遍历所有文件
fileNameStr = trainingFileList[i]
#提取文件名称
fileStr = fileNameStr.split('.')[0]
#提取数字标识
classNumStr = int(fileStr.split('_')[0])
#数字9记为-1类
if classNumStr == 9: hwLabels.append(-1)
#其他数字记为+1类
else: hwLabels.append(1)
#提取图像向量，填充数据集
trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr))
#返回数据集和数据标签
return trainingMat, hwLabels

def testDigits(kTup = ('rbf',10)):
"""
Function：   手写数字分类函数

Input：      kTup：核函数采用径向基函数

Output： 输出打印信息
"""
#导入数据集
dataArr, labelArr = loadImages('trainingDigits')
#调用Platt SMO算法
b, alphas = smoPK(dataArr, labelArr, 200, 0.0001, 10000, kTup)
#初始化数据矩阵和标签向量
datMat = mat(dataArr); labelMat = mat(labelArr).transpose()
#记录支持向量序号
svInd = nonzero(alphas.A > 0)[0]
#读取支持向量
sVs = datMat[svInd]
#读取支持向量对应标签
labelSV = labelMat[svInd]
#输出打印信息
print("there are %d Support Vectors" % shape(sVs)[0])
#获取数据集行列值
m, n = shape(datMat)
#初始化误差计数
errorCount = 0
#遍历每一行，利用核函数对训练集进行分类
for i in range(m):
#利用核函数转换数据
kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
#仅用支持向量预测分类
predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
#预测分类结果与标签不符则错误计数加一
if sign(predict)!=sign(labelArr[i]): errorCount += 1
#打印输出分类错误率
print "the training error rate is: %f" % (float(errorCount)/m)
#导入测试数据集
dataArr,labelArr = loadImages('testDigits')
#初始化误差计数
errorCount = 0
#初始化数据矩阵和标签向量
datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
#获取数据集行列值
m,n = shape(datMat)
#遍历每一行，利用核函数对测试集进行分类
for i in range(m):
#利用核函数转换数据
kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
#仅用支持向量预测分类
predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
#预测分类结果与标签不符则错误计数加一
if sign(predict)!=sign(labelArr[i]): errorCount += 1
#打印输出分类错误率
print "the test error rate is: %f" % (float(errorCount)/m)

  上面的大部分代码我们都已经很熟悉了，这里就不再赘述，直接进行测试。

>>> reload(svmMLiA)
>>> svmMLiA.testDigits(('rbf', 20))
L==H
fullSet, iter: 6 i: 398, pairs changed 0
j not moving enough
fullSet, iter: 6 i: 399, pairs changed 0
fullSet, iter: 6 i: 400, pairs changed 0
j not moving enough
fullSet, iter: 6 i: 401, pairs changed 0
iteration number: 7
there are 58 Support Vectors
the training error rate is: 0.000000
the test error rate is: 0.021505

  尝试不同的σ，并尝试了线性核函数，总结得到如图结果：

Kernel,settings	Training error (%)	Test error (%)	Support vectors
RBF,0.1	0	52	402
RBF,5	0	3.2	402
RBF,10	0	0.5	99
RBF,50	0.2	2.2	41
RBF,100	1.5	4.3	26
Linear	2.7	2.2	28

  结果表明，当径向基函数的参数σ取10左右时，就可以得到最小的测试错误率。同时，最小的训练错误率，并不对应于最小的支持向量数目。另外，线性核函数的效果并不是特别糟糕，可以以牺牲线性核函数的错误率来换取分类速度的提高。

总结

  支持向量机是一种分类器。之所以称为“机”是因为它会产生一个二值决策结果，即它是一种决策“机”。支持向量机的泛化错误率较低，具有良好的学习能力，且学到的结果具有很好的推广性。这些优点使得支持向量机十分流行，有些人认为他是监督学习中最好的定式算法。

系列教程持续发布中，欢迎订阅、关注、收藏、评论、点赞哦～～(￣▽￣～)～

完的汪(∪｡∪)｡｡｡zzz