N_iter = 1000;
Record_fitness = zeros(N_iter, 1);

for iter = 1:N_iter
    theta = struct('rho', rho, 'v', v, 'mu', mu);
    [K, grads] = get_K(x, p, theta, k0, c);
    
    Ky = K + sigma2_n*eye(N);
    Ky_inv_y = Ky\y; 
    
    ell = -y'*Ky_inv_y - log(det(Ky));
    Record_fitness(iter) = real(ell);
    
    d_ell_K = conj(Ky_inv_y)*(Ky_inv_y.') - inv(Ky);
    gmus = [grads.mu];
    d_ell_mu = 2*real(gmus*d_ell_K(:)); 
    vnt = d_ell_mu;
    
    s = 1;  a = 0.5;
    muN = mu + a*s*vnt;
    thetaN = struct('rho', rho, 'v', v, 'muN', muN);
    [K_N, gradsN] = get_K(x, p, thetaN, k0, c);
    KyN = K_N + sigma2_n*eye(N);
    Ky_inv_y_N = KyN\y; 
    d_ell_K_N = conj(Ky_inv_y_N)*(Ky_inv_y_N.') - inv(KyN);
    gmusN = [gradsN.muN];
    d_ell_mu_N = 2*real(gmusN*d_ell_K_N(:)); 
    vntN = d_ell_mu_N;
    vntN = vntN+(vntN'*(vntN))/(vnt'*vnt)*vnt;           % Fletcher-Reeves NCG
   % vntN = vntN+(vntN'*(vntN-vnt))/(vnt'*vnt)*vnt;      % Polak-Ribiere NCG
     
    mu = mu + vntN;
end