monqpcinfty.m

基于核分析的多类分类器

function [xnew, lambda, ind] = monqp0(H,c,A,b,l,verbose,X,ps,xinit)
% function [xnew, lambda, ind] = monqp0(H,c,A,b,l,verbose,X,ps)
%
% min 1/2  x' H x - c' x                
%  x 
% contrainte   A' x = b                     y'*x = 0 
% 
% et         0 <= x 
%
% Cr閐its : J.P. Yvon (INSA de Rennes) pour l'algorithme
%           O. Bodard (Clermont Ferrand) pour le debugage on line
%
%  S CANU - [e-mail]
%  Mai 2001

 
 
%-------------------------------------------------------------------------- 
%                                verifications 
%-------------------------------------------------------------------------- 
[n,d] = size(H); 
[nl,nc] = size(c); 
[nlc,ncc] = size(A); 
[nlb,ncb] = size(b); 
if d ~= n 
	error('H must be a squre matrix n by n'); 
end 
if nl ~= n 
	error('H and c must have the same number of row'); 
end 
 
if nlc ~= n 
	error('H and A must have the same number of row'); 
end 
if nc ~= 1 
	error('c must be a row vector'); 
end 
if ncb ~= 1 
	error('b must be a row vector'); 
end 
if ncc ~= nlb  
	error('A'' and b  must have the same number of row'); 
end 
 
if nargin < 5		% default value for the regularisation parameter 
	l = 0;				 
end; 
 
if nargin < 6		% default value for the display parameter 
	verbose = 0;				 
end; 

if exist('xinit','var')~=1
    xinit=[];
end;
 
 
fid = 1; %default value, curent matlab window 
%-------------------------------------------------------------------------- 
 

%-------------------------------------------------------------------------- 
%                       I N I T I A L I S A T I O N  
%-------------------------------------------------------------------------- 
 
OO = zeros(ncc);
H = H+l*eye(length(H)); % preconditionning

xnew = -1;ness = 0;
if isempty(xinit)
while (min(xnew) < 0 ) & ness < 100
   ind = randperm(n);
   ind = sort(ind(1:ncc)');
   aux=[H(ind,ind) A(ind,:);A(ind,:)' OO];
   aux= aux+l*eye(size(aux));
   newsol = aux\[c(ind) ; b];
   xnew = newsol(1:length(ind));
   ness = ness+1;
end;
x = xnew;
lambda = newsol(length(ind)+1:length(ind)+ncc);
else  % start with a predefined x

%     ind=find(xinit);
%        
%    aux=[H(ind,ind) A(ind,:);A(ind,:)' OO];
%     aux= aux+l*eye(size(aux));
%     newsol = aux\[c(ind) ; b];
%     xnew = newsol(1:length(ind));
%     lambda = newsol(length(ind)+1:length(ind)+ncc);
%     
%     x = xnew;
    
 
    ind=find(xinit>0);    
    x = xinit(ind);
    lambda=ones(ncc,1);
    
    
end;


%ind = 1;
%x = C/2; xnew = x;
%keyboard;
Hini = H;Aini=A;cini=c;bini=b;

eyen = l*eye(n);

%-------------------------------------------------------------------------- 
%                            M A I N   L O O P 
%-------------------------------------------------------------------------- 

Jold = 10000000000000000000;  
C = Jold; % for the cost function
if verbose ~= 0 
  disp('      Cost     Delta Cost  #support  #up saturate'); 
  nbverbose = 0; 
end 

STOP = 0;
while STOP ~= 1



indd = zeros(n,1);indd(ind)=ind;nsup=length(ind);indsuptot=[];
[J,yx] = cout(H,x,b,C,indd,c,A,b,lambda); 
if verbose ~= 0 
  nbverbose = nbverbose+1; 
  if nbverbose == 20 
  disp('      Cost     Delta Cost  #support  #up saturate'); 
  nbverbose = 0; 
end 
if Jold == 0 
   fprintf(fid,'| %11.4e | %8.4f | %6.0f | %6.0f |\n',[J (Jold-J) nsup length(indsuptot)]); 
   elseif  (Jold-J)>0
fprintf(fid,'| %11.4e | %8.4f | %6.0f | %6.0f |\n',[J min((Jold-J)/abs(Jold),99.9999) nsup length(indsuptot)]); 
   else
fprintf(fid,'| %11.4e | %8.4f | %6.0f | %6.0f | bad mouve \n',[J max((Jold-J)/abs(Jold),-99.9999) nsup length(indsuptot)]); 
end 
end 
Jold = J; 
 
% nouvele resolution du syst鑝e

% newsol = [H(ind,ind)+l*eye(length(ind)) A(ind,:);A(ind,:)' OO]\[c(ind) ; b];
% xnew = newsol(1:length(ind));
% lambda = newsol(length(ind)+1:length(newsol));

% %keyboard
% U = chol(H(ind,ind));      % CHOLESKI INSIDE!
% At=A(ind,:)';
% Ae=A(ind,:);
% ce = c(ind);
% Ut = U';
% M = At*(U\(Ut\Ae));
% d = U\(Ut\ce);
% d = (At*d - b);    % second membre  (homage au gars OlgaZZZZZZ qui nous a bien aid