## 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', 0); randn('seed', 0);

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 = 330, time = 4.12e-02, solution relative error = 7.79e-07
```

## 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 = 86, time = 1.06e-02, solution relative error = 2.83e-07
```

## 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 = 122, time = 1.46e-02, solution relative error = 9.81e-07
```

## 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 = 118, time = 1.50e-02, solution relative error = 2.60e-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');
```