机器学习之旅---SVM分类器

本次内容主要讲解什么是支持向量，SVM分类是如何推导的，最小序列SMO算法部分推导。

最后给出线性和非线性2分类问题的smo算法matlab实现代码。

一、什么是支持向量机(Support Vector Machine)

本节内容部分翻译Opencv教程：

http://docs.opencv.org/doc/tutorials/ml/introduction_to_svm/introduction_to_svm.html#introductiontosvms

 

SVM是一个由分类超平面定义的判别分类器。也就是说给定一组带标签的训练样本，算法将会输出一个最优超平面对新样本(测试样本)进行分类。

什么样的超平面才是最优的？我们考虑如下场景：对于一个由二维坐标点构成的线性可分集合，如何找到一条直线将两类坐标分隔开？

      上图中，你可以看到存在多条直线将两类坐标分开。它们中有没有最好的？我们可以直观地定义如下规则：如果一条分割的直线离坐标点太近，那它就不是最佳的。因为它会对噪声敏感，不能正确的推广。因此，我们的目标是找到一条分割线，它要离所有的样本点都尽可能的远。

       SVM算法就是找一个超平面，并且它到离他最近的训练样本的距离要最大。即最优分割超平面最大化训练样本边界。

原谅哥的翻译水平，看着自己也别扭，重新整理下：

对线性可分集，总能找到使样本正确划分的分界面，而且有无穷多个，哪个是最优？ 必须寻找一种最优的分界准则，是两类模式分开的间隔最大。

在介绍如何推到之前，了解一个定义：

分隔超平面：上述将数据集分割开来的直线叫做分隔超平面。

超平面：如果数据集是N维的，那么就需要N-1维的某对象来对数据进行分割。该对象叫做超平面，也就是分类的决策边界。

间隔：

一个点到分割面的距离，称为点相对于分割面的距离。

数据集中所有的点到分割面的最小间隔的2倍，称为分类器或数据集的间隔。

最大间隔：SVM分类器是要找最大的数据集间隔。

支持向量：离分割超平面最近的那些点。

SVM算法其实就是最大化支持向量到分割面的距离。另外，目前上面讲的都是针对线性可分的数据集。非线性的数据集，需要核函数转换空间，才具有非线性数据处理能力。

最优分类超平面只由少数支持向量决定，问题具有稀疏性。

 

二、推导

下面讲解下SVM的数学原理，其实就是复述我之前写过的东西。

http://blog.csdn.net/jinshengtao/article/details/17636953

（这里yi(wTx+b)称为点到分割面的函数间隔。yi(wTx+b)/|w|称为点到分割面的几何间隔）

            现在的目标就是找出分类器定义中的w和b。为此，我们必须找到具有最小间隔的数据点，而这些数据点也就是前面提到的支持向量。一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化。这就可以写作：

          直接求解上面的公式很难。并且到目前为止，我们知道所有数都不那么干净，不能100%线性可分。我们可以通过引入松弛变量，来允许有些数据点可以处于分割面的错误的一侧。如下面几幅图所示：

 εi是松弛变量，常数C代表惩罚系数，即如果某个x是属于某一类，但是它偏离了该类，跑到边界上后者其他类的地方去了，C越大表明越不想放弃这个点，边界就会缩小，错分点更少，但是过拟合情况会比较严重。

为此，我们引入拉格朗日乘子，用条件极值求解最优分界面，构造拉格朗日函数，也是之前公式(3)的对偶形式：

也就是说，只要求的拉格朗日乘子ai，我们就能得到分类边界。

三、序列最小化方法SMO

    C-SVM是一个凸二次规划问题，SMO算法能快速解决SVM的二次规划问题，把它分解成一些子问题来解决。在每一步中，SMO算法仅选择两个拉格朗日乘子进行优化，然后再更新SVM以反应新的优化值。

   对于C-SVM目标函数和优化过程中必须遵循的约束条件如下：

     我们需要一种算法来求解最优化问题。（在以前的实现，我用了matlab提供的求条件极值的函数，http://blog.csdn.net/jinshengtao/article/details/17675653）。

现在打算自己用smo算法实现，下面就重点讲下smo算法。

       C-SVM的KTT条件为：

(KTT条件是指在满足一些有规则的条件下, 一个非线性规划(Nonlinear Programming)问题能有最优化解法的一个必要和充分条件.即最优解都需要满足KTT条件)。现在这里对于C-SVM的KTT条件为：

其中，b*由上面的公式计算得到。

SMO算法包括两个步骤：一是用分析的方法解一个简单的优化问题；二是选择待优化的拉格朗日乘子的策略。

