optimal_phase_shift.asv

function [v, R, history] = optimal_phase_shift(h_d, H_r, h_v, error_th,...
    sigma_n_sqr, p_sig, iteration_count, L, early_stop)
Phi = H_r * diag(h_v);
A = Phi'*Phi;
b = Phi'*h_d;
discrete_levels = linspace(0, 2*pi, L);
N = size(H_r,2);
I = ones(N,1);
I = randi([1 L],N,1);
% theta = zeros(N,1);
for n=1:N
    theta(n) = discrete_levels(I(n));
end
v = exp(1i * theta);
arg_max_theta = theta;
obj_value = achievable_rate(h_d + H_r*diag(v)*h_v, sigma_n_sqr,p_sig);
max_obj_value = obj_value;
history = zeros(iteration_count+1,2);
history(1,1) = obj_value;
history(1,2) = max_obj_value;
for iter = 1:iteration_count
    for n=1:N
        kappa_n = b(n);
        for l=1:N
            if l~=n
                kappa_n = kappa_n + A(n,l)*v(l);
            end
        end
        min_diff = 4*pi;
        for ang=discrete_levels
            diff = abs(angle(exp(1i*ang))-angle(kappa_n));
            if diff < min_diff
                min_diff = diff;
                theta(n) = ang;
            end
        end
    end
    v = exp(1i * theta);
    prev_obj_value = max_obj_value;
    obj_value=achievable_rate(h_d + H_r*diag(v)*h_v, sigma_n_sqr,p_sig);
    history(iter + 1,1) = obj_value;
    if max_obj_value < obj_value
        max_obj_value = obj_value;
        arg_max_theta = theta;
    end
    history(iter + 1,2) = max_obj_value;
    if early_stop==1 && abs(max_obj_value - prev_obj_value) < error_th
        break;
    end
end
R = max_obj_value;
v = exp(1i * arg_max_theta);
end