Blahut–Arimoto Algorithms for Inner and Outer Bounds on Capacity Regions of Broadcast Channels ^†

Abstract

The celebrated Blahut–Arimoto algorithm computes the capacity of a discrete memoryless point-to-point channel by alternately maximizing the objective function of a maximization problem. This algorithm has been applied to degraded broadcast channels, in which the supporting hyperplanes of the capacity region are again cast as maximization problems. In this work, we consider general broadcast channels and extend this algorithm to compute inner and outer bounds on the capacity regions. Our main contributions are as follows: first, we show that the optimization problems are max–min problems and that the exchange of minimum and maximum holds; second, we design Blahut–Arimoto algorithms for the maximization part and gradient descent algorithms for the minimization part; third, we provide convergence analysis for both parts. Numerical experiments validate the effectiveness of our algorithms.

1. Introduction

In 1972, Cover [1] introduced the two-receiver discrete memoryless broadcast channel $p(y,z|x)$ to model a system of downlink communication in which X is the sender and $(Y,Z)$ are the receivers. In the same paper, he proposed a coding scheme which resulted in the superposition coding inner bound (SCIB). It turns out that the SCIB is indeed the capacity region for two-receiver broadcast channels in which the receivers are comparable in the following partial orders: degraded [2], less noisy [3], and more capable [3]. However, for a general broadcast channel the single-letter capacity region remains open.

To characterize the capacity region of a broadcast channel, a standard approach is to show that one inner bound matches another outer bound. Currently, the best inner bound for general broadcast channels is Marton’s inner bound (MIB) [4], while the UV outer bound (UVOB) [5] was the best outer bound until recently, when a better one called J version outer bound was proposed in [6].

The evaluation of inner and outer bounds is critical in the following aspects: (1) the evaluation of an inner bound usually results in an optimal input distribution which can help in the design of practical coding schemes; (2) the identification of the capacity region of a particular broadcast channel through the comparison of one inner bound and one outer bound relies on the evaluation of these two bounds; and (3) when claiming to establish a new bound, it is necessary to show that the new bound strictly improves on old ones through the evaluation of bounds on a particular channel.

Remark 1.

This is the full version of conference papers accepted by ISIT 2022 and 2023 [7,8].

However, this evaluation is usually difficult due to its non-convexity [9]. To alleviate this issue, there exist a number of generic optimization algorithms, such as interior point [10], active set [10], and sequential quadratic programming [11]. However, efficient algorithms should use the domain knowledge of information theory as well; from this viewpoint, we consider the Blahut–Arimoto (BA) algorithm, which is specially customized for information theory.

The original BA algorithm was independently developed by Blahut [12] and Arimoto [13] to calculate the channel capacity $C={max}_{q\left(x\right)}I(X;Y)$ for a general point-to-point channel $p\left(y\right|x)$. The algorithm transforms the original maximization problem into an alternating maximization problem:

$$\begin{array}{c}\hfill \underset{q\left(x\right)}{max}\sum _{x,y}q\left(x\right)p\left(y\right|x)ln\frac{q\left(x\right|y)}{q\left(x\right)},\to \underset{q\left(x\right),Q\left(x\right|y)}{max}\sum _{x,y}q\left(x\right)p\left(y\right|x)ln\frac{Q\left(x\right|y)}{q\left(x\right)}\end{array}$$

where the updating formulae are explicit within each iteration.

There have been numerous extensions of the BA algorithm to various scenarios in information theory. For example, [14] applied the BA algorithm to compute the sum rate of the multiple-access channel. Later, using the idea from the BA algorithm, the whole capacity region of the multiple-access channel was formulated in [15] as a rank-one constrained problem and solved by relaxation methods. It is beyond the scope of this paper to list all of these references. Instead, we discuss those papers closely related to computing the bounds on capacity regions of broadcast channels.

In [16], the authors considered the capacity region of a degraded broadcast channel $p(y,z|x)$, where receiver Z is a degraded version of Y.

In this scenario, the capacity region of the rate pairs $(R_Y,R_Z)$ is known, and can be achieved by the simplified version of superposition coding. The supporting hyperplanes can be characterized as

$$\begin{array}{c}\hfill \theta R_Y+(1-\theta )R_Z=\underset{q(u,x)}{max}\theta I(X;Y|U)+(1-\theta )I(U;Z).\end{array}$$

Using a similar idea to that of the BA algorithm, the authors designed an algorithm to alternatively maximize the objective function.

The method in [16] is directly applicable to less noisy broadcast channels, as the characterization of the capacity region is the same as that of the degraded case. However, this equivalence no longer holds for the more capable case, as this time the value of the supporting hyperplane $\theta R_Y+(1-\theta )R_Z$ is characterized as a max–min optimization problem (e.g., see Equation (14)). As a mater of fact, the supporting hyperplanes of the above-mentioned bounds, that is, SCIB, MIB, and UVOB, are all of the max–min form. The main issue is that the minimization part is inside the maximization part, which prevents the application of the BA algorithm to the whole problem.

The algorithms for calculating inner bounds and outer bounds for general broadcast channels are very limited. The authors of [17] considered MIB (see Section 3.2) and designed a BA algorithm to compute the sum rate $R_Y+R_Z$ of the simplified version, where the auxiliary random variable $W=\mathsf{\varnothing }$. The objective function,

$$\begin{array}{c}\hfill \sum _{u,v,x,y,z}q(u,v)q\left(x\right|u,v)p(y,z|x)ln\frac{Q\left(u\right|y\left)Q\right(v\left|z\right)}{q(u,v)}\end{array}$$

is convex in $q\left(x\right|u,v)$, which means that the maximum input X is a function of $(U,V)$. Noticing this, the authors performed optimization over all fixed mappings $q\left(x\right|u,v)$. However, discarding W can result in a strictly smaller sum rate [18], making it is necessary to consider the complete version of MIB.

In this paper, we seek to design BA algorithms for general broadcast channels in order to compute the following inner and outer bounds: SCIB, MIB, and UVOB. The key difference here is that the optimization problems are max–min problems, rather than only containing a maximization part. In

Table 1. Comparison of typical scenarios related to the BA algorithm.

Channel	Reference	Objective	Form	Algorithm
point-to-point	[12,13]	capacity	max	BA
multiple access	[14]	sum-rate	max	BA
multiple access	[15]	inner/outer bounds	max	relaxation
degraded broadcast	[16]	capacity region	max	BA
general broadcast	this paper	inner/outer bounds	max–min	BA + gradient

, we provide an intuitive comparison of related references.

The notation we use is as follows. p denotes a fixed (conditional) probability distribution such as $p(y,z|x)$, while q and Q are used for (conditional) probabilities that are changeable. Calligraphic letters such as $\mathcal{S}$ are used to denote sets. The use of square brackets in the function $f\left[g\right]$ means that f is specified by the variable g; $\overline{\theta }$ denotes $\overline{\theta }:=1-\theta$; and unless otherwise specified, we use the natural logarithm. To make the mathematical expressions more concise, we use the following abbreviations of the Kullback–Leibler divergences:

$$\begin{array}{cc}\hfill D_Y\left(p\right|\left|q\right):=& \sum _{y}p\left(y\right)ln\frac{p\left(y\right)}{q\left(y\right)},\hfill \\ \hfill D_{Y|x}\left(p\right|\left|q\right):=& \sum _{y}p\left(y\right|x)ln\frac{p\left(y\right|x)}{q\left(y\right|x)},\hfill \\ \hfill D_{Y|X}\left(p\right|\left|q\right):=& \sum _{x}p\left(x\right)·D_{Y|x}\left(p\right|\left|q\right).\hfill \end{array}$$

The organization is as follows. First, in Section 2 we introduce the necessary background on the BA algorithm and its extension in [16]. Then, in Section 3 we extend the BA algorithm to the evaluation of SCIB, MIB, and UVOB. Convergence analyses of these algorithms are presented in Section 4. Finally, in Section 5 we perform numerical experiments to validate the effectiveness and efficiency of our algorithms.

2. Mathematical Background of Blahut–Arimoto Algorithms

We first introduce the standard BA algorithm in Section 2.1, as we will rely on several of its properties later. Then, in Section 2.2 we discuss why the method in [16] cannot be applied to general broadcast channels.

2.1. Blahut–Arimoto Algorithm for Point-to-Point Channel

For a point-to-point channel $p\left(y\right|x)$, the capacity C is the maximum of the mutual information $C={max}_{q\left(x\right)}I(X;Y)={max}_{q\left(x\right)}C\left(q\right)$, where

$$\begin{array}{c}\hfill C\left(q\right)=H\left(X\right)-H\left(X\right|Y)=\sum _{x,y}q\left(x\right)p\left(y\right|x)ln\frac{q\left(x\right|y)}{q\left(x\right)}.\end{array}$$

(1)

By replacing $q\left(x\right|y)$ with a free variable $Q\left(x\right|y)$, the BA algorithm performs alternating maximization ${max}_{q,Q}C(q,Q)$, where

$$\begin{array}{c}\hfill C(q,Q)=\sum _{x,y}q\left(x\right)p\left(y\right|x)ln\frac{Q\left(x\right|y)}{q\left(x\right)}.\end{array}$$

(2)

Notice here that we abuse the notation of $C(·)$, which should not cause confusion in general.

The above objective function can be reformulated as follows:

$$\begin{array}{cc}\hfill C(q,Q)& =\sum _{x}q\left(x\right)(d\left[Q\right]\left(x\right)-lnq\left(x\right)),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \mathrm{where}d\left[Q\right]\left(x\right)& =\sum _{y}p\left(y\right|x)lnQ\left(x\right|y).\hfill \end{array}$$

(4)

We call this the basic form. For different scenarios, BA algorithms mainly differ in the distribution $q(·)$ and function $d\left[Q\right](·)$.

The following theorem (see proof in Appendix A) provides the explicit formulae for the maximum Q given q (denoted as $Q\left[q\right]$) and maximum q given Q (denoted as $q\left[Q\right]$).

Theorem 1.

The following properties hold for the problem ${max}_{q,Q}C(q,Q)$.

• Given a fixed q, $C(q,Q)$ is concave in Q, and the maximum point $Q\left[q\right]$ is induced by the input and the channel,

  $$\begin{array}{c}\hfill Q\left[q\right]\left(x\right|y)=q\left(x\right|y)=\frac{q\left(x\right)p\left(y\right|x)}{{\sum }_{x^{\prime }}q\left(x^{\prime }\right)p\left(y\right|x^{\prime })}.\end{array}$$

  (5)
  Further, the function values satisfy

  $$\begin{array}{cc}\hfill C\left(q\right)& =C(q,Q[q\left]\right)\ge C(q,Q).\hfill \end{array}$$

  (6)
• Given a fixed Q, $C(q,Q)$ is concave in q, and the maximum point $q\left[Q\right]$ is obtained by the Lagrangian,

  $$\begin{array}{c}\hfill q\left[Q\right]=\frac{exp\left\{d\right[Q\left]\right(x\left)\right\}}{{\sum }_{x^{\prime }}exp\left\{d\left[Q\right]\left(x^{\prime }\right)\right\}}.\end{array}$$

  (7)
  Further, evaluation of the function value results in

  $$\begin{array}{cc}\hfill C\left(q\right[Q],Q)& =ln\sum _{x^{\prime }}exp\left\{d\left[Q\right]\left(x^{\prime }\right)\right\}=d\left[Q\right]\left(x\right)-lnq\left[Q\right]\left(x\right),\forall x.\hfill \end{array}$$

  (8)

Starting from an initial $q_0\left(x\right)>0$, the BA algorithm performs alternating maximization and produces a sequence of points

$$\begin{array}{c}\hfill q_0\to Q_1\to q_1\to Q_2\to \dots \to q_n\to Q_{n+1}\to \dots \end{array}$$

where

$$\begin{array}{cc}\hfill Q_n& =Q\left[q_{n-1}\right],\hfill \\ \hfill q_n& =q\left[Q_n\right]\hfill \end{array}$$

according to Equations (5) and (7), respectively.

The criterion for stopping the iterations is based on the following result [12,13].

Proposition 1

(Proposition 1 in [13], Theorem 2 in [12]). It is the case that $q_{\ast }$ maximizes $C\left(q\right)$ if and only if the following holds for $Q\left[q_{\ast }\right]$ and some scalar D:

