标签:
SMO(Sequentil Minimal Optimization)算法在支持向量机中用来求解对偶问题,即
在这个问题中,变量是拉格朗日乘子
SMO算法是一个启发式的算法,它的基本思路是:如果所有变量的解都满足KKT条件,即:
由于每次只选择
最后可以得到
这里
因为还有限制条件,所以还要对这个解进行剪辑
然后依次更新
SMO称第一个变量的选择过程为外层循环,外层循环要选择一个违反KKT条件的变量,具体来说,若
又因为
也就是说要满足KKT条件要满足
同样的推导过程可以得到
若变量违反了KKT条件,我们就选择它为第一个变量,
在检验过程中,外层循环首先在所有变量中遍历,遇到违反KKT条件的变量就接着选择第二个变量,然后在整个集合上遍历一次后,
然后在所有非边界变量上遍历所谓非边界变量,就是指满足
的变量,为什么要这么遍历呢,因为随着多次子优化过程,边界变量倾向于留在边界,而非边界变量倾向于波动,这一步启发式的选择算法是基于节省时间考虑的,并且算法会一直在非边界变量集合上遍历,直到所有非边界变量都满足KKT条件(self-consistent)
随后算法继续在整个集合上遍历寻找违反KKT条件的变量作为优化的第一个变量,要注意的是,算法在整个集合上最多只连续遍历一次,但在非边界变量集合上可能连续遍历多次
选择第一个变量的代码如下:
def smo(oS, maxIter,):
iter = 0
entireSet = True
numChanged = 0
while (iter < maxIter) and ((numChanged > 0) or (entireSet)):
numChanged == 0
iter += 1
if entireSet:
for i in range(oS.m):
numChanged += innerL(oS, i)
else:
nonBounds = nonzero((oS.alphas.A > 0) * (oS.alphas.A < oS.C))[0]
for i in nonBounds:
numChanged += innerL(oS, i)
if entireSet:
entireSet = False
elif numChanged == 0:
entireSet = True
maxIter是最大迭代次数
上面说了第一个变量的选择方法,然后是第二个变量的选择方法
SMO称选择第二个变量的过程为内层循环,第二个变量的选择标准是希望
Platt原文说的是计算Kernel函数太耗时了,所以就以
至此两个要优化的变量就选择了出来
但是这里还有一个小问题,如果
具体方法是,首先在非边界元素上遍历,注意这里要随机选择位置,这是为了防止从头开始的话算法偏好靠前的元素,如果非边界元素上找不到使得
代码如下:
# 如果两个选择的两个输入向量的X相同,那么eta为0,此时要另外选择一个j
def secondChoiceJ(oS, i):
# 首先先从非边界元素中寻找是否eta不为零的j
nonBounds = nonzero((oS.alphas.A > 0) * (oS.alphas.A < oS.C))[0]
m = len(nonBounds)
st = int(random.uniform(0, m))
for i in range(m):
j = nonBounds[(st+i) % m]
if (j == i): continue
kIi = oS.K[i, i]
kIj = oS.K[i, j]
kJj = oS.K[j, j]
eta = kIi + kJj - 2*kIj
if (eta > 0):
return j
# 如果非边界找不到,那么到边界上找J
bounds = nonzero((oS.alphas.A == 0) + (oS.alphas.A == oS.C))[0]
m = len(bounds)
st = int(random.uniform(0, m))
for i in range(m):
j = bounds[(st+i) % m]
if (j == i): continue
kIi = oS.K[i, i]
kIj = oS.K[i, j]
kJj = oS.K[j, j]
eta = kIi + kJj - 2*kIj
if (eta > 0):
return j
return -1
预测函数为
如果是线性核函数的话,可以预先处理出
下面是Python代码,可以使用线性核或者高斯核函数
from numpy import *
# SVM算法
# 预处理数据集
def loadData(fileName):
dataSet = []
labels = []
fr = open(fileName)
for line in fr.readlines():
lineArr = line.strip().split(‘\t‘)
n = len(lineArr)
fltLine = list(map(float, lineArr[0:n-1]))
dataSet.append(fltLine)
labels.append(float(lineArr[-1]))
return dataSet, labels
# SMO算法
# 核函数
def kernelTrans(X, A, kTup):
m, n = X.shape
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] = deltaRow * (deltaRow.T)
K = exp(K/(-1*kTup[1]**2))
else: raise NameError(‘Unrecognizable Kernel‘)
return K
# 建立一个数据结构用来保存SMO算法中的变量
class optStruct:
def __init__(self, dataSet, labels, C, eps, kTup):
self.X = dataSet
self.Y = labels.T
self.C = C
self.eps = eps
self.m, self.n = shape(dataSet)
self.alphas = mat(zeros((self.m, 1)))
self.b = 0
self.eCache = mat(zeros((self.m, 2)))
self.K = mat(zeros((self.m, self.m)))
for i in range(self.m):
self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)
# 计算f(x_k)-y_k,E是函数f(x)对于输入x的预测值和真实值的偏差
def calcEk(oS, k):
# multiply表示两个矩阵对应元素之间相乘
fXk = multiply(oS.alphas, oS.Y).T*oS.K[:, k] + oS.b
return fXk - oS.Y[k]
# 更新E_k
def updateEk(oS, k):
Ek = calcEk(oS, k)
oS.eCache[k] = [1, Ek]
# 如果当前只有一个有效的E,则随机选择第二个alpha
def selectJrand(i, m):
j = i
while j == i:
j = int(random.uniform(0, m))
return j
# 内层启发式选择第二个alpha
def selectJ(oS, i, Ei):
maxK = -1
maxDelta = -1
Ej = 0
oS.eCache[i] = [1, Ei]
# .A返回矩阵的一个二维视图,不做任何拷贝
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(Ek-Ei)
if deltaE > maxDelta:
maxDelta = deltaE
maxK = k
Ej = Ek
return maxK, Ej
else:
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej
# 如果两个选择的两个输入向量的X相同,那么eta为0,此时要另外选择一个j
def secondChoiceJ(oS, i):
# 首先先从非边界元素中寻找是否eta不为零的j
nonBounds = nonzero((oS.alphas.A > 0) * (oS.alphas.A < oS.C))[0]
m = len(nonBounds)
st = int(random.uniform(0, m))
for i in range(m):
j = nonBounds[(st+i) % m]
if (j == i): continue
kIi = oS.K[i, i]
kIj = oS.K[i, j]
kJj = oS.K[j, j]
eta = kIi + kJj - 2*kIj
if (eta > 0):
return j
# 如果非边界找不到,那么到边界上找J
bounds = nonzero((oS.alphas.A == 0) + (oS.alphas.A == oS.C))[0]
m = len(bounds)
st = int(random.uniform(0, m))
for i in range(m):
j = bounds[(st+i) % m]
if (j == i): continue
kIi = oS.K[i, i]
kIj = oS.K[i, j]
kJj = oS.K[j, j]
eta = kIi + kJj - 2*kIj
if (eta > 0):
return j
return -1
# SMO算法中的内层循环
def innerL(oS, i):
Ei = calcEk(oS, i)
if ((oS.alphas[i] < oS.C) and (oS.Y[i]*Ei < -oS.eps)) or ((oS.alphas[i] > 0) and (oS.Y[i]*Ei > oS.eps)):
j, Ej = selectJ(oS, i, Ei)
alphaIold = oS.alphas[i].copy()
alphaJold = oS.alphas[j].copy()
s = oS.Y[i] * oS.Y[j]
if s < 0:
L = max(0, alphaJold-alphaIold)
H = min(oS.C, oS.C+alphaJold-alphaIold)
else:
L = max(0, alphaJold+alphaIold-oS.C)
H = min(oS.C, alphaJold+alphaIold)
if (L == H):
return 0
# eta
kIi = oS.K[i, i]
kIj = oS.K[i, j]
kJj = oS.K[j, j]
eta = kIi + kJj - 2*kIj
# 如果eta非正,那么需要重新选取一个j
if (eta <= 0):
return 0
j = secondChoiceJ(oS, i)
if j < 0: return 0
# 此时需要重新计算j的信息
alphaJold = oS.alphas[j].copy()
s = oS.Y[i] * oS.Y[j]
if s < 0:
L = max(0, alphaJold-alphaIold)
H = min(oS.C, oS.C+alphaJold-alphaIold)
else:
L = max(0, alphaJold+alphaIold-oS.C)
H = min(oS.C, alphaJold+alphaIold)
if (L == H):
return 0
kIj = oS.K[i, j]
kJj = oS.K[j, j]
eta = kIi + kJj - 2*kIj
aJ = alphaJold + oS.Y[j]*(Ei-Ej)/eta
## 剪辑alpha_j
if aJ > H: aJ = H
elif aJ < L: aJ = L
oS.alphas[j] = aJ
# 如果改进量不大,那么拒绝这次修改
if (abs(oS.alphas[j]-alphaJold) < oS.eps):
return 0
oS.alphas[i] = alphaIold + s*(alphaJold-oS.alphas[j])
bi = -Ei - oS.Y[i]*kIi*(oS.alphas[i]-alphaIold) - oS.Y[j]*kIj*(oS.alphas[j]-alphaJold) + oS.b
bj = -Ej - oS.Y[i]*kIj*(oS.alphas[i]-alphaIold) - oS.Y[j]*kJj*(oS.alphas[j]-alphaJold) + oS.b
if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
oS.b = bi
elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
oS.b = bj
else:
oS.b = (bi+bj)/2.0
updateEk(oS, i)
updateEk(oS, j)
return 1
return 0
# SMO算法中的外层循环
def smo(oS, maxIter, kTup=(‘lin‘, 0)):
iter = 0
entireSet = True
numChanged = 0
while (iter < maxIter) and ((numChanged > 0) or (entireSet)):
numChanged == 0
iter += 1
if entireSet:
for i in range(oS.m):
numChanged += innerL(oS, i)
else:
nonBounds = nonzero((oS.alphas.A > 0) * (oS.alphas.A < oS.C))[0]
for i in nonBounds:
numChanged += innerL(oS, i)
if entireSet:
entireSet = False
elif numChanged == 0:
entireSet = True
def testRbf(k1=1.3):
dataSet, labels = loadData(‘testSetRBF.txt‘)
dataSet = mat(dataSet)
labels = mat(labels)
oS = optStruct(dataSet, labels, 200, 0.0001, (‘rbf‘, k1))
smo(oS, 100, (‘rbf‘, k1))
# 构造支持向量矩阵
svInd = nonzero(oS.alphas.A > 0)[0]
sVs = oS.X[svInd]
labelSV = oS.Y[svInd]
print (‘there are %d support vectors‘ % len(svInd))
# 样本集出错概率
errorCount = 0
for i in range(oS.m):
kernelEval = kernelTrans(sVs, oS.X[i, :], (‘rbf‘, k1))
predict = kernelEval.T * multiply(labelSV, oS.alphas[svInd]) + oS.b
if sign(predict) != sign(oS.Y[i]):
errorCount += 1
print (‘Training error rate is: %f‘ % (float(errorCount)/oS.m))
# 新的测试集上的错误率
dataArr, labelArr = loadData(‘testSetRBF2.txt‘)
errorCount = 0
dataMat = mat(dataArr)
labelMat = mat(labelArr).T
m, n = shape(dataMat)
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], (‘rbf‘, k1))
predict = kernelEval.T * multiply(labelSV, oS.alphas[svInd]) + oS.b
if sign(predict) != sign(labelMat[i]):
errorCount += 1
print (‘test error rate is %f‘ % (float(errorCount)/m))
testRbf()
标签:
原文地址:http://blog.csdn.net/u014664226/article/details/51706921