经典SVM之SMO算法实现

最新推荐文章于 2021-12-10 20:28:39 发布

YouthDance

最新推荐文章于 2021-12-10 20:28:39 发布

阅读量1w

点赞数 7

分类专栏： _____MachineLearning_____

本文链接：https://blog.csdn.net/zhongshijunacm/article/details/52006552

版权

_____MachineLearning_____ 专栏收录该内容

11 篇文章 0 订阅

订阅专栏

经典SVM之SMO算法实现

一、浅谈理论

（1）原始问题部分

对于理论不做过深的解释和讨论。原因有两个。第一，我也不会。第二，写公式符号太头疼!所以，只是简单给出一些最重要的公式或者程序直接使用的公式并做适当的解释。现在就给出SVM最核心最重要的一个公式：

线性支持向量机原始最优化问题：

非线性支持向量机学习算法：

【1】选取适当的核函数K(x,z)和适当的参数C，构造并求解最优化问题

（1.1）

（1.2）

求得最优解

（1.3）

【2】选择

的一个正分量0 < αj < C,计算

（1.4）

【3】构造决策函数：

（1.5）

当K(x,z)是正定核函数时，（1.1）~（1.5）是凸二次规划问题，解是存在的。

（2）问题求解部分

对于凸二次规划问题的求解有很多方法，比较流行的一种方法是用SMO(sequence minimal optimization)算法。

【1】SMO算法介绍

SMO算法是一种启发式算法，其基本思路是：如果所有变量的解都满足此最优化问题的KKT条件（1.1~1.5），那么这个最优化问题的解就得到了。因为KKT条件是该最优化问题的充分必要条件。否则选择两个变量，固定其他变量针对这两个变量构建一个二次规划问题。

【2】两个变量二次规划的求解方法

处理边界L，H的值：

时，

（1.6）

时，

（1.7）

求解决策误差Ei:

（1.8）

当i=1,2时，Ei为函数g(x)对输入xi的预测值与真实值yi之差。

对于求解未剪辑前的第二个变量的公式：

（1.9）

经剪辑后的

的解是：

（2.0）

由

求得

是

（2.1）

【3】计算阀值b

（2.2）

【3】总结分析

（1）第1个变量的选择

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

sigmai = 0 <==>yi * g(xi) >= 1

0 < sigmai < C <==> yi * g(xi) = 1

sigmai = C <==> yi * g(xi) <= 1

其中，g(xi) = ∑[j=1,..,N] sigmaj * yj * K(xi, xj) + b

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

（2）、第2个变量的选择

SMO称选择第2个变量的过程为内层循环。假设在外层循环中已经找到了第1个变量sigma1，现在要在内层循环中找到第2个变量sigma2.第2个变量选择的标准是希望能使sigma2有最够大的变化。

sigma2[new] 是依赖|E1 - E2|的，为了加快计算速度，一种简单地做法是选择sigma2，使其对应的|E1 - E2|最大。

还有一个公式，在线性核时候会用到。计算w*

（2.3）

二、代码实现

（1）数据结构

class optStruct:
	#@dataMatIn: 数据集,type: mat
	#@classLabels: 类标签,type: mat
	#@C：自设调节参数
	#@toler: 自设容错大小
	#@kTup: 核函数类型
	def __init__(self, dataMatIn, classLabels, C, toler, kTup):
		self.X = dataMatIn
		self.labelMat = classLabels
		self.C = C
		self.toler = toler
		self.m = shape(dataMatIn)[0]
		self.alphas = mat(zeros((self.m, 1)))   #初始化一个m的列向量α
		self.b = 0
		self.eCache = mat(zeros((self.m, 2)))	#误差(Ei)缓存
		self.K = mat(zeros((self.m, self.m)))   #初始化一个存储核函数值得m*m维的K
		for i in range(self.m):					#获得核函数的值K(X,xi)
			self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)

（2）核转换函数（高斯核）

#@kTup: lin: 线性 
#       rbf: 高斯核 公式：exp{(xj - xi)^2 / (2 * δ^2)} | j = 1,...,N
#       δ：有用户自设给出kTup[1]   
def kernelTrans(X, A, 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] = exp(K / (-1 * kTup[1] ** 2))
	else: raise NameError('Houston we have a promble -- That kernel is not recoginzed')
	return K

