SVM——Support Vector Machines
SVM１
SVM２
以上两个链接详细的证明了SVM的理论基础
每一步都很详细，不过数学基础差的同学可以直接跳过
毕竟。。。。。看不懂
但要想真正掌握这个算法，还是要沉下心来跟着他的推导一步步来
因为这个算法的实现就是整个推导过程的实现
拿出一整天来，跟着他的推导一步步跟
你的理解会更深刻
这是整个推导总结出来的最终步骤，接下来也是根据这个步骤来进行
在实现代码之前，需要几个辅助函数
from numpy import *
from time import sleep
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 #we want to select any J not equal to i
while (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
这几步很好理解，仔细看过推导的人都能很快明白这几步是干啥的
以下是简化版的SVM
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
#转换为numpy的mat存储
dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
#初始化b参数，统计dataMatrix的维度
b = 0; m,n = np.shape(dataMatrix)
#初始化alpha参数，设为0
alphas = np.mat(np.zeros((m,1)))
#初始化迭代次数
iter_num = 0
#最多迭代matIter次
while (iter_num < maxIter):
alphaPairsChanged = 0
for i in range(m):
#步骤1：计算误差Ei
fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
Ei = fXi - float(labelMat[i])
#优化alpha，更设定一定的容错率。
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
#随机选择另一个与alpha_i成对优化的alpha_j
j = selectJrand(i,m)
#步骤1：计算误差Ej
fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
Ej = fXj - float(labelMat[j])
#保存更新前的aplpha值，使用深拷贝
alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
#步骤2：计算上下界L和H
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
#步骤3：计算eta
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
#步骤4：更新alpha_j
alphas[j] -= labelMat[j]*(Ei - Ej)/eta
#步骤5：修剪alpha_j
alphas[j] = clipAlpha(alphas[j],H,L)
if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小"); continue
#步骤6：更新alpha_i
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
#步骤7：更新b_1和b_2
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
#步骤8：根据b_1和b_2更新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("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num,i,alphaPairsChanged))
#更新迭代次数
if (alphaPairsChanged == 0): iter_num += 1
else: iter_num = 0
print("迭代次数: %d" % iter_num)
return b,alphas
这是简化版的SVM，其中有几个continue可能是不好理解的部分，这也是推导中未曾提及的
１.Ｌ＝＝H，即上下界一样，不符合实际即continue重新开始循环
2.eta>=0,此处eta并不是推导过程中的η，观察他的表达式知道他其实是-η，所以此处条件即为η<=0,对η表达式研究过可知，他不可能为负数，因为xTx即为向量的长度，令x1的长度为a，x2的长度为b，则η为a方+b方-2abcosα，此处α为两向量夹角，此式子恒>=0,当他等于0的时候，即两向量同向，那么也就没有必要修改更新误差的必要了，所以continue重新开始下一轮循环
3.abs(alphas[i]-alphas[j])<0.00001,此时二者近似近似相等，即基本同向没有修改的必要了，所以continue
下面这两个公式想必已经不再陌生：
在实现SMO算法的时候，先计算η，再更新a_j。为了加快第二个α_j乘子的迭代速度，需要让直线的斜率增大，对于α_j的更新公式，其中η值没有什么文章可做，于是只能令：
因此，我们可以明确自己的优化方法了：
最外层循环，首先在样本中选择违反KKT条件的一个乘子作为最外层循环，然后用”启发式选择”选择另外一个乘子并进行这两个乘子的优化
在非边界乘子中寻找使得|E_i - E_j|最大的样本
如果没有找到，则从整个样本中随机选择一个样本
接下来，让我们看看完整版SMO算法如何实现。
首先需要再加入一些辅助函数
class optStruct:
def __init__(self,dataMatIn, classLabels, C, toler, kTup): # Initialize the structure with the parameters
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))) #first column is valid flag
self.K = mat(zeros((self.m,self.m)))
for i in range(self.m):
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
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 selectJ(i, oS, Ei): #this is the second choice -heurstic, and calcs Ej
maxK = -1; maxDeltaE = 0; Ej = 0
oS.eCache[i] = [1,Ei] #set valid #choose the alpha that gives the maximum delta E
validEcacheList = nonzero(oS.eCache[:,0].A)[0]
if (len(validEcacheList)) > 1:
for k in validEcacheList: #loop through valid Ecache values and find the one that maximizes delta E
if k == i: continue #don't calc for i, waste of time
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if (deltaE > maxDeltaE):
maxK = k; maxDeltaE = deltaE; Ej = Ek
return maxK, Ej
else: #in this case (first time around) we don't have any valid eCache values
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej
def updateEk(oS, k):#after any alpha has changed update the new value in the cache
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
此外仍需要刚才添加的辅助函数
def selectJrand(i,m):
j=i #we want to select any J not equal to i
while (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
以下是完整版Platt SMO核心代码
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) #this has been changed from selectJrand
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,i] - oS.K[j,j] #changed for kernel
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) #added this for the Ecache
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])#update i by the same amount as j
updateEk(oS, i) #added this for the Ecache #the update is in the oppostie direction
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,j]- 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
以上是内循环，就是只针对某一个i寻找最佳j与之配对，return的值只有0或1，当为1的时候就是已成功修改了alpha对的值，当为0的时候就是这个i有毒，修改失败。
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #full Platt SMO
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: #go over all
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:#go over non-bound (railed) alphas
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
if entireSet: entireSet = False #toggle entire set loop
elif (alphaPairsChanged == 0): entireSet = True
print "iteration number: %d" % iter
return oS.b,oS.alphas
以上是外循环框架，其中用到了内循环。
其中entireSet用的很奇妙，虽然不明白怎么想到的，怎么创造出来这么一个中间参数，但是仔细体会就会发现很奇妙，此外循环条件也发生了改变。
当迭代次数超过指定的最大值，或者遍历整个集合都未对任意alpha对进行修改时就退出循环。
此外还有就是
nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
这行代码很奇妙，他把非边界的样本筛选了出来，从而通过entireSet只需要对非边界点进行修改即可
解释是：由于非边界alpha数量很少(大部分alpha=0)，并且非边界alpha上更新更有效，原启发式对这一部分的处理合情合理完满
此次完整版的Platt SMO代码与原先的简化版的改进总结：
在实际运行中，发现原文中第二个启发式，即确定拉格朗日乘子alpha1的过程中，存在一些cpu算力浪费的现象：
1.在函数examineExample()中，对于使|E1-E2|最大的alpha1不能有效更新的情况，原启发式将会从随机起点处遍历所有的非边界alpha，逐一对其进行可能的更新，由于非边界alpha数量很少(大部分alpha=0)，并且非边界alpha上更新更有效，原启发式对这一部分的处理合情合理完满
2.但是，在遍历所有非边界alpha后，仍然不能有效更新的情况，或者非边界alpha数量小于2，原启发式将会从随机起点处遍历所有可能的alpha，逐一对其进行可能的更新。这种策略实际运行中效率十分低下，主要原因是遍历所有的alpha使之与alpha2形成一对拉格朗日乘子，每次计算的alpha2的更新值a2，都会非常靠近alpha2，更新量太小根本不值得更新，但是这个过程却要花费算力去计算a2，因此效率很低
3.一个合理的想法是，将从随机起点处遍历所有可能的alpha变换为：
在所有可能的alpha中随机选择一个作为alpha1，对乘子对(alpha1,alpha2)尝试更新only once，若更新不成功，就直接跳出examineExample()，以期待一个更优秀的alpha2
这样就可以将大量的算力用在真正大步长更新的乘子对上，会显著加快收敛速度
此外简单说一下核函数的技巧
核函数的意义在于把低维数据映射到高维平面里从而更易于找到超平面，这对我们处理非线性问题起到了巨大的作用
对于这种数据分布，我们很难在二维平面里面找到一条线把数据分开，所以我们用核函数把他映射到更高维度中从而可以分离数据
采用映射(x,y)->(x,y,x^2＋ y^2)就可以得到以下的样子
显然我们把二维数据变成了三维，然后瞬间就能想到一个超平面将他们分割开，从而达到区分数据的目的
这是核函数的大致思想，详情看本blog最上面的第二个链接，里面有更详细的解答。
最后鸣谢Jack-Cui大V的资源分享
扫一扫
2209
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
