机器学习 SVM

libSVM源码解读（4）

Solve类

Posted by Welt Xing on July 15, 2021

Solver类引入

我们这里会研究libSVM中的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：一个存储优化问题的解对应参数的结构体：

    struct SolutionInfo {
    double obj; // 目标函数最优值
    double rho; // 对应ρ
    double upper_bound_p;
    double upper_bound_n;
    double r;  // for Solver_NU
 };

  • active_size：实际参加矩阵运算的样本数，我们在这里详细介绍了关于active_size变量的含义；

  • y：类别向量，$y\in{-1,+1}$；

  • G：优化函数的梯度向量；

  • 枚举变量enum { LOWER_BOUND, UPPER_BOUND, FREE }：分别表示对应的$\alpha_i=0$，$\alpha_i=C$和$\alpha_i\in(0,C)$；

  • alpha_status：上面对应枚举值的数组；

  • alpha：也就是$\alpha$向量；

  • Q：也就是前面提到的$Q$矩阵，核函数和Solver相互结合，可以产生多种SVC,SVR；

  • QD：Q的对角元素构成的数组；

  • eps：判断迭代是否终止的阈值$\epsilon$；

  • Cp, Cn，表示Positive Class和Negative Class，用于多分类问题；

  • p：对应优化目标函数$f=\dfrac12\alpha^\top Q\alpha+p^\top\alpha$中的向量$p$；

  • active_set：也就是参加矩阵运算的样本数的序列集合，长度为active_size；

  • G_bar：对应论文中的$\bar G$，详细解析笔者写在这里；

  • l：样本数；

  • unshrinked：表示是否还没使用Shrinking技巧.

然后是成员函数，构造和析构函数跳过：

  • get_C：返回对应样本的种类；设置不同的Cp和Cn是为了处理数据的不平衡，类似于“再缩放”（?）；
  • update_alpha_status：用于在每轮迭代后对$\alpha_i$的状态更新，状态对应前面的枚举值；
  • 三个is_*函数，用于判断$\alpha_i$状态的封装函数；
  • swap_index：交换两样本产生的内容；
  • reconstruct_gradient：梯度的增量式更新，数学推导和代码讲解放在这里；
  • select_working_set：选择工作集，数学推导和代码讲解放在这里；
  • calculate_rho：计算参数$\rho$，相应的也就是在计算$b$；
  • do_shrinking：采取收缩算法，数学推导放在这里；
  • be_shrunk：判断参数$\alpha_i$是否收敛到最优；
  • Solve：核心函数，用于求解问题。

我们接下来对上面一些复杂或没有提到的函数进行解读。

$\rho$的计算

可以证明，$C-$SVC和$\varepsilon$-SVR中的$b$与One-class SVM中的$\rho$等价，且满足

\[b+\rho=0\]

对于一个支持向量$\pmb x_i$，当我们得到最优的$\pmb\alpha$后，有

\[b=y_i-\sum_{j=1}^ly_i\alpha_jK(\pmb x_j,\pmb x_i)\]

采用更鲁棒的方法，对所有支持向量对应的$b$取平均值：

\[\begin{aligned} b^*&=\dfrac{1}{\vert S\vert }\sum_{i\in S}(y_i-\sum_{j=1}^ly_i\alpha_jK(\pmb x_j,\pmb x_i))\\ &=\dfrac{1}{\vert S\vert}\sum_{i\in S}\big(-y_i\nabla f(\alpha)_i\big)\\ \end{aligned}\]

其中$S$为支持向量的下标集。由$b$与$\rho$的关系，我们也可得到

\[\rho=\dfrac{1}{\vert S\vert}\sum_{i\in S}y_i\nabla f(\alpha)_i\]

此时应该有

\[-M(\alpha)\leqslant\rho\leqslant-m(\alpha)\]

而不幸的是如果$S=\emptyset$，那么我们令

\[r_1=\dfrac{1}{2}\bigg(\max_{\alpha_i=C,y_i=1}\nabla f(\alpha)_i+\min_{\alpha_i=0,y_i=1}\nabla f(\alpha)_i\bigg)\\ r_2=\dfrac{1}{2}\bigg(\max_{\alpha_i=C,y_i=-1}\nabla f(\alpha)_i+\min_{\alpha_i=0,y_i=-1}\nabla f(\alpha)_i\bigg)\]

从而设

\[\rho=\dfrac{r_1+r_2}{2},-b=\dfrac{r_1-r_2}{2}\]

我们来看calculate_rho的代码：

double Solver::calculate_rho() {
    double r;
    int nr_free = 0;
    double ub = INF, lb = -INF, sum_free = 0;
    for (int i = 0; i < active_size; i++) {
        double yG = y[i] * G[i];

        if (is_upper_bound(i)) {
            if (y[i] == -1)
                ub = min(ub, yG);
            else
                lb = max(lb, yG);
        } else if (is_lower_bound(i)) {
            if (y[i] == +1)
                ub = min(ub, yG);
            else
                lb = max(lb, yG);
        } else {
            ++nr_free;
            sum_free += yG;
        }
    }

    if (nr_free > 0)
        r = sum_free / nr_free;
    else
        r = (ub + lb) / 2;

    return r;
}