(1)   简单的优化问题----两个拉格朗日乘子优化问题的解

    为了求解只有两个乘子的优化问题【因为(xi,yi)都已知】，SMO方法首先计算它的约束，然后再解带有约束的最小化问题。对于两个乘子，其二维约束很容易表达：边界约束使得乘子在方框内，而线性等式约束使得乘子在对角线上，这里yi∈{1,-1}，

    现在，在一条线段上求目标函数的极值，相当于一个一维的极值问题。我们可以把a1用a2表示，对a2求无条件极值，如果目标函数是严格凹的，最小值就一定在这一极值点(极值点在区间内)或在区间端点(极值点在区间外)。a2确定下来后，a1也确定了。因此，我们要先找到a2的优化区间，再在这个区间中对a2求最小值。

根据上面两幅图可以知道，a2优化区间：

当y1和y2异号时：【L、H的计算是由a2-a1正负未知所导致的】

当y1和y2同号时：

(b)当L的二阶导数为0，即训练样本中出现相同的特征时，我们计算目标函数是单调的。最小值在线段两个端点上的取值，我们将a2=L和a2=H分别代入目标函数，这样拉格朗日乘子就修正到目标函数较小的端点上，比较ΨL和ΨH就可以求得Ψ的最小值。最终得到a2的值，接着一串值都解出来了。

至于偏置值b，每一步都要更新，因为前面的KKT条件指出了ai和yiui关系，而ui和b有关，在每一步计算出ai后，根据KKT条件来调整b，分如下几种情况

(2)   选择策略----选取两个拉格朗日乘子的启发规则

事实上即使我们不采用任何找点法，只是按顺序抽取ai，aj的所有组合进行优化，目标函数也会不断下降。直到任一对ai，aj都不能继续优化，目标函数就会收敛到极小值。我们采取某种找点方法只是为了使算法收敛得更快。

这种试探法先选择最有可能需要优化的a2，再针对这样的a2选择最有可能取得较大修正步长的a1。这样，我们在程序中使用两个层次的循环：

外层循环首先遍历所有样本查找违反KTT规则的乘子进行优化，当完成一次遍历后外层循环再遍历那些非边界乘子的样本(0<a<C)，挑选那些违反KTT条件的样本进行优化。外层循环重复执行，直到所有的非边界样本在一定范围内均满足KKT条件。

    在选定第一个拉格朗日乘子ai后，内层循环会通过最大化步长的方式来挑选第二个拉格朗日乘子，即最大化|Ei-Ej|,当Ei为正时最小化Ej，当为负Ei时最大化Ej.

下面给出matlab代码实现 ：

本次内容主要讲解什么是支持向量，SVM分类是如何推导的，最小序列SMO算法部分推导。

最后给出线性和非线性2分类问题的smo算法matlab实现代码。

一、什么是支持向量机(Support Vector Machine)

本节内容部分翻译Opencv教程：

http://docs.opencv.org/doc/tutorials/ml/introduction_to_svm/introduction_to_svm.html#introductiontosvms

 

SVM是一个由分类超平面定义的判别分类器。也就是说给定一组带标签的训练样本，算法将会输出一个最优超平面对新样本(测试样本)进行分类。

什么样的超平面才是最优的？我们考虑如下场景：对于一个由二维坐标点构成的线性可分集合，如何找到一条直线将两类坐标分隔开？

      上图中，你可以看到存在多条直线将两类坐标分开。它们中有没有最好的？我们可以直观地定义如下规则：如果一条分割的直线离坐标点太近，那它就不是最佳的。因为它会对噪声敏感，不能正确的推广。因此，我们的目标是找到一条分割线，它要离所有的样本点都尽可能的远。

       SVM算法就是找一个超平面，并且它到离他最近的训练样本的距离要最大。即最优分割超平面最大化训练样本边界。

原谅哥的翻译水平，看着自己也别扭，重新整理下：

对线性可分集，总能找到使样本正确划分的分界面，而且有无穷多个，哪个是最优？ 必须寻找一种最优的分界准则，是两类模式分开的间隔最大。

在介绍如何推到之前，了解一个定义：

分隔超平面：上述将数据集分割开来的直线叫做分隔超平面。

超平面：如果数据集是N维的，那么就需要N-1维的某对象来对数据进行分割。该对象叫做超平面，也就是分类的决策边界。

间隔：

一个点到分割面的距离，称为点相对于分割面的距离。

数据集中所有的点到分割面的最小间隔的2倍，称为分类器或数据集的间隔。

最大间隔：SVM分类器是要找最大的数据集间隔。

支持向量：离分割超平面最近的那些点。

SVM算法其实就是最大化支持向量到分割面的距离。另外，目前上面讲的都是针对线性可分的数据集。非线性的数据集，需要核函数转换空间，才具有非线性数据处理能力。

最优分类超平面只由少数支持向量决定，问题具有稀疏性。

 

二、推导

下面讲解下SVM的数学原理，其实就是复述我之前写过的东西。

http://blog.csdn.net/jinshengtao/article/details/17636953

