解密SVM系列（四）：SVM非线性分类原理实验

前面几节我们讨论了SVM原理、求解线性分类下SVM的SMO方法。本节将分析SVM处理非线性分类的相关问题。

一般的非线性分类如下左所示（后面我们将实战下面这种情况）：

可以看到在原始空间中你想用一个直线分类面划分开来是不可能了，除非圆。而当你把数据点映射一下成右图所示的情况后，现在数据点明显看上去是线性可分的，那么在这个空间上的数据点我们再用前面的SVM算法去处理，就可以得到每个数据点的分类情况了，而这个分类情况也是我们在低维空间的情况。也就是说，单纯的SVM是不能处理非线性问题的，说白了只能处理线性问题，但是来了非线性样本怎么办呢？我们是在样本上做的文章，我把非线性样本变成线性样本，再去把变化后的线性样本拿去分类，经过这么一圈，就达到了把非线性样本分开的目的，所以只看开头和结尾的话发现，SVM竟然可以分非线性问题，其实呢还是分的线性问题。

现在的问题是如何找到这个映射关系对吧，就比如上面那个情况，我们可以人为计算出这种映射，比如一个样本点是用坐标表示的(x1,x2),它有个类标签，假设为1，那么把这个点映射到三维中变成(x21,2√x1x2,x22)。对每个点我都这么去映射，假设一个原始点样本集是这样的：

然后按照上面那个公式去把每个点映射成3维坐标点后，画出来是这样的：

可以看到是线性可分的吧，如果还看不清把视角换个角度（右视图）：

现在能看清楚了吧。

那这是二维的点到三维，映射的关系就是上面的那个关系，那如果是三维到四维，四维到N维呢？这个关系你还想去找吗？理论上是找的到的，但是实际上人工去找你怎么去找？你怎么知道数据的映射关系是这样的是那样的？不可能知道。然而我们真的需要找到这种关系吗？答案是不需要的，返回去看看前三节的关于SVM的理论部分可以看到，无论是计算α呀，还是b呀等等，只要涉及到原始数据点的，都是以内积的形式出来的，也就是说是一个点的向量与另一个点的向量相乘的，向量内积出来是一个值。

就拿α的更新来说，如下：

αnew2=αold2−y2(E1−E2)η其中Ei=ui−yiη=2K(x1,x2)−K(x1,x1)−K(x2,x2)

可以看到α的更新用到原始数据点的就在η的计算上。而η的计算可以看到K表示的是内积，就是两个点的内积或者自己内积。假设现在两个点的坐标非别是X1=(x1,y1),X2=(x2,y2)，那么
K(X1,X2)=XT1∗X2=x1x2+y1y2=C1
其他依次类推，K出来的是一个数比如C1。那么假设这个样本点都映射到三维以后了，每个样本点是一个1*3的向量，那么计算到这一块来了，就变成了下面这样：K(X1,X2)=XT1∗X2=(x1,y1,z1)T∗(x2,y2,z2)=x1y1+x2y2+x3y3=C2
最后也是得到一个值比如C2。既然SVM里面所有涉及到原始数据的地方都是以向量的形式出现的，那么我们还需要管它的映射关系吗？因为它也不需要你去计算说具体到比如说三维以后，三维里面的三个坐标值究竟是多少，他需要的是内积以后的一个结果值。那么好办了，我就假设有一个黑匣子，输入原始数据维度下的两个坐标向量，然后经过黑匣子这么一圈，出来一个值，这个值我们就认为是高维度下的值。而黑匣子的潜在意义就相当于一个高维映射器一样。更重要的是我们并不需要知道黑匣子究竟是怎么映射的，只需要知道它的低纬度下的形式就可以了。常用的黑匣子就是径向基函数，而这个黑匣子在数学上就叫做核函数，例如径向基函数的外在形式如下所示：

K(X1,X2)=exp(−||X1−X2||2σ)

