% =================================================================================
%  Matrix Algebraic Pursuit algorithms for low rank + sparse factorization - Demo
% =================================================================================
% Matrix Algebraic Pursuit (ALPS) algorithms are accelerated hard thresholding 
% methods on low rank recovery of linear inverse systems. Here, we consider the extension
% on low rank + sparse matrix recovery from an incomplete set of observations. In 
% particular, let
%             A   : m x n -> p x 1 linear mapping (p < mn)
%             L*  : m x n rank-k matrix component
%             M*  : m x n s-sparse matrix component with entries with
%                   arbitrary energy
%             X*  : m x n data matrix such that X* = L* + M*
%             e   : p x 1 additive noise vector
% then, y = A(X*) + e is the undersampled p x 1 measurement vector. Matrix ALPS 
% solve the following minimization problem 
%          minimize ||y - A(L + M)||_F^2    subject to rank(L) <= k. 
%                                                      ||M||_0 <= s.
%
% 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, 2011.
% ================================================================================
% 08/12/2011, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
% ================================================================================

clear all
close all
clc

addpath ALPS/;
addpath operators/;
addpath utility/;
addpath PROPACK/;

%% General parameters
n = 1000;                       % Number of columns in (L* + M*)
m = 600;                       % Number of rows in (L* + M*)
k = 80;                         % Rank of L*
s = ceil(.05*m*n);             % Sparsity level of M*
sigma = 0*10^-4;

a = -50;                       % Define range of magnitude of the entries in M*
b = 50;

% ALPS parameters
params.ALPSiters = 500;         % Maximum number of iterations
params.tol = 1e-5;              % Convergence tolerance
params.xpath = 1;               % Keep history log
params.svdMode = 'propack';     % SVD method - default mode for Matrix ALPS II - other options: 'svds', 'svd'
params.cg_tol = 1e-10;          % Conjugate gradients tolerance
params.cg_maxiter = 500;        % Maximum number of conjugate gradient iterations
params.svdApprox = 0;           % Set to 1 for column subset selection - really slow...
params.power = 2;               % Number of iterations in subspace range finder

maxits = 0;
MC = 1;                         % Monte Carlo iterations

A = @(z) A_id(z);               % Linear operator for Robust PCA
At = @(z) At_id(z, m, n);   

