支持向量机的smo算法（MATLAB code）

% function [alpha,bias] = smo(X, y, C, tol)

function model = smo(X, y, C, tol)

% SMO: SMO algorithm for SVM
%
%Implementation of the Sequential Minimal Optimization (SMO)
%training algorithm for Vapnik's Support Vector Machine (SVM)
%
% This is a modified code from Gavin Cawley's MATLAB Support
% Vector Machine Toolbox
% (c) September 2000.
%
% Diego Andres Alvarez.
%
% USAGE: [alpha,bias] = smo(K, y, C, tol)
%
% INPUT:
%
%   K: n x n kernel matrix
%   y: 1 x n vector of labels, -1 or 1
%   C: a regularization parameter such that 0 <= alpha_i <= C/n
%   tol: tolerance for terminating criterion
%
% OUTPUT:
%
%   alpha: 1 x n lagrange multiplier coefficient
%   bias: scalar bias (offset) term

% Input/output arguments modified by JooSeuk Kim and Clayton Scott, 2007

global SMO;

y = y';

ntp = size(X,1);
%recompute C
% C = C/ntp;

%initialize
ii0 = find(y == -1);
ii1 = find(y == 1);

i0 = ii0(1);
i1 = ii1(1);

alpha_init = zeros(ntp, 1);
alpha_init(i0) = C;
alpha_init(i1) = C;
bias_init = C*(X(i0,:)*X(i1,:)' -X(i0,:)*X(i1,:)') + 1;

%Inicializando las variables
SMO.epsilon = 10^(-6); SMO.tolerance = tol;
SMO.y = y'; SMO.C = C;
SMO.alpha = alpha_init; SMO.bias = bias_init;
SMO.ntp = ntp; %number of training points

%CACHES:
SMO.Kcache = X*X'; %kernel evaluations
SMO.error = zeros(SMO.ntp,1); %error

numChanged = 0; examineAll = 1;

%When all data were examined and no changes done the loop reachs its
%end. Otherwise, loops with all data and likely support vector are
%alternated until all support vector be found.
while ((numChanged > 0) || examineAll)
numChanged = 0;
if examineAll
%Loop sobre todos los puntos
for i = 1:ntp
numChanged = numChanged + examineExample(i);
end;
else
%Loop sobre KKT points
for i = 1:ntp
%Solo los puntos que violan las condiciones KKT
if (SMO.alpha(i)>SMO.epsilon) && (SMO.alpha(i)<(SMO.C-SMO.epsilon))
numChanged = numChanged + examineExample(i);
end;
end;
end;

if (examineAll == 1)
examineAll = 0;
elseif (numChanged == 0)
examineAll = 1;
end;
end;
alpha = SMO.alpha';
alpha(alpha < SMO.epsilon) = 0;
alpha(alpha > C-SMO.epsilon) = C;
bias = -SMO.bias;

model.w = (y.*alpha)* X; %%%%%%%%%%%%%%%%%%%%%%
model.b = bias;

return;

function RESULT = fwd(n)
global SMO;
LN = length(n);
RESULT = -SMO.bias + sum(repmat(SMO.y,1,LN) .* repmat(SMO.alpha,1,LN) .* SMO.Kcache(:,n))';
return;

function RESULT = examineExample(i2)
%First heuristic selects i2 and asks to examineExample to find a
%second point (i1) in order to do an optimization step with two
%Lagrange multipliers

global SMO;
alpha2 = SMO.alpha(i2); y2 = SMO.y(i2);

if ((alpha2 > SMO.epsilon) && (alpha2 < (SMO.C-SMO.epsilon)))
e2 = SMO.error(i2);
else
e2 = fwd(i2) - y2;
end;

% r2 < 0 if point i2 is placed between margin (-1)-(+1)
% Otherwise r2 is > 0. r2 = f2*y2-1

