本文主要是介绍SMO算法。对偶问题请转至上篇。
SVM的目的就是通过二次规划算法求解下式的最值。
但是，该问题的规模正比于训练样本数。数据样本越大，计算成本越大。因此，提出高效算法SMO。
SMO(Sequential Minimal Optimization)算法是提供了解决二次规划问题的优化算法。
有关SMO的历史，自行百度。
本文重点参考周志华老师的西瓜书和《Machine Learning In Action》
个人的理解是，如果单单看任何一本书可能都会有点迷惑，最好是相互参照，这样一来一往不仅加深印象，同时加深了理解。
"Talk is cheap, show me the code."
直接上代码吧。
#数据导入模块
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): #生成随机数 m为最大的值
j=i
while (j==i): #避免生成的随机数与导入的数相同
j = int(random.uniform(0,m))
return j
def clipAlpha(aj,H,L): #固定ai之外的其他变量，smo每次选择两个变量，对其范围进行约束
if aj >H:
aj =H
if aj <L:
aj =L
return aj
在alpha的选择上，我们为什么要选择两个alpha，alpha[i] , alpha[j]？
根据我们得到的约束条件，
Assume 我们提出alpha1
我们得到，
如果我们固定除了alpha1以外的参数，实际上，alpha1也变成了常数。
因此，我们需要选取两个alpha。
到这里我们发现固定除了alpha1、alpha2的其他alpha，alpha、alpha2变成了相互影响的二维变量。
我们的目的就是通过优化符合区间的alpha从而计算出b，到最后得到最优分割平面方程。
以上是三个辅助模块，接下来是核心模块。
def smoSimple(dataMatIn, classLabels, C, toler, maxIter)
# dataMatin数据矩阵， classLabels标签矩阵，C常数作用后面会提到，toler容错值，maxIter最大迭代数
dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
b=0; m,n = shape(dataMatrix) #初始化b
alphas = mat(zeros((m,1))) #alpha初始化0矩阵
iter =0
while (iter < maxIter):
alphaPairsChanged = 0
for i in range(m):
fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
#f(x)=w^T*x+b w=\sum alpha_i * label_i * x_i
Ei = fXi - float(labelMat[i]) #计算预测值和真实结果的误差
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei>toler) and (alphas[i]>0)):
#我们选取误差超过容错值的alpha值进行优化，同时，对于alpha=C/0的，已经在边界，因此不对这些进行优化
j = selectJrand(i,m) # alpha_i 确定了后，确定alpha_j，采用随机
fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T))+b
Ej = fXj - float(labelMat[j]) #计算预测值和真实结果的误差
alphaIold = alphas[i].copy(); #备份alpha优化前的值
alphaJold = alphas[j].copy(); #同上
if (labelMat[i] != labelMat[j]): #涉及到给alpha_i、alpha_j的值域，该内容请看前篇。
L = max(0, alphas[j] - alphas[i])
H = min(C, C+alphas[j] - alphas[j])
else:
L = max(0, alphas[j]+alphas[i]-C)
H = min(C, C+alphas[j]+alphas[j])
if L==H:
print("L==H"); #无效，直接跳出循环
continue
eta = 2.0 * dataMatrix[i,:] * dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T-dataMatrix[j,:]*dataMatrix[j,:].T #衡量x_i和x_j相似性的计算
if eta >= 0:
print("eta>=0)";
continue
alphas[j] -= labelMat[j]*(Ei - Ej)/eta #对alpha_j进行优化
alphas[j] = clipAlpha(alphas[j], H,L) #保证alpha_j正确范围内
if (abs(alphas[j] - alphaJold) <0.00001): #检验alpha_j变化浮动，如果浮动不大认为alpha_j优化结束
print("j not moving enough");
continue
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold-alphas[j]) #对alpha_i进行优化
#接下来计算b的阈值
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
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("iter: %d i:%d, pairs changed %d" %(iter, i, alphaPairsChanged)
if (alphaPairsChanged == 0):
iter += 1
else:
iter = 0
print("iteration number:%d" % iter)
return(b,alphas)
