# Computation of Channel Capacity Based on Self-Concordant Functions

Authors: Da-gang Tian; Yi-qun Huang
Journal: Journal of Electrical and Computer Engineering
 (2012)

---

## Abstract

The computation of channel capacity is a classical issue in information theory. We prove that algorithms based on self-concordant functions can be used to deal with such issues, especially when constrains are included. A new algorithm to compute the channel capacity per unit cost is proposed. The same view is suited to the computation of maximum entropy. All the algorithms are of polynomial time.


---

## Body

## 1. Introduction

The computation of channel capacity is a classical issue in information theory. Much attention has been paid to it [1–3]. Despite some progress [3], the Arimoto-Blahut algorithm [4–7] is still the most used method. This method has been developed into a general algorithm, alternative minimization method [8], and has found applications in various fields [9–11]. Recently, the computation of channel capacity with constraints, such as the computation of capacity-cost function, is becoming increasingly important [5, 12–15]. The Arimoto-Blahut algorithm is very limited in dealing with such issues. Due to lack of appropriate means to eliminate the Lagrange multipliers, the method in [5], which is based on Arimoto-Blahut algorithm, cannot obtain the optimal point.The computation of channel capacity has always been an optimization problem, but no one optimization algorithm, such as the gradient method and Newton’s method, and so forth, has become the mainstream method. One reason is that so far no one optimization algorithm can be as effective as the Arimoto-Blahut algorithm. However, with advances in optimization theory and with the increasing importance of constrained channel capacity issues, the situation is changing.In this paper, some new algorithms based on self-concordant function are proposed. Besides the effectiveness—it is of polynomial time—the algorithm can compute constrained channel capacity easily. In Section2, we introduced the definition of the self-concordant function and the main properties, as well as some algorithms for convex optimization problems. In Section 3, we give some new algorithms for channel capacity with or without constraints and for channel capacity per unit cost. Because it is based on the self-concordant function, these algorithms are of polynomial time. Some conclusions are given in Section 4.


## 2. Self-Concordant Function and Convex Programming

For convenience, in this section, we make a brief introduction to self-concordant functions, including the definition, main propositions, and algorithms for convex optimization.


### 2.1. Self-Concordant Function

Self-concordant function theory, published by Nesterov and Nemirovskii in 1994 [16, 17], is of landmark importance for establishing general interior point polynomial-time algorithms for convex programming. The theory reveals the main property of the interior point method for linear programming initiated by Karmarkar in 1984 [18]: self-concordance of the barrier function.Suppose thatf:D→R is three times continuously differentiable and strictly convex, where the domain D is an open convex subset in Rn.Definition 1 (see [19]).
Letn=1. A convex function f:R→R is self-concordant if
(1)|f′′′(x)|≤2f′′(x)3/2
for all x∈D.Definition 2.
A functionf:Rn→R is self-concordant if it is self-concordant along every line in its domain, that is, if the function g(t)=f(x+tv) is a self-concordant function of t for all x∈D and for all v∈Rn.
Since linear and convex quadratic functions have zero third derivatives, they are trivial self-concordant functions. The most common nontrivial self-concordant function is the logarithmic function −logx.
We can get more self-concordant functions due to the following simple and useful combination rules [19].Proposition 3.
Iff1, f2 are self-concordant on D and a1≥1, a2≥1, then a1f1+a2f2 is self-concordant on D.Proposition 4.
Iff is self-concordant on D⊂Rn, and A:Rm→Rn is a linear map, b∈Rn, then f(Ax+b) is self-concordant on {x∈Rm∣Ax+b∈D}.Proposition 5.
Iffi is self-concordant on Di⊂Rni,i=1,2,…,m. Let D=D1×⋯×Dm. Then, f(x1,…,xm)=f1(x1)+⋯+fm(xm) is self-concordant on D.


### 2.2. Newton’s Method for Self-Concordant Functions

