function [tau, x, R_tau, tau_opt] = gradient_descent(tau, lambda, k, R_tau, R_aver, gamma)
% This function updates the value of lambda by using gradient descent

x = 10 * ones(1, k);
x_new = lambda; % x_new = initial value of lambda
fun_g_new = sum(x_new .* R_tau) - sum(x_new) * R_aver;
fun_g = inf;
epsilon = 1e-2;
step = 0;

while abs(fun_g_new - fun_g) > epsilon * abs(fun_g_new) 
    x = x_new;
    step = step + 1;
    if rem(step, 10) ==0
        tau = tau /2;
    end
    fun_g = fun_g_new;
 
    g = R_tau - R_aver; % the gradient of g(lambda) over lambda
    x_new = x - tau * g; % update the value of x_new 
    x_new = max(x_new, 1e-3);
    % based on newly obtained value of lambda (here is x_new), update the
    % optimal valuse of tau, R(tau), and R_aver, which are related to the
    % gradient of g(lamda) in x_new
    [R_tau, tau_opt] = tau_opt_cal(k, x_new, gamma);
    fun_g_new = sum(x_new .* R_tau) - sum(x_new) * R_aver;
end
