mulassoqp.m

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

function beta = MULASSOqp(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_LASSOqp(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',20000);
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('-------------+--------------+--------------+--------------+--------------+--------------+');


iter=0;
t1=cputime;    
actp=1:prod(size(betap));
actn=actp;

YYt = Y*Y';


SST=trace(YYt)

YXt=Y*X';

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;
cost=inf;
act=1:length(beta);

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);   
   
    bpXXt=zeros(size(beta));
    bnXXt=bpXXt;
    bnXXt(act)=XXt(act,act)*betan(act)';
    betap=update_nqp(betap,XXtp,XXtn,-YXt-bnXXt+lambda,act);
    bpXXt(act)=XXt(act,act)*betap(act)';
    betan=update_nqp(betan,XXtp,XXtn,YXt-bpXXt+lambda,act);
    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;

% -------------------------------------------------------------------------
function v=update_nqp(v,Ap,An,b,act)

v(act)=v(act).*(-b(act)+sqrt(b(act).*b(act)+4*(v(act)*Ap(act,act)).*(v(act)*An(act,act))))./(2*v(act)*Ap(act,act)+eps);


% -------------------------------------------------------------------------
% 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