function X_est = global_omp_lowcom(A, Y, tol_coef)
% Solves Y = AX for a sparse matrix X.
% 
% Note that here X has no specific structured sparsity.

[M, N] = size(A);
[~, L] = size(Y);

X_est = zeros(N, L);

R = Y;
tol = tol_coef * norm(vec(Y));
col_norm = vecnorm(A);
A = A./ col_norm;

support_sets = cell(1, L);     % Stores the indices of atoms in the support
inv_Gram_matrices = cell(1, L); % Stores inv(A_sub' * A_sub)
A_corr_Y_vectors = cell(1, L);  % Stores A_sub' * Y(:,j)

k = 0;

while norm(R, 'fro') > tol
    k = k + 1;
    % Compute correlation of residual with dictionary atoms
    P = A' * R;

    % Find the (atom, measurement) pair with the max correlation
    [~, linear_idx] = max(abs(P(:)));
    [atom_idx, y_idx] = ind2sub(size(P), linear_idx);
    d_new = A(:, atom_idx);
    S_j = support_sets{y_idx};

    % Append the new atom
    S_j_new = [S_j; atom_idx];
    support_sets{y_idx} = S_j_new;
    
    if isempty(S_j)
        inv_G_new = 1;
        A_corr_Y_new = d_new' * Y(:, y_idx);
        
    else
        inv_G_old = inv_Gram_matrices{y_idx};
        A_corr_Y_old = A_corr_Y_vectors{y_idx};
        D = A(:, S_j); % The previous dictionary for this column
        
        u = D' * d_new;
        s = d_new' * d_new;
        
        % Schur complement
        schur_comp = s - u' * inv_G_old * u;

        if abs(schur_comp) < 1e-4
            fprintf('Warning: New atom is nearly linearly dependent.');
            break;
        end

        inv_s = 1 / schur_comp;
        
        % Use the blockwise matrix inversion formula to update the inverse Gram matrix
        inv_G_new = [inv_G_old + inv_G_old * u * inv_s * u' * inv_G_old, -inv_G_old * u * inv_s;
                     -inv_s * u' * inv_G_old,                            inv_s];
        
        A_corr_Y_new = [A_corr_Y_old; d_new' * Y(:, y_idx)];
    end

    inv_Gram_matrices{y_idx} = inv_G_new;
    A_corr_Y_vectors{y_idx} = A_corr_Y_new;

    x_solution = inv_G_new * A_corr_Y_new;
    X_est(S_j_new, y_idx) = x_solution;

    % Update residual
    R = Y - A * X_est;
end
X_est = X_est./col_norm';

end