（这里yi(wTx+b)称为点到分割面的函数间隔。yi(wTx+b)/|w|称为点到分割面的几何间隔）

            现在的目标就是找出分类器定义中的w和b。为此，我们必须找到具有最小间隔的数据点，而这些数据点也就是前面提到的支持向量。一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化。这就可以写作：

          直接求解上面的公式很难。并且到目前为止，我们知道所有数都不那么干净，不能100%线性可分。我们可以通过引入松弛变量，来允许有些数据点可以处于分割面的错误的一侧。如下面几幅图所示：

 εi是松弛变量，常数C代表惩罚系数，即如果某个x是属于某一类，但是它偏离了该类，跑到边界上后者其他类的地方去了，C越大表明越不想放弃这个点，边界就会缩小，错分点更少，但是过拟合情况会比较严重。

为此，我们引入拉格朗日乘子，用条件极值求解最优分界面，构造拉格朗日函数，也是之前公式(3)的对偶形式：

也就是说，只要求的拉格朗日乘子ai，我们就能得到分类边界。

三、序列最小化方法SMO

    C-SVM是一个凸二次规划问题，SMO算法能快速解决SVM的二次规划问题，把它分解成一些子问题来解决。在每一步中，SMO算法仅选择两个拉格朗日乘子进行优化，然后再更新SVM以反应新的优化值。

   对于C-SVM目标函数和优化过程中必须遵循的约束条件如下：

     我们需要一种算法来求解最优化问题。（在以前的实现，我用了matlab提供的求条件极值的函数，http://blog.csdn.net/jinshengtao/article/details/17675653）。

现在打算自己用smo算法实现，下面就重点讲下smo算法。

       C-SVM的KTT条件为：

(KTT条件是指在满足一些有规则的条件下, 一个非线性规划(Nonlinear Programming)问题能有最优化解法的一个必要和充分条件.即最优解都需要满足KTT条件)。现在这里对于C-SVM的KTT条件为：

其中，b*由上面的公式计算得到。

SMO算法包括两个步骤：一是用分析的方法解一个简单的优化问题；二是选择待优化的拉格朗日乘子的策略。

(1)   简单的优化问题----两个拉格朗日乘子优化问题的解

    为了求解只有两个乘子的优化问题【因为(xi,yi)都已知】，SMO方法首先计算它的约束，然后再解带有约束的最小化问题。对于两个乘子，其二维约束很容易表达：边界约束使得乘子在方框内，而线性等式约束使得乘子在对角线上，这里yi∈{1,-1}，

    现在，在一条线段上求目标函数的极值，相当于一个一维的极值问题。我们可以把a1用a2表示，对a2求无条件极值，如果目标函数是严格凹的，最小值就一定在这一极值点(极值点在区间内)或在区间端点(极值点在区间外)。a2确定下来后，a1也确定了。因此，我们要先找到a2的优化区间，再在这个区间中对a2求最小值。

根据上面两幅图可以知道，a2优化区间：

当y1和y2异号时：【L、H的计算是由a2-a1正负未知所导致的】

当y1和y2同号时：

(b)当L的二阶导数为0，即训练样本中出现相同的特征时，我们计算目标函数是单调的。最小值在线段两个端点上的取值，我们将a2=L和a2=H分别代入目标函数，这样拉格朗日乘子就修正到目标函数较小的端点上，比较ΨL和ΨH就可以求得Ψ的最小值。最终得到a2的值，接着一串值都解出来了。

至于偏置值b，每一步都要更新，因为前面的KKT条件指出了ai和yiui关系，而ui和b有关，在每一步计算出ai后，根据KKT条件来调整b，分如下几种情况

(2)   选择策略----选取两个拉格朗日乘子的启发规则

事实上即使我们不采用任何找点法，只是按顺序抽取ai，aj的所有组合进行优化，目标函数也会不断下降。直到任一对ai，aj都不能继续优化，目标函数就会收敛到极小值。我们采取某种找点方法只是为了使算法收敛得更快。

这种试探法先选择最有可能需要优化的a2，再针对这样的a2选择最有可能取得较大修正步长的a1。这样，我们在程序中使用两个层次的循环：

外层循环首先遍历所有样本查找违反KTT规则的乘子进行优化，当完成一次遍历后外层循环再遍历那些非边界乘子的样本(0<a<C)，挑选那些违反KTT条件的样本进行优化。外层循环重复执行，直到所有的非边界样本在一定范围内均满足KKT条件。

    在选定第一个拉格朗日乘子ai后，内层循环会通过最大化步长的方式来挑选第二个拉格朗日乘子，即最大化|Ei-Ej|,当Ei为正时最小化Ej，当为负Ei时最大化Ej.

下面给出matlab代码实现 ：

1.线性可分简单smo

[plain] view plain copy