r2 = e2*y2;
%KKT conditions:
% r2>0 and alpha2==0 (well classified)
% r2==0 and 0% r2<0 and alpha2==C (support vectors between margins)
%
% Test the KKT conditions for the current i2 point.
%
% If a point is well classified its alpha must be 0 or if
% it is out of its margin its alpha must be C. If it is at margin
% its alpha must be between 0%take action only if i2 violates Karush-Kuhn-Tucker conditions
if ((r2 < -SMO.tolerance) && (alpha2 < (SMO.C-SMO.epsilon))) || ...
((r2 > SMO.tolerance) && (alpha2 > SMO.epsilon))
% If it doens't violate KKT conditions then exit, otherwise continue.

%Try i2 by three ways; if successful, then immediately return 1;
RESULT = 1;
% First the routine tries to find an i1 lagrange multiplier that
% maximizes the measure |E1-E2|. As large this value is as bigger
% the dual objective function becames.
% In this first test, only support vectors will be tested.

POS = find((SMO.alpha > SMO.epsilon) & (SMO.alpha < (SMO.C-SMO.epsilon)));
[MAX,i1] = max(abs(e2 - SMO.error(POS)));
if ~isempty(i1)
if takeStep(i1, i2, e2), return;
end;
end;

%The second heuristic choose any Lagrange Multiplier that is a SV and tries to optimize
for i1 = randperm(SMO.ntp)
if (SMO.alpha(i1) > SMO.epsilon) & (SMO.alpha(i1) < (SMO.C-SMO.epsilon))
%if a good i1 is found, optimise
if takeStep(i1, i2, e2), return;
end;
end
end

%if both heuristc above fail, iterate over all data set
for i1 = randperm(SMO.ntp)
if ~((SMO.alpha(i1) > SMO.epsilon) & (SMO.alpha(i1) < (SMO.C-SMO.epsilon)))
if takeStep(i1, i2, e2), return;
end;
end
end;
end;

%no progress possible
RESULT = 0;
return;

function RESULT = takeStep(i1, i2, e2)
% for a pair of alpha indexes, verify if it is possible to execute
% the optimisation described by Platt.

global SMO;
RESULT = 0;
if (i1 == i2), return;
end;

% compute upper and lower constraints, L and H, on multiplier a2
alpha1 = SMO.alpha(i1); alpha2 = SMO.alpha(i2);
y1 = SMO.y(i1); y2 = SMO.y(i2);
C = SMO.C; K = SMO.Kcache;

s = y1*y2;
if (y1 ~= y2)
L = max(0, alpha2-alpha1); H = min(C, alpha2-alpha1+C);
else
L = max(0, alpha1+alpha2-C); H = min(C, alpha1+alpha2);
end;

if (L == H), return;
end;

if (alpha1 > SMO.epsilon) & (alpha1 < (C-SMO.epsilon))
e1 = SMO.error(i1);
else
e1 = fwd(i1) - y1;
end;

%if (alpha2 > SMO.epsilon) & (alpha2 < (C-SMO.epsilon))
% e2 = SMO.error(i2);
%else
% e2 = fwd(i2) - y2;
%end;

%compute eta
k11 = K(i1,i1); k12 = K(i1,i2); k22 = K(i2,i2);
eta = 2.0*k12-k11-k22;

%recompute Lagrange multiplier for pattern i2
if (eta < 0.0)
a2 = alpha2 - y2*(e1 - e2)/eta;

%constrain a2 to lie between L and H
if (a2 < L)
a2 = L;
elseif (a2 > H)
a2 = H;
end;
else
%When eta is not negative, the objective function W should be
%evaluated at each end of the line segment. Only those terms in the
%objective function that depend on alpha2 need be evaluated...

ind = find(SMO.alpha>0);

aa2 = L; aa1 = alpha1 + s*(alpha2-aa2);

Lobj = aa1 + aa2 + sum((-y1*aa1/2).*SMO.y(ind).*K(ind,i1) + (-y2*aa2/2).*SMO.y(ind).*K(ind,i2));

aa2 = H; aa1 = alpha1 + s*(alpha2-aa2);
Hobj = aa1 + aa2 + sum((-y1*aa1/2).*SMO.y(ind).*K(ind,i1) + (-y2*aa2/2).*SMO.y(ind).*K(ind,i2));

