《机器学习实战》 - 支持向量机(SVM)

1. 简介

支持向量机(Support Vector Machines,SVM)

SVM有很多实现,这里仅关注其中最流行实现:序列最小优化(Sequential Minimal Optimization,SMO)。

之后,将介绍使用核函数(kernel)方式将SVM扩展到更多数据集上。

最后,回顾手写识别,考察能否通过SVM提高识别效果。

SMO: 一种求解SVM二次规划的算法

2. 基于最大间隔 分割数据

SVM:

优点:

  1. 泛化错误率低
  2. 计算开销不大
  3. 结果易解释

缺点:

  1. 对 参数调节 和 核函数选择 敏感
  2. 原始分类器 不加修改 仅 适用于处理二类问题

适用数据类型:

  1. 数值型
  2. 标称型

线性可分(linearly separable) 概念:

问:能否画一直线将原形点和方形点分开?

图:4个线性不可分的数据集

下图 将数据集分隔开来的直线称为 分隔超平面(separating hyperplane)

由于数据点在二维平面,因此此时 分割超平面 只是一直线,若数据集是三维,则分割数据的是一个平面,更高维依此类推。

超平面 即 分类的决策边界。

若 数据点 离决策边界 越远,那么预测结果越可信

图B、C、D三直线均能将数据分隔开,但哪一条最好?

是否 最小化数据点到分隔超平面的平均距离来求最佳曲线?

我们希望 找到离分隔超平面最近的点,确保它们离分隔面的距离尽可能远

点到分隔面的距离---间隔(margin)

我们希望间隔尽可能大,因为若我们犯错或者在有限数据上训练分类器的话,我们希望分类器尽可能健壮。

两个间隔的概念:

一个是点到分隔面的距离,称为点相对于分隔面的间隔;

另一个是数据集中所有点到分隔面的最小间隔的2倍,称为分类器或数据集的间隔。

一般论文书籍中所提到的“间隔”多指后者。

SVM分类器是要找最大的数据集间隔。书中没有特意区分上述两个概念,请根据上下文理解

支持向量:离分隔超平面最近的那些点 组成的向量

接下来 最大化支持向量 到 分隔面 的距离

3. 寻找最大 间隔

分隔超平面形式: \[ wx + b \] 计算 点 A到分隔超平面的距离,即 点到分隔面的法线或垂线 的长度,该值为: \[ \frac {|w^T A + b|} {|w|} \] 这里常数 b 类似于 Logistic回归中 截距 \(w_0\)

这里的向量 w 和常数 b 一起描述了所给数据的分隔线或超平面

图:A到分隔平面的距离 就是 该点到分隔面的法线长度

3.1 分类器求解 的优化问题

由于这里的约束条件都是基于数据点的,因此我们就可以将超平面写成数据点的形式。于是,优化目标函数最后可以写成:

其约束条件为:

至此,一切都很完美,但是这里有个假设:数据必须100%线性可分。目前为止,我们知道几乎所有数据都不不可能那么“干净”。这时我们就可以通过引入所谓松弛变量(slack variable),来允许有些数据点可以处于分隔面的错误一侧。这样我们的优化目标就能保持仍然不变,但是此时新的约束条件则变为:

这里的常数C用于控制“最大化间隔”和“保证大部分点的函数间隔小于1.0”这两个目标的权重。

3.2 SVM应用的一般框架

SVM的一般流程

  1. 收集数据:可以使用任意方法。
  2. 准备数据:需要数值型数据。
  3. 分析数据:有助于可视化分隔超平面。
  4. 训练算法:SVM的大部分时间都源自训练,该过程主要实现两个参数的调优。
  5. 测试算法:十分简单的计算过程就可以实现。
  6. 使用算法:几乎所有分类问题都可以使用SVM,值得一提的是,SVM本身是一个二类分类器,对多类问题应用SVM需要对代码做一些修改。

4. SMO高效优化算法

4.1 Platt的SMO算法

1996年,John Platt发布了一个称为SMO的强大算法,用于训练SVM

SMO表示序列最小优化(Sequential Minimal Optimization)

SMO算法 将大优化问题分解为多个小优化问题,对它们进行顺序求解与作为整体求解 结果一致,但时间短很多。

SMO算法目标:

求出一系列 alpha 和 b,一旦求出这些 alpha,就很容易计算出权重向量w ,并得到分隔超平面

SMO算法原理:

每次循环中选择两个 alpha 进行优化处理。

一旦找到 一对合适 alpha,就增大其中一个 同时 减小另一个

“合适”:两个 alpha 必须符合一定条件:

  1. 两个 alpha 必须在 间隔边界之外
  2. 两个 alpha 还没有进行区间化处理 或 不在边界上

4.2 应用简化版SMO算法处理小规模数据集

简化版代码虽然量少但执行速度慢。

Platt SMO算法中的外循环确定要优化的最佳alpha对。

而简化版却会跳过这一部分,首先在数据集上遍历每一个alpha,然后在剩下的alpha集合中随机选择另一个alpha,从而构建alpha对。

