function [metrics] = HiHTP(y, A, s, sigma, params)
    
    N = params.M;
    n = params.D;
    matrix_dims = [n, N];
    % Initialization
    xk = zeros(N*n, 1);
    support = [];
    tic;
    for iter = 1:params.maxIt
        
        support_old = support;
        % 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);

        % Identify the support
        support = find(x_thresh ~= 0);
        
        % Solve least squares
        [n_indices, N_indices] = ind2sub(matrix_dims, support);
        A_sub = kr(A{2}(:,N_indices), A{1}(:,n_indices));           
        z_sub = pinv(A_sub) * y;

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


function K = kr(A, B)
[m,k] = size(A);
[n,k2] = size(B);
assert(k==k2,'Both factors must have same #columns');
K = zeros(m*n,k);
for j = 1:k
    K(:,j) = kron( A(:,j), B(:,j) );
end
end
