机器学习实战支持向量机SVM 代码解析

《机器学习实战》用代码实现了算法，理解源代码更有助于我们掌握算法，但是比较适合有一定基础的小伙伴。svm这章代码看起来风轻云淡，实则对于新手来说有（shi）点(fen)晦涩，必须先搞清楚svm原理和具体推导优化步骤。这里推荐一个知乎的回答，支持向量机 SVM（非常详细）。这篇文章只推到了优化目标公式，最后的优化过程用的是序列最小最优化（SMO）算法。具体过程参见SMO算法详解

下面是github 上写了注释的代码svm支持向量机

smoSimple(dataMatIn, classLabels, C, toler, maxIter):
"""smoSimple
Args:dataMatIn    特征集合classLabels  类别标签C   松弛变量(常量值)，允许有些数据点可以处于分隔面的错误一侧。控制最大化间隔和保证大部分的函数间隔小于1.0这两个目标的权重。可以通过调节该参数达到不同的结果。toler   容错率（是指在某个体系中能减小一些因素或选择对某个系统产生不稳定的概率。）maxIter 退出前最大的循环次数
Returns:b       模型的常量值alphas  拉格朗日乘子
"""
dataMatrix = mat(dataMatIn)
# 矩阵转置 和 .T 一样的功能
labelMat = mat(classLabels).transpose()
m, n = shape(dataMatrix)# 初始化 b和alphas(alpha有点类似权重值。)
b = 0
alphas = mat(zeros((m, 1)))# 没有任何alpha改变的情况下遍历数据的次数
iter = 0
while (iter < maxIter):# w = calcWs(alphas, dataMatIn, classLabels)# print("w:", w)# 记录alpha是否已经进行优化，每次循环时设为0，然后再对整个集合顺序遍历alphaPairsChanged = 0for i in range(m):# print 'alphas=', alphas# print 'labelMat=', labelMat# print 'multiply(alphas, labelMat)=', multiply(alphas, labelMat)# 我们预测的类别 y[i] = w^Tx[i]+b; 其中因为 w = Σ(1~n) a[n]*label[n]*x[n]fXi = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[i, :].T)) + b# 预测结果与真实结果比对，计算误差EiEi = fXi - float(labelMat[i])# 约束条件 (KKT条件是解决最优化问题的时用到的一种方法。我们这里提到的最优化问题通常是指对于给定的某一函数，求其在指定作用域上的全局最小值)# 0<=alphas[i]<=C，但由于0和C是边界值，我们无法进行优化，因为需要增加一个alphas和降低一个alphas。# 表示发生错误的概率：labelMat[i]*Ei 如果超出了 toler， 才需要优化。至于正负号，我们考虑绝对值就对了。'''# 检验训练样本(xi, yi)是否满足KKT条件yi*f(i) >= 1 and alpha = 0 (outside the boundary)yi*f(i) == 1 and 0<alpha< C (on the boundary)yi*f(i) <= 1 and alpha = C (between the boundary)'''if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):

对于条件((labelMat[i]*Ei < -toler) and (alphas[i] < C))：
αi<C,即对于边界上或边界外的点，
当yi=1时,labelMat[i]*Ei=yi(wxi+b-yi)=f(i)-1<-toler，f(i)<1-toler；当yi=-1时，f(i)>-1+toler
也就是边界以外的点xi经过wxi+b分类后得到的结果超出了设置的容错率，即分类错误，需要继续优化
同理，对于条件((labelMat[i]*Ei > toler) and (alphas[i] > 0))：
当yi=1时，labelMat[i]*Ei=yi(f(i)-yi)=f(i)-1>toler,即f(i)>toler+1；当yi=-1时，f(i)<-1-toler,说明样本i是边界外的点。那么αi，也就是样本i对优化目标限制条件的贡献应该为零，而条件alphas[i]>0却使样本i对限制条件的贡献大于零，说明αi是不对的，需要继续优化。

插播一条对αi的理解：SVM就是一个不等式约束的优化问题，约束条件就是每个正样本到分类面（超平面）的（带符号的）距离大于+1/ω，每个负样本到超平面的距离小于-1/ω，优化目标就是1/ω。通俗来说就是把每个样本都分类正确的情况下，要使正样本和负样本到超平面的距离最大，距离越大就越安全。要想把所有样本都分类正确其实不需要把所有样本到超平面的距离都求出来然后让他们大于1/ω，只需要让那些离超平面最近的点（也就是支持向量）到超平面的距离大于1/ω就行了。也就是说边界以外的样本对超平面的优化是没有约束的。
回到我们构造的Lagrange函数:
　　　　　　　　　L(x,α)=f(x)+∑αihi(x)
其中αi是样本i对优化目标f(x)约束的权重，边界外的点对超平面没有约束，所以αi等于零。对于下图中的坏点，αi就代表无法正确分类的样本对目标函数的惩罚，应该取一个较大的值，假如我们设定它为C。而对于刚好位于边界上的点（支持向量），他们的权重应该大于零而小于C。这其实就是所谓的KKT条件了。
我们看到αi对样本i的各个维度的值其实是没有调节作用的，它是通过调整每个支持向量在限制条件中的权重，而每个支持向量的不同维度特征值又是不同的，从而影响了ω在各个维度的值，最终影响超平面。

