perform_homotopy.m~

signal procesing toolbox

function [X,lambda_list,sparsity_list] = perform_homotopy(D,y,options)

% perform_homotopy - compute the L1 regularization path
%
%   X = perform_homotopy(D,y,options);
%
%   Copyright (c) 2008 Gabriel Peyre

n = size(D,1);
p = size(D,2);

niter = getoptions(options, 'niter', min(,n/2);

X = zeros(p,niter);
xbp = zeros(p,1);
niter = round(1.5*M);
sparsity_list = [];
lambda_list = [];
for i=1:niter
    % correlation
    C = D'*(y-D*xbp);
    lambda = max(abs(C));
    % support for update
    S = find( abs(abs(C/lambda)-1)<1e-9);
    I = ones(p,1); I(S)=0;
    Sc = find(I);
    % update direction
    d = zeros(p,1);
    d(S) = (D(:,S)'*D(:,S)) \ sign( C(S) );
    % useful vector
    v = D(:,S)*d(S);
    % Compute minimum |gamma| so that situation 1) is in force.
    w = ( lambda-C(Sc) ) ./ ( 1 - D(:,Sc)'*v );
    gamma1 = min(w(w>0));
    % Compute minimum |gamma| so that situation 2) is in force.
    w = ( lambda+C(Sc) ) ./ ( 1 + D(:,Sc)'*v );
    gamma2 = min(w(w>0));
    % Compute minimum |gamma| so that situation 3) is in force.
    w = -xbp(S)./d(S);
    gamma3 = min(w(w>0));
    % any condition is in force
    gamma = min([gamma1 gamma2 gamma3]);
    % new solution
    xbp = xbp + gamma*d;
    % record sparsity and lambda
    sparsity_list(i) = sum( abs(xbp)>1e-9 );
    lambda_list(i) = lambda;
end
% Display the evolution of the sparsity as a function of |lambda|
clf;
plot(lambda_list,sparsity_list, '.-'); axis('tight');
set_graphic_sizes([], 20);
set_label('lambda', 'sparsity');