转自：http://www.cnblogs.com/yuwenchao/archive/2011/10/28/csvc.html
本文讨论对于一般的分类问题，线性支持向量分类机(C-SVC)的实现原理及MATLAB代码实现。
由于需要处理线性不可分问题，我们没法找到一个超平面可以完全正确的分化训练集，因此需要“软化”一些条件。由标准的支持向量机(SVM)的最大间隔法所导出的约束条件可以进行如下的“软化”：
yi((w·xi)+b) ≥ 1-ξi , i = 1, 2, ... , n
当ξi足够大时，训练点(yi, xi)总是可以满足条件的，所以我们不能让ξi太大，因此得给ξi一个惩罚系数，所以C-SVC的原问题可以归纳如下：
minω,b,ξ 1/2||ω||2 + CΣξi
s.t. yi((w·xi)+b) ≥ 1-ξi , i = 1, 2, ... , n
ξi ≥ 0, i = 1, 2, ... , n
根据原问题，我们可以得到如下的对偶问题：
minα 1/2αTHα-eTα
s.t. αTy = 0
0 ≤ α ≤ C
其中，H = ΣΣyiyj(xi·xj), i = 1, 2, ... , n, j = 1, 2, ... , n，α = (α1, α2, ... , αn)，y = (y1, y2, ... , yn)。
对于得到的对偶问题，我们可以用函数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的值。由于在求对偶问题的过程中，我们知道ω = Σαiyixi , i = 1, 2, ... , n。所以b = yj-ωTxj。在这里需要注意的是α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. 对于线性不可分问题的最大间隔法有理论基础吗？因为此时已经没有间隔了...
