支持向量机smo算法C语言,SVM支持向量机（四）R语言实现、SMO算法

一、求解支持向量机。

上篇笔记讲到，如何求解拉格朗日乘子向量。基本的想法就是，每次选出两个乘子，对其他的乘子赋值，此时，只剩两个乘子。问题变成了一个两元一次方程和求二元函数最小值的问题。如果乘子可以更新(既违反了KKT条件)，则把其中一个乘子用令一个乘子代替，带入到二元函数中，再求函数取最小值时(通过公式可以看出这是一个开口向上的抛物线)，未知数的值。重复上面的过程直到所有的乘子都稳定下来，不再发生变化。此时问题求解成功。

记要更新的乘子为ai,aj，aj的更新过程如下

其中E1，E2是我们分类的预测值和实际值的差，既wx+b - y 。

K为核函数。 aj每次更新时等于更新前的值加上y2(E1-E2)/n ，这个公式实际上就是上面说的，对开口向上的抛物线求导，令导数等于0求得的解。不过这个公式不太好推导，下面参考的博客里有详细的推导过程。这种写法同时体现了坐标下降的思路，和梯度下降有点像，但梯度下降每次向梯度方向下降固定的步长，而坐标下降一次性达到局部最优的位置。但是，a2new只是抛物线顶点，还不是真正的解。因为还有限制条件，将限制条件画到图上，就是一个正方形，而ai,aj的关系是一条直线，直线截距可正可负，所以要分情况看。如果根据抛物线导数求解的a2new没有落在平行线上和正方形内时，应取直线与正方形的边界值。根据aj,ai的符号是否相同，限制条件的边界值也有所不同：

ai,aj异号

ai,aj同号

根据上面的公式，每次可以更新一个乘子，由于两个乘子是线性关系，所以另外一个乘子往相反的方向更新即可。由于更新乘子的过程中，使用到了b，所以每次更新完，还得把b也更新一下。b的更新公式为

当某个乘子a满足 0 < a < C 时，这个乘子所对应的样本就是支持向量，因此有wx+b=1或者wx+b=-1, 因为w可以由乘子计算出来，后面的正负一就是样本的分类，所以，b可以求出来，推导之后b的计算方式也可以写成上面的更新公式。如果两个乘子都不满足0 < a < C ,也就是说两个样本点都不是支持向量，此时，b不能精确到底是多少，但可以肯定在b1和b2之间，这里一般写成两者平均数。

上面的公式展示了每一次是如何更新一对乘子，并更新b的。

在简单求解的过程中，可以随机选取一些乘子来更新。如果随机选取到的乘子都不需要再更新，既所有的乘子都满足KKT条件，当这样的情况发生到一定次数时，停止迭代，求解完成。

下面我们用R语言实现上面算法，测试数据集为《机器学习实战》第6章的简单算法测试数据。原始数据集图：

使用如下代码，运行后，画出分隔平面

src 

names(src) 

#作原始数据图像

library(ggplot2)

qplot(x,y,data=src,geom="point",xlab="x",ylab="y",color=label,position="jitter")

#挑选一个随机数，当选择ai后，随机选取一个aj

selectRandom 

j 

while(j == i){

j 

}

return(j)

}

#选择下一个拉格朗日乘子，如果顶点在限制条件外，应取边界点

nextAlpha 

#cat("aj,H,L is ",aj,H,L,"\n")

if(aj > H){

aj 

}

if(aj 

aj 

}

return(aj)

}

#定义一个核函数，为方便，先实现内积

kernel 

if(type == "linear"){

return(sum(vector1*vector2))

}

}

#创建核矩阵，避免重复计算

createKernelMatrix 

n 

km 

for(i in 1:n){

for(j in i:n){

value 

km[i,j] 

km[j,i] 

}

}

return(km)

}

#简单支持向量机算法

simpleSVM 

#初始化数据

dataSet 

label 

aSet 

b 

kernelMatrix 

#开始迭代

iter 

while( iter 

changedCount 

for(i in 1:length(label)){

j 

x1 

x2 

E1 

#是否可以优化？ 违反了KKT条件的都可以优化

if(((E1*label[i]  miss) && (aSet[i] > 0))){

E2 

H 

L 

if(label[i]*label[j] > 0){

L 

H 

}else{

L 

H 

}

if(L==H){

cat("L=H continue \n")

next

}

nita 

if(nita == 0){

cat("nita is 0 \n")

next

}

newAj 

oldAj 

aSet[j] 

oldAi 

aSet[i] 

#开始更新b

if(aSet[i] > 0 && aSet[i] 

b 

label[j]*(aSet[j]-oldAj)*kernelMatrix[i,j]

}else{

if(aSet[j] > 0 && aSet[j] 

b 

label[j]*(aSet[j]-oldAj)*kernelMatrix[j,j]

}else{

b1 

label[j]*(aSet[j]-oldAj)*kernelMatrix[i,j]

b2 

label[j]*(aSet[j]-oldAj)*kernelMatrix[j,j]

b 

}

}

changedCount 

}

}

if(changedCount == 0){

iter 

}else{

iter 

}

}

temp 

w1 

w2 

w 

return(list("w"=w,"b"=b,"a"=aSet))

}

model 

src$calY 

p 

p 

p

感觉还行，不过代码里很多地方写死了，只能用于2维数据，为了不写的太复杂，自己都搞不清，就这样了。

二、SMO算法

在上面的实现中，我们顺序选取ai然后随机选取aj。这个过程浪费了非常多计算量。从经验来看，绝大多数的样本点都不是支持向量，当优化到该样本时，乘子会变成0。而乘子变成0或C之后基本不会再发生变化，反而那些优化后仍然处于0和C之间的乘子，往往需要不断的优化。SMO算法遵循“启发式”选择方法，既优先优化那些可以一次性优化到位的乘子。然后再优化其他的乘子。SMO通过两层循环来挑选要优化的乘子。

SMO称选择第一个变量的过程为外层循环。外层训练在训练样本中选取违法KKT条件最严重的样本点。并将其对应的变量作为第一个变量。 该检验是在ε范围内进行的。在检验过程中，外层循环首先遍历所有满足条件0

优先选择遍历非边界数据样本，因为非边界数据样本更有可能需要调整，边界数据样本常常不能得到进一步调整而留在边界上。由于大部分数据样本都很明显不可能是支持向量，因此对应的α乘子一旦取得零值就无需再调整。遍历非边界数据样本并选出他们当中违反KKT 条件为止。当某一次遍历发现没有非边界数据样本得到调整时，遍历所有数据样本，以检验是否整个集合都满足KKT条件。如果整个集合的检验中又有数据样本被 进一步进化，则有必要再遍历非边界数据样本。这样，不停地在遍历所有数据样本和遍历非边界数据样本之间切换，直到整个样本集合都满足KKT条件为止。以上 用KKT条件对数据样本所做的检验都以达到一定精度ε就可以停止为条件。如果要求十分精确的输出算法，则往往不能很快收敛。

对整个数据集的遍历扫描相当容易，而实现对非边界αi的扫描时，首先需要将所有非边界样本的αi值(也就是满足0

在选择第一个αi后，算法会通过一个内循环来选择第二个αj值。因为第二个乘子的迭代步长大致正比于|Ei-Ej|，所以我们需要选择能够最大化|Ei-Ej|的第二个乘子(选择最大化迭代步长的第二个乘子)。在这里，为了节省计算时间，我们建立一个全局的缓存用于保存所有样本的误差值，而不用每次选择的时候就重新计算。我们从中选择使得步长最大或者|Ei-Ej|最大的αj。

后期实现。