SVM之序列最小最优化算法（SMO算法）

Leave a reply

SVM回顾

支持向量机（SVM）的一大特点是最大化间距（max margin）。对于如上图的二分类问题，虽然有很多线可以将左右两部分分开，但是只有中间的红线效果是最好的，因为它的可活动范围（margin）是最大的，从直观上来说很好理解。

对于线性二分类问题，假设分类面为

\begin{equation} u=\vec w \cdot \vec x-b \end{equation}

则margin为

\begin{equation} m=\frac{1}{||w||_2} \end{equation}

根据max margin规则和约束条件，得到如下优化问题，我们要求的就是参数\vec w和b：

\begin{equation} \min\limits_{\vec w,b}\frac{1}{2}||\vec w||^2 \quad\text{subject to}\quad y_i(\vec w\cdot \vec x_i-b) \geq 1, \forall i,\end{equation}

对于正样本，类标号y_i为+1，反之则为-1。根据拉格朗日对偶，(3)可以转换为如下的二次规划（QP）问题，其中\alpha_i为拉格朗日乘子。

\begin{equation} \min\limits_{\vec \alpha}\Psi(\vec\alpha)=\min\limits_{\vec \alpha}\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^Ny_iy_j(\vec x_i\cdot\vec x_j)\alpha_i\alpha_j-\sum_{i=1}^N\alpha_i,\end{equation}

其中N为样本数量。上式还需满足如下两个约束条件：

\begin{equation} \alpha_i\geq 0, \forall i,\end{equation}

\begin{equation} \sum_{i=1}^Ny_i\alpha_i=0.\end{equation}

一旦求解出所有的拉格朗日乘子，则我们可以通过如下的公式得到分类面参数\vec w和b。

\begin{equation}\vec w=\sum_{i=1}^Ny_i\alpha_i\vec x_i,\quad b=\vec w\cdot\vec x_k-y_k\quad\text{for some}\quad\alpha_k>0.\end{equation}

当然并不是所有的数据都可以完美的线性划分，可能有少量数据就是混在对方阵营，这时可以通过引入松弛变量\xi_i得到软间隔形式的SVM：

\begin{equation} \min\limits_{\vec w,b,\vec\xi}\frac{1}{2}||\vec w||^2+C\sum_{i=1}^N\xi_i \quad\text{subject to}\quad y_i(\vec w\cdot \vec x_i-b) \geq 1-\xi_i, \forall i,\end{equation}

其中的\xi_i为松弛变量，能假装把错的样本分对，C对max margin和margin failures的trades off。对于这个新的优化问题，约束变成了一个box constraint：

\begin{equation}0\leq\alpha_i \leq C,\forall i.\end{equation}

而松弛变量\xi_i不再出现在对偶公式中了。

对于线性不可分的数据，可以用核函数K将其投影到高维空间，这样就可分了，由此得到一般的分类面公式：

\begin{equation}u=\sum_{j=1}^Ny_j\alpha_jK(\vec x_j,\vec x)-b,\end{equation}

终极优化问题就变成了下面这个样子：

\begin{equation} \min\limits_{\vec \alpha}\Psi(\vec\alpha)=\min\limits_{\vec \alpha}\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^Ny_iy_jK(\vec x_i,\vec x_j)\alpha_i\alpha_j-\sum_{i=1}^N\alpha_i,\\0\leq\alpha_i \leq C,\forall i,\\ \sum_{i=1}^Ny_i\alpha_i=0.\end{equation}

要满足的KKT条件为：

\begin{equation}\alpha_i=0\Leftrightarrow y_iu_i\geq1,\\0<\alpha_i<C\Leftrightarrow y_iu_i=1,\\ \alpha_i=C\Leftrightarrow y_iu_i\leq 1.\end{equation}

SMO算法

为了求解式(11)，SMO算法通过启发式方法选择两个\alpha_i,\alpha_j当变量，固定其他\alpha_k，然后用解析的方法求解两个变量的二次规划问题。关于这两个变量的解应该更接近原始的二次规划问题，更重要的是，这时子问题可以通过解析方法求解，避免了矩阵运算，大大提高了整个算法的计算速度。如此，SMO算法将原问题不断分解为子问题并对子问题求解，进而达到求解原问题的目的。

