% ====================================================================================
%  Matrix Algebraic Pursuit algorithms for low rank + sparse factorization - PDE 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.
% ===================================================================================
% 25/03/2012, by Anastasios Kyrillidis. anastasios.kyrillidis@epfl.ch, EPFL.
% ===================================================================================

clear all
close all
clc

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

N = 300;
x = linspace(-1,1,N);
[X,Y] = meshgrid(x,x);
 
% This matrix has a sv decay, but robust PCA can do better in total storage.
% A = exp(-sqrt(X.*X+Y.*Y)) + exp(-sqrt((X-0.5).*(X-0.5)+(Y-0.2).*(Y-0.2)));
 
%% Non-smooth function, C1 along the diagonal
C1 = abs(X-Y);
 
% Non-smooth function, infinite on the diagonal
% C1 = 1./abs(X-Y+1e-5);

%% Problem parameters
m = N;
n = N;

k = 18;
s = 19908;

params.ALPSiters = 3000;          % Maximum number of iterations
params.tol = 1e-10;   
params.xpath = 1;                 % Keep history log
params.svdMode = 'propack';       % SVD method
params.cg_tol = 1e-12;            % Conjugate gradients tolerance
params.cg_maxiter = 3000;         % Maximum number of conjugate gradient iterations
params.svdApprox = 0;
params.power= 2;
params.tau = 0;

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

Z = zeros(size(C1));
maxits = 0; 
i = 1;

%% Matrix ALPS I for Low rank + Sparse
t0 = clock;
[L1, M1, numiter1, ~, ~, LM_path1] = matrixALPSI_LS(y, A, At, m, n, k, s, params, C1, Z); 
toc_time = etime(clock, t0);
disp('==================================================================');
str = sprintf('Matrix ALPS I for Low rank + Sparse terminated:');
disp(str)
str = sprintf('Total error norm: %f ', norm(C1 - (M1 + L1), 'fro')/norm(C1, 'fro'));
disp(str);

if (numiter1 > maxits)
    maxits = numiter1;
end;

time(1, i) = toc_time;
num_iter(1, i) = numiter1;
LM_path(1, i, 1:numiter1) = LM_path1;
total_error(1, i) = norm(C1 - (M1 + L1), 'fro')/norm(C1, 'fro');

%% Accelerated Matrix ALPS II with constant tau
t0 = clock;
[L2, M2, numiter2, ~, ~, LM_path2] = accelerated_matrixALPSII_constant_tau(Y', k, s, n, m, params, C1, Z); 
toc_time = etime(clock, t0);
disp('==================================================================');
str = sprintf('Accelerated Matrix ALPS II with memory - constant tau - terminated:');
disp(str)
str = sprintf('Total error norm: %f ', norm(C1 - (M2' + L2'), 'fro')/norm(C1, 'fro'));
disp(str);

if (numiter2 > maxits)
    maxits = numiter2;
end;

time(2, i) = toc_time;
num_iter(2, i) = numiter2;
LM_path(2, i, 1:numiter2) = LM_path2;
total_error(2, i) = norm(C1 - (M2' + L2'), 'fro')/norm(C1, 'fro');

%% Accelerated Matrix ALPS II with adaptive tau
t0 = clock;
[L3, M3, numiter3, ~, ~, LM_path3] = accelerated_matrixALPSII_adaptive_tau(Y', k, s, n, m, params, C1, Z); 
toc_time = etime(clock, t0);
disp('==================================================================');
str = sprintf('Accelerated Matrix ALPS II with memory - adaptive tau - terminated:');
disp(str)
str = sprintf('Total error norm: %f ', norm(C1 - (M3' + L3'), 'fro')/norm(C1, 'fro'));
disp(str);

if (numiter3 > maxits)
    maxits = numiter3;
end;

time(3, i) = toc_time;
num_iter(3, i) = numiter3;
LM_path(3, i, 1:numiter3) = LM_path3;
total_error(3, i) = norm(C1 - (M3' + L3'), 'fro')/norm(C1, 'fro');    


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;
