线性支持向量分类机及其实现

转自：http://www.cnblogs.com/yuwenchao/archive/2011/10/28/csvc.html

本文讨论对于一般的分类问题，线性支持向量分类机(C-SVC)的实现原理及MATLAB代码实现。

　　由于需要处理线性不可分问题，我们没法找到一个超平面可以完全正确的分化训练集，因此需要“软化”一些条件。由标准的支持向量机(SVM)的最大间隔法所导出的约束条件可以进行如下的“软化”：

　　    y_i((w·x_i)+b) ≥ 1-ξ_{i ,} i = 1, 2, ... , n

　　当ξ_i足够大时，训练点(y_i, x_i)总是可以满足条件的，所以我们不能让ξ_i太大，因此得给ξ_i一个惩罚系数，所以C-SVC的原问题可以归纳如下：

  min_ω,b,ξ  1/2||ω||² + CΣξ_i

　　　　　　 s.t.  y_i((w·x_i)+b) ≥ 1-ξ_{i ,} i = 1, 2, ... , n

                   ξ_i ≥ 0_, i = 1, 2, ... , n　　  　　

　　根据原问题，我们可以得到如下的对偶问题：

    min_α  1/2α^THα-e^Tα　　   　　

  s.t.  α^Ty = 0　　　　　　　

    0 ≤ α ≤ C　　　　

　　其中，H = ΣΣy_iy_j(x_i·x_j), i = 1, 2, ... , n, j = 1, 2, ... , n，α = (α₁, α₂, ... , α_n)，y = (y₁, y₂, ... , y_n)。

　　对于得到的对偶问题，我们可以用函数quadprog()对这个凸二次规划问题进行相应的求解，所得到的结果即为α的解。代码如下：

 1 % 选定适当的惩罚系数
 2 C = 1.95;
 3 
 4 % 实验数据heartData的存储格式为第一列是正负类的标记，其他列为属性
 5 % 根据具体情况对y和X进行赋值
 6 y = heartData(:,1);
 7 X = heartData(:,2:end);
 8 
 9 % 记录数据的规模：样本数*属性维数
10 numbers = size(X);
11 
12 % 求维数为“样本数*样本数”的矩阵
13 H = zeros(numbers(1),numbers(1));
14 for i = 1:numbers(1)
15     for j = 1:numbers(1)
16         H(i,j) = y(i)*y(j)*X(i,:)*X(j,:)';
17     end
18 end
19 
20 % 以下是几个约束：
21 % Aeq 线性等式约束的矩阵
22 % beq 线性等式约束的向量
23 % f 向量
24 % lb 存储下界的向量
25 % ub 存储上界的向量
26 Aeq = y';
27 beq = 0;
28 f = -ones(numbers(1), 1);
29 lb = zeros(numbers(1), 1);
30 ub = C * ones(numbers(1), 1);
31  
32 % 计算拉格朗日乘子，返回向量x，即拉格朗日乘子
33 [x,fval,exitflag,output,lambda] = quadprog(H,f,[],[],Aeq,beq,lb, ub, [], optimset('LargeScale','off','MaxIter', 500));
34  
35 % 将x中小于一定阀值的数置为零（机器误差）
36 for j = 1:numbers(1)
37     if(x(j) < 0.001)
38         x(j) = 0;
39     end
40 end

　　接下来我们需要计算b的值。由于在求对偶问题的过程中，我们知道ω = Σα_iy_ix_i , i = 1, 2, ... , n。所以b = y_j-ω^Tx_j。在这里需要注意的是α_j取(0, C)之间的数的时候，所得到的b是唯一的(具体证明在以后的随笔中会介绍)。因此我们可以得到划分的超平面 ω^Tx+b = 0。所以决策函数为 f(x) = sgn( ω^Tx+b)。求b以及决策函数的代码如下：

 1 % 计算向量w的值
 2 w = zeros(1,numbers(2));
 3 for i = 1:numbers(1)
 4     w = w + y(i)*x(i)*X(i,:);
 5 end
 6 
 7 % 计算b的值（这里计算每一个y对应的b的值）
 8 b = zeros(numbers(1), 1);
 9 for i = 1:numbers(1)
10     if(x(i) > 0 && x(i) < C)
11         b(i) = y(i) - w*X(i,:)';
12     end
13 end
14 % 将b中小于一定阀值的数置为零（机器误差）
15 for j = 1:numbers(1)
16     if(abs(b(j)) < 0.001)
17         b(j) = 0;
18     end
19 end
20 % 由于计算的原因，b不一定总是相等，所以这里取b的值为最大和最小值的均值
21 bNonZero = find( b ~= 0 );
22 bval = 1/2*(max(b(bNonZero)) + min(b(bNonZero)));
23  
24 %% 测试：通过已经求出来的判别函数对数据进行分类的错误率分析
25 
26 % 记录分类正确的数目
27 correctNumbers = 0;
28 % 记录通过判别函数计算的结果
29 result = zeros(numbers(1),1);
30 for i = 1:numbers(1)
31     result(i) = w*X(i,:)' + bval;
32     if(((result(i) > 0) && (y(i) == 1))||((result(i) < 0)&&(y(i) == -1)))
33         correctNumbers = correctNumbers + 1;
34     end
35 end
36 % 输出正确率
37 disp (correctNumbers/numbers(1))

　　上面代码中的测试正确率是这样进行计算的，将得到的决策函数在原始的训练集中跑一遍，因为是线性不可分的，所以不可能完全划分正确，我们将正确的数目除以总的样本数就得到一个简单正确率的计算，需要说明的是，要测试决策函数的好坏，判断C的选取，应该使用LOOCV等方法。

问题与思考：

1. 我曾采用的二维测试数据集为：

y = [1;1;1;1;1;-1;-1;-1;-1;-1];

X1 = [1.5 15;0.2 11;0.3 17;4 1;2 20;-10 -21;-2 -1;-15 -11;-1.4 -1;-2 -6.20];

X2 = [1.5 15;0.2 11;0.3 17;4 1;2 20;-10 -21;-2 -1;0.7 13;-1.4 -1;-2 -6.20];

其中X1是线性可分的情况、X2是线性不可分的情况。

对于数据集(y,X1)得到的结果如下图所示，有两个支持向量（第四个和第八个）。此时b不为零的取值为-0.4188和-0.4252。

对于数据集(y,X2)得到的结果如下图所示，有五个支持向量（第二、三、四、八和第九个）。此时b不为零的取值均为-0.4435。

所以我特别好奇为什么在前一个数据集上的测试中，得到的b不唯一。带着这个问题，我将不同的b的值都画出了相应的直线，得到了几乎重合的直线，因此我的解释是机器误差产生的，但是不知道这样解释靠不靠谱。这个问题还会继续去看看有没有办法解决。

2. 对于线性不可分问题的最大间隔法有理论基础吗？因为此时已经没有间隔了...