σ是需要预先设定的参数。至于这个黑匣子把初始数据映射到多少维了，谁知道呢，既然是黑匣子，那就是看不到的，上帝给了人类这么一个黑匣子就已经很够意思了。可以看到的是原始数据结果黑匣子算了以后，出来就是一个值了，而这个值就认为是高维度下的数据通过内积计算而来的值。当然上帝还留了一个窗户，就是σ，相传σ选取的越小，数据映射的维度越大，小到一定程度，维度空间大到无穷维。反之越大，映射的维度空间就越小，但是会不会小到低于原始空间维度呢？谁知道了，然而通过实验我发现，大到一定程度，样本点分的乱七八糟，并且σ正好在一定范围的时候效果非常好，这个范围既不是极小的范围，也不是极大的范围，那这暗示了什么呢？也就是说非线性原始样本是有一个属于他自己的最佳高维空间的，大了小了似乎都不好。

好了既然黑匣子是藏着的，那也就只能说这么多了。有趣的是上帝给的这个黑匣子不止一个，有好几个，只是上面的那个普遍效果更好而已。基于此，那么对于上节的SMO算法，如果拿来求解非线性数据的话，我们只需要将其中对应的内积部分改成核函数的形式即可。一个数据核函数程序如下：

function result = Kernel(data1,data2,sigma)% data里面每一行数据是一个样本(的行向量)[m1,~] = size(data1);[m2,~] = size(data2);result = zeros(m1,m2);for i = 1:m1    for j = 1:m2        result(i,j) = exp(-norm(data1(i,:)-data2(j,:))/(2*sigma^2));    endend

有了此核函数，我们用上节的随机遍历α的方式（这个函数代码少一点）来实验一下非线性样本，非线性样本如下：

然后把主程序对应的部分用上述核函数代替：

%%% * svm 简单算法设计%%% 加载数据% * 最终data格式：m*n，m样本数，n维度% * label:m*1  标签必须为-1与1这两类clcclear% close alldata = load('data_test1.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 = 80;% 核函数的参数sigma = 4;%%while iter < max_iter    alpha_change = 0;    for i = 1:num_data        %输出目标值        pre_Li = (alphas.*label)'*Kernel(data,data(i,:),sigma) + 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)'*Kernel(data,data(j,:),sigma) + 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*Kernel(data(i,:),data(j,:),sigma)- ...                Kernel(data(i,:),data(i,:),sigma)...                - Kernel(data(j,:),data(j,:),sigma);            %保存旧值            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)*...                Kernel(data(i,:),data(i,:),sigma) - label(j)*...                (alphas(j)-alphasJ_old)*Kernel(data(i,:),data(j,:),sigma);            b2 = b - Ej - label(i)*(alphas(i)-alphasI_old)*...                Kernel(data(i,:),data(j,:),sigma)- label(j)*...                (alphas(j)-alphasJ_old)*Kernel(data(j,:),data(j,:),sigma);            %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)'*Kernel(data,data,sigma) + b;predict = sign(predict);%% 显示结果figure;index1 = find(predict==-1);data1 = (data(index1,:))';plot(data1(1,:),data1(2,:),'+r');hold onindex2 = find(predict==1);data2 = (data(index2,:))';plot(data2(1,:),data2(2,:),'*');hold ondataw = (data(index_sup,:))';plot(dataw(1,:),dataw(2,:),'og','LineWidth',2);title(['核函数参数sigma = ',num2str(sigma)]);

下面是几个不同参数下的结果：




可以看到σ到4以后就分不出来了。绿色的为支持向量，可以看到在σ在0.6到1之间是最少的，结果应该也是最好的。至此SMO实验非线性样本完毕。

当今学者已经有非常多的人研究SVM算法，同时开发了许多开源的程序，这些程序都是经过不断优化的，性能比起我们这里自己编的来说要好得多，所以在实际应用中通常都是用他们无私贡献的软件包。一个典型的软件包就是台湾一个教授团队的LIBSVM软件包，那么你是否想一窥其用法，看看它的性能如何呢？请看下节matlab下LIBSVM的简单使用。