这里有一点相当重要,就是我们要同时改变两个alpha。之所以这样做是因为我们有一个约束条件: \[ \sum \alpha \cdot label^{(i)} = 0 \] 由于改变 一个 alpha 可能会导致 该约束条件失效,因此 我们总是同时改变两个 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
import numpy as np
from time import sleep

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):
"""
其中i是第一个alpha的下标,m是所有alpha的数目。
只要函数值不等于输入值i,函数就会进行随机选择。
"""
# 我想要选择任何一个 不等于 i 的 j
j = i
while(j==i):
j = int(random.uniform(0,m))

return j


# 用于在数值太大时对其进行调整
def clipAlpha(aj, H, L):
"""
用于调整 大于H 或 小于L 的alpha值
"""
if aj > H:
aj = H
if L > aj:
aj = L

return aj
1
2
3
4
5
dataArr,labelArr = loadDataSet('testSet.txt')

print(labelArr)

# [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0]

下面,SMO算法的第一个版本

SMO函数的伪代码:

1
2
3
4
5
6
7
8
创建一个alpha向量并将其初始化为0向量
当迭代次数小于最大迭代次数时(外循环)
对数据集中的每个数据向量(内循环):
如果该数据向量可以被优化:
随机选择另外一个数据向量
同时优化这两个向量
如果两个向量都不能被优化,退出内循环
如果所有向量都没被优化,增加迭代数目,继续下一次循环
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
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
"""
数据集、类别标签、常数C、容错率、退出前最大的循环次数

"""

dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
b = 0; m,n = np.shape(dataMatrix)
alphas = np.mat(np.zeros((m,1)))
iter = 0
while (iter < maxIter):
alphaPairsChanged = 0
for i in range(m):
fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
Ei = fXi - float(labelMat[i]) #if checks if an example violates KKT conditions
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
j = selectJrand(i,m)
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();
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 = 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] -= labelMat[j]*(Ei - Ej)/eta
alphas[j] = clipAlpha(alphas[j],H,L)
if (abs(alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); continue
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
#the update is in the oppostie direction
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
print("iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
if (alphaPairsChanged == 0): iter += 1
else: iter = 0
print("iteration number: %d" % iter)
return b,alphas

np.multiply()函数

数组和矩阵对应位置相乘,输出与相乘数组/矩阵的大小一致

1
2
3
4
5
b,alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)

print(b)
print('-'*100)
print(alphas)
1
2
3
4
# 仅仅包含大于0的值
print(alphas[alphas>0])

# [[0.12754069 0.24131263 0.36885332]]
1
2
3
4
# 支持向量的个数
np.shape(alphas[alphas>0])

# (1, 3)
1
2
3
4
5
6
7
8
# 了解哪些数据点是支持向量
for i in range(100):
if alphas[i]>0.0:
print(dataArr[i],labelArr[i])

# [4.658191, 3.507396] -1.0
# [3.457096, -0.082216] -1.0
# [6.080573, 0.418886] 1.0

图:在示例数据集上运行简化版SMO算法后得到的结果,包括画圈的支持向量与分隔超平面

5. 利用完整Platt SMO算法加速优化

在这两个版本中,实现alpha的更改和代数运算的优化环节一模一样。

在优化过程中,唯一的不同就是选择alpha的方式。完整版的Platt SMO算法应用了一些能够提速的启发方法。

​ Platt SMO算法是通过一个外循环来选择第一个alpha值的,并且其选择过程会在两种方式之间进行交替:一种方式是在所有数据集上进行单遍扫描,另一种方式则是在非边界alpha中实现单遍扫描。而所谓非边界alpha指的就是那些不等于边界0或C的alpha值。对整个数据集的扫描相当容易,而实现非边界alpha值的扫描时,首先需要建立这些alpha值的列表,然后再对这个表进行遍历。同时,该步骤会跳过那些已知的不会改变的alpha值。

​ 在选择第一个alpha值后,算法会通过一个内循环来选择第二个alpha值。在优化过程中,会通过最大化步长的方式来获得第二个alpha值。在简化版SMO算法中,我们会在选择j之后计算错误率Ej。但在这里,我们会建立一个全局的缓存用于保存误差值,并从中选择使得步长或者说Ei-Ej最大的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
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):
"""
计算E值并返回
"""
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
"""
用于选择第二个alpha或者说内循环的alpha值
"""
# 内循环中的启发式方法
maxK = -1; maxDeltaE = 0; Ej = 0
# eCache的第一列给出的是eCache是否有效的标志位,而第二列给出的是实际的E值
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)
# (以下两行)选择具有最大步长的j
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]

完整版Platt SMO的支持函数

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
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):
"""
计算E值并返回
"""
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
"""
用于选择第二个alpha或者说内循环的alpha值
"""
# 内循环中的启发式方法
maxK = -1; maxDeltaE = 0; Ej = 0
# eCache的第一列给出的是eCache是否有效的标志位,而第二列给出的是实际的E值
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)
# (以下两行)选择具有最大步长的j
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]

print('完成')

接下来将简单介绍一下用于寻找决策边界的优化例程。

完整Platt SMO算法中的优化例程

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
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 = np.max(0, oS.alphas[j] - oS.alphas[i])
H = np.min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
else:
L = np.max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
H = np.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

print('完成')

完整版Platt SMO的外循环代码

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

print('完成')

Q&A

补充

参考

感谢帮助!

  • 《机器学习实战》[美] Peter Harrington