mutual_mimo_deri_symm_barr.m

实现mimo预编码环境中容量的所有计算,包括不同的准则

function [fgrad_mean, fhess_mean, fmutual_mean] = ...
    mutual_mimo_deri_symm_barr(barr_t, M, N, H0, Rt_sqrt, gamma, Q_new, ... 
                               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
%   Q_new: the previous transmit covariance matrix 
%   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_new)] - logdet(Q_new)/barr_t
%     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 variables
fgrad = zeros(Nx,1);
fhess = zeros(Nx,Nx);
fmutual = 0;

% calculating the gradient and the Hessian by sample mean
% approximation. this is the most computationally intensive part.
for k=1:NSAMPLE
  
  H_k = H0 + H(:,:,k)*Rt_sqrt;
  S_k = gamma*H_k'*H_k;
  P = I + S_k*Q_new;
  
  % 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
  fgradZ_2 = diag(Z12(N+1:Nr,N+1:Nr));
  fgrad = fgrad - [diag(real(Z12)); imag(fgradZ_2)];

  W12 = transpose(Z12).*Z12;
  W11 = conj(Z22).*Z11;
    
  % hessian
  hess_11 = (real(W11)+real(W12));
  hess_22 = (real(W11)-real(W12));
  hess_12 = (imag(W12)-imag(W11));
  hess_21 = (imag(W12)+imag(W11));
  
  fhess = fhess + [hess_11 hess_12(:,N+1:Nr); hess_21(N+1:Nr,:) ...
                   hess_22(N+1:Nr,N+1:Nr)];

end

% add in the gradient, hessian and function values of the barrier part
R = inv(Q_new);
R12 = R(ind1, ind2);
R11 = R(ind1, ind1);
R22 = R(ind2, ind2);
fgradR_2 = diag(R12(N+1:Nr,N+1:Nr));
fgrad_mean = fgrad/NSAMPLE - ([diag(real(R12)); imag(fgradR_2)]/barr_t);

T12 = transpose(R12).*R12;
T11 = conj(R22).*R11;    
hess_11 = (real(T11)+real(T12));
hess_22 = (real(T11)-real(T12));
hess_12 = (imag(T12)-imag(T11));
hess_21 = (imag(T12)+imag(T11));
fhess_mean = fhess/NSAMPLE + ([hess_11 hess_12(:,N+1:Nr); hess_21(N+1:Nr,:) ...
                    hess_22(N+1:Nr,N+1:Nr)]/barr_t);  

fmutual_mean = fmutual/NSAMPLE - log(real(det(Q_new)))/barr_t;

% check that the hessian is non-negative
% D_temp = eig(fhess_mean);
% if any(D_temp < 0)
%   disp('Error! Negative Hessian');
%   D_temp
% end