monqpold1.m

surpport vector machine,matlab

function [xnew, lambda, pos] = monqp(H,c,A,b,C,l,verbose,X,ps,xinit)
% function [xnew, lambda, pos] = monqp(H,c,A,b,C,l,verbose,X,ps,xinit)
%
% min 1/2  x' H x - c' x                 
%  x 
% contrainte   A' x = b                    
% 
% et         0 <= x_i  <= C_i
%
% 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  length(C)~=nl
    C=ones(nl,1)*C;  % vectorizing the UpperBound Constraint
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;     

%keyboard


if isempty(xinit)
    while (( sum (xnew < 0)) >0 | ( sum(xnew >C) > 0))  & ness < 100   %% 03/01/2002
        ind = randperm(n);
        ind = sort(ind(1:ncc)');
        aux=[H(ind,ind) A(ind,:);A(ind,:)' OO];
        aux= aux+l*eye(size(aux));
        if rcond(aux)>1e-5
            newsol = aux\[c(ind) ; b];
            xnew = newsol(1:length(ind));
            ness = ness+1; 
        else
            ness=101;
            break;
        end;
    end;
    %keyboard
    if ness < 100
        x = xnew;
        lambda = newsol(length(ind)+1:length(ind)+ncc);
    else
      %  keyboard
        ind = 1;
        x = C(ind)/2.*ones(size(ind,1),1);                               %% 03/01/2002
        %lambda = newsol(length(ind)+1:length(ind)+ncc);
        lambda=ones(ncc,1);    
       
    end
    indsuptot = [];	 
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;
    %  keyboard
    
    indsuptot=find(xinit==C);
    indsup=indsuptot;
    ind=find(xinit>0 & xinit ~=C);    
    x = xinit(ind);
    lambda=ones(ncc,1);
    
end;
%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);indd(indsuptot)=-1;
    [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)]); 
            %keyboard    
        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));
    
    ce = c(ind);
    be = b;
    
%     if (isempty(indsuptot)==0) % if NOT empty
%         keyboard
%         %     ce= ce - C*(sum(H(ind,indsuptot),2)+sum(H(indsuptot,ind),1)');%    min     x'Hx + c'x
%         ce= ce - C*(sum(H(ind,indsuptot),2));%                              min     x'Hx + c'x
%         be= be - C*sum(A(indsuptot,:),1)';                            %    with    Ax=b
%     end 
    

%03/01/2003
    if (isempty(indsuptot)==0) % if NOT empty
     %   keyboard
        %     ce= ce - C*(sum(H(ind,indsuptot),2)+sum(H(indsuptot,ind),1)');%    min     x'Hx + c'x
        ce= ce - min(C)*(sum(H(ind,indsuptot),2));%                              min     x'Hx + c'x
        %be= be - min(C)*sum(A(indsuptot,:),1)';                            %    with    Ax=b
        %be=be;
    end 
    


 %   keyboard
    auxH=H(ind,ind);
    [U,testchol] = chol(auxH);
    At=A(ind,:)';
    Ae=A(ind,:);
    
    if testchol==0
        Ut = U';

        M = At*(U\(Ut\Ae));
        d = U\(Ut\ce);
        d = (At*d - be);   % second membre  (homage au gars OlgaZZZZZZ qui nous a bien aid