lasso_box.m

function [z, history] = lasso_box(A, b, lambda, rho, alpha)
% lasso  Solve lasso problem via ADMM
%
% [z, history] = lasso_box(A, b, lambda, rho, alpha);
% 
% Solves the following problem via ADMM:
%
%   minimize 1/2*|| Ax - b ||_2^2 + \lambda || x ||_1
%	subject to: x belong to [0,1]
%
% 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).
%

t_start = tic;

%% Global constants and defaults

QUIET    = 0;
MAX_ITER = 1000;
ABSTOL   = 1e-10;
RELTOL   = 1e-10;

%% Data preprocessing

[m, n] = size(A);

% save a matrix-vector multiply
Atb = A'*b;

%% ADMM solver

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

% cache the factorization
[L U] = factor(A, rho);

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

for k = 1:MAX_ITER
    
    % x-update
%     q = Atb + rho*(z - u);    % temporary value
    q = Atb + rho*(0.5 * (z - u) +  0.5 * (z1 - u1));
    if( m >= n )    % if skinny
       x = U \ (L \ q);
    else            % if fat
       x = q/rho - (A'*(U \ ( L \ (A*q) )))/rho^2;
    end

    % z-update with relaxation
    zold = z;
    x_hat = alpha*x + (1 - alpha)*zold;
    z = shrinkage(x_hat + u, lambda/rho);
    
    zold1 = z1;
    z1 = box0_1(x_hat+u1);    
    
    % u-update
    u = u + (x_hat - z);
    u1= u1 + (x_hat - z1);
    
    % diagnostics, reporting, termination checks
    history.objval(k)  = objective(A, b, lambda, x, z);

    history.r_norm(k)  = norm(x - z);
    history.s_norm(k)  = norm(-rho*(z - zold));
    history.r1_norm(k) = norm(x - z1);
    history.s1_norm(k) = norm(-rho*(z1 - zold1));
    
    history.eps_pri(k) = sqrt(n)*ABSTOL + RELTOL*max(norm(x), norm(-z));
    history.eps_dual(k)= sqrt(n)*ABSTOL + RELTOL*norm(rho*u);
    history.eps_pri1(k)= sqrt(n)*ABSTOL + RELTOL*max(norm(x), norm(-z1));
    history.eps_dual1(k)=sqrt(n)*ABSTOL + RELTOL*norm(rho*u1);
    
    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.r1_norm(k),history.eps_pri1(k), ...
            history.s1_norm(k),history.eps_dual1(k), ...
            history.objval(k));
    end

    if (history.r_norm(k) < history.eps_pri(k) && ...
       history.s_norm(k) < history.eps_dual(k)&& ...
       history.r1_norm(k) < history.eps_pri1(k)&& ...
       history.s1_norm(k) < history.eps_dual1(k))
         break;
    end

end

if ~QUIET
    toc(t_start);
end

end

function p = objective(A, b, lambda, x, z)
    p = ( 1/2*sum((A*x - b).^2) + lambda*norm(z,1) );
end

function z = shrinkage(x, kappa)
    z = max( 0, x - kappa ) - max( 0, -x - kappa );
end

function [L U] = factor(A, rho)
    [m, n] = size(A);
    if ( m >= n )    % if skinny
       L = chol( A'*A + rho*speye(n), 'lower' );
    else            % if fat
       L = chol( speye(m) + 1/rho*(A*A'), 'lower' );
    end
    
    % force matlab to recognize the upper / lower triangular structure
    L = sparse(L);
    U = sparse(L');
end

function z_hat = box0_1(v)
    z_hat = max(0,v);
    z_hat = min(z_hat,1);
end