osunaint.m

支持向量机 这是我的老师编写的matlab源文件 希望对大家有用

function [nsv, alpha, b0] = osuna(X,Y,trainsize,transize,ker,C)
%the Improved Osuna SVC Support Vector Classification
%
%  [nsv, alpha, b0] = osuna(X,Y,trainsize,transize,ker,C)
%
%  Parameters: X         - Training inputs
%              Y         - Training targets
%              trainsize - the size of one training
%              transize  - the number of the exchange examples
%              ker       - kernel function
%              C         - upper bound (non-separable case)
%              nsv       - number of support vectors
%              alpha     - Lagrange Multipliers
%              b0        - bias term
% 
%  Author: Zhou weida ([e-mail])

  if (nargin <4 | nargin>6) % check correct number of arguments
    help osuna
 else
    st = cputime;

    fprintf('Osnua Support Vector Classification\n')
    fprintf('_____________________________\n')
    
        if (nargin<6) C=Inf;, end
    if (nargin<5) ker='linear';, end

    % tolerance for Support Vector Detection
    epsilon = svtol(C);
    kktepsilon=epsilon;
    %construct the samples index randomly for the first training
    numsam=size(Y,1);
    xindex=uint16([round((numsam/2-1)*rand(trainsize/2,1)+0.5);...
       round(numsam/2)+round((numsam/2-1)*rand(trainsize/2,1)+0.5)]);
    xindex=uint16(sort(double(xindex)));
 
    %avoid the same train samples
    for k=1:trainsize-1
       if xindex(k,1)>=xindex(k+1,1)
          xindex(k+1,1)=uint16(double(xindex(k,1))+1);
       end
    end
    
    %initialize the lagrangre multiplier for the whole samples
    globalalpha=zeros(numsam,1);
    
    count=0;%count the number of loop
    ok=0;   %the flag to stop the loop 
    disp('number of samples against kkt:');

    
    while(ok==0)
       count=count+1;
       ok=1;
       
    %construct the (trainsize)sub-samples from X Y for training
       trainX=X(xindex,:);
       trainY=Y(xindex,:); 
       trainflag=uint8(zeros(numsam,1));      %the flag for the samples be selected for training 
       trainflag(xindex,1)=uint8(1);
       remaindindex=find(trainflag==0);%index for the unselected samples
       remaindindex=uint16(remaindindex);
    
      
       % Construct the Kernel matrix
       %fprintf('Constructing ...\n');
       H = zeros(trainsize,trainsize);  
       for i=1:trainsize
          for j=1:trainsize
              H(i,j) = double(trainY(i))*double(trainY(j))*svkernel(ker,double(trainX(i,:)),double(trainX(j,:)));
          end
       end
       % Add small amount of zero order regularisation to 
       % avoid problems when Hessian is badly conditioned.     
       %H0=H;
       H = H+1e-10*eye(size(H));

    
       %Construct Qbn;
       %add the information into the training object function(linear item)
       numn=size(remaindindex,1);
       Qbn=zeros(trainsize,1);
       for k=1:trainsize
          for l=1:numn
             if (globalalpha(remaindindex(l,1)) > epsilon)
                    Qbn(k,1)=Qbn(k,1)+globalalpha(remaindindex(l,1))*double(Y(remaindindex(l,1)))*...
                        svkernel(ker,double(X(xindex(k,1),:)),double(X(remaindindex(l,1),:)));
              end
           end
           Qbn(k,1)=Qbn(k,1)*double(Y(xindex(k,1),1));
       end
       c = -ones(trainsize,1)+Qbn;  

    
       % Set up the parameters for the Optimisation problem

       vlb = zeros(trainsize,1);      % Set the bounds: alphas >= 0
       vub = C*ones(trainsize,1);     %                 alphas <= C
       x0 = zeros(trainsize,1);       % The starting point is [0 0 0   0]
       neqcstr =1;     % Set the number of equality constraints (1 or 0)  
       A = double(Y(xindex,1)');, b = -1*(globalalpha(remaindindex,1)'*...
                  double(Y(remaindindex,1)));     % Set the constraint Ax = b
      
       % Solve the Optimisation Problem
    
                  
        [alpha lambda how] = qp(H, c, A, b, vlb, vub, x0, neqcstr);
       globalalpha(xindex,1)=alpha;

           
       % Compute the number of Support Vectors
       svi = find( globalalpha > epsilon);
       svi=uint16(svi);
       nsv = length(svi);
       
       % Implicit bias, b0
       b0 = 0;

      Hb=[];
      svii = find(globalalpha > epsilon & globalalpha < (C - epsilon));
      svii=uint16(svii);
         for i=1:length(svii)
            for j=1:length(svi)
               Hb(i,j) =double(Y(svii(i,1),1))*double(Y(svi(j,1),1))...
                  *svkernel(ker,double(X(svii(i,1),:)),double(X(svi(j,1),:)));
            end
         end
         if length(svii) > 0
            b0 =  (1/length(svii))*sum(double(Y(svii,1)) - Hb*globalalpha(svi).*double(Y(svii,1)));
         else 
            fprintf('No support vectors on margin - cannot compute bias.\n');
         end
  
       remaindX=X(remaindindex,:);
       
       %compute the output for the unselected samples
       galphaindex=find(globalalpha>epsilon);
       galphaindex=uint16(galphaindex);
       fremaindX=soutput(double(X(galphaindex,:)),double(Y(galphaindex,1)),globalalpha(galphaindex,1),ker,b0,double(remaindX));
       
       %the flag for the unselcted samples
       %rmdflag(:,1)   -  index.
       %rmdflag(:,2)   -  the kind of against KKT conditions 
       %               -  0:no against KKT
       %               -  1:alpha=0,against KKT,
       %               -  2:0<alpha<C,against KKT,
       %               -  3:alpha=C,against KKT.
       %rmdflag(:,3)   -  the measure for against KKT conditions =1-f(xi)*yi 
       rmdflag=zeros(numn,3);
       rmdflag(:,1)=double(remaindindex);
       for k=1:numn
          if globalalpha(rmdflag(k,1),1)<=kktepsilon
             if double(Y(rmdflag(k,1),1))*fremaindX(k,1)<1-kktepsilon
                ok=0;
                rmdflag(k,2)=1;
                rmdflag(k,3)=1-double(Y(rmdflag(k,1),1))*fremaindX(k,1);
             end
          end
          if (globalalpha(rmdflag(k,1),1)>=kktepsilon &...
             globalalpha(rmdflag(k,1),1)<C-kktepsilon)
          if (double(Y(rmdflag(k,1),1))*fremaindX(k,1))<(1-kktepsilon) ...
                | (double(Y(rmdflag(k,1),1))*fremaindX(k,1))>(1+kktepsilon)
                ok=0;
                rmdflag(k,2)=2;
                rmdflag(k,3)=abs(1-double(Y(rmdflag(k,1),1))*fremaindX(k,1));
             end
          end
          if globalalpha(rmdflag(k,1),1)>=C-kktepsilon
             if double(Y(rmdflag(k,1),1))*fremaindX(k,1)>1
                ok=0;
                rmdflag(k,2)=3;
                rmdflag(k,3)=double(Y(rmdflag(k,1),1))*fremaindX(k,1)-1;
             end
          end
       end
       
       % construct the exchange samples.
       
       alpha1=alpha;%save the alpha for training samples
       zeroalpha=find(alpha<kktepsilon);
       zeroalpha=uint16(zeroalpha);
       agtkkt=find(rmdflag(:,2)>0);
       agtkkt=uint16(agtkkt);
       countzero=size(zeroalpha,1);
       countagtkkt=size(agtkkt,1);
       disp(countagtkkt);
       if countzero<countagtkkt
          if countzero<=transize
             tsize=transize;
          else
             tsize=countzero-1;
          end
       else
          tsize=countagtkkt;
       end
       tranindex=uint16(zeros(tsize,2));   %index of the exchange samples 
                                   %tranindex(:,1) -  the index of out samples,
                                   %tranindex(:,2) -  the index of in samples.
       for k=1:tsize
           [talpha1,tindex1]=min(alpha);     %select the out samples whose alpha is minimun.
           alpha(tindex1,1)=C;  
           tranindex(k,1)=xindex(tindex1,1);
           [talpha2,tindex2]=max(rmdflag(:,3));%select the in samples whose the againtKKT is maxim.
           rmdflag(tindex2,3)=0;
           tranindex(k,2)=rmdflag(tindex2,1);
           xindex(tindex1,1)=uint16(rmdflag(tindex2,1));%implement the exchange .
       end
        %disp(count);
        %disp(xindex);
       xindex=uint16(sort(double(xindex)));
    end  
    
    
    
    
    %w2 = alpha'*H*alpha;
    %fprintf('|w0|^2    : %f\n',w2);
    %fprintf('Margin    : %f\n',2/sqrt(w2));
    fprintf('Sum alpha : %f\n',sum(alpha));
    fprintf('Execution time: %4.1f seconds\n',cputime - st);

   
     svi = find( globalalpha > epsilon);
     nsv = length(svi);
     alpha=globalalpha; 
     %figure(2);
     %clf;
     %svcplot(X,Y,ker,globalalpha,b0);

       
end