function [metrics] = HiIHT(y, s, sigma, params)
    N = params.M;
    n = params.D;
    % Initialization
    xk = zeros(N*n, 1);
    tic;
    for iter = 1:params.maxIt
        
        xk_old = xk;
        % Residual
        Ax = compute_kron_Ax(xk, params);
        r  = y - Ax;        
        g  = compute_kron_Atx(r, params);   

        % Threshold
        temp = xk + g/(params.M * params.N);
        x_thresh = hierarchicalThresholdOperator(temp, N, n, s, sigma);

        xk = x_thresh;

        if norm(xk-xk_old)/norm(xk) < 1e-5
            xk = x_thresh;
            break;
        end
    end
    time = toc;
    metrics={xk;
            time
            };
end