代码的行为和算法一致：遍历active_set，同时计算sum_free和约束上下界，为两种情况做准备。

Shrinking相关

我们直接看源码：

void Solver::do_shrinking() {
    int i;
    double Gmax1 = -INF;  // max { -y_i * grad(f)_i | i in I_up(\alpha) }
    double Gmax2 = -INF;  // max { y_i * grad(f)_i | i in I_low(\alpha) }

    // find maximal violating pair first
    for (i = 0; i < active_size; i++) {
        if (y[i] == +1) {
            if (!is_upper_bound(i)) {
                if (-G[i] >= Gmax1) Gmax1 = -G[i];
            }
            if (!is_lower_bound(i)) {
                if (G[i] >= Gmax2) Gmax2 = G[i];
            }
        } else {
            if (!is_upper_bound(i)) {
                if (-G[i] >= Gmax2) Gmax2 = -G[i];
            }
            if (!is_lower_bound(i)) {
                if (G[i] >= Gmax1) Gmax1 = G[i];
            }
        }
    }
	...
}

注意到这里是找最大违反对，其实是求解$m(\alpha)$和$M(\alpha)$。

	...
    if (unshrink == false && Gmax1 + Gmax2 <= eps * 10) {
        unshrink = true;
        reconstruct_gradient();
        active_size = l;
        info("*");
    }
	...

这一部分在论文里也有提及：为了防止收缩策略幅度太大，于是利用收敛阈值$\varepsilon$对其限制：如果

\[m(\alpha)\leqslant M(\alpha^k)+10\varepsilon\]

就进行收缩。最后一个for循环则是将所有已经被收缩的变量放到后面，为访问提供便利，在论文中这一操作被称为“Index Rearrangement”：

	...
	for (i = 0; i < active_size; i++)
        if (be_shrunk(i, Gmax1, Gmax2)) {
            active_size--;
            while (active_size > i) {
                if (!be_shrunk(active_size, Gmax1, Gmax2)) {
                    swap_index(i, active_size);
                    break;
                }
                active_size--;
            }
        }
}

Solve函数

Solve函数是Solver类甚至是整个libSVM的核心函数，负责求解二次规划问题。我们先来看函数声明：

void Solver::Solve(int l, const QMatrix &Q, const double *p_, const schar *y_,
                   double *alpha_, double Cp, double Cn, double eps,
                   SolutionInfo *si, int shrinking);

函数参数都是我们前面提到过的，它们可以刻画出一个二次规划问题。我们将它的实现拆开来讲：

初始化

先是一些数据和参数的拷贝和赋值：

this->l = l;
this->Q = &Q;
QD = Q.get_QD();
clone(p, p_, l);
clone(y, y_, l);
clone(alpha, alpha_, l);
this->Cp = Cp;
this->Cn = Cn;
this->eps = eps;
unshrink = false;

然后初始化所有的$\alpha_i$都是active的，在之后的迭代中随着shrink而减少。

{
	active_set = new int[l];
	for (int i = 0; i < l; i++) 
        active_set[i] = i;
	active_size = l;
}

初始化梯度，包括G和G_bar：

{
    G = new double[l];
    G_bar = new double[l];
    int i;
    for (i = 0; i < l; i++) {
        G[i] = p[i];
        G_bar[i] = 0;
    }
    for (i = 0; i < l; i++)
        if (!is_lower_bound(i)) {
            const Qfloat *Q_i = Q.get_Q(i, l);
            double alpha_i = alpha[i];
            int j;
            for (j = 0; j < l; j++) 
                G[j] += alpha_i * Q_i[j];
            if (is_upper_bound(i))
                for (j = 0; j < l; j++) 
                    G_bar[j] += get_C(i) * Q_i[j];
        }
}

由

\[\begin{aligned} f(\alpha)&=\dfrac12\alpha^\top Q\alpha-p^\top\alpha\\ \nabla f(\alpha)_i&=\sum_{j=1}^lQ_{ij}\alpha_i-p_i \end{aligned}\]

在第一个循环中，对每个梯度分量先加上$p_i$；此外，由于所有变量都是active的，G_bar初始化为0。第二个循环显然是上式的求和过程。注意到为了减少无意义的循环，当$\alpha_i$为0时就不再计算；如果发现存在变成inactive的$\alpha_i$，就对G_bar进行更新；

优化

初始化完成后就可以对问题进行优化，为了防止不收敛，我们得确定一个最大循环次数：

int max_iter = max(10000000, l > INT_MAX / 100 ? INT_MAX : 100 * l);

