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

　　本文讨论对于一般的分类问题，线性支持向量分类机(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()对这个凸二次规划问题进行相应的求解，所得到的结果即为α的解。代码如下：

% 选定适当的惩罚系数
C = 1.95;

% 实验数据heartData的存储格式为第一列是正负类的标记，其他列为属性
% 根据具体情况对y和X进行赋值
y = heartData(:,1);
X = heartData(:,2:end);

% 记录数据的规模：样本数*属性维数
numbers = size(X);

% 求维数为“样本数*样本数”的矩阵
H = zeros(numbers(1),numbers(1));
for i = 1:numbers(1)
for j = 1:numbers(1)
H(i,j) = y(i)*y(j)*X(i,:)*X(j,:)';
end
end

% 以下是几个约束：
% Aeq 线性等式约束的矩阵
% beq 线性等式约束的向量
% f 向量
% lb 存储下界的向量
% ub 存储上界的向量
Aeq = y';
beq = 0;
f = -ones(numbers(1), 1);
lb = zeros(numbers(1), 1);
ub = C * ones(numbers(1), 1);

% 计算拉格朗日乘子，返回向量x，即拉格朗日乘子
[x,fval,exitflag,output,lambda] = quadprog(H,f,[],[],Aeq,beq,lb, ub, [], optimset('LargeScale','off','MaxIter', 500));

% 将x中小于一定阀值的数置为零（机器误差）
for j = 1:numbers(1)
if(x(j) < 0.001)
x(j) = 0;
end
end

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

% 计算向量w的值
w = zeros(1,numbers(2));
for i = 1:numbers(1)
w = w + y(i)*x(i)*X(i,:);
end

% 计算b的值（这里计算每一个y对应的b的值）
b = zeros(numbers(1), 1);
for i = 1:numbers(1)
if(x(i) > 0 && x(i) < C)
b(i) = y(i) - w*X(i,:)';
end
end
% 将b中小于一定阀值的数置为零（机器误差）
for j = 1:numbers(1)
if(abs(b(j)) < 0.001)
b(j) = 0;
end
end
% 由于计算的原因，b不一定总是相等，所以这里取b的值为最大和最小值的均值
bNonZero = find( b ~= 0 );
bval = 1/2*(max(b(bNonZero)) + min(b(bNonZero)));

%% 测试：通过已经求出来的判别函数对数据进行分类的错误率分析

% 记录分类正确的数目
correctNumbers = 0;
% 记录通过判别函数计算的结果
result = zeros(numbers(1),1);
for i = 1:numbers(1)
result(i) = w*X(i,:)' + bval;
if(((result(i) > 0) && (y(i) == 1))||((result(i) < 0)&&(y(i) == -1)))
correctNumbers = correctNumbers + 1;
end
end
% 输出正确率
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. 对于线性不可分问题的最大间隔法有理论基础吗？因为此时已经没有间隔了...