function [A_m,mu_m,Iteration_Err] = optimize_a_mu(hp_samples,D,w,rho,num_iter_1,num_iter_2,thold,real_A)
    % Input the training samples, parameters, number of iterations 
    % and the real A for iterative learning.
    A_m = rand(D,D);          % Start with a random matrix to make optimization
    A_m = 0.8 * A_m./max(abs(eig(A_m)));  % Similar with the real A, scale the spectral radius of A_m to 0.8
    mu_m = rand(D, 1)./D;     % Start with a random mu value to make optimization as well 
    UL = zeros(D,D);          % Matrix U's Parameter of Lowrank Property
    ZL = zeros(D,D);          % Matrix Z's Parameter of Lowrank Property
    US = zeros(D,D);          % Matrix U's Parameter of Sparse Property
    ZS = zeros(D,D);          % Matrix U's Parameter of Sparse Property

    for i = 0:num_iter_1      % Outer while loop for algorithm 1 in the paper
		rho = rho * 1.1;      % Iteratively change the value of rho each time
        for j = 0:num_iter_2  % Inner while loop for algorithm 1 in the paper
            % Main Algorithm mentioned above
            fprintf('No. %d outter while iteration | No. %d inner while iteration', i+1,j+1);
			[A_m,mu_m,RelErr] = update_a_mu(A_m,mu_m,hp_samples,UL,ZL,US,ZS,w,rho,D,real_A);
            Iteration_Err(j+1,i+1) = RelErr;   % Capture the error during each iteration
        end
        [s,v,d] = svd(A_m + UL);     % Update ZL via SVD (Singular Value Decomposition)
        v = v - thold/rho;
        v(v < 0) = 0;                % Get ride of negatives
        ZL = s * v * d';             % Update ZL (Z for LowRank)
        UL = UL + (A_m - ZL);        % Update UL (U for LowRank)

       tmp = (abs(A_m + US)-thold/rho);
       tmp(tmp < 0) = 0;
       ZS = (sign(A_m + US).*tmp);  % Update ZS (Z for Sparse)
       US = US + (A_m-ZS);          % Update US (U for Sparse)
    end
end