1. function [b,alphas] = smoSimple(data, class, C, toler, maxIter)  
2. b = 0;  
3. [m,n] = size(data);  
4. alphas = zeros(m,1);  
5. iter=0;  
6. while (iter < maxIter)  
7.     alphasChanges = 0;  
8.     for k=1:1:m  
9.         fxk = (alphas .* class)' * data * data(k,:)' + b;   % f = wx+b  
10.         ek = fxk - class(k);  
11.         if (((ek*class(k) < -toler) && (alphas(k) < C)) || ((ek*class(k) > toler) && (alphas(k) > 0)))  
12.             j = selectJrand(k,m);  
13.             fxj = (alphas .* class)' * data * data(j,:)' + b;   % f = wx+b  
14.             ej = fxj - class(j);  
15.               
16.             temp_k = alphas(k);  
17.             temp_j = alphas(j);  
18.             if(class(k) ~= class(j))  
19.                 L = max(0, alphas(j) - alphas(k));  
20.                 H = min(C, C + alphas(j) - alphas(k));  
21.             else  
22.                 L = max(0, alphas(k) + alphas(j) - C);  
23.                 H = min(C, alphas(k) + alphas(j));  
24.             end  
25.             if L == H  
26.                 continue;  
27.             end  
28.             eta = 2.0 * data(k,:) * data(j,:)' - data(k,:) * data(k,:)' - data(j,:) * data(j,:)';  
29.             if eta >= 0  
30.                 continue;  
31.             end  
32.             alphas(j) = alphas(j) - class(j) * (ek - ej) / eta;  
33.             alphas(j) = clipalpha(alphas(j), H, L);  
34.               
35.             if(abs(alphas(j) - temp_j) < 0.00001)  
36.                 continue;  
37.             end  
38.               
39.             alphas(k) = alphas(k) + class(k) * class(j) * (temp_j - alphas(j));  
40.             b1 = b - ek - class(k) * (alphas(k) - temp_k) * data(k,:) * data(k,:)' - class(j) * (alphas(j) - temp_j) * data(k,:) * data(j,:)';  
41.             b2 = b - ej - class(k) * (alphas(k) - temp_k) * data(k,:) * data(j,:)' - class(j) * (alphas(j) - temp_j) * data(j,:) * data(j,:)';   
42.               
43.             if (alphas(k) > 0 && alphas(k) < C)  
44.                 b = b1;  
45.             elseif(alphas(j) > 0 && alphas(j) < C)  
46.                 b = b2;  
47.             else  
48.                 b = (b1 + b2)/2;  
49.             end  
50.             alphasChanges = alphasChanges + 1;  
51.         end   
52.     end  
53.     if alphasChanges == 0  
54.         iter = iter + 1;  
55.     else  
56.         iter = 0;  
57.     end  
58. end  
59. end  
60.   
61. function index = selectJrand(k,m)  
62.     index = k;  
63.     while(index == k)  
64.         index = randi([1,m],1,1);    
65.     end  
66. end  
67.   
68. function res = clipalpha(a, H, L)  
69.     if a > H  
70.         a = H;  
71.     end  
72.       
73.     if a < L  
74.         a = L;  
75.     end  
76.     res = a;  
77. end  

[plain] view plain copy

1. clc;  
2. clear;  
3.   
4. load Data  
5.   
6. [r,c] = size(Data);  
7. Test = Data(:,1:2);  
8. Label = Data(:,3);  
9.   
10. [b, alphas] = smoSimple(Test, Label, 0.6, 0.001, 40);  
11.    
12. %%画图  
13. figure(1)  
14. axis([-2 12 -8 6])  
15. for k = 1:1:r  
16.     hold on  
17.     if Data(k,3) == 1  
18.         plot(Data(k,1),Data(k,2),'r+');  
19.     else  
20.         plot(Data(k,1),Data(k,2),'b*');  
21.     end  
22. end  
23.   
24. %画支持向量及分割面  
25. %result=[];  
26. for k=1:1:r  
27.     if alphas(k)~= 0  
28.         hold on  
29.         %result =[result;alphas(k)];  
30.         QX = plot(Data(k,1:1),Data(k,2:2),'Ok','MarkerSize',12);  
31.         set(QX,'LineWidth',2.0);  
32.     end  
33. end  
34. W=(alphas.*Label)'*Data(:,1:2);  
35. y=(-W(1).* Data(:,1:1)-b) ./W(2);  
36. plot(Data(:,1:1),y);  

结果：

上述代码运行事件较长，此处仅是100个点的小规模数据集，对于更大的数据集收敛时间更长。

2.完整的smo算法

[plain] view plain copy