Consider the convex optimization problem(2)min⁡Df(x),
where D⊂Rn is a bounded, closed, convex subset with a nonempty interior. Let f(x) be a convex function.Damped Newton Algorithm
Given an initial pointx0∈int(D), with tolerance ε>0.Repeat
(1)
Compute the Newton step(3)Δxnt=-∇2f(x)-1∇f(x),
and Newton decrement
(4)λ(f,x)=(∇f(x)T∇2f(x)-1∇f(x))1/2.
(2) Stopping criterion.Quit if λ(f,x)2/2≤ε.
(3) Ifλ(f,x)>1/4, then
(5)x=x+11+λ(f,x)Δxnt,
else(6)x=x+Δxnt.
Iff(x) is a self-concordant function on D, then the damped Newton algorithm above possesses the following good properties:(i)
allx∈int(D) if x0∈int(D).(ii)
ifλ(f,x)>0.25, then f(xi)-f(xi+1)≥0.026.(iii)
if at some iterationi we have λ(f,xi)≤0.25, then we are in the region of quadratic convergence of the method, that is, for every j≥i,
(7)λ(f,xj+1)≤2λ2(f,xj)≤(12)λ(f,xj),f(xj)-p*≤λ2(f,xj)2(1-λ(f,xj)).(iv)
the number of Newton steps required to find a pointx∈int(D) with λ(f,xi)≤ε<1 is bounded by(8)C((f(x(0))-p*)+log2log⁡2(1ε))
for some constant C.
For the convex optimization problem with equality constraintsAx=b, Newton’s method and the bound (8) are still valid as long as the initial point x0∈int(D) satisfies Ax0=b, and the Newton step Δxnt and decrement λ(f,x) are computed as follows:
(9)[∇2f(x)ATA0][Δxntw]=[-∇f(x)0],λ(f,x)=(ΔxntT∇2f(x)Δxnt)1/2.
For general convex optimization problems are(10)min⁡f0(x),s.tfi(x)≤0,i=1,…,m,Ax=b,
where f0,…,fm:Rn→R are convex and twice continuously differentiable, and A is a p×n matrix with rank (A)=p<n. By means of the logarithmic barrier function, (10) can be formulated approximately as an equality constrained problem as follows:
(11)min⁡tf0(x)+ϕ(x),s.tAx=b,
where
(12)ϕ(x)=-∑i=1mlog⁡⁡(-fi(x))
is called the logarithmic barrier function for (10). Let x*(t) be the optimal point of (11) and let p* be the optimal value of (10), then f0(x*(t))-p*≤m/t. Therefore, we can simply take t=m/ε and solve the equality constrained problem (11) using Newton’s method to get an ε-solution of (10), and if the objective function of (11) is self-concordant, then the number of Newton’s iterate is bounded by (8). However, this method usually requires a large t, which may bring about a large f(x(0))-p* and a large number of Newton’s iterate, thus it is rarely used. A commonly used method is the path-following method as follows.Path-Following Method
Given strictly feasiblex,t=t0>0,μ>1, with tolerance ε>0.Repeat
(1)
Centering step.Computex*(t) by minimizing tf0(t)+φ, subject to Ax=b, starting at x.(2)
Update.x=x*(t).(3)
Stopping criterion.Quit if m/t<ε.(4)
Increaset. t=μt.If for anyt>0,f(x)=tf0(x)+φ(x) is self-concordant, then the total number of Newton steps in the path-following method, not counting the initial centering step, is(13)N=[log⁡(m/(t0ε))log⁡μ]×(m(μ-1-log⁡⁡μ)+log2log2(1ε)).The term in brackets is the iterate number oft from t0 to m/ε, and the term in parentheses is the iterate number of Newton’s method per centering step.


## 2.1. Self-Concordant Function

Self-concordant function theory, published by Nesterov and Nemirovskii in 1994 [16, 17], is of landmark importance for establishing general interior point polynomial-time algorithms for convex programming. The theory reveals the main property of the interior point method for linear programming initiated by Karmarkar in 1984 [18]: self-concordance of the barrier function.Suppose thatf:D→R is three times continuously differentiable and strictly convex, where the domain D is an open convex subset in Rn.Definition 1 (see [19]).
Letn=1. A convex function f:R→R is self-concordant if
(1)|f′′′(x)|≤2f′′(x)3/2
for all x∈D.Definition 2.
A functionf:Rn→R is self-concordant if it is self-concordant along every line in its domain, that is, if the function g(t)=f(x+tv) is a self-concordant function of t for all x∈D and for all v∈Rn.
Since linear and convex quadratic functions have zero third derivatives, they are trivial self-concordant functions. The most common nontrivial self-concordant function is the logarithmic function −logx.
We can get more self-concordant functions due to the following simple and useful combination rules [19].Proposition 3.
Iff1, f2 are self-concordant on D and a1≥1, a2≥1, then a1f1+a2f2 is self-concordant on D.Proposition 4.
Iff is self-concordant on D⊂Rn, and A:Rm→Rn is a linear map, b∈Rn, then f(Ax+b) is self-concordant on {x∈Rm∣Ax+b∈D}.Proposition 5.
Iffi is self-concordant on Di⊂Rni,i=1,2,…,m. Let D=D1×⋯×Dm. Then, f(x1,…,xm)=f1(x1)+⋯+fm(xm) is self-concordant on D.


## 2.2. Newton’s Method for Self-Concordant Functions

