function [L_hat, M_hat, numiter, L_path, M_path, LM_path] = matrixALPSI_LS(y, A, At, m, n, k, s, params, L, M)
% =========================================================================
%   Matrix ALPS I for Low rank + Sparse recovery algorithm - Beta Version
% =========================================================================
% Matrix Algebraic Pursuit (ALPS) algorithm with memory acceleration. 
% 
% Detailed discussion on the algorithm can be found in 
% [1] "Matrix ALPS: Accelerated Low rank and Sparse Matrix Reconstruction", 
% written by Anastasios Kyrillidis and Volkan Cevher, Technical Report, 2012.
% =========================================================================
% INPUT ARGUMENTS:
% y                         p x 1 undersampled measurement vector.
% A, At                     Linear operator and its adjoint:
%                           A: m x n |-> p x 1 mapping.
%                           At: p x 1 |-> m x n mappting
% k                         Rank prior information of L* or desired
%                           rank of computed solution.
% s                         Sparsity prior information of M* or desired
%                           sparsity of computed solution
% m, n                      Dimensions
% params                    Structure of parameters. These are:
%    tau,...                Momentum step size selection. Default value:
%                           tau = 0 for optimum step size selection. 
%    tol,...                Early stopping tolerance. Default value: tol =
%                           1-e5
%    ALPSiters,...          Maximum number of algorithm iterations. Default
%                           value: 300. 
%    cg_maxiter,...         Maximum iterations for Conjugate-Gradients method.
%    cg_tol,...             Tolerance variable for Conjugate-Gradients method.
%    xpath,...              Set history log to null.
%    svdApprox,...          Set to nonzero value in (0,1) for SVD approximation
% =========================================================================
% 
% OUTPUT ARGUMENTS:
% L_hat                     m x n recovered rank-k matrix.
% M_hat                     m x n recovered s-sparse matrix.
% numiter                   Number of iterations executed.
% X_path                    Keeps a series of computed m x n low-rank matrices 
%                           until the end of the iterative process. In
%                           case this information is not needed, please set
%                           params.xpath = 0
% =========================================================================
% 09/12/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
% =========================================================================
% cgsolve.m is written by Justin Romberg, Caltech, Oct. 2005.
%                         Email: jrom@acm.caltech.edu
% =========================================================================
% This work was supported in part by the European Commission under Grant 
% MIRG-268398 and ERC (Future Proof project). VC also would like 
% to acknowledge Rice University for his Faculty Fellowship.
% =========================================================================

%% Initialize to zero matrix
if (params.xpath == 1)
    L_path = zeros(1, params.ALPSiters);
    M_path = zeros(1, params.ALPSiters);
    LM_path = zeros(1, params.ALPSiters);
end;

length_y = length(y);

% Low rank part variables
L_cur = zeros(m, n);

QL_cur = zeros(m, n);
QM_cur = zeros(m, n);

ALM_cur = zeros(length_y, 1);
Ucur = [];
ULcur = [];

% Sparse part variables
M_cur = zeros(m, n);
M_prev = zeros(m, n);
Mcur = [];
UMcur = [];
complementary_Mi = ones(m,n);

options.tol = 10^-3;

I = eye(m,m);

i = 1;