（3）计算误差值Ei

#根据公式【4】，计算出Ei
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

（4）根据（1.9）选择第二个变量α

#随机选择第一个α
def selectJrand(i, m):
	j = i
	while (j == i):
		j = int(random.uniform(0, m))
	return j
	
#根据最大步长选择第二个变量	
def selectJ(i, oS, Ei):
	maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1, Ei]        		 #Ei保存的两维中，第一维表示Ei的有效性，第二维才是真实值
	validEcachelist = nonzero(oS.eCache[:, 0].A)[0]
	if (len(validEcachelist)) > 1:
		for k in validEcachelist:
			if k == i: continue
			Ek = calcEk(oS, k)
			deltaE = abs(Ei - Ek)
			if (deltaE > maxDeltaE):      #选择具有最大步长(max|Ei - Ej|)的j,
				maxK = k; maxDeltaE = deltaE; Ej = Ek
		return maxK, Ej
	else:                                 #随机选择第一个α
		j = selectJrand(i, oS.m)
		Ej = calcEk(oS, j)
	    return j, Ej

（5）跟新Ei值

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

（6）Platt SMO算法中的优化例程(选择第一个变量α)

#调整大于H或小于L的α值
def clipAlpha(aj, H, L):
	if aj > H:
		aj = H
	if aj < L:
	    aj = L
	return aj
	
#选择第一个变量α
def innerL(i, oS):
	Ei = calcEk(oS, i)
	if ((oS.labelMat[i] * Ei < -oS.toler) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i] * Ei > oS.toler) and (oS.alphas[i] > 0)):
		j, Ej = selectJ(i, oS, Ei)      									 #第二个α选择的启发式方法
		alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy()  	 #进行深复制处理
		
		#根据公式[6]
		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
		
		#根据公式[5]
		eta = 2.0 *  oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
		if eta >= 0; print "eta >= 0"; return 0                             #??为什么eta>=0就可表示剪辑过？
		oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
		oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
		updatEk(oS, j)														#跟新误差缓存
		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])
		updatEk(oS, i)
		
		#根据公式[3]
		b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, i] - oS.labelMat[j] * (oS.labelMat[j] - alphaJold) * oS.K[i, j]
		b2 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, j] - oS.labelMat[j] * (oS.labelMat[j] - alphaJold) * oS.K[j, j]
	    
		#根据SMO中优化的部分选择b的公式
		if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
		elif (o < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
		else: oS.b = (b1 + b2) / 2.0
		return 1
	else: return 0

（7）Platt SMO的外循环代码

#@maxIter:程序最大的循环次数
#其他变量的含义已经在数据结构中介绍了，这里就不在赘诉
	
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin', ):
	oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
	iter = 0
	entireSet = True; alphaPairsChanged = 0
	while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
		alphaPairsChanged = 0
		if entireSet:															#遍历所有值
			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:																	#遍历非边界值（任意一对α值发生改变，那么就会返回1）
			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
		iter += 1
		if entireSet: entireSet = False
		elif (alphaPairsChanged == 0): entireSet = True
		print "iteration number: %d" % iter
	return oS.b, oS.alphas

（7）计算w值

#@param: alphas, dataAr, classLabels --> Type: List		
def calWs(alphas, dataArr, classLabels):
	X = mat(dataArr); labelMat = mat(classLabels).transpose()
	m, n = shape(X)
	w = zeros((n, 1))
	for i in range(m):
		w += multiply(alphas[i] * labelMat[i], X[i, :].T)
	return w

（8）根据决策函数验证结果

def testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('testSetRBF.txt')
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd] #get matrix of only support vectors
    labelSV = labelMat[svInd];
    print "there are %d Support Vectors" % shape(sVs)[0]
    m,n = shape(datMat)
    errorCount = 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]): errorCount += 1
    print "the training error rate is: %f" % (float(errorCount)/m)
    dataArr,labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    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]): errorCount += 1    
    print "the test error rate is: %f" % (float(errorCount)/m)