function [L_hat, M_hat, numiter, L_path, M_path, LM_path] = accelerated_matrixALPSII_constant_tau(Y, k, s, m, n, params, L, M)
% =====================================================================================
%   Accelerated Matrix ALPS I for Low rank + Sparse recovery algorithm - Beta Version
% =====================================================================================
% Accelerated Matrix Algebraic Pursuit (ALPS) algorithm with memory
% acceleration for low rank + sparse matrix recovery in the robust PCA
% problem:
%
%                          Y = L* + M* + e
% 
% In this memory scheme, we assume constant mommentum step size selection:
%                       
%                             tau = 1/4
%
% 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                         m x n measurement matrix.
% 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;

%initialization of L and S
res = Y;
L_cur = zeros(m, n);
M_cur = zeros(m, n);
QM = zeros(m, n);

L_prev = zeros(m, n);
M_prev = zeros(m, n);
i = 1;

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;
    
    %% Update of L
    Y2 = randn(n, k);
    for j = 1:params.power+1
        Y1 = res*Y2;
        Y2 = res'*Y1;
    end
    [Q, ~] = qr(Y2,0);
    L_cur = (res*Q)*Q';
    
    QL = L_cur + 1/4*(L_cur - L_prev);
    
    T = res - QL + QM;
    res = QL;
    
    %% Update of S
    [~, idx] = sort(abs(T(:)),'descend');
    M_cur = zeros(m, n);
    M_cur(idx(1:s)) = T(idx(1:s));
        
    QM = M_cur + 1/4*(M_cur - M_prev);
    
    T = T - QM;
    % 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;
    
    res = res + T;
    
    i = i + 1;
    L_prev = L_cur;
    M_prev = M_cur;
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;