lar.m

基于核分析的多类分类器

function [solution, solution_OLS] = LAR(x,y, AllBorne, Limites, lambda, verbose, xappforplot)
%
% USAGE :
%  [solution, solution_OLS] = LAR(x,y, AllBorne, Limites, lambda, verbose, xappforplot)
%
%  Looks for a regression solution : y_hat = x Beta +b
%  by solving : min ||y_hat-y||^2 
%               sum(abs(Beta)) < t
%
%  The algorithm is based on the LARS of Efron et al.
%  Stopping criterion and multiple kernel formulation are based on
%  V. Guigue's paper untitled Kernel Basis Pursuit (google it!).
%
% ARGUMENTS :
%   verbose    : plot (>0) or no (==0) lfigures.
%               figures are numbered from verbose
%   x          : input data
%   y          : output data
%   AllBorne :  This paraameter is a structure that deals with the stopping
%               criterion. One can combine several kinds of stopping
%               criterion for the processing of the regularization path.
%
%                    Allborne{i}.type : type of stopping criterion 
%                    Allborne{i}.borne : bound value on this stopping
%                                       criterion
%                                     
%                   Here we give the possible type of stopping criterion
%
%                   'sumBeta'  ->  sum(abs(Beta)) < borne_Beta
%                       *  Allborne.borne : vector of bounds. The algorithm
%                       will output as many solutions as bounds.
%                       
%                       (usual bound as defined in the LARS algorithm)
%
%                   'pcSV'     ->  borne_Beta = percent of selected support
%                               vectors
%                       *  Allborne.borne : vector of bounds. The algorithm
%                       will output as many solutions as bounds.

%                   'nbSV'     ->  borne_Beta = nb of support vector
%                       *  Allborne.borne : vector of bounds. The algorithm
%                       will output as many solutions as bounds.

%                   'trapscale'-> scale or kernel representing noise 
%                       * structure.indTS : index of the kernel trapscale
%                       in the multiple kernel setting. It should be one of
%                       the last.                   
%                       *  Allborne.borne : vector of bounds. The algorithm
%                       will output as many solutions as bounds. (counts
%                       the number of time the trapscale has been
%                       selected.)
%
%                   'RVtrap'   ->  Automatically add some random noise as
%                   information sources.
%                       *  Allborne.borne : vector of bounds. The algorithm
%                       will output as many solutions as bounds. (counts
%                       the number of time the RV has been selected).
%
%                   'countback' -> allowed backward movement
%                       *  Allborne.borne : vector of bounds. The algorithm
%                       will output as many solutions as bounds.
%
%   Limites         * structure.nbmaxMoveBack : maximal number of backward
%                     movement for the LASSO.
%                   * structure.nbmaxSV : maximal number of SV allowed.
%                       this is useful in a Large Scale problem in order to
%                       avoid memory crash.
%
%   lambda : regularisation term
%
% OUTPUT:
%   solution  : all the solutions in a structures, the beta are only the
%                   non-zeros one
%                   solution{i}.b   
%                   solution{i}.Beta
%                   solution{i}.indxsup
%                   solution{i}.type % de solution !
%                   solution{i}.borne
%
%   solution_OLS : same solution  but the last step is based on the OLS not
%   on LARS.
%
%





% V. Guigue, A. Rakotomamonjy 12/2004


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PRE TEST

if nargin<6
    verbose = 0;
end

nbvar    =  size(x,2);
nbptsapp =  size(x,1);
% if sum(std(x)) > nbvar+0.5 | sum(std(x)) < nbvar-0.5
%     %error('LARS method should be used on normalized data\n');   
%     fprintf('LARS method should be used on normalized data\n');
% end

 if sum(mean(x))>1e-3;
     fprintf('LARS method should be used on normalized data\n');
 end

nbexptot = 0;
nbborne  = length(AllBorne);
for i=1:nbborne,
    nbexptot = nbexptot + length(AllBorne{i}.borne);
    
    if strcmp(AllBorne{i}.type,'trapscale') % transformation de l'indice de debut TS
        AllBorne{i}.indTS = (AllBorne{i}.indTS-1)*nbptsapp;   
    end
    
