function [metrics] = HiHTP(A, y, N, n, s, sigma, maxIter)
    % Initialization
    col_norm = vecnorm(A);
    A = A./col_norm;

    xk = zeros(N*n, 1);
    support = [];
    time = 0;
    for iter = 1:maxIter
        tic;
        support_old = support;
        % Residual
        r   = y - A*xk;        
        g   = A'*r;            

        % Threshold
        temp = xk + g;
        x_thresh = hierarchicalThresholdOperator(temp, N, n, s, sigma);

        % Identify the support
        support = find(x_thresh ~= 0);
        
        % Solve least squares
        A_sub = A(:, support);           
        z_sub = pinv(A_sub)*y;

        % Construct new
        x_new = zeros(N*n, 1);
        x_new(support) = z_sub;
        
        time = time + toc;
        if isequal(support_old,support)
            xk = x_new;
            break;
        end
        
        xk = x_new;
    end

    x_new = x_new ./ col_norm';
    metrics={x_new;
            time
            };
end
