LibSVM是是臺灣林智仁(Chih-Jen Lin)教授2001年開發的一套支援向量機的庫，可以很方便的對資料做分類或迴歸。由於LibSVM程式小，運用靈活，輸入引數少，並且是開源的，易於擴充套件，因此成為目前國內應用最多的SVM的庫，同時sklearn.svm也是使用的該庫。

網路上對於LibSVM原始碼的講解有很多，但個人感覺絕大多數講解的不夠清晰，很多都是貼個理論公式圖片再粘段程式碼就一帶而過。並且網路上基本都是對SVC的講解，SVR部分幾乎沒有提及（雖然SVR只是SVC的擴充套件）。因此本篇博文將系統地講解LibSVM中SVR訓練與預測部分的原始碼（想學習SVC的同學同樣適用）。

LibSVM整體流程

train：

//根據svm_type的不同進行初始化
svm_train()
    //根據svm_type的不同調用不同的分類迴歸訓練函式
    svm_train_one()
        //針對epsilon-SVR這一型別進行模型引數初始化
        solve_epsilon_svr()
            //使用SMO演算法求解對偶問題（二次優化問題）
            Solver::Solve()
                    //每隔若干次迭代進行一次shrinking，對樣本集進行縮減降低計算成本
                    Solver::
 
do_shrinking()
                    //若滿足停止條件則進行梯度重建並跳出迴圈
                    Solver::reconstruct_gradient()
                    //選擇出當前最大違反對i,j
                    Solver::select_working_set()
                //計算引數優化後的rho
                Solver::caculate_rho()
            //得到優化後的alpha和SolutionInfo物件si
        //得到優化後的alpha和SolutionInfo物件si
 

    //得到decision_function物件f
//得到svm_model物件model

predict

//根據svm_type的不同開闢dec_value空間
svm_predict()
    //根據svm_type的不同調用k_function函式
    svm_predict_values()
        //根據kernel_type的不同計算k(i,j)
        Kernel::k_function()
        //得到k(x_train[i],x_test[j])
    //得到預測值y_pre[j]
//得到預測值y_pre[j]

svm.h檔案解析

svm_node

//儲存一個樣本（假設為樣本i）的一個特徵
struct svm_node{
    int index;   //樣本i的特徵索引值，最後一個為-1
    double value;   //樣本i第index個特徵的值，最後一個為NULL
};

如：x[i]={0.23,1.2,3.5,1.5}
則需使用五個svm_node來表示x[i]，具體對映如下：

svm_problem

//儲存參加運算的所有樣本資料（訓練集）
struct svm_problem{
        int l;    //樣本總數
        double *y;    //樣本輸出值（所屬類別）
        struct svm_node **x;    //樣本輸入值
};

下圖中最右邊的長條格同上表，儲存了三維資料。

**svm_problem中的y與類Solver中的y並不完全一樣!!!**對於一般SVC而言可以看出一樣的，其值為-1與+1，對於多分類而言svm_problem.y[i]可以是1、2、3等，而多類計算其實是二分類的組合，因此在二分類中y[i]依然等於+1與-1.更特殊的，在SVR中，svm_problem的y[i]等於其目標值，如：11.234、56.24、5.23等，在計算時svm_problem.y[i]整合到了Solver.p[i]與Solver.p[i+svm_problem.l]中（具體的問題後續章節再詳細解釋），而在Solver.y[i]依然為+1和-1.

svm_parameter

//svm_type和svm_type可能取值
enum { C\_SVC, NU\_SVC, ONE\_CLASS, EPSILON\_SVR, NU\_SVR };/* svm_type */
enum { LINEAR, POLY, RBF, SIGMOID }; /* kernel_type */

//svm模型訓練引數
struct svm_parameter
    {
        int svm_type;
        int kernel_type;
        int degree; /* for poly */
        double gamma;   /* for poly/rbf/sigmoid */
        double coef0;   /* for poly/sigmoid */

        /* these are for training only */
        double cache_size; /* in MB */
        double eps; /* stopping criteria */
        double C;   /* for C_SVC, EPSILON_SVR and NU_SVR */
        int nr_weight;      /* for C_SVC */
        int *weight_label;  /* for C_SVC */
        double* weight;     /* for C_SVC */
        double nu;  /* for NU_SVC, ONE_CLASS, and NU_SVR */
        double p;   /* for EPSILON_SVR */
        int shrinking;  /* use the shrinking heuristics */
        int probability; /* do probability estimates */
    };

LibSVM中的核函式如下：

各引數解釋如下：

svm_model

