《机器学习实战》之支持向量机(4)核函数及其实现



前言:

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



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

核函数简介

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

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

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

径向基函数

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

$$
k(x,y) = exp(\frac{ {-\begin{Vmatrix} x-y \end{Vmatrix} }^2}{2 \sigma^2})
$$

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

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

核函数实现

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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

  其他函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
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及类别标签值求积。特别需要注意观察的是,我们是如何做到只需要支持向量数据就可以进行分类的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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))

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
>>> 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)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
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)

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

1
2
3
4
5
6
7
8
9
10
11
12
13
>>> 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

0%