经典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)