最小化的目标函数、优化过程中必须要遵循的额约束条件。不久之前,人们使用二次规划求解工具来解决上述最优化问题,这种工具是一种用于在线性约束下优化具有多个变量的二次目标函数的软件,而这些二次规划求解工具需要强大的计算能力支撑,另外在实现上也十分复杂。所有需要做的围绕优化的事情就是训练分类器,一旦得到alpha的最优值,我们就得到了分隔超平面并能够将之用于数据分类。
platt的SMO算法
SMO表示序列最小优化。platt的SMO算法时将大优化问题分解为多个小优化问题来求解的。这些小优化问题往往很容易求解,并且对它们进行顺序求解的结果与将它们作为整体来求解的结果是完全一致的。在结果完全相同的同时,SMO算法的求解时间短很多。
SMO算法的目标是求出一系列的alpha和b,一旦求出了这些alpha,就很容易计算出权重向量w并得到分隔超平面。
SMO算法的工作原理是:每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha,那么就增大其中一个同时减少另一个。这里所谓的合适是指两个alpha必须要符合一定的条件,条件之一就是这两个alpha必须要在间隔边界之外,第二个条件则是这两个alpha还没有进行过区间化处理或者不在边界上。
应用简化版SMO算法处理小规模数据集
首先在数据集上遍历每一个alpha,然后在剩下的alpha集合中随机选择另一个alpha,从而构建alpha对。这里有一点非常重要,就是我们要同时改变两个alpha,之所以这样做是因为我们有一个约束条件:
由于改变一个alpha可能会导致改约束条件失效,因此我们总是同时改变两个alpha。
构建一个辅助函数,用于在某个区间范围内随机选择一个整数,同时,我们还需要另一个辅助函数,用于在数值太大时对其进行调整:
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):
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
运行验证:
dataArr,labelArr=loadDataSet('testSet.txt')
print(labelArr)
这类可以看出来,类别标签都是-1和1,而不是0和1.
下面是SMO第一个版本的伪代码:
创建一个alpha向量并将其初始化为0向量
当迭代次数小于最大迭代次数时(外循环):
对数据集中的每个数据向量(内循环):
如果该数据向量可以被优化:
随机选择另外一个数据向量
同时优化这两个向量
如果两个向量都不能被优化,退出内循环
如果所有向量都没被优化,增加迭代数目,继续下一次循环
下面是实现过程:
def smoSimple(dataMatIn,classLabels,C,toler,maxIter):
#输入参数分别为:数据集、类别标签、常数C、容错率、退出前最大的循环次数
#转换为mat矩阵,得到列向量
dataMatrix=mat(dataMatIn)
labelMat=mat(classLabels).transpose()
b=0
m,n=shape(dataMatrix)
#alpha列矩阵,初始化为0
alphas=mat(zeros((m,1)))
#遍历数据集的次数,当iter=maxIter时,函数结束运行并退出
iter=0
while (iter<maxIter):
#用来记录alpha是否已经进行优化。
alphaPairsChanges=0
for i in range(m):
fXi=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T))+b
#误差
Ei=fXi-float(labelMat[i])
#如果误差很大,可以对该数据实例对应的alpha值进行优化
#另外,如果alpha小于0或者大于C将被调整为0或C,一旦他们等于这两个值,就不值得强化了
if ((labelMat[i]*Ei<-toler) and (alphas[i]<C)) or ((labelMat[i]*Ei>toler) and (alphas[i]>0)):
#如果alpha可以更改,进入优化过程
j=selectJrand(i,m)
#随机选择第二个alpha
fXj=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T))+b
# 继续计算误差值
Ej=fXj-float(labelMat[j])
#通过copy()实现
alphaIold=alphas[i].copy()
alphaJold=alphas[j].copy()
#保证alpha在0与C之间
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为最优修改量,如果eta=0,需要退出循环的当前迭代过程。
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]=alphas[j]-labelMat[j]*(Ei-Ej)/eta
alphas[j]=clipAlpha(alphas[j],H,L)
if (abs(alphas[j]-alphaJold)<=0.00001):
print('J没有足够移动')
continue
#对i进行修改,修改量与j相同,但方向相反
alphas[i]=alphas[i]+labelMat[j]*labelMat[i]*(alphaJold-alphas[j])
#设置常数项
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
alphaPairsChanges=alphaPairsChanges+1
if alphaPairsChanges==0:
iter=iter+1
else:
iter=0
return b,alphas
实际运行:
dataArr,labelArr=loadDataSet('testSet.txt')
# print(labelArr)
b,alphas=smoSimple(dataArr,labelArr,0.6,0.001,40)
print(b)
print(alphas[alphas>0])
alphas[alphas>0]命令是数组过滤的一个实例,而且它只对NumPy类型有用,却并不适用于Python中的正则表,如果输入alpha>0,那么就会得到一个布尔数组,并且在不等式成立的情况下,其对应值是正确的。于是,在将该布尔数组应用到原始的矩阵当中时,就会得到一个NumPy矩阵,并且其中矩阵仅仅包含大于0的值。
为了得到支持向量的个数,输入:
print(shape(alphas[alphas>0]))
for i in range(100):
if alphas[i]>0.0:
print(dataArr[i],labelArr[i])
在原始数据集上对这些支持向量画圈后的结果如图: