mulasso.m

lasso算法是stanford的tibshirani提出的一种变量子集选择算法。

function beta = MULASSO(Y,X,lambda,varargin)

% The LASSO based on Multiplicative updates, i.e. solving the problem
%
% beta = argmin{ min 0.5*||Y-beta X||_2^2+lambda*|beta|_1 }
%
% The algorithm is described in: 
% Morten M鴕up and Line Harder Clemmensen "Multiplicative updates for the
% LASSO" submitted to MLSP2007. Please refer to this article when
% publishing results obtained from the algorithm.
%
%
% Usage: 
%   beta=MU_LASSO(Y,X,lambda,opts)
%
% Input:
%   Y                   I x N matrix of observations
%   X                   M x N matrix of measured variables
%   lambda              L1 regularization penalty
%   opts.
%           beta        I x M matrix of initial guess on beta
%           maxiter     maximum number of iterations (default 15000)
%           alpha       Stepsize of multiplicative update can be either set to 0.5 or 1 (default 0.5)
%           conv_crit   convergence criterion (default 1e-8)
%
% Output:
%   beta            I x M esimated regressor
%
% Copyright (c) 2007  Morten M?rup, Line Harder Clemmensen and the Technical
% University of Denmark 


if nargin>=4, opts = varargin{1}; else opts = struct; end
conv_crit=mgetopt(opts,'conv_crit',1e-8);
maxiter=mgetopt(opts,'maxiter',30000);
alpha=mgetopt(opts,'alpha',0.5);
if alpha~=0.5 & alpha~=1
   disp('Illegal stepsize, stepsize (alpha) reset to 0.5');
   alpha=0.5;
end
%random initialisation
betap = max(abs(Y(:)))*(0.5+0.5*rand(size(Y,1),size(X,1)));
betan=betap+1e-2;
beta=betap-betan;
disp([' '])
disp(['Sparse Least Squares Regression using L1 regularisation'])
disp([' '])
disp(['To stop algorithm press control C'])
disp([' ']);
dheader = sprintf('%12s | %12s | %12s | %12s ','Iteration','dbeta/|beta|',' Time(s)   ','Cost','dcost','Var expl');
dline = sprintf('-------------+--------------+--------------+--------------+--------------+--------------+');

cost=inf;
SST=norm(Y,'fro')^2;
iter=0;
t1=cputime;    
actp=1:prod(size(betap));
actn=actp;

YYt = Y*Y';
YXt=Y*X';
YXtp=zeros(size(YXt));
YXtn=YXtp;
pind=find(YXt>0);
nind=find(YXt<0);
YXtp(pind)=YXt(pind);
YXtn(nind)=-YXt(nind);

XXt=X*X';
XXtp=zeros(size(XXt));
XXtn=XXtp;
pind=find(XXt>0);
nind=find(XXt<0);
XXtp(pind)=XXt(pind);
XXtn(nind)=-XXt(nind);
dbeta=inf;

while  iter<maxiter & dbeta>conv_crit
    if mod(iter,500)==0
         disp(dline); disp(dheader); disp(dline);
    end
    
    iter=iter+1;
    % find active set
    actoldp=actp;
    actoldn=actn;
    betaoldp=betap;
    betaoldn=betan;
    betaold=beta;
    actpt =find(betap(actp)/(norm(betap(actp),'fro')+eps)>1e-6);
    actp=actp(actpt);
    actnt =find(betan(actn)/(norm(betan(actn),'fro')+eps)>1e-6);
    actn=actn(actnt);
    remp=setdiff(actoldp,actp);
    betap(remp)=0;
    remn=setdiff(actoldn,actn);
    betan(remn)=0;    
    act=union(actp,actn);   
    bpXXtn(act)=betap(act)*XXtn(act,act);
    bpXXtp(act)=betap(act)*XXtp(act,act);
    bnXXtn(act)=betan(act)*XXtn(act,act);
    bnXXtp(act)=betan(act)*XXtp(act,act);
    
    gradp=((YXtp(actp)+bpXXtn(actp)+bnXXtp(actp)+eps)./(YXtn(actp)+bpXXtp(actp)+bnXXtn(actp)+lambda+eps));
    gradn=((YXtn(actn)+bpXXtp(actn)+bnXXtn(actn)+eps)./(YXtp(actn)+bpXXtn(actn)+bnXXtp(actn)+lambda+eps));
    if alpha==0.5
        betap(actp) = betap(actp).*sqrt(gradp);
        betan(actn) = betan(actn).*sqrt(gradn);
    else
        betap(actp) = betap(actp).*gradp;
        betan(actn) = betan(actn).*gradn;
    end

    
    beta(act)=betap(act)-betan(act);
    dbeta=beta(act)-betaold(act);
    dbeta=dbeta/(norm(betaold,'fro')+eps);
    tf=abs(dbeta)<1e-6;
    dbeta=norm(dbeta,'fro');

%    if rem(iter,50)==0   
%        told=t1;
%        t1=cputime;
%        cost_old=cost;
%        SSE=SST -2*trace(YXt*beta')+trace(beta*XXt*beta');
%        cost=0.5*SSE+lambda*sum(betap+betan);
%        dcost=cost_old-cost;
%        disp(sprintf('%12.0f | %12.4f | %12.4f | %12.4f | %12.4f| %12.4f',iter,dbeta,t1-told,cost,dcost,((SST-SSE)/SST)));
%   end   
    gzero=beta(act)>0;
    ind=find(gzero & tf);
    betap(act(ind))=beta(act(ind));
    betan(act(ind))=0;
    ind=find(~gzero & tf);
    betan(act(ind))=-beta(act(ind));
    betap(act(ind))=0;    
end
beta=betap-betan;

% -------------------------------------------------------------------------
% Parser for optional arguments
function var = mgetopt(opts, varname, default, varargin)
if isfield(opts, varname)
    var = getfield(opts, varname); 
else
    var = default;
end
for narg = 1:2:nargin-4
    cmd = varargin{narg};
    arg = varargin{narg+1};
    switch cmd
        case 'instrset',
            if ~any(strcmp(arg, var))
                fprintf(['Wrong argument %s = ''%s'' - ', ...
                    'Using default : %s = ''%s''\n'], ...
                    varname, var, varname, default);
                var = default;
            end
        otherwise,
            error('Wrong option: %s.', cmd);
    end
end