Contents

% Comparison of different linearized Bregman approaches on recovering a sparse solution x from Ax = b
% The nonzero entries of x have iid Gaussian random values
% A is a partial DCT matrix

clear

Generate problem data

rand('seed', 2012); randn('seed', 2012);

m = 250; n = 500; % matrix dimension m-by-n
k = 25; % sparsity

A = dctmtx(n); % dct matrix
A = A([1 randperm(n,m-1)+1],:); % randomly select m rows but always include row 1
x_ref = zeros(n,1); % true vector
x_ref(randsample(n,k)) = randn(k,1); % Gaussian random values
b = A*x_ref; % finish generating equations Ax = b

set parameters

alpha = 5*norm(x_ref,inf); % don't need to be exact, roughly 1 - 10 times norm(x_ref,inf) is fine
opts.tol = 1e-6;   % stop once norm(Ax-b)<tol*norm(b)
opts.maxit = 2000; % run maximally 1000 iterations
opts.x_ref = x_ref;

LBreg: fixed stepsize ---

opts.stepsize = 2/alpha/normest(A*A.',1e-2); % roughly 2/alpha/norm(A)^2
t0 = tic;
[x,out] = lbreg_fixedstep(A,b,alpha,opts);
time = toc(t0);
fprintf('iter = %d, time = %4.2e, ', out.iter, time);
fprintf('solution relative error = %4.2e\n', norm(x - x_ref)/norm(x_ref));
opts = rmfield(opts, 'stepsize');

iter = 268, time = 3.37e-02, solution relative error = 1.29e-06

LBreg: Barzilai-Borwein and non-montone line search ---

opts.stepsize = 2/alpha/normest(A*A.',1e-2); % roughly 2/alpha/norm(A)^2
t0 = tic;
[x1,out1] = lbreg_bbls(A,b,alpha,opts);
time = toc(t0);
fprintf('iter = %d, time = %4.2e, ', out1.iter, time);
fprintf('solution relative error = %4.2e\n', norm(x1 - x_ref)/norm(x_ref));
opts = rmfield(opts, 'stepsize');

iter = 63, time = 7.99e-03, solution relative error = 1.14e-06

LBreg: accelerated, no restart ---

t0 = tic;
[x2,out2] = lbreg_accelerated(A,b,alpha,opts);
time = toc(t0);
fprintf('iter = %d, time = %4.2e, ', out2.iter, time);
fprintf('solution relative error = %4.2e\n', norm(x2 - x_ref)/norm(x_ref));

iter = 108, time = 1.43e-02, solution relative error = 1.13e-06

LBreg: skip: monotonic gradient scheme ---

opts.reset = 6;
t0 = tic;
[x3,out3] = lbreg_accel_w_reset(A,b,alpha,opts); % no restart
time = toc(t0);
fprintf('iter = %d, time = %4.2e, ', out3.iter, time);
fprintf('solution relative error = %4.2e\n', norm(x3 - x_ref)/norm(x_ref));
opts = rmfield(opts, 'reset');

iter = 86, time = 1.08e-02, solution relative error = 7.84e-07

Reporting

figure;
plot(1:out.iter,  out.hist_obj,  'k-', ...
     1:out1.iter, out1.hist_obj, 'b--', ...
     1:out2.iter, out2.hist_obj, 'g--', ...
     1:out3.iter, out3.hist_obj, 'r--', ...
     'LineWidth', 2);
legend('fixed stepsize', 'BB+line search', 'Nesterov accel', 'Nesterov+skip','Location','SouthEast');
title('Dual objective')
xlabel('iteration'); ylabel('dual objective');

figure;
semilogy(1:out.iter,  out.hist_err/norm(x_ref),  'k-', ...
         1:out1.iter, out1.hist_err/norm(x_ref), 'b--', ...
         1:out2.iter, out2.hist_err/norm(x_ref), 'g--', ...
         1:out3.iter, out3.hist_err/norm(x_ref), 'r--', ...
         'LineWidth', 2);
legend('fixed stepsize', 'BB+line search', 'Nesterov accel', 'Nesterov+skip','Location','NorthEast');
title('Primal solution relative error')
xlabel('iteration'); ylabel('||x - x_{ref}||_2/||x_{ref}||_2');

Download this m-file

Published with MATLAB® R2012b