1. function [b, res_alphas] = smoP(data, class, C, toler, maxIter)  
2.     [m,n] = size(data);  
3.     iter = 0;  
4.     entireSet = 1;  
5.     alphaPairsChanged = 0;  
6.     oS = init(data,class,C,toler,m);  
7.           
8.     while(((iter<maxIter)&&(alphaPairsChanged > 0)) || (entireSet == 1))  
9.         alphaPairsChanged = 0;  
10.         if entireSet == 1  
11.             for k = 1:1:m  
12.                 [ret, oS] = innerL(k, oS);  
13.                 alphaPairsChanged = alphaPairsChanged + ret;  
14.             end  
15.             iter = iter + 1;  
16.         else  
17.             nonBoundIs = [];  
18.             for k = 1:1:m  
19.                if ((oS.alphas(k) < C) && (oS.alphas(k) > 0))  
20.                    nonBoundIs = [nonBoundIs k];  
21.                end  
22.             end  
23.             [r,c] = size(nonBoundIs);  
24.             for k = 1:1:c  
25.                 index = nonBoundIs(k);  
26.                 [ret, oS] = innerL(index, oS);  
27.                 alphaPairsChanged = alphaPairsChanged + ret;  
28.             end  
29.             iter = iter + 1;  
30.         end  
31.         if entireSet == 1  
32.             entireSet = 0;  
33.         elseif alphaPairsChanged == 0  
34.             entireSet = 1;  
35.         end  
36.     end  
37.     b = oS.b;  
38.     res_alphas = oS.alphas;  
39. end  
40.   
41. function oS = init(data,class,C,toler,m)  
42.     alphas = zeros(m,1);  
43.     eCache = zeros(m,2);  
44.     b = 0;  
45.       
46.     oS.data = data;  
47.     oS.class = class;  
48.     oS.C = C;  
49.     oS.toler = toler;  
50.     oS.m = m;  
51.     oS.alphas = alphas;  
52.     oS.b = b;  
53.     oS.eCache = eCache;  
54.       
55. end  
56.   
57. function [ret,oS] = innerL(k, oS)  
58.     Ei = calcEk(oS, k);  
59.     if(((oS.class(k)*Ei < -oS.toler) && (oS.alphas(k) < oS.C)) || ((oS.class(k)*Ei > oS.toler) && (oS.alphas(k) > 0)))  
60.         [j, Ej] = selectJ(k, oS, Ei);  
61.         temp_k = oS.alphas(k);  
62.         temp_j = oS.alphas(j);  
63.           
64.         if oS.class(k) ~= oS.class(j)  
65.             L = max(0, oS.alphas(j) - oS.alphas(k));  
66.             H = min(oS.C, oS.C +oS.alphas(j) - oS.alphas(k));  
67.         else  
68.             L = max(0, oS.alphas(j) + oS.alphas(k) - oS.C);  
69.             H = min(oS.C, oS.alphas(j) + oS.alphas(k));  
70.         end  
71.         if L == H  
72.             ret = 0;  
73.             return;  
74.         end  
75.         eta = 2.0 * oS.data(k,:) * oS.data(j,:)' - oS.data(k,:) * oS.data(k,:)' - oS.data(j,:) * oS.data(j,:)';  
76.         if eta >=0   
77.             ret = 0;  
78.             return;  
79.         end  
80.         oS.alphas(j) = oS.alphas(j) - oS.class(j) * (Ei - Ej) / eta;  
81.         oS.alphas(j) = clipalpha(oS.alphas(j), H, L);  
82.           
83.         %update Ek  
84.         Et = calcEk(oS, j);  
85.         oS.eCache(j,:) = [1 Et];  
86.           
87.         if(abs(oS.alphas(j) - temp_j) < 0.00001)  
88.             ret = 0;  
89.             return;  
90.         end  
91.           
92.         oS.alphas(k) =   oS.alphas(k) + oS.class(j)*oS.class(k)*(temp_j - oS.alphas(j));  
93.         Et = calcEk(oS, k);  
94.         oS.eCache(k,:) = [1 Et];  
95.           
96.         b1 = oS.b - Ei - oS.class(k) * (oS.alphas(k) - temp_k) * oS.data(k,:) * oS.data(k,:)' - oS.class(j) * (oS.alphas(j) - temp_j) * oS.data(k,:) * oS.data(j,:)';  
97.         b2 = oS.b - Ej - oS.class(k) * (oS.alphas(k) - temp_k) * oS.data(k,:) * oS.data(j,:)' - oS.class(j) * (oS.alphas(j) - temp_j) * oS.data(j,:) * oS.data(j,:)';   
98.           
99.         if (oS.alphas(k)>0) && (oS.alphas(k)<oS.C)  
100.             oS.b = b1;  
101.         elseif(oS.alphas(j)>0) && (oS.alphas(j)<oS.C)  
102.             oS.b = b2;  
103.         else  
104.             oS.b = (b1+b2)/2.0;  
105.         end  
106.         ret = 1;  
107.         return;  
108.     else  
109.         ret = 0;  
110.         return;  
111.     end  
112. end  
113.   
114. function Ek = calcEk(oS, k)  
115.     fXk = (oS.alphas .* oS.class)' * oS.data * oS.data(k,:)' + oS.b;  
116.     Ek = fXk - oS.class(k);  
117. end  
118.   
119. function [j, Ej] = selectJ(k, oS, Ei)  
120.     maxK = -1;  
121.     maxDeltaE = 0;   
122.     Ej = 0;  
123.     oS.eCache(k,:) =[1 Ei];  
124.     validEcacheList = [];  
125.       
126.     for l = 1:1:oS.m  
127.         if oS.eCache(l,1:1) ~= 0  
128.             validEcacheList = [validEcacheList l];  
129.         end  
130.     end  
131.     [r, c] = size(validEcacheList);  
132.     if c > 1  
133.         for l=1:1:c  
134.             index = validEcacheList(l)  
135.             if index == k  
136.                 continue;  
137.             end  
138.             Ek = calcEk(oS,index);  
139.             deltaE = abs(Ei - Ek);  
140.             if(deltaE > maxDeltaE)  
141.                 maxK = index;  
142.                 maxDeltaE = deltaE;  
143.                 Ej = Ek;  
144.             end  
145.         end  
146.         j = maxK;  
147.     else  
148.         j = selectJrand(k, oS.m);  
149.         Ej = calcEk(oS, j);  
150.     end  
151. end  
152.   
153. function index = selectJrand(k,m)  
154.     index = k;  
155.     while(index == k)  
156.         index = randi([1,m],1,1);    
157.     end  
158. end  
159.   
160. function res = clipalpha(a, H, L)  
161.     if a > H  
162.         a = H;  
163.     end  
164.       
165.     if a < L  
166.         a = L;  
167.     end  
168.     res = a;  
169. end  

