前言
在支持向量机这篇文章中详细叙述了支持向量机的理论推导部分，其中在实现支持向量机的部分，我们说到支持向量机最后就是求解一个二次规划问题，在那里我们使用的是MATLAB中求解二次规划的函数进行求解，这次将会采用一个新的算法来求解支持向量机的二次规划问题。
SMO算法介绍
SMO算法相比常规求解二次规划问题的算法，是一种非常快速的方法。在支持向量机的对偶问题中，有：
SMO算法实现
1.变量的选择方法
在学习SMO算法时，我是对照着《统计学习方法》这本书结合代码来学习的，在这里其实已经是第二遍了，第一遍时我只编出了代码，并没有给出详细的解释，这次打算从0实现一个SMO，顺便介绍我实现代码的思路。
变量的选择方法是SMO算法最巧妙的地方，在SMO算法中，如果对所有的N个变量进行优化，那种速度太慢，SMO采用的方法就是选择两个α进行优化，在选择的两个变量中，至少有一个变量是违反了KKT条件的。而这两个变量是一一进行选择的，不是一起把两个的变量都选择，而是一个一个选择，我们把第一个变量的选择称为外层循环，第二个变量的选择称为内层循环。
1.1 第一个变量的选择
在选择变量之前，我们需要对N个α赋初值，只有赋完初值后才能选择，初始化所有的α均为0，SMO在选择第一个α时主要看样本点中违反KKT条件最严重的点，让其作为第一个变量，如何判断样本点违反KKT条件的程度呢？
在KKT条件中有下面的式子：
1.2.第二个变量的选择
当第一个变量选取之后，我们就可以选取第二个变量了，如果说第一个变量选取的标准是违反KKT条件最大的点，那么第二个变量选取的标准就是能够使α2有很大的变化，所谓α变化最大，就是指在选择α2之后进行优化，我们新的 α2new改变幅度要大。至于如何判断α的变化幅度有多大，需要了解新的α2new是如何产生的，由于变量的解析方法在下一节介绍，这里只简单引入一下，那就是通过比较∣E1−E2∣的值来判断改变的幅度大小，而E又是什么呢？E的计算公式为：
SMO算法的实现
SMO的理论部分终于结束，我相信绝大多数人都不一定会喜欢这么多的数学公式的，但是没办法，想要实现算法必须得了解每一道公式得含义，既然已经学习完理论部分，我们就可以上手编程，使用的软件为MATLAB。整理编程主要就三大版块：选择变量，更新变量，更新参数
代码如下：
clc;
clear;
tic
%% 数据提取
data = [
27 53 -1
49 37 -1
56 39 -1
28 60 -1
68 75 1
57 69 1
64 62 1
77 68 1
70 54 1
56 63 1
25 41 -1
66 34 1
55 79 1
77 31 -1
46 66 1
30 23 -1
21 45 -1
68 42 -1
43 43 -1
56 59 1
79 68 1
60 34 -1
49 32 -1
80 79 1
77 46 1
26 66 1
29 29 -1
77 34 1
20 71 -1
49 25 -1
58 65 1
33 57 -1
31 79 1
20 78 1
77 37 -1
73 34 -1
60 26 -1
77 66 1
71 75 1
35 36 -1
49 61 1
26 37 -1
42 73 1
36 50 -1
66 73 1
71 43 1
33 62 1
43 41 -1
42 29 -1
58 20 -1
];
Alpha = [];
% save('data.mat') % 保存数据，不需要运行
X = data(:,1:2); % 提取特征向量
y = data(:,3); % 提取标签
% data_1 =[];data_2 = [];
% for i = 1:length(data)
% if data(i,3)==1
% data_1 = [data_1;data(i,1:2)];
% else
% data_2 = [data_2;data(i,1:2)];
% end
% end
% scatter(data_1(:,1),data_1(:,2),'*r')
% hold on
% scatter(data_2(:,1),data_2(:,2),'*b')
%% SMO序列最小优化算法
b = 0; % 初始化b
[m,~] = size(X); % 获得特征向量的大小
alphas = zeros(m,1); % 构造初始解
iter = 0; % 初始化迭代次数
err = 0.01; % 这是误差允许范围
Iter_max = 300; % 最大迭代次数
E = zeros(m,1); % 初始化存储矩阵
i_isok =zeros(m,1); % 记录无法选择的变量
j_isok = zeros(m,1); % 记录无法选择的变量
MSE = zeros(Iter_max,1);
C = 0.01; % 惩罚参数
while (iter <= Iter_max)
disp(['当前迭代的次数为：',num2str(iter)])
% 选择第一个变量
violatedMax = 0; % 定义中间变量
I = 0; % 初始化第一个变量的下标
for i =1:m % 遍历循环，选择第一个变量
G_X = (X(i,:)*X')*(alphas .* y) + b; % 先计算G_X的值
E(i) = G_X - y(i); % 计算E的值
% 判断KKT条件
% 如果违反第一个条件
if (alphas(i) == 0 && y(i) * G_X < 1 - err)
err_i = abs(1 - err - y(i) * G_X); % 计算误差的绝对值大小
% 如果误差比最大误差还大并且在i的记录中没有当前i时i，选择此i为第一个优化变量
if (err_i > violatedMax && i_isok(i) == 0)
violatedMax = err_i; % 更新最大误差
I = i; % 更新当前i为第一个变量
end
end
% 如果违反第二个条件
if (alphas(i) > 0 && alphas(i) < C && (y(i) * G_X < 1 - err) || (y(i) * G_X > 1 + err) )
err_i = max(abs(1 - err - y(i) * G_X),abs(1 + err - y(i) * G_X)); % 计算此时的误差绝对值大小
% 如果误差比最大误差还大并且在i的记录中没有当前i时i，选择此i为第一个优化变量
if (err_i > violatedMax && i_isok(i) == 0)
violatedMax = err_i; % 更新最大误差
I = i; % 更新当前i为第一个变量
end
end
% 如果违反第三个条件
if (alphas(i) == C && y(i) * G_X > 1 + err)
err_i = abs(1 + err - y(i) * G_X); % 计算误差的绝对值大小
% 如果误差比最大误差还大并且在i的记录中没有当前i时i，选择此i为第一个优化变量
if (err_i > violatedMax && i_isok(i) == 0)
violatedMax = err_i; % 更新最大误差
I = i; % 更新当前i为第一个变量
end
end
% 第一个变量选择的结果是要么I = i,要么I = -1;
end
% 找到第一个变量后，下面找第二个变量
J = 0; % 初始化第二个变量的标记
if (I~=0) % 如果第一个变量已经被选
% 第一种情况，详细参考《统计学习方法》147页
if E(I) > 0
E2_value = max(E); % 初始化E2的中间值为E的最大，后续逐渐找到最小值
for j = 1:m
% 如果当前的j满足条件，则选择当前的j
if (E(j) < E2_value && j ~= I && j_isok(j)==0)
E2_value = E(j); % 将更小的值赋给E2_value
J = j; % 并将选择次j为第二个变量
end
end
end
% 第二种情况
if E(I) < 0
E2_value = min(E); % 初始化E2_value的值
for j = 1:m
% 如果当前的j满足条件，则选择当前的j
if (E(j) > E2_value && j ~= I && j_isok(j)==0)
E2_value = E(j); % 将更大的值赋给E2_value
J = j; % 并将选择次j为第二个变量
end
end
end
% 如果第一个变量没有选，第二个变量也没有选，则说明循环结束
if I==0 && J==0
break;
end
% 在此I的情况下未找到合适的J，需要重新寻找I
if I ~= 0 && J == 0
j_isok = zeros(m,1); % 清空禁搜表
i_isok(I) = 1; % 说明这个I的情况下没有合适的j值，需要重新挑选
continue
end
% 下面计算变量的更新方法
alphaI_old = alphas(I);
alphaJ_old = alphas(J);
% 计算H和L，具体参考《统计学习方法》的145页
if (y(I) ~= y(J))
L = max(0,alphaJ_old-alphaI_old);
H = min(C,C + alphaJ_old-alphaI_old);
else
L = max(0,alphaJ_old + alphaI_old - C);
H = min(C,alphaJ_old + alphaI_old);
end
if L == H
j_isok(J) = 1; % 并将当前的第二个变量放进禁搜表中
continue
end
% 如果当前的第二个变量已经选中，则进行后续的计算
eta = X(I,:) * X(I,:)' + X(J,:) * X(J,:)' - 2 * X(I,:) * X(J,:)';
alphasJ_unCut = alphaJ_old + y(J) * (E(I) - E(J)) / eta; % 计算未剪辑的解
% 这里是通过将剪辑后的解与剪辑前的解做对比，如果变化小，则重新选择第二个变量
if (abs(clipAlpha(alphasJ_unCut,H,L) - alphaJ_old) < 0.0001)
j_isok(J) = 1;
continue;
end
% 更新I，J
alphas(J) = clipAlpha(alphasJ_unCut,H,L);
alphas(I) = alphaI_old + y(I) * y(J) * (alphaJ_old - alphas(J));
% 更新b1
b1 = b - E(I) - y(I) * X(I,:) * X(I,:)' * (alphas(I) - alphaI_old) - y(J) * X(J,:) * X(I,:)'* (alphas(J) - alphaJ_old);
% 更新b2
b2 = b - E(J) - y(I) * X(I,:) * X(J,:)' * (alphas(I) - alphaI_old) - y(J) * X(J,:) * X(J,:)' * (alphas(J) - alphaJ_old);
% 更新总的b
if (0 < alphas(I) && C > alphas(I))
b = b1;
elseif (0 < alphas(J) && C > alphas(J))
b = b2;
else
b = (b1 + b2) / 2;
end
end
% 如果所有的I都是满足KKT条件的，就结束循环
iter = iter + 1; % 迭代次数加一
i_isok = zeros(m,1); % 清空禁选表
j_isok = zeros(m,1); % 清空禁选表
% draw(alphas,X,y)
% pause(0.1)
% hold off
MSE(iter) = Distance(alphas,X,y); % 计算迭代误差
Alpha = [Alpha,alphas];
end修正函数：
function [alpha] = clipAlpha(alpha,H,L)
%CLIPALPHA
% 对alpha的值进行剪枝
if alpha > H
alpha = H;
end
if L > alpha
alpha = L;
end
end绘图函数：
为了便于以后理解，上面大部分代码都已经加好了注释，并且就是按照文章实现的思路进行的。
其中部分环节没有和SMO算法叙述的完全相同，因为在第二个变量选择时，我是直接遍历所有的样本然后找到 ∣E1−E2∣最大的值作为第二个变量的，没有考虑目标函数下降的问题，实验表明这样也是可以的。
评论（0）
您还未登录，请登录后发表或查看评论