Consider the convex optimization problem(2)min⁡Df(x),
where D⊂Rn is a bounded, closed, convex subset with a nonempty interior. Let f(x) be a convex function.Damped Newton Algorithm
Given an initial pointx0∈int(D), with tolerance ε>0.Repeat
(1)
Compute the Newton step(3)Δxnt=-∇2f(x)-1∇f(x),
and Newton decrement
(4)λ(f,x)=(∇f(x)T∇2f(x)-1∇f(x))1/2.
(2) Stopping criterion.Quit if λ(f,x)2/2≤ε.
(3) Ifλ(f,x)>1/4, then
(5)x=x+11+λ(f,x)Δxnt,
else(6)x=x+Δxnt.
Iff(x) is a self-concordant function on D, then the damped Newton algorithm above possesses the following good properties:(i)
allx∈int(D) if x0∈int(D).(ii)
ifλ(f,x)>0.25, then f(xi)-f(xi+1)≥0.026.(iii)
if at some iterationi we have λ(f,xi)≤0.25, then we are in the region of quadratic convergence of the method, that is, for every j≥i,
(7)λ(f,xj+1)≤2λ2(f,xj)≤(12)λ(f,xj),f(xj)-p*≤λ2(f,xj)2(1-λ(f,xj)).(iv)
the number of Newton steps required to find a pointx∈int(D) with λ(f,xi)≤ε<1 is bounded by(8)C((f(x(0))-p*)+log2log⁡2(1ε))
for some constant C.
For the convex optimization problem with equality constraintsAx=b, Newton’s method and the bound (8) are still valid as long as the initial point x0∈int(D) satisfies Ax0=b, and the Newton step Δxnt and decrement λ(f,x) are computed as follows:
(9)[∇2f(x)ATA0][Δxntw]=[-∇f(x)0],λ(f,x)=(ΔxntT∇2f(x)Δxnt)1/2.
For general convex optimization problems are(10)min⁡f0(x),s.tfi(x)≤0,i=1,…,m,Ax=b,
where f0,…,fm:Rn→R are convex and twice continuously differentiable, and A is a p×n matrix with rank (A)=p<n. By means of the logarithmic barrier function, (10) can be formulated approximately as an equality constrained problem as follows:
(11)min⁡tf0(x)+ϕ(x),s.tAx=b,
where
(12)ϕ(x)=-∑i=1mlog⁡⁡(-fi(x))
is called the logarithmic barrier function for (10). Let x*(t) be the optimal point of (11) and let p* be the optimal value of (10), then f0(x*(t))-p*≤m/t. Therefore, we can simply take t=m/ε and solve the equality constrained problem (11) using Newton’s method to get an ε-solution of (10), and if the objective function of (11) is self-concordant, then the number of Newton’s iterate is bounded by (8). However, this method usually requires a large t, which may bring about a large f(x(0))-p* and a large number of Newton’s iterate, thus it is rarely used. A commonly used method is the path-following method as follows.Path-Following Method
Given strictly feasiblex,t=t0>0,μ>1, with tolerance ε>0.Repeat
(1)
Centering step.Computex*(t) by minimizing tf0(t)+φ, subject to Ax=b, starting at x.(2)
Update.x=x*(t).(3)
Stopping criterion.Quit if m/t<ε.(4)
Increaset. t=μt.If for anyt>0,f(x)=tf0(x)+φ(x) is self-concordant, then the total number of Newton steps in the path-following method, not counting the initial centering step, is(13)N=[log⁡(m/(t0ε))log⁡μ]×(m(μ-1-log⁡⁡μ)+log2log2(1ε)).The term in brackets is the iterate number oft from t0 to m/ε, and the term in parentheses is the iterate number of Newton’s method per centering step.


## 3. Algorithms Based on Self-Concordant Functions for Channel Capacity

In this section, it is proved that the channel capacity, with or without constrains, can be computed by the path-following method in polynomial time. Meanwhile, we prove that channel capacity per unit cost is a single-peak function of the expected cost.


### 3.1. Channel Capacity without Constrains

LetX={x1,…,xm} and Y={y1,…,yn} be the input and output alphabets, respectively. Let p(x)=(p(x1),…,p(xm)) and q(y)=(q(y1),…,q(yn)) be the distributions of X and Y, respectively. Let T=(p(yj∣xi))m×n=(pij)m×nbe a transition matrix from X to Y. Hence,(14)q(yj)=p(x1)p1j+⋯+p(xm)pmj,j=1,…,n.IfI(X;Y) is the mutual information between X and Y, then(15)I(X;Y)=∑i,jp(xi)p(yj∣xi)log⁡p(yj∣xi)q(yj)(16)=∑i,jp(xi)p(yj∣xi)log⁡p(xi∣yj)p(xi).Channel capacity is defined as follows:(17)C=max⁡p(x)I(X;Y).
Equation (16) is a function of p(x) and p(x|y) and we denote it by I(p(x);p(x|y)).The Arimoto-Blahut algorithm utilized (16) as follows. In (16), both p(x) and p(x|y) are unknown, but given p(x), (16) can get its maximum at(18)p(xi∣yj)=p(xi)p(yj∣xi)∑k=1mp(xk)p(yj∣xk),i=1,…,m.j=1,…,n,
and given p(x|y), (16) can get its maximum at(19)p(xi)=∏jp(xi∣yj)p(yj∣xi)∑i=1m∏jp(xi∣yj)p(yj∣xi),i=1,…,m.
If p0(x)=(1/m,…,1/m), then the channel capacity can be approximated by computing p0(x|y),p1(x),p1(x|y),… through (18) and (19), alternatively. Let C(N+1,N) be the mutual information I(pN+1(x);pN(x|y)), then,(20)C-C(N+1,N)≤log⁡mN.In order to get an algorithm based on the self-concordant function, we utilize (15). Let(21)ci=-∑j=1np(yj∣xi)log⁡p(yj∣xi).
By (15), the channel capacity can be expressed as follows:(22)min⁡∑i=1mcipi+∑j=1nqjlog⁡⁡qj,s.tqj=∑i=1mpip(yj∣xi),j=1,…,n,p1+⋯+pm=1,p1>0,…,pm>0.
It is a standard convex optimization problem. Unfortunately, the objective function of (22) is not self-concordant, even when the logarithmic barrier −log pi is added.It seems that we should addqj>0 as constraints for (22) since there are logqj in the objective function of (22). But the qj are not independent variables. It can be computed by (14) and maintain positive as long as all the pi>0,i=1,2,…,m, hold. Nevertheless, we still add the logarithmic barrier −log qj in the objective of (22) to get a self-concordant objective function.Fort>0, consider the convex optimization problem as follows:(23)min⁡t∑i=1mcipi+t∑j=1nqjlog⁡qj-∑i=1mlog⁡pi-∑j=1nlog⁡qj,s.tqj=∑i=1mpip(yj∣xi),j=1,…,n,p1+⋯+pm=1.
In order to show the self-concordance of the objective function of (23), we need Proposition 6.Proposition 6.
For anyt>0,f(x)=txlogx-logx is a self-concordant function on x>0.Proof.
We have(24)f′(x)=tlog⁡x+t-1x,f′′(x)=tx+1x2,|f′′′(x)|=tx2+2x3.
Thus
(25)|f′′′(x)|f′′(x)3/2=(tx+2)(tx+1)3/2≤2.By Propositions3–6, the objective function of (23) is self-concordant, so we can solve it by the path following method, and the number of Newton iterations is bounded by (13).Equation (13) shows that solving (23) by the path following method is a polynomial time algorithm. In addition, a main advantage is that the algorithm can deal with constrains very easily.


