irlssvmf.m

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

function [pall,alphaold,lho]=irlssvmf(X,y,alpha,pars)

K=length(unique(y));c=K-1;
[N,d]=size(X);
lho=-Inf;

R=eye(c*d);
bufferlhold=-Inf*ones(pars.deltak,1);
stop1=Inf;
itt=0;
%initialise probs
if sum(alpha.^2)==0 || length(alpha)~=c*N,
  alpha=zeros(c*N,1);
  pall=(1/K)*ones(N,c);  
  if (pars.info==1), disp('reset1');end
else
disp('controle op correctheid')
  alphaz=zeros(c*N,1);
  alphamatz=reshape(alphaz,N,c); 
  pallz=(1/K)*ones(N,c); 
  ff = zeros(N*c,1);
  %lhnz=lhlambda(y,wmatz,pallz,pars);
  lhnz = lhlambdaAlpha(y,alphamatz,pallz,ff,pars)
  alphamat=reshape(alpha,d,c); 
  %pall=probs(wmat,X);
  [pall,ff]=alphaprobs(X,y,alpha,pars,X);
  lhnp=lhlambda(y,alphamat,pall,ff,pars);
  if lhnp<lhnz,
      alpha=alphaz;
      pall=pallz;
      if (pars.info==1), disp('reset2');end
  end
end

alphaold=zeros(N*c,1); 
alphaomat=zeros(N,(K-1)); 
eta=1;ib=1;con=1;
while (itt<pars.maxitt002) && (pars.epsil<stop1) && con, %WHILE 1
  itt=itt+1;
  
  if pars.redds == 1 ,
    %exclude datapoints with no informative value
    %indexi=selectInput(X,y,pall,pars);
    disp('pars.redds not implemented for normale KLOGREG');
  end  
  %use all points
  indexi=1:size(X,1);
  N=length(indexi);

  ytemp=y(indexi);
  pallt=pall(indexi,:);
  
  zr=zeros(c*N,1);Yr=[];
  weightr=sparse(zeros(c*N,c*N));
  for kk=1:N,
     Blocp=pallt(kk,:)';
     BlocW=-Blocp*Blocp';
     BlocW=BlocW+diag(sparse(Blocp)); % 1 block of weight matrix\
     yl=zeros(K-1,1);
     if ytemp(kk)<K
         yl(ytemp(kk))=1;
     end
     Yr=[Yr ; yl];
     weightr(c*(kk-1)+(1:c),c*(kk-1)+(1:c))=BlocW;
     zr((kk-1)*c+(1:c))=BlocW \ (yl - Blocp); 
   end;

   Xr=[];
   Xr=sparse(zeros(N*c,c*d)); % matrix 
   for cc=1:c,
     row=(0:1:N-1)*c+cc;
     col=d*(cc-1)+(1:1:d);
     Xr(row,col)=X(indexi,:);
   end

   omegaSmall = kernel_matrix(Xr,pars.ker,pars.sig2);
   
   %for kernel matrix
   FZ=zeros(N*c,c*(d)); % matrix 
   for cc=1:c,
     row=(0:1:N-1)*c+cc;
     col=d*(cc-1)+(1:1:d);
     FZ(row,col)=ones(N,d);
   end
   FZ=(FZ*FZ')/d;
   
   omegaSmall=FZ.*omegaSmall; %force kernelmatrix to have zeros on predefined places

   Hr = omegaSmall + (N*pars.lambda)*inv(weightr);

   [crap,ff]=alphaprobs(X,y,alphaomat,pars,X);

  
   ff=reshape(ff',N*c,1); 
   alphanew=Hr \ (zr+ff); 
   alphanmat=reshape(alphanew,K-1,N)'; % matrix notation
   
   [palln,ffn]=alphaprobs(X,y,alphanmat,pars,X);

   lhn = lhlambdaAlpha(y,alphanmat,palln,ffn,pars);
   
   if lhn>=lho,
    alphaomat=alphanmat;
    alphaold=alphanew;
    lho=lhn;  
    pall=palln; 
    stop1= abs(bufferlhold(end) - lhn ) / abs(lhn);  
    bufferlhold = [lhn; bufferlhold(1:end-1)];
    if (pars.info==1), disp(lhn); end
   else
       con =0; 
   end
 
end % WHILE 1