$$\begin{array}{c}\hfill d\left[Q\left[q_{\ast }\right]\right]\left(x\right)-ln{q}_{\ast }\left(x\right)\left\{\begin{array}{cc}=D,\hfill & q_{\ast }\left(x\right)>0,\hfill \\ \le D,\hfill & q_{\ast }\left(x\right)=0.\hfill \end{array}\right.\end{array}$$

Remark 2.

It should be mentioned that in order to avoid infinity minus infinity for those $q_{\ast }\left(x\right)=0$ in the above proposition, the following equivalent formulae can be used:

$$\begin{array}{cc}\hfill d\left[Q\right[q\left]\right]\left(x\right)-lnq\left(x\right)& =\sum _{y}p\left(y\right|x)ln\frac{q\left(x\right|y)}{q\left(x\right)}=\sum _{y}p\left(y\right|x)ln\frac{p\left(y\right|x)}{q\left(y\right)}.\hfill \end{array}$$

(9)

This is not be an issue in the BA algorithm, as $q_n\left(x\right)>0$ according to Equation (7).

Thus, at the end of the n-th step, if the difference

$$\begin{array}{cc}\hfill & \underset{x}{max}\{d\left[Q\left[q_{n-1}\right]\right]\left(x\right)-ln{q}_{n-1}\left(x\right)\}-C(q_n,Q_n)\hfill \\ \hfill & =\underset{x}{max}\{d\left[Q_n\right]\left(x\right)-ln{q}_{n-1}\left(x\right)\}-ln\sum _{x^{\prime }}exp\left\{d\left[Q_n\right]\left(x^{\prime }\right)\right\}\hfill \end{array}$$

is small enough then the iteration is stopped.

Summarizing the above details, we arrive at the BA algorithm depicted in Algorithm 1. The convergence of the resulted sequence of $C(q_n,Q_n)$ is characterized in the following theorem (see proof in Appendix A).

Theorem 2

(Theorem 1 in [13], Theorem 3 in [12]). If $p_0>0$, then the value $C(q_n,Q_n)$ converges monotonically from below to the capacity C.

Algorithm 1: Computing channel capacity
Input: $p\left(y\right|x)$, maximum iterations N, threshold $ϵ>0$;
Initialization: $q_0\left(x\right)>0$, ${ϵ}_0>ϵ$, $n=0$;
while $n<N$ and ${ϵ}_n>ϵ$ do
		$\phantom{(}$ $\phantom{(}$ $n\leftarrow n+1$;
		$Q_n=Q\left[q_{n-1}\right]$ using Equation (5);
		$q_n=q\left[Q_n\right]$ using Equation (7);
		$C(q_n,Q_n)=ln{\sum }_{x}exp\left\{d\left[Q_n\right]\left(x\right)\right\}$ using Equations (8) and (4);
		${ϵ}_n={max}_x\{d\left[Q_n\right]\left(x\right)-ln{q}_{n-1}\left(x\right)\}-C(q_n,Q_n)$ using Equation (4);
end
Output: $q_n\left(x\right)$, $Q_n\left(x\right|y)$, $C(q_n,Q_n)$

2.2. Blahut–Arimoto Algorithm for Degraded Broadcast Channel

In [16], the authors considered the capacity region of the degraded broadcast channel. The original objective function for the value of the supporting hyperplane $\theta R_Y+\overline{\theta }{R}_Z$ (where $\theta <\frac{1}{2}$) is:

$$\begin{array}{cc}\hfill & F\left(q\right(u,x\left)\right)\hfill \\ \hfill & =\theta I(X;Y|U)+\overline{\theta }I(U;Z)\hfill \\ \hfill & =\overline{\theta }(I(X;Y|U)+I(U;Z))-(\overline{\theta }-\theta )I(X;Y|U)\hfill \end{array}$$

(10)

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}=\overline{\theta }(H(U,X)-H\left(X\right|Y,U)-H\left(U\right|Z))-(\overline{\theta }-\theta )(H\left(Y\right|U)-H\left(Y\right|X)).\hfill \end{array}$$

(11)

Similar to the BA algorithm, the new objective function is

$$\begin{array}{cc}\hfill F(q,Q)& =\overline{\theta }·\sum _{u,x}q(u,x)(d\left[Q\right](u,x)-lnq(u,x)),\hfill \end{array}$$

(12)

where

$$\begin{array}{cc}\hfill d\left[Q\right](u,x)& =\sum _{y,z}p(y,z|x)\left(lnQ\left(x\right|y,u)Q\left(u\right|z)+\frac{\overline{\theta }-\theta }{\overline{\theta }}ln\frac{Q\left(y\right|u)}{p\left(y\right|x)}\right).\hfill \end{array}$$

(13)

Then, the authors designed an extended BA algorithm that alternately maximizes $F(q,Q)$ and analysed the convergence.

However, the above method does not generalize to allow for evaluating the capacity bounds of general broadcast channels. The main reason can be summarized in just one sentence: the minimum of the expectation is, in general, greater than the expectation of the minimum. Taking SCIB as an example, using the representation $\mathcal{A}(U,X)$ in Lemma 1, the supporting hyperplane is

$$\begin{array}{cc}\hfill \theta R_Y+\overline{\theta }{R}_Z& =\underset{q(u,x)}{max}\theta I(X;Y|U)+\overline{\theta }·min\{I(U;Z),I(U;Y)\}\hfill \\ \hfill & =\underset{q(u,x)}{max}(\overline{\theta }(H(U,X)-H\left(X\right|Y,U)+min\{-H\left(U\right|Z),-H\left(U\right|Y)\})\hfill \\ \hfill & -(\overline{\theta }-\theta )(H\left(Y\right|U)-H\left(Y\right|X))).\hfill \end{array}$$

A direct extension might try to reformulate the last expression in the form of Equation (12) by letting

$$\begin{array}{cc}\hfill d\left[Q\right](u,x)& =\sum _{y,z}p(y,z|x)\left(lnQ\left(x\right|y,u)+min\{lnQ\left(u\right|z),lnQ\left(u\right|y)\}+\frac{\overline{\theta }-\theta }{\overline{\theta }}ln\frac{Q\left(y\right|u)}{p\left(y\right|x)}\right),\hfill \end{array}$$

however, this is not an equivalent reformulation, as

$$\begin{array}{c}\hfill min\{-H\left(U\right|Z),-H\left(U\right|Y)\}\ge \sum _{u,x}q(u,x)·\sum _{y,z}p(y,z|x)·min\{lnQ\left(u\right|z),lnQ\left(u\right|y)\}.\end{array}$$

In the following section, we first use the fact that $min\{f,g\}={min}_{\alpha \in [0,1]}\{\alpha f+\overline{\alpha }g\}$ to parameterize the minimum, then show the max–min exchanges so that we can apply the BA algorithm in the maximization part.

3. Blahut–Arimoto Algorithms for Capacity Bounds of Broadcast Channels

In this section, we first introduce two inner bounds and one outer bound on the capacity region of the broadcast channel. We characterize their supporting hyperplanes as max–min problems and show that the maximum and minimum can be exchanged. Then, we design BA algorithms for the maximization parts and gradient descent algorithms for the minimization parts.

3.1. Superposition Coding Inner Bound

The superposition coding inner bound was proposed by Cover in [1], which corresponds to region $\mathcal{B}$ in the following lemma (see proof in Appendix B). This region actually has three equivalent characterizations.

Lemma 1

(folklore). The following regions $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$ are equivalent characterizations of the superposition coding inner bound:

$$\begin{array}{cc}\hfill \mathcal{A}& :=\bigcup _{q(u,x)}\mathcal{A}(U,X):=\left\{\begin{array}{cc}\hfill R_Y& \le I(X;Y|U),\hfill \\ \hfill R_Z& \le min\left\{I\right(U;Z),I(U;Y\left)\right\},\hfill \end{array}\right.\hfill \\ \hfill \mathcal{B}& :=\bigcup _{q(u,x)}\mathcal{B}(U,X):=\left\{\begin{array}{cc}\hfill R_Y& \le I(X;Y|U),\hfill \\ \hfill R_Z& \le I(U;Z),\hfill \\ \hfill R_Y+R_Z& \le I(X;Y),\hfill \end{array}\right.\hfill \\ \hfill \mathcal{C}& :=\bigcup _{q(u,x)}\mathcal{C}(U,X):=\left\{\begin{array}{cc}\hfill R_Z& \le I(U;Z),\hfill \\ \hfill R_Y+R_Z& \le min\left\{I\right(X;Y\left|U\right)+I(U;Z),I(X;Y)\}.\hfill \end{array}\right.\hfill \end{array}$$

To characterize the supporting hyperplane $\theta R_Y+\overline{\theta }{R}_Z=F$, we choose to use the representation $\mathcal{A}(U,X)$. It is clear that

$$\begin{array}{cc}\hfill F& =\underset{q(u,x)}{max}\theta I(X;Y|U)+\overline{\theta }·min\{I(U;Z),I(U;Y)\}\hfill \\ \hfill & =\underset{q(u,x)}{max}\theta I(X;Y|U)+\overline{\theta }·\underset{\alpha \in [0,1]}{min}\left(\alpha I(U;Z)+\overline{\alpha }I(U;Y)\right)\hfill \\ \hfill & =\underset{q(u,x)}{max}\underset{\alpha \in [0,1]}{min}\theta I(X;Y|U)+\overline{\theta }\alpha I(U;Z)+\overline{\theta }\overline{\alpha }I(U;Y).\hfill \end{array}$$

(14)

Notice here that we cannot use the BA algorithm directly, as there is a minimum inside the maximum. If we are able to swap the orders of the maximum and minimum, then we can adopt the BA algorithm in the maximization part.

To show this kind of exchange, we first introduce a Terkelsen-type min–max result in the following lemma.

Lemma 2

(Corollary 2 in Appendix A of [19]). Let ${\Lambda }_d$ be the d-dimensional simplex, i.e., ${\lambda }_i\ge 0$ and ${\sum }_{i=1}^d{\lambda }_i=1$, let $\mathcal{P}$ be a set of probability distributions $p\left(u\right)$, and let $T_i\left(p\left(u\right)\right)$, $i=1,..,d$ be a set of functions such that the set $\mathcal{T}$ defined by

$$\begin{array}{cc}\hfill \mathcal{T}& =\{(a_1,...,a_d)\in {\mathbb{R}}^d:a_i\le T_i\left(p\left(u\right)\right)\mathrm{for}\mathrm{some}p\left(u\right)\in \mathcal{P}\}\hfill \end{array}$$

is a convex set; then,

$$\underset{p\left(u\right)\in \mathcal{P}}{sup}\underset{\lambda \in {\mathsf{\Lambda }}_d}{min}\sum _i{\lambda }_i{T}_i\left(p\left(u\right)\right)=\underset{\lambda \in {\mathsf{\Lambda }}_d}{min}\underset{p\left(u\right)\in \mathcal{P}}{sup}\sum _i{\lambda }_i{T}_i\left(p\left(u\right)\right).$$

With this lemma, we are to establish the following theorem (proof in Appendix B).

Theorem 3.

The supporting hyperplane $\theta R_Y+\overline{\theta }{R}_Z=F$ of the superposition coding inner bound is as follows: if $\theta \in [\frac{1}{2},1]$, then $F={max}_{q\left(x\right)}\theta I(X;Y)$; otherwise, if $\theta \in [0,\frac{1}{2})$ then

$$\begin{array}{cc}\hfill F=min\{& \underset{\alpha \le \frac{1-2\theta }{1-\theta }}{min}\underset{q\left(x\right)}{max}\overline{\theta }\overline{\alpha }I(X;Y)+\overline{\theta }\alpha I(X;Z),\hfill \\ \hfill & \underset{\alpha >\frac{1-2\theta }{1-\theta }}{min}\underset{q(u,x)}{max}\theta I(X;Y|U)+\overline{\theta }\overline{\alpha }I(U;Y)+\overline{\theta }\alpha I(U;Z)\}.\hfill \end{array}$$

Further, it suffices to consider the cardinality size: $\left|\mathcal{U}\right|\le \left|\mathcal{X}\right|$.

For the maximization part and the nontrivial case where $\theta \in (0,\frac{1}{2})$, following the above theorem, two types of BA algorithms can be designed according to the value of $\alpha$.

When $\alpha \in (\frac{1-2\theta }{1-\theta },1]$, the original objective function is

$$\begin{array}{cc}\hfill F(\alpha ,q)& =\theta I(X;Y|U)+\overline{\theta }\overline{\alpha }I(U;Y)+\overline{\theta }\alpha I(U;Z)\hfill \\ \hfill & =\overline{\theta }(I(X;Y|U)+\overline{\alpha }I(U;Y)+\alpha I(U;Z))-(\overline{\theta }-\theta )I(X;Y|U)\hfill \\ \hfill & =\overline{\theta }(H(U,X)-H\left(X\right|Y,U)-\overline{\alpha }H\left(U\right|Y)-\alpha H\left(U\right|Z))\hfill \end{array}$$

(15)

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+(\overline{\theta }-\theta )(H\left(Y\right|X)-H\left(Y\right|U))\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\overline{\theta }·\sum _{u,x,y,z}q(u,x)p(y,z|x)\left(ln\frac{q\left(x\right|y,u)}{q(u,x)}q{\left(u\right|y)}^{\overline{\alpha }}q{\left(u\right|z)}^{\alpha }+\frac{\overline{\theta }-\theta }{\overline{\theta }}ln\frac{q\left(y\right|u)}{p\left(y\right|x)}\right).\hfill \end{array}$$

(16)

By replacing the conditional qs with free variables Qs, we have the new objective function

$$\begin{array}{cc}\hfill F(\alpha ,q,Q)& =\overline{\theta }·\sum _{u,x}q(u,x)(d\left[Q\right](u,x)-lnq(u,x)),\hfill \end{array}$$

(17)

where

$$\begin{array}{cc}\hfill d\left[Q\right](u,x)& =\sum _{y,z}p(y,z|x)\left(lnQ\left(x\right|y,u)Q{\left(u\right|y)}^{\overline{\alpha }}Q{\left(u\right|z)}^{\alpha }+\frac{\overline{\theta }-\theta }{\overline{\theta }}ln\frac{Q\left(y\right|u)}{p\left(y\right|x)}\right).\hfill \end{array}$$

(18)

When $\alpha \in [0,\frac{1-2\theta }{1-\theta }]$, similar to Equation (4), the new objective function is

$$\begin{array}{cc}\hfill F(\alpha ,q,Q)& =\overline{\theta }·\sum _{x}q\left(x\right)(d\left[Q\right]\left(x\right)-lnq\left(x\right)),\hfill \end{array}$$

(19)

where

$$\begin{array}{cc}\hfill d\left[Q\right]\left(x\right)& =\sum _{y,z}p(y,z|x)(\overline{\alpha }lnQ\left(x\right|y)+\alpha lnQ\left(x\right|z)).\hfill \end{array}$$

(20)

For the minimization part, it is possible to use the optimal q and Qs obtained in the maximization part to update $\alpha$. Because the values of the optimal q and Qs may vary greatly when $\alpha$ changes, we propose changing $\alpha$ locally in the neighbourhood. A candidate approach is to use the gradient descent method, as follows:

$$\begin{array}{c}\hfill {\alpha }_{k+1}={\alpha }_k-{\tau }_k·\frac{\partial }{\partial \alpha }F({\alpha }_k,q),\end{array}$$

(21)

where

$$\begin{array}{c}\hfill \frac{\partial F(\alpha ,q)}{\partial \alpha }=\left\{\begin{array}{cc}\overline{\theta }·(I(X;Z)-I(X;Y)),\hfill & \mathrm{if}0\le \alpha \le \frac{1-2\theta }{1-\theta },\hfill \\ \overline{\theta }·(I(U;Z)-I(U;Y)),\hfill & \mathrm{if}\frac{1-2\theta }{1-\theta }<\alpha \le 1.\hfill \end{array}\right.\end{array}$$

(22)

If the change in $\alpha$ is sufficiently small, it can be assumed that the optimization with respect to $\alpha$ converges and then stop the iteration.

We summarize the above procedures in Algorithm 2. Note that the updating rules for the q and Qs depend on the interval in which the value of $\alpha$ falls.

Algorithm 2: Computing the superposition coding inner bound for $\theta \in (0,\frac{1}{2})$
Input: $p(y,z|x)$, maximum iterations K, N, thresholds $\eta$, $ϵ>0$, step size $\tau >0$;
Initialization: ${\alpha }_0\in (0,1)$, $q_0(u,x)>0$, ${\eta }_{\alpha }>\eta$, $k=0$;
while $k<K$ and ${\eta }_{\alpha }>\eta$ do
		initialize ${ϵ}_q>ϵ$, $n=0$;
		while $n<N$ and ${ϵ}_q>ϵ$ do
				$n\leftarrow n+1$;
				$Q_n=Q\left[q_{n-1}\right]$ using Equation (5) similarly;
				$q_n=q\left[Q_n\right]$ using Equation (7) similarly;
				$F({\alpha }_k,q_n,Q_n)=\overline{\theta }·ln{\sum }_{u,x}exp\left\{d\left[Q_n\right]\right\}$ using Equation (20) or (18);
				${ϵ}_q=\overline{\theta }·max\{d\left[Q_n\right]-ln{q}_{n-1}\}-F({\alpha }_k,q_n,Q_n)$;
		end
		$\phantom{(}$ $k\leftarrow k+1$;
		calculate ${\alpha }_k$ using Equations (21) and (22);
		${\alpha }_k\leftarrow min\{1,max\{0,{\alpha }_k\}\}$;
		${\eta }_{\alpha }=|{\alpha }_k-{\alpha }_{k-1}|$;
		$q_0\leftarrow q_n$;
end
Output: ${\alpha }_k$, $q_n(u,x)$, $Q_n$, $F({\alpha }_{k-1},q_n,Q_n)$

Remark 3.

Given ${\alpha }_k$, according to Equation (8),

$$\begin{array}{cc}\hfill F({\alpha }_k,q_n,Q_n)& =\overline{\theta }·ln\sum _{u,x}exp\left\{d\left[Q_n\right](u,x)\right\},\hfill \\ \hfill \Rightarrow exp\left\{\frac{1}{\overline{\theta }}F({\alpha }_k,q_n,Q_n)\right\}& =\sum _{u,x}exp\left\{d\left[Q_n\right](u,x)\right\}.\hfill \end{array}$$

According to Equation (18) or (20), ${\sum }_{u,x}exp\left\{d\left[Q_n\right]\right\}$ equals a sum of exponents, as it is a function of ${\alpha }_k$. This is a simple function; thus, we might wonder whether we can minimize $F({\alpha }_k,q_n,Q_n)$ to update ${\alpha }_k$. It turns out that this kind of global updating rule can result in an oscillating effect, as can be observed from Figure 1 in [7]. The main reason for this is that $q\left[Q\right]$ depends locally on α; therefore, it is not suitable to update α globally.

3.2. Marton’s Inner Bound

Marton’s inner bound [4] refers to the union over $q(u,v,w,x)$ (such that $(U,V,W)\to X\to (Y,Z)$ is Markov) of the non-negative rate pairs $(R_Y,R_Z)$ satisfying

$$\begin{array}{cc}\hfill R_Y& \le I(U,W;Y),\hfill \\ \hfill R_Z& \le I(V,W;Z),\hfill \\ \hfill R_Y+R_Z& \le min\left\{I\right(W;Y),I(W;Z\left)\right\}+I(U;Y|W)+I(V;Z\left|W\right)-I(U;V|W).\hfill \end{array}$$

For general broadcast channels, this is the most well known inner bound.

In the following, we characterize the supporting hyperplane $\theta R_Y+\overline{\theta }{R}_Z=M$ of MIB. Because the expressions in MIB have symmetry in Y and Z, without loss of generality we can assume that $\theta \le \overline{\theta }$, i.e., $\theta \in [0,\frac{1}{2}]$. According to [20], the supporting hyperplane is stated in the following lemma.

Lemma 3.

(Equations (2) and (5) in [20]). The supporting hyperplane $\theta R_Y+\overline{\theta }{R}_Z=M$ of Marton’s inner bound, where $\theta \in [0,\frac{1}{2}]$, is

$$\begin{array}{cc}\hfill M& =\underset{q(u,v,w,x)}{max}\underset{\alpha \in [0,1]}{min}M(\alpha ,q(u,v,w,x))\hfill \\ \hfill & =\underset{\alpha \in [0,1]}{min}\underset{q(u,v,w,x)}{max}M(\alpha ,q(u,v,w,x)),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & M(\alpha ,q)=(\overline{\theta }-\alpha \theta )I(W;Z)+\alpha \theta I(W;Y)+\overline{\theta }I(V;Z|W)+\theta I(U;Y|W)-\theta I(U;V|W).\hfill \end{array}$$

(23)

Further, it suffices to consider the following cardinalities: $\left|\mathcal{U}\right|$, $\left|\mathcal{V}\right|\le \left|\mathcal{X}\right|$, and $\left|\mathcal{W}\right|\le \left|\mathcal{X}\right|+4$.

To compute the value of this supporting hyperplane, we can reformulate $M(\alpha ,q)$ as follows:

$$\begin{array}{cc}\hfill M(\alpha ,q)& =\overline{\theta }H(U,V,W,X)-(\overline{\theta }-\alpha \theta )H\left(W\right|Z)-\alpha \theta H\left(W\right|Y)-\overline{\theta }H\left(V\right|W,Z)\hfill \\ \hfill & -\theta H\left(U\right|W,Y)-(\overline{\theta }-\theta )H\left(U\right|V,W)-\overline{\theta }H\left(X\right|U,V,W).\hfill \end{array}$$

(24)

Then, the objective function $M(\alpha ,q,Q)$ can be expressed as

$$\begin{array}{cc}\hfill M(\alpha ,q,Q)& =\overline{\theta }\sum _{u,v,w,x}q(u,v,w,x)(d\left[Q\right](u,v,w,x)-lnq(u,v,w,x)),\hfill \end{array}$$

(25)

where

$$\begin{array}{cc}\hfill d\left[Q\right]& =\sum _{y,z}p(y,z|x)(\frac{\overline{\theta }-\alpha \theta }{\overline{\theta }}lnQ\left(w\right|z)+\frac{\alpha \theta }{\overline{\theta }}lnQ\left(w\right|y)+lnQ\left(v\right|w,z)\hfill \\ \hfill & +\frac{\theta }{\overline{\theta }}lnQ\left(u\right|w,y)+\frac{\overline{\theta }-\theta }{\overline{\theta }}lnQ\left(u\right|v,w)+lnQ\left(x\right|u,v,w)).\hfill \end{array}$$

(26)

For minimization over $\alpha$, similar to Section 3.1, we update $\alpha$ along the gradient

$$\begin{array}{c}\hfill {\alpha }_{k+1}={\alpha }_k-{\tau }_k·\frac{\partial }{\partial \alpha }M({\alpha }_k,q)={\alpha }_k-{\tau }_k·\theta ·(I(W;Y)-I(W;Z)).\end{array}$$

(27)

Similar to Algorithm 2, we summarize the algorithm for MIB in Algorithm 3.

Algorithm 3: Computing Marton’s inner bound for $\theta \in (0,\frac{1}{2}]$
Input: $p(y,z|x)$, maximum iterations K, N, thresholds $\eta$, $ϵ>0$, step size $\tau >0$;
Initialization: ${\alpha }_0\in (0,1)$, $q_0(u,v,w,x)>0$, ${\eta }_{\alpha }>\eta$, $k=0$;
while $k<K$ and ${\eta }_{\alpha }>\eta$ do
		initialize ${ϵ}_q>ϵ$, $n=0$;
		while $n<N$ and ${ϵ}_q>ϵ$ do
				$n\leftarrow n+1$;
				$Q_n=Q\left[q_{n-1}\right]$ using Equation (5) similarly;
				$q_n=q\left[Q_n\right]$ using Equation (7) similarly;
				$M({\alpha }_k,q_n,Q_n)=\overline{\theta }·ln{\sum }_{u,v,w,x}exp\left\{d\left[Q_n\right]\right\}$ using Equation (26);
				${ϵ}_q=\overline{\theta }·max\{d\left[Q_n\right]-ln{q}_{n-1}\}-M({\alpha }_k,q_n,Q_n)$;
		end
		$k\leftarrow k+1$;
		calculate ${\alpha }_k$ using Equation (27);
		${\alpha }_k\leftarrow min\{1,max\{0,{\alpha }_k\}\}$;
		${\eta }_{\alpha }=|{\alpha }_k-{\alpha }_{k-1}|$;
		$q_0\leftarrow q_n$;
end
Output: ${\alpha }_k$, $q_n(u,v,w,x)$, $Q_n$, $M({\alpha }_{k-1},q_n,Q_n)$

3.3. UV Outer Bound

The UV outer bound [5] refers to the union over $q(u,v,x)$ of non-negative rate pairs $(R_Y,R_Z)$ satisfying

$$\begin{array}{cc}\hfill R_Y& \le I(U;Y),\hfill \\ \hfill R_Z& \le I(V;Z),\hfill \\ \hfill R_Y+R_Z& \le I(U;Y)+I(X;Z|U),\hfill \\ \hfill R_Y+R_Z& \le I(V;Z)+I(X;Y|V).\hfill \end{array}$$

For general broadcast channels, this was the best outer bound until [6] strictly improved upon it over an erasure Blackwell channel. The following theorem (proof in Appendix B) characterizes the supporting hyperplanes.

Theorem 4

(Claim 2 and Remark 1 in [21]). The supporting hyperplane of the UV outer bound is

$$\begin{array}{cc}\hfill \theta R_Y+\overline{\theta }{R}_Z& =\underset{\alpha ,\beta }{min}\underset{q(u,v,x)}{max}\overline{\theta }\alpha I(V;Z)+\overline{\theta }\overline{\alpha }I(X;Z|U)+\theta \beta I(U;Y)+\theta \overline{\beta }I(X;Y|V),\hfill \end{array}$$

(28)

where $\alpha ,\beta \in [0,1]$ satisfy $\overline{\theta }\alpha +\theta \beta \ge max\{\theta ,\overline{\theta }\}$. Further, it suffices to consider the cardinality sizes $\left|\mathcal{U}\right|,\left|\mathcal{V}\right|\le \left|\mathcal{X}\right|$.

The original objective function can be reformulated as

$$\begin{array}{cc}\hfill G(\alpha ,\beta ,q)& =\overline{\theta }\alpha (I(X;Y|V)+I(V;Z))-(\overline{\theta }\alpha -\theta \overline{\beta })I(X;Y|V)\hfill \\ \hfill & +\theta \beta (I(X;Z|U)+I(U;Y))-(\theta \beta -\overline{\theta }\overline{\alpha })I(X;Z|U).\hfill \end{array}$$

The right-hand side contains two parts, both of which are similar to Equation (10), i.e., the objective function of the degraded broadcast channel. It seems workable to apply the BA algorithm twice, however, it should be noted that these two parts are coupled by the same $q\left(x\right)$.

Observe that the first part depends only on $q(v,x)$, while the other depends on $q(u,x)$. It suffices to consider the subset of distributions such that $q(u,v,x)=q\left(x\right)q\left(v\right|x\left)q\right(u\left|x\right)$. Thus, it is natural to decouple these two parts by fixing $q\left(x\right)$ and applying the BA algorithm separately to $q\left(v\right|x)$ and $q\left(u\right|x)$. After some manipulations, we have

$$\begin{array}{cc}\hfill & G(\alpha ,\beta ,q,Q)\hfill \\ \hfill & =(\overline{\theta }\alpha +\theta \beta )H\left(X\right)\hfill \\ \hfill & +\overline{\theta }\alpha (H\left(V\right|X)-H\left(X\right|Y,V)-H\left(V\right|Z))-(\overline{\theta }\alpha -\theta \overline{\beta })(H\left(Y\right|V)-H\left(Y\right|X))\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & +\theta \beta (H\left(U\right|X)-H\left(X\right|Z,U)-H\left(U\right|Y))-(\theta \beta -\overline{\theta }\overline{\alpha })(H\left(Z\right|U)-H\left(Z\right|X))\hfill \\ \hfill & =-(\overline{\theta }\alpha +\theta \beta )\sum _{x}q\left(x\right)lnq\left(x\right)\hfill \\ \hfill & +\overline{\theta }\alpha \sum _{x,v}q\left(x\right)q\left(v\right|x)(d_1\left[Q\right](v,x)-lnq\left(v\right|x))\hfill \\ \hfill & +\theta \beta \sum _{x,u}q\left(x\right)q\left(u\right|x)(d_2\left[Q\right](u,x)-lnq\left(u\right|x)).\hfill \end{array}$$

(30)

The functions $d_1$ and $d_2$ in the above are

$$\begin{array}{cc}\hfill d_1\left[Q\right](v,x)& =\sum _{y,z}p(y,z|x)\left(lnQ\left(x\right|y,v)Q\left(v\right|z)+\frac{\overline{\theta }\alpha -\theta \overline{\beta }}{\overline{\theta }\alpha }ln\frac{Q\left(y\right|v)}{p\left(y\right|x)}\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill d_2\left[Q\right](u,x)& =\sum _{y,z}p(y,z|x)\left(lnQ\left(x\right|z,u)Q\left(u\right|y)+\frac{\theta \beta -\overline{\theta }\overline{\alpha }}{\theta \beta }ln\frac{Q\left(z\right|u)}{p\left(z\right|x)}\right).\hfill \end{array}$$

(32)

For fixed $q\left(x\right)q\left(v\right|x\left)q\right(u\left|x\right)$, according to Equation (5), the optimal Qs are induced $Q\left[q\right]$s. For fixed Qs, according to Equation  (7), for each x we have

$$\begin{array}{cc}\hfill q\left[Q\right]\left(v\right|x)& =\frac{exp\left\{d_1\left[Q\right](v,x)\right\}}{{\sum }_{v^{\prime }}exp\left\{d_1\left[Q\right](v^{\prime },x)\right\}},\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill q\left[Q\right]\left(u\right|x)& =\frac{exp\left\{d_2\left[Q\right](u,x)\right\}}{{\sum }_{u^{\prime }}exp\left\{d_2\left[Q\right](u^{\prime },x)\right\}}.\hfill \end{array}$$

(34)

The value of the objective function is

$$\begin{array}{cc}\hfill & G(\alpha ,\beta ,q\left(x\right),q\left[Q\right]\left(v\right|x),q\left[Q\right]\left(u\right|x),Q)=(\overline{\theta }\alpha +\theta \beta )\sum _{x}q\left(x\right)(d\left[Q\right]\left(x\right)-lnq\left(x\right)),\hfill \end{array}$$

(35)

where

$$\begin{array}{c}\hfill d\left[Q\right]\left(x\right)=\frac{\overline{\theta }\alpha }{\overline{\theta }\alpha +\theta \beta }ln\sum _{v}exp\left\{d_1\left[Q\right](v,x)\right\}+\frac{\theta \beta }{\overline{\theta }\alpha +\theta \beta }ln\sum _{u}exp\left\{d_2\left[Q\right](u,x)\right\}.\end{array}$$

(36)

Again, according to Equation (7) the optimal $q\left[Q\right]\left(x\right)$ and corresponding function value are

$$\begin{array}{cc}\hfill q\left[Q\right]\left(x\right)& =\frac{exp\left\{d\right[Q\left]\right(x\left)\right\}}{{\sum }_{x^{\prime }}exp\left\{d\left[Q\right]\left(x^{\prime }\right)\right\}},\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill G(\alpha ,\beta ,q[Q],Q)& =(\overline{\theta }\alpha +\theta \beta )ln\sum _{x}exp\left\{d\left[Q\right]\left(x\right)\right\}.\hfill \end{array}$$

(38)

For minimization over $(\alpha ,\beta )$, similar to Section 3.1, we update $(\alpha ,\beta )$ along the gradient:

$$\begin{array}{c}\hfill {\alpha }_{k+1}={\alpha }_k-{\tau }_k·\frac{\partial }{\partial \alpha }G({\alpha }_k,{\beta }_k,q)={\alpha }_k-{\tau }_k·\overline{\theta }·(I(V;Z)-I(X;Z|U)),\end{array}$$

(39)

$$\begin{array}{c}\hfill {\beta }_{k+1}={\beta }_k-{\tau }_k·\frac{\partial }{\partial \beta }G({\alpha }_k,{\beta }_k,q)={\beta }_k-{\tau }_k·\theta ·(I(U;Y)-I(X;Y|V)).\end{array}$$

(40)

Here, it should be mentioned that $(\alpha ,\beta )$ must satisfy the constraint $\overline{\theta }\alpha +\theta \beta \ge max\{\theta ,\overline{\theta }\}$. Thus, if the resulting $({\alpha }_{k+1},{\beta }_{k+1})$ violate this constraint, then we need to scale $\overline{\theta }{\alpha }_{k+1}+\theta {\beta }_{k+1}$ up to be (at least) equal to $max\{\theta ,\overline{\theta }\}$. One way to accomplish this is to use the equality to make $\beta$ dependent on $\alpha$, in which case the gradient descent update becomes ${\alpha }_{k+1}={\alpha }_k-{\tau }_k·d_k$, ${\beta }_{k+1}={\beta }_k-{\tau }_k·(-\overline{\theta }{d}_k/\theta )$, where

$$\begin{array}{c}\hfill d_k=\frac{\partial }{\partial \alpha }G({\alpha }_k,{\beta }_k,q)-\frac{\overline{\theta }}{\theta }·\frac{\partial }{\partial \beta }G({\alpha }_k,{\beta }_k,q).\end{array}$$

Similar to Algorithm 2, we summarize the algorithm for UVOB in Algorithm 4.

Algorithm 4: Computing the UV outer bound
Input: $p(y,z|x)$, maximum iterations K, N, thresholds $\eta$, $ϵ>0$, step size $\tau >0$;
Initialization: ${\alpha }_0,{\beta }_0\in (0,1)$, $q_0(u,v,x)>0$, $\overline{\theta }{\alpha }_0+\theta {\beta }_0\ge max\{\theta ,\overline{\theta }\}$, ${\eta }_{\alpha },{\eta }_{\beta }>\eta$,
$k=0$;
while $k<K$ and $max\{{\eta }_{\alpha },{\eta }_{\beta }\}>\eta$ do
		initialize ${ϵ}_q>ϵ$, $n=0$;
		while $n<N$ and ${ϵ}_q>ϵ$ do
				$n\leftarrow n+1$;
				$Q_n=Q\left[q_{n-1}\right]$ using Equation (5) similarly;
				$q_n=q\left[Q_n\right]$ using Equations (33), (34) and (37);
				$G({\alpha }_k,{\beta }_k,q_n,Q_n)=(\overline{\theta }\alpha +\theta \beta )·ln{\sum }_{x}exp\left\{d\left[Q_n\right]\right\}$ using Equation (36);
				${ϵ}_q=(\overline{\theta }\alpha +\theta \beta )·max\{d\left[Q_n\right]\left(x\right)-ln{q}_{n-1}\left(x\right)\}-G({\alpha }_k,{\beta }_k,q_n,Q_n)$;
		end
		$k\leftarrow k+1$;
		calculate ${\alpha }_k$ and ${\beta }_k$ using Equations (39) and (40);
		${\alpha }_k\leftarrow min\{1,max\{0,{\alpha }_k\}\}$;
		${\beta }_k\leftarrow min\{1,max\{0,{\beta }_k\}\}$;
		if $\overline{\theta }{\alpha }_k+\theta {\beta }_k<max\{\theta ,\overline{\theta }\}$, scale up to equality;
		${\eta }_{\alpha }=|{\alpha }_k-{\alpha }_{k-1}|$
		${\eta }_{\beta }=|{\beta }_k-{\beta }_{k-1}|$;
		$q_0\leftarrow q_n$;
end
Output: ${\alpha }_k$, ${\beta }_k$, $q_n(u,v,x)$, $Q_n$, $G({\alpha }_{k-1},{\beta }_{k-1},q_n,Q_n)$

4. Convergence Analysis

Here, we aim to show that certain convergence results hold if $q_n$ lies in a proper convex set which contains the global maximizer $q_{\ast }$. For this purpose, we first introduce the first-order characterization of a concave function.

Lemma 4

(Lemma 3 in [22]). Given a convex set $\mathcal{S}$, a differentiable function f is concave in $\mathcal{S}$ if and only if, for all $x,y\in \mathcal{S}$,

$$\begin{array}{c}\hfill f\left(y\right)-f\left(x\right)-{(y-x)}^T{\nabla f|}_x\le 0.\end{array}$$

(41)

Similar to [22], we use the superlevel set to construct the convex set $\mathcal{S}$. Let ${\mathcal{S}}_F(\alpha ,k)$ be the superlevel set of the objective function $F(\alpha ,q)$ of SCIB:

$$\begin{array}{c}\hfill {\mathcal{S}}_F(\alpha ,k):=\left\{q\right|F(\alpha ,q)\ge k\}.\end{array}$$

(42)

For a fixed k, it is possible for ${\mathcal{S}}_F(\alpha ,k)$ to contain more than one connected set. For $q\in {\mathcal{S}}_F(\alpha ,k)$, we denote the connected set that contains q as ${\mathcal{T}}_F(\alpha ,k,q)$.

Similarly, for MIB and UVOB we define the corresponding (connected) superlevel sets: ${\mathcal{S}}_M(\alpha ,k)$, ${\mathcal{T}}_M(\alpha ,k,q)$, ${\mathcal{S}}_G(\alpha ,\beta ,k)$, ${\mathcal{T}}_G(\alpha ,\beta ,k,q)$. Note that k here should not be confused with the notation indicating the number of iterations in the algorithms.

4.1. Superposition Coding Inner Bound

According to Theorem 3, the expression of the objective function $F(\alpha ,q)$ of SCIB depends on the value of $\alpha$. Without loss of generality, we can consider the objective function depicted in Equation (16). An equivalent condition for $F(\alpha ,q)$ to be concave is provided in the following lemma.

Lemma 5.

Given a convex set $\mathcal{S}$ with a distribution $q(u,x)$, then $F(\alpha ,q)$ as depicted in Equation (16) is concave in $\mathcal{S}$ if and only if, for all $q_1$, $q_2\in \mathcal{S}$, we have

$$\begin{array}{c}\hfill \overline{\theta }(-D_{UX}+D_{X|YU}+\overline{\alpha }{D}_{U|Y}+\alpha D_{U|Z})+(\overline{\theta }-\theta )D_{Y|U}\le 0,\end{array}$$

where $D_{A|B}$ denotes $D_{A|B}(q_2\left|\right|q_1)$.

The following lemma shows that $q_{n+1}$ lies in the same connected superlevel set as that of $q_n$. The proof (see Appendix C) is similar to that for Lemma 4 in [16].

Lemma 6.

In Algorithm 2, if $q_n(u,x)\in {\mathcal{S}}_F(\alpha ,k)$, then $q_{n+1}\in {\mathcal{T}}_F(\alpha ,k,q_n)$.

Fixing $\alpha$ and letting $q_{\ast }(u,x)$ be the maximizer, the following theorem states that the function values $F(\alpha ,q_n,Q_n)$ converge. The proof (see Appendix C) is similar to that of Theorem 2.

Theorem 5.

If $q_{\ast },q_0\in {\mathcal{T}}_F(\alpha ,k,\tilde{q})$ for some k and $\tilde{q}$, and if $F(\alpha ,q)$ is concave in ${\mathcal{T}}_F(\alpha ,k,\tilde{q})$, then the sequence $F(\alpha ,q_n,Q_n)$ generated by Algorithm 2 converges monotonically from below to $F(\alpha ,q_{\ast })$.

The following corollary is implied by the proof of Theorem 5.

Corollary 1.

If $q_{\ast },q_0\in {\mathcal{T}}_F(\alpha ,k,\tilde{q})$ for some k and $\tilde{q}$, and if $F(\alpha ,q)$ is concave in ${\mathcal{T}}_F(\alpha ,k,\tilde{q})$, then

$$\begin{array}{cc}\hfill & F(\alpha ,q_{\ast })-F(\alpha ,q_N,Q_N)\le \frac{\overline{\theta }}{N}·D_{UX}(q_{\ast }\left|\right|q_0).\hfill \end{array}$$

The above analyses deal with $F(\alpha ,q_n,Q_n)$ for a fixed $\alpha$. When ${\alpha }_m$ changes to ${\alpha }_{m+1}$, the estimation for the one-step change in the function value is presented in the following proposition (see proof in Appendix C).

Proposition 2.

Given ${\alpha }_m$, suppose that Algorithm 2 converges to the optimal variables ${\tilde{q}}_{\ast }$ and ${\tilde{Q}}_{\ast }$ such that ${\tilde{q}}_{\ast }=q\left[{\tilde{Q}}_{\ast }\right]$ and ${\tilde{Q}}_{\ast }=Q\left[{\tilde{q}}_{\ast }\right]$. Letting ${\alpha }_{m+1}$ be updated using Equation (22) and letting $q_0={\tilde{q}}_{\ast }$ be the initial point for the next round, we have

$$\begin{array}{c}\hfill F({\alpha }_{m+1},q_1,Q_1)\approx F({\alpha }_m,{\tilde{q}}_{\ast })-\frac{{({\alpha }_{m+1}-{\alpha }_m)}^2}{{\tau }_m}.\end{array}$$

4.2. Marton’s Inner Bound

Next, we present the convergence results of the BA algorithm for MIB. The proofs are omitted, as they are similar to those for SCIB.

Lemma 7.

Given a convex set $\mathcal{S}$ with a distribution $q(u,v,w,x)$, $M(\alpha ,q)$ as depicted in Equation (24) is concave in $\mathcal{S}$ if and only if, for all $q_1$, $q_2\in \mathcal{S}$, we have

$$\begin{array}{cc}\hfill & -\overline{\theta }{D}_{UVWX}+(\overline{\theta }-\alpha \theta )D_{W|Z}+\alpha \theta D_{W|Y}+\overline{\theta }{D}_{V|WZ}\hfill \\ \hfill & +\theta D_{U|WY}+(\overline{\theta }-\theta )D_{U|VW}+\overline{\theta }{D}_{X|UVW}\le 0,\hfill \end{array}$$

where $D_{A|B}$ denotes $D_{A|B}(q_2\left|\right|q_1)$.

Lemma 8.

In Algorithm 3, if $q_n(u,v,w,x)\in {\mathcal{S}}_M(\alpha ,k)$, then $q_{n+1}\in {\mathcal{T}}_M(\alpha ,k,q_n)$.

Fixing $\alpha$ and letting $q_{\ast }(u,v,w,x)$ be the maximizer, the following theorem states that the function values $M(\alpha ,q_n,Q_n)$ converge.

Theorem 6.

If $q_{\ast },q_0\in {\mathcal{T}}_M(\alpha ,k,\tilde{q})$ for some k and $\tilde{q}$ and if $M(\alpha ,q)$ is concave in ${\mathcal{T}}_M(\alpha ,k,\tilde{q})$, then the sequence $M(\alpha ,q_n,Q_n)$ generated by Algorithm 3 converges monotonically from below to $M(\alpha ,q_{\ast })$.

The following corollary is implied by the proof of Theorem 6.

Corollary 2.

If $q_{\ast },q_0\in {\mathcal{T}}_M(\alpha ,k,\tilde{q})$ for some k and $\tilde{q}$, and if $M(\alpha ,q)$ is concave in ${\mathcal{T}}_M(\alpha ,k,\tilde{q})$, then

$$\begin{array}{cc}\hfill & M(\alpha ,q_{\ast })-M(\alpha ,q_N,Q_N)\le \frac{\overline{\theta }}{N}·D_{UVWX}(q_{\ast }\left|\right|q_0).\hfill \end{array}$$

The estimation for the one-step change in the function value for MIB is presented in the following proposition.

Proposition 3.

Given ${\alpha }_m$, suppose that Algorithm 3 converges to the optimal variables ${\tilde{q}}_{\ast }$ and ${\tilde{Q}}_{\ast }$ such that ${\tilde{q}}_{\ast }=q\left[{\tilde{Q}}_{\ast }\right]$ and ${\tilde{Q}}_{\ast }=Q\left[{\tilde{q}}_{\ast }\right]$. Letting ${\alpha }_{m+1}$ be updated using Equation (27) and letting $q_0={\tilde{q}}_{\ast }$ be the initial point for the next round, we have

$$\begin{array}{c}\hfill M({\alpha }_{m+1},q_1,Q_1)\approx M({\alpha }_m,{\tilde{q}}_{\ast })-\frac{{({\alpha }_{m+1}-{\alpha }_m)}^2}{{\tau }_m}.\end{array}$$

4.3. UV Outer Bound

Now, we present the convergence results of of the BA algorithm for UVOB. The proofs are again omitted, as they are similar to those of SCIB.

Lemma 9.

Given a convex set $\mathcal{S}$ of distribution $q(u,v,x)$, $G(\alpha ,\beta ,q)$ as depicted in Equation (29) is concave in $\mathcal{S}$ if and only if, for all $q_1$, $q_2\in \mathcal{S}$, we have

$$\begin{array}{cc}\hfill -(\overline{\theta }\alpha +\theta \beta )D_X& +\overline{\theta }\alpha (-D_{V|X}+D_{X|YV}+D_{V|Z})+(\overline{\theta }\alpha -\theta \overline{\beta })D_{Y|V}\hfill \\ \hfill & +\theta \beta (-D_{U|X}+D_{X|ZU}+D_{U|Y})+(\theta \beta -\overline{\theta }\overline{\alpha })D_{Z|U}\le 0,\hfill \end{array}$$

where $D_{A|B}$ denotes $D_{A|B}(q_2\left|\right|q_1)$.

Lemma 10.

In Algorithm 4, if $q_n(u,v,x)\in {\mathcal{S}}_G(\alpha ,\beta ,k)$, then $q_{n+1}\in {\mathcal{T}}_G(\alpha ,\beta ,k,q_n)$.

Fixing $(\alpha ,\beta )$ and letting $q_{\ast }(u,v,x)$ be the maximizer, the following theorem states that the function values $G(\alpha ,\beta ,q_n,Q_n)$ converge.

Theorem 7.

If $q_{\ast },q_0\in {\mathcal{T}}_G(\alpha ,\beta ,k,\tilde{q})$ for some k and $\tilde{q}$, and if $G(\alpha ,\beta ,q)$ is concave in ${\mathcal{T}}_G(\alpha ,\beta ,k,\tilde{q})$, then the sequence $G(\alpha ,\beta ,q_n,Q_n)$ generated by Algorithm 4 converges monotonically from below to $G(\alpha ,\beta ,q_{\ast })$.

The following corollary is implied by the proof of Theorem 7.

Corollary 3.

If $q_{\ast },q_0\in {\mathcal{T}}_G(\alpha ,\beta ,k,\tilde{q})$ for some k and $\tilde{q}$, and if $G(\alpha ,\beta ,q)$ is concave in ${\mathcal{T}}_G(\alpha ,\beta ,k,\tilde{q})$, then

$$\begin{array}{cc}\hfill & G(\alpha ,\beta ,q_{\ast })-G(\alpha ,\beta ,q_N,Q_N)\le \frac{\overline{\theta }\alpha +\theta \beta }{N}·D_{UVX}(q_{\ast }\left|\right|q_0).\hfill \end{array}$$

The estimation for the one-step change in the function value for UVOB is presented in the following proposition.

Proposition 4.

Given $({\alpha }_m,{\beta }_m)$, suppose that Algorithm 4 converges to the optimal variables ${\tilde{q}}_{\ast }$ and ${\tilde{Q}}_{\ast }$ such that ${\tilde{q}}_{\ast }=q\left[{\tilde{Q}}_{\ast }\right]$ and ${\tilde{Q}}_{\ast }=Q\left[{\tilde{q}}_{\ast }\right]$. Letting $({\alpha }_{m+1},{\beta }_{m+1})$ be updated using Equations (39) and (40) and letting $q_0={\tilde{q}}_{\ast }$ be the initial point for the next round, we have

$$\begin{array}{c}\hfill G({\alpha }_{m+1},{\beta }_{m+1},q_1,Q_1)\approx G({\alpha }_m,{\beta }_m,{\tilde{q}}_{\ast })-\frac{{({\alpha }_{m+1}-{\alpha }_m)}^2}{{\tau }_m}-\frac{{({\beta }_{m+1}-{\beta }_m)}^2}{{\tau }_m}.\end{array}$$

(43)

5. Numerical Results

We take the binary skew-symmetric broadcast channel $p(y,z|x)$ as the test channel. The conditional probability matrices are

$$\begin{array}{c}\hfill P_{Y|X}=\left[\begin{array}{cc}1& 0\\ 0.5& 0.5\end{array}\right],P_{Z|X}=\left[\begin{array}{cc}0.5& 0.5\\ 0& 1\end{array}\right].\end{array}$$

This is perhaps the simplest broadcast channel for which the capacity region is still unknown.

This broadcast channel plays a very important role in research on capacity bounds. It was first studied in [23] to show that the time-sharing random variable is useful for the Cover–van der Meulen inner bound [24,25]. Later, [26,27,28] demonstrated that the sum rate of the UVOB for this broadcast channel is strictly larger than that of the MIB, showing for the first time that at least one of these two bounds are suboptimal.

Our algorithms are important in at least the following sense: supposing that it is not known whether the MIB matches the UVOB (or the other two bounds for a new scenario) and we want to check this; we can perform an exhaustive search on channel matrices of size $2\times 2$ (or of higher dimensions) to check whether they match. According to the results shown below in Section 5.5, this does not take very much time compared with generic algorithms.

In the following, we apply the algorithms to compute the value of the supporting hyperplane $\theta R_Y+\overline{\theta }{R}_Z$, where $\theta =0.4$. The initial values of $\alpha$ and $\beta$ are ${\alpha }_0={\beta }_0=0.7$. This set of parameters is feasible for UVOB, as $\overline{\theta }\alpha +\theta \beta =0.7>max\{\theta ,\overline{\theta }\}$.

We demonstrate the algorithms in the following aspects: (1) the maximization part; (2) the minimization part; (3) the change from the maximum part to the minimization part; (4) the superlevel set; and (5) comparison with generic non-convex algorithms.

5.1. Maximization Part

In this part, we fix $\alpha$ and $\beta$ to the initial values and let the BA algorithms iterate for $N=200$ times. The results are presented in Figure 1. Because this is the maximization part, the function values increase as the iterations proceed. It is clear that the function values behave properly for fixed $\alpha$ and $\beta$.

5.2. Minimization Part

In this part, we start with the initial ${\alpha }_0$ and ${\beta }_0$, then let the algorithms iterate for $K=200$ times. The results for $({\alpha }_k,{\beta }_k)$ are presented in Figure 2. Because this is the minimization part, the function values decrease as the iterations proceed. It is clear that ${\alpha }_k$ in SCIB and MIB gradually changes as k grows. For UVOB, it is necessary to ensure that $\overline{\theta }\alpha +\theta \beta \ge max\{\theta ,\overline{\theta }\}$. When the updated $({\alpha }_{k+1},{\beta }_{k+1})$ makes $\overline{\theta }\alpha +\theta \beta$ fall below this value, it becomes necessary to scale it back. This happens approximately starting from $k=5$.

5.3. Change from Maximization to Minimization

In this part, we consider UVOB and let K in the algorithm be $K=100$. Figure 3 plots the following three values in Equation (43):

$$\begin{array}{cc}\hfill & G({\alpha }_k,{\beta }_k,{\tilde{q}}_{\ast }),\hfill \\ \hfill & LHS:=G({\alpha }_{k+1},{\beta }_{k+1},q_1,Q_1),\hfill \\ \hfill & RHS:=G({\alpha }_k,{\beta }_k,{\tilde{q}}_{\ast })-\frac{{({\alpha }_{k+1}-{\alpha }_k)}^2}{{\tau }_k}-\frac{{({\beta }_{k+1}-{\beta }_k)}^2}{{\tau }_k}.\hfill \end{array}$$

As the algorithm iterates, the estimate in Equation (43) becomes more and more accurate, as $exp\left\{x\right\}\approx 1+x$ and $ln(1+x)\approx x$ for small x.

5.4. Superlevel Set

To visualize the convergence of $q_n$ and its relation with the superlevel set, we take SCIB as an example and fix $q\left(u\right)$ such that $q\left(x\right|u)$ has two free variables. We reformulate the objective function of SCIB depicted in Equation (16) as follows:

$$\begin{array}{cc}\hfill \tilde{F}(\alpha ,q\left(x\right|u),Q)& =\overline{\theta }H\left(U\right)+\overline{\theta }(H\left(X\right|U)-H\left(X\right|Y,U)-\overline{\alpha }H\left(U\right|Y)-\alpha H\left(U\right|Z))\hfill \\ \hfill & -(\overline{\theta }-\theta )(H\left(Y\right|U)-H\left(Y\right|X))\hfill \\ \hfill & =\overline{\theta }H\left(U\right)+\overline{\theta }\sum _{x,u}q\left(u\right)q\left(x\right|u)(d\left[Q\right](u,x)-lnq\left(x\right|u)),\hfill \end{array}$$

where

$$\begin{array}{c}\hfill d\left[Q\right](u,x)=\sum _{y,z}p(y,z|x)\left(lnQ\left(x\right|y,u)Q{\left(u\right|y)}^{\overline{\alpha }}Q{\left(u\right|z)}^{\alpha }+\frac{\overline{\theta }-\theta }{\overline{\theta }}ln\frac{Q\left(y\right|u)}{p\left(y\right|x)}\right).\end{array}$$

In particular, we fix ${\alpha }_0=0.7$ and $P_U=(0.3,0.7)$, then use the algorithm to find the values of $P(X=0|U=0)$ and $P(X=0|U=1)$. The results are shown in Figure 4. In this case, $q_n$ for large enough n lies in the concave part of the superlevel set, meaning that the algorithm converges. Here, it should be mentioned that it is possible that the algorithm may not converge to the optimal point for some initial $q_0$s that do not lie in the concave part.

5.5. Comparison with Generic Non-Convex Algorithms

Here, we compare our algorithms with the following generic algorithms implemented using the “fmincon” MATLAB function: interior-point, active-set, and sequential quadratic programming (sqp). For simplicity, we only compare the sum rate of MIB, for which the optimal value is 0.2506717… nats (0.3616428… bits). The optimization problem for computing the sum rate is

$$\begin{array}{c}\hfill \underset{q(u,v,w,x)}{max}min\{I(W;Y),I(W;Z)\}+I(U;Y|W)+I(V;Z|W)-I(U;V|W).\end{array}$$

According to Lemma 3, the cardinality size is $\left|\mathcal{U}\right|·\left|\mathcal{V}\right|·\left|\mathcal{W}\right|·\left|\mathcal{X}\right|={\left|\mathcal{X}\right|}^3\left(\right|\mathcal{X}|+4)=48$.

Notice that we do not carry out a comparison with the method in [16], as it cannot be applied to cases where there is a minimum. For scenarios in which [16] can be used, our algorithms degenerate to the method in [16].

The initial point of $q(u,v,w,x)$ is randomly generated for all the algorithms.

Table 2. Comparison with generic non-convex algorithms on BSSC.

Method	Time (Seconds)	Sum-Rate of MIB (Nats)
interior-point	513.82	0.25060…
active-set	2438.57	0.25061…
sqp	1.7621	0.25067…
this paper	0.0629	0.25067…

lists the experimental results. For the first three algorithms, a randomly picked starting point usually does not provide a good enough result. Thus, we ran the first three algorithms multiple times until the best function value hit 0.2506 in order to test their effectiveness. It is clear from the table that only sqp can be considered comparable to our algorithms.

For further comparison with sqp, we randomly generated broadcast channels with cardinalities of $\left|\mathcal{X}\right|=3$, 4, 5, 6, and $\left|\mathcal{Y}\right|=\left|\mathcal{Z}\right|=\left|\mathcal{X}\right|$. The corresponding dimensions are ${\left|\mathcal{X}\right|}^3\left(\right|\mathcal{X}|+4)=189$, 512, 1125, 2160. Because the optimal sum rate is not yet known, we ran sqp once to record the running time. The results in

Table 3. Comparison with sqp on random channels with alphabet sizes $\left|\mathcal{X}\right|=3,4,5,6$.

Method	Time (Seconds)				Sum-Rate of MIB (Nats)
	$\left|\mathcal{X}\right|$ = 3	$\left|\mathcal{X}\right|$ = 4	|$\mathcal{X}|$ = 5	$\left|\mathcal{X}\right|$ = 6	$\left|\mathcal{X}\right|$ = 3	$\left|\mathcal{X}\right|$ = 4	$\left|\mathcal{X}\right|$ = 5	$\left|\mathcal{X}\right|$ = 6
sqp	2.6342	22.44	168.12	1065.51	0.1840	0.2348	0.1983	0.2351
this paper	0.0771	0.1031	0.1450	0.2086	0.1863	0.2375	0.2019	0.2423

suggest that our algorithms are highly scalable. This meets our expectation, as the updating formulae in Equations (5) and (7) are all explicit and can be computed rapidly.

6. Discussion and Conclusions

6.1. Initial Points of Algorithms

Taking MIB as an example, we next discuss how to choose the initial points. When there is no prior knowledge on the optimization problem, the initial point is usually generated randomly. In this paper, Theorem 6 and Lemma 7 provide some guidance on the choice of the initial point $p_0(u,v,w,x)$. A possibly workable method is to randomly generate an initial point and slightly perturb it to check whether these two points satisfy the inequality in Lemma 7. If the answer is no, then it is possible that the objective function is not concave in the neighbourhood of this point, and we continue to generate new initial points.

For the initial point ${\alpha }_0$, because it lies in $[0,1]$ it is affordable to perform a grid search, especially when $\left|\mathcal{X}\right|$ is small. For example, we can take 0.1 as the equal space and try each ${\alpha }_0\in \{0,0.1,0.2,\dots ,1\}$. This approach can to some extent help us avoid becoming stuck in local extreme points.

6.2. J Version Outer Bound

As mentioned earlier, the best general outer bound is the J version outer bound proposed in [6]. However, the evaluation of this outer bound turns out to be even harder, as there are additional constraints on the free variables and the auxiliary channel with the joint distribution

$$\begin{array}{c}\hfill q\left(x\right)q(u,v,w|x)q(\tilde{u},\tilde{v},\tilde{w}|x)q(\widehat{u},\widehat{v},\widehat{w}|x)p(y,z|x)T\left(j\right|x,y,z).\end{array}$$

These constraints are presented in Equations (18a)–(18c) and (19a)–(19c) in [6]. Taking Equations (18a) and (19a) as an example,

$$\begin{array}{cc}\hfill \left(18a\right):& I(\tilde{W};Z)-I(\tilde{W};J)+I(\widehat{W};J)-I(\widehat{W};Y)=I(W;Z)-I(W;Y),\hfill \\ \hfill \left(19a\right):& 0\le I(X;Z|\tilde{U},\tilde{W})-I(X;J|\tilde{U},\tilde{W})\le I(\tilde{V};Z|\tilde{W})-I(\tilde{V};J|\tilde{W}).\hfill \end{array}$$

Direct application of Equation (7) does not yield an updated $q\left[Q\right]$ guaranteed to satisfy these constraints; thus, the design of BA algorithms for the J version outer bound should carefully address this kind of problem. We leave this for future research.

Finally, to conclude our paper, the extension of the BA algorithm to inner and outer bounds for general broadcast channels encounters max–min problems. We have shown that the max–min order can be changed to min–max. Based on this observation, we have designed BA algorithms for the maximization parts and gradient descent algorithms for the minimization parts, then performed convergence analysis and numerical experiments to support our analysis. We have compared our algorithms to the following generic non-convex algorithms: interior-point, active-set, and sequential quadratic programming. The results show that our algorithms are both effective and efficient.