if (Lobj>Hobj+SMO.epsilon)
a2 = H;
elseif (Lobj<Hobj-SMO.epsilon)
a2 = L;
else
a2 = alpha2;
end;
end;

if (abs(a2-alpha2) < SMO.epsilon*(a2+alpha2+SMO.epsilon))
return;
end;

% recompute Lagrange multiplier for pattern i1
a1 = alpha1 + s*(alpha2-a2);

w1 = y1*(a1 - alpha1); w2 = y2*(a2 - alpha2);

%update threshold to reflect change in Lagrange multipliers
b1 = SMO.bias + e1 + w1*k11 + w2*k12;
bold = SMO.bias;

if (a1>SMO.epsilon) & (a1<(C-SMO.epsilon))
SMO.bias = b1;
else
b2 = SMO.bias + e2 + w1*k12 + w2*k22;
if (a2>SMO.epsilon) & (a2<(C-SMO.epsilon))
SMO.bias = b2;
else
SMO.bias = (b1 + b2)/2;
end;
end;

% update error cache using new Lagrange multipliers
SMO.error = SMO.error + w1*K(:,i1) + w2*K(:,i2) + bold - SMO.bias;
SMO.error(i1) = 0.0; SMO.error(i2) = 0.0;

% update vector of Lagrange multipliers
SMO.alpha(i1) = a1; SMO.alpha(i2) = a2;

%report progress made
RESULT = 1;
return;

clear
X = []; Y=[];
figure;
% Initialize training data to empty; will get points from user
% Obtain points froom the user:
trainPoints=X;
trainLabels=Y;
clf;
axis([-5 5 -5 5]);
if isempty(trainPoints)
% Define the symbols and colors we'll use in the plots later
symbols = {'o','x'};
classvals = [-1 1];
trainLabels=[];
hold on; % Allow for overwriting existing plots
xlim([-5 5]); ylim([-5 5]);

for c = 1:2
title(sprintf('Click to create points from class %d. Press enter when finished.', c));
[x, y] = getpts;

plot(x,y,symbols{c},'LineWidth', 2, 'Color', 'black');

% Grow the data and label matrices
trainPoints = vertcat(trainPoints, [x y]);
trainLabels = vertcat(trainLabels, repmat(classvals(c), numel(x), 1));
end

end

% C = 10;tol = 0.001;
% par = SMOforSVM(trainPoints, trainLabels , C, tol );
% p=length(par.b); m=size(trainPoints,2);
%  if m==2
% %     for i=1:p
% %         plot(X(lc(i)-l(i)+1:lc(i),1),X(lc(i)-l(i)+1:lc(i),2),'bo')
% %         hold on
% %     end
%     k = -par.w(1)/par.w(2);
%     b0 = - par.b/par.w(2);
%     bdown=(-par.b-1)/par.w(2);
%     bup=(-par.b+1)/par.w(2);
%     for i=1:p
%         hold on
%         h = refline(k,b0(i));
%         set(h, 'Color', 'r')
%         hdown=refline(k,bdown(i));
%         set(hdown, 'Color', 'b')
%         hup=refline(k,bup(i));
%         set(hup, 'Color', 'b')
%     end
%  end
% xlim([-5 5]); ylim([-5 5]);
%
% pause

C = 10;tol = 0.001;
par = smo(trainPoints, trainLabels, C, tol);
p=length(par.b); m=size(trainPoints,2);
if m==2
%     for i=1:p
%         plot(X(lc(i)-l(i)+1:lc(i),1),X(lc(i)-l(i)+1:lc(i),2),'bo')
%         hold on
%     end
k = -par.w(1)/par.w(2);
b0 = - par.b/par.w(2);
bdown=(-par.b-1)/par.w(2);
bup=(-par.b+1)/par.w(2);
for i=1:p
hold on
h = refline(k,b0(i));
set(h, 'Color', 'r')
hdown=refline(k,bdown(i));
set(hdown, 'Color', 'b')
hup=refline(k,bup(i));
set(hup, 'Color', 'b')
end
end
xlim([-5 5]); ylim([-5 5]);