//儲存訓練後的模型引數
struct svm_model{
        struct svm_parameter param; /* parameter */
        int nr_class;       /* number of classes, = 2 in regression/one class svm */
        int l;          /* total #SV */
        struct svm_node **SV;       /* SVs (SV[l]) */
        double **sv_coef;   /* coefficients for SVs in decision functions (sv_coef[k-1][l]) */
        double *rho;        /* constants in decision functions (rho[k*(k-1)/2]) */
        double *probA;      /* pariwise probability information */
        double *probB;
        int *sv_indices;        /* sv_indices[0,...,nSV-1] are values in [1,...,num_traning_data] to indicate SVs in the training set */

                                /* for classification only */

        int *label;     /* label of each class (label[k]) */
        int *nSV;       /* number of SVs for each class (nSV[k]) */
                        /* nSV[0] + nSV[1] + ... + nSV[k-1] = l */
                        /* XXX */
        int free_sv;        /* 1 if svm_model is created by svm_load_model*/
                            /* 0 if svm_model is created by svm_train */
};

各引數解釋如下：

svm.cpp檔案解析

下圖為svm.cpp中的類繼承和組合情況（實現表示繼承關係，虛線表示組合關係）：

Cache類主要負責運算所涉及的記憶體的管理，包括申請、釋放等。本篇博文主要講解SVM求解過程，對於Cache類將不予解析。

Kernel類

class Kernel : public QMatrix {
public:
    Kernel(int l, svm_node * const * x, const svm_parameter& param);
    virtual ~Kernel();

    static double k_function(const svm_node *x, const svm_node *y,
        const svm_parameter& param);
    virtual Qfloat *get_Q(int column, int len) const = 0;
    virtual double *get_QD() const = 0;
    virtual void swap_index(int i, int j) const // no so const...
    {
        swap(x[i], x[j]);
        if (x_square) swap(x_square[i], x_square[j]);
    }
protected:

    double (Kernel::*kernel_function)(int i, int j) const;

private:
    const svm_node **x;
    double *x_square;

    // svm_parameter
    const int kernel_type;
    const int degree;
    const double gamma;
    const double coef0;

    static double dot(const svm_node *px, const svm_node *py);
    double kernel_linear(int i, int j) const
    {
        return dot(x[i], x[j]);
    }
    double kernel_poly(int i, int j) const
    {
        return powi(gamma*dot(x[i], x[j]) + coef0, degree);
    }
    double kernel_rbf(int i, int j) const
    {
        return exp(-gamma * (x_square[i] + x_square[j] - 2 * dot(x[i], x[j])));
    }
    double kernel_sigmoid(int i, int j) const
    {
        return tanh(gamma*dot(x[i], x[j]) + coef0);
    }
    double kernel_precomputed(int i, int j) const
    {
        return x[i][(int)(x[j][0].value)].value;
    }
};

成員變數

成員函式

Kernel(int l, svm_node * const * x, const svm_parameter& param);

建構函式。初始化類中的部分常量、指定核函式、克隆樣本資料。

Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param)
    :kernel_type(param.kernel_type), degree(param.degree),
    gamma(param.gamma), coef0(param.coef0)
{
    switch (kernel_type)    //根據kernel_type的不同定義相應的函式kernel_function()
    {
    case LINEAR:
        kernel_function = &Kernel::kernel_linear;
        break;
    case POLY:
        kernel_function = &Kernel::kernel_poly;
        break;
    case RBF:
        kernel_function = &Kernel::kernel_rbf;
        break;
    case SIGMOID:
        kernel_function = &Kernel::kernel_sigmoid;
        break;
    case PRECOMPUTED:
        kernel_function = &Kernel::kernel_precomputed;
        break;
    }

    clone(x, x_, l);

    if (kernel_type == RBF)    //如果使用RBF 核函式，則計算x_sqare[i]，即x[i]^T*x[i]
    {
        x_square = new double[l];
        for (int i = 0; i<l; i++)
            x_square[i] = dot(x[i], x[i]);
    }
    else
        x_square = 0;
}

static double dot(const svm_node *px, const svm_node *py);

點乘函式，點乘兩個樣本資料，按svm_node 中index (一般為特徵)進行運算，一般來說，index為1，2，…直到-1。返回點乘總和。例如：x1={1,2,3} ,x2={4,5,6}總和為sum=1*4+2*5+3*6;在svm_node[3]中儲存index=-1時，停止計算。

