解密SVM系列（三）：SMO算法原理与实战求解

上节我们讨论到解SVM问题最终演化为求下列带约束条件的问题：

minW(α)=12(∑i,j=1Nαiyiαjyjxi∗xj)−∑i=1Nαis.t.0≤αi≤C∑i=1Nαiyi=0min\quad W(\alpha)=\dfrac{1}{2}(\sum_{i,j=1}^{N}\alpha_iy_i\alpha_jy_jx_i*x_j)-\sum_{i=1}^{N}\alpha_i\\s.t. \quad 0\le \alpha_i\le C\\ \quad \quad \quad \quad \sum_{i=1}^{N}\alpha_iy_i=0

问题的解就是找到一组αi=(α1,α2,...,αn)\alpha_i = (\alpha_1,\alpha_2,...,\alpha_n)使得W最小。

现在我们来看看最初的约束条件是什么样子的：

1−yi(Wxi+b)≤0 \quad 1- y_i(Wx_i+b) \le 0

这是最初的一堆约束条件吧，现在有多少个约束条件就会有多少个αi\alpha_i。那么KKT条件的形成就是让αi(1−yi(Wxi+b))=0 \alpha_i(1- y_i(Wx_i+b) )=0。

我们知道αi≥0\alpha_i\ge0,而后面那个小于等于0，所以他们中间至少有一个为0（至于为什么要这么做，第一节讨论过）。再简单说说原因，假设现在的分类问题如下：




某一次迭代后，分类面为粗蓝线所示，上下距离为1的分界面如细蓝线所示，而理想的分界面如紫虚线所示。那么我们想想，要想把粗蓝线变化到紫虚线，在这一次是哪些个点在起作用？很显然是界于细蓝线边上以及它们之间的所有样本点在起作用吧，而对于那些在细蓝线之外的点，比如正类的四个圈和反类的三个叉，它们在这一次的分类中就已经分对了，那还考虑它们干什么？所以这一次就不用考虑这些分对了的点。那么我们用数学公式可以看到，对于在这一次就分对了的点，它们满足什么关系，显然yi(Wxi+b)>1y_i(Wx_i+b) > 1,然后还得满足αi(1−yi(Wxi+b))=0 \alpha_i(1- y_i(Wx_i+b) )=0，那么显然它们的αi=0\alpha_i=0。对于那些在边界内的点，显然yi(Wxi+b)≤1y_i(Wx_i+b) \le 1，而这些点我们说是要为下一次达到更好的解做贡献的，那么我们就取这些约束条件的极限情况，也就是yi(Wxi+b)=1y_i(Wx_i+b) = 1，在这些极限约束条件下，我们就会得到一组新的权值W与b，也就是改善后的解。那么既然这些点的yi(Wxi+b)=1y_i(Wx_i+b) = 1，那它对应的αi\alpha_i就可以不为0了，至于是多少，那就看这些点具体属于分界面内的什么位置了，偏离的越狠的点，我想它对应的αi\alpha_i就越大，这样才能把这个偏得非常狠的点给拉回来，或者说使其在下一次的解中更靠近正确的分类面。

好了这就是KKT为什么要这么做的原因。那么整理一下，最终带松弛变量的KKT条件就如下所示：

(1)αi=0⇒yiui≥1(分对了)(2)0≤αi≤C⇒yiui=1(边界上，支持向量)(3)αi=C⇒yiui≤1(边界之间)(1) \alpha_i=0\quad\Rightarrow\quad y_iu_i\ge1\quad (分对了)\\(2)0\le\alpha_i\le C\quad\Rightarrow\quad y_iu_i=1\quad (边界上，支持向量)\\(3)\alpha_i=C\quad\Rightarrow\quad y_iu_i\le1\quad (边界之间)

那么满足KKT条件的，我们说如果一个点满足KKT条件，那么它就不需要调整，一旦不满足，就需要调整。由上可知，不满足KKT条件的也有三种情况：

