clear all; close all
randn('seed',2020);rand('seed',2020); 
%% setup
n = 50; r = 5;

outlier_ratio = 0.3;
m = ceil(outlier_ratio*n^2);
Uo = randn(n,r);
L = Uo*Uo';

E = zeros(n,n);

indx = randperm(n^2);

E(indx(1:m)) = 10*randn(m,1);

X = L + E;
L_fro = norm(L,'fro');

%%
maxiter = 3e3;
%% exact-parameterization L1
U_0 = 0.001*randn(n,r);
mu_0 = .002;
U = U_0;
%dist_over_constant(1) = norm(U*U' - L,'fro')/L_fro;
%L1residual_exact_constant(1) = sum(sum(abs(U*U' - X)));
for i = 1:maxiter
    grad = (sign(U*U'-X))*U + (sign(U*U'-X))'*U;
    mu = mu_0;
    U = U - mu*grad;
    dist_exact_constant(i) = norm(U*U' - L,'fro')/L_fro;
    L1residual_exact_constant(i) = sum(sum(abs(U*U' - X)));
end

%% over-parameterization L1
U_0 = 0.001*randn(n,n);
mu_0 = .002;
U = U_0;
%dist_over_constant(1) = norm(U*U' - L,'fro')/L_fro;
%L1residual_over_constant(1) = sum(sum(abs(U*U' - X)));
for i = 1:maxiter
    grad = (sign(U*U'-X))*U + (sign(U*U'-X))'*U;
    mu = mu_0;
    U = U - mu*grad;
    dist_over_constant(i) = norm(U*U' - L,'fro')/L_fro;
    L1residual_over_constant(i) = sum(sum(abs(U*U' - X)));
end


%% double over-parameterization
mu_U = .001; mu_S = mu_U*sqrt(n); 
Sp_0 = 0.001*randn(n,n); Sm_0 = Sp_0;
U = U_0; Sp = Sp_0; Sm = Sm_0;
%dist_double_over(1) = norm(U*U' - X,'fro')/L_fro;
%L1residual_double_over(1) = sum(sum(abs(U*U' - L)));
for i = 1:maxiter
    
    % compute gradient
    temp = U*U' + Sp.*Sp - Sm.*Sm - X;
    gradU = (temp + temp')*U;
    gradSp = 2*temp.*Sp;
    gradSm = -2*temp.*Sm;
     
    %gradient descent
    U = U - mu_U*gradU;
    Sp = Sp - mu_S*gradSp;
    Sm = Sm - mu_S*gradSm;
    
    % compute distance for the low-rank matrix
    dist_double_over(i) = norm(U*U' - L,'fro')/L_fro;
    L1residual_double_over(i) = sum(sum(abs(U*U' - X)));
end

%% plot the results
fontsize = 26;
plotStyle = {'b:','k--','r-','m-','g:'};
figure
semilogy(1:maxiter,dist_exact_constant,plotStyle{1},'linewidth',3); hold on
%legendInfo{1} = ['$\|\textup{exact-parameterization}\|_1$']; hold on
semilogy(1:maxiter,dist_over_constant,plotStyle{2},'linewidth',3);
semilogy(1:maxiter,dist_double_over,plotStyle{3},'linewidth',3);

legendInfo{1} = ['MF-$\ell_1 (r''=r)$'];
legendInfo{2} = ['MF-$\ell_1 (r'' = n)$'];
legendInfo{3} = ['Ours $(r'' = n)$'];
legend(legendInfo,'Interpreter','LaTex','Location','Best')
xlabel('Iteration','FontSize',fontsize,'FontName','Times New Roman');
%xlim([1 maxiter])
ylabel('$\|{\bf X}_\star - {\bf U\bf U}^\top\|_F/\|{\bf X}_\star\|_F$','FontSize',fontsize,'FontName','Times New Roman','Interpreter','LaTex');
set(gca,'YDir','normal')
set(gca, ...
    'LineWidth' , 2                     , ...
    'FontSize'  , fontsize              , ...
    'FontName'  , 'Times New Roman'         );
set(gcf, 'Color', 'white');
set(gcf, 'Position', [100, 100, 600, 400]);
grid on