function result = iht(A, y, k)

    [~, n] = size(A);

    x_hat = zeros(n, 1);

    max_iter = 200;
    const = 1.9 / (norm(A)^2);
    
    tic;

    for iter = 1:max_iter
        
        x_prev = x_hat;

        residual = y - A * x_hat;
        x_hat = x_hat + const * A' * residual;

        [~, indices] = sort(abs(x_hat), 'descend');
        
        x_hat(indices(k+1:end)) = 0;

        if norm(x_hat - x_prev) / norm(x_hat) < 1e-6
            break;
        end
    end

    time_total = toc;
    x_re = x_hat;

    result = cell(2,1);
    result{1} = x_re;
    result{2} = time_total;

end