double Kernel::dot(const svm_node *px, const svm_node *py)
{
    double sum = 0;
    while (px->index != -1 && py->index != -1)
    {
        if (px->index == py->index)
        {
            sum += px->value * py->value;
            ++px;
            ++py;
        }
        else
        {
            if (px->index > py->index)
                ++py;
            else
                ++px;
        }
    }
    return sum;
}

static double k_function(const svm_node *x, const svm_node *y, const svm_parameter& param);

功能類似kernel_function,不過kerel_function用於訓練，k_function用於預測。

double Kernel::k_function(const svm_node *x, const svm_node *y,
    const svm_parameter& param)    //輸入資料為兩個資料樣本，其中一個為訓練樣本一個為測試樣本
{
    switch (param.kernel_type)
    {
    case LINEAR:
        return dot(x, y);
    case POLY:
        return powi(param.gamma*dot(x, y) + param.coef0, param.degree);
    case RBF:
    {
        double sum = 0;
        while (x->index != -1 && y->index != -1)
        {
            if (x->index == y->index)
            {
                double d = x->value - y->value;
                sum += d * d;
                ++x;
                ++y;
            }
            else
            {
                if (x->index > y->index)
                {
                    sum += y->value * y->value;
                    ++y;
                }
                else
                {
                    sum += x->value * x->value;
                    ++x;
                }
            }
        }

        while (x->index != -1)
        {
            sum += x->value * x->value;
            ++x;
        }

        while (y->index != -1)
        {
            sum += y->value * y->value;
            ++y;
        }

        return exp(-param.gamma*sum);
    }
    case SIGMOID:
        return tanh(param.gamma*dot(x, y) + param.coef0);
    case PRECOMPUTED:  //x: test (validation), y: SV
        return x[(int)(y->value)].value;
    default:
        return 0;  // Unreachable
    }
}

其中RBF部分很有講究。因為儲存時，0值不保留。如果所有0值都保留，第一個while就可以都做完了；如果第一個while做不完，在x，y中任意一個出現index＝-1，第一個while就停止，剩下的程式碼中兩個while只會有一個工作，該迴圈直接把剩下的計算做完。

virtual Qfloat *get_Q(int column, int len);

純虛擬函式，將來在子類中實現(如class SVR_Q)，計算Q值。相當重要的函式。

virtual Qfloat *get_Q(int column, int len) const = 0;

virtual void swap_index(int i, int j);

虛擬函式，x[i]和x[j]中所儲存指標的內容。如果x_square不為空，則交換相應的內容。

virtual void swap_index(int i, int j) const // no so const...
    {
        swap(x[i], x[j]);
        if (x_square) swap(x_square[i], x_square[j]);
    }

virtual double *get_QD();

純虛擬函式，將來在子類中實現(如class SVR_Q),計算Q[i,i]值。

virtual Qfloat *get_Q(int column, int len) const = 0;

double (Kernel::*kernel_function)(int i, int j)；

函式指標，根據相應的核函式型別，來決定所使用的函式。在計算矩陣Q時使用。

double (Kernel::*kernel_function)(int i, int j) const;

Solver類

class Solver {
public:
    Solver() {};
    virtual ~Solver() {};

    struct SolutionInfo {
        double obj;
        double rho;
        double upper_bound_p;
        double upper_bound_n;
        double r;   // for Solver_NU
    };

    void Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
        double *alpha_, double Cp, double Cn, double eps,
        SolutionInfo* si, int shrinking);
protected:
    int active_size;
    schar *y;
    double *G;      // gradient of objective function
    enum { LOWER_BOUND, UPPER_BOUND, FREE };
    char *alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE
    double *alpha;
    const QMatrix *Q;
    const double *QD;
    double eps;
    double Cp, Cn;
    double *p;
    int *active_set;
    double *G_bar;      // gradient, if we treat free variables as 0
    int l;
    bool unshrink;  // XXX

    double get_C(int i)
    {
        return (y[i] > 0) ? Cp : Cn;
    }
    void update_alpha_status(int i)
    {
        if (alpha[i] >= get_C(i))
            alpha_status[i] = UPPER_BOUND;
        else if (alpha[i] <= 0)
            alpha_status[i] = LOWER_BOUND;
        else alpha_status[i] = FREE;
    }
    bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
    bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
    bool is_free(int i) { return alpha_status[i] == FREE; }
    void swap_index(int i, int j);
    void reconstruct_gradient();
    virtual int select_working_set(int &i, int &j);
    virtual double calculate_rho();
    virtual void do_shrinking();
private:
    bool be_shrunk(int i, double Gmax1, double Gmax2);
};

成員變數

結構體SolutionInfo為求解優化中的引數資訊。

各引數解釋如下：