Python学习-机器学习实战-ch06 支持向量机

支持向量机是我学习机器学习时期最重要的分类算法（没有之一），但是每次都似懂若无的感觉。其中用到了间隔最大化、对偶问题求解等等步骤。其中，支持向量机的学习问题可以形式化为求解凸二次规划问题，有许多最优化算法可以用于这个问题求解。当样本容量很大时，这些方法往往变得非常低效。本章还涉及支持向量机的最优化求解方法。

==================================================================

序列最小最优算法（sequential minimal optimization,SMO）是一种启发式方法，其基本思路是：如果所有变量的解都满足此最优化问题的KKT条件，那么这个最优化问题的解就得到了。因为KKT条件是该优化问题的充分必要条件。否则，选择两个变量，固定其他变量，针对这两个变量构建一个二次规划问题，这个二次规划问题关于这两个变量的解应该更接近原始二次规划问题的解，因为这会使得原始二次规划问题的目标函数值变得更小。

整个SMO算法包括两个部分：求解两个变量二次规划的解析方法和选择变量的启发式方法。

——李航《统计学习方法》

1.简化版的SMO

原版SMO算法需要在外循环确定要优化的alpha对。简化版的SMO首先在数据集上遍历每一个alpha，然后在剩下的alpha中随机选择另一个alpha，构建alpha对。

+++++++++++++++++++++++++++++++++++++++++

简化版的SMO的伪代码如下：

创建一个alpha向量并将其初始化为0向量

当迭代次数下去最大迭代次数时（外循环）：

对数据集中的每个数据向量（内循环）：

如果该数据向量可以被优化：

随机选择另外一个数据向量

同时优化这两个向量

如果这两个向量都不能被优化，退出内循环

如果所有向量都没被优化，增加迭代数目，继续下一次循环

+++++++++++++++++++++++++++++++++

辅助函数：

from numpy import *

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的总个数
j=i
while(j==i):
#j是第二个alpha的下标
#如果j与i相同，则重新选取一个
j=int(random.uniform(0,m))
return j

def clipAlpha(aj,H,L):
#设置上下界限
if aj>H:
aj=H
if L>aj:
aj=L
return aj

简化版的SMO算法：

<pre name="code" class="python">def smoSimple(dataMatIn,classlabels,C,toler,maxIter):
#函数输入：数据、标签集、常数、容错率、最大循环次数
dataMatrix=mat(dataMatIn)
labelMat=mat(classlabels).transpose()
b=0
m,n=shape(dataMatrix)
#m表示样本个数，n表示特征维度
alphas=mat(zeros((m,1)))
iter=0
while(iter<maxIter):
#外循环
alphaPairsChanged=0
for i in range(m):
fXi=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T))+b
#multiply是numpy的乘法函数，.T是转置
#第i个样本对应的判别结果
#(dataMatrix*dataMatrix[i,:].T)是一个核函数计算
Ei=fXi-float(labelMat[i])
#误差，预测结果与真实类别值
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
#计算第二个alpha的值
Ej=fXj-float(labelMat[j])
#得到两个样本对应的两个alpha对应的误差值
alphaIold=alphas[i].copy()
alphaJold=alphas[j].copy()
#存储原本的alpha值
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
#计算上下阈值
#针对y1,y2的值相同与否，上下值也不同
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
#更新alpha
alphas[j]=clipAlpha(alphas[j],H,L)
if (abs(alphas[j]-alphaJold)<0.00001):
print('j is not moving enough')
continue
alphas[i]+=labelMat[j]*labelMat[i]*(alphaJold-alphas[j])
#由一个alpha确定另一个alpha
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-Ei-labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T-labelMat[j]*(alphas[j]-alphaJold)\
*dataMatrix[j,:]*dataMatrix[j,:].T
#更新两个b值
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

补充说明：

eta的公式是：



alpha2的未处理值为：



alpha的截断处理：



由alpha2得出alpha1：



子问题有两个变量，一个是违反KKT条件最严重的哪那一个，另一个由约束条件自动确定。

SMO称选择第一个变量的过程为外层循环。外循环在训练样本中选取违反KKT条件最严重的样本点，并将其对应的变量作为第一个变量。具体地，检验训练样本点是否满足KKT条件，即



外层循环首先遍历所有满足0到C范围内的alpha，即在间隔边界上的支持向量点，检验它们是否满足KKT条件。如果这些样本点都满足KKT条件，那么遍历整个训练集，检验它们是否满足KKT条件。