[plain] view plain copy

1. clc;  
2. clear;  
3.   
4. load Data  
5.   
6. [r,c] = size(Data);  
7. Test = Data(:,1:2);  
8. Label = Data(:,3);  
9.   
10. [b, alphas] = smoP(Test, Label, 0.6, 0.001, 40);  
11.   
12. %%画图  
13. figure(1)  
14. axis([-2 12 -8 6])  
15. for k = 1:1:r  
16.     hold on  
17.     if Data(k,3) == 1  
18.         plot(Data(k,1),Data(k,2),'r+');  
19.     else  
20.         plot(Data(k,1),Data(k,2),'b*');  
21.     end  
22. end  
23.   
24. %画支持向量及分割面  
25. %result=[];  
26. for k=1:1:r  
27.     if alphas(k)~= 0  
28.         hold on  
29.         %result =[result;alphas(k)];  
30.         QX = plot(Data(k,1:1),Data(k,2:2),'Ok','MarkerSize',12);  
31.         set(QX,'LineWidth',2.0);  
32.     end  
33. end  
34. W=(alphas.*Label)'*Data(:,1:2);  
35. y=(-W(1).* Data(:,1:1)-b) ./W(2);  
36. plot(Data(:,1:1),y);  

运行结果：

与第一个代码唯一的不同就是选择alphas的方式。第一个代码覆盖了所有数据集。常数C=0.6，一方面要保证所有样例的间隔不小于1.0，另一方面又要使得分类间隔尽可能大，并且要在 这两方面平衡。如果C很大，那么分类器将力图通过分割超平面对所有的样例都正确分类。小圆点标注的是支持向量。如果数据集非线性可分，支持向量 会在超平面附近聚集成团

四、非线性可分问题

对于上图，在二维平面中很难用直线分割，但是这里明显存在着两类数据。接下来，我们就使用一种称为核函数的工具将数据转化成易于分类器理解的形式。

1.      利用核函数将数据映射到高维空间

对于上图而言，如果只在x和y构成的坐标系中插入直线进行分类的话，我们不会得到理想的结果。我们可以对数据进行转换从而得到某些新的变量来表示数据。在这种情况下，我们就更容易得到大于零或小于零的测试结果。数学家们将数据从一个特征空间转换到另一特征空间的过程称为特征空间映射，通常我们将低维特征空间映射到高维特征空间。下面举个例子来形象地理解核函数：

我们把横轴上端点a和b之间红色部分里的所有点定为正类，两边的黑色部分里的点定为负类。试问能找到一个线性函数把两类正确分开么？不能，因为二维空间里的线性函数就是指直线，显然找不到符合条件的直线。但我们可以找到一条曲线，例如下面这一条：

    显然通过点在这条曲线的上方还是下方就可以判断点所属的类别（你在横轴上随便找一点，算算这一点的函数值，会发现负类的点函数值一定比0大，而正类的一定比0小）。这条曲线就是我们熟知的二次曲线。

上述过程即完成了一维空间向二维空间的映射。

对于SVM分类问题，所有的运算都可以写成内积形式(点积)，我们把内积运算替换成核函数，即可完成特征映射。核函数主要有：

l  多项式核

l  傅立叶核

l  B样条核

l  Sigmod核

l  高斯径向基核

核函数并不仅仅应用于支持向量机，很多其他机器学习算法也要用到。下面就介绍高斯径向基核函数。

径向基函数是一种采用向量作为自变量的函数，能够基于向量距离输出一个标量，具体数学公式：

     其中，σ是用户定义的用于确定到达率或者说是函数值跌落到0的速度参数。这个高斯核函数将数据从其特征空间映射到更高维的空间，具体说来这里是映射到一个无穷维的空间。我们不用确切地理解数据是如何表现的。

