lar.m

基于核分析的多类分类器

%         subplot(3,1,3)
%         plot(xappforplot,R,'b');
%         h= title('Residu');
%         set(h,'fontsize',16);
% 
%         figure(verbose+2)
%         imagesc(reshape(c(1:end-1,:),nbptsapp,nbvar/nbptsapp));
%         h= title('Correlation residu/K_i(x_j,\cdot)');
%         set(h,'fontsize',16);
        
        drawnow;
%         pause
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CALCUL DU PAS
    
    S        = sign(c);   
    
    % METHODE CHOL
    
        Xsup     = x(:,indxsup) * diag(S(indxsup));
        

            regterm    = lambda;
    
        Gsup     = Xsup'*Xsup + eye(length(indxsup))*regterm;           
        
%   CHOL        
%         Rchol    = chol(Gsup);
%         bchol    = Rchol'\ones(size(Gsup,1),1);
%         ychol    = Rchol\bchol;
%         Asup     = (sum(ychol))^(-0.5);
%         Wsup     = Asup*ychol;
    %     
% INVERSION BRUTALE !
    
    invGsup  = inv(Gsup);    % INVERSION MAT !
    Asup     = (sum(sum(invGsup)))^(-0.5);
    Wsup     = sum(Asup*invGsup,2);
    
