function x = Lp_proj(a,lambda,q)

% compute the Lq regularized projection problem 
% min \|x-a\|^2/2 + lambda*\|x\|^q
%
% modified on March 6, 2017: changed the while loop

b = abs(a);
n = max(size(b));
fmin = b.^2/2;
x = zeros(n,1);
xl = ((1/(q*(1-q))./lambda)).^(1/(q-2));
gxl = xl - b + q*(lambda.*(xl.^(q-1)));
I = find(gxl <0);
tmpa = b(I);
%xl = xl(I);
xu = b(I);
tmpl = lambda(I);

%tmpx = (xl+xu)/2;
tmpx = xu;
gx = tmpx - tmpa + q*(tmpl.*(tmpx.^(q-1)));
hx = 1 + (q*(q-1))*(tmpl.*(tmpx.^(q-2)));

while max(abs(gx)) > 1e-12
    tmpx = tmpx - gx./hx;
    gx = tmpx - tmpa + q*(tmpl.*(tmpx.^(q-1)));
    hx = 1 + (q*(q-1))*(tmpl.*(tmpx.^(q-2)));
end 
ftmpx = (tmpx-tmpa).^2/2 + tmpl.*(abs(tmpx).^q);
tmpf = fmin(I);
J = find (ftmpx < tmpf);
x(I(J)) = tmpx(J);

x = x.*sign(a);
return;