### 3.2. Channel Capacity with Constrains and Channel Capacity Per Unit Cost

Consider the channel capacity with constraints as follows:(26)max⁡I(X;Y),s.t∑j=1maijpj≤Ei,i=1,…,r,p1+⋯+pm=1,p1≥0,…,pm≥0.
For r=2, [5] gave some algorithms based on the Arimoto-Blahut method. Because one cannot eliminate the Lagrange multiplier, the algorithms in [5] cannot get the optimal solution for (26).By adding logarithmic barriers for the inequality constraints in the objective of (26) as follows, we can get a polynomial time algorithm for (26):(27)min⁡t∑i=1mcipi+t∑j=1n∑i=1mpip(yj∣xi)log⁡∑i=1mpip(yj∣xi)-∑i=1mlog⁡pi-∑j=1nlog⁡∑i=1mpip(yj∣xi)-∑k=1rlog⁡(Ek-∑i=1mpiari),s.tp1+⋯+pm=1.
The objective function of (27) is self-concordant, so we can solve (27) by the path-following method.Channel capacity per unit cost is one of the typical problems of channel capacity with constraints [15, 20–22]. By [15], channel capacity per unit cost can be computed by(28)C=sup⁡β>0C(β)β,
where C(β) is the capacity-cost function, which is the solution of the following constrained problem [12]:(29)max⁡I(X;Y),s.t∑j=1majpj≤β,p1+⋯+pm=1,p1≥0,…,pm≥0.
Equation (29) is a special case of (26) (r=1). Here we discuss the channel capacity per unit cost in a more general manner.Theorem 7.
LetD⊆Rn be a convex set. Let F(x) be a concave function on D and G(x) a convex function on D. Let a=inf{G(x)|x∈D},b=sup{G(x)|x∈D}, and suppose a≥0. Let C(β)=sup{F(x)|G(x)≤β}. Then,(1)
C(β)is a concave increasing function,(2)
let(30)ρ(β)=C(β)β,
then ρ(β) is either a unimodal function or a monotone function. 												(31)sup⁡DF(x)G(x)=sup⁡a<β<bC(β)β.Proof.
(1) Letβ1>β2, then {x|G(x)≤β2}⊆{x|G(x)≤β1}. Therefore, C(β2)≤C(β1), For any ε>0, there are x1 such that G(x1)≤β1 and C(β1)=F(x1)+ε, and there are x2 such that G(x2)≤β2 and C(β2)=F(x2)+ε. Since G(x) is convex and F(x) is concave, we get
(32)G(αx1+(1-α)x2)≤αG(x1)+(1-α)G(x2)≤αβ1+(1-α)β2,αC(β1)+(1-α)C(β2)=αF(x1)+(1-α)F(x2)+ε≤F(αx1+(1-α)x2)+ε≤C(αβ1+(1-α)β2)+ε.
Since ε is arbitrary, C(β) is concave.
(2) Without loss of generality, we can assumeC(β) is differentiable. If not, we can approach it by a differentiable function to any precision [23]. Suppose there is a β0∈(a,b), such that
(33)ρ′(β)|β=β0=0,thatisC′(β0)β0-C(β0)=0.
Letβ1>β0; since C(β) is concave, we have C′(β1)≤C′(β0). Let
(34)ρ′(β)|β=β0=0,thatisC′(β0)β0-C(β0)=0.
It is the equation of a straight line that has slope C′(β1) and passes through the (β0,C(β0)). Since C(β) is increasing and concave, we get y1(β1)≤C(β1) and
(35)C′(β1)-C(β1)β1≤C′(β1)-y1(β1)β1=C′(β1)β0β1-C(β0)β1≤C′(β0)β0β1-C(β0)β1=0,
that is, ρ′(β1)≤0.
Letβ2<β0, since C(β) is concave, C′(β2)≥C′(β0). Let
(36)y0(β)=C′(β0)(β-β0)+C(β0).
It is the equation of a straight line that has slope C′(β0) and passes through (β0,C(β0)). Since C(β) is increasing and concave, we get y0(β2)≥C(β2) and
(37)C′(β2)-C(β2)β2≥C′(β2)-y0(β2)β2≥C′(β0)-y0(β2)β2=C′(β0)β0β2-C(β0)β2=0,
that is, ρ′(β2)≥0.
Therefore,β0 is a maximum point of function ρ(β)=C(β)/β.
(3) If(38)sup⁡DF(x)G(x)=+∞,
then for any N>0, there are x0, such that
(39)F(x0)G(x0)>N.
Let G(x0)=β0, then C(β0)=sup{F(x)|G(x)≤β0}≥F(x0), so
(40)sup⁡a<β<bC(β)β≥C(β0)β0≥F(x0)G(x0)>N.
Since N is arbitrary, we get
(41)sup⁡a<β<bC(β)β=+∞.
In turn, if
(42)sup⁡a<β<bC(β)β=+∞,
then for any N>0, there are β0, a<β0<b, such that
(43)C(β0)β0>N.
Since
(44)sup⁡DF(x)G(x)≥sup⁡G(x)≤β0F(x)G(x)≥sup⁡G(x)≤β0F(x)β0=C(β0)β0>N
so
(45)sup⁡DF(x)G(x)=+∞.
Therefore, if
(46)sup⁡DF(x)G(x)<+∞
is valid, then so is
(47)sup⁡a<β<bC(β)β<+∞
and vice versa.
Suppose (46) is valid, for any ε>0, there are x0∈D, such that
(48)F(x0)G(x0)+ε=sup⁡DF(x)G(x).
Let G(x0)=β0, then
(49)sup⁡DF(x)G(x)=F(x0)G(x0)+ε≤C(β0)β0+ε≤sup⁡a<β<bC(β)β+ε.
Since ε is arbitrary, we get
(50)sup⁡DF(x)G(x)≤sup⁡a<β<bC(β)β.
On the other hand, since (47) is valid too, for any ε>0, there are a<β0<b, such that
(51)C(β0)β0+ε=sup⁡a<β<bC(β)β.
Let x0∈D, such that G(x0)≤β0 and F(x0)+εβ0=sup{F(x)|G(x)≤β0}=C(β0), then
(52)sup⁡DF(x)G(x)≥sup⁡G(x)≤β0F(x)G(x)≥F(x0)β0=sup⁡a<β<bC(β)β-2ε.
Since ε is arbitrary, we get
(53)sup⁡DF(x)G(x)≥sup⁡a<β<bC(β)β.IfF(x)=I(X,Y) is the average mutual information defined by (15) and G(x)=E[b(X)], then (3) of Theorem 7 is Theorem  2 in [15].In addition to average mutual information, entropy is another frequently used function in applications [24]. It is obvious that Theorem 7 is valid when F(x) is an entropy function (as a function of the distribution of input symbols). It is useful to point out that Theorem 7 is still valid when G(x) is a positive definite quadratic function [25].Making use of (2) of Theorem7, we can write an algorithm to compute the sup(F(x)/G(x))as follows. The algorithm is based on the 0.618 method used to locate the optimal point of a single peak function.Algorithm 8 (0.618 method).
Leta=inf{G(x)|x∈D},b=sup{G(x)|x∈D}, and assume a≥0.
(1) Given search interval[a,b], and tolerance ε>0. compute λ and μ as follows:
λ=a+0.382(b-a),μ=a+0.618(b-a); evaluate C(λ)/λ and C(μ)/μ. Let k=1.
(2) IfC(λ)/λ>C(μ)/μ, then go to 3, else go to 4.
(3) Ifb-λ≤ε, then output μ and C(μ)/μ, else, let a=λ,λ=μ,C(λ)/λ=C(μ)/μ,μ=a+0.618(b-a); evaluate C(μ)/μ, go to 5.
(4) Ifμ-a≤ε, then output λ and C(λ)/λ, else, let b=μ,μ=λ,C(μ)/μ=C(λ)/λ,λ=a+0.382(b-a); evaluate C(λ)/λ, go to 5.
(5)k=k+1, go to 2.
Even ifa=0, that is when there are zero cost input symbols, the algorithm is still valid.


