function [L,S,fvals,resids] = pcp_admm(M,lambda,rho)
%
% [L,S] = pcp_admm(M)
% [L,S] = pcp_admm(M,lambda)
% [L,S] = pcp_admm(M,lambda,rho)
% 
% Computes a solution to the Principal Component Pursuit (PCP) problem
%   minimize:   ||L||_* + lambda*||S||_1
%   subject to: L + S == M
% via the Alternating Direction Method of Multipliers (ADMM) method
% 
% Original Paper:
% ---------------
% Emmanuel J. Candes, Xiaodong Li, Yi Ma, and John Wright. 
% Robust principal component analysis. Technical report,
% Stanford University, 2009.
%
% Input: 
%   M      = m-by-n data matrix
%   lambda = regularization parameter for original problem
%            default lambda = 1/sqrt(max(m,n))
%   rho    = starting penalty parameter for ADMM
%            default rho = 1
%            this script automatically varies the penalty parameter
%            and stops on a relative criterion
%
% Output:
%   L = m-by-n, low-rank
%   S = m-by-n, sparse

[m,n] = size(M);

% default parameter values
if nargin <= 2
    rho = 1;
end
if nargin <= 1
    lambda = 1/sqrt(max(m,n));
end

% algorithm parameters and initial point
MAX_ITER = 5000;
TOL = 1e-2;
L = ones(m,n);
S = M - L;
W = ones(m,n)/rho;

fvals = [];
resids = [];

rhonext = rho;

fprintf('Running ADMM...\n');
fprintf('Iteration |    Objective | normf(L+S-M) |   Dual resid |     Time (s)\n');

for k=1:MAX_ITER
    titerstart = tic;
    
    % 1. prox of the nuclear norm
    [U,Sig,V] = svd(M-S-W,0);
    Lnext = U * diag(soft_thresh(diag(Sig), 1/rho)) * V';
    
    % 2. prox of the l1 norm
    Snext = soft_thresh(M-Lnext-W, lambda/rho);
    
    % 3. residual running sum
    Wnext = W + (Lnext + Snext - M);
    
    % 4. evaluate primal and dual residuals
    % r^{k+1} = L^{k+1} + S^{k+1} - M
    % s^{k+1} = rho*(S^{k+1}-S^{k})
    primal_resid = norm(Lnext + Snext - M,'fro');
    dual_resid = rho*norm(Snext - S,'fro');
    
    % 5. automatically vary the penalty parameter
    if primal_resid > 10*dual_resid
        rhonext = 2*rho;
        Wnext = Wnext/2;
        fprintf('--------> Increasing penalty to rho=%e\n', rhonext);
    elseif dual_resid > 10*primal_resid
        rhonext = rho/2;
        Wnext = 2*Wnext;
        fprintf('--------> Decreasing penalty to rho=%e\n', rhonext);
    else
        rhonext = rho;
    end
    
    % 6. update iterates
    L = Lnext;
    S = Snext;
    W = Wnext;
    rho = rhonext;
    
    % stopping criterion
    if primal_resid <= TOL && dual_resid <= TOL
        fprintf('Converged to TOL=%e in %d iterations\n', TOL, k);
        break;
    end
    
    % save the objective and feasibility residuals
    if mod(k, 1) == 0
        fval = norm_nuc(L) + lambda*norm(S(:),1);
        fvals = [fvals, fval];
        resids = [resids, [primal_resid; dual_resid]];
    end
    
    % show progress
    titerstop = toc(titerstart);
    if mod(k, 1) == 0
        fprintf('%9d | %e | %e | %e | %12f\n', ...
            k, fval, primal_resid, dual_resid, titerstop);
    end
end

if k == MAX_ITER
    fprintf('Failed to converge in MAX_ITER=%d iterations\n', MAX_ITER);
end


end

function result = soft_thresh(X, kappa)
% computes the shrinkage S_kappa(X)
assert(kappa > 0, 'soft thresholding on negative kappa');
result = max(X - kappa,0) - max(-X - kappa,0);
end