svmclassnpa.m

SVM algorithm in matlab.

function  [xsup,alpha,b,pos]=svmclassnpa(x,y,C,kernel,kerneloption,verbose);

% USAGE 
%   [xsup,alpha,b,pos]=svmclassnpa(x,y,C,kernel,kerneloption,verbose); 
%
%
%    Main ROUTINE For Nearest Point Algorithm
%
%   This function solve the SVM Classifier problem by 
%   transforming the soft-margin problem into a hard margin problem
%   by means of a slight modification of the kernel and the
%   introduction of a quadratic penalization term.
%
%   The problem is equivalent  to look for the nearest points of 
%   two convex hulls. This is solved by the point of minimum norm is the 
%   Minkowski set (U-V) where U is the class 1 set and V the class 2 set 
%   
%   This method is deeply described is 
%   
%   Tech Rep : A Fast Iterative NPA for SVM Classifier
%   S. Keerthi et al
% 

%   Last Modification : 31/03/2001  A. R
%


global  n;
global  m;
% global  C;
% global  kernel;
% global  kerneloption;
% global  x;
% global  y;
global  deluu;
global  delvv;
global  deluv;
global  delzz;
global  eu;
global  ev;
global  alind;
global  beta;




%  Initialization


eps1=0.001; %type I loop
eps=0.001; %type II loop

m=length(y);
beta=zeros(1,m);
eu=zeros(7,m);
ev=zeros(1,m);

indpos=find(y==1);
indmoins=find(y==-1);
i=indpos(1);
j=indmoins(1);
beta(i)=1;
beta(j)=1;
alind(i)=1;
alind(j)=1;
nsupp=2;
deluu=svmkernel(x(i,:),kernel,kerneloption,x(i,:))+1/C;
delvv=svmkernel(x(j,:),kernel,kerneloption,x(j,:))+1/C;
deluv=svmkernel(x(i,:),kernel,kerneloption,x(j,:));
delzz=deluu+delvv-2*deluv;
eu(i)=deluu;
eu(j)=deluv;
ev(i)=deluv;
ev(j)=delvv;
Examine_NonSV=1;SV_Optimality=1;NonSV_Optimality=0;toohigheps=1;
Num_Type4_Updates=0;Num_Type2_Updates=0;

%-------------------------------------------------------------------%
%                                                                   %
%                               MAIN LOOP                           %
%                                                                   %
%-------------------------------------------------------------------%

