完整版和简化版的smo不通之处在于alpha的选择方式上,alpha的更改和代数运算和简化版本,完整版的alpha的选择方式采用了启发式进行选择,前面文章已经详解了,这里再简单的叙述一下:
Platt SMO算法是通过一个外循环来选择第一个alpha值的,并且其选择过程会在两种方式之间进行交替即:一种方式是在所有数据集上进行单遍扫描,另一种方式则是在非边界alpha中实现单遍扫描,而所谓非边界alpha指的就是那些不等于边界0或者C的alpha值。对整个数据集扫描相当容易,而实现非边界alpha值的扫描时,首先需要建立这些alpha值的列表,然后在对这个表进行遍历,同时该步骤跳过那些已知不会改变的alpha值。
在选择第一个alpha值后,算法会通过一个内循环来选择第二个alpha值,在优化过程中,会通过最大步长的方式来获的第二个alpha值,在简易的smo中,我会在选择j后计算误差率E,但这里会建立一个全局的缓存用于保存误差值,并从中选择使得步长或者说Ei-Ej最大的alpha值,为什么是这样呢?看下式:
从上式可以清楚的看到的步长取决于和,而又无法人为改变,所以能改变的只有,使其差值越大,说明前进步长越大,下面给出代码,代码均已详细注释:
# 完整版smo算法实现# 下面的函数就是参数方面的存储定义# 简化版的smo都是在函数体内定义,这里单独拿出来,不同的是处理E值,增加了ecache项class optStruct:def __init__(self, dataMatIn, classLabels, C, toler):# dataMatIn、classLabels为输入数据和标签,C为松弛因子,toler为容忍度,简单来说就是alpha多大才认为不在变化了# kTup为核函数的信息self.X = dataMatInself.labelMat = classLabelsself.C = Cself.tol = tolerself.m = np.shape(dataMatIn)[0]self.alphas = np.mat(np.zeros((self.m, 1))) # 给alpha占位self.b = 0self.eCache = np.mat(np.zeros((self.m, 2)))# 首先这是一个矩阵,两列,第一列是表面E是否有效的标志位,第二列给出的是实际的E值# self.K = np.mat(np.zeros((self.m, self.m))) # 核函数方面的参数#for i in range(self.m):# self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)# 其实就是f(x) = ∑alpha_i*y_i*<x_i,x> + b,前面理论也详细讲了这一点,就是预测x的类别,只是单独拿出来了# 然后计算 差值 Ekdef calcEk(oS, k):fXk = float(np.multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.bEk = fXk - float(oS.labelMat[k])return Ek# 用来选择第二个alpha参数,或者说内循环的alpha值,其目标是选择合适的alpha值使其步长最大def selectJ(i, oS, Ei): # this is the second choice -heurstic, and calcs EjmaxK = -1 # 定义使的计算E的差值最大maxDeltaE = 0 # E值的最大变化量Ej = 0 # 使的差值最大的那个EjoS.eCache[i] = [1, Ei] # 这是把计算出来的差值记录并标为有效validEcacheList = np.nonzero(oS.eCache[:, 0].A)[0] # 刚开始有点不理解.A是什么意思,查了一下,解释一下# 情况下面的示例,其中.A的意思是把矩阵转换成数组,因为oS.eCache[:, 0].A,取的是一列的数据E值有效位,转换后还是一列# 而nonzero(oS.eCache[:, 0].A)[0]返回的是非零的目录列表值并取出赋给validEcacheList,共下面处理if (len(validEcacheList)) > 1:for k in validEcacheList: # 索引validEcacheList的标号,求最大差值if k == i: continue # 如果相等,说明选取两个alpha一样,直接跳出本次循环,进行下一次循环Ek = calcEk(oS, k) # 计算EkdeltaE = abs(Ei - Ek) # 求两个Ek的差值,选择最大的步长if (deltaE > maxDeltaE): # 比较,选择最大的差值maxK = k;maxDeltaE = deltaE;Ej = Ekreturn maxK, Ejelse:# 这个是为第一次准备的,不满足上面的条件即(len(validEcacheList)) > 1,因为第一次的E都为0,# 所以长度都为0,不满足条件,来到这里,进行随机抽取alpha 的 ijj = selectJrand(i, oS.m)Ej = calcEk(oS, j)return j, Ej'''In [12]: a = [[0,0,3],[0,5,0],[7,0,9]]In [13]: b = np.mat(a)In [14]: c = b.AIn [15]: print('type(a): ',type(a),'type(b): ',type(b),'type(c): ',type(c))type(a): <class 'list'> type(b): <class 'numpy.matrixlib.defmatrix.matrix'> type(c): <class 'numpy.ndarray'>In [16]: bOut[16]:matrix([[0, 0, 3],[0, 5, 0],[7, 0, 9]])In [17]: cOut[17]:array([[0, 0, 3],[0, 5, 0],[7, 0, 9]])In [18]: v =np.nonzero(c)[0]In [19]: vOut[19]: array([0, 1, 2, 2], dtype=int64)'''# 每一次循环计算的E都会缓存,即用于更新def updateEk(oS, k): # after any alpha has changed update the new value in the cacheEk = calcEk(oS, k)oS.eCache[k] = [1, Ek]# 完整版的SMO优化例程,和前面的简化版类似def innerL(i, oS):Ei = calcEk(oS, i) # 计算第一个alpha的差值# 寻找非边界的alpha,为寻找第二个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)):j,Ej = selectJ(i, oS, Ei) # 挑选第二个alpha,并得到Ej# 保存旧的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, 大家还记得上篇的参数更新里的: K_11 + K_22 - 2K_12吗,其实就是这个eta,有点不同,# 是因为假设条件不一样,但是原理相同eta = 2.0 * oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].Tif eta >= 0:print("eta>=0")return 0# 更新alphaoS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/etaoS.alphas[j] = clipAlpha(oS.alphas[j],H,L)# 把选出的第二个alpha对应的E添加到EcacheupdateEk(oS, j) #added this for the Ecache# 比较第二alpha的值变化是否在容忍度以内,如果在这属于不在变化的alpha了,则返回,反之继续执行if (abs(oS.alphas[j] - alphaJold) < 0.00001):print("j not moving enough")return 0# 更新另一个alpha,即第一个alphaoS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j# 同样添加到缓存区updateEk(oS, i) #added this for the Ecache#the update is in the oppostie direction# 计算b值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,:].Tb2 = 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# 根据条件更新bif (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2else: oS.b = (b1 + b2)/2.0return 1else:return 0# 外循环即第一个alpha的选择def smoPK(dataMatIn, classLabels, C, toler, maxIter): #full Platt SMO# 把数据传入数据结构中oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler)iter = 0entireSet = TruealphaPairsChanged = 0# 执行while的条件是1.迭代次数达到最大,2.alpha值存在变化,或者整个开始启动即entireSet = Truewhile (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):alphaPairsChanged = 0if entireSet: #go over allfor i in range(oS.m): # 遍历所有的alpha值alphaPairsChanged += innerL(i,oS) # 从上面我们知道进入innerL是第一个alpha选取就为i# 第二个跟剧最大的E差值进行选择也是在innerl完成选择# 如果寻找满足条件则返回1,反之返回0,所以alphaPairsChanged是记录变化的alpha的变化次数print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))iter += 1 # 整个遍历完以后alpha更新结束为一次迭代else:# 遍历非边界的alpha值,也就是不在0或者C上的值nonBoundIs = np.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 += 1if entireSet:entireSet = False # 开始时先遍历整个alpha,然后迭代完成一次以后,开始转向非边界进行遍历elif (alphaPairsChanged == 0):entireSet = True # 如果遍历完所有的非边界都没有可更新的alpha,则继续转为全局遍历print("iteration number: %d" % iter)return oS.b,oS.alphas# 计算w,大家还记的求w的式子吗?# W = ∑alpha_i*y_i*x_idef calcWs(alphas, dataArr, classLabels):X = np.mat(dataArr);labelMat = np.mat(classLabels).transpose()m, n = np.shape(X)w = np.zeros((n, 1))for i in range(m):w += np.multiply(alphas[i] * labelMat[i], X[i, :].T)# 虽然该过程遍历所有的样本点,但是通过前面的原理详解可知,大部分的alpha值为0,只有支持向量的值不为零# 所以这个计算还是很简单的return w
下面给出测试代码:
b, alphas = SMO_simplify.smoPK(dataarr, labelarr, 0.6, 0.001, 40)print('b: ',b)print('alphas: ',alphas[alphas>0])ws = SMO_simplify.calcWs(alphas, dataarr, labelarr)print('ws: ',ws)# 对数据进行预测datmat = np.mat(dataarr)print('datmat[0]*np.mat(ws) + b = ', datmat[0]*np.mat(ws) + b)print('labelarr[0]',labelarr[0])err=0# 计算全体数据的错误率for i in range(len(dataarr)):y = datmat[i]*np.mat(ws) + bif y > 0:y = 1else:y = -1if y != labelarr[i]:err += 1;print('err: ', err)print('错误率为: ', err/len(dataarr))
测试 结果:
non-bound, iter: 1 i:46, pairs changed 0non-bound, iter: 1 i:55, pairs changed 0non-bound, iter: 1 i:94, pairs changed 0iteration number: 2j not moving enoughfullSet, iter: 2 i:0, pairs changed 0fullSet, iter: 2 i:1, pairs changed 0fullSet, iter: 2 i:2, pairs changed 0j not moving enoughb: [[-2.89901748]]alphas: [[0.06961952 0.0169055 0.0169055 0.0272699 0.04522972 0.02726990.0243898 0.06140181 0.06140181]]ws: [[ 0.65307162][-0.17196128]]datmat[0]*np.mat(ws) + b = [[-0.92555695]]labelarr[0] -1.0err: 0错误率为: 0.0
简单测试100%正确,这是线性分类,如果是非线性呢?下一节将详细讨论核函数实战代码