SVM回顾
支持向量机(SVM)的一大特点是最大化间距(max margin)。对于如上图的二分类问题,虽然有很多线可以将左右两部分分开,但是只有中间的红线效果是最好的,因为它的可活动范围(margin)是最大的,从直观上来说很好理解。
对于线性二分类问题,假设分类面为
则margin为
根据max margin规则和约束条件,得到如下优化问题,我们要求的就是参数和
:
对于正样本,类标号为+1,反之则为-1。根据拉格朗日对偶,(3)可以转换为如下的二次规划(QP)问题,其中
为拉格朗日乘子。
其中N为样本数量。上式还需满足如下两个约束条件:
一旦求解出所有的拉格朗日乘子,则我们可以通过如下的公式得到分类面参数和
。
当然并不是所有的数据都可以完美的线性划分,可能有少量数据就是混在对方阵营,这时可以通过引入松弛变量得到软间隔形式的SVM:
其中的为松弛变量,能假装把错的样本分对,
对max margin和margin failures的trades off。对于这个新的优化问题,约束变成了一个box constraint:
而松弛变量不再出现在对偶公式中了。
对于线性不可分的数据,可以用核函数将其投影到高维空间,这样就可分了,由此得到一般的分类面公式:
终极优化问题就变成了下面这个样子:
要满足的KKT条件为:
SMO算法
为了求解式(11),SMO算法通过启发式方法选择两个当变量,固定其他
,然后用解析的方法求解两个变量的二次规划问题。关于这两个变量的解应该更接近原始的二次规划问题,更重要的是,这时子问题可以通过解析方法求解,避免了矩阵运算,大大提高了整个算法的计算速度。如此,SMO算法将原问题不断分解为子问题并对子问题求解,进而达到求解原问题的目的。
不失一般性,假设选择的两个拉格朗日乘子是,因为
和变量
是一一对应的,所以也可以说选择了变量
。固定其他变量
此时,式(11)的优化问题转换为如下的优化问题:
满足
因为只有是变量,其他
是固定的(可以认为是常数),所以式(13)省略了不含
的常数项。又因为
,
有限制(9),所以
被限制在
的盒子里,且位于对角线上,如图Figure 1.
由于需满足(9),所以最优值
的取值范围必须满足条件
其中L与H是所在的对角线段端点的界。如果
(Figure 1.左图),则
如果(Figure 1.右图),则
令
因为,所以
,式(14)同乘以
,用
表示
并带入式(13),得到只含一个变量
的表达式,该表达式对
求偏导并等于0,求出
的更新公式为:
其中是第
个样本的训练误差。经过截断之后的
为:
将的结果代入式(14),得到
的更新公式:
当时,由式(12)中间的条件可以得到阈值
的更新公式:
同理,当时,得到阈值
的更新公式:
综合以上,阈值的更新公式为:
至此,如果我们选定了两个变量,则可以通过公式(19)和(21)来更新它们,并通过公式(24)更新阈值
,通过多次迭代逼近最优结果。但是怎样选择
呢,SMO通过启发式方法选择!
对于第一个样本,我们选择违反KKT条件式(12)(最严重)的样本
;由式(19)可知,
是依赖于
的,为了加快计算,一种简单的做法是选择
,使其对应的
最大。
至此,我们有了启发式选择两个变量的方法,以及求解这两个变量二次规划的解析方法,下面介绍SMO算法的具体实现。
MATLAB代码实现
算法伪代码如下:
输入:训练数据,其中
,
,
,精度为
。
输出:拉格朗日乘子和阈值
的近似解。
- 初始化
- 选一个违反KKT条件的变量
- 选一个使得
最大的变量
- 根据公式(19)~(21)更新变量
- 根据公式(22)~(24)更新阈值
- 如果不满足结束条件,则转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