(1)yiui≤1但是αi<C(2)yiui≥1但是αi>0(3)yiui≤1但是αi=0或αi=C(1) y_iu_i\le1\quad但是\quad\alpha_i0 \\(3)y_iu_i\le1\quad 但是\quad \alpha_i=0或\alpha_i=C

这三种条件下的α\alpha都需要调整。那么怎么调整呢？比如说某一次迭代完后发现有10个点不满足，也就是有10个α\alpha需要调整，那么10个α\alpha在一起，你怎么知道去增加或者减少哪一个或者哪几个α\alpha呢？这个时候SMO采取的策略是选择这10个中的两个α\alpha出来，就假设为α1,α2\alpha_1,\alpha_2吧，调整它们。看看前面有一个条件∑Ni=1αiyi=0\sum_{i=1}^{N}\alpha_iy_i=0，所以在调整前后满足：

αnew1y1+αnew2y2=αold1y1+αold2y2=ϵ\alpha_1^{new}y_1+\alpha_2^{new}y_2=\alpha_1^{old}y_1+\alpha_2^{old}y_2=\epsilon

至于ϵ\epsilon是多少，管它多少，又不用它。也就是说我们只要知道了α1,α2\alpha_1,\alpha_2中一个调整以后的值，根据条件另一个值不用算就知道。那么怎么求呢？假设我们来求α2\alpha_2吧。观察上面那个等式，y1,y2y1,y2是标签，要么1要么-1。而两个α>=0\alpha>=0。所以新的α\alpha是有范围的。假设现在y1=y2=1或−1y1=y2=1或-1，先=1吧，那么αnew1+αnew2=αold1+αold2=ϵ\alpha_1^{new}+\alpha_2^{new}=\alpha_1^{old}+\alpha_2^{old}=\epsilon

因为αnew1\alpha_1^{new}是在0-C之间，所以假设αnew1=0\alpha_1^{new}=0,那么αnew2\alpha_2^{new}可以取到最大值为ϵ \epsilon,也就是αold1+αold2\alpha_1^{old}+\alpha_2^{old}。而αnew2\alpha_2^{new}又不能大于C，所以其最大取值为min(C,αold1+αold2)min(C,\alpha_1^{old}+\alpha_2^{old})。同理如果αnew1=C\alpha_1^{new}=C，那么αnew2\alpha_2^{new}可以取到最小值为ϵ−C \epsilon-C,也就是αold1+αold2−C\alpha_1^{old}+\alpha_2^{old}-C,而αnew2\alpha_2^{new}最小不能小于0，那玩意这个值比0小怎么办？所以αnew2\alpha_2^{new}的下限值就为max(0,αold1+αold2−C)max(0,\alpha_1^{old}+\alpha_2^{old}-C)。这是相等取1，那么相等取-1呢？同理吧，不过是

−αnew1−αnew2=−αold1−αold2=ϵ-\alpha_1^{new}-\alpha_2^{new}=-\alpha_1^{old}-\alpha_2^{old}=\epsilon两边乘以-1，就是αnew1+αnew2=αold1+αold2=−ϵ\alpha_1^{new}+\alpha_2^{new}=\alpha_1^{old}+\alpha_2^{old}=-\epsilon，因为你也不知道−ϵ-\epsilon是多少，不过一个字母而已，那么用ϵ\epsilon代替−ϵ-\epsilon吧，接下来的讨论过程一模一样。从而属于同类别的时候αnew2\alpha_2^{new}上下限就有了。同理去计算下不同类（1，-1）的时候的上下限。最终也能得到一个类似的结果，这里就省了，给出最后的结果：

ify1≠y2L=max(0,αold2−αold1),H=min(C,C+αold2−αold1)ify1=y2L=max(0,αold2+αold1−C),H=min(C,αold2+αold1)if\quad y_1\neq y_2\quad L=max(0,\alpha_2^{old}-\alpha_1^{old}),H=min(C,C+\alpha_2^{old}-\alpha_1^{old})\\if\quad y_1=y_2\quad L=max(0,\alpha_2^{old}+\alpha_1^{old}-C),H=min(C,\alpha_2^{old}+\alpha_1^{old})