不失一般性，假设选择的两个拉格朗日乘子是\alpha_1,\alpha_2，因为\alpha_i和变量x_i是一一对应的，所以也可以说选择了变量x_1,x_2。固定其他变量\alpha_i (i=3,4,...,N)

此时，式(11)的优化问题转换为如下的优化问题：

\begin{equation} \min\limits_{\alpha_1,\alpha_2}W(\alpha_1,\alpha_2)=\frac{1}{2}K_{11}\alpha_1^2+\frac{1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1\alpha_2\\-(\alpha_1+\alpha_2)+y_1\alpha_1\sum_{i=3}^Ny_i\alpha_iK_{i1}+y_2\alpha_2\sum_{i=3}^Ny_i\alpha_iK_{i2}\end{equation}

满足

\begin{equation}\alpha_1y_1+\alpha_2y_2=-\sum_{i=3}^Ny_i\alpha_i=\zeta\end{equation}

因为只有\alpha_1,\alpha_2是变量，其他\alpha_i (i=3,4,...,N)是固定的（可以认为是常数），所以式(13)省略了不含\alpha_1,\alpha_2的常数项。又因为y_i=\pm 1，\alpha_i有限制(9)，所以\alpha_1,\alpha_2被限制在[0,C]\times[0,C]的盒子里，且位于对角线上，如图Figure 1.

SMO-1

由于\alpha_2^{new}需满足(9)，所以最优值\alpha_2^{new}的取值范围必须满足条件

\begin{equation}L\leq\alpha_2^{new}\leq H\end{equation}

其中L与H是\alpha_2^{new}所在的对角线段端点的界。如果y_1\neq y_2（Figure 1.左图），则

\begin{equation}L=max(0,\alpha_2^{old}-\alpha_1^{old}),\quad H=min(C,C+\alpha_2^{old}-\alpha_1^{old}).\end{equation}

如果y_1= y_2（Figure 1.右图），则

\begin{equation}L=max(0,\alpha_2^{old}+\alpha_1^{old}-C),\quad H=min(C,\alpha_2^{old}+\alpha_1^{old}).\end{equation}

令

\begin{equation}\eta=K(\vec x_1,\vec x_1)+K(\vec x_2,\vec x_2)-2K(\vec x_1,\vec x_2).\end{equation}

因为y_i=\pm 1，所以y_i^2=1，式(14)同乘以y_1，用\alpha_2表示\alpha_1并带入式(13)，得到只含一个变量\alpha_2的表达式，该表达式对\alpha_2求偏导并等于0，求出\alpha_2的更新公式为：

\begin{equation}\alpha_2^{new}=\alpha_2^{old}+\frac{y_2(E_1-E_2)}{\eta},\end{equation}

其中E_i=u_i-y_i是第i个样本的训练误差。经过截断之后的\alpha_2^{new}为：

\begin{equation}\alpha_2^{new,clipped}=\begin{cases} H, & \mbox{if }\alpha_2^{new}\geq H; \\ \alpha_2^{new}, & \mbox{if }L<\alpha_2^{new}<H;\\ L, & \mbox{if }\alpha_2^{new}\leq L.\end{cases}\end{equation}

将\alpha_2^{new,clipped}的结果代入式(14)，得到\alpha_1的更新公式：

\begin{equation}\alpha_1^{new}=\alpha_1^{old}+y_1y_2(\alpha_2^{old}-\alpha_2^{new,clipped}).\end{equation}

当0<\alpha_1^{new}<C时，由式(12)中间的条件可以得到阈值b_1的更新公式：

\begin{equation}b_1^{new}=b^{old}-E_1-y_1K_{11}(\alpha_1^{new}-\alpha_1^{old})-y_2K_{21}(\alpha_2^{new,clipped}-\alpha_2^{old}),\end{equation}

同理，当0<\alpha_2^{new,clipped}<C时，得到阈值b_2的更新公式：

\begin{equation}b_2^{new}=b^{old}-E_2-y_1K_{12}(\alpha_1^{new}-\alpha_1^{old})-y_2K_{22}(\alpha_2^{new,clipped}-\alpha_2^{old}).\end{equation}

综合以上，阈值b的更新公式为：

\begin{equation}b^{new}=\begin{cases} b_1^{new}, & \mbox{if }0<\alpha_1^{new}<C; \\ b_2^{new}, & \mbox{if }0<\alpha_2^{new,clipped}<C;\\ (b_1^{new}+b_2^{new})/2, & \mbox{otherwise}.\end{cases}\end{equation}

至此，如果我们选定了两个变量\alpha_2,\alpha_1，则可以通过公式(19)和(21)来更新它们，并通过公式(24)更新阈值b，通过多次迭代逼近最优结果。但是怎样选择\alpha_2,\alpha_1呢，SMO通过启发式方法选择！

对于第一个样本\alpha_2，我们选择违反KKT条件式(12)（最严重）的样本\alpha_2；由式(19)可知，\alpha_2^{new}是依赖于|E_1-E_2|的，为了加快计算，一种简单的做法是选择\alpha_1，使其对应的|E_1-E_2|最大。

至此，我们有了启发式选择两个变量\alpha_2,\alpha_1的方法，以及求解这两个变量二次规划的解析方法，下面介绍SMO算法的具体实现。

MATLAB代码实现

算法伪代码如下：

输入：训练数据T=\{(\vec x_1,y_1),(\vec x_2,y_2),...,(\vec x_N,y_N)\}，其中\vec x_i \in R^n，y_i\in\{-1,+1\}，i=1,2,...,N，精度为\epsilon。

输出：拉格朗日乘子\vec\alpha=\{\alpha_1,\alpha_2,...,\alpha_N\}和阈值b的近似解。

  1. 初始化\vec \alpha=\vec 0, b=0
  2. 选一个违反KKT条件的变量\alpha_2
  3. 选一个使得|E_1-E_2|最大的变量\alpha_1
  4. 根据公式(19)~(21)更新变量\alpha_2,\alpha_1
  5. 根据公式(22)~(24)更新阈值b
  6. 如果不满足结束条件，则转2；否则算法结束

算法的结束条件是：如果所有变量都已经检查过，且没有变量可以进一步优化时，算法结束。

下面是SMO算法的MATLAB简化实现，省略了论文[1]的公式(19)。代码主流程参考论文[1]中的伪代码，部分实现细节参考[2]。

本代码和MATLAB自带的seqminopt.m算法接口相同，可直接替换使用。下载代码。

function [alphas, offset] = bitjoy_seqminopt(data, targetLabels, boxConstraints, ...
    kernelFunc, smoOptions)

% http://bitjoy.net/2016/05/02/svm-smo-algorithm/
% 参考文献:
% [1] A Fast Algorithm for Training  Support Vector Machines,
%     http://research.microsoft.com/pubs/69644/tr-98-14.pdf
% [2] CSDN, http://blog.csdn.net/techq/article/details/6171688

N = length(data); % 样本数量
alphas = zeros([N,1]); % 所有拉格朗日乘子初始化为0
offset = 0.0; % 偏移量b也初始化为0

numChanged = 0; % 拉格朗日乘子改变的个数
examineAll = 1; % 是否需要检查所有拉格朗日乘子

% 主流程
while numChanged > 0 || examineAll
    numChanged = 0;
    if examineAll == 1
        for i = 1 : N
            numChanged = numChanged + examineExample(i);
        end
    else
        for i = 1 : N
            if alphas(i) ~= 0 && alphas(i) ~= boxConstraints(i)
                numChanged = numChanged + examineExample(i);
            end
        end
    end
    if examineAll == 1
        examineAll = 0;
    elseif numChanged == 0
        examineAll = 1;
    end
end

    % 计算当前参数下第idx个样本的函数输出
    function [svm_o] = wxpb(idx)
        svm_o = 0.0;
        for j = 1 : N
            svm_o = svm_o + alphas(j) * targetLabels(j) * kernelFunc(data(j,:),data(idx,:));
        end
        svm_o = svm_o + offset;
    end

    % 选定第二个变量i2后，根据max|E1-E2|的启发式规则，选择i1；
    % 如果没有满足条件的i1，返回-1.
    function [i1] = selectSecondChoice(i2, E2)
        i1 = -1;
        maxDelta = -1;
        for j = 1 : N
            if j ~= i2 && alphas(j) ~= 0 && alphas(j) ~= boxConstraints(j)
                Ej = wxpb(j) - targetLabels(j);
                if abs(E2 - Ej) > maxDelta
                    i1 = j;
                    maxDelta = abs(E2 - Ej);
                end
            end
        end
    end

    % 根据选定的两个变量i1,i2，代入更新公式计算；
    % 最后还更新了偏移量offset，也就是y=wx+b中的b。
    function [flag] = takeStep(i1, i2)
        alpha1 = alphas(i1);
        y1 =  targetLabels(i1);
        E1 = wxpb(i1) - y1;
        alpha2 = alphas(i2);
        y2 =  targetLabels(i2);
        E2 = wxpb(i2) - y2;
        s = y1 * y2;
        if y1 ~= y2
            L = max(0, alpha2 - alpha1);
            H = min(boxConstraints(i2), boxConstraints(i1) + alpha2 - alpha1);
        else
            L = max(0, alpha2 + alpha1 - boxConstraints(i1));
            H = min(boxConstraints(i2), alpha2 + alpha1);
        end
        if L == H
            flag = 0;
            return;
        end
        k11 = kernelFunc(data(i1,:),data(i1,:));
        k12 = kernelFunc(data(i1,:),data(i2,:));
        k22 = kernelFunc(data(i2,:),data(i2,:));
        eta = k11 + k22 - 2 * k12;
        if eta > 0
            a2 = alpha2 + y2 * (E1 - E2) / eta;
            %截断
            if a2 < L
                a2 = L;
            elseif a2 > H
                a2 = H;
            end
        else % 并未处理该case，论文[1]公式(19)有更详细的方法
            flag = 0;
            return;
        end
        if abs(a2 - alpha2) < eps * (a2 + alpha2 + eps)
            flag = 0;
            return;
        end
        a1 = alpha1 + s * (alpha2 - a2);
        alphas(i1) = a1;
        alphas(i2) = a2;
        
        % 更新offset
        b1 = offset - E1 - y1 * k11 * (a1 - alpha1) - y2 * k12 * (a2 - alpha2);
        b2 = offset - E2 - y1 * k12 * (a1 - alpha1) - y2 * k22 * (a2 - alpha2);
        if a1 > 0 && a1 < boxConstraints(i1)
            offset = b1;
        elseif a2 > 0 && a2 < boxConstraints(i2)
            offset = b2;
        else
            offset = (b1 + b2) / 2;
        end
        flag = 1;
    end

    % 检查样本。
    % 首先判断i2是否满足KKT条件，如果不满足，
    % 则根据启发式规则再选择i1样本，
    % 然后更新i1和i2的拉格朗日乘子。
    function [flag] = examineExample(i2)
        y2 =  targetLabels(i2);
        alpha2 = alphas(i2);
        E2 = wxpb(i2) - y2;
        r2 = E2 * y2;
        if (r2 < -smoOptions.TolKKT && alpha2 < boxConstraints(i2)) || (r2 > smoOptions.TolKKT && alpha2 > 0)
            i1 = selectSecondChoice(i2, E2);
            if i1 == -1
                i1 = floor(1 + rand() * N); % 随机选一个i1
                while i1 == i2
                    i1 = floor(1 + rand() * N);
                end
                flag = takeStep(i1,i2);
            else
                flag = takeStep(i1,i2);
            end
        else
            flag = 0;
        end
    end
end

参考资料

[1]. J.C. Platt: A Fast Algorithm for Training Support Vector Machines

[2]. CSDN, techq, SVM算法实现（一）

[3]. 国科大叶齐祥老师机器学习方法与应用课程资料

[4]. 支持向量机（SVM）的详细推导过程及注解（一）

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.