function X = mmv_htp(A, Y, k)
maxItr = 100;
col_norm = vecnorm(A);
A = A./col_norm;

[~,n] = size(A);
[~,L] = size(Y);
X = zeros(n,L);
i = 0;
support = [];

while i < maxItr 
    
    support_old = support;
    r = Y - A*X;
    g = A'*r;

    temp = X + g;
    [~, idx] = sort(sum(abs(temp).^2,2), 'descend');
    support = idx(1:k);                 % 1×k

    Asub = A(:, support);                            % m×k
    Z = pinv(Asub) * Y;                              % k×L

    X_new = zeros(n, L);

    X_new(support, :) = Z;

    if isequal(support_old,support)
        X = X_new;
        break;
    end

    X = X_new;
    i = i + 1;
end

X = X./col_norm';

end