function [X,err,fval_ls,iter,time] = LRS_sparse(A,k,s,lambda,opts)

% This function computes a stationary point for the following model
%
% min 0.5*||X - A||^2 + delta_{S}(X) + dist^2(X, Pi_k)/(2*lambda)
% where S = {X: ||vec(X)||_0 <= s} and Pi_k = {X: rank(X) <= k}
% Notice that in this case, the constraint {max_{i,j}|X_ij| <= tau} is not needed actually,
% since f(X) = 0.5*||X - A||_F^2 is level bounded.

% It uses proximal DCA with nonmonotone line search

%
% Input
%
% A           - noise matrix
% k           - rank bound
% s           - sparsity bound
% lambda      - smoothing parameter
% opts        - contains the following structure variables, default values are in brackets:
%               X0      - starting point [0]
%               tol     - tolerance [1e-6]
%               freq    - frequency for displaying output [5000]
%               maxiter - maximum number of iterations [inf]
%               disp    - display [1]: 1 - display per freq iterations
%                                  0 - no display
%               nonM       - memory in nonmonotone linesearch [5]
%
% Output
%
% X       - approximate stationary point
% err     - feasibility measurement: dist(X, Pi_k)
% fval_ls - ls value 0.5*||X - A||^2
% iter    - number of iterations
% time    - cputime

% further optimized on May 16, 2018

tstart = tic;
[m,n] = size(A); eta = 2; sigma = 1e-4; step_init = 1;

%% Parameters and defaults
if isfield(opts,'X0'),      X0      = opts.X0;      else X0 = zeros(m,n);end
if isfield(opts,'tol'),     tol     = opts.tol;     else tol = 1e-6;     end
if isfield(opts,'freq'),    freq    = opts.freq;    else freq = 5000;    end
if isfield(opts,'maxiter'), maxiter = opts.maxiter; else maxiter = inf;  end
if isfield(opts,'disp'),    disp    = opts.disp;    else disp = 1;       end
if isfield(opts,'nonM'),    nonM    = opts.nonM;    else nonM = 5;       end

X = X0;
lambdat = 1 + 1/lambda;
%% Compute function value and gradient
fval_ls = 0.5*norm(X - A,'fro')^2;
[U,S,V] = svd(X,'econ'); 
s0 = diag(S);
err = norm(s0(k+1:end));
fnew = fval_ls + err^2/(2*lambda);

g1 = lambdat*X - A;
U = U(:,1:k); S2 = S(1:k,1:k); V = V(:,1:k);   % projection to Pi_k
g2 = (U*S2*V')/lambda;
%% Set_up for nonmonotone line search
f1 = -inf(nonM,1);
f1(1) = fnew;
fmax = max(f1);

%% Main iteration
for iter = 1: maxiter
    fold = fnew;
    
    stepsize = step_init;
    Xnew = X - stepsize*(g1 - g2);
    tmp = Xnew(:); [~, ind] = sort(abs(tmp),'ascend'); tmp(ind(1:m*n-s)) = 0; Xnew = reshape(tmp,m,n); % projection to S
    % Compute function value
    fval_ls = 0.5*norm(Xnew - A,'fro')^2;
    [U,S,V] = svd(Xnew,'econ'); 
    s0 = diag(S);
    err = norm(s0(k+1:end));
    fnew = fval_ls + err^2/(2*lambda);
    count = 0;
    while fnew > fmax - sigma/2*norm(Xnew - X,'fro')^2
        stepsize = stepsize/eta;
        Xnew = X - stepsize*(g1 - g2);
        tmp = Xnew(:); [~, ind] = sort(abs(tmp),'ascend'); tmp(ind(1:m*n-s)) = 0; Xnew = reshape(tmp,m,n); % projection to S
        % Compute function value
        fval_ls = 0.5*norm(Xnew - A,'fro')^2;
        [U,S,V] = svd(Xnew,'econ');
        s0 = diag(S);
        err = norm(s0(k+1:end));
        fnew = fval_ls + err^2/(2*lambda);
        count = count + 1;
    end
    
    %% Display
     if mod(iter,freq) == 0 && disp == 1
        fprintf(' iter = %5d, fval = %8.6f,  fval_ls = %8.6f, ls_used = %3d\n', iter, fnew, fval_ls, count)
    end

    
    %% Check for termination
    if norm(Xnew - X,'fro')/max(1,norm(Xnew,'fro')) < tol*stepsize || abs(fnew - fold)/max(1,abs(fnew)) < 1e-12 
        if disp == 1
            fprintf(' TERMINATION \n')
            fprintf(' iter = %5d, fval = %8.6f, fval_ls = %8.6f, err = %2.1e\n', iter, fnew, fval_ls, err)
            if abs(fnew - fold)/max(1,abs(fnew)) < 1e-12
                fprintf(' Terminate by lack of progress in function value\n')
            end
        end
        X = Xnew;
        break
    end
    
    %% Barzilai-Borwein stepsize
    step_init = max(min(1/lambdat,1e8),1e-8);
    
    % update iterate, gradient and function value
    X = Xnew;
    g1 = lambdat*Xnew - A;
    U = U(:,1:k); S2 = S(1:k,1:k); V = V(:,1:k);  % projection to Pi_k
    g2 = (U*S2*V')/lambda;
    f1(mod(iter,nonM)+1,1) = fnew;
    fmax = max(f1);
end
time = toc(tstart);