function [fgrad_mean, fhess_mean, fmutual_mean] = ...
    mutual_mimo_deri_symm(M, N, H0, Rt_sqrt, gamma, Q, NSAMPLE, ...
                          ind1, ind2, I, Nx, Nr);

%------------------------------------------------------------
% Function to calculate the gradient and Hessian for general Hermitian
% covariance matrix case.
%
% Author: Mai Vu
% Date: 11/30/2003
% Revised: 04/12/2004
%------------------------------------------------------------

%------------------------------------------------------------
% INPUTS
%   M: number of receive antennas
%   N: number of transmit antennas
%   H0: channel mean
%   Rt_sqrt: square root of transmit correlation
%   gamma: SNR
%   lambda: the previous transmit covariance matrix eigenvalues
%   NSAMPLE: number of samples used in Monte-Carlo simulation
%   ind1, ind2: reorganizing indexes
%   I: identity matrix of size NxN
%   Nx: number of complex unknowns
%   Nr: number of real unknowns
%
% THE PROBLEM
%
%   The unknown covariance matrix Q is a Hermitian matrix of size N^2. The
%   number of complex unknowns is Nx = N*(N+1)/2, which is equivalent to a
%   number of real unknowns Nr = N^2. The reason that Nx and Nr are inputs
%   to the function is to speed up the run time as this function is called
%   often, hence saving the time to calculate Nx and Nr each time the
%   function is called.
%   
%   Objective function
%     f = -E[logdet(I + SQ)]
%     where  S = H'*H
%            H = H0 + Hw*Rt_sqrt
%   This function calculates the gradient and Hessian of f wrt the real and
%   imaginary entries of Q.
%
% OUTPUTS
%   fgrad_mean: average gradient
%   fhess_mean: average Hessian
%   fmutual_mean: average mutual information
%
%------------------------------------------------------------

% generate NSAMPLE of MxN Gaussian random matrices
H = (randn(M,N,NSAMPLE) + j*randn(M,N,NSAMPLE))/sqrt(2);

% create the variables
fgrad = zeros(Nx,1);
fhess = zeros(Nx,Nx);
fmutual = 0;

for k=1:NSAMPLE
  
  H_k = H0 + H(:,:,k)*Rt_sqrt;
  S_k = gamma*H_k'*H_k;
  P = I + S_k*Q;
  
  % function value
  fmutual = fmutual - real(log(det(P)));
  Y = inv(P);    
  Z = Y*S_k;
  
  Z12 = Z(ind1, ind2);
  Z11 = Z(ind1, ind1);
  Z22 = Z(ind2, ind2);
  
  % gradient
  fgrad_2 = diag(Z12(N+1:Nr,N+1:Nr));
  fgrad = fgrad - [diag(real(Z12)); (imag(fgrad_2))]; % *2
  
  W12 = transpose(Z12).*Z12;
  W11 = conj(Z22).*Z11;

  % hessian
  hess_11 = (real(W11)+real(W12));  % *2
  hess_22 = (real(W11)-real(W12));  % *2
  hess_12 = (imag(W12)-imag(W11));  % *2
  hess_21 = (imag(W12)+imag(W11));  % *2
  
  fhess = fhess + [hess_11 hess_12(:,N+1:Nr); hess_21(N+1:Nr,:) ...
                   hess_22(N+1:Nr,N+1:Nr)];
end

% take sample mean
fgrad_mean = fgrad/NSAMPLE;
fhess_mean = fhess/NSAMPLE;  
fmutual_mean = fmutual/NSAMPLE;