到这只是给出了αnew2\alpha_2^{new}的范围，那么它怎么解呢？解这个复杂一点，这里引用牛人博客里的证明。然后只给出一个解的思想。首先我们只想考虑α1,α2\alpha_1,\alpha_2，而原问题：

W(α)=12(∑Ni,j=1αiyiαjyjxi∗xj)−∑Ni=1αiW(\alpha)=\dfrac{1}{2}(\sum_{i,j=1}^{N}\alpha_iy_i\alpha_jy_jx_i*x_j)-\sum_{i=1}^{N}\alpha_i

里面有所有的α\alpha，这里把这个式子乘开，把含有α1,α2\alpha_1,\alpha_2都单独拿出来，其他的作为一堆，就变成下面这样：W(α)=12K11α21+12K22α22+y1y2K12α1α2+y1α1v1+y2α2v2−α1−α2+W(α3...αn),W(\alpha) = \dfrac{1}{2}K_{11}\alpha_1^2+\dfrac{1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1\alpha_2+y_1\alpha_1v_1+y2\alpha_2v_2\\-\alpha_1-\alpha_2+W(\alpha_3...\alpha_n),

v是一个与α1,α2,y1,y2\alpha_1,\alpha_2,y_1,y_2有关的长式子，K是<x1∗x2.>内积。可以看到后面一堆与α1,α2\alpha_1,\alpha_2就没有关系。然后因为αnew1y1+αnew2y2=αold1y1+αold2y2\alpha_1^{new}y_1+\alpha_2^{new}y_2=\alpha_1^{old}y_1+\alpha_2^{old}y_2这个关系，又可以把α1\alpha_1给消掉是不是，这样新的W前面一部分只与α2\alpha_2有关，后面一部分因为不含α1,α2\alpha_1,\alpha_2所以与之没关系。而里面的αold1,αold2\alpha_1^{old},\alpha_2^{old}是上一次迭代的结果，是知道的。这样这个式子对α2\alpha_2求导再等于0，就可以解出α2\alpha_2了，应该是αnew2\alpha_2^{new}。（上面那个索引有详细的推导过程）。那么思路就是这样的，最终得到的结果为：

αnew2=αold2−y2(E1−E2)η其中Ei=ui−yiη=2K(x1,x2)−K(x1,x1)−K(x2,x2)\alpha_2^{new}=\alpha_2^{old}-\dfrac{y_2(E_1-E_2)}{\eta}\\其中\quad\quad E_i=u_i-y_i\\\eta=2K(x_1,x_2)-K(x_1,x_1)-K(x_2,x_2)

好简单的式子，然后看看αnew2\alpha_2^{new}的大小是否符合上面求出来的范围，超出了将其限制在范围内。

有了αnew2\alpha_2^{new}，再根据

αnew1y1+αnew2y2=αold1y1+αold2y2\alpha_1^{new}y_1+\alpha_2^{new}y_2=\alpha_1^{old}y_1+\alpha_2^{old}y_2就可以计算出αnew1\alpha_1^{new}了。如下：

αnew1=αold1+y1y2(αold2−αnew2)\alpha_1^{new}=\alpha_1^{old}+y_1y_2(\alpha_2^{old}-\alpha_2^{new})

因为αnew2\alpha_2^{new}已经限制范围，而这个范围就是已经认为αnew1\alpha_1^{new}的范围属于0-C下而来的，所以αnew1\alpha_1^{new}的范围一定是在0-C之间的，对于αnew1\alpha_1^{new}就不需要限制范围了。

好了有了αnew1,αnew2\alpha_1^{new},\alpha_2^{new},问题不就解决了吗？那么这一次更新以后的权值W怎么算呢？用这个式子w=∑Ni=1αiyixiw=\sum_{i=1}^{N}\alpha_iy_ix_i 就可以算出来这一次的W了，那么b呢？我们发现b还没有算呀，它怎么更新？回到问题，α1\alpha_1由αold1\alpha_1^{old}变化到了αnew1\alpha_1^{new}，那么再回到我们最最初的KKT条件，如果要改变的α1\alpha_1不为0，那么怎么样？回去看看是不是它对应的另一个乘子yi(w∗x+b)−1=0y_i(w*x+b)-1=0。w刚刚说了可以用那个求和表示，那么对于bold、bnewb_{old}、b_{new}，分别列出这么一个yi(w∗x+b)−1=0y_i(w*x+b)-1=0等式。然后把bnewb_{new}对应等式中相同的部分用boldb_{old}对应的等式里面的东西替换掉，比如说里面有一个求和，拆开后是从α3\alpha_3以后的求和(因为α1,α2\alpha_1,\alpha_2要用),这个求和在前后是一样的替换掉。那么一顿替换化简以后对应α1\alpha_1的就会有一个b1newb1_{new}，同理对于α2\alpha_2的就会有一个b2newb2_{new}，他们的最终结果如下：

bnew1=bold−E1−y1(αnew1−αold1)K(x1,x1)−y2(αnew2−αold2)K(x1,x2)bnew2=bold−E2−y1(αnew1−αold1)K(x1,x2)−y2(αnew2−αold2)K(x2,x2)b_1^{new}=b^{old}-E_1-y_1(\alpha_1^{new}-\alpha_1^{old})K(x_1,x_1)-y_2(\alpha_2^{new}-\\\alpha_2^{old})K(x_1,x_2)\\b_2^{new}=b^{old}-E_2-y_1(\alpha_1^{new}-\alpha_1^{old})K(x_1,x_2)-y_2(\alpha_2^{new}-\\\alpha_2^{old})K(x_2,x_2)

那么每次更新会出来两个b，选择哪个呢？谁准就选谁，那么怎么判断准不准呢？就是变换后的哪个α\alpha在0-C之间，也就是在分界线的边界上，也就是属于支持向量了，那么谁就准。于是乎就给一个判断吧：

b=⎧⎩⎨⎪⎪b1,b2,(b1+b2)/2if0≤αnew1≤Cif0≤αnew2≤Cothers b= \begin{cases}b_1, & {if \quad 0\le \alpha_1^{new}\le C} \\ b_2, & {if \quad 0\le \alpha_2^{new}\le C}\\ (b_1+b_2)/2 &others \end{cases}

至此我们可以说，简单的，线性的，带有松弛条件（可以容错的）的整个SMO算法就完了，剩下的就是循环，选择两个α\alpha，看是否需要更新（如果不满足KKT条件），不需要再选，需要就更新。一直到程序循环很多次了都没有选择到两个不满足KKT条件的点，也就是所有的点都满足KKT了，那么就大功告成了。

当然了，这里面还有些问题就是如何去优化这些步骤，最明显的就是怎么去选择这两个α\alpha，程序越到后期，你会发现只有那么几个点不满足KKT条件，这个时候如果你再去随机选择两个点的α\alpha，那么它是满足的，就不更新，循环，这样一直盲目的找呀找，程序的效率明显就下来了。当然这在后面是有解决办法的。

先不管那么多，就先让他盲目盲目的找吧，设置一个代数，盲目到一定代数停止就ok了，下面就来一个盲目找α\alpha的matlab程序，看看我们的SMO算法如何。

上程序之前先规整一下步骤吧：

（1）根据αold1,αold2\alpha_1^{old},\alpha_2^{old}和上面的公式找到αnew2\alpha_2^{new}的上下限。

（2）计算一个二阶导数，也就是上面的

η=2K(x1,x2)−K(x1,x1)−K(x2,x2)\eta=2K(x_1,x_2)-K(x_1,x_1)-K(x_2,x_2)，K表示里面两个向量的内积。

（3）根据公式更新αnew2\alpha_2^{new},矫正它的范围，再更新αnew1\alpha_1^{new}

（4）按照公式计算两个b，再更新b。

我的样本是这样的：




程序如下：

%%
% * svm 简单算法设计
%
%% 加载数据
% * 最终data格式：m*n，m样本数，n维度
% * label:m*1  标签必须为-1与1这两类
clc
clear
close all
data = load('data_test2.mat');
data = data.data;
train_data = data(1:end-1,:)';
label = data(end,:)';
[num_data,d] = size(train_data);
data = train_data;
%% 定义向量机参数
alphas = zeros(num_data,1);
% 系数
b = 0;
% 松弛变量影响因子
C = 0.6;
iter = 0;
max_iter = 40;
%%
while iter < max_iter
    alpha_change = 0;
    for i = 1:num_data
        %输出目标值
        pre_Li = (alphas.*label)'*(data*data(i,:)') + b;
        %样本i误差
        Ei = pre_Li - label(i);
        % 满足KKT条件
        if (label(i)*Ei<-0.001 && alphas(i)<C)||(label(i)*Ei>0.001 && alphas(i)>0)
           % 选择一个和 i 不相同的待改变的alpha(2)--alpha(j)
            j = randi(num_data,1);  
            if j == i
                temp = 1;
                while temp
                    j = randi(num_data,1);
                    if j ~= i
                        temp = 0;
                    end
                end
            end
            % 样本j的输出值
            pre_Lj = (alphas.*label)'*(data*data(j,:)') + b;
            %样本j误差
            Ej = pre_Lj - label(j);
            %更新上下限
            if label(i) ~= label(j) %类标签相同
                L = max(0,alphas(j) - alphas(i));
                H = min(C,C + alphas(j) - alphas(i));
            else
                L = max(0,alphas(j) + alphas(i) -C);
                H = min(C,alphas(j) + alphas(i));
            end
            if L==H  %上下限一样结束本次循环
                continue;end
            %计算eta
            eta = 2*data(i,:)*data(j,:)'- data(i,:)*data(i,:)' - ...
                data(j,:)*data(j,:)';
            %保存旧值
            alphasI_old = alphas(i);
            alphasJ_old = alphas(j);
            %更新alpha(2)，也就是alpha(j)
            alphas(j) = alphas(j) - label(j)*(Ei-Ej)/eta;
            %限制范围
            if alphas(j) > H
                alphas(j) = H;
            elseif alphas(j) < L
                    alphas(j) = L;
            end
            %如果alpha(j)没怎么改变，结束本次循环
            if abs(alphas(j) - alphasJ_old)<1e-4
                continue;end
            %更新alpha(1)，也就是alpha(i)
            alphas(i) = alphas(i) + label(i)*label(j)*(alphasJ_old-alphas(j));
            %更新系数b
            b1 = b - Ei - label(i)*(alphas(i)-alphasI_old)*data(i,:)*data(i,:)'-...
                label(j)*(alphas(j)-alphasJ_old)*data(i,:)*data(j,:)';
            b2 = b - Ej - label(i)*(alphas(i)-alphasI_old)*data(i,:)*data(j,:)'-...
                label(j)*(alphas(j)-alphasJ_old)*data(j,:)*data(j,:)';
            %b的几种选择机制
            if alphas(i)>0 && alphas(i)<C
                b = b1;
            elseif alphas(j)>0 && alphas(j)<C
                b = b2;
            else
                b = (b1+b2)/2;
            end
            %确定更新了，记录一次
            alpha_change = alpha_change + 1;
        end
    end
    % 没有实行alpha交换，迭代加1
    if alpha_change == 0
        iter = iter + 1;
    %实行了交换，迭代清0
    else
        iter = 0;
    end
    disp(['iter ================== ',num2str(iter)]);
