osuna.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 <2 | nargin>6) % check correct number of arguments
    help osuna
 else
    fprintf('Osnua Support Vector Classification\n')
    fprintf('_____________________________\n')
        if (nargin<4) C=Inf;, end
    if (nargin<3) 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=[round((numsam/2-1)*rand(trainsize/2,1)+0.5);...
       round(numsam/2)+round((numsam/2-1)*rand(trainsize/2,1)+0.5)];
    xindex=sort(xindex);
 
    %avoid the same train samples
    for k=1:trainsize-1
       if xindex(k,1)>=xindex(k+1,1)
          xindex(k+1,1)=xindex(k,1)+1;
       end
    end
    
    %initialize the lagrangre multiplier for the whole samples
    globalalpha=zeros(numsam,1);
    
    st = cputime;
    count=0;%count the number of loop
    ok=0;   %the flag to stop the loop 
    disp('number of samples against kkt:\n');

    
    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=zeros(numsam,1);      %the flag for the samples be selected for training 
       trainflag(xindex,1)=1;
       remaindindex=find(trainflag==0);%index for the unselected samples
    
      
       % Construct the Kernel matrix
       %fprintf('Constructing ...\n');
       H = zeros(trainsize,trainsize);  
       for i=1:trainsize
          for j=1:trainsize
              H(i,j) = trainY(i)*trainY(j)*svkernel(ker,trainX(i,:),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))*Y(remaindindex(l,1))*...
                        svkernel(ker,X(xindex(k,1),:),X(remaindindex(l,1),:));
              end
           end
           Qbn(k,1)=Qbn(k,1)*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 = Y(xindex,1)';, b = -1*(globalalpha(remaindindex,1)'*...
                  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);
       nsv = length(svi);
       
       % Implicit bias, b0
       b0 = 0;

      Hb=[];
         svii = find(globalalpha > epsilon & globalalpha < (C - epsilon));
         for i=1:length(svii)
            for j=1:length(svi)
               Hb(i,j) = Y(svii(i,1),1)*Y(svi(j,1),1)*svkernel(ker,X(svii(i,1),:),X(svi(j,1),:));
            end
         end
         if length(svii) > 0
            b0 =  (1/length(svii))*sum(Y(svii,1) - Hb*globalalpha(svi).*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);
       fremaindX=soutput(X(galphaindex,:),Y(galphaindex,1),globalalpha(galphaindex,1),ker,b0,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)=remaindindex;
       for k=1:numn
          if globalalpha(rmdflag(k,1),1)<=kktepsilon
             if Y(rmdflag(k,1),1)*fremaindX(k,1)<1-kktepsilon
                ok=0;
                rmdflag(k,2)=1;
                rmdflag(k,3)=1-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 (Y(rmdflag(k,1),1)*fremaindX(k,1))<(1-kktepsilon) ...
                | (Y(rmdflag(k,1),1)*fremaindX(k,1))>(1+kktepsilon)
                ok=0;
                rmdflag(k,2)=2;
                rmdflag(k,3)=abs(1-Y(rmdflag(k,1),1)*fremaindX(k,1));
             end
          end
          if globalalpha(rmdflag(k,1),1)>=C-kktepsilon
             if Y(rmdflag(k,1),1)*fremaindX(k,1)>1
                ok=0;
                rmdflag(k,2)=3;
                rmdflag(k,3)=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);
       agtkkt=find(rmdflag(:,2)>0);
       countzero=size(zeroalpha,1);
       countagtkkt=size(agtkkt,1)
       if countzero<countagtkkt
          if countzero<=transize
             tsize=transize;
          else
             tsize=countzero-1;
          end
       else
          tsize=countagtkkt;
       end
       tranindex=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)=rmdflag(tindex2,1);%implement the exchange .
       end
        %disp(count);
        %disp(xindex);
       xindex=sort(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