end

% fabrication d'un ensemble de point aleatoire :
RVset = randn(nbptsapp,15)';


nbexpcur = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AJOUT DE B NORMALISE

x = [x ones(size(x,1),1)/sqrt(nbptsapp)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMTRES PAR DEFAUT


test_OLS = 1;
eps = 1e-12;

if isfield(Limites,'nbmaxMoveBack') == 0
    Limites.nbmaxMoveBack = 1e8;
end
if isfield(Limites,'nbmaxSV') == 0
    Limites.nbmaxSV = 1e8;
end

%lambda
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INITIALISATION
y_hat     = zeros(size(y));

% points courrants
indxsup      = [];
gamma_sup    = [];
touslesind   = 1:size(x,2);
nbdim        = size(x,2);
% ind_nonxsup  = touslesind;
Beta         = zeros(1,nbdim);
arret         = 0;
move_backward = 0;


compt_trap     = 0;
compt_backward = 0;
compt_borne    = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INITIALISATION MEMOIRES
c_mem           = [];
Beta_mem       = [zeros(1,nbdim)];
cout_tot_mem   = [];
sum_Beta_mem   = [];

iter           = 0;
ptsadd         = [];
ptsret         = [];
scale_mem      = [];
stdR_mem       = [];
R_im1          = 100000;

if verbose ~=0
    gcf=figure(verbose);
    set(gcf, 'Position', [100 500  1000 400]);
    %keyboard
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BOUCLES

while arret == 0,
        
    ind_nonxsup          = setdiff(touslesind,indxsup);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CORRELATION VARIABLE RESIDU
     R                    = (y - y_hat);
    % R                       = max(0,1-(y.*y_hat));
    % R                       = R.*sign(y);
    
    %keyboard
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % VERIFICATION DE CVG 
%     if sum(abs(R)) > R_im1 + lambda %cvg non assure
%         fprintf('Cvg non assuree !\n');
%         break;
%     end
%     R_im1                = sum(abs(R));
    
    c                    = x'*R;
    
    [max_c,indmax_c]     = max(abs(c(ind_nonxsup))); % point de plus fort cout
    indmax_c             = ind_nonxsup(indmax_c);
    
    if move_backward == 0
        indxsup              = [indxsup; indmax_c];   
        ptsadd               = [ptsadd; iter indmax_c];
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  Existe-t-il une VA qui est plus correlee ?
    c_RV                 = abs(RVset*R); % std = 1, mean = 0
    Compt_RVpluscor      = length(find(c_RV > max_c));

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % MEMOIRE ET TRACE
    
    cout_tot_mem = [cout_tot_mem sum(R.^2)];
    sum_Beta_mem = [sum_Beta_mem sum(abs(Beta))];
    
    if verbose ~=0
        figure(verbose)  
        nbfig = 3;
        subplot(1,nbfig,1);
        
        h= plot(cout_tot_mem); 
        set(h,'linewidth',2);
        h= title('Cost (MSE)');
        set(h,'fontsize',16);
        
        subplot(1,nbfig,2);
        %figure(verbose+1)  
        h= plot(sum_Beta_mem,'r');
        set(h,'linewidth',2);
        h= title('\Sigma_i |\beta_i|');
        set(h,'fontsize',16);
        subplot(1,nbfig,3);
        %figure(verbose+2)  
        if length(ptsadd) > 0
            h= plot(ceil(ptsadd(:,2)/size(x,1)));
            set(h,'linewidth',2);
        end
        h= title('Scales of chosen points');
        set(h,'fontsize',16);
        
%         figure(verbose+1)
%         subplot(3,1,1)
%         plot(xappforplot,y,'b');
%         h= title('cos(exp(\omega t))');
%         set(h,'fontsize',16);
%         subplot(3,1,2)
%         plot(xappforplot,y_hat,'b');
%         h= title('\hat{y}');
%         set(h,'fontsize',16);