第2个遍历作为内层循环。假设子外层循环已经找到第1个遍历，则在内层循环找第2个变量。第2个变量的标准是使第2个变量有足够大的变化。alpha2是依赖于E1-E2的。为了加快速度，一种简单的方法是选择alpha2，使E1-E2最大。这是因为alpha1确定的情况下，E1也就确定了。如果E1是正的，选择最小的Ei作为E2。如果E1是负的，则选择最大的Ei作为E2.【《统计学习》】



【《机器学习实战》】

通过一个外循环来选择第一个alpha值，并且在其选择过程中会在一下两种方式之间交替进行：

1、在所有数据集上进行单遍扫描

2、在非边界alpha中实现单遍扫描，非边界alpha是指那些不等于边界0或C的alpha值。

选择完第一个alpha值有后，算法通过一个内循环来选择第二个alpha值，通过最大化步长的方式获得。

class optStruct:
def __init__(self,dataMatIn,classLabels,C,toler):
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)))
#用来缓存误差
#是两列的，第一列表示是否有效，第二列是实际的E值

def calcEk(oS,k):
#用来计算误差的函数
fXk=float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T))+oS.b
Ek=fXk-float(oS.labelMat[k])
return Ek

def selectJ(i,oS,Ei):
#在确定好第一个alpha的情况下，确定第二个
#求找最大步长，E1-E2
maxK=-1
maxDeltaE=0
Ej=0
oS.eCache[i]=[1,Ei]
#将Ei设置为有效
validEcacheList=nonzero(oS.eCache[:,0].A)[0]
#nonzero返回一个列表，这个列表中包含以输入列表为目录的列标识
#返回非零E值所对应的alpha值
#因为在eCache的第一列代表是否有效，非0代表有效
if(len(validEcacheList))>1:
for k in validEcacheList:
if k==i:continue
Ek=calcEk(oS,k)
deltaE=abs(Ei-Ek)
if (deltaE>maxDeltaE):
maxK=k
maxDeltaE=deltaE
Ej=Ek
return maxK,Ej
else:
#如果都不满足要求，直接随机选一个
j=selectJrand(i,oS.m)
Ej=calcEk(oS,j)
return j,Ej

def updateEk(oS,k):
Ek=oS.calcEk(oS,k)
oS.eCache[k]=[1,Ek]

建立一个数据结构，以及辅助函数。在之前简易SMO的基础上进行修改。主要是在选择第二个alpha的函数有了更大修改。

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)
#启发式方法选择第二个alpha
alphaIold=oS.alphas[i].copy()
alphaJold=oS.alphas[j].copy()
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[j,:].T-oS.X[j,:]*oS.X[j,:].T
#计算eta
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)
#更新对应的误差值
if(abs(oS.alphas[j]-alphaJold)<0.00001):
print('j is not moving enough')
return 0
oS.alphas[i]+=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j])
#由第二个alpha求解另一个alpha
updateEk(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]-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
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

innerL函数是一个更新版，可以与smoSimple函数进行比较。主要有两点不同：1、使用了数据结构了对参数进行传递。2、修改了选择第二alpha的方式。

接下来是外循环函数代码。

def smoP(dataMatIn,classLabels,C,toler,maxIter,kTup=('lin',0)):
oS=optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
#初始化数据结构
iter=0
entireSet=True
alphaPairChanged=0
while(iter<maxIter) and ((alphaPairChanged>0)or(entireSet)):
#循环条件：1、迭代次数少于最大迭代数；2、遍历着数据集对alpha进行了改变
alphaPairChanged=0
if entireSet:
for i in range(oS.m):
#oS.m表示数据的个数
alphaPairChanged+=innerL(i,oS)
#此处i对数据集进行遍历，InnerL选择第二个alpha，如果有改变返回1，否则返回0
print('fullSet,iter: %d i: %d, Pairs changed: %d' %(iter,i,alphaPairChanged))
iter+=1
else:
nonBoundIs=nonzero((oS.alphas.A>0)*(oS.alphas.A<C))[0]
#找所有非零的非边界值
for i in nonBoundIs:
alphaPairChanged+=innerL(i,oS)
print('non-bound,iter: %d i: %d, Pairs changed: %d' %(iter,i,alphaPairChanged))
iter+=1
if entireSet:entireSet=False
elif (alphaPairChanged==0):entireSet=True
print('iteration number: %d' %iter)
return oS.b,oS.alphas

此处entireSet相当于一个开关，要思考什么时候外循环会停止。

两种条件停止：
1、迭代次数大于maxIter

2、(alphaPairChanged>0)or(entireSet))的结果为0。即要嘛没有alpha对发生改变，要嘛entireSet等于false。

