机器学习-python通过序列最小优化算法(SMO)方法编写支持向量机(SVM)

代码及数据集下载：SVM

学习SVM，然后想自己手写一个，看了一些材料后发现了最流行的SMO，下面总结如下。SMO算法的思想很简单，它将大优化的问题分解成多个小优化的问题。这些小问题往往比较容易求解，并且对他们进行顺序求解的结果与将他们作为整体来求解的结果完全一致。在结果完全一致的同时，SMO的求解时间短很多。在深入SMO算法之前，我们先来了解下坐标下降这个算法，SMO其实基于这种简单的思想的。SMO在论文《Sequential Minimal Optimization: A Fast Algorithm for TrainingSupport Vector Machines》中提出，感兴趣的可以向我要论文。
SVM的对偶约束问题为：

KKT条件为：

---

这里要介绍一个坐标下降（上升）法的概念：

minαW(α1,α2,...,αm)

如果要解如上优化问题，通常采用梯度下降法（若max，则为上升法），每次迭代调整每个αi变量。而坐标下降法为，每次只调整一个α值，其他α在调整过程中保持不变。

for until convergence:    for i in range m:        alpha_i := argmin W(alpha_1,...,alpha_m)

一个二维例子可以很好的解释：

黄色箭头是最快的梯度方向，即为牛顿优化法，但是高维情况计算较慢（矩阵求逆复杂）。下面看坐标下降法，先固定x，让y沿梯度方向到B，再固定y，x沿梯度到C。同理收敛到相同的F点。在维数较高时，坐标下降法更加高效。

---

下面介绍SMO算法：
SMO算法继承了坐标下降法的思想，但是为了满足KKT的等式约束，因此每次选取两个αi值。
一、寻找两个αi的步骤如下:
1.外循环遍历所有αi，选取非边界样本，即0<α<C的点，选择最背离KKT条件的点，因为非边界点才有调整的空间，直到优化完所有的样本点都不违背KKT。选取最大背离为α1。
2.内循环遍历所有非边界样本，选取更新步长最大的αi作为α2。步长与|Ei−E2|成正比，E=u−y，u为SVM输出，y为理想输出。如果E1为正，选最小的E2，如果E1为负，选最大的E2。

Ei=∑j=1N(yjαj<xj,xi>)+b−yi

二、由于不等式的约束，则α1与α2的调成范围可以被限制在一个方框内，又由于等式约束，α1与α2的调成范围被限制在一条直线上。

因此当y1=y2时，范围可如下计算：

当y1!=y2时，范围可如下计算：

三、参数调整方法
1.α调整方法




其中s=y1y2
2.b调整方法

b1=b−E1−y1(αnew1−αold1)<x1,x1>−y2(αnew2−αold2)<x2,x1>b2=b−E2−y1(αnew1−αold1)<x1,x2>−y2(αnew2−αold2)<x2,x2>

如果α1与α2都为非边界样本，则b1=b2。如果其中一个为边界样本，则选取另一个非边界样本对应的b。若都为边界样本，则选取b=b1+b22。
四、流程用伪代码表示如下

for 迭代最大次数：    for 所有数据：        在0<alpha<C,中找到偏离KKT条件最多的alpha1        for 非边界样本(0<alpha<C)：            找到步长最大的alpha2，使得|E1-E2|最大            检查是否满足容忍度                确定L与H                确定eta                计算新的alpha2                限制alpha2的范围                计算新的alpha1                计算b1与b2，选取非边界样本确定的b，若无选(b1+b2)/2    若没有参数更新，则迭代次数+1，否则为0        

---

python简单实现

import numpy as npdef loadDataSet(fileName):      dataMat = []     labelMat = []    with open(fileName) as f:        for line in f.readlines():            line = line.strip().split()            dataMat.append([float(line[0]),float(line[1])])            labelMat.append(int(line[2]))    return dataMat,labelMatdef selectJrand(i,m):    j = i    while(j == i):        j = int(np.random.uniform(0,m))    return jdef clipAlpha(aj,H,L):    if aj > H:        aj = H    if aj < L:        aj = L    return ajdef smoSimple(dataMatIn,classLabels,C,toler,maxIter):    dataMatrix = np.mat(dataMatIn)    labelMat = np.mat(classLabels).transpose()    b = 0    m,n = np.shape(dataMatrix)    alphas = np.mat(np.zeros((m,1)))    iters = 0    while (iters < maxIter):        alphaPairsChanged = 0        for i in range(m):            fxi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b            Ei = fxi - float(labelMat[i])            if( ((labelMat[i]*Ei)<-toler) and (alphas[i]<C) ) or ((labelMat[i]*Ei>toler) and(alphas[i] > 0)):                j = selectJrand(i,m)                fxj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b                Ej = fxj - float(labelMat[j])                alphaIold = alphas[i].copy()                alphaJold = alphas[j].copy()                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                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                alphas[j] -= labelMat[j]*(Ei - Ej)/eta                alphas[j] = clipAlpha(alphas[j],H,L)                if (abs(alphas[j]-alphaJold)<0.00001):                    print('J is not move')                    continue                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold-alphas[j])                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:%d'%(iters,i,alphaPairsChanged))        if(alphaPairsChanged == 0):            iters += 1        else:            iters = 0        print('iteration number :%d'%iters)    return b,alphasif __name__ == '__main__':    data,label = loadDataSet('testSet.txt')    b,alpha = smoSimple(data,label,0.6,0.001,40)    print(b,alpha)

注：编程技巧
1.向量对应元素相乘

aOut[13]: matrix([[2],        [2],        [2]])bOut[14]: matrix([[1],        [2],        [3]])np.multiply(a,b)Out[16]: matrix([[2],        [4],        [6]])

2.过滤numpy的元素

alphas[alphas>0]