for i = 1:MC

    %% Generate Data
    L = randn(m, k)*randn(k, n);    % Left and right component such that rank(L) = k

    M = zeros(size(L));             % Sparse compoment with s non-zero entries
    indices = randperm(m*n); 
    indices = indices(1:s);
    r = a + (b-a).*rand(100,1);
    M(indices) = a + (b-a).*rand(size(indices));

    %% Generate noise and measurements
    e = randn(m*n,1);
    e = sigma*e/norm(e);

    Y = L + M + reshape(e, size(L));
    y = A(L + M) + e;
    
    %% Matrix ALPS I for Low rank + Sparse
    t0 = clock;
    [L1, M1, numiter1, L_path1, M_path1, LM_path1] = matrixALPSI_LS(y, A, At, m, n, k, s, params, L, M); 
    toc_time = etime(clock, t0);
    disp('==================================================================');
    str = sprintf('Matrix ALPS I for Low rank + Sparse terminated:');
    disp(str)
    str = sprintf('Low rank matrix error norm: %f ', norm(L-L1, 'fro')/norm(L, 'fro'));
    disp(str);
    str = sprintf('Sparse matrix error norm: %f ', norm(M-M1, 'fro')/norm(M, 'fro'));
    disp(str);
    str = sprintf('Total error norm: %f ', norm((M + L) - (M1 + L1), 'fro')/norm(M + L, 'fro'));
    disp(str);
    
    if (numiter1 > maxits)
        maxits = numiter1;
    end;
   
    time(1, i) = toc_time;
    num_iter(1, i) = numiter1;
    L_path(1, i, 1:numiter1) = L_path1;
    M_path(1, i, 1:numiter1) = M_path1;
    LM_path(1, i, 1:numiter1) = LM_path1;
    L_error(1, i) = norm(L-L1, 'fro')/norm(L, 'fro');
    M_error(1, i) = norm(M-M1, 'fro')/norm(M, 'fro');
    total_error(1, i) = norm((M + L) - (M1 + L1), 'fro')/norm(M + L, 'fro');
    
    %% Accelerated Matrix ALPS II with constant tau
    t0 = clock;
    [L2, M2, numiter2, L_path2, M_path2, LM_path2] = accelerated_matrixALPSII_constant_tau(Y', k, s, n, m, params, L', M'); 
    toc_time = etime(clock, t0);
    disp('==================================================================');
    str = sprintf('Accelerated Matrix ALPS II with memory - constant tau - terminated:');
    disp(str)
    str = sprintf('Low rank matrix error norm: %f ', norm(L-L2', 'fro')/norm(L, 'fro'));
    disp(str);
    str = sprintf('Sparse matrix error norm: %f ', norm(M-M2', 'fro')/norm(M, 'fro'));
    disp(str);
    str = sprintf('Total error norm: %f ', norm((M + L) - (M2' + L2'), 'fro')/norm(M + L, 'fro'));
    disp(str);
    
    if (numiter2 > maxits)
        maxits = numiter2;
    end;
       
    time(2, i) = toc_time;
    num_iter(2, i) = numiter2;
    L_path(2, i, 1:numiter2) = L_path2;
    M_path(2, i, 1:numiter2) = M_path2;
    LM_path(2, i, 1:numiter2) = LM_path2;
    L_error(2, i) = norm(L-L2', 'fro')/norm(L, 'fro');
    M_error(2, i) = norm(M-M2', 'fro')/norm(M, 'fro');
    total_error(2, i) = norm((M + L) - (M2' + L2'), 'fro')/norm(M + L, 'fro');
    
    %% Accelerated Matrix ALPS II with adaptive tau
    t0 = clock;
    [L3, M3, numiter3, L_path3, M_path3, LM_path3] = accelerated_matrixALPSII_adaptive_tau(Y', k, s, n, m, params, L', M'); 
    toc_time = etime(clock, t0);
    disp('==================================================================');
    str = sprintf('Accelerated Matrix ALPS II with memory - adaptive tau - terminated:');
    disp(str)
    str = sprintf('Low rank matrix error norm: %f ', norm(L-L3', 'fro')/norm(L, 'fro'));
    disp(str);
    str = sprintf('Sparse matrix error norm: %f ', norm(M-M3', 'fro')/norm(M, 'fro'));
    disp(str);
    str = sprintf('Total error norm: %f ', norm((M + L) - (M3' + L3'), 'fro')/norm(M + L, 'fro'));
    disp(str);
    
    if (numiter3 > maxits)
        maxits = numiter3;
    end;
       
    time(3, i) = toc_time;
    num_iter(3, i) = numiter3;
    L_path(3, i, 1:numiter3) = L_path3;
    M_path(3, i, 1:numiter3) = M_path3;
    LM_path(3, i, 1:numiter3) = LM_path3;
    L_error(3, i) = norm(L-L3', 'fro')/norm(L, 'fro');
    M_error(3, i) = norm(M-M3', 'fro')/norm(M, 'fro');
    total_error(3, i) = norm((M + L) - (M3' + L3'), 'fro')/norm(M + L, 'fro');    
end;

set(0, 'DefaultAxesFontSize', 16);
figure;
semilogy(squeeze(median(L_path(1,:,:),2)),'k','LineWidth',4); hold on
semilogy(squeeze(median(L_path(2,:,:),2)),'g','LineWidth',4);  hold on
semilogy(squeeze(median(L_path(3,:,:),2)),'r--','LineWidth',4); hold on
grid on;
xlabel('\# of iterations','FontSize', 18, 'Interpreter', 'latex')
ylabel(strcat(['$\|\mathbf{X}(i) $','- ','$ \mathbf{X}^\ast\|_F$']),'Interpreter', 'latex','FontSize', 18)
axis([1 maxits + 20 10^-6 1]);
semilogy((norm(e))*ones(1,maxits + 20),'k-.','LineWidth',2);

legend(strcat('Matrix ALPS I [', num2str(median(time(1,:))),']'), ...        
       strcat('Accelerated Matrix ALPS II - constant tau - [', num2str(median(time(2,:))),']'), ...;
       strcat('Accelerated Matrix ALPS II - adaptive tau - [', num2str(median(time(3,:))),']'));   
h  = legend;
set(h, 'interpreter', 'latex','FontSize', 18);
title('Low rank matrix component estimation');
shg;

set(0, 'DefaultAxesFontSize', 16);
figure;
semilogy(squeeze(median(M_path(1,:,:),2)),'k','LineWidth',4); hold on
semilogy(squeeze(median(M_path(2,:,:),2)),'g','LineWidth',4);  hold on
semilogy(squeeze(median(M_path(3,:,:),2)),'r--','LineWidth',4); hold on
grid on;
xlabel('\# of iterations','FontSize', 18, 'Interpreter', 'latex')
ylabel(strcat(['$\|\mathbf{X}(i) $','- ','$ \mathbf{X}^\ast\|_F$']),'Interpreter', 'latex','FontSize', 18)
axis([1 maxits + 20 10^-6 1]);
semilogy((norm(e))*ones(1,maxits + 20),'k-.','LineWidth',2);

legend(strcat('Matrix ALPS I [', num2str(median(time(1,:))),']'), ...        
       strcat('Accelerated Matrix ALPS II - constant tau - [', num2str(median(time(2,:))),']'), ...;
       strcat('Accelerated Matrix ALPS II - adaptive tau - [', num2str(median(time(3,:))),']'));   
h  = legend;
set(h, 'interpreter', 'latex','FontSize', 18);
title('Sparse matrix component estimation');
shg;

set(0, 'DefaultAxesFontSize', 16);
figure;
semilogy(squeeze(median(LM_path(1,:,:),2)),'k','LineWidth',4); hold on
semilogy(squeeze(median(LM_path(2,:,:),2)),'g','LineWidth',4);  hold on
semilogy(squeeze(median(LM_path(3,:,:),2)),'r--','LineWidth',4); hold on
grid on;
xlabel('\# of iterations','FontSize', 18, 'Interpreter', 'latex')
ylabel(strcat(['$\|\mathbf{X}(i) $','- ','$ \mathbf{X}^\ast\|_F$']),'Interpreter', 'latex','FontSize', 18)
axis([1 maxits + 20 10^-6 1]);
semilogy((norm(e))*ones(1,maxits + 20),'k-.','LineWidth',2);

legend(strcat('Matrix ALPS I [', num2str(median(time(1,:))),']'), ...        
       strcat('Accelerated Matrix ALPS II - constant tau - [', num2str(median(time(2,:))),']'), ...;
       strcat('Accelerated Matrix ALPS II - adaptive tau - [', num2str(median(time(3,:))),']'));   
h  = legend;
set(h, 'interpreter', 'latex','FontSize', 18);
title('Low rank + sparse matrix estimation');
shg;
