Solving rate distortion function for a discrete memoryless source using CVX


#1

I want to solve the rate distortion function of a discrete memoryless source which translates to the minimization of mutual information subject to a distortion constraint. This is a classical convex optimization problem in information theory where mutual information is minimized over a conditional distribution whereas a source distribution is kept fixed. This problem is the ``dual’’ of the channel capacity convex optimization problem in information theory where one is interested to maximize with respect to an input probability distribution and the conditional distribution is kept fixed. That problem is given as an example in the illustrative list of examples in CVX guide. Can this be solved as well since it is a convex optimization problem? I wrote the following but I get an error which I paste in the very end of my script.

% We consider a discrete memoryless source, with input
% alphabet X(t) \in {1,…,n}, and output alphabet Y(t) \in {1,…,m}, for t = 1,2,…
% The relation between the input and output is given statistically:
% p_ij = Prob(Y(t)=i|X(t)=j), i=1,…,m, j=1,…,n
% The matrix P is called the channel transition matrix.
% The rate distortion R is given by
% R = min{ I(X;Y) | P >= 0, sum(P(:,1)) = 1, sum(P(:,2)) = 1, sum(rhoPx) <= D},
% I(X;Y) is the mutual information between X and Y, and it can be shown
% that: I(X;Y) = c’x - sum_{i=1}^m y_ilog_2(y_i)
% where c_j = sum_{i=1}^m p_ijlog_2(p_ij), j=1,…,m
%and
% P >= 0 means that the matrix we are searching for has positive entries,
% sum(P(:,1)) = 1, sum(P(:,2)) = 1 mean that the matrix is stochastic
% sum(rho * P
x) <= D is the distortion constraint set

%===========================================================
% Initialization
s=0;
rng(s,‘v5uniform’);

%============================================================
%
%Initial data
%
%============================================================
n = 2; %dim X
m = 2; %dim Y
D=2; %Distortion level

rho=[0 1;
1 0]; %distortion function

if (0)

x=rand(n,1); %non-uniform vector source
x = x./repmat(sum(x),n,1); %normalize source vector x

elseif (1)

x=ones(n,1);
x = x./repmat(sum(x),n,1); %generate uniform vector source

end

%Rate Distortion Function via CVX
cvx_begin
variable P(n,n)
c = sum(P.*log ( P) )’;
y = P *x;
minimize(sum(entr(y))/log(2)+(c’*x)/log(2))
P >= 0;
sum(P(:,1)) == 1;
sum(P(:,2)) == 1;
sum(rho *y)<=D
cvx_end
R = cvx_optval;

% Result
disp([‘The rate distortion function is: ’ num2str ( R) ’ bits.’])

%===================ERROR==========================================

Error using .* (line 173)
Disciplined convex programming error:
Cannot perform the operation: {real affine} .* {concave}

Error in cvx_rate_distortion_memoryless_source (line 45)
c = sum(P.*log ( P))’;

So apparently the problem boils down to how one can write mutual information in the objective function of the cvx if I understand correctly. But the question is how.


(Mark L. Stone) #2

help entr

entr Scalar entropy.
entr(X) returns an array of the same size as X with the unnormalized
entropy function applied to each element:
{ -X.*LOG(X) if X > 0,
entr(X) = { 0 if X == 0,
{ -Inf otherwise.
If X is a vector representing a discrete probability distribution, then
SUM(entr(X)) returns its entropy.

Disciplined convex programming information:
    entr(X) is concave and nonmonotonic in X. Thus when used in CVX
    expressions, X must be real and affine. Its use will effectively 
    constrain X to be nonnegative: there is no need to add an
    additional X >= 0 to your model in order to enforce this.

You should read http://cvxr.com/cvx/doc/funcref.html#new-functions as well as http://cvxr.com/cvx/doc/funcref.html#built-in-functions to find out about other CVX functions.