lssvcbs.m

The goal of SPID is to provide the user with tools capable to simulate, preprocess, process and clas

function  [alphaTot,b,sig2, gam_opt, mu_opt, zeta_opt, L3_opt,neffPar]=lssvcBaySps(ptrn, ttrn, knlType, maxstep, dispOpt, sigs)
%function  [alphaTot,b,sig2, gam_opt, mu_opt, zeta_opt]=lssvcBaySps(ptrn, ttrn, knlType, maxstep, dispOpt)
% knlType: 'rbf', or 'lin'
% maxstep: # steps for sparse approximation (i.e. the pruning of data points with negative alpha value)
% dispOpt<0 : not display any msg on sigma and gamma selelction

if nargin<3,
    knlType='rbf';
end,

if nargin<4,
   maxstep=20;
end,

if nargin<5,
    dispOpt=0;
end,

if nargin<6,
    sigs=[];
end
 
switch knlType,
case 'rbf',
   knl=1;
   if isempty(sigs),
        nlsig=4; %sigs=[];
        sig_opt=1; 
    else
        nlsig=1; 
        sig_opt=sigs(1); 
    end
   gam_opt=1;
case 'lin',
   knl=2;
   nlsig=1; sigs=[1];
   sig_opt=1; gam_opt=1;
case 'pol',
    knl=2;
    if isempty(sigs),
        sigs=[1:4];
    end
    nlsig=1; gam_opt=1;

case 'cwrbf',
    knl=-1; gebr_kernel='RBFaddkernel';
   if isempty(sigs),
        nlsig=4; 
        sig_opt=1; 
    else
        nlsig=1; 
        sig_opt=sigs(1); 
    end
   gam_opt=1;
    
case 'cwlin'
    knl=-2; gebr_kernel='LINaddkernel';
    nlsig=1; sigs=[1];

case 'cwpol',
     knl=-2; gebr_kernel='LINaddkernel';
    if isempty(sigs),
        sigs=[1:4];
    end
    nlsig=1; gam_opt=1;
end,


ntot=size(ptrn,1);
if dispOpt<0,
    disp2=-1;
else
    disp2=0;
end,

st=cputime;
global logKK0 KK0,
alphaTot=ones(ntot,1);
for istep=1:abs(maxstep),
   
    idxM=find(alphaTot>0); 
    if dispOpt>=0,
    disp(['step ', num2str(istep), ' ; #data: ', num2str(length(idxM))]),
    end,
    Xs=ptrn(idxM,:); Ys=ttrn(idxM);
    
   [sig_opt, gam_opt, L3_opt, neffPar, mu_opt, zeta_opt]=...
      optSigGamLssvmBay(Xs, Ys, sigs, nlsig, disp2, knl);
  
%gam0=3; sig0=0.1;
%global gradType, gradType.gam='on'; gradType.sig='on';
%[sig_opt, gam_opt, L3_opt, neffPar, mu, zeta]  = optim_gamsigBayLSSVMard(gam0, sig0, Xs, Ys,knl,10, 1)
  if dispOpt>0, disp([sig_opt, gam_opt, L3_opt, neffPar, mu_opt, zeta_opt]), end,
  if istep==1,
      neffPold=neffPar;
  else
      rneff=neffPar/neffPold*ntot/length(Ys);
      if dispOpt>0,
      disp(['rneff(neffPnew/neffPold*Ntot/Npruned):', num2str(rneff)]),
      end,
%      if  rneff< 0.1 |rneff>5,
      if  rneff< 0.001 |rneff>5,
          if dispOpt>0,          
              disp(['Pruning end (rneff<<)!']),
          end
          break,
      end,
  end,
  
   pK=logKK0;
   if knl==1,
      sig2=sig_opt*sig_opt; 
      KK=exp(pK/sig2);
      K=Ys*Ys'.*KK;
  elseif knl==2,
      sig2=sig_opt;
      KK=pK.^sig2;
      K=Ys*Ys'.*KK;
  elseif knl<0,
      if knl==1,
          sig2=sig_opt*sig_opt;
      else
          sig2=sig_opt;
      end
      logKK0=[]; pK=[];
      for i=1:size(Xs,1),
          for j=1:i,
              KK(i,j) = eval([gebr_kernel,'(Xs(i,:),Xs(j,:),sig_opt)']);    
              KK(j,i) = KK(i,j);              
          end
      end
      K=Ys*Ys'.*KK;
   end,  
   if istep==1,
%       save MATkkk KK, 
       KK0=pK;
   end,
   clear KK,
   
    Nm=length(Ys);
    Ic=eye(Nm); vec1=ones(Nm,1);
    %Igam=1/gam_opt*ones(Nm,1); Igam(find(Ys>0))=1e-5; M=[0 Ys';Ys (K+diag(Igam))];
	M=[0 Ys';Ys (K+Ic/gam_opt)];
    alphab=M\[0;vec1];
    alpha=alphab(2:length(alphab));
    b=alphab(1);
    
    alpha0=alpha; Km=K; epsa=0;
    for jstep=1:maxstep,
        if isempty(find(alpha0<=epsa)),
            break,
        end,
        idxM0=find(alpha>epsa); 
        Ym=Ys(idxM0);
        nsv=length(idxM0);
    	Ic=eye(nsv); vec1=ones(nsv,1);
        Km=K(idxM0, idxM0);
    	M=[0 Ym';Ym (Km+Ic/gam_opt)];
       alphab0=M\[0;vec1];
       alpha0=alphab0(2:length(alphab0));  
       b=alphab0(1);
   
       alpha(idxM0)=alpha0;
    end, % for jstep
      
    if isempty(find(alpha<=epsa)),
       alphaTot(idxM)=alpha; 
       if dispOpt>0,
       disp('Pruning end!'),
       end
       break,
    else
       if maxstep>0
       alpha(find(alpha<=epsa))=0; 
       end
       alphaTot(idxM)=alpha;
    end     
end, % for istep

if dispOpt>0,
    fprintf('Execution time: %4.1f seconds\n',cputime - st);
    w2 = alpha0'*Km*alpha0;
    fprintf('|w0|^2    : %f\n',w2);
    fprintf('Margin    : %f\n',2/sqrt(w2));
    fprintf('Sum alpha : %f\n',sum(alpha0));
end,
return,