function [fval,x_new,iter] = pro_gra_extra_log( A, b, nmA, lambda, tol, maxiter, re_freq, epsilon, dis_freq)

fprintf('\n Start proximal gradient algorithm with extrapolation\n')

% set matrix handles
Ahandle = @(x) A*x;
AThandle = @(x) A'*x;

% initialization
[~, n] = size(A);
L = nmA;

theta = 1;
theta0 = 1;
x_old = zeros(n, 1);
x_new = x_old;

num_restart = 0;
iter = 0;

while iter < maxiter
    
    % gradient of log function
    tmp = lambda * x_new./(epsilon*(epsilon + abs(x_new)));
    % extrapolation
    y = x_new + (theta0/theta - 1/theta) * (x_new - x_old);
    
    % proximal gradient
    y_tmp = y - (1/L)*(AThandle(Ahandle(y)-b) - tmp);
    x_old = x_new;
    x_new = soft_thresh(y_tmp,lambda/L/epsilon);
    
    % check for termination
    if norm(x_new - x_old)/max(1,norm(x_new)) < tol
        iter = iter + 1;
        z = Ahandle(x_new) - b;
        fval = 1/2*norm(z)^(2) + lambda*sum(log(1+abs(x_new)/epsilon));
        fprintf(' Termination: it_num = %d, restart_num = %d\n', iter, num_restart);
        break
    end
    
    
    % Update theta
    theta0 = theta;
    theta = (1 + sqrt(1+4*theta^2))/2;
    
    if re_freq < inf
        if (iter>0 && (mod(iter, re_freq) == 0 || (y-x_new)'*(x_new-x_old) > 0) )
            num_restart = num_restart + 1;
            theta0 = 1;
            theta = 1;
        end
    end
    
    if mod(iter, dis_freq) == 0
        z = Ahandle(x_new) - b;
        fval = 1/2*norm(z)^(2) + lambda*sum(log(1+abs(x_new)/epsilon));
        fprintf(' iter %d, fval %g, norm_suc_x %g, alpha %g\n', iter, fval, norm(x_new - x_old), (theta/theta0 - 1/theta));
    end
    
    iter = iter + 1;
    
end











