套索识别算法

function [history, stopReason,k] = lars(yin, xin, XTX, type, stopCriterion, regularization, trace, quiet) % %{ % This program implements "LARS" algorithm and its lasso modification % introduced by Efron et. al. 2003. Read the paper to understand codes % of this function. Each line of this file has corresponding equation % number in Efron et. al. 2003 for reader's convenience. % % % *** CAUTION % history(1).mu_OLS contains original 'yin' to provide convinience in % writing user defined stop criterion function. Actual mu_OLS of the first % step should be just mean of yin which is a simple array of copy of % history(1).mu. So, if history(1).mu_OLS contains that information, it is % redundant. In fact, this is also contained in history(1).b, which is a % bias of the output. Therefore, to provide more information to user who % want to write his/her own stop criterion function, history(1).mu_OLS % contains yin. % % % Example 1: moderate size x. % stopCrioterion = {}; % stopCrioterion{1,1} = 'maxKernels'; % stopCrioterion{1,2} = 100; % % XTX = lars_getXTX(x_original); % this takes long time. % sol = lars(y, x, XTX,'lasso', stopCriterion); % % Example 2: very small size x, or a really really big size x % stopCrioterion = {}; % stopCrioterion{1,1} = 'maxKernels'; % stopCrioterion{1,2} = 100; % % sol = lars(y, x,XTX, 'lasso', stopCriterion); % % Note: % Users can add any kind of stop criterion by editing % the corresponding portion of this file. See the code % for existing examples. % % Note 2: % This m-file does not implement routine for missing data. % %} % global USING_CLUSTER; % global RESOLUTION_OF_LARS; global REGULARIZATION_FACTOR; RESOLUTION_OF_LARS=0.0001; lars_init(); regularization_factor = REGULARIZATION_FACTOR; % Tikhonov regularization factor (or the ridge regression factor) % -> This should be small enough in this case to get % a reasonable pseudoinverse. stopReason = {}; %%% Check parameters if isempty(yin) || isempty(xin) warning('\nInput or Output has zero length.\n'); history.active_set = []; stopReason{1} = 'Parameter error'; stopReason{2} = 0; return; end if size(yin,1) ~= size(xin,1) warning('\nSize of y does not match to that of x.\n'); history.active_set = []; stopReason{1} = 'Parameter error'; stopReason{2} = 0; return; end if ~strcmp(type, 'lasso') & ~strcmp(type, 'lars') & ~strcmp(type, 'forward_stepwise') warning('\nUnknown type of regression.\n'); history.active_set = []; stopReason{1} = 'Parameter error'; stopReason{2} = 0; return; end if strcmp(type, 'forward_stepwise') warning('\nForward_stepwise is not implemented.\n'); history.active_set = []; stopReason{1} = 'Parameter error'; stopReason{2} = 0; return; end if exist('regularization','var') && ~isempty(regularization) regularization = 10; else regularization = 0; end if ~exist('trace','var') || isempty(trace) trace=0; end if ~exist('quiet','var') || isempty(quiet) quiet=0; elseif quiet==1 trace=0; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Data preparation % Program automatically centers and standardizes predictors. if ~exist('XTX','var') XTX=[]; end no_xtx = 0; if ~isempty(XTX) if ~quiet && trace >=0 fprintf('\nLars is using the provided xtx.\n'); end elseif size(xin,2)^2 > 10^6 if ~quiet && trace >=0 fprintf('Too large matrix (size(x,2)^2 > 10^6). lars will not pre-calculate xtx.\n'); end no_xtx = 1; XTX = lars_getXTX(xin,no_xtx); else % fprintf('\nCalculating xtx.\n'); XTX = lars_getXTX(xin); end x = XTX.x; % normalized xin mx = XTX.mx; % mean xin sx = XTX.sx; % length of each column of xin ignores = XTX.ignores; % indices for constant terms all_candidate = XTX.all_candidate;% indices for all possible columns if ~no_xtx xtx = XTX.xtx; % xtx matrix dup_columns = XTX.dup_columns; % duplicated columns which will be automatically ignored. end my = mean(yin); y = yin-my; n = size(x,1); % # of samples m = size(x,2); % # of predictors % Now, we can determine the maximum number of kernels. %maxKernels = min(maxKernels, min(size(xin,1)-1, length(all_candidate))); %maxKernels = min(maxKernels, min(rank(xin), length(all_candidate))); existMaxKernels = 0; existMSE = 0; for is = 1:size(stopCriterion,1) if strcmp(stopCriterion{is,1},'maxKernels') existMaxKernels = 1; stopCriterion{is,2} = min(stopCriterion{is,2}, min(size(xin,1)-1, length(all_candidate))); % stopCriterion{is,2} = min(stopCriterion{is,2}, min(rank(xin), length(all_candidate))); if stopCriterion{is,2}<1 warning('Max Kernel is less than 1. It must be larger than 0.\n'); stopCriterion{is,2} = 1; end end if strcmp(stopCriterion{is,1},'MSE') existMSE = 1; if stopCriterion{is,2}<1.0e-10 warning('Maximum MSE is too small. Automatically set to 1.0e-10\n'); stopCriterion{is,2} = 1.0e-10; end end end if ~existMaxKernels is = size(stopCriterion,1); stopCriterion{is+1,1} = 'maxKernels'; % stopCriterion{is+1,2} = min(size(xin,1)-1, length(all_candidate)); % Stop when size of active set is data.maxKernels. stopCriterion{is+1,2} = min(rank(xin), length(all_candidate)); % Stop when size of active set is data.maxKernels. end if ~existMSE is = size(stopCriterion,1); stopCriterion{is+1,1} = 'MSE'; stopCriterion{is+1,2} = 1.0e-10; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialization active = []; % active set inactive = all_candidate; % inactive set mu_a = zeros(n,1); % current estimate (eq. 2.8) mu_a_plus = 0; % next estimate (eq. 2.12) mu_a_OLS = 0; % OLS estimate (eq. 2.19) beta = zeros(1,size(x,2)); beta_new = beta; beta_OLS = beta; history.active_set = []; history.add = []; history.drop = []; history.beta_norm = []; history.beta = []; history.b = my; history.mu = my; history.beta_OLS_norm = []; history.beta_OLS = []; history.b_OLS = my; history.mu_OLS = my*ones(size(yin)); history.MSE = sum(y.^2)/length(y); history.R_square = 0; history.resolution_warning = []; if var(yin)==0 stopReason{1} = 'zeroVarY'; stopReason{2} = var(yin); return; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Main loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c = 0; % correlation vector C_max = max(abs(c)); C_max_ind = []; C_max_ind_pl = []; drop = []; % used for 'lasso' k = 1; % iteration index %if ~quiet && trace >= 0 % fprintf('Active predictors / total ::::: Current iteration\n '); %end while 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%% Exit Criterions %%%%%%%%%%%%%%%%%%