libSVM源码分析

最新推荐文章于 2020-12-16 16:23:45 发布

le_zhou

最新推荐文章于 2020-12-16 16:23:45 发布

阅读量1.3k

点赞数

文章标签： svm 梯度算法计算机

转载请注明原载地址：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,&param);

->统计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(&param);

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个类的组合方式为：

[cpp] view plain copy

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、计算梯度

主要代码如下：

[cpp] view plain copy

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中α₁由α₂表示，即：，结合上图由解析几何可得α₂的取值范围：

（9.2）

经过一系列变换，可以得到的α₂更新值α₂^new：

（9.3）

结合9.2和9.3式得到α₂^new最终表达式：

（9.4）

得到α₂^new后，就可以由9.1式求α₁^new。

libSVM实现：

具体操作的时候，把选择后的序号i和j代替这里的2和1就可以了。当然，编程时，这些公式还是太抽象。对于9.2式，还需要具体细分。比如，对于y₁y₂=-1时的L = max(0,α₂- α₁)，是0大还α₂- α₁是大的问题。总共需要分8种情况。至于程序中在一个分支中给α₁^new、α₂^new同时赋值，是因为两者之间存在的关系：

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