机器学习之四：支持向量机，推导及实现
一、支持向量机（SVM）
1、哪个模型更优
先来讨论一下下面二分类问题。
数据样本如下图所示：
其可能的划分边界可能是如下图：
图中的蓝、红两种分类边界，都能达到分类的效果，但是哪一个效果更优？
评判一个模型的优劣，最主要的，是其对样本的泛化能力。
若有一个新样本，其二维特征位置，如下图蓝色方块所示：
两个模型，所得到的结果，显示是各不相同的，但从聚集的情况下来，该样本更偏向于“圆”的分类。
也就是说，蓝色边界所表示的模型，更有代表性。而它，也是SVM算法，所得到的模型。
2、支持向量机
支持向量机，是一个二分类模型，其基本思想是：找到集合边缘上的若干数据（称为支持向量），用这些点找出一个平面（称为决策面），使得支持向量到该平面的距离最长。
以上面的样本为例，SVM以下图圈起来的三个样本为支持向量，根据这三个样本，找到分类的边界——即中间蓝线表示的平面。
综上所述，虽然还未涉及SVM的具体运算，但可以初步看出，支持向量机有如下优点：
(1) 计算量少。只使用了支持向量进行计算，而其他距离边界很远的点，其特征很明显，不会引起歧义，不用参与计算。
(2) 泛化能力好。因其定义的边界，距离正负样本距离是最长的，具备良好的区分效果。
二、SVM推导
根据数据集的属性，可以将这个问题分为两个层次：
-
线性可分：数据分布简单，如上面的例子。可以找到一个超平面，直接在原始空间中将数据进行切分。
-
线性不可分：数据分布复杂，不存在这样的超平面，则通过将原始空间映射到一个高维空间，在高维空间对数据进行划分。
1、线性可分决策面
线性可分是一个最简单的分类问题，下面做推导。
1.1、首先是样本集
\[
{(x_1,y_1),(x_2,y_2),…,(x_n,y_n)}
\]
其中，\(y\) 的取值为 {+1,−1}。
1.2、决策面方程
从第一节知道，分类的手段是取到一个符合条件的决策平面。
平面可以用方程表示：
\[
w^Tx + b = 0
\]
那么问题的解决，则转化为找到符合条件的 \(w\) 和 \(b\)，条件是：使得数据集的边缘若干点，到这个平面的距离是最长的。
该平台将样本切分为两个部分，故：
\(y = -1\) 的样本点可表示为：\(w^Tx + b < 0\)
\(y = +1\) 的样本点可表示为：\(w^Tx + b > 0\)
那么支持向量所在的平面可以表示为：\(w^Tx + b = \pm A\)
在后面的优化中，其实\(A\)的值为多少，并不影响结果，为了方便计算，令 \(A = 1\)，即：
\[w^Tx + b = \pm 1\]
如下图所示：
1.3、最长距离
由几何知识可知，点到平面的距离可以表示为：
\[
\gamma = \frac{w^Tx+b}{{||w||}}
\]
此距离称为几何距离，存在正负。
而算法的目标，就是找到一个集合 \(x\) (支持向量)，使得 \(\gamma\) 的值最小。因 \(y\) 的取值为{+1,−1}，所以可以乘上一个 \(y\)，使该距离恒为正（前提该样本点 \(x\) 被正确分类的情况下）。
故有最大距离为
\[
max({\frac{y(w^Tx+b)}{{||w||}}})
\]
结上一节的假设，支持向量所在方程：
\[
w^Tx + b = \pm 1
\]
故而
\[
y(w^Tx+b) = 1
\]
则最大距离简化为
\[
max({\frac{1}{{||w||}}})
\]
即，等同于求解下式的最小化值
\[
min({\frac{1}{{2}}}||w||^2)
\]
这里增加了一个1/2系数、和一个平方，其实是为了方便求导。一求导两者就相消了。
这个式子有还有一些限制条件，完整的写下来，应该是这样的：
\[
min({\frac{1}{{2}}}||w||^2),s.t,y(w^Tx+b)\geq1
\]
s.t.后面的限制条件可以看做是一个凸多面体，我们要做的就是在这个凸多面体中找到最优解。
1.4、求最优解
上面优化问题，可以用拉格朗日乘子法去解，使用了KKT条件的理论，这里直接作出这个式子的拉格朗日目标函数：
\[
L(w,b,a) =\frac{1}{{2}}||w||^2-\sum_{i=1}^n \alpha_i(y_i(w^Tx+b)-1)
\]
目标是让 $L(ω,b,α) $ 针对 $ α $ 达到最大值。如何求解？
L是关于 \(w,b,α\) 三个变量的函数，要得到使L最大值的 \(α\)，需进行两步操作：
-
首先需要先排除掉 \(w\) 和 \(b\) 的影响，让\(L\)关于\(w、b\) 最小化，即不管 \(w、b\) 如何改变，\(L\) 都不会再变小
-
接着再让 \(L\) 关于 \(α\) 取最大值
(1) \(L\) 关于 \(w、b\) 最小化
分别令 \(L\) 关于 \(w、b\) 的偏导数为0，得到关于原问题的一个表达式：
令：
\[
\frac{\partial L}{\partial w} = 0
\]
得到：
\[
w=\sum_{i=1}^n\alpha_iy_ix_i
\]
令：
\[
\frac{\partial L}{\partial b} = 0
\]
得到：
\[
\sum_{i=1}^n\alpha_iy_i = 0
\]
将上两式带回 $L(w,b,a) $得到对偶问题的表达式
\[
L(w,b,a) =\frac{1}{{2}}\sum_{i=1}^n\alpha_i – \frac{1}{{2}}\sum_{i,j=1}^n\alpha_i\alpha_jy_iy_jx^T_ix_j
\]
于是新的目标问题及限制条件为（对偶问题）:
\[
max_a(\sum_{i=1}^n\alpha_i – \frac{1}{{2}}\sum_{i,j=1}^n\alpha_i\alpha_jy_iy_jx^T_ix_j)
\]
\[
s.t.,\alpha \geq0,i=1,2,…,n
\]
\[
\sum_{i=1}^n\alpha_iy_i = 0
\]
这个就是我们需要最终优化的式子，只关于 \(α\) 向量的式子。
(2) L关于 \(α\) 的最大值
上式最终的对偶问题，是一个凸二次规划问题，理论上用任何一个解决凸二次规划的软件包都可以解决。一般使用SMO算法。
SMO基本思想是，不断执行如下两个步骤，直至收敛：
-
选取一对参数\((\alpha_i,\alpha_j)\)
-
固定 \(\alpha\) 向量的其他参数，将\((\alpha_i,\alpha_j)\)代入上述表达式进行求最优解获得更新后的\((\alpha_i,\alpha_j)\)
具体解法，属于另一个层面的东西，就不在此处细述。
总之，通过SMO算法，能得到L的最大值，相当于得到 \(w\) 的值，进而可以算出 \(b\) ，从而得出需要的边界。
2、线性不可分决策面
上面讨论了线性可分的数据集的处理方式。但是，实际应用中的数据样本，可能更多的是线性不可分的，即不能找到一个可以将数据分类的超平台。
这种情况下，一般使用核函数将原始空间映射到一个高维空间，在高维高间对数据进行划分。理论上只要维度足够高，那么总能做到分类。
2.1 决策面方程
在线性可分的基础上，将样本 \(x\)，进行一次变换，得到\(\phi(x)\)
超平台变为：
\[
w^T\phi(x) + b = 0
\]
2.2 求最优解
整个推理过程，与线性可分的基本一样，唯一不同的，是将各个公式化中的 \(x\)，换成 \(\phi(x)\)。
即得到对遇问题：
\[
max_a(\sum_{i=1}^n\alpha_i – \frac{1}{{2}}\sum_{i,j=1}^n\alpha_i\alpha_jy_iy_j\phi(x_i)^T\phi(x_j))
\]
\[
s.t.,\alpha \geq0,i=1,2,…,n
\]
\[
\sum_{i=1}^n\alpha_iy_i = 0
\]
令：
\[
\phi(x_i)^T\phi(x_j) = K(x_i,x_j)
\]
这个式子所做的事情就是将线性的空间映射到高维的空间，K函数有很多种，下面是比较典型的两种：
- 多项式核
\[
K(x_i,x_j) = (x_i.x_j+1)^d
\]
- 高斯核
\[
K(x_i,x_j) = exp(-\frac{(x_i-x_j)^2}{2\sigma^2})
\]
高斯核是将原始空间映射为无穷维空间。
最后一样能通过各类软件包，例如SMO，实现求解。
三、程序实现
高斯核函数
function sim = gaussianKernel(x1, x2, sigma)
x1 = x1(:); x2 = x2(:);
sim = exp(-(x1 - x2)' * (x1 - x2) / (2 * sigma ^2));
endSVM训练函数
function [model] = svmTrain(X, Y, C, kernelFunction, ...
tol, max_passes)
if ~exist('tol', 'var') || isempty(tol)
tol = 1e-3;
end
if ~exist('max_passes', 'var') || isempty(max_passes)
max_passes = 5;
end
% Data parameters
m = size(X, 1);
n = size(X, 2);
% Map 0 to -1
Y(Y==0) = -1;
% Variables
alphas = zeros(m, 1);
b = 0;
E = zeros(m, 1);
passes = 0;
eta = 0;
L = 0;
H = 0;
% Pre-compute the Kernel Matrix since our dataset is small
% (in practice, optimized SVM packages that handle large datasets
% gracefully will _not_ do this)
%
% We have implemented optimized vectorized version of the Kernels here so
% that the svm training will run faster.
if strcmp(func2str(kernelFunction), 'linearKernel')
% Vectorized computation for the Linear Kernel
% This is equivalent to computing the kernel on every pair of examples
K = X*X';
elseif strfind(func2str(kernelFunction), 'gaussianKernel')
% Vectorized RBF Kernel
% This is equivalent to computing the kernel on every pair of examples
X2 = sum(X.^2, 2);
K = bsxfun(@plus, X2, bsxfun(@plus, X2', - 2 * (X * X')));
K = kernelFunction(1, 0) .^ K;
else
% Pre-compute the Kernel Matrix
% The following can be slow due to the lack of vectorization
K = zeros(m);
for i = 1:m
for j = i:m
K(i,j) = kernelFunction(X(i,:)', X(j,:)');
K(j,i) = K(i,j); %the matrix is symmetric
end
end
end
% Train
fprintf('\nTraining ...');
dots = 12;
while passes < max_passes,
num_changed_alphas = 0;
for i = 1:m,
% Calculate Ei = f(x(i)) - y(i) using (2).
% E(i) = b + sum (X(i, :) * (repmat(alphas.*Y,1,n).*X)') - Y(i);
E(i) = b + sum (alphas.*Y.*K(:,i)) - Y(i);
if ((Y(i)*E(i) < -tol && alphas(i) < C) || (Y(i)*E(i) > tol && alphas(i) > 0)),
% In practice, there are many heuristics one can use to select
% the i and j. In this simplified code, we select them randomly.
j = ceil(m * rand());
while j == i, % Make sure i \neq j
j = ceil(m * rand());
end
% Calculate Ej = f(x(j)) - y(j) using (2).
E(j) = b + sum (alphas.*Y.*K(:,j)) - Y(j);
% Save old alphas
alpha_i_old = alphas(i);
alpha_j_old = alphas(j);
% Compute L and H by (10) or (11).
if (Y(i) == Y(j)),
L = max(0, alphas(j) + alphas(i) - C);
H = min(C, alphas(j) + alphas(i));
else
L = max(0, alphas(j) - alphas(i));
H = min(C, C + alphas(j) - alphas(i));
end
if (L == H),
% continue to next i.
continue;
end
% Compute eta by (14).
eta = 2 * K(i,j) - K(i,i) - K(j,j);
if (eta >= 0),
% continue to next i.
continue;
end
% Compute and clip new value for alpha j using (12) and (15).
alphas(j) = alphas(j) - (Y(j) * (E(i) - E(j))) / eta;
% Clip
alphas(j) = min (H, alphas(j));
alphas(j) = max (L, alphas(j));
% Check if change in alpha is significant
if (abs(alphas(j) - alpha_j_old) < tol),
% continue to next i.
% replace anyway
alphas(j) = alpha_j_old;
continue;
end
% Determine value for alpha i using (16).
alphas(i) = alphas(i) + Y(i)*Y(j)*(alpha_j_old - alphas(j));
% Compute b1 and b2 using (17) and (18) respectively.
b1 = b - E(i) ...
- Y(i) * (alphas(i) - alpha_i_old) * K(i,j)' ...
- Y(j) * (alphas(j) - alpha_j_old) * K(i,j)';
b2 = b - E(j) ...
- Y(i) * (alphas(i) - alpha_i_old) * K(i,j)' ...
- Y(j) * (alphas(j) - alpha_j_old) * K(j,j)';
% Compute b by (19).
if (0 < alphas(i) && alphas(i) < C),
b = b1;
elseif (0 < alphas(j) && alphas(j) < C),
b = b2;
else
b = (b1+b2)/2;
end
num_changed_alphas = num_changed_alphas + 1;
end
end
if (num_changed_alphas == 0),
passes = passes + 1;
else
passes = 0;
end
fprintf('.');
dots = dots + 1;
if dots > 78
dots = 0;
fprintf('\n');
end
if exist('OCTAVE_VERSION')
fflush(stdout);
end
end
fprintf(' Done! \n\n');
% Save the model
idx = alphas > 0;
model.X= X(idx,:);
model.y= Y(idx);
model.kernelFunction = kernelFunction;
model.b= b;
model.alphas= alphas(idx);
model.w = ((alphas.*Y)'*X)';
end
绘制边界
function visualizeBoundary(X, y, model, varargin)
plotData(X, y)
% Make classification predictions over a grid of values
x1plot = linspace(min(X(:,1)), max(X(:,1)), 100)';
x2plot = linspace(min(X(:,2)), max(X(:,2)), 100)';
[X1, X2] = meshgrid(x1plot, x2plot);
vals = zeros(size(X1));
for i = 1:size(X1, 2)
this_X = [X1(:, i), X2(:, i)];
vals(:, i) = svmPredict(model, this_X);
end
% Plot the SVM boundary
hold on
contour(X1, X2, vals, [0.5 0.5], 'r','LineWidth',2);
hold off;
end
主调用方法
load('ex6data2.mat');
% SVM Parameters
C = 1; sigma = 0.1;
model= svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));
visualizeBoundary(X, y, model);