Contents

% An example of recovery a sparse solution x from Ax = b using linearized Bregman iteration with Nesterov's acceleration and different reset schemes
% The nonzero entries of x have iid Gaussian random values

clear

Generate problem data

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

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

A = randn(m,n); % random matrix
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

Solve problem

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 = 1000; % run maximally 1000 iterations
opts.x_ref = x_ref;

% --- accelerated, no restart ---
opts.reset = 0;
t0 = tic;
[x,out] = lbreg_accel_w_reset(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));

% --- restart: monotonic dual objective ---
opts.reset = 1;
t0 = tic;
[x1,out1] = lbreg_accel_w_reset(A,b,alpha,opts); % no restart
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));

% --- skip: monotonic residual size ---
opts.reset = 2;
t0 = tic;
[x2,out2] = lbreg_accel_w_reset(A,b,alpha,opts); % no restart
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));

% --- restart: monotonic gradient scheme ---
opts.reset = 5;
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));

% --- skip: monotonic gradient scheme ---
opts.reset = 6;
t0 = tic;
[x4,out4] = lbreg_accel_w_reset(A,b,alpha,opts); % no restart
time = toc(t0);
fprintf('iter = %d, time = %4.2e, ', out4.iter, time);
fprintf('solution relative error = %4.2e\n', norm(x4 - x_ref)/norm(x_ref));

iter = 260, time = 2.57e-02, solution relative error = 1.42e-06
iter = 183, time = 2.58e-02, solution relative error = 1.04e-06
iter = 140, time = 2.12e-02, solution relative error = 5.61e-07
iter = 179, time = 2.46e-02, solution relative error = 9.19e-07
iter = 141, time = 1.89e-02, solution relative error = 8.43e-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:', ...
     1:out4.iter, out4.hist_obj, 'm--', ...
     'LineWidth', 2);
legend('no restart', 'restart: monotonic obj', 'skip: monotonic obj', 'restart: grad scheme', 'skip: grad scheme', '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:', ...
         1:out4.iter, out4.hist_err/norm(x_ref), 'm--', ...
         'LineWidth', 2);
legend('no restart', 'restart: monotonic obj', 'skip: monotonic obj', 'restart: grad scheme', 'skip: grad scheme', 'Location','SouthWest');
title('Primal solution relative error')
xlabel('iteration'); ylabel('||x - x_{ref}||_2/||x_{ref}||_2');

Download this m-file

Published with MATLAB® R2012b