从零开始学习SVM（四）---SMO算法线性分类与代码精解析

从零开始学习SVM（三）—SMO算法原理了解SMO算法来龙去脉之后,我们终于到了实践这步了。只有在学习理论基础之后并实践才能更好的理解这个精妙的算法。我建议先搞懂SMO算法原理，也就是我上一篇博客，弄懂了各个公式的来龙去脉后，撸起码来才能得心应手。首先

selectJrand

函数是用来随机选择另外一个alpha,在numpy方法里面

np.random.uniform

，这个uniform是选择另外一个alpha采用均匀分布的方式选取，也就说0到m任意选取一个值，并且概率一致

def selectJrand(i,m):
'''
随机选择j 另外一个aplha
'''
j=i #we want to select any J not equal to i
while (j==i):
j = int(np.random.uniform(0,m))#概率相同的情况下，随机选取
return j

因为alpha的取值的有上下限的，所以我们需要对超过范围的alpha进去截断操作

def clipAlpha(alpha_j,H,L):
'''
更新的alpha要满足kkt条件，如果超出取值范围就，截断，把值设置为上下限
'''
if alpha_j>H:
alpha_j=H
if alpha_j<L:
alpha_j=L
return alpha_j

这是完整的简单的SMO算法代码

def smoSimple(X,Y,C,toler,maxIter):
X=np.mat(X)
Y=np.mat(Y).transpose()
b=0
m,n=X.shape
#每个样本都具有一个alpha,所以给各个样本初始化一个alpha
alphas=np.mat(np.zeros([m,1]))
iterNum=0
while(iterNum<maxIter):
alpha_pairs_changed=0
for i in range(m):
#超平面函数公式
fxi=float(np.multiply(alphas,Y).T*(X*X[i,:].T))+b
Ei=fxi-float(Y[i])
'''
这里不管是正间隔还是负间隔都会被测试，并且alpha的取值范围0到c，否则
它们将会被赋值在间隔上，也就没有优化的意义了
'''
if ((Y[i]*Ei<-toler)  and (alphas[i]<C))or ((Y[i]*Ei>toler) and (alphas[i]>0)) :
j = selectJrand(i,m)
#                随机抽取一个样本
fxj=float(np.multiply(alphas,Y).T*(X*X[j,:].T))+b
#得到偏差
Ej=fxj-float(Y[j])
#现在是对象了，所以得用深拷贝
alphaIold=alphas[i].copy()
alphaJold=alphas[j].copy()
#考虑两者同号或者异号
if Y[i]!=Y[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[i]+alphas[j])
if H==L:#如果上下限相同，则是直线上的一个点
print('H==L')
continue
eta=2.0*X[i,:]*X[j,:].T-X[i,:]*X[i, :].T-X[j,:]*X[j,:].T
if eta>=0:#eta大于零则没意义了，相同就说明是同一个样本
continue
alphas[j]-=Y[j]*(Ei-Ej)/eta
#处理上下限
alphas[j]=clipAlpha(alphas[j],H,L)
#如果alpha_j 相对于以前仅仅变化一点，则不更新，重新选择一对
if abs(alphas[j]-alphaJold)<0.00001:
print('alpha2的更新程度小，没必要更新这对alpha')
continue
alphas[i]+=Y[j]*Y[i]*(alphaJold-alphas[j])
#更新b
b1=b-Ei-Y[i]*(alphas[i]-alphaIold)*X[i,:]*X[i,:].T-Y[j]*(alphas[j]-alphaJold)*X[i,:]*X[j,:].T
b2=b-Ej-Y[i]*(alphas[i]-alphaIold)*X[i,:]*X[j,:].T-Y[j]*(alphas[j]-alphaJold)*X[j,:]*X[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.
#标记一下，改变了一对alpha
alpha_pairs_changed+=1
print("iter: %d i:%d, pairs changed %d" % (iterNum,i,alpha_pairs_changed))
if alpha_pairs_changed==0:
iterNum+=1
else:
iterNum=0
print('iteration number: %d'%iterNum)
return b,alphas

不要被这大片的代码吓到，现在我们一步步解析这些代码：首先得到样本与对应的标签矩阵，np.mat是转换成矩阵的方式，用这个方法好处就是直接用*就可以表达出内积而不用使用np.dot方法。每个样本都对应一个alpha,所以我们应该初始化几个alpha呢？没错就是m个，因此我们初始化一个m*1的alpha零矩阵：

alphas=np.mat(np.zeros([m,1]))

X=np.mat(X)
Y=np.mat(Y).transpose()
b=0
m,n=X.shape
#每个样本都具有一个alpha,所以给各个样本初始化一个alpha
alphas=np.mat(np.zeros([m,1]))
iterNum=0

fxi=float(np.multiply(alphas,Y).T*(X*X[i,:].T))+b

这是什么？我相信你不会忘记超平面的函数，把W替换为alpha的方式就得到了上述代码。其实这个都好理解，最难理解的是

if ((Y[i]*Ei<-toler)  and (alphas[i]<C))or ((Y[i]*Ei>toler) and (alphas[i]>0)) :

这个条件没少让我下功夫，现在我们来彻底弄懂这个条件，相信会让你更加深刻的理解KKT条件，二话不说上KKT条件



函数中的toler就是我们第二节的松弛变量也叫容忍系数，对应ζ. 把

Y[i]*Ei

公式化，并展开得yi(f(xi)−yi)>ζ展开yif(xi)−1−ζ>0；没道理啊怎么会得到这么个玩意，这是哪里的KKT条件？且慢，透过现象看本质。其实这从一开始就错了，好好思考一下什么样的a不满足KKT条件呢？当然是分错的样本不满足啊，那就是yif(xi)−1+ζ<0了，这样就算错了，对。两边同时除以yi,当yi=1时候得yi(f(xi)−1yi)<ζyi=>yi(f(xi)−yi)<ζ=>yiEi<−toler

同理可得，当yi=−1,则：yi(f(xi)−1yi)>ζyi=>yi(f(xi)−yi)>ζ=>yiEi>toler

对于这些条件，alpha为0或者C就不用考虑了，这些样本将在边界上，而且alpha不会被改变，所以只要考虑在（0，C）这个区间内。alpha在初始化的时候为0，并且后面更新alpha对的时候，会最后调用截断alpha函数，所以所有的alpha都会在[0,C]区间内。所以上面的两个条件可以随意搭配α>0,α<C这两个条件。我为了验证这个想法，我把这两个条件全部去掉，重新运行模型，完全正常




运行次数不同得到的支持向量的个数也不同，这是正常情况。

#更新b
b1=b-Ei-Y[i]*(alphas[i]-alphaIold)*X[i,:]*X[i,:].T-Y[j]*(alphas[j]-alphaJold)*X[i,:]*X[j,:].T
b2=b-Ej-Y[i]*(alphas[i]-alphaIold)*X[i,:]*X[j,:].T-Y[j]*(alphas[j]-alphaJold)*X[j,:]*X[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.

上面这段来自于下面这个公式：



这个公式，至于选择哪一个b为最终的b，当然是哪个alpha满足条件用哪一个，在这一点b的选择略显得随意。本以为实践应该很简单，但是真正去做的时候却发现有很多小细节，比较零散的细节我都注释到代码中了，查了很久的资料但是网上基本流传的都是这段代码且没有详细的解释，书中解释的也很泛，索性写一篇比较low的解释博客吧。希望能帮到一些跟我一样初入门的童鞋，少一点痛苦多一点捷径。虽然我的粉丝少的可怜，基本都是威逼利诱下才关注我的师弟师妹，在此我还是打一份广告吧，欢迎关注我的微信公共号，哈哈！！