irwls.m

支持向量机的一个训练软件

function [nsv,al3,bi,T,Lp]=irwls(x,y,ker,C,par,tol);
%
%
%  This function solves the Support Vector Machine for pattern recognition.
%
%     [nsv,al3,bi,T,Lp]=irwls(x,y,ker,C,par,tol);
%
% This function solves the SVM using the IRWLS procedure. 
% Parameters:
%
%              x: matrix of input vectors. Each row of x1 represents a d-dimensional vector. 
%                 The number of rows is the number of training samples.
%              y: column vector of lables (-1 or +1). Each row represent the lable of the sample in x.
%              ker: the kernel to be employed. It has to be a string out of: 'linear', 'poly_h', 'poly_i' and 'rbf'.
%                Type help kernel for further information.
%              C: is the penalty factor of the SVM soft-margin formulation.
%              par: is the parameter of the used kernel. type help kernel for further information.
%              tol: tolerance for the value of Lp, the primal functional, to stop the algorithm.
%                The default value if 10^-5. (optional)
%
% Outputs:
%
%              nsv: number of support vector.
%              al3: the value of the alphas times the lables. It is the value of alpha for the positive samples.
%                and its negative value for the negatives samples.
%              bi: the value of the bias.
%              T: running time.
%              Lp: A vector that show the value of the primal \|w\|^2+C\sum_{i=0}^n \xi_i in each iteration
%               of the algorithm.
%
% Examples:
%
%       with linear kernel:
%               x=[randn(1,50)-3 randn(1,50)+3;randn(1,100)]';
%               y=[-ones(1,50) ones(1,50)]';
%               [nsv,alpha,bias]=irwls(x,y,'linear',10);
%               svcplot(x,y,'linear',alpha,bias);
%
%               W=x'*alpha;
%               X_test=[randn(1,500)-3 randn(1,500)+3;randn(1,1000)]';
%               Y_test=[-ones(1,500) ones(1,500)]';
%               Output=X_test*W+bias;
%               number_error=sum(abs(Y_test-sign(Output)))/2;
%
%               
%       with nonlinear (RBF) kernel:
%               x=[randn(1,50)-1 randn(1,50)+1;randn(1,100)]';
%               y=[-ones(1,50) ones(1,50)]';
%               [nsv,alpha,bias]=irwls(x,y,'rbf',10,4);
%               svcplot(x,y,'rbf',alpha,bias,4);
%
%               X_test=[randn(1,500)-1 randn(1,500)+1;randn(1,1000)]';
%               Y_test=[-ones(1,500) ones(1,500)]';
%               Output=kernel('rbf',X_test,x(find(alpha),:),4)*alpha(find(alpha))+bias;
%               number_error=sum(abs(Y_test-sign(Output)))/2/;
%
% Author:   Fernando Perez-Cruz ([e-mail])
% Version:  1.0
% Date:     18th March 2002.
%
%
T=clock;
if (nargin<4 | nargin>6)
   help irwls
