svmnlspex01.m

数据挖掘的新方法－支持向量机书中算法例子

%SVMnLSPex01.m
%
%Two Dimension SVM Problem, Two Class  Un-separable Situation
%
%Method from Thorsten Joachims:
%"Making Large-Scale SVM Learning Practical", function (4)/(5)/(6)
%	
%  Objective:   min "f(A)=-sum(A)+(A'*Q*A)/2" ,							function (4);
%					 here:	Q = diag(Y)*K(X,X')*diag(Y)  and in matrix Q any component Q(ij)=yi*yj*xi*xj';
%					 K(x,y) is Kernel Function, it define a opration relation between two vectors x and y;
%  Subject to:  sum{A'*Y}=0 ,													function (5);
%					 0 <= ai <= C ,												function (6);
%					 here yi={1,-1} indicates a sample xi belongs to either class of two ;
%The optimizing variables is "Lagrange Multipliers": A=[a1,a2,...,am],m is the number of total samples.
%
%For Support Vector Machines, when two class is linear separable, it is originally a Quadratic Program problem:
%
%	Objective:		Min "||W||"
%	Subject to:		yi(W*xi + b) - 1 >= 0, i = 1,2,...,m
%
% Here, X = [xi] is Learning Samples and yi = {-1, 1} correspondes each of xi to indicate which class it is.
%
%

% The Linear-Separable Situation:
%
% When the optimizing problem is replaced by function (4) and (5)/(6), the variables is altered as A = [ai].
% From the QP problem (4),(5) and (6) we can get all of the A = [ai],and the original object function,can be  
% gotten indirectly from  W = Sum(ai*yi*xi).
% In the other hand, ai = 0 means that the sample xi is a SV (Support Vector). And then, we can get each SV by
% check ai,in succession, substitute arbitrary SV (xi) we have gotten to function "yi(W*xi + b) - 1 >= 0",to get
% the intercept b. 
% Thus the classification function f(x) = W*x + b has been obtained. Gaving an arbitrary sample xs, if f(xs)>0 xs
% is decided to class A otherwise if f(xs)<0 xs should be class B. Take f(x) = 0, we can get a classifying line.

% The Nonlinear Separable Situation:
%
% We call the space where original problem lies "L" space and assume there is a mapping function F map "L" to a high
% space "H". Then each sample/vector xi in L corresponds a vector Xi in H and Xi = F(xi). 
% Assuming two linear un-separable set AL and BL in L, their mapped sets AH and BH in H can be linear separated.
% So, if we know the mapping function F we can still use the former method to get the classification function f(X)in H.
% But unfortunately, finding mapping F is difficult. In the Other hand, in high space H the opration about samples is
% only Vector Dot Multiplication F(x)*F(y) (in function 4), so we can define an operation K(x,y) between two L space
% samples x and y equivalent to a Dot Multiplying in H space, that is F(x)*F(y)= K(x,y). 
%
% By this, the object function in H space is:
%  Objective:  min "-sum(A)+(A'*Q*A)/2" ,							function (4);
%		  here:  Q = diag(Y)*K(X,X')*diag(Y);
% K(x,y) is so-call Kernel Function. Now, the A = [ai] is still Lagrange Multipliers but its counterpart sample is
% a vector F(xi) in H and from each ai = 0 we can find the H space SV F(xi,ai = 0). 
%
% Nevertheless, by hypothesis for any sample x0
% the Classifying Function in H is line:  				f(x0) = W*F(x0) + b;								(a)
% and the constraint by any sample x0 is:					yi*(W*F(x0) + b) - 1 >=0;						(b)
% W in H space should be:										W = sum(ai*yi*F(xi));							(c)								
% However the Dot Multiplying W*F(x0) can wrote as: 	W*F(x0) = sum(ai*yi*K(xi,x0));				(d)
%                                     
% Finally, we can check A = [ai] zero to get SV F(xj), substitute xj into (b) to obtaine b, use b and expression (d) 
% complete Classifying Function (a). Giving a sample x0, it should A when f(x0) > 0 otherwise class B when f(x0) < 0.
% Giving a continued series of x0 we can draw a classifying line.
%
% In this example, we use polynomial Kernel Function.
%

function mm=SVMnLSex01(mm) %mm -- select sample number

echo off;
close all;
fclose all;

%Generate Two Set of Data each set with m samples
	m1=mm;					% The number of class A which yi = 1 ;
   m2=mm;					% The number of class B which yi = -1;
   n = 2;					% The dimension of sample vector, that is xi = [xi(1), xi(2)];
   al=pi/8;					%
   
   x=rand(m1,2);
	X1(1:m1,1)=2*x(1:m1,1);
	X1(1:m1,2)=x(1:m1,2)+2.5+sin(x(1:m1,1)*pi*2.7);
	%X1(1:m1,1)=x(1:m1,1)*cos(al)-x(1:m1,2)*sin(al);
	%X1(1:m1,2)=x(1:m1,2)*cos(al)+x(1:m1,1)*sin(al);

	x=rand(m2,2);
	X2(1:m2,1)=2*x(1:m2,1);
	X2(1:m2,2)=x(1:m2,2)+1+sin(x(1:m2,1)*pi*2.7);
	%X2(1:m2,1)=x(1:m2,1)*cos(al)-x(1:m2,2)*sin(al);
	%X2(1:m2,2)=x(1:m2,2)*cos(al)+x(1:m2,1)*sin(al);
   
   clear x;
   
   X=[X1;X2];							%X1 is Positive Samples Set and X2 is Negative Samples Set;
   m=m1+m2;								%Total number of samples in matrix X(m*n);
     
   figure(1);
  	whitebg(1,[0,0.1,0]);
   plot(X(1:m1,1),X(1:m1,2),'m+');
   hold on
   plot(X(m1+1:m,1),X(m1+1:m,2),'c*');
   
   Xa = min(0,min(X(1:m,1)));		%Minimum of X axes of X
   Xb = max(X(1:m,1));				%Maximum of X axes of X
   Ya = min(0,min(X(1:m,2)));		%Minimum of Y axes of Y
   Yb = max(X(1:m,2));				%Maximum of Y axes of Y

	axis([Xa,Xb,min(X(1:m,2)),max(X(1:m,2))])
   
   asw = 0;
   while (asw~=1)&(asw~=2)
      asw=input('Continue or no? (1/2)');
   end
   if asw==1
   	%
		%
		%Transmiting Samples to Object and Constraint Function, in SVMnLSex01FUN.m,
		fid=fopen('a.dat','w');
      fwrite(fid,[m1,m2,n],'float');
      for k=1:n
         fwrite(fid,X(1:m,k),'float');
      end
      fclose(fid);
		%
		%Keep to MATLAB prescribe,Here
		%x0 is start point to optimizing;
      %Set options(13)=1 indicate the first constraint function contained in SVMnLSPex01FUN.m is a equality,
      %and others are "<=0",Reference to MATLAB function "CONSTR". Of cause, here we have only one constraint;
      %
   	x0 = 0.5*ones(m,1);	%Initial point of optimizing					
      options(13) = 1;
      
	   C = 500;					
      VLB = zeros(m,1);			%	The lower-bound of Lagrange Multiplier ai
      VUB = C*ones(m,1);		%	The up-bound of Lagrange Multiplier ai
%
% About The up-bound of Lagrange Multiplier:  
% A Multiplier a(i) isn't zero indicate a corresponding constranit in the original optimizing problem is active.
% That is the constraint yi*(W*xi +b )-1>=0 is active, on the other words, sample (xi) is a SVs.
% The other Multiplier aj = 0 corresponding to the samples who aren't SV.
%
% But, too small up-bound imposed to Multipliers will cause a Mutiplier who corrrespond un-active constraint 
% cann't descent to zeros in optimizing procedure. 
% However, we decide if a sample is SV just only dependent on the fact that its Multiplier is zero or no,so
% this will cause difficult when distinguish SVs from samples.
%
     
[A,options,L] = constr('SVMnLSPex01FUN', x0, options, VLB, VUB);

      %
		%Return parameter: 
      % A=[ai] optimizing variables, Lagrange Multipliyer to original optimizing problem correspond to each samples;
      % Each component of A, ai, corresponds to a sample xi which formed a constraint in original optimizing problem:
      % Min ||W||, Subject to yi(W*xi + b) - 1>=0, i = 1,2,...,m;
      % 
      %
      %Un-separable Situation
      
	   Y=[ones(m1,1);-ones(m2,1)];	%The classfication value for both Positive and Negative Samples Set
		W = sum((diag(A)*diag(Y))*X);	%Normal vector of the superplane which divides two classes with largest margin.

		[z,I] = max(A(1:m));				%The Support Vectors xj are corresponding to a(j) <> 0;
      XI = X(I,1:n)';					%Pick out a arbitrary Support Vector to determine the "b" in classifying function
      YI = Y(I);
      
      %
      %Assume have gotten a SV x0, then
      %y0*(W*x0 +b )-1 = y0*(sum(ai*yi*K(xi,x0)) + b)-1 >= 0,  b = y0-sum(ai*yi*K(xi,x0)) = y0-(diag(A)*Y)'*K(X,x0);
      %
      
      b = YI - (diag(A)*Y)'*K(X,XI);
      
      Xd = [Xa:0.05:Xb+0.025];		%X axes partition 
      Yk = [Ya:0.005:Yb];				%For a fine curve you can set small interval but will need more time
      
      for k = 1:length(Xd),							%Giving a X axes point
         Xy = [Xd(k)*ones(1,length(Yk));Yk];		%Construct matrix depict a comlun line cross X 
         fx = (diag(A)*Y)'*K(X,Xy) + b;			
         [l,p] = min(abs(fx));
         Yd(k) = Yk(p);
         if max(fx)<=0,
            	Yd(k) = max(Yk);
         	else if min(fx)>=0,
                  Yd(k) = min(Yk);
               end
            end
      end
      
      plot(Xd,Yd,'y')
         
   	for n=1:m										% ai = 0 corresponds a SV
 	     if A(n) > 1e-6,
            plot(X(n,1),X(n,2),'wo');
 	     end
	   end
      
      delete('a.dat');

end


function Kernel = K(U,E)									%Kernel function
      Kernel = (U*E+1).^4;