end
%% 计算权值W
W = (alphas.*label)'*data;
%记录支持向量位置
index_sup = find(alphas ~= 0);
%计算预测结果
predict = (alphas.*label)'*(data*data') + b;
predict = sign(predict);
%% 显示结果
figure;
index1 = find(predict==-1);
data1 = (data(index1,:))';
plot(data1(1,:),data1(2,:),'+r');
hold on
index2 = find(predict==1);
data2 = (data(index2,:))';
plot(data2(1,:),data2(2,:),'*');
hold on
dataw = (data(index_sup,:))';
plot(dataw(1,:),dataw(2,:),'og','LineWidth',2);
% 画出分界面，以及b上下正负1的分界面
hold on
k = -W(1)/W(2);
x = 1:0.1:5;
y = k*x + b;
plot(x,y,x,y-1,'r--',x,y+1,'r--');
title(['松弛变量范围C = ',num2str(C)]);

程序中设置了松弛变量前的系数C是事先规定的，表明松弛变量项的影响程度大小。下面是几个不同C下的结果：










这是80个样本点，matlab下还是挺快2秒左右就好了。上图中，把真实分界面，上下范围为1的界面，以及那些α\alpha不为0的点（绿色标出）都画了出来，可以看到，C越大，距离越起作用，那么落在分界线之间的点就越少。同时可以看到，三种情况下，真实的分界面（蓝色）都可以将两种样本完全分开（我的样本并没有重叠，也就是完全是可分的）。