另外画了张图来帮助理解

            # 如果满足优化的条件，我们就随机选取非i的一个点，进行优化比较j = selectJrand(i, m)# 预测j的结果fXj = float(multiply(alphas, labelMat).T*(dataMatrix*dataMatrix[j, :].T)) + bEj = fXj - float(labelMat[j])alphaIold = alphas[i].copy()alphaJold = alphas[j].copy()# L和H用于将alphas[j]调整到0-C之间。如果L==H，就不做任何改变，直接执行continue语句# labelMat[i] != labelMat[j] 表示异侧，就相减，否则是同侧，就相加。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是alphas[j]的最优修改量，如果eta==0，需要退出for循环的当前迭代过程# 参考《统计学习方法》李航-P125~P128<序列最小最优化算法>eta = 2.0 * dataMatrix[i, :]*dataMatrix[j, :].T - dataMatrix[i, :]*dataMatrix[i, :].T - dataMatrix[j, :]*dataMatrix[j, :].Tif eta >= 0:print("eta>=0")continue# 计算出一个新的alphas[j]值alphas[j] -= labelMat[j]*(Ei - Ej)/eta# 并使用辅助函数，以及L和H对其进行调整alphas[j] = clipAlpha(alphas[j], H, L)# 检查alpha[j]是否只是轻微的改变，如果是的话，就退出for循环。if (abs(alphas[j] - alphaJold) < 0.00001):print("j not moving enough")continue# 然后alphas[i]和alphas[j]同样进行改变，虽然改变的大小一样，但是改变的方向正好相反alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])# 在对alpha[i], alpha[j] 进行优化之后，给这两个alpha值设置一个常数b。# w= Σ[1~n] ai*yi*xi => b = yj- Σ[1~n] ai*yi(xi*xj)# 所以：  b1 - b = (y1-y) - Σ[1~n] yi*(a1-a)*(xi*x1)# 为什么减2遍？ 因为是 减去Σ[1~n]，正好2个变量i和j，所以减2遍b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i, :]*dataMatrix[i, :].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i, :]*dataMatrix[j, :].Tb2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i, :]*dataMatrix[j, :].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j, :]*dataMatrix[j, :].Tif (0 < alphas[i]) and (C > alphas[i]):b = b1elif (0 < alphas[j]) and (C > alphas[j]):b = b2else:b = (b1 + b2)/2.0

这里b的更新注释讲得比较简略，也没有说明后面的if条件是怎么来的。b的更新需要使得优化后的α满足KKT条件，更详尽的推导可以参考知乎b值更新公式推导

            alphaPairsChanged += 1print("iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))# 在for循环外，检查alpha值是否做了更新，如果更新则将iter设为0后继续运行程序# 直到更新完毕后，iter次循环无变化，才退出循环。if (alphaPairsChanged == 0):iter += 1

什么情况下alphaPairsChanged 会等于零呢：
遍历所有样本，每个样本都满足以下条件之一
1. 满足KKT条件
2. 满足第一个continue条件，L=H。新的α取值范围上限等于下限，无法优化
3. 满足第二个continue条件，eata=0。最佳优化值为零
当所有样本要么满足KKT条件，不满足KKT条件的点也没有了优化的空间，自然可以进行下一次循环了。

    else:iter = 0print("iteration number: %d" % iter)
return b, alphas

看完这个之后应该对SVM有了比较全面的了解了，推荐再去看一下吴恩达老师对SVM的讲解，吴恩达机器学习。虽然讲的是一个东西，但是用的是完全不同的思路。首先借助逻辑回归的代价函数构造了一个满足SVM模型的代价函数，然后就直接求代价函数最小值就好了。是不是很神奇，其实是跟拉格朗日乘子法殊途同归了。优化目标和约束条件都是一样的，Andrew还以向量内积的角度解释了svm最后为什么会找到最大间隔的超平面，妙啊。
　　看看不同的思路，融会贯通，不禁感叹前辈们都太聪明了吧！