%%%%%%%%%%%% Gaussian SVM %%%%%%%%%%%%%%%%%%%%
clear
close all
%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%
NNN=50; %%% Training Sample number
BETA=5; %%%%% 1 5 25 125 Parameter of Gaussian Kernel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NOISEY=0.0; %%% NOISE LEVEL in Training Samples
CCC=1000; %%% PARAMETR for soft margin
CYCLE=20000; %%% Cycle of optimization process
DISPLAY=20; %%% Diplay Number
INTERVAL=(CYCLE/20); %%% Drawing cycle
ETA=0.01; %%% Optimization coefficient
LLL=1; %%% Parameter for Condition dot(YYY,Alpha)=0
SUPPORT=0.01; %%% Judgement of support vector
RAND_0=100; %%% random seed
Complexity=3*pi; %%% compexity of the true distribution
%%%%%%%%%%%%%%%%%%%%%%%%%% Ture parameter %%%%%%%%%%%%%%%%%%%%
A10=1;
A20=1;
B00=-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Make Inputs and Outputs %%%%%%%%%%%%%%%%%%%%%%%%%%%
XXX=zeros(2,NNN); %%% (2,NNN) Inputs of samples
XXX(1,:)=rand(1,NNN);
XXX(2,:)=rand(1,NNN);
HHH=0.5*A10*(ones(1,NNN)+sin(Complexity*XXX(1,:)))+A20*XXX(2,:)+B00+NOISEY*(1-2*rand(1,NNN));
YYY=sign(HHH); %%% (1,NNN) Outputs of samples
ZZZ=(YYY+1)/2;
for i=1:1:NNN
subplot(2,2,1); plot(XXX(1,i),XXX(2,i),'+','Color',[ZZZ(i) 0 1-ZZZ(i)]); hold on
end
title('Blue -1 : Red : +1');
drawnow;
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Kernel%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% KXX=XXX'*XXX; %%% (NNN,NNN) (x,x)
KXX=zeros(NNN,NNN);
for i=1:1:NNN
for j=1:1:NNN
DXX=XXX(:,i)-XXX(:,j);
KXX(i,j)=DXX'*DXX;
end
end
KXX=exp(-BETA*KXX); %%% (NNN,NNN)
KYY=YYY'*YYY; %%% (NNN,NNN) yy
KYYXX=KYY.*KXX; %%% (NNN,NNN) yy(x,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Initialize Optimization%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CCCS=CCC*ones(1,NNN);
MATYX=KYYXX+eye(NNN)+LLL*YYY'*YYY;
AT=MATYX\ones(NNN,1);
Alpha=AT';
AlphaA=zeros(1,NNN);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Dual Optimization %%%%%%%%%%%%%%%%%%%%%%%%
CCCS=CCC*ones(1,NNN);
MATYX=KYYXX+eye(NNN)+LLL*YYY'*YYY;
AT=MATYX\ones(NNN,1);
Alpha=AT';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Optimization%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Dual Parameters are optimized by steepest descent %%%%%%%%%%%%%%
III=zeros(1,DISPLAY);
EEE=zeros(1,DISPLAY);
COND=zeros(1,DISPLAY);
for i=1:1:CYCLE
ETA2=ETA*(1-i/CYCLE);
TTT=dot(YYY,Alpha);
AB=Alpha+ETA2*(1-Alpha*KYYXX'-LLL*TTT*YYY); %%% Steepest Descent
AB=AB-dot(YYY,AB)*YYY/NNN; %%%
AB=(AB+abs(AB))/2; %%% make ( AB >= 0 )
AB=(AB+CCCS-abs(AB-CCCS))/2; %%% make (AB <= CCCS )
Alpha=AB; %%% Update Alpha
if(mod(i,INTERVAL)==0)
III(i/INTERVAL)=i;
EEE(i/INTERVAL)=dot(Alpha,ones(1,NNN))-0.5*Alpha*KYYXX*Alpha';
COND(i/INTERVAL)=0.5*LLL*TTT*TTT;
end
end
subplot(2,2,2);
plot(III,EEE,'r-');
title('Horizontal: cycle, Red: Dual Loss');
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% WWW and BBB are calculated %%%%%%%%%%%%%%%%%%%%
fprintf('Optimization Completed \n');
TTT=dot(YYY,Alpha);
CONDITION=TTT*TTT;
fprintf('CONDITION=%e should be sufficiently small. \n',CONDITION); %%%%% dot(YYY,Alpha)=0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Result Drawing %%%%%%
COUNTER=0;
BBB=0;
for i=1:1:NNN
if(Alpha(i)>SUPPORT)
%%% Alpha(i) %%% SUPPORT VECTOR
subplot(2,2,3); plot(XXX(1,i),XXX(2,i),'o','Color',[0 1 0]); hold on
UUU=0;
for j=1:1:NNN
if(Alpha(j)>SUPPORT)
UUU=UUU+Alpha(j)*YYY(j)*KXX(j,i);
end
end
BBB=BBB+YYY(1,i)-UUU;
COUNTER=COUNTER+1;
end
end
BBB=BBB/COUNTER;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Display Results %%%%%%%%%%%%%%
for i=1:1:NNN
subplot(2,2,3);
plot(XXX(1,i),XXX(2,i),'+','Color',[ZZZ(i) 0 1-ZZZ(i)]); hold on
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[X10, X20]=ndgrid(0:0.05:1,0:0.05:1);
Z00=0.5*A10*(ones(21,21)+sin(Complexity*X10))+A20*X20+B00;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ZZZ=zeros(21,21);
for p=1:1:21
for q=1:1:21
for i=1:1:NNN
if(Alpha(i)>SUPPORT)
%%%if(i==10)
D10=XXX(1,i)-X10(p,q);
D20=XXX(2,i)-X20(p,q);
ZZZ(p,q)=ZZZ(p,q)+Alpha(i)*YYY(i)*exp(-BETA*(D10*D10+D20*D20));
end
end
ZZZ(p,q)=ZZZ(p,q)+BBB;
end
end
subplot(2,2,3);
contour(X10,X20,Z00,[0 0],'r'); hold on
contour(X10,X20,ZZZ,[0 0],'k'); hold on
title('Red: True, Black: Estimated, o: Support Vectors');
hold off
subplot(2,2,4); bar(Alpha);
title('Optimized Values of Dual Variables');
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Generalization Error
Correctans=0;
for p=1:1:21
for q=1:1:21
H00=0.5*A10*(1+sin(Complexity*X10(p,q)))+A20*X20(p,q)+B00;
Y00=sign(H00); %%%
if(Y00*ZZZ(p,q)>0)
Correctans = Correctans + 1;
end
end
end
fprintf('Number of Support Vectors = %g. \n',COUNTER);
fprintf('Generalization: Recognition Rate = %f. \n',Correctans/441);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%