好了，这就是随机选取α\alpha的实验，第一个α\alpha是按顺序遍历所有的α\alpha，第二个α\alpha是在剩下的α\alpha中在随机选一个。当第二个α\alpha选了iter次还没有发现不满足KKT条件的，就退出循环。同时程序中的KKT条件略有不同，不过是一样的。下面介绍如何进行启发式的选取α\alpha呢？我们分析分析，比如上一次的一些点的α\alpha在0-C之间，也就是这些点不满足条件需要调整，那么一次循环后，他调整了一点，在下一次中这些点是不是还是更有可能不满足条件，因为每一次的调整比较不可能完全。而那些在上一次本身满足条件的点，那么在下一次后其实还是更有可能满足条件的。所以在启发式的寻找α\alpha过程中，我们并不是遍历所有的点的α\alpha，而是遍历那些在0-C之间的α\alpha，而0-C反应到点上就是那些属于边界之间的点是不是。当某次遍历在0-C之间找不到α\alpha了，那么我们再去整体遍历一次，这样就又会出现属于边界之间α\alpha了，然后再去遍历这些α\alpha，如此循环。那么在遍历属于边界之间α\alpha的时候，因为是需要选两个α\alpha的，第一个可以随便选，那第二个呢？这里在用一个启发式的思想，第1个α\alpha选择后，其对应的点与实际标签是不是有一个误差，属于边界之间α\alpha的所以点每个点都会有一个自己的误差，这个时候选择剩下的点与第一个α\alpha点产生误差之差最大的那个点。