### 3.3. Maximum Entropy

A similar but simpler problem is the computation of maximum entropy. The typical form is as follows:(54)min⁡∑i=1mpilog⁡piπi,s.t∑j=1maijpj=Ei,i=1,…,r,p1+⋯+pm=1,p1≥0,…,pm≥0,
where πi are known, and Ei can be computed by sample data, i=1,…,r. The GIS algorithm [26–28] is one of the most common algorithms for (54), which is only suited to linear constraints. In spite of that there is a good discussion on the computation in [27]; however, we cannot ensure that GIS is a polynomial time algorithm.In [25] a maximum entropy problem with nonlinear constraints is considered. It is as follows:(55)min⁡∑i=1mpilog⁡piπi,s.t∑i=1mpiai-∑i,jσijpipj≥b,p1+⋯+pm=1,p1≥0,…,pm≥0,
where (σij)m×m is a positive definite matrix.We can solve (54) and (55) in polynomial time. In fact, for (54), we add logarithmic barriers -logpi in its objective function, for (55), we add logarithmic barrier -logpi and(56)-log⁡(∑i=1mpiai-∑i,jσijpipj-b)
in its objective function, and the problems can be transformed into (3). It is good to know that the function in (56) is self-concordant. Therefore, the path following method is valid.


## 3.1. Channel Capacity without Constrains

