    function x = shrinkFuncComp(y,alpha,qNorm, method)
    % solve f(x)=y for x yields shrinkage function

    f = @(x) x+alpha.*qNorm.*sign(x).*abs(x).^(qNorm-1);
    df = @(x) 1+alpha.*qNorm.*(qNorm-1).*abs(x).^(qNorm-2);

    maxN = 10000;
    tol = 1E-7;
    Y = y(:);

    tic;

    if strcmp(method, 'newton')
        x = Y;
        diff = f(x)-Y;
        dx = 1;
        n = 0;

        % use dx as condition due to high slope at x=0
        while n<=maxN && max(abs(dx))>tol
            n = n+1;
            tmp = x;

            % newton step, but do not let x change its sign
            % (f is strictly increasing + f(0) = 0)
            x(x>0) = max(tol*x(x>0), x(x>0) - diff(x>0)./df(x(x>0)));
            x(x<0) = min(tol*x(x<0), x(x<0) - diff(x<0)./df(x(x<0)));

            diff = f(x)-Y;
            dx = x - tmp;
        end
    else % bisection
        x = zeros(size(Y));

        lowerBounds = min(x, Y);
        upperBounds = max(x, Y);

        x = (lowerBounds+upperBounds)/2;
        diff = f(x)-Y;
        dx = 1;
        n = 0;

        while n<=maxN && max(abs(dx))>tol
            n = n+1;
            tmp = x;

            upperBounds(diff>0) = x(diff>0);
            lowerBounds(diff<0) = x(diff<0);

            x = (lowerBounds+upperBounds)/2;
            diff = f(x)-Y;
            dx = x - tmp;
        end
    end

    t = toc;

    if n>maxN
        fprintf('WARNING: shrinkFuncComp.m : maximum number of iterations for inversion reached!\n');
        fprintf('         n = %i, ||diff|| = %4.2e, ||dx|| = %4.2e, t = %4.4fs\n', n, max(abs(diff)), max(abs(dx)),t);
    end

    x = reshape(x,size(y));
    end