while (SV_Optimality==0 | NonSV_Optimality==0 ) 
    ind=find(beta<=0);
    alind(ind)=zeros(size(ind));
    beta(ind)=zeros(size(ind));
    if Examine_NonSV==1             % TYPE 1 LOOOP
        Num_Type1_Updates=0;
        NonSV_Optimality=1;
        Examine_NonSV=0;
        ind_NonSV=find(alind==0);
        
        for k=(ind_NonSV)
            ind_I=intersect(find(y==1),find(alind==1));       % Positive class
            ind_J=intersect(find(y==-1),find(alind==1));      % Negative class  
           % keyboard
            eu(k)=sum(beta(ind_I).*(svmkernel(x(k,:),kernel,kerneloption,x(ind_I,:))+(ind_I==k)./C));
            ev(k)=sum(beta(ind_J).*(svmkernel(x(k,:),kernel,kerneloption,x(ind_J,:))+(ind_J==k)./C));
            zdk=eu(k)-ev(k);
            zdu=deluu-deluv;
            zdv=deluv-delvv;  
            if ( (y(k)==1 & (zdu-zdk)>=0.5*eps1*delzz ) | (y(k)==-1 & (zdk-zdv)>=0.5*eps1*delzz) )
                NonSV_Optimality=0;
                beta(k)=0;
                alind(k)=1;
                success=Take_Step(k,x,y,C,kernel,kerneloption);
                if delzz<0         
                    keyboard
                end;
                if success==0
                    alind(k)=0;
                else
                    Num_Type1_Updates=Num_Type1_Updates+1;
                end;
            end;
            
        end;
        
        
    else
        Examine_NonSV=1;            % Support Vector Processing....  
        nsupp_old=nsupp;
        nsupp=sum(alind==1);
        if abs(nsupp_old-nsupp) >=0.05*nsupp
            Max_Updates=m;
        else
            Max_Updates=10*m;
        end;
        Loop_Completed=0;
        Num_Type2_Updates=0;
        nbloop=0;
        while(Loop_Completed==0)
            nbloop=nbloop+1;
            zdu=deluu-deluv;
            zdv=deluv-delvv;
            ind_I=intersect(find(y==1),find(alind==1));       % Positive SV
            ind_J=intersect(find(y==-1),find(alind==1));      % Negative SV
            zdi=eu(ind_I)-ev(ind_I);
            [maxi imax]=max(zdu-zdi);
            imax=ind_I(imax);
            zdj=eu(ind_J)-ev(ind_J);
            [maxi jmax]=max(zdj-zdv);
            jmax=ind_J(jmax);
            gu=zdu-eu(imax)+ev(imax);
            gv=eu(jmax)-ev(jmax)-zdv;
            if (gu<=0.5*eps*delzz) & (gv <=0.8*eps*delzz)
                Loop_Completed=1;
                SV_Optimality=1; 
            else
                if (gu >=gv) 
                    kmax=imax;
                else
                    kmax=jmax;
                end;
                %[success,cas]=Take_Step(kmax);
                success=Take_Step(kmax,x,y,C,kernel,kerneloption);
                if delzz<0          
                    keyboard
                end;
                if success==0
                    Loop_Completed=1;
                    SV_Optimality=0;
                else
                    
                    Num_Type2_Updates=Num_Type2_Updates+1;
                end;
                if Num_Type2_Updates>Max_Updates
                    Loop_Completed=1;
                    SV_Optimality=0;
                end;
            end;
            
        end;
    end;
    
    if SV_Optimality & NonSV_Optimality  
        fprintf('Optimality criteria Satisfied...\n');
 
        pos=find(alind==1);
        xsup=x(pos,:);
        ysup=y(pos,:);
        if ~isempty(imax) & ~isempty(imax)
            h_U = ev(imax) - eu(imax);
            h_V = eu(jmax) - ev(jmax);
            gamma = 2.0 / (-h_U - h_V);
            b  = (h_V - h_U) / (h_V + h_U);
            alpha=beta(pos).*gamma.*y(pos)';
            alpha=alpha(:);
        end;
        return
    end;
    
    if Num_Type1_Updates==0 & Num_Type2_Updates==0 
        fprintf('Algorithm has not converged...\n');
        pos=find(alind==1);
        xsup=x(pos,:);
        ysup=y(pos,:);
        if ~isempty(imax) & ~isempty(imax)
            h_U = ev(imax) - eu(imax);
            h_V = eu(jmax) - ev(jmax);
            gamma = 2.0 / (-h_U - h_V);
            b  =(h_V - h_U) / (h_V + h_U);
            alpha=beta(pos).*gamma.*y(pos)';
            alpha=alpha(:);
        end;
        return
    end;
    
    
    if ~isempty(verbose) & verbose==1
    fprintf(' ||z||= %f SV_Optimality:%d  NonSV_Optimality:%d\n',delzz,SV_Optimality, NonSV_Optimality); 
end;
    
    
end;