【这里扯一下我的同学，他的论文《基于矩阵运算的单隐层Madaline网络批量学习》，人家提出数据往低维空间映射，比较神奇哈】

最终的分类平面：（推导参考：http://blog.csdn.net/wangran51/article/details/7354915http://blog.csdn.net/wangran51/article/details/7354915）

代码：

[plain] view plain copy

1. function [b, res_alphas] = rbf_smoP(data, class, C, toler, maxIter, k1)  
2.     [m,n] = size(data);  
3.     iter = 0;  
4.     entireSet = 1;  
5.     alphaPairsChanged = 0;  
6.     oS = init(data, class, C, toler, m, k1);  
7.           
8.     while(((iter<maxIter)&&(alphaPairsChanged > 0)) || (entireSet == 1))  
9.         alphaPairsChanged = 0;  
10.         if entireSet == 1  
11.             for k = 1:1:m  
12.                 [ret, oS] = innerL(k, oS);  
13.                 alphaPairsChanged = alphaPairsChanged + ret;  
14.             end  
15.             iter = iter + 1;  
16.         else  
17.             nonBoundIs = [];  
18.             for k = 1:1:m  
19.                if ((oS.alphas(k) < C) && (oS.alphas(k) > 0))  
20.                    nonBoundIs = [nonBoundIs k];  
21.                end  
22.             end  
23.             [r,c] = size(nonBoundIs);  
24.             for k = 1:1:c  
25.                 index = nonBoundIs(k);  
26.                 [ret, oS] = innerL(index, oS);  
27.                 alphaPairsChanged = alphaPairsChanged + ret;  
28.             end  
29.             iter = iter + 1;  
30.         end  
31.         if entireSet == 1  
32.             entireSet = 0;  
33.         elseif alphaPairsChanged == 0  
34.             entireSet = 1;  
35.         end  
36.     end  
37.     b = oS.b;  
38.     res_alphas = oS.alphas;  
39. end  
40.   
41. function K = kernelTrans(X, A, k1)  
42.     [m, n] = size(X);  
43.     K = zeros(m,1);  
44.     for j = 1:1:m  
45.         deltaRow = X(j,:) - A;  
46.         K(j) = deltaRow * deltaRow';  
47.     end  
48.     K = exp(K./(-2*k1));  
49. end  
50.   
51. function oS = init(data,class,C,toler,m,k1)  
52.     alphas = zeros(m,1);  
53.     eCache = zeros(m,2);  
54.     b = 0;  
55.       
56.     oS.data = data;  
57.     oS.class = class;  
58.     oS.C = C;  
59.     oS.toler = toler;  
60.     oS.m = m;  
61.     oS.alphas = alphas;  
62.     oS.b = b;  
63.     oS.eCache = eCache;  
64.     oS.K = zeros(m,m);  
65.     for j = 1:1:m  
66.         oS.K(:,j) = kernelTrans(oS.data,oS.data(j,:),k1);  
67.     end  
68. end  
69.   
70. function [ret,oS] = innerL(k, oS)  
71.     Ei = calcEk(oS, k);  
72.     if(((oS.class(k)*Ei < -oS.toler) && (oS.alphas(k) < oS.C)) || ((oS.class(k)*Ei > oS.toler) && (oS.alphas(k) > 0)))  
73.         [j, Ej] = selectJ(k, oS, Ei);  
74.         temp_k = oS.alphas(k);  
75.         temp_j = oS.alphas(j);  
76.           
77.         if oS.class(k) ~= oS.class(j)  
78.             L = max(0, oS.alphas(j) - oS.alphas(k));  
79.             H = min(oS.C, oS.C +oS.alphas(j) - oS.alphas(k));  
80.         else  
81.             L = max(0, oS.alphas(j) + oS.alphas(k) - oS.C);  
82.             H = min(oS.C, oS.alphas(j) + oS.alphas(k));  
83.         end  
84.         if L == H  
85.             ret = 0;  
86.             return;  
87.         end  
88.         eta = 2.0 * oS.K(k,j) - oS.K(k,k) - oS.K(j,j);  
89.         if eta >=0   
90.             ret = 0;  
91.             return;  
92.         end  
93.         oS.alphas(j) = oS.alphas(j) - oS.class(j) * (Ei - Ej) / eta;  
94.         oS.alphas(j) = clipalpha(oS.alphas(j), H, L);  
95.           
96.         %update Ek  
97.         Et = calcEk(oS, j);  
98.         oS.eCache(j,:) = [1 Et];  
99.           
100.         if(abs(oS.alphas(j) - temp_j) < 0.00001)  
101.             ret = 0;  
102.             return;  
103.         end  
104.           
105.         oS.alphas(k) =   oS.alphas(k) + oS.class(j)*oS.class(k)*(temp_j - oS.alphas(j));  
106.         Et = calcEk(oS, k);  
107.         oS.eCache(k,:) = [1 Et];  
108.           
109.         b1 = oS.b - Ei - oS.class(k) * (oS.alphas(k) - temp_k) * oS.K(k,k) - oS.class(j) * (oS.alphas(j) - temp_j) * oS.K(k,j);  
110.         b2 = oS.b - Ej - oS.class(k) * (oS.alphas(k) - temp_k) * oS.K(k,j) - oS.class(j) * (oS.alphas(j) - temp_j) * oS.K(j,j);   
111.           
112.         if (oS.alphas(k)>0) && (oS.alphas(k)<oS.C)  
113.             oS.b = b1;  
114.         elseif(oS.alphas(j)>0) && (oS.alphas(j)<oS.C)  
115.             oS.b = b2;  
116.         else  
117.             oS.b = (b1+b2)/2.0;  
118.         end  
119.         ret = 1;  
120.         return;  
121.     else  
122.         ret = 0;  
123.         return;  
124.     end  
125. end  
126.   
127. function Ek = calcEk(oS, k)  
128.     fXk = (oS.alphas .* oS.class)' * oS.K(:,k) + oS.b;  
129.     Ek = fXk - oS.class(k);  
130. end  
131.   
132. function [j, Ej] = selectJ(k, oS, Ei)  
133.     maxK = -1;  
134.     maxDeltaE = 0;   
135.     Ej = 0;  
136.     oS.eCache(k,:) =[1 Ei];  
137.     validEcacheList = [];  
138.       
139.     for l = 1:1:oS.m  
140.         if oS.eCache(l,1:1) ~= 0  
141.             validEcacheList = [validEcacheList l];  
142.         end  
143.     end  
144.     [r, c] = size(validEcacheList);  
145.     if c > 1  
146.         for l=1:1:c  
147.             index = validEcacheList(l);  
148.             if index == k  
149.                 continue;  
150.             end  
151.             Ek = calcEk(oS,index);  
152.             deltaE = abs(Ei - Ek);  
153.             if(deltaE > maxDeltaE)  
154.                 maxK = index;  
155.                 maxDeltaE = deltaE;  
156.                 Ej = Ek;  
157.             end  
158.         end  
159.         j = maxK;  
160.     else  
161.         j = selectJrand(k, oS.m);  
162.         Ej = calcEk(oS, j);  
163.     end  
164. end  
165.   
166. function index = selectJrand(k,m)  
167.     index = k;  
168.     while(index == k)  
169.         index = randi([1,m],1,1);    
170.     end  
171. end  
172.   
173. function res = clipalpha(a, H, L)  
174.     if a > H  
175.         a = H;  
176.     end  
177.       
178.     if a < L  
179.         a = L;  
180.     end  
181.     res = a;  
182. end  

