function [A,mu,RelErr] = update_a_mu(A_m,mu_m,hp_samples,UL,ZL,US,ZS,w,rho,D,real_A)
    num_samples = length(hp_samples);
    %following the Algorthm 1: ADM4
    mu_numerator = zeros(D, 1);%init [D,1] all zeros; D -> dim of HP
    C = zeros(size(A_m));           % Initiation of the Matrix C (Equation 10 in the paper)
    A_Step = zeros(size(A_m)) + 2 * rho;
    B = zeros(size(A_m)) + rho*(UL-ZL) + rho*(US-ZS);% Add rho to B and its initialization (Equation 10 in the paper)
    for s = 1:num_samples           % Iterate through all the samples, 50,000
        timestamp = hp_samples(s).timestamp;% A vector, let's say 100 time stamps
        event = hp_samples(s).event;% A vector, each value a dim, same size as Time
        tc = timestamp(end);        % End time
        nc = length(event);
        dt = tc - timestamp;        % dt for all every time step
        for i = 1:nc                % Iterate through all the timestamps and events in ONE sample
            ui = event(i);          % Scalar
            ti = timestamp(i);      % Scalar
            int_g = kernel_int_g(dt, w);
            B(ui,:) = B(ui,:) + double(A_m(ui,:)>0).*repmat(int_g(i),[1,D]);
            pii = [];               % Initiation of the pii parameter is the Q objective function (Equation 8 in the paper)
            pij = [];               % Initiation of the pij parameter is the Q objective function (Equation 8 in the paper)
            ag_arr = [];            % Sum Auiuj*g(ti-tj)
            if i>1
                tj = timestamp(1:i-1);                        % Vector
                uj = event(1:i-1);                            % Vector
                ag_arr = A_m(uj, ui) .* kernel_g(ti - tj, w)';% Vector
            end          
            pii = mu_m(ui)./(mu_m(ui) + sum(ag_arr));    % Calculation of pii in Equation 8 in the paper
            if i>1
                pij = ag_arr./(mu_m(ui) + sum(ag_arr));  % Calculation of pij in Equation 8 in the paper
                if ~isempty(pij) && sum(pij(:))>0
                    for j = 1:length(uj)
                        uuj = uj(j);
                        C(uuj,ui) = C(uuj,ui) - pij(j,:);% Calculation of matrix C in Equation 10 in the paper
                    end
                end
            end
            mu_numerator(ui) = mu_numerator(ui) + pii;
        end                                           % for i = 1:nc (number of sequence in a sample)
    end                                               % for s = 1:num_samples (number of samples)
    mu = mu_numerator./(zeros(D, 1) + tc);            % Equation 9 in the paper
    
    A = (-B + sqrt(B.^2 - 4*A_Step.*C))./(2*A_Step);  % Equation 10 in the paper
                                                      % Given AmatA = zeros(size(A_m)) + 2 * rho;
    RelErr = real_err(real_A,A);                      % Storage of Iteration Relative Error
    fprintf('non-zero in mu= %d\n',nnz(mu));          % Testing of the sparse consition for mu
    fprintf('non-zero in C= %d\n',nnz(C));            % Testing of the sparse consition for C
    fprintf('non-zero in B = %d, non-zero in sqrt = %d\n',nnz(B),nnz(sqrt(B.^2 - 4*A_Step.*C)));
    fprintf('real error = %.4f, correlation = %.4f, #non-zero in A = %d\n\n',real_err(real_A,A),sum(abs(corr(real_A(:),A(:)))),nnz(A));
    A(isnan(A))=0;                                    % Refer to Equation 10 in the paper
    A(isinf(A))=0;
    %(A==0) = 0.00001;
end                                                   %end function