程序如下：

%%
% * svm 简单算法设计 --启发式选择
%
%% 加载数据
% * 最终data格式：m*n，m样本数，n维度
% * label:m*1  标签必须为-1与1这两类
clc
clear
close all
data = load('data_test2.mat');
data = data.data;
train_data = data(1:end-1,:)';
label = data(end,:)';
[num_data,d] = size(train_data);
data = train_data;
%% 定义向量机参数
alphas = zeros(num_data,1);
b = 0;
error = zeros(num_data,2);
tol = 0.001;
C = 0.6;
iter = 0;
max_iter = 40;
%%
alpha_change = 0;
entireSet = 1;%作为一个标记看是选择全遍历还是部分遍历
while (iter < max_iter) && ((alpha_change > 0) || entireSet)
    alpha_change = 0;
    %% -----------全遍历样本-------------------------
    if entireSet 
        for i = 1:num_data
            Ei = calEk(data,alphas,label,b,i);%计算误差
            if (label(i)*Ei<-0.001 && alphas(i)<C)||...
                    (label(i)*Ei>0.001 && alphas(i)>0)
                %选择下一个alphas
                [j,Ej] = select(i,data,num_data,alphas,label,b,C,Ei,entireSet);
                alpha_I_old = alphas(i);
                alpha_J_old = alphas(j);
                if label(i) ~= label(j)
                    L = max(0,alphas(j) - alphas(i));
                    H = min(C,C + alphas(j) - alphas(i));
                else
                    L = max(0,alphas(j) + alphas(i) -C);
                    H = min(C,alphas(j) + alphas(i));
                end
                if L==H
                    continue;end
                eta = 2*data(i,:)*data(j,:)'- data(i,:)*...
                    data(i,:)' - data(j,:)*data(j,:)';
                if eta >= 0 
                    continue;end
                alphas(j) = alphas(j) - label(j)*(Ei-Ej)/eta;
                %限制范围
                if alphas(j) > H
                    alphas(j) = H;
                elseif alphas(j) < L
                    alphas(j) = L;
                end
                if abs(alphas(j) - alpha_J_old) < 1e-4
                    continue;end
                alphas(i) = alphas(i) + label(i)*...
                    label(j)*(alpha_J_old-alphas(j));
                b1 = b - Ei - label(i)*(alphas(i)-alpha_I_old)*...
                    data(i,:)*data(i,:)'- label(j)*...
                    (alphas(j)-alpha_J_old)*data(i,:)*data(j,:)';
                b2 = b - Ej - label(i)*(alphas(i)-alpha_I_old)*...
                    data(i,:)*data(j,:)'- label(j)*...
                    (alphas(j)-alpha_J_old)*data(j,:)*data(j,:)';
                if (alphas(i) > 0) && (alphas(i) < C)
                    b = b1;
                elseif (alphas(j) > 0) && (alphas(j) < C)
                    b = b2;
                else
                    b = (b1+b2)/2;
                end
                alpha_change = alpha_change + 1;
            end
        end
         iter = iter + 1;
   %% --------------部分遍历(alphas=0~C)的样本--------------------------
    else
        index = find(alphas>0 & alphas < C);
        for ii = 1:length(index)
            i = index(ii);
            Ei = calEk(data,alphas,label,b,i);%计算误差
            if (label(i)*Ei<-0.001 && alphas(i)<C)||...
                    (label(i)*Ei>0.001 && alphas(i)>0)
                %选择下一个样本
                [j,Ej] = select(i,data,num_data,alphas,label,b,C,Ei,entireSet);
                alpha_I_old = alphas(i);
                alpha_J_old = alphas(j);
                if label(i) ~= label(j)
                    L = max(0,alphas(j) - alphas(i));
                    H = min(C,C + alphas(j) - alphas(i));
                else
                    L = max(0,alphas(j) + alphas(i) -C);
                    H = min(C,alphas(j) + alphas(i));
                end
                if L==H
                    continue;end
                eta = 2*data(i,:)*data(j,:)'- data(i,:)*...
                    data(i,:)' - data(j,:)*data(j,:)';
                if eta >= 0
                    continue;end
                alphas(j) = alphas(j) - label(j)*(Ei-Ej)/eta;  
                %限制范围
                if alphas(j) > H
                    alphas(j) = H;
                elseif alphas(j) < L
                    alphas(j) = L;
                end
                if abs(alphas(j) - alpha_J_old) < 1e-4
                    continue;end
                alphas(i) = alphas(i) + label(i)*...
                    label(j)*(alpha_J_old-alphas(j));
                b1 = b - Ei - label(i)*(alphas(i)-alpha_I_old)*...
                    data(i,:)*data(i,:)'- label(j)*...
                    (alphas(j)-alpha_J_old)*data(i,:)*data(j,:)';
                b2 = b - Ej - label(i)*(alphas(i)-alpha_I_old)*...
                    data(i,:)*data(j,:)'- label(j)*...
                    (alphas(j)-alpha_J_old)*data(j,:)*data(j,:)';
                if (alphas(i) > 0) && (alphas(i) < C)
                    b = b1;
                elseif (alphas(j) > 0) && (alphas(j) < C)
                    b = b2;
                else
                    b = (b1+b2)/2;
                end
                alpha_change = alpha_change + 1;
            end
        end
        iter = iter + 1;
    end
    %% --------------------------------
    if entireSet %第一次全遍历了，下一次就变成部分遍历
        entireSet = 0;
    elseif alpha_change == 0 
        %如果部分遍历所有都没有找到需要交换的alpha，再改为全遍历
        entireSet = 1;
    end
    disp(['iter ================== ',num2str(iter)]);    
