grabcode("http://web.stanford.edu/~boyd/papers/admm/basis_pursuit/basis_pursuit.html")

grabcode("http://web.stanford.edu/~boyd/papers/admm/basis_pursuit/basis_pursuit_example.html")

function [z, history] = basis_pursuit(A, b, rho, alpha)
% basis_pursuit  Solve basis pursuit via ADMM
%
% [x, history] = basis_pursuit(A, b, rho, alpha)
% 
% Solves the following problem via ADMM:
% 
%   minimize     ||x||_1
%   subject to   Ax = b
%
% The solution is returned in the vector x.
%
% history is a structure that contains the objective value, the primal and 
% dual residual norms, and the tolerances for the primal and dual residual 
% norms at each iteration.
% 
% rho is the augmented Lagrangian parameter. 
%
% alpha is the over-relaxation parameter (typical values for alpha are 
% between 1.0 and 1.8).
%
% More information can be found in the paper linked at:
% http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html
%

t_start = tic;

%% Global constants and defaults

QUIET    = 0;
MAX_ITER = 1000;
ABSTOL   = 1e-4; % 绝对误差容限
RELTOL   = 1e-2; % 相对误差容限

%% Data preprocessing

[~, n] = size(A);

%% ADMM solver

x = zeros(n,1);
z = zeros(n,1);
u = zeros(n,1);

if ~QUIET
    fprintf(‘%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n‘, ‘iter‘, ...
      ‘r norm‘, ‘eps pri‘, ‘s norm‘, ‘eps dual‘, ‘objective‘);
end

% precompute static variables for x-update (projection on to Ax=b)
AAt = A*A‘;
P = eye(n) - A‘ * (AAt \ A);
q = A‘ * (AAt \ b);

for k = 1:MAX_ITER
    % x-update
    x = P*(z - u) + q;

    % z-update with relaxation
    zold = z;
    x_hat = alpha*x + (1 - alpha)*zold; % 带松弛的更新方法
    z = shrinkage(x_hat + u, 1/rho);

    u = u + (x_hat - z);

    % diagnostics, reporting, termination checks
    history.objval(k)  = objective(A, b, x);

    history.r_norm(k)  = norm(x - z); % primal residual
    history.s_norm(k)  = norm(-rho*(z - zold)); % dual residual
    
    history.eps_pri(k) = sqrt(n)*ABSTOL + RELTOL*max(norm(x), norm(-z)); % primal feasibility tolerance
    history.eps_dual(k)= sqrt(n)*ABSTOL + RELTOL*norm(rho*u); %  dual feasibility tolerance

    if ~QUIET
        fprintf(‘%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n‘, k, ...
            history.r_norm(k), history.eps_pri(k), ...
            history.s_norm(k), history.eps_dual(k), history.objval(k));
    end

    if (history.r_norm(k) < history.eps_pri(k) && ...
       history.s_norm(k) < history.eps_dual(k))
         break;
    end
end

if ~QUIET
    toc(t_start);
end

end

function obj = objective(A, b, x)
    obj = norm(x,1);
end

function y = shrinkage(a, kappa)
    y = max(0, a-kappa) - max(0, -a-kappa);
end

% Basis pursuit example

%% Generate problem data
rand(‘seed‘, 0);
randn(‘seed‘, 0);

n = 30; 
m = 10;
A = randn(m,n);

x = sprandn(n, 1, 0.1*n);
b = A*x;

xtrue = x;

%% Solve problem

[x, history] = basis_pursuit(A, b, 1.0, 1.0);

%% Reporting
K = length(history.objval);                                                                                                        

h = figure;
plot(1:K, history.objval, ‘k‘, ‘MarkerSize‘, 10, ‘LineWidth‘, 2); 
ylabel(‘f(x^k) + g(z^k)‘); xlabel(‘iter (k)‘);

g = figure;
subplot(2,1,1);                                                                                                                    
semilogy(1:K, max(1e-8, history.r_norm), ‘k‘, ...
    1:K, history.eps_pri, ‘k--‘,  ‘LineWidth‘, 2); % 原始问题的残差
ylabel(‘||r||_2‘); 

subplot(2,1,2);                                                                                                                    
semilogy(1:K, max(1e-8, history.s_norm), ‘k‘, ...
    1:K, history.eps_dual, ‘k--‘, ‘LineWidth‘, 2);   % 对偶问题的n

ylabel(‘||s||_2‘); xlabel(‘iter (k)‘);