转载请注明原载地址:http://blog.csdn.net/xinhanggebuguake/article/details/8705648
在此之前,上海交大模式分析与机器智能实验室对2.6版本的svm.cpp做了部分注解,《LibSVM学习(四)——逐步深入LibSVM》也介绍了libSVM的思路,很精彩。而我写这篇博客更侧重与理解算法流程与具体代码的结合点。(环境:LibSVM2.6 C-SVC型SVM RBF核函数)
函数调用流程:
svm-train.c
main()
parse_command_line();//解析命令行,将数据读入param,并获取input file、model file
read_problem();//读取input file中的数据到prob中。
(do_cross_validation();//该函数将试验所有的核函数,根据交叉验证选择最好的核函数)
svm_train(&prob,¶m);
->统计classes的数量以及每个classes下样本数量
->把相同类别的训练数据分组,每个分组开始的索引记录在start数组里。
->计算每个类别的惩罚因子C
->训练k*(k-1)/2个分类器模型
->svm_train_one();
->solve_c_svc();
->s.Solve();
->初始化alpha_status、active_set、active_size
->求梯度
迭代优化: ->do_shrinking(); //把数据分成active_size和active_size-L的部分集中排序。
->select_working_set(); //选择两个样本
->更新alpha[i]和alpha[j]的值
->更新G和G_Bar
->calculate_rho();//计算b值
->计算目标值
svm_save_model(model_file_name,model);
svm_destroy_model(model);
svm_destroy_param(¶m);
1、read_problem()
prob.y //记录每行样本所属类别
prob.x //指针数组(L维),每个指针指向x_sapce(实际存储特征词)的一维
x_space//实际的存储结构,记录所有样本的特征词(L*(k+1)个),可以形象化为L维,虽然每一维的长度可能不同。
prob.y、prob.x与x_space的关系如下图所示:
2、统计classes数据
使用以下变量遍历所有样本,统计数据。
label[i] //记录类别
count[i] //记录类别中样本的数量
index[i] //记录位置为i的样本的类别
nr_class //索引类别的数目
3、训练数据分组
训练数据进行分组时使用到了以下数据:
int *start = Malloc(int,nr_class);
svm_node **x = Malloc(svm_node *,l);
两者之间通过index进行过渡,因为index记录了位置i的样本的类别,每个类别在start中只有一个位置,即该类别在x中的起始的索引。x是各类的排列顺序是按照原始样本中各类出现的先后顺序排列的,prob中则是原始样本的label序号排列,而start中记录的是各类的起始序号,而这个序号是在x的序号。
4、训练k*(k-1)/2个分类器模型
svm对于多类别的分类方法有多种,但将实现分为两个过程:训练阶段,判别阶段。
1)1-V-R方式
对于k类问题,把其中某一类的n个训练样本视为一类,所有其他类别归为另一类,因此共有k个分类器。最后,判别使用竞争方式,也就是哪个类得票多就属于那个类。
2)1-V-1方式
即one-against-one方式。该方法把其中的任意两类构造一个分类器,共有(k-1)×k/2个分类器。最后判别也采用竞争方式。
3)1-V-1在libSVM中的实现
LibSVM采用的是1-V-1方式,因为这种方式思路简单,并且许多实践证实效果比1-V-R方式要好。该方法在训练阶段采用1-V-1方式,而判别阶段采用一种两向有向无环图的方式。
训练阶段:
上图是一个5类1-V-1组合的示意图,红色是0类和其他类的组合,紫色是1类和剩余类的组合,绿色是2类与右端两类的组合,蓝色只有3和4的组合。因此,对于nr_class个类的组合方式为:
- for(i = 0; i < nr_class; i ++)
- {
- for(j = i+1; i < nr_class; j ++)
- {
- 类 i –V – 类 j
- }
- }
判别阶段:
在对一篇文章进行分类之前,我们先按照下面图的样子来组织分类器(如你所见,这是一个有向无环图,因此这种方法也叫做DAG-SVM)
在分类时,我们就可以先问分类器“1对5”(意思是它能够回答“是第1类还是第5类”),如果它回答5,我们就往左走,再问“2对5”这个分类器,如果它还说是“5”,我们就继续往左走,这样一直问下去,就可以得到分类结果。
5、计算梯度
主要代码如下:
- G[i] = b[i];
- G_bar[i] = 0;
- Qfloat *Q_i = Q.get_Q(i,l);
- for(i=0;i<l;i++)
- {
- for(j=0;j<l;j++)
- G[j] += alpha_i*Q_i[j];
- for(j=0;j<l;j++)
- G_bar[j] += get_C(i) * Q_i[j];
- }
首先,Q.get_Q(i,l)返回data,而
data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));
翻译成公式,即:
所以,以上计算梯度的代码翻译成公式,则:
G为:
(5.1)
G_bar为:
(5.2)
6、数据选择select_working_set(i,j)
理论依据:
对于样本数量比较多的时候(几千个),SVM所需要的内存是计算机所不能承受的。目前,对于这个问题的解决方法主要有两种:块算法和分解算法。这里,libSVM采用的是分解算法中的SMO(串行最小化)方法,其每次训练都只选择两个样本。我们不对SMO做具体的讨论,要想深入了解可以查阅相关的资料,这里只谈谈和程序有关的知识。
一般SVM的对偶问题为:
S.t. (6.1)
SVM收敛的充分必要条件是KKT条件,其表现为:
(6.2)
由6.1式求导可得:
(6.3)
进一步推导可知:
(6.4)
也就是说,只要所有的样本都满足6.4式,那么得到解就是最优值。因此,在每轮训练中,每次只要选择两个样本(序号为i和j),是最违反KKT条件(也就是6.4式)的样本,就能保证其他样本也满足KKT条件。序号i和j的选择方式如下:
(6.5)
libSVM实现:
由公式5.1和6.5可知,select_working_set的过程,只跟G_bar和C有关,所以根据is_lower_bound与is_upper_bound判断C的范围,再根据y[i],可以将公式6.5分为8个分支。循环遍历所有样本,就能查找到最违反KTT条件的样本的index。
7、数据缩放do_shrinking()
上面说到SVM用到的内存巨大,另一个缺陷就是计算速度,因为数据大了,计算量也就大,很显然计算速度就会下降。因此,一个好的方式就是在计算过程中逐步去掉不参与计算的数据。因为,实践证明,在训练过程中,alpha[i]一旦达到边界(alpha[i]=0或者alpha[i]=C),alpha[i]值就不会变,随着训练的进行,参与运算的样本会越来越少,SVM最终结果的支持向量(0<alpha[i]<C)往往占很少部分。
LibSVM采用的策略是在计算过程中,检测active_size中的alpha[i]值,如果alpha[i]到了边界,那么就应该把相应的样本去掉(变成inactived),并放到栈的尾部,从而逐步缩小active_size的大小。
8、迭代优化停止准则
LibSVM程序中,停止准则蕴含在了函数select_working_set(i,j)返回值中。也就是,当找不到符合6.5式的样本时,那么理论上就达到了最优解。但是,实际编程时,由于KKT条件还是蛮苛刻的,要进行适当的放松。令:
(8.1)
由6.4式可知,当所有样本都满足KKT条件时,gi ≤ -gj
加一个适当的宽松范围ε,也就是程序中的eps,默认为0.001,那么最终的停止准则为:
gi ≤ -gj +ε → gi + gj ≤ε
9、因子α的更新
理论依据:
由于SMO每次都只选择2个样本,那么4.1式的等式约束可以转化为直线约束:
(9.1)
转化为图形表示为:
把式9.1中α1由α2 表示,即:,结合上图由解析几何可得α2的取值范围:
(9.2)
经过一系列变换,可以得到的α2更新值α2new:
(9.3)
结合9.2和9.3式得到α2new最终表达式:
(9.4)
得到α2new后,就可以由9.1式求α1new。
libSVM实现:
具体操作的时候,把选择后的序号i和j代替这里的2和1就可以了。当然,编程时,这些公式还是太抽象。对于9.2式,还需要具体细分。比如,对于y1y2=-1时的L = max(0,α2- α1),是0大还α2- α1是大的问题。总共需要分8种情况。至于程序中在一个分支中给α1new、α2new同时赋值,是因为两者之间存在的关系:
diff = alpha[i] - alpha[j];
依据公式9.4,最内层对alpha[i](或alpha[j])判断可以得出alpha[i] (或alpha[j])的值,代入以上公式可得另外一个的值。
10、更新G和G_Bar
根据的变化更新G(i),更新alpha_status和较简单,根据alpha状态前后是否有变化,适当更新,更新的内容参考公式5.2
11、截距b的计算
b计算的基本公式为:
(11.1)
理论上,b的值是不定的。当程序达到最优后,只要用任意一个标准支持向量机(0<alpha[i]<C)的样本带入11.1式,得到的b值都是可以的。目前,求b的方法也有很多种。在libSVM中,分别对y=+1和y=-1的两类所有支持向量求b,然后取平均值:
12、计算目标函数值
因为目标值的计算公式为:1/2*alpha*Sigma (G[i]+b[i]),
G[i]转换为公式为alpha_i*Q_i[j]+b[i]
由于在传递给Solve函数的minus_ones将所有值赋为-1,所以b[i]=-1,以上公式就转换为
1/2*alpha*Sigma (Q_i[j]) + 1/2*alpha*2*(-1);
上面的公式不正是我们的目标函数吗。所以可以理解libSVM中的实现。
总结:由于刚接触相关方面的知识,疏漏之处在所难免,希望各位高手能不吝赐教!
参考资料:
http://blog.csdn.net/xinhanggebuguake/article/details/8705631
http://www.blogjava.net/zhenandaci/archive/2009/03/26/262113.html
http://www.cnblogs.com/biyeymyhjob/archive/2012/07/17/2591592.html
http://blog.csdn.net/flydreamGG/article/details/4470121
libsvm-2[2].8程序代码导读 刘国平
序列最小化方法 罗林开