% METHODE SANS INVERSION -> uniquement pour les bases orthogonales
% 
%     Wsup       = [R'*Xsup]'; % projection du residu sur les variables actives
%     %Usup       = Xsup*Wsup;
%     %Usup       = Usup/sqrt(Usup'*Usup); % normalisation    
%     b          = x'*x*W;

    % FIN COMMUNE
%     
    Usup     = Xsup*Wsup;
    b        = x'*Usup;
    W          = zeros(nbdim,1);
    W(indxsup) = Wsup;
        
    ind_nonxsup          = setdiff(touslesind,indxsup);
    
    
    if length(ind_nonxsup) ~= 0
        gamma_out_xsup_moins = (max_c - c(ind_nonxsup))./(Asup - b(ind_nonxsup));
        gamma_out_xsup_plus  = (max_c + c(ind_nonxsup))./(Asup + b(ind_nonxsup));
        all_gamma            = [gamma_out_xsup_moins; gamma_out_xsup_plus];
        all_gamma            = all_gamma(find(all_gamma>0)); % uniquement les valeurs positives des gammas
        gamma                = min(all_gamma);
        if isempty(gamma)
            gamma = eps;
        end
    else
        gamma                = max_c/Asup;
    end
    
    move_backward = 0;
    % calcul de beta

    Beta_tp1             = Beta + gamma*S'.*W';  
    
    % gamma_OLS            = max_c/Asup;
     gamma_OLS            = max(c(indxsup)./b(indxsup));
    % gamma_OLS            = (gamma + 3*min(c(indxsup)./b(indxsup)))/4;
    if isempty(gamma_OLS)
        gamma_OLS = eps;
    end
    Beta_OLS             = Beta + gamma_OLS*S'.*W';
    
    inddiffsign = [];
    inddiff0             = find(Beta > eps | Beta < -eps);
    inddiffsign          = find(sign(Beta(inddiff0)) ~= sign(Beta_tp1(inddiff0)));
    inddiffsign          = inddiff0(inddiffsign);
    
    if length(inddiffsign) == 0
        Beta          = Beta_tp1;
    elseif length(inddiffsign) == 1
        if verbose ~= 0
            fprintf('move backward (%d) \n', length(indxsup));
        end
        
        indtodel      = inddiffsign(1);        
        gamma         = Beta(indtodel)/(W(indtodel)'*-S(indtodel));
        
        Beta          = Beta + gamma*S'.*W';        
        indtodel      = find(indxsup == indtodel);
        move_backward = 1;
        compt_backward = compt_backward+1;
    else     
        if verbose ~= 0
            fprintf('move backward mult (%d) \n', length(indxsup));
        end
        
        gamma         = Beta(inddiffsign)./(W(inddiffsign)'.*-S(inddiffsign)');
        [gamma,ind]   = min(gamma);
        indtodel      = inddiffsign(ind);
        
        Beta          = Beta + gamma*S'.*W';       
        indtodel      = find(indxsup == indtodel);
        move_backward = 1;
        compt_backward = compt_backward+1;
    end
    
    
    % stockage gamma
    gamma_sup            = [gamma_sup gamma];
    
    % calcul de y_hat
    y_hat                = x*Beta';
    
    y_hat_OLS            = x*Beta_OLS';
    
    Beta_mem   = [Beta_mem; Beta];
    sumabsBeta = sum(abs(Beta_mem),2);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % LASSO (SUITE)
    
    % enlever le vecteur support correspondant
    if move_backward == 1     
        ptsret                 = [ptsret; iter indxsup(indtodel)];
        indxsup(indtodel)      = [];
        gamma_sup(indtodel)    = [];
        sum_Beta_mem(indtodel) = [];
        cout_tot_mem(indtodel) = [];
        ind_nonxsup            = setdiff(touslesind,indxsup);
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CRITERES ARRET
    
    if length(indxsup) > Limites.nbmaxSV | compt_backward >Limites.nbmaxMoveBack 
        % il faut tout arreter !
        fprintf('Bornes de securite depassees.\n');
        break;
    end
    
    if move_backward ~= 1  % ne pas prendre de solution sous-optimale
        for i=1:nbborne
            RETENIR = 0;    
            if strcmp(AllBorne{i}.type,'sumBeta' ) & length(AllBorne{i}.borne)>0
                if sum(abs(Beta)) > AllBorne{i}.borne(1)  
                    RETENIR = 1;
                end
                
            elseif strcmp(AllBorne{i}.type,'pcSV' ) & length(AllBorne{i}.borne)>0
                if length(indxsup)/nbdim > AllBorne{i}.borne(1)
                    RETENIR = 1;
                end
                
            elseif strcmp(AllBorne{i}.type,'nbSV') & length(AllBorne{i}.borne)>0
                if length(indxsup) >= AllBorne{i}.borne(1)
                    RETENIR = 1;            
                end
                
            elseif strcmp(AllBorne{i}.type,'trapscale')  & length(AllBorne{i}.borne)>0
                if (indxsup(end) > AllBorne{i}.indTS & indxsup(end) < size(x,2)) % le point appartient a une echelle interdite
                    compt_trap = length(find(indxsup > AllBorne{i}.indTS & indxsup < size(x,2))); % combien ?
                    
                    if compt_trap > AllBorne{i}.borne(1) 
                        RETENIR = 1;    
                    end
                end
                
            elseif strcmp(AllBorne{i}.type,'RVtrap') & length(AllBorne{i}.borne)>0
                if Compt_RVpluscor > AllBorne{i}.borne(1)
                    RETENIR = 1;   
                end
                
            elseif strcmp(AllBorne{i}.type,'countback') & length(AllBorne{i}.borne)>0
                if compt_backward > AllBorne{i}.borne(1) 
                    RETENIR = 1;    
                end    
                
            end
            
            
            if RETENIR == 1
                solution{nbexpcur}.indxsup = indxsup;
                solution{nbexpcur}.Beta    = Beta(indxsup);
                solution{nbexpcur}.type    = AllBorne{i}.type;
                if strcmp(AllBorne{i}.type,'trapscale'),
                    solution{nbexpcur}.borne   = AllBorne{i}.borne(1)*0.01 + AllBorne{i}.indTS;% recomposition !
                else
                    solution{nbexpcur}.borne   = AllBorne{i}.borne(1);
                end
                
                if test_OLS == 1
                    solution_OLS{nbexpcur}.indxsup = indxsup;
                    solution_OLS{nbexpcur}.Beta    = Beta_OLS(indxsup);
                    solution_OLS{nbexpcur}.type    = AllBorne{i}.type;% '_OLS'];
                    if strcmp(AllBorne{i}.type,'trapscale'),
                        solution_OLS{nbexpcur}.borne   = AllBorne{i}.borne(1)*0.01 + AllBorne{i}.indTS; % recomposition !
                    else
                        solution_OLS{nbexpcur}.borne   = AllBorne{i}.borne(1);
                    end
                        
                end
                
                
                %fprintf('%4d/%4d : %15s (%f) # %d\n',nbexpcur,nbexptot,solution{nbexpcur}.type,solution{nbexpcur}.borne,length(indxsup));
                fprintf('\b\b\b\b\b\b\b\b\b%4d/%4d',nbexpcur,nbexptot);
                
                AllBorne{i}.borne(1)     = []; % elimination !
                nbexpcur                 = nbexpcur +1;     
                %            size(indxsup)
            end
        end
    end
    if length(indxsup) == nbdim 
        fprintf('Regression totale !\n');
        break;
    end
    
    
    
    iter = iter+1;
    
    if nbexpcur-1==nbexptot, % le compteur est en avance !
        break;
    end
    % pause
end
fprintf('\n');

nbexpcur    = nbexpcur-1;

for i=1:nbexpcur
    
    indb = find(solution{i}.indxsup == size(x,2));
    if length(indb) ~= 0
        solution{i}.b = solution{i}.Beta(indb) / sqrt(nbptsapp); % "denormalisation"
        solution{i}.Beta(indb) = [];
        solution{i}.indxsup(indb) = [];    
    else
        solution{i}.b = 0;
    end
    
    if test_OLS == 1
        if length(indb) ~= 0
            solution_OLS{i}.b = solution_OLS{i}.Beta(indb) / sqrt(nbptsapp); % "denormalisation"
            solution_OLS{i}.Beta(indb) = [];
            solution_OLS{i}.indxsup(indb) = [];    
        else
            solution_OLS{i}.b = 0;
        end
    end
end