else
   if(nargin==4 & ker(1)=='p')
       disp('The used kernel needs an input parameter. We have set it to 2.');
       par=2;
   elseif(nargin==4 & ker(1)=='r')
       fprintf(1,'The used kernel needs an input parameter. We have set it to %1.3f\n',sqrt(size(x,2)));
       par=sqrt(size(x,2));
   elseif(nargin==4)
       par=0;
   end
   if(nargin<6)
      tol=10^-5;
   end         
   K=10^6;
   N=size(x,1);
   ns=-1;
   k=2;
   hacer=1;
   bi=0;
   
   H=kernel(ker,x,x,par);
   
   i1p=1:2:N;
   i1n=2:2:N;
   i1=[i1p i1n]';   
   
   a=zeros(N,1);
   a(i1)=C;
   
   al3=zeros(N,1);
   bi=0;
   i2p=[];
   i2n=[];
   
   Lp(1)=N*C;
   e=y;
   while(hacer)
      al3_a=al3;
      bi_a=bi;
      al3=zeros(N,1);
      al3(i2p)=C;
      al3(i2n)=-C;
      if(length(i1))
         Xi=inv([H(i1,i1)+diag(1./(a(i1))) ones(length(i1),1);ones(1,length(i1)) 0]);
         aux=Xi*[y(i1)+C*H(i1,i2n)*ones(length(i2n),1)-C*H(i1,i2p)*ones(length(i2p),1);C*length(i2n)-C*length(i2p)];
         al3(i1)=aux(1:length(i1));
         bi=aux(length(i1)+1);
      end
      e_a=e;
      
      I=find(al3-al3_a);      
      if(length(I))
         e=e_a-H(:,I)*(al3(I)-al3_a(I))-(bi-bi_a);
      end      
      
      Lp(k)=al3'*H*al3/2+C*sum(e(find(e.*y>0)).*y(find(e.*y>0)));
 
      if((abs(Lp(k)-Lp(k-1))/Lp(k-1))<tol & k>(ns+1))
         nsv=length(i1)+length(i2n)+length(i2p);
         hacer=0;                  
      else                     
          if((Lp(k)-Lp(k-1))/Lp(k-1)>tol)
              LL=Lp(k-1);
              HH=Lp(k);
              posL=0;
              posH=1;
              while(posL==0)
                  jj=(posH+posL)/2;
                  al3b=al3*jj+al3_a*(1-jj);          
                  bib=bi*jj+bi_a*(1-jj);
                  I=find(al3b-al3_a);
                  if(length(I))
                      eb=e_a-H(:,I)*(al3b(I)-al3_a(I))-(bib-bi_a);              
                  end
                  if(HH>LL)
                      HHa=HH;
                      posHa=posH;
                      HH=al3b'*H*al3b/2+C*sum(eb(find(eb.*y>0)).*y(find(eb.*y>0)));
                      posH=jj;
                  else
                      LL=al3b'*H*al3b/2+C*sum(eb(find(eb.*y>0)).*y(find(eb.*y>0)));
                      posL=jj;
                  end
              end
              jj=posH+posL;
              al3b=al3*jj+al3_a*(1-jj);          
              bib=bi*jj+bi_a*(1-jj);
              I=find(al3b-al3_a);
              if(length(I))
                  eb=e_a-H(:,I)*(al3b(I)-al3_a(I))-(bib-bi_a);              
              end
              valor=al3b'*H*al3b/2+C*sum(eb(find(eb.*y>0)).*y(find(eb.*y>0)));
              
              
              posL=posH;
              LL=HH;
              posH=jj;
              HH=valor;
              if(valor>Lp(k-1))
                  while(abs(valor-Lp(k-1))/Lp(k-1)>10^-3)                      
                      jj=(posH+posL)/2;
                      al3b=al3*jj+al3_a*(1-jj);          
                      bib=bi*jj+bi_a*(1-jj);
                      I=find(al3b-al3_a);
                      if(length(I))
                          eb=e_a-H(:,I)*(al3b(I)-al3_a(I))-(bib-bi_a); 
                      end
                      valor=al3b'*H*al3b/2+C*sum(eb(find(eb.*y>0)).*y(find(eb.*y>0)));                                            
                      if(valor>Lp(k-1))                    
                          HH=valor;
                          posH=jj;
                      else
                          LL=valor;
                          posL=jj;
                      end
                  end
              end                                   
              pos=jj;
              
              al3=al3*pos+al3_a*(1-pos);
              bi=bi*pos+bi_a*(1-pos);
              I=find(al3-al3_a);      
              if(length(I))
                  e=e_a-H(:,I)*(al3(I)-al3_a(I))-(bi-bi_a);
              end  
              Lp(k)=al3'*H*al3/2+C*sum(e(find(e.*y>0)).*y(find(e.*y>0)));
 
              if(pos<10^-10)
                  hacer=0;
                  nsv=length(i1)+length(i2n)+length(i2p);
              end
          end
                     
          i1p=find(e>=0 & y==1);
          i1n=find(e<=0 & y==-1);
          
         a=zeros(N,1);
         a(i1p)=C./e(i1p);
         a(i1n)=-C./e(i1n);
         a(find(a>K))=K;

         i2p=find(al3>=(C-C/100) & al3<=(C+C/100) & al3_a>=(C-C/100) & al3_a<=(C+C/100) & e>0 & y==1 & 15*e>14*e_a);
         if(length(setdiff(i1p,i2p)))
             i1p=setdiff(i1p',i2p)';
         else
             [val,pos]=min(e(i2p));
             i2p=setdiff(i2p,i2p(pos));
             i1p=setdiff(i1p',i2p)';
         end
         
         i2n=find(al3>=(-C-C/100) & al3<=(-C+C/100) & al3_a>=(-C-C/100) & al3_a<=(-C+C/100) & e<0 & y==-1 & 15*e<14*e_a);
         if(length(setdiff(i1n,i2n)))
             i1n=setdiff(i1n',i2n)';
         else
             [val,pos]=max(e(i2n));
             i2n=setdiff(i2n,i2n(pos));
             i1n=setdiff(i1n',i2n)';
         end
         
         i1=[i1p;i1n];
         
         k=k+1;
      end   
   end
end
T=etime(clock,T);