一 前言
上篇文章讲解的是线性SVM的推导过程以及简化版SMO算法的代码实现。本篇文章将讲解SMO算法的优化方法以及非线性SVM。
二 SMO算法优化
在几百个点组成的小规模数据集上,简化版SMO算法的运行是没有什么问题的,但是在更大的数据集上的运行速度就会变慢。简化版SMO算法的第二个α的选择是随机的,针对这一问题,我们可以使用启发式选择第二个α值,来达到优化效果。
1 启发选择方式
下面这两个公式想必已经不再陌生:
在实现SMO算法的时候,先计算η,再更新a_j。为了加快第二个α_j乘子的迭代速度,需要让直线的斜率增大,对于α_j的更新公式,其中η值没有什么文章可做,于是只能令:
因此,我们可以明确自己的优化方法了:
- 最外层循环,首先在样本中选择违反KKT条件的一个乘子作为最外层循环,然后用”启发式选择”选择另外一个乘子并进行这两个乘子的优化
- 在非边界乘子中寻找使得|E_i - E_j|最大的样本
- 如果没有找到,则从整个样本中随机选择一个样本
1)第一个变量αi的选择:
SMO称选择第一个变量的过程为外层循环。外层训练在训练样本中选取违法KKT条件最严重的样本点。并将其对应的变量作为第一个变量。具体的,检验训练样本(xi, yi)是否满足KKT条件,也就是:
该检验是在ε范围内进行的。在检验过程中,外层循环首先遍历所有满足条件0<αj<C的样本点,即在间隔边界上的支持向量点,检验他们是否满足KKT条件,然后选择违反KKT条件最严重的αi。如果这些样本点都满足KKT条件,那么遍历整个训练集,检验他们是否满足KKT条件,然后选择违反KKT条件最严重的αi。
优先选择遍历非边界数据样本,因为非边界数据样本更有可能需要调整,边界数据样本常常不能得到进一步调整而留在边界上。由于大部分数据样本都很明显不可能是支持向量,因此对应的α乘子一旦取得零值就无需再调整。遍历非边界数据样本并选出他们当中违反KKT 条件为止。当某一次遍历发现没有非边界数据样本得到调整时,遍历所有数据样本,以检验是否整个集合都满足KKT条件。如果整个集合的检验中又有数据样本被进一步进化,则有必要再遍历非边界数据样本。这样,不停地在遍历所有数据样本和遍历非边界数据样本之间切换,直到整个样本集合都满足KKT条件为止。以上用KKT条件对数据样本所做的检验都以达到一定精度ε就可以停止为条件。如果要求十分精确的输出算法,则往往不能很快收敛。
对整个数据集的遍历扫描相当容易,而实现对非边界αi的扫描时,首先需要将所有非边界样本的αi值(也就是满足0<αi<C)保存到新的一个列表中,然后再对其进行遍历。同时,该步骤跳过那些已知的不会改变的αi值。
2)第二个变量αj的选择:
在选择第一个αi后,算法会通过一个内循环来选择第二个αj值。因为第二个乘子的迭代步长大致正比于|Ei-Ej|,所以我们需要选择能够最大化|Ei-Ej|的第二个乘子(选择最大化迭代步长的第二个乘子)。在这里,为了节省计算时间,我们建立一个全局的缓存用于保存所有样本的误差值,而不用每次选择的时候就重新计算。我们从中选择使得步长最大或者|Ei-Ej|最大的αj
接下来,让我们看看完整版SMO算法如何实现。
2 完整版SMO算法
完整版Platt SMO算法是通过一个外循环来选择违反KKT条件的一个乘子,并且其选择过程会在这两种方式之间进行交替:
- 在所有数据集上进行单遍扫描
- 在非边界α中实现单遍扫描
非边界α指的就是那些不等于边界0或C的α值,并且跳过那些已知的不会改变的α值。所以我们要先建立这些α的列表,用于才能出α的更新状态。
在选择第一个α值后,算法会通过”启发选择方式”选择第二个α值。
3 编写代码
我们首先构建一个仅包含init方法的optStruct类,将其作为一个数据结构来使用,方便我们对于重要数据的维护。代码思路和之前的简化版SMO算法是相似的,不同之处在于增加了优化方法,如果上篇文章已经看懂,我想这个代码会很好理解。创建一个svm-smo.py文件,编写代码如下:
# -*-coding:utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import random
"""
Author:
Jack Cui
Blog:
http://blog.csdn.net/c406495762
Zhihu:
https://www.zhihu.com/people/Jack--Cui/
Modify:
2017-10-03
"""
class optStruct:
"""
数据结构,维护所有需要操作的值
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
"""
def __init__(self, dataMatIn, classLabels, C, toler):
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))) #根据矩阵行数初始化alpha参数为0
self.b = 0 #初始化b参数为0
self.eCache = np.mat(np.zeros((self.m,2))) #根据矩阵行数初始化虎误差缓存,第一列为是否有效的标志位,第二列为实际的误差E的值。
def loadDataSet(fileName):
"""
读取数据
Parameters:
fileName - 文件名
Returns:
dataMat - 数据矩阵
labelMat - 数据标签
"""
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 calcEk(oS, k):
"""
计算误差
Parameters:
oS - 数据结构
k - 标号为k的数据
Returns:
Ek - 标号为k的数据误差
"""
fXk = float(np.multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T) + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek
def selectJrand(i, m):
"""
函数说明:随机选择alpha_j的索引值
Parameters:
i - alpha_i的索引值
m - alpha参数个数
Returns:
j - alpha_j的索引值
"""
j = i #选择一个不等于i的j
while (j == i):
j = int(random.uniform(0, m))
return j
def selectJ(i, oS, Ei):
"""
内循环启发方式2
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Ei - 标号为i的数据误差
Returns:
j, maxK - 标号为j或maxK的数据的索引值
Ej - 标号为j的数据误差
"""
maxK = -1; maxDeltaE = 0; Ej = 0 #初始化
oS.eCache[i] = [1,Ei] #根据Ei更新误差缓存
validEcacheList = np.nonzero(oS.eCache[:,0].A)[0] #返回误差不为0的数据的索引值
if (len(validEcacheList)) > 1: #有不为0的误差
for k in validEcacheList: #遍历,找到最大的Ek
if k == i: continue #不计算i,浪费时间
Ek = calcEk(oS, k) #计算Ek
deltaE = abs(Ei - Ek) #计算|Ei-Ek|
if (deltaE > maxDeltaE): #找到maxDeltaE
maxK = k; maxDeltaE = deltaE; Ej = Ek
return maxK, Ej #返回maxK,Ej
else: #没有不为0的误差
j = selectJrand(i, oS.m) #随机选择alpha_j的索引值
Ej = calcEk(oS, j) #计算Ej
return j, Ej #j,Ej
def updateEk(oS, k):
"""
计算Ek,并更新误差缓存
Parameters:
oS - 数据结构
k - 标号为k的数据的索引值
Returns:
无
"""
Ek = calcEk(oS, k) #计算Ek
oS.eCache[k] = [1,Ek] #更新误差缓存
def clipAlpha(aj,H,L):
"""
修剪alpha_j
Parameters:
aj - alpha_j的值
H - alpha上限
L - alpha下限
Returns:
aj - 修剪后的alpah_j的值
"""
if aj > H:
aj = H
if L > aj:
aj = L
return aj
def innerL(i, oS):
"""
优化的SMO算法
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Returns:
1 - 有任意一对alpha值发生变化
0 - 没有任意一对alpha值发生变化或变化太小
"""
#步骤1:计算误差Ei
Ei = calcEk(oS, i)
#优化alpha,设定一定的容错率。
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)):
#使用内循环启发方式2选择alpha_j,并计算Ej
j,Ej = selectJ(i, oS, Ei)
#保存更新前的aplpha值,使用深拷贝
alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
#步骤2:计算上下界L和H
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
#步骤3:计算eta
eta = 2.0 * oS.X[i,:] * oS.X[j,:].T - oS.X[i,:] * oS.X[i,:].T - oS.X[j,:] * oS.X[j,:].T
if eta >= 0:
print("eta>=0")
return 0
#步骤4:更新alpha_j
oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta
#步骤5:修剪alpha_j
oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
#更新Ej至误差缓存
updateEk(oS, j)
if (abs(oS.alphas[j] - alphaJold) < 0.00001):
print("alpha_j变化太小")
return 0
#步骤6:更新alpha_i
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
#更新Ei至误差缓存
updateEk(oS, i)
#步骤7:更新b_1和b_2
b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*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]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
#步骤8:根据b_1和b_2更新b
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):
"""
完整的线性SMO算法
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
Returns:
oS.b - SMO算法计算的b
oS.alphas - SMO算法计算的alphas
"""
oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler) #初始化数据结构
iter = 0 #初始化当前迭代次数
entireSet = True; alphaPairsChanged = 0
while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): #遍历整个数据集都alpha也没有更新或者超过最大迭代次数,则退出循环
alphaPairsChanged = 0
if entireSet: #遍历整个数据集
for i in range(oS.m):
alphaPairsChanged += innerL(i,oS) #使用优化的SMO算法
print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
else: #遍历非边界值
nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] #遍历不在边界0和C的alpha
for i in nonBoundIs:
alphaPairsChanged += innerL(i,oS)
print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
if entireSet: #遍历一次后改为非边界遍历
entireSet = False
elif (alphaPairsChanged == 0): #如果alpha没有更新,计算全样本遍历
entireSet = True
print("迭代次数: %d" % iter)
return oS.b,oS.alphas #返回SMO算法计算的b和alphas
def showClassifer(dataMat, classLabels, w, b):
"""
分类结果可视化
Parameters:
dataMat - 数据矩阵
w - 直线法向量
b - 直线解决
Returns:
无
"""
#绘制样本点
data_plus = [] #正样本
data_minus = [] #负样本
for i in range(len(dataMat)):
if classLabels[i] > 0:
data_plus.append(dataMat[i])
else:
data_minus.append(dataMat[i])
data_plus_np = np.array(data_plus) #转换为numpy矩阵
data_minus_np = np.array(data_minus) #转换为numpy矩阵
plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7) #正样本散点图
plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) #负样本散点图
#绘制直线
x1 = max(dataMat)[0]
x2 = min(dataMat)[0]
a1, a2 = w
b = float(b)
a1 = float(a1[0])
a2 = float(a2[0])
y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
plt.plot([x1, x2], [y1, y2])
#找出支持向量点
for i, alpha in enumerate(alphas):
if