[plain] view plain copy

1. clc;  
2. clear;  
3.   
4. load NData  
5. load NTest  
6.   
7. Data = ndata;  
8. Data_Test = ntest;  
9. [r,c] = size(Data);  
10. Test = Data(:,1:2);  
11. Label = Data(:,3);  
12.   
13. [b, alphas] = rbf_smoP(Test, Label, 200, 0.0001, 1000,1.3);  
14.   
15. %%画图  
16. figure(1)  
17. axis([-1.5 1.5 -1.5 1.5])  
18. for k = 1:1:r  
19.     hold on  
20.     if Data(k,3) == 1  
21.         plot(Data(k,1),Data(k,2),'r+');  
22.     else  
23.         plot(Data(k,1),Data(k,2),'b*');  
24.     end  
25. end  
26. %%画支持向量  
27. support_vector = [];  
28. lable_sv = [];  
29. alphas_sv = [];  
30. for k=1:1:r  
31.     if alphas(k)~= 0  
32.         hold on  
33.         support_vector = [support_vector; Test(k,1:2)];  
34.         lable_sv = [lable_sv Label(k)];  
35.         alphas_sv = [alphas_sv alphas(k)];  
36.         %result =[result;alphas(k)];  
37.         QX = plot(Data(k,1:1),Data(k,2:2),'Ok','MarkerSize',12);  
38.         set(QX,'LineWidth',2.0);  
39.     end  
40. end  
41. %%预测  
42. temp = lable_sv .* alphas_sv;  
43. [m, n] = size(Data_Test);  
44. errorCount = 0;  
45. for k = 1:1:m  
46.     value = kernelTrans(support_vector, Data_Test(k,1:2),1.3);  
47.     predict = temp * value + b;  
48.     if predict > 0  
49.         predict = 1;  
50.     else  
51.         predict = -1;  
52.     end  
53.     if predict ~= Data_Test(k,3:3)  
54.         errorCount = errorCount + 1;  
55.     end  
56. end  
57. errorCount  

运行结果：

支持向量围绕超平面成团了。。。

预测结果，错分类2，效果不错。

代码地址：

http://download.csdn.net/detail/jinshengtao/8134089

原址

http://blog.csdn.net/jinshengtao/article/details/40900865