libSVM源碼分析

轉載請註明原載地址：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);

一、read_problem()

prob.y //記錄每行樣本所屬類別

prob.x //指針數組(L維)，每一個指針指向x_sapce(實際存儲特徵詞)的一維

x_space//實際的存儲結構，記錄全部樣本的特徵詞（L*(k+1)個），能夠形象化爲L維，雖然每一維的長度可能不一樣。

prob.y、prob.x與x_space的關係以下圖所示：

二、統計classes數據

使用如下變量遍歷全部樣本，統計數據。

label[i] //記錄類別

count[i] //記錄類別中樣本的數量

index[i] //記錄位置爲i的樣本的類別

nr_class //索引類別的數目

三、訓練數據分組

訓練數據進行分組時使用到了如下數據：

int *start = Malloc(int,nr_class);

svm_node **x = Malloc(svm_node *,l);

二者之間經過index進行過渡，由於index記錄了位置i的樣本的類別，每一個類別在start中只有一個位置，即該類別在x中的起始的索引。x是各種的排列順序是按照原始樣本中各種出現的前後順序排列的，prob中則是原始樣本的label序號排列，而start中記錄的是各種的起始序號，而這個序號是在x的序號。

四、訓練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」，咱們就繼續往左走，這樣一直問下去，就能夠獲得分類結果。

五、計算梯度

主要代碼以下：

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)   

六、數據選擇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。

七、數據縮放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的大小。

八、迭代優化中止準則

    LibSVM程序中，中止準則蘊含在了函數select_working_set(i,j)返回值中。也就是，當找不到符合6.5式的樣本時，那麼理論上就達到了最優解。可是，實際編程時，因爲KKT條件仍是蠻苛刻的，要進行適當的放鬆。令： 

                                    （8.1）

由6.4式可知，當全部樣本都知足KKT條件時，gi ≤ -gj

加一個適當的寬鬆範圍ε，也就是程序中的eps，默認爲0.001，那麼最終的中止準則爲：

                     gi ≤ -gj +ε  →    gi + gj ≤ε

九、因子α的更新

理論依據：

因爲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])的值，代入以上公式可得另一個的值。

十、更新G和G_Bar

根據的變化更新G(i)，更新alpha_status和較簡單，根據alpha狀態先後是否有變化，適當更新，更新的內容參考公式5.2

十一、截距b的計算

b計算的基本公式爲：

            （11.1）

理論上，b的值是不定的。當程序達到最優後，只要用任意一個標準支持向量機（0<alpha[i]<C）的樣本帶入11.1式，獲得的b值都是能夠的。目前，求b的方法也有不少種。在libSVM中，分別對y=+1和y=-1的兩類全部支持向量求b，而後取平均值：

                                 

十二、計算目標函數值

由於目標值的計算公式爲：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程序代碼導讀       劉國平

序列最小化方法                 羅林開