end
%% 计算权值W
W = (alphas.*label)'*data;
%记录支持向量位置
index_sup = find(alphas ~= 0);
%计算预测结果
predict = (alphas.*label)'*(data*data') + b;
predict = sign(predict);
%% 显示结果
figure;
index1 = find(predict==-1);
data1 = (data(index1,:))';
plot(data1(1,:),data1(2,:),'+r');
hold on
index2 = find(predict==1);
data2 = (data(index2,:))';
plot(data2(1,:),data2(2,:),'*');
hold on
dataw = (data(index_sup,:))';
plot(dataw(1,:),dataw(2,:),'og','LineWidth',2);
% 画出分界面，以及b上下正负1的分界面
hold on
k = -W(1)/W(2);
x = 1:0.1:5;
y = k*x + b;
plot(x,y,x,y-1,'r--',x,y+1,'r--');
title(['松弛变量范围C = ',num2str(C)]);

其中的子函数，一个是计算误差函数，一个是选择函数如下：

function Ek = calEk(data,alphas,label,b,k)
pre_Li = (alphas.*label)'*(data*data(k,:)') + b;
Ek = pre_Li - label(k);

function [J,Ej] = select(i,data,num_data,alphas,label,b,C,Ei,choose)
maxDeltaE = 0;maxJ = -1;
if choose == 1 %全遍历---随机选择alphas
    j = randi(num_data ,1);
    if j == i
        temp = 1;
        while temp
            j = randi(num_data,1);
            if j ~= i
                temp = 0;
            end
        end
    end
    J = j;
    Ej = calEk(data,alphas,label,b,J);
