function x_rel = classicSBL(A,Y,noise_var)

A_ori = A;
N = size(A,2);
L = size(Y,2);
gammac = 1*ones(N,1); % initialize the prior

keep_listc = 1:N;

normc = 1;
itrc = 1;

r_rele = 1e-3;
thres = 1e-4;
R_max = 200;

while(itrc < R_max + 1 && normc > thres) % do the iteration

    gammal = zeros(N,1);
    gammal(keep_listc) = gammac;
    gammac_old = gammal;
    
    if min(abs(gammac)) < r_rele*max(abs(gammac))%
        indexc = find(gammac > r_rele*max(abs(gammac))); % should keep
        gammac = gammac(indexc);
        A = A(:,indexc);
        keep_listc = keep_listc(indexc);
    end
    gammad = gammac;

    % update the posterior
    Gammac = diag(gammac);
    invPhic = A'* (noise_var * eye(size(A,1)) + A * Gammac * A')^-1 * A;
    Sigma_x = Gammac - Gammac * invPhic * Gammac;
    x_re = noise_var^-1 * Sigma_x * A' * Y;
    
    % updating gamma
    Sigma_x_diag = real(diag(Sigma_x));
    gammac = Sigma_x_diag + sum(abs(x_re).^2,2)/L;
    
    % updating noise variance
    temp = sum(ones(length(gammad),1) - gammad.^-1 .* Sigma_x_diag);
    noise_var = (norm(Y - A*x_re,'fro')^2/L + noise_var*temp)/size(Y,1);

    
    gammal = zeros(N,1);
    gammal(keep_listc) = gammac;

    normc = norm(gammal-gammac_old)/norm(gammac_old);
    norm(gammal-gammac_old)/norm(gammac_old);
    itrc = itrc + 1; 
end
% re-estimate
Gammac = diag(gammal);
invPhic = A_ori'* (noise_var * eye(size(A_ori,1)) + A_ori * Gammac * A_ori')^-1 * A_ori;
Sigma_x = Gammac - Gammac * invPhic * Gammac;
x_rel = noise_var^-1 * Sigma_x * A_ori' * Y;

end
