function x = ls_log_sub(y,step,theta,lambda,n)

% This solves the subproblem
% min 0.5/step*||x - y||^2 + lambda*sum(log(abs(x) + theta))

% initialization
x = zeros(n,1); 
absy = abs(y);

% determine when f' = 0
delta = (theta - absy).^2 - 4*(lambda*step - theta*absy);
index1 = (delta >= 0);

if isempty(index1) == 0
    % locate all possible nonzero stationary points
    ytmp = absy(index1);
    t1 = max((ytmp - theta + sqrt(delta(index1)))/2,0);
    
    % compare function values; indices restricted to index1
    fun1 = 1/2*(t1 - ytmp).^2 + (lambda*step)*log(1 + t1/theta) - 1/2*ytmp.^2;
    inde1 = (fun1 < 0);
    
    % put entries back to the big/original vector
    index1(index1) = inde1;
    x(index1) = t1(inde1);
    x = sign(y).*x;
end

end
