monqp.m

支持向量机的核心程序

function [xnew, lambda, pos,mu] = 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 <6		% default value for the regularisation parameter 
    l = 0;				 
end; 

if nargin <7		% default value for the display parameter 
    verbose = 0;				 
end; 
if  size(C,1)~=nl
 %   keyboard
    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*ones(size(C));ness = 0;     
ind=1;

%keyboard

if isempty(xinit)

    while (( sum (xnew < 0)) >0 | ( sum(xnew >C(ind)) > 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-12
            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;


%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)]); 
            
        end 
    end 
    
    if isempty(ind)
        xnew=0;
        lambda=0;
        pos=1;
        return
    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

       %
       %  on remplace tout les x qui sont indsuptot par C
       %
        
       Cmat= ones(length(ind),1)*(C(indsuptot)'); 
       if size(ce)~= size(sum( Cmat.*H(ind,indsuptot),2))
            keyboard;
        end;
       ce= ce - sum( Cmat.*H(ind,indsuptot),2);%                              min     x'Hx + c'x
       Cmat= C(indsuptot)*ones(1,size(A,2)); 
       be= be - sum(Cmat.*A(indsuptot,:),1)';                            %    with    Ax=b

    end 
    


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

        M = At*(U\(U'\Ae));
        d = U\(U'\ce);
        d = (At*d - be);   % second membre  (homage au gars OlgaZZZZZZ qui nous a bien aid?
        % On appelle le gc fabuleux de mister OlgaZZZ Hoduc
        % lambdastart = rand([m,1]);
        % [lambda] = gradconj(lambdastart,M*M,M*d,MaxIterZZZ,tol,1000);
        if rcond(M)<l
            M=M+l*eye(size(M));
        end;
        lambda = M\d;
        % On reconstitue le vecteur en r閟olvant le syst鑝e lin閍ire Hx = c-Alambda
       % size(lambda)
       % size(Ae)
 %      keyboard
%        xnew = U\(Ut\(ce-Ae*lambda));
        xnew=auxH\(ce-Ae*lambda);
    else
        
        % if non definite positive;
        
        
        auxM=At*(auxH\Ae);
        M=auxM'*auxM;
        d=auxH\ce;
        d=At*d-be;
        if rcond(M)<l
            M=M+l*eye(size(M));
        end;
        lambda=M\d;
        xnew=auxH\(ce-Ae*lambda);
    end;
    
    %-------------------------------------------------------------------------------------------------------
  
    [minxnew minpos]=min(xnew);
    
    %if (min(xnew) < 0  | max(xnew) > C) &  length(ind) > ncc
    
    if (sum(xnew< 0) > 0 | sum(xnew > C(ind)) > 0) &  length(ind) > ncc  %% 03/01/2002
    
    % keyboard
    % on projette dans le domaine admissible et on sature la contrainte associ閑
        %-------------------------------------------------------------------------
        d = xnew - x;                % projection direction
        indad = find(xnew < 0); 
        indsup = find(xnew > C(ind) );                      %% 03/01/2002
        [tI indmin] = min(-x(indad)./d(indad));
        [tS indS] = min((C(ind(indsup))-x(indsup))./d(indsup)); % pas de descente 
        if isempty(tI) , tI = tS + 1; end; 
        if isempty(tS) , tS = tI + 1; end; 
        t = min(tI,tS); 
        
        x = x + t*d;                 % projection into the admissible set
        
        if t == tI 
            ind(indad(indmin)) = [];            % the associate variable is set to 0
            x(indad(indmin)) = [];              % the associate variable is set to 0
        else
            indsuptot = [indsuptot ; ind(indsup(indS))]; 
            ind(indsup(indS))= []; 
            x(indsup(indS))= []; 
        end
        
    else
        xt = zeros(n,1);           % keyboard;
        xt(ind) = xnew;              % keyboard;
        xt(indsuptot) = C(indsuptot); indold=ind;             % keyboard;                                     %% 03/01/2002
        
        mu = H*xt - c + A*lambda;    % calcul des multiplicateurs de lagrange associ閑s aux contraintes
        
        indsat = 1:n;                           % on ne regarde que les contraintes satur閑s
        indsat([ind ; indsuptot]) = [];
        
        [mm mpos]=min(mu(indsat));
        [mmS mposS]=min(-mu(indsuptot));  %keyboard
        
        if (~isempty(mm) & (mm < -sqrt(eps)))  | (~isempty(mmS) & (mmS <  -sqrt(eps))) 
            if isempty(indsuptot) | (mm < mmS)             % il faut rajouter une variable
                ind = sort([ind ; indsat(mpos)]);        % on elimine une contrainte de type x=0
                x = xt(ind);
            else
                ind = sort([ind ; indsuptot(mposS)]);  % on elimine la contrainte sup si necessaire
                x = xt(ind);                           % on elimine une contrainte de type x=C
                indsuptot(mposS)=[];
            end
        else
 
            STOP = 1;                % AILLLET !
            pos=sort([ind ; indsuptot]);
            xt = zeros(n,1);             % keyboard;
            xt(ind) = xnew;              % keyboard;
            xt(indsuptot) = C(indsuptot);           % keyboard;
            indout = 1:n;
            indout(pos)=[];
            xnew = xt;
            xnew(indout)=[];
              
        end
        
    end
end

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