LetX={x1,…,xm} and Y={y1,…,yn} be the input and output alphabets, respectively. Let p(x)=(p(x1),…,p(xm)) and q(y)=(q(y1),…,q(yn)) be the distributions of X and Y, respectively. Let T=(p(yj∣xi))m×n=(pij)m×nbe a transition matrix from X to Y. Hence,(14)q(yj)=p(x1)p1j+⋯+p(xm)pmj,j=1,…,n.IfI(X;Y) is the mutual information between X and Y, then(15)I(X;Y)=∑i,jp(xi)p(yj∣xi)log⁡p(yj∣xi)q(yj)(16)=∑i,jp(xi)p(yj∣xi)log⁡p(xi∣yj)p(xi).Channel capacity is defined as follows:(17)C=max⁡p(x)I(X;Y).
Equation (16) is a function of p(x) and p(x|y) and we denote it by I(p(x);p(x|y)).The Arimoto-Blahut algorithm utilized (16) as follows. In (16), both p(x) and p(x|y) are unknown, but given p(x), (16) can get its maximum at(18)p(xi∣yj)=p(xi)p(yj∣xi)∑k=1mp(xk)p(yj∣xk),i=1,…,m.j=1,…,n,
and given p(x|y), (16) can get its maximum at(19)p(xi)=∏jp(xi∣yj)p(yj∣xi)∑i=1m∏jp(xi∣yj)p(yj∣xi),i=1,…,m.
If p0(x)=(1/m,…,1/m), then the channel capacity can be approximated by computing p0(x|y),p1(x),p1(x|y),… through (18) and (19), alternatively. Let C(N+1,N) be the mutual information I(pN+1(x);pN(x|y)), then,(20)C-C(N+1,N)≤log⁡mN.In order to get an algorithm based on the self-concordant function, we utilize (15). Let(21)ci=-∑j=1np(yj∣xi)log⁡p(yj∣xi).
By (15), the channel capacity can be expressed as follows:(22)min⁡∑i=1mcipi+∑j=1nqjlog⁡⁡qj,s.tqj=∑i=1mpip(yj∣xi),j=1,…,n,p1+⋯+pm=1,p1>0,…,pm>0.
It is a standard convex optimization problem. Unfortunately, the objective function of (22) is not self-concordant, even when the logarithmic barrier −log pi is added.It seems that we should addqj>0 as constraints for (22) since there are logqj in the objective function of (22). But the qj are not independent variables. It can be computed by (14) and maintain positive as long as all the pi>0,i=1,2,…,m, hold. Nevertheless, we still add the logarithmic barrier −log qj in the objective of (22) to get a self-concordant objective function.Fort>0, consider the convex optimization problem as follows:(23)min⁡t∑i=1mcipi+t∑j=1nqjlog⁡qj-∑i=1mlog⁡pi-∑j=1nlog⁡qj,s.tqj=∑i=1mpip(yj∣xi),j=1,…,n,p1+⋯+pm=1.
In order to show the self-concordance of the objective function of (23), we need Proposition 6.Proposition 6.
For anyt>0,f(x)=txlogx-logx is a self-concordant function on x>0.Proof.
We have(24)f′(x)=tlog⁡x+t-1x,f′′(x)=tx+1x2,|f′′′(x)|=tx2+2x3.
Thus
(25)|f′′′(x)|f′′(x)3/2=(tx+2)(tx+1)3/2≤2.By Propositions3–6, the objective function of (23) is self-concordant, so we can solve it by the path following method, and the number of Newton iterations is bounded by (13).Equation (13) shows that solving (23) by the path following method is a polynomial time algorithm. In addition, a main advantage is that the algorithm can deal with constrains very easily.


## 3.2. Channel Capacity with Constrains and Channel Capacity Per Unit Cost

