统计学习第七章-SMO算法

SOM基本原理

基本约束为：

KKT条件为：

得到更新公式为

α由于受到边界条件约束，则可以得到如下的结论：


以上我们可以计算α2了，通过α2我们可以计算α1

b的更新公式也可以得到：

详细情况见SMO算法

SOM完整算法

有了完整的更新公式和算法，我们可以很容易的写出程序。

#找一个结构体，用于存储数据
class optStruct:
	def _init_(self, 数据集, 类别标签, 常数C, 容错率, 最大循环次数):
		self.X = 数据集
		self.lableMat = 类别标签
		self.C = 常数C
		self.tol = 容错率
		#返回数据集的行数目
		self.m = 数据集.shape[0]
		#α为m行1列的矩阵
		self.alphas = mat(zeros((self.m, 1)))
		#初始化b=0
		self.b = 0
		#误差为m行2列的矩阵，一行用来判定是不是为0，一行存新的Ei
		self.eCache = mat(zeros((self.m, 2)))

	def selectrand(i,m)
		#用于随机选择一个第二α，在无法比较的时候
		#这里首先赋值，是j=i，那么下面的条件必然生效，选择出一个不等于i的值，如果直接用选择，可能会选择出i
		j = i
		while(j == i):
			j = int(random.uniform(0,m))
		return j


	def calcEk(oS, k):
		#用于计算Ek = g(x)-yi
		#对应函数位置相乘multipy()
		Fxk = float(multipy(oS.alphas,oS.lableMat).T*(oS.X*oS.X[k,:].T)+oS.b)
		Ek = Fxk- float(oS.lableMat[k])
		#选择第二个循环的α


	def selectJ(i, oS, Ei):
	#初始化最大k，初始化最大间隔maxDeltaE，最大误差0
		maxK = -1; maxDeltaE = 0; Ej =0
		#将数据保存
		oS.eCache[i] = [1,Ei]
		#.A为转换矩阵数据类型，
		#α的更新函数严格依赖于|E1-E2|，提取出那些E不为0的，如果为0则保持不变，说明有很多解都是等效的，进行下一步比较
		validEcacheList = nonzero(oS.eCache[:，0].A)[0]
		#如果只有一个不为0，说明剩下的都为0，则跳转到else，随机选择一个
		if(len(validEcacheList)) > 1:
			#对不为0的值，取一个遍历
			for k in validEcacheList:
			#如果这个点是本身就跳过
				if k == i:
					continue
				#计算对每个点的Ek，计算对这个点的最大间隔
				Ek = calcEk(oS, k)
				deltaE = abs(Ei-Ek)
				#最大间隔肯定比0大，所以必然会有一个点返回，但实际上可能比0的效果差，这里就暂时算作还可以吧
				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


	#函数跟新缓存并放置到eCache，第一行为1，表示存在这个数
	def updata(oS, k):
		Ek = calcEk(oS, k)
		oS.eCache[k] = [1,Ek]


	#简单的选择函数
	def clipAlpha(aj,H,L):
		if aj > H:
			aj = H
		if aj < L:
			aj = L
		return aj


#函数的作用是输入αi选择一个最优的αj，然后进行一次更新
def	innerL(i, oS)
	#计算Ek,为Ej选择做准备
	Ei = calcEk(oS, i)
	#判断是否在精度范围可优化
	if ((oS.lableMat[i]*Ei < -oS.tol)and(oS.alphas[i<os.C)or(oS.lableMat[i]*Ei > oS.tol)and(oS.alphas[i]>0)):
		j,Ej = selectJ(i, oS, Ei)
		#将i，j进行复制为i_old，j_old
		alphaIold = oS.alphas[i].copy
		alphaJold = oS.alphas[j].copy
		#y1不等于y2计算边界,计算H，L
		if (oS.lableMat[i] != oS.lableMat[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])	
		#如果相等，则说明解为0，返回解		
		if L == H:
			print("L==H")
			return 0
		#计算λ，
		eta = 2*oS.K[i,:]*oS.K[j,:].T-oS.K[i,:]*oS.K[i,:].T-oS.K[j,:]*oS.K[j,:].T
		#判断是否大于0 ，我觉得不可能大于0.等于0是可能的
		if eta >= 0:
			print("eta = 0")
			return 0
		#判断边界条件，更新αj
		oS.alphas[j] = oS.alphas[j]-oS.lableMat[j]*(Ei-Ej)/eta
		oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
		update(oS, j)
		#判断αj的最佳更新边界
		if (abs(oS.alphas[j]-alphaJold) < 0.0001):
			print("行吧，不更新了吧")
			return 0
		#判断边界条件，更新αi
		oS.alphas[i] = oS.lableMat[i]*oS.lableMat[j]*(oS.alphas[j]-alphaJold)
		update(oS, i)
		#更新b1 ，b2
		b1 = oS.b - Ei - oS.lableMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,:]*oS.K[i,:].T-oS.lableMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,:]*oS.K[i,:].T
		b2 = oS.b - Ej - oS.lableMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,:]*oS.K[j,:].T-oS.lableMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,:]*oS.K[j,:].T
		#对b的取值做判断
		if (0 < oS.alphas[i])and(oS.alphas[i] < oS.C)
			oS.b = b1
		elif (0 < oS.alphas[j])and(oS.alphas[j] < oS.C)
			oS.b = b2
		else
			oS.b = (b2+b1)/2
		return 1
	else
		return 0

def SMO(数据集, 类别标签, 常数C, 容错率, 最大循环次数, kTup = ('lin', 0)):
#初始化oS数据，iter产生迭代器
	oS = optStruct(mat(数据集), mat(类别标签), 常数C, 容错率, 最大循环次数)
	iter = 0
	#entireSet用于转换判断
	entireSet = True
	#记录每个α的跟新次数
	alphaPairsChanged = 0
	#如果α不变了，优化达到最大次数了，则停止优划
	while((iter < 最大循环次数) and (alphaPairsChanged > 0) or entireSet):
		#每次遍历重新统计优化次数
		alphaPairsChanged = 0
		#entireSet为真，则进行内部条件的遍历，先对所有的α进行优划
		if entireSet:
		#对M个α进行遍历，然后进行优划
			for i in range(oS.m):
				alphaPairsChanged = alphaPairsChanged + innerL(i, oS)
				print("迭代器：%d，αi：%d，迭代次数：%d"%(iter,i,alphaPairsChanged))
				iter = iter +1
		else:#对支持向量α进行优划，去掉那些不满足条件的α
			nonBoundIs = nonzero((oS.alphas.A>0)*(oS.alphas.A<C))[0]
			for i in nonBoundIs:
				alphaPairsChanged = alphaPairsChanged + innerL(i, oS)
				print("迭代器：%d，αi：%d，迭代次数：%d"%(iter,i,alphaPairsChanged))
				iter = iter +1
		#对所有的α遍历完了之后，将entireSet变为False，只对满足条件的点进行优化
		if entireSet:
			entireSet = False
		## 如果内部点不能再优化了，重新对所有点进行优化
		elif alphaPairsChanged == 0:
			entireSet = True
	#对所有点都优化不动了，则alphaPairsChanged = 0，当 entireSet = False的时候，跳出循环
	#当达到最大次数时，当 entireSet = False的时候，跳出循环
	print('算完了')
	return oS.b, os.alphas

对于核函数，采用映射的方式，将上文中的核内积替换就行，或者在根本上取更新K。