ADMM_SPCP.m

function [out] = ADMM_SPCP(A,options)
% minimize(0.5*square_pos(norm(E,'fro')) + lambda_S*norm(S(:),1) + lambda_L*norm_nuc(L))

% subject to
% E + S + L == A;
% based on the ADMM code in https://web.stanford.edu/~boyd/papers/prox_algs/matrix_decomp.html
%
% PARAMETRIZATION NOTE (relation to the SPCP papers and to the paper's Table 1)
% -----------------------------------------------------------------------------
% The literature (Zhou et al., "Stable Principal Component Pursuit", ISIT 2010)
% writes SPCP with the nuclear norm normalized to coefficient 1:
%
%       minimize ||L||_* + lambda*||S||_1 + (mu/2)*||E||_F^2   s.t. A = L+S+E,
%
% with the theoretically motivated choices lambda = 1/sqrt(max(m,n)) and
% mu = 1/(sqrt(2*max(m,n))*sigma), where sigma is the std of the dense noise.
% This file minimizes the equivalent scaled form (multiply the objective above
% by 1/mu = lambda_L):
%
%       minimize 0.5*||E||_F^2 + lambda_S*||S||_1 + lambda_L*||L||_*,
%
% so the two parametrizations map to each other via
%
%       lambda_L = 1/mu       = sqrt(2*max(m,n))*sigma,
%       lambda_S = lambda/mu  = sqrt(2)*sigma,
%       lambda   = lambda_S/lambda_L = 1/sqrt(max(m,n))   (the sqrt(2)'s and
%                                                          sigma cancel in the ratio).
%
% Hence the sqrt(2) factors appear in BOTH lambda_L and lambda_S here, while in
% the (lambda, mu) parametrization the factor 2 survives only in mu. The two
% choices are identical problems; only the normalization differs.
%
% [S+L_Embedding project modification] If options.L_true and options.S_true are
% supplied, each output entry gains two extra fields, L_true_relative_error and
% S_true_relative_error (plus a monitoring-excluded .time), logging the
% per-iteration relative error of the L / S iterates to the ground truth. This
% is diagnostic only and does not change the algorithm; omit the options to get
% the original behavior.

if ~exist('options','var')
    options = populate_vars([],A);
else
    options = populate_vars(options,A);
end

verbose = options.verbose;
lambda_S = options.lambda_S;
lambda_L = options.lambda_L;
MAX_ITER = options.max_itr;
ABSTOL = options.abs_tol;
RELTOL = options.rel_tol;


[m,n] = size(A);
out = struct([]);

% [S+L_Embedding project addition] Ground-truth error monitoring. When both
% options.L_true and options.S_true are supplied, log the per-iteration relative
% errors of the low-rank (L) and sparse (S) iterates to them as extra fields on
% the output struct (see header). Diagnostic only; does not change ADMM.
track_true_error = isfield(options, 'L_true') && isfield(options, 'S_true');
if track_true_error
    L_true      = options.L_true;
    S_true      = options.S_true;
    norm_L_true = norm(L_true, 'fro');
    norm_S_true = norm(S_true, 'fro');
    monitor_time_sum = 0;   % monitoring time, excluded from out(k).time
end

N = 3;
tic;

lambda = 1;
rho = 1/lambda;

E = zeros(m,n);
S = zeros(m,n);
L = zeros(m,n);
z   = zeros(m,N*n);
U   = zeros(m,n);

fprintf('\n%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n', 'iter', ...
    'r norm', 'eps pri', 's norm', 'eps dual', 'objective');

for k = 1:MAX_ITER

    B = avg(E, S, L) - A./N + U;

    % x-update
    E = (1/(1+lambda))*(E - B);
    S = prox_l1(S - B, lambda*lambda_S);
    L = prox_matrix(L - B, lambda*lambda_L, @prox_l1);
    x = [E S L];
    zold = z;
    z = x + repmat(-avg(E, S, L) + A./N, 1, N);

    % u-update
    U = B;

    % diagnostics, reporting, termination checks
    h.objval(k)   = objective(E, lambda_S, S, lambda_L, L);
    h.r_norm(k)   = norm(x - z,'fro');
    h.s_norm(k)   = norm(-rho*(z - zold),'fro');
    h.eps_pri(k)  = sqrt(m*n*N)*ABSTOL + RELTOL*max(norm(x,'fro'), norm(-z,'fro'));
    h.eps_dual(k) = sqrt(m*n*N)*ABSTOL + RELTOL*sqrt(N)*norm(rho*U,'fro');

    %if k == 1 || mod(k,2) == 0
    if verbose
        fprintf('%4d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n', k, ...
            h.r_norm(k), h.eps_pri(k), h.s_norm(k), h.eps_dual(k), h.objval(k));
    end
    %end

    % [S+L_Embedding project addition] Record this iteration's relative errors to
    % the true low-rank (L) and sparse (S) matrices as extra fields on the output
    % struct. Timed with its own tic so the monitoring cost is excluded from
    % out(k).time (same convention as R_Trust / R3_Trust / VBRPCA / GoDec).
    if track_true_error
        monitor_id = tic;
        out(k).L_true_relative_error = norm(L - L_true, 'fro') / norm_L_true;
        out(k).S_true_relative_error = norm(S - S_true, 'fro') / norm_S_true;
        monitor_time_sum = monitor_time_sum + toc(monitor_id);
        out(k).time = toc - monitor_time_sum;
    end

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

end
out(k).L = L;
out(k).S = S;
out(k).E = A - out(k).L - out(k).S;

end

function x = avg(varargin)
N = length(varargin);
x = 0;
for k = 1:N
    x = x + varargin{k};
end
x = x/N;
end

function p = objective(E, lambda_S, S, lambda_L, L)
p = norm(E,'fro').^2 + lambda_S*norm(S(:),1) + lambda_L*norm(svd(L),1);
end

function x = prox_l1(v, lambda)
% PROX_L1    The proximal operator of the l1 norm.
%
%   prox_l1(v,lambda) is the proximal operator of the l1 norm
%   with parameter lambda.

x = max(0, v - lambda) - max(0, -v - lambda);
end


function x = prox_matrix(v, lambda, prox_f)
% PROX_MATRIX    The proximal operator of a matrix function.
%
%   Suppose F is a orthogonally invariant matrix function such that
%   F(X) = f(s(X)), where s is the singular value map and f is some
%   absolutely symmetric function. Then
%
%     X = prox_matrix(V,lambda,prox_f)
%
%   evaluates the proximal operator of F via the proximal operator
%   of f. Here, it must be possible to evaluate prox_f as prox_f(v,lambda).
%
%   For example,
%
%     prox_matrix(V,lambda,prox_l1)
%
%   evaluates the proximal operator of the nuclear norm at V
%   (i.e., the singular value thresholding operator).

[U,S,V] = svd(v,'econ');
x = diag(prox_f(diag(S), lambda));
x = U*x;
x = x*V';



[U,s,V] = svd(v,'econ','vector');

s = prox_f(s, lambda);

idx = s ~= 0;

x2 = (U(:,idx) .* s(idx).') * V(:,idx)';
%x = U*diag(prox_f(diag(S), lambda))*V';
end

%% Function to populate options.
function [options] = populate_vars(options,A)

if isempty(options)
    options.verbose  = 1;
    options.lambda_S = 0.15*norm(A(:),inf);
    options.lambda_L = 0.15*norm(A);
    options.max_itr  = 100;
    options.abs_tol  = 1e-4;
    options.rel_tol  = 1e-2;
else

    if ~isfield(options, 'verbose')
        options.verbose = 1;
    end

    if ~isfield(options, 'lambda_S')
        options.lambda_S = 0.15*norm(A(:),inf);
    end

    if ~isfield(options, 'lambda_L')
        options.lambda_L = 0.15*norm(A);
    end

    if ~isfield(options, 'max_itr')
        options.max_itr  = 100;
    end

    if ~isfield(options, 'abs_tol')
        options.abs_tol  = 1e-4;
    end

    if ~isfield(options, 'rel_tol')
        options.rel_tol  = 1e-2;
    end

end
end