else %部分遍历--启发式的选择alphas
    index = find(alphas>0 & alphas < C);
    for k = 1:length(index)
        if i == index(k)
            continue;
        end
        temp_e = calEk(data,alphas,label,b,k);
        deltaE = abs(Ei - temp_e); %选择与Ei误差最大的alphas
        if deltaE > maxDeltaE
            maxJ = k;
            maxDeltaE = deltaE;
            Ej = temp_e;
        end
    end
    J = maxJ;
end

至此算是完了，试验了一下，两者的效果其实差不多（反而随机选取的效果更好一点，感觉是因为随机保证了更多的可能，毕竟随机选择包括了你的特殊选择，但是特殊选择到后期是特殊不起来的，反而随机会把那些差一点的选择出来），但是速度上当样本小的时候，基本上差不多，但是当样本大的时候，启发式的特殊选择明显占优势了。我试验了400个样本点的情况，随机选择10多秒把，而启发式2,3秒就好了。可见效果差不多的情况下，启发式选择是首要选择。

至此两种方式下的方法都实验完了。那么我们看到，前面(三节)所讲的一切以及实验，分类的样本都是线性样本，那么如果来了非线性样本该如何呢？而SVM的强大之处更在于对非线性样本的准确划分。那么前面的理论对于非线性样本是否适用？我们又该如何处理非线性样本呢？请看下节SVM非线性样本的分类。