Consider the channel capacity with constraints as follows:(26)max⁡I(X;Y),s.t∑j=1maijpj≤Ei,i=1,…,r,p1+⋯+pm=1,p1≥0,…,pm≥0.
For r=2, [5] gave some algorithms based on the Arimoto-Blahut method. Because one cannot eliminate the Lagrange multiplier, the algorithms in [5] cannot get the optimal solution for (26).By adding logarithmic barriers for the inequality constraints in the objective of (26) as follows, we can get a polynomial time algorithm for (26):(27)min⁡t∑i=1mcipi+t∑j=1n∑i=1mpip(yj∣xi)log⁡∑i=1mpip(yj∣xi)-∑i=1mlog⁡pi-∑j=1nlog⁡∑i=1mpip(yj∣xi)-∑k=1rlog⁡(Ek-∑i=1mpiari),s.tp1+⋯+pm=1.
The objective function of (27) is self-concordant, so we can solve (27) by the path-following method.Channel capacity per unit cost is one of the typical problems of channel capacity with constraints [15, 20–22]. By [15], channel capacity per unit cost can be computed by(28)C=sup⁡β>0C(β)β,
where C(β) is the capacity-cost function, which is the solution of the following constrained problem [12]:(29)max⁡I(X;Y),s.t∑j=1majpj≤β,p1+⋯+pm=1,p1≥0,…,pm≥0.
Equation (29) is a special case of (26) (r=1). Here we discuss the channel capacity per unit cost in a more general manner.Theorem 7.
LetD⊆Rn be a convex set. Let F(x) be a concave function on D and G(x) a convex function on D. Let a=inf{G(x)|x∈D},b=sup{G(x)|x∈D}, and suppose a≥0. Let C(β)=sup{F(x)|G(x)≤β}. Then,(1)
C(β)is a concave increasing function,(2)
let(30)ρ(β)=C(β)β,
then ρ(β) is either a unimodal function or a monotone function. 												(31)sup⁡DF(x)G(x)=sup⁡a<β<bC(β)β.Proof.
(1) Letβ1>β2, then {x|G(x)≤β2}⊆{x|G(x)≤β1}. Therefore, C(β2)≤C(β1), For any ε>0, there are x1 such that G(x1)≤β1 and C(β1)=F(x1)+ε, and there are x2 such that G(x2)≤β2 and C(β2)=F(x2)+ε. Since G(x) is convex and F(x) is concave, we get
(32)G(αx1+(1-α)x2)≤αG(x1)+(1-α)G(x2)≤αβ1+(1-α)β2,αC(β1)+(1-α)C(β2)=αF(x1)+(1-α)F(x2)+ε≤F(αx1+(1-α)x2)+ε≤C(αβ1+(1-α)β2)+ε.
Since ε is arbitrary, C(β) is concave.
(2) Without loss of generality, we can assumeC(β) is differentiable. If not, we can approach it by a differentiable function to any precision [23]. Suppose there is a β0∈(a,b), such that
(33)ρ′(β)|β=β0=0,thatisC′(β0)β0-C(β0)=0.
Letβ1>β0; since C(β) is concave, we have C′(β1)≤C′(β0). Let
(34)ρ′(β)|β=β0=0,thatisC′(β0)β0-C(β0)=0.
It is the equation of a straight line that has slope C′(β1) and passes through the (β0,C(β0)). Since C(β) is increasing and concave, we get y1(β1)≤C(β1) and
(35)C′(β1)-C(β1)β1≤C′(β1)-y1(β1)β1=C′(β1)β0β1-C(β0)β1≤C′(β0)β0β1-C(β0)β1=0,
that is, ρ′(β1)≤0.
Letβ2<β0, since C(β) is concave, C′(β2)≥C′(β0). Let
(36)y0(β)=C′(β0)(β-β0)+C(β0).
It is the equation of a straight line that has slope C′(β0) and passes through (β0,C(β0)). Since C(β) is increasing and concave, we get y0(β2)≥C(β2) and
(37)C′(β2)-C(β2)β2≥C′(β2)-y0(β2)β2≥C′(β0)-y0(β2)β2=C′(β0)β0β2-C(β0)β2=0,
that is, ρ′(β2)≥0.
Therefore,β0 is a maximum point of function ρ(β)=C(β)/β.
(3) If(38)sup⁡DF(x)G(x)=+∞,
then for any N>0, there are x0, such that
(39)F(x0)G(x0)>N.
Let G(x0)=β0, then C(β0)=sup{F(x)|G(x)≤β0}≥F(x0), so
(40)sup⁡a<β<bC(β)β≥C(β0)β0≥F(x0)G(x0)>N.
Since N is arbitrary, we get
(41)sup⁡a<β<bC(β)β=+∞.
In turn, if
(42)sup⁡a<β<bC(β)β=+∞,
then for any N>0, there are β0, a<β0<b, such that
(43)C(β0)β0>N.
Since
(44)sup⁡DF(x)G(x)≥sup⁡G(x)≤β0F(x)G(x)≥sup⁡G(x)≤β0F(x)β0=C(β0)β0>N
so
(45)sup⁡DF(x)G(x)=+∞.
Therefore, if
(46)sup⁡DF(x)G(x)<+∞
is valid, then so is
(47)sup⁡a<β<bC(β)β<+∞
and vice versa.
Suppose (46) is valid, for any ε>0, there are x0∈D, such that
(48)F(x0)G(x0)+ε=sup⁡DF(x)G(x).
Let G(x0)=β0, then
(49)sup⁡DF(x)G(x)=F(x0)G(x0)+ε≤C(β0)β0+ε≤sup⁡a<β<bC(β)β+ε.
Since ε is arbitrary, we get
(50)sup⁡DF(x)G(x)≤sup⁡a<β<bC(β)β.
On the other hand, since (47) is valid too, for any ε>0, there are a<β0<b, such that
(51)C(β0)β0+ε=sup⁡a<β<bC(β)β.
Let x0∈D, such that G(x0)≤β0 and F(x0)+εβ0=sup{F(x)|G(x)≤β0}=C(β0), then
(52)sup⁡DF(x)G(x)≥sup⁡G(x)≤β0F(x)G(x)≥F(x0)β0=sup⁡a<β<bC(β)β-2ε.
Since ε is arbitrary, we get
(53)sup⁡DF(x)G(x)≥sup⁡a<β<bC(β)β.IfF(x)=I(X,Y) is the average mutual information defined by (15) and G(x)=E[b(X)], then (3) of Theorem 7 is Theorem  2 in [15].In addition to average mutual information, entropy is another frequently used function in applications [24]. It is obvious that Theorem 7 is valid when F(x) is an entropy function (as a function of the distribution of input symbols). It is useful to point out that Theorem 7 is still valid when G(x) is a positive definite quadratic function [25].Making use of (2) of Theorem7, we can write an algorithm to compute the sup(F(x)/G(x))as follows. The algorithm is based on the 0.618 method used to locate the optimal point of a single peak function.Algorithm 8 (0.618 method).
Leta=inf{G(x)|x∈D},b=sup{G(x)|x∈D}, and assume a≥0.
(1) Given search interval[a,b], and tolerance ε>0. compute λ and μ as follows:
λ=a+0.382(b-a),μ=a+0.618(b-a); evaluate C(λ)/λ and C(μ)/μ. Let k=1.
(2) IfC(λ)/λ>C(μ)/μ, then go to 3, else go to 4.
(3) Ifb-λ≤ε, then output μ and C(μ)/μ, else, let a=λ,λ=μ,C(λ)/λ=C(μ)/μ,μ=a+0.618(b-a); evaluate C(μ)/μ, go to 5.
(4) Ifμ-a≤ε, then output λ and C(λ)/λ, else, let b=μ,μ=λ,C(μ)/μ=C(λ)/λ,λ=a+0.382(b-a); evaluate C(λ)/λ, go to 5.
(5)k=k+1, go to 2.
Even ifa=0, that is when there are zero cost input symbols, the algorithm is still valid.


