function [metrics] = shtp(A, y, N, s)
maxIter = 100;

% Normalization
col_norm = vecnorm(A);
A = A./col_norm;

% Initialization
xk = zeros(prod(N), 1);  % x^0 = 0
support = [];

tic;
for iter = 1:maxIter
    support_old = support;
    r   = y - A*xk;        % residual
    g   = A'*r;            % gradient

    temp = xk + g;
    x_thresh = KroSupThresholdOperator(temp, N, s);

    support = find(x_thresh ~= 0);

    A_sub = A(:, support);           
    z_sub = pinv(A_sub)*y;

    x_new = zeros(prod(N), 1);
    x_new(support) = z_sub;

    if isequal(support_old,support)
        xk = x_new;
        break;
    end

    xk = x_new;
end

time_total = toc;

xk = xk ./ col_norm';
metrics={'vector',xk;
    'time',time_total
    };
end