但我目前并没有找到文献资料能说明这样的设计的原因是什么，姑且理解为一种启发式方法。

正式优化之前，libSVM会每隔一定的循环次数就进行shrink操作：

if (--counter == 0) {
	counter = min(l, 1000);
	if (shrinking) 
        do_shrinking();
	info(".");
}

然后是工作集选取：

int i, j;
if (select_working_set(i, j) != 0) {
	// reconstruct the whole gradient
	reconstruct_gradient();
	// reset active set size and check
    active_size = l;
    info("*");
	if (select_working_set(i, j) != 0)
		break;
	else
		counter = 1;  // do shrinking next iteration
}

如果选取工作集失败则会进入这个if语句，此时会重构梯度，将所有的$\alpha$都设置为active，然后重新选择，如果再一次选不出来，说明这个问题已经收敛了（?），退出迭代；否则就继续迭代，而且必然会进行一次Shrink，因为我们将所有拉格朗日乘子都激活了。

接着是对选定的$\alpha_i$和$\alpha_j$进行更新，这里涉及到更新和剪辑操作，要对$y_i=y_j$和$y_i\neq y_j$进行分类讨论：

const Qfloat *Q_i = Q.get_Q(i, active_size);
const Qfloat *Q_j = Q.get_Q(j, active_size);

double C_i = get_C(i);
double C_j = get_C(j);

double old_alpha_i = alpha[i];
double old_alpha_j = alpha[j];

if (y[i] != y[j]) {
    double quad_coef = QD[i] + QD[j] + 2 * Q_i[j];
    if (quad_coef <= 0) quad_coef = TAU;
    double delta = (-G[i] - G[j]) / quad_coef;
    double diff = alpha[i] - alpha[j];
    alpha[i] += delta;
    alpha[j] += delta;

    if (diff > 0) {
        if (alpha[j] < 0) {
            alpha[j] = 0;
            alpha[i] = diff;
        }
    } else {
        if (alpha[i] < 0) {
            alpha[i] = 0;
            alpha[j] = -diff;
        }
    }
    if (diff > C_i - C_j) {
        if (alpha[i] > C_i) {
            alpha[i] = C_i;
            alpha[j] = C_i - diff;
        }
    } else {
        if (alpha[j] > C_j) {
            alpha[j] = C_j;
            alpha[i] = C_j + diff;
        }
    }
} else {
    double quad_coef = QD[i] + QD[j] - 2 * Q_i[j];
    if (quad_coef <= 0) quad_coef = TAU;
    double delta = (G[i] - G[j]) / quad_coef;
    double sum = alpha[i] + alpha[j];
    alpha[i] -= delta;
    alpha[j] += delta;

    if (sum > C_i) {
        if (alpha[i] > C_i) {
            alpha[i] = C_i;
            alpha[j] = sum - C_i;
        }
    } else {
        if (alpha[j] < 0) {
            alpha[j] = 0;
            alpha[i] = sum;
        }
    }
    if (sum > C_j) {
        if (alpha[j] > C_j) {
            alpha[j] = C_j;
            alpha[i] = sum - C_j;
        }
    } else {
        if (alpha[i] < 0) {
            alpha[i] = 0;
            alpha[j] = sum;
        }
    }
}

具体的算法细节被放置在这里。

辅助变量的更新

优化完成后需要对辅助变量G和G_bar进行更新，准备下一次迭代，我们在梯度重构中讲解了这部分代码。

迭代超限的处理

如果到了指定迭代次数但仍未收敛，我们会重构一次梯度以计算此时的优化目标函数值，同时输出警告信息：

if (iter >= max_iter) {
	if (active_size < l) {
	// reconstruct the whole gradient to calculate objective value
	reconstruct_gradient();
	active_size = l;
	info("*");
	}
	fprintf(stderr, "\nWARNING: reaching max number of iterations\n");
}

收尾工作

在迭代结束后，便是最优解和参数的处理工作：

// calculate rho
si->rho = calculate_rho();

// calculate objective value
{
	double v = 0;
	int i;
	for (i = 0; i < l; i++) 
        v += alpha[i] * (G[i] + p[i]);
    
	si->obj = v / 2;
}

// put back the solution
{
	for (int i = 0; i < l; i++) 
        alpha_[active_set[i]] = alpha[i];
}

包括$\rho$，目标函数值的代入，以及参数解的赋值。别忘了申请空间的释放：

delete[] p;
delete[] y;
delete[] alpha;
delete[] alpha_status;
delete[] active_set;
delete[] G;
delete[] G_bar;

至此，我们对完成了libSVM的求解器的分析。

CATALOG
    FEATURED TAGS
    C/C++ Java Python 机器学习 SVM 集成学习 优化算法 神经网络 文献解读 数学 深度学习 卷积神经网络 自动微分 异常检测 数字信号处理 LIBLINEAR 在线学习 图神经网络