%% SpaRCS with ALPS and Matrix ALPS
while (i <= params.ALPSiters)
    if (params.xpath == 1)
        LM_path(1,i) = norm(L_cur + M_cur - (L + M),'fro')/norm((L + M), 'fro');
        L_path(1,i) = norm(L_cur - L,'fro')/norm(L, 'fro');
        M_path(1,i) = norm(M_cur - M,'fro')/norm(M, 'fro');
    end;
        
    % Compute the residual      
    res = y - A(QL_cur + QM_cur);
    
    % Compute the gradient
    grad = At(res);   
    
    %% Low rank matrix part
                
    % Active subspace expansion step: Si_L (Di_L := Pk(P_{\mathcal{L}_i}^\bot grad))
    if (i == 1)
        [Uout, ~, ~] = lansvd(grad, k, 'L',options);
    else
        [Uout, ~, ~] = lansvd(ortho_UQL_i*grad, k, 'L',options); 
    end;
    
    Si_L = [Ucur Uout];
    
    % Error norm reduction via gradient descent
    proj_grad = Si_L*(Si_L'*grad);    
    mu = norm( proj_grad,'fro')^2/norm( A(proj_grad),2)^2;    
    Vi_L = QL_cur + (mu)*proj_grad;
    
    % Best rank-k subspace selection
    [UWi_L, SWi_L, VWi_L] = lansvd(Vi_L, k, 'L');  
    Wi_L = UWi_L*(SWi_L*VWi_L');
        
    % Debias via gradient descent
    res = y - A(Wi_L) - A(M_cur);        
    grad = At(res);    
    proj_grad = UWi_L*(UWi_L'*grad);
    xi = norm( proj_grad,'fro')^2/norm( A(proj_grad),2)^2;    
    L_prev = L_cur;
    L_cur = Wi_L + (xi)*proj_grad;
    ULprev = ULcur;
    ULcur = UWi_L;        
    
    Ucur = [ULcur ULprev];
    ortho_UQL_i = I - Ucur*Ucur';
    
    % Update current estimates ALM_cur, ALM_prev, QLM_cur
    ALM_prev = ALM_cur;
    ALM_cur = A(L_cur + M_prev);
    ALM = ALM_cur - ALM_prev;
    tau_LM = ((y - ALM_cur)'*(ALM))/norm(ALM, 2)^2; 
    QL_cur = L_cur  + tau_LM*(L_cur - L_prev);
    
    %% Sparse matrix part

    % Compute the residual
    res = y - A(QL_cur + QM_cur);
    
    % Compute the gradient
    grad = At(res);

    % Active support expansion step: Si_M (Di_M := Ps(P_{\mathcal{M}_i}^\bot grad))
    complementary_Mi(Mcur) = 0;
    [~, ind] = sort(abs(grad(:)).*complementary_Mi(:), 'descend');
    complementary_Mi(Mcur) = 1;
    Si_M = union(Mcur, ind(1:s));

    % Error norm reduction via gradient descent
    proj_grad = grad(Si_M);
    mu = norm(proj_grad, 2)^2 / norm(A_I(A, proj_grad, Si_M, m*n),2)^2;
    Vi_M = QM_cur(Si_M) + mu*proj_grad;

    % Best s-sparse support selection
    [~, ind] = sort(abs(Vi_M), 'descend');
    UWi_M = Si_M(ind(1:s));
    Wi_M = zeros(m, n); 
    Wi_M(UWi_M) = Vi_M(ind(1:s));

    % Debias via gradient descent
    res = y - A(Wi_M) - A(L_cur);  
    grad = At(res);
    proj_grad = grad(UWi_M);
    xi = norm(proj_grad, 2)^2/norm( A_I(A, proj_grad, UWi_M, m*n) ,2)^2;
    M_temp = Wi_M(UWi_M) + xi*proj_grad;
    M_prev = M_cur;
    M_cur = zeros(m,n); M_cur(UWi_M) = M_temp;
    UMprev = UMcur;
    UMcur = UWi_M;
    Mcur = union(UMprev, UMcur);
    
    % Update current estimates AM_cur, AM_prev, QM_cur
    ALM_prev = ALM_cur;
    ALM_cur = A(L_cur + M_cur);
    ALM = ALM_cur - ALM_prev;
    tau_LM = ((y - ALM_cur)'*(ALM))/norm(ALM, 2)^2;
    QM_cur = M_cur + tau_LM*(M_cur - M_prev);
    
    % Test stopping criterion
    if (i > 1) && (norm((L_cur + M_cur) - (L_prev + M_prev), 'fro') < params.tol*norm((L_cur + M_cur), 'fro'))
        break;
    end;
    i = i + 1;          
            
end;

L_hat = L_cur;
M_hat = M_cur; 
numiter = i;

if (params.xpath == 1)
    if (i > params.ALPSiters)
        LM_path = LM_path(1,1:numiter-1);
        L_path = L_path(1,1:numiter-1);
        M_path = M_path(1,1:numiter-1);
        numiter = numiter - 1;
    else
        LM_path = LM_path(1,1:numiter);
        L_path = L_path(1,1:numiter);
        M_path = M_path(1,1:numiter);
    end;
else LM_path = []; L_path = []; M_path = [];
end;