## 3.3. Maximum Entropy

A similar but simpler problem is the computation of maximum entropy. The typical form is as follows:(54)min⁡∑i=1mpilog⁡piπi,s.t∑j=1maijpj=Ei,i=1,…,r,p1+⋯+pm=1,p1≥0,…,pm≥0,
where πi are known, and Ei can be computed by sample data, i=1,…,r. The GIS algorithm [26–28] is one of the most common algorithms for (54), which is only suited to linear constraints. In spite of that there is a good discussion on the computation in [27]; however, we cannot ensure that GIS is a polynomial time algorithm.In [25] a maximum entropy problem with nonlinear constraints is considered. It is as follows:(55)min⁡∑i=1mpilog⁡piπi,s.t∑i=1mpiai-∑i,jσijpipj≥b,p1+⋯+pm=1,p1≥0,…,pm≥0,
where (σij)m×m is a positive definite matrix.We can solve (54) and (55) in polynomial time. In fact, for (54), we add logarithmic barriers -logpi in its objective function, for (55), we add logarithmic barrier -logpi and(56)-log⁡(∑i=1mpiai-∑i,jσijpipj-b)
in its objective function, and the problems can be transformed into (3). It is good to know that the function in (56) is self-concordant. Therefore, the path following method is valid.


## 4. An Example

Let transmission matrixP be as follows: (57)P=(0.90.010.010.010.010.010.010.010.010.0201000000000.010.020.010.90.010.010.010.010.010.010.500.500000000.10.10.10.10.10.10.10.10.10.100.10.100.10.100.10.50000000.40.10.100.40.10.20.30.4000000000.200.200.20.200.200.10.500.200.10.100).Letε=0.00001. Suppose that the initial distribution is discrete uniform, with penalty factor μ=2 and initial penalty parameter t0=10. Iterating 70 times by the path-following method, we get channel capacity C=1.370656 and the optimal distribution:(0.0893518,0.2414400,0.1743530,0.1548483,0.0000003,0.0811754,0.1946234,0.0000004,0.0642042,0.0000080).With the Arimoto-Blahut algorithm, iterate 180 times, we get the same results.Letb=(1,2,5,7,6,8,3,9,3,4)T be cost vector of the input symbols. Let G(p)=bT·p be the expect cost of the input symbols. The curve of C(β)/β is shown in Figure 1; By Algorithm 8, the channel capacity per unit cost is 0.5801, and it is attained at β=1.8325. The search interval is [1, 5.8]. The number of iterations of the path-following is 17. The number of iterations of the 0.618 method is 19. The number of iterations of Newton’s algorithm is 5.Figure 1
Curve ofC(β)/β, where G(x) is the expected cost of the input symbol.LetG(x)=b′x+x′Dx be a positive definite quadratic function, where(58)D=(4.27933.63903.32993.15052.59393.29773.17773.01313.22074.71073.63904.45363.68853.66442.25663.77253.39123.16152.16363.73113.32993.68854.84362.88522.96073.99941.95092.65992.76744.19903.15053.66442.88524.10042.27732.89203.84542.61062.67863.30132.59392.25662.96072.27733.63622.27642.04941.38582.49843.08553.29773.77253.99942.89202.27645.41582.92413.21133.05803.08603.17773.39121.95093.84542.04942.92416.02923.86793.23173.45823.01313.16152.65992.61061.38583.21133.86794.92803.42043.52313.22072.16362.76742.67862.49843.05803.23173.42045.18563.65434.71073.73114.19903.30133.08553.08603.45823.52313.65436.5953).The curve of C(β)/β is shown in Figure 2. The maximum is 0.2156, and it is attained at β=5.9093. The number of iterations of the path-following is 17. The number of iterations of the 0.618 method is 19. The number of iterations of Newton’s algorithm is 4.Figure 2
Curve ofC(β)/β, where G(x) is a positive definite quadratic function.


## 5. Conclusion

By means of self-concordant function theory, the computation of channel capacity, especially, when there are constraints and when constraints are nonlinear, becomes very simple. When the numerator is a general concave function and the denominator is a general convex function, the formula about channel capacity per unit cost is still valid. Furthermore, the functionC(β)/β is single peak, hence we can get some new algorithms for channel capacity per unit cost.



---

*Source: 318946-2012-02-26.xml*