entireSet==False的情况按照if语句来看也有两种，一是对所有的数据执行完遍历。（这里没懂）

if entireSet:entireSet=False

还有就是对所有的非边界值进行遍历，发生了alpha修改。

def calcWs(alphas,dataArr,classLabel):
X=mat(dataArr)
labelMat=mat(classLabel).transpose()
#获得数据和类别标签
m,n=shape(X)
#m表示数据个数，n表示数据维度
w=zeros((n,1))
#初始化权重向量
for i in range(m):
w+=multiply(alphas[i]*labelMat[i],X[i,:].T)
return w

通过alpha求解权重w，就可以得到分类器了。



2、核技巧

对于更为复杂的非线性情况，上述方法是没法解决的。可以通过数据的特征空间变换，从低维度空间映射到高维度空间，使得数据更容易被分类。

书中介绍了一种常用的核函数——径向基函数。径向基函数是一个采用向量作为自变量的函数，能够基于向量距离运算输出一个标量。

径向基函数的高斯版本：



上式中的Y容易造成误导，并不是表示类别标签，而是其中的x,y表示两个向量。

def kernelTrans(X,A,kTup):
#kTup是一个包含核函数信息的元祖
m,n=shape(X)
K=mat(zeros((m,1)))
if kTup[0]=='lin':K=X*A.T
#最简单的内积方式
elif kTup[0]=='rbf':
for j in range(m):
deltaRow=X[j,:]-A
K[j]=deltaRow*deltaRow.T
K=exp(K/(-1*kTup[1]**2))
else:
raise NameError('Houston We have a Problem -- That Kernel is not recognized')
return K

这是一个核函数转换的函数。

kTup包含核函数的信息，第一个参数表示核函数类型，第二个是核函数的参数。上面的函数提供了两种核函数方法。

1、内积，计算所有数据集和数据集中的一行之间的展开。

2、for循环对每个元素计算高斯函数的值，for循环后，将计算过程应用的整个向量中去。

deltaRow表示的是|x-y|

class optStruct:
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)

相比原本的数据结构，增加了核函数的信息。

K是先构造矩阵，然后通过循环将数据进行转换填充进去。

然后将InnerL和calcEk函数里头的内积替换成核函数转换后的结果。

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)
#启发式方法选择第二个alpha
alphaIold=oS.alphas[i].copy()
alphaJold=oS.alphas[j].copy()
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.K[i,j]-oS.K[i,j].T-oS.K[j,j]
#计算eta
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)
#更新对应的误差值
if(abs(oS.alphas[j]-alphaJold)<0.00001):
print('j is not moving enough')
return 0
oS.alphas[i]+=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j])
#由第二个alpha求解另一个alpha
updateEk(oS,i)
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,i]\
-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

def calcEk(oS,k):
#用来计算误差的函数
fXk=float(multiply(oS.alphas,oS.labelMat).T*(oS.K[:,k]))+oS.b
Ek=fXk-float(oS.labelMat[k])
return Ek

编写一个测试用的函数，如下：

def testRbf(k1=1.3):
dataArr,labelArr=loadDataSet('testSetRBF.txt')
b,alphas=smoP(dataArr,labelArr,200,0.0001,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)
errCount=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]):
errCount+=1
print('the training error rate is %f' %(float(errCount)/m))
dataArr,labelArr=loadDataSet('testSetRBF2.txt')
errCount=0
datMat=mat(dataArr)
labelMat=mat(labelArr).transpose()
m,n=shape(dataArr)
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]):
errCount+=1
print('the training error rate is %f' %(float(errCount)/m))

其实，还有一个地方要修改，就是smoP函数的optStruct构造时需要加入kTup参数，同时将函数输入的内积默认给去掉。

def loadImages(dirName):
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])
if classNumStr==9:hwLabels.append(-1)
else:hwLabels.append(1)
trainingMat[i,:]=img2vector('%s/%s' % (dirName,fileNameStr))
return trainingMat,hwLabels

def testDigits(kTup=('rbf',10)):
dataArr,labelArr=loadImages('trainingDigits')
b,alphas=smoP(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)
errCount=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]):
errCount+=1
print('the test error rate is %f' %(float(errCount)/m))
dataArr,labelArr=loadImages('testDigits')
errCount=0
datMat=mat(dataArr)
labelMat=mat(labelArr).transpose()
m,n=shape(dataArr)
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]):
errCount+=1
print('the test error rate is %f' %(float(errCount)/m))

此过程中将数字9作为正类，其余数字作为反类。

该测试函数与之前无异，只是多了将图像转换成数据的时间，运算时间较久。