SVM回顾
支持向量机(SVM)的一大特点是最大化间距(max margin)。对于如上图的二分类问题,虽然有很多线可以将左右两部分分开,但是只有中间的红线效果是最好的,因为它的可活动范围(margin)是最大的,从直观上来说很好理解。
对于线性二分类问题,假设分类面为
则margin为
根据max margin规则和约束条件,得到如下优化问题,我们要求的就是参数$\vec w$和$b$:
对于正样本,类标号$y_i$为+1,反之则为-1。根据拉格朗日对偶,(3)可以转换为如下的二次规划(QP)问题,其中$\alpha_i$为拉格朗日乘子。
其中N为样本数量。上式还需满足如下两个约束条件:
一旦求解出所有的拉格朗日乘子,则我们可以通过如下的公式得到分类面参数$\vec w$和$b$。
当然并不是所有的数据都可以完美的线性划分,可能有少量数据就是混在对方阵营,这时可以通过引入松弛变量$\xi_i$得到软间隔形式的SVM:
其中的$\xi_i$为松弛变量,能假装把错的样本分对,$C$对max margin和margin failures的trades off。对于这个新的优化问题,约束变成了一个box constraint:
而松弛变量$\xi_i$不再出现在对偶公式中了。
对于线性不可分的数据,可以用核函数$K$将其投影到高维空间,这样就可分了,由此得到一般的分类面公式:
终极优化问题就变成了下面这个样子:
要满足的KKT条件为:
SMO算法
为了求解式(11),SMO算法通过启发式方法选择两个$\alpha_i,\alpha_j$当变量,固定其他$\alpha_k$,然后用解析的方法求解两个变量的二次规划问题。关于这两个变量的解应该更接近原始的二次规划问题,更重要的是,这时子问题可以通过解析方法求解,避免了矩阵运算,大大提高了整个算法的计算速度。如此,SMO算法将原问题不断分解为子问题并对子问题求解,进而达到求解原问题的目的。
不失一般性,假设选择的两个拉格朗日乘子是$\alpha_1,\alpha_2$,因为$\alpha_i$和变量$x_i$是一一对应的,所以也可以说选择了变量$x_1,x_2$。固定其他变量$\alpha_i (i=3,4,...,N)$
此时,式(11)的优化问题转换为如下的优化问题:
满足
因为只有$\alpha_1,\alpha_2$是变量,其他$\alpha_i (i=3,4,...,N)$是固定的(可以认为是常数),所以式(13)省略了不含$\alpha_1,\alpha_2$的常数项。又因为$y_i=\pm 1$,$\alpha_i$有限制(9),所以$\alpha_1,\alpha_2$被限制在$[0,C]\times[0,C]$的盒子里,且位于对角线上,如图Figure 1.
由于$\alpha_2^{new}$需满足(9),所以最优值$\alpha_2^{new}$的取值范围必须满足条件
其中L与H是$\alpha_2^{new}$所在的对角线段端点的界。如果$y_1\neq y_2$(Figure 1.左图),则
如果$y_1= y_2$(Figure 1.右图),则
令
因为$y_i=\pm 1$,所以$y_i^2=1$,式(14)同乘以$y_1$,用$\alpha_2$表示$\alpha_1$并代入式(13),得到只含一个变量$\alpha_2$的表达式,该表达式对$\alpha_2$求偏导并等于0,求出$\alpha_2$的更新公式为:
其中$E_i=u_i-y_i$是第$i$个样本的训练误差。经过截断之后的$\alpha_2^{new}$为:
将$\alpha_2^{new,clipped}$的结果代入式(14),得到$\alpha_1$的更新公式:
当$0<\alpha_1^{new}<C$时,由式(12)中间的条件可以得到阈值$b_1$的更新公式:
同理,当$0<\alpha_2^{new,clipped}<C$时,得到阈值$b_2$的更新公式:
综合以上,阈值$b$的更新公式为:
至此,如果我们选定了两个变量$\alpha_2,\alpha_1$,则可以通过公式(19)和(21)来更新它们,并通过公式(24)更新阈值$b$,通过多次迭代逼近最优结果。但是怎样选择$\alpha_2,\alpha_1$呢,SMO通过启发式方法选择!
对于第一个样本$\alpha_2$,我们选择违反KKT条件式(12)(最严重)的样本$\alpha_2$;由式(19)可知,$\alpha_2^{new}$是依赖于$|E_1-E_2|$的,为了加快计算,一种简单的做法是选择$\alpha_1$,使其对应的$|E_1-E_2|$最大。
至此,我们有了启发式选择两个变量$\alpha_2,\alpha_1$的方法,以及求解这两个变量二次规划的解析方法,下面介绍SMO算法的具体实现。
MATLAB代码实现
算法伪代码如下:
输入:训练数据$T=\{(\vec x_1,y_1),(\vec x_2,y_2),...,(\vec x_N,y_N)\}$,其中$\vec x_i \in R^n$,$y_i\in\{-1,+1\}$,$i=1,2,...,N$,精度为$\epsilon$。
输出:拉格朗日乘子$\vec\alpha=\{\alpha_1,\alpha_2,...,\alpha_N\}$和阈值$b$的近似解。
- 初始化$\vec \alpha=\vec 0, b=0$
- 选一个违反KKT条件的变量$\alpha_2$
- 选一个使得$|E_1-E_2|$最大的变量$\alpha_1$
- 根据公式(19)~(21)更新变量$\alpha_2,\alpha_1$
- 根据公式(22)~(24)更新阈值$b$
- 如果不满足结束条件,则转2;否则算法结束
算法的结束条件是:如果所有变量都已经检查过,且没有变量可以进一步优化时,算法结束。
下面是SMO算法的MATLAB简化实现,省略了论文[1]的公式(19)。代码主流程参考论文[1]中的伪代码,部分实现细节参考[2]。
本代码和MATLAB自带的seqminopt.m算法接口相同,可直接替换使用。下载代码。
function [alphas, offset] = bitjoy_seqminopt(data, targetLabels, boxConstraints, ...
kernelFunc, smoOptions)
% http://bitjoy.net/2016/05/02/svm-smo-algorithm/
% 参考文献:
% [1] A Fast Algorithm for Training Support Vector Machines,
% http://research.microsoft.com/pubs/69644/tr-98-14.pdf
% [2] CSDN, http://blog.csdn.net/techq/article/details/6171688
N = length(data); % 样本数量
alphas = zeros([N,1]); % 所有拉格朗日乘子初始化为0
offset = 0.0; % 偏移量b也初始化为0
numChanged = 0; % 拉格朗日乘子改变的个数
examineAll = 1; % 是否需要检查所有拉格朗日乘子
% 主流程
while numChanged > 0 || examineAll
numChanged = 0;
if examineAll == 1
for i = 1 : N
numChanged = numChanged + examineExample(i);
end
else
for i = 1 : N
if alphas(i) ~= 0 && alphas(i) ~= boxConstraints(i)
numChanged = numChanged + examineExample(i);
end
end
end
if examineAll == 1
examineAll = 0;
elseif numChanged == 0
examineAll = 1;
end
end
% 计算当前参数下第idx个样本的函数输出
function [svm_o] = wxpb(idx)
svm_o = 0.0;
for j = 1 : N
svm_o = svm_o + alphas(j) * targetLabels(j) * kernelFunc(data(j,:),data(idx,:));
end
svm_o = svm_o + offset;
end
% 选定第二个变量i2后,根据max|E1-E2|的启发式规则,选择i1;
% 如果没有满足条件的i1,返回-1.
function [i1] = selectSecondChoice(i2, E2)
i1 = -1;
maxDelta = -1;
for j = 1 : N
if j ~= i2 && alphas(j) ~= 0 && alphas(j) ~= boxConstraints(j)
Ej = wxpb(j) - targetLabels(j);
if abs(E2 - Ej) > maxDelta
i1 = j;
maxDelta = abs(E2 - Ej);
end
end
end
end
% 根据选定的两个变量i1,i2,代入更新公式计算;
% 最后还更新了偏移量offset,也就是y=wx+b中的b。
function [flag] = takeStep(i1, i2)
alpha1 = alphas(i1);
y1 = targetLabels(i1);
E1 = wxpb(i1) - y1;
alpha2 = alphas(i2);
y2 = targetLabels(i2);
E2 = wxpb(i2) - y2;
s = y1 * y2;
if y1 ~= y2
L = max(0, alpha2 - alpha1);
H = min(boxConstraints(i2), boxConstraints(i1) + alpha2 - alpha1);
else
L = max(0, alpha2 + alpha1 - boxConstraints(i1));
H = min(boxConstraints(i2), alpha2 + alpha1);
end
if L == H
flag = 0;
return;
end
k11 = kernelFunc(data(i1,:),data(i1,:));
k12 = kernelFunc(data(i1,:),data(i2,:));
k22 = kernelFunc(data(i2,:),data(i2,:));
eta = k11 + k22 - 2 * k12;
if eta > 0
a2 = alpha2 + y2 * (E1 - E2) / eta;
%截断
if a2 < L
a2 = L;
elseif a2 > H
a2 = H;
end
else % 并未处理该case,论文[1]公式(19)有更详细的方法
flag = 0;
return;
end
if abs(a2 - alpha2) < eps * (a2 + alpha2 + eps)
flag = 0;
return;
end
a1 = alpha1 + s * (alpha2 - a2);
alphas(i1) = a1;
alphas(i2) = a2;
% 更新offset
b1 = offset - E1 - y1 * k11 * (a1 - alpha1) - y2 * k12 * (a2 - alpha2);
b2 = offset - E2 - y1 * k12 * (a1 - alpha1) - y2 * k22 * (a2 - alpha2);
if a1 > 0 && a1 < boxConstraints(i1)
offset = b1;
elseif a2 > 0 && a2 < boxConstraints(i2)
offset = b2;
else
offset = (b1 + b2) / 2;
end
flag = 1;
end
% 检查样本。
% 首先判断i2是否满足KKT条件,如果不满足,
% 则根据启发式规则再选择i1样本,
% 然后更新i1和i2的拉格朗日乘子。
function [flag] = examineExample(i2)
y2 = targetLabels(i2);
alpha2 = alphas(i2);
E2 = wxpb(i2) - y2;
r2 = E2 * y2;
if (r2 < -smoOptions.TolKKT && alpha2 < boxConstraints(i2)) || (r2 > smoOptions.TolKKT && alpha2 > 0)
i1 = selectSecondChoice(i2, E2);
if i1 == -1
i1 = floor(1 + rand() * N); % 随机选一个i1
while i1 == i2
i1 = floor(1 + rand() * N);
end
flag = takeStep(i1,i2);
else
flag = takeStep(i1,i2);
end
else
flag = 0;
end
end
end
参考资料
[1]. J.C. Platt: A Fast Algorithm for Training Support Vector Machines