function [x] = lp_prox_mapping(y, lambda, p)
% compute the proximal mapping of lp regularization
% min f(X) := lambda*||x||_p^p + 1/2 * ||x - y||^2 (0 < p <= 1)

iter = 0;
if p == 1
    x = sign(y) .* max(abs(y) - lambda, 0);
else
    b = abs(y); x = zeros(size(b));
    tau = (2*lambda*(1 - p))^(1/(2-p)) + lambda*p*(2*lambda*(1 - p))^((p-1)/(2-p));
    % Using Newton method to find the minimizer in (xl, +\infty), where x1 = (2*p*(1 - p)*lambda)^(1/(2-p))
    index = find(b > tau);
    tmpx = b(index); tmpY = b(index);
    gx = tmpx - tmpY + lambda*p*tmpx.^(p-1);
    while max(abs(gx)) > 1e-12
        hx = 1 + lambda*p*(p-1)*tmpx.^(p-2);
        tmpx = tmpx - gx./hx;
        gx = tmpx - tmpY + lambda*p*tmpx.^(p-1);
        iter = iter + 1;
    end
    x(index) = tmpx.*sign(y(index));
end