pos=find(alind==1);
xsup=x(pos,:);
ysup=y(pos,:);
h_U = ev(imax) - eu(imax);
h_T = eu(jmax) - ev(jmax);
gamma = 2.0 / (-h_U - h_V);
b  =(h_V - h_U) / (h_N + h_U);
alpha=(beta(pos).*gamma.*y(pos)');
alpha=alpha(:);



%---------------------------------------------------------------------%
%                                                                     %  
%                           Function Take_Step                        %  
%                                                                     %  
%---------------------------------------------------------------------%
function [success,cas]=Take_Step(kmax,x,y,C,kernel,kerneloption)

global   n;
global  m;
% global C;
% global kernel;
% global kerneloption;
% global  x;
% global y;
global   deluu;
global   delvv;
global   deluv;
global   delzz;
global   eu;
global   ev;
global   alind;
global   beta;


% compute kmin%
success=0;
zdu=deluu-deluv;
zdv=deluv-delvv;
ind_betapos=find(beta>=0);
ind_I=intersect(find(y==1),find(alind==1));      % Positive SV
ind_J=intersect(find(y==-1),find(alind==1));     % Negative SV   
ind_Ibetapos=intersect(ind_I,ind_betapos);
ind_Jbetapos=intersect(ind_J,ind_betapos);
[minii,imin]=min(zdu-eu(ind_Ibetapos)+ev(ind_Ibetapos));
imin=ind_Ibetapos(imin);
[minij,jmin]=min(-zdv+eu(ind_Jbetapos)-ev(ind_Jbetapos));
jmin=ind_Jbetapos(jmin);
if (minii)<(minij)
    kmin=imin;
else
    kmin=jmin;
end;








% Do Modified Gilbert Step if beta(kmin)=1

if (beta(kmin)>=1 ) | kmin==0        
    cas=5;
    success=Gilbert_Step(kmax,x,y,C,kernel,kerneloption);
    return;
end;
mu=beta(kmin)/(1-beta(kmin));
%kmin
%kmax
if isempty(kmin)
    kmin=kmax;
end;
%kmin
%kmax
%size(x)
psi = svmkernel(x(kmin,:),kernel,kerneloption,x(kmax,:))+(kmin==kmax)./C;
psi1= svmkernel(x(kmax,:),kernel,kerneloption,x(kmax,:))+(kmax==kmax)./C;
psi2= svmkernel(x(kmin,:),kernel,kerneloption,x(kmin,:))+(kmin==kmin)./C;
eumax = eu(kmax);
evmax = ev(kmax);
eumin = eu(kmin);
evmin = ev(kmin);

if y(kmax)==1 & y(kmin)==1
    cas=1;
    imax=kmax; imin=kmin;
    d11=delzz;
    t1=deluu - eumin - deluv + evmin;
    d22 = delzz + mu*mu*(deluu + psi2 - 2.0*eumin) + 2.0*mu*t1;
    d33 = psi1 + delvv - 2.0*evmax;
    d12 = delzz + mu*t1;
    d13 = eumax - deluv - evmax + delvv;
    d23 = d13 + mu*(eumax - deluv - psi + evmin);
    [d,la1,la2,la3,flag]=triangle(d11,d22,d33,d12,d13,d23);
    if d >=delzz | d < 0
        success=0;
        return;
    end;
    success=1;
    delzz=d;
    la1b = la1 + la2 + la2*mu;
    la2b = -la2*mu;
    la3b = la3;
    deluu = la1b*la1b*deluu + la2b*la2b*psi2 + la3b*la3b*psi1 + 2.0*la1b*la2b*eumin ;
    deluu = deluu + 2.0*la1b*la3b*eumax +2.0*la2b*la3b*psi;
    deluv = la1b*deluv + la2b*evmin + la3b*evmax;
    ind_SV=find(alind==1);
    kkmin=svmkernel(x(imin,:),kernel,kerneloption,x(ind_SV,:))+(ind_SV==imin)./C;
    kkmax=svmkernel(x(imax,:),kernel,kerneloption,x(ind_SV,:))+(ind_SV==imax)./C;
    eu(ind_SV)=la1b*eu(ind_SV)+la2b*kkmin+la3b*kkmax;
    ind_I=intersect(find(y==1),find(alind==1));
    beta(ind_I)=la1b*beta(ind_I);
    beta(imax)=beta(imax)+la3b;
    if (flag ==2 | flag== 5 | flag == 3) 
        alind(imin)=0;
        beta(imin)=0;
    else
        beta(imin) = beta(imin) + la2b;
    end;
end;
if y(kmax)==-1 & y(kmin)==-1
    cas=2; 
    jmax=kmax;jmin=kmin;
    d11 = delzz;
    t1 = delvv - evmin - deluv + eumin; 
    d22 = delzz + mu*mu*(delvv + psi2 - 2.0*evmin) + 2.0*mu*t1;
    d33 = psi1 + deluu - 2.0*eumax;
    d12 = delzz + mu*t1;
    d13 = evmax - deluv - eumax + deluu;
    d23 = d13 + mu*(eumin - deluv - psi + evmax);
    [d,la1,la2,la3,flag]=triangle(d11,d22,d33,d12,d13,d23);
    if d >=delzz | d < 0
        success=0;
        return;
    end;
    success = 1;
    delzz = d;
    la1b = la1 + la2 + la2*mu;
    la2b = -la2*mu;
    la3b = la3;
    delvv = la1b*la1b*delvv + la2b*la2b*psi2 + la3b*la3b*psi1 +   2.0*la1b*la2b*evmin;
    delvv=delvv+ 2.0*la1b*la3b*evmax +  2.0*la2b*la3b*psi;
    deluv = la1b*deluv + la2b*eumin + la3b*eumax;
    ind_SV=find(alind==1);
    kkmin=svmkernel(x(jmin,:),kernel,kerneloption,x(ind_SV,:))+(ind_SV==jmin)./C;
    kkmax=svmkernel(x(jmax,:),kernel,kerneloption,x(ind_SV,:))+(ind_SV==jmax)./C;
    ev(ind_SV)=la1b*ev(ind_SV)+la2b*kkmin+la3b*kkmax;
    ind_J=intersect(find(y==-1),find(alind==1));
    beta(ind_J)=la1b*beta(ind_J);
    beta(jmax)=beta(jmax)+la3b;
    if (flag ==2 | flag== 5 | flag == 3) 
        alind(jmin)=0;
        beta(jmin)=0;
    else
        beta(jmin) = beta(jmin) + la2b;
    end;
end;
if y(kmax)==1 & y(kmin)==-1
    cas=3; 
    imax=kmax;jmin=kmin;