Minimum--Entropy Couplings and their Applications

TL;DR

This paper addresses the problem of finding joint distributions with minimum entropy for coupled discrete variables, providing efficient approximation algorithms with bounded additive gaps, and extending the approach to multiple variables and other entropy measures.

Contribution

The paper introduces the first efficient approximation algorithm for minimum-entropy couplings with bounded additive error, extending to multiple variables and alternative entropy measures.

Findings

Provides an algorithm with at most 1-bit excess entropy for two variables.

Extends to $k$ variables with at most $ ext{log }k$ bits excess.

Discusses applications and entropy measure extensions.

Abstract

Given two discrete random variables $X$ and $Y,$ with probability distributions ${\bf p}=(p_1, \ldots , p_n)$ and ${\bf q}=(q_1, \ldots , q_m)$, respectively, denote by ${\cal C}({\bf p}, {\bf q})$ the set of all couplings of ${\bf p}$ and ${\bf q}$, that is, the set of all bivariate probability distributions that have ${\bf p}$ and ${\bf q}$ as marginals. In this paper, we study the problem of finding a joint probability distribution in ${\cal C}({\bf p}, {\bf q})$ of \emph{minimum entropy} (equivalently, a coupling that \emph{maximizes} the mutual information between $X$ and $Y$), and we discuss several situations where the need for this kind of optimization naturally arises. Since the optimization problem is known to be NP-hard, we give an efficient algorithm to find a joint probability distribution in ${\cal C}({\bf p}, {\bf q})$ with entropy exceeding the minimum possible at most by {1 bit}, thus providing an approximation algorithm with an additive gap of at most 1 bit. Leveraging on this algorithm, we extend our result to the problem of finding a minimum--entropy joint distribution of arbitrary $k\geq 2$ discrete random variables $X_1, \ldots , X_k$, consistent with the known $k$ marginal distributions of the individual random variables $X_1, \ldots , X_k$. In this case, our algorithm has an { additive gap of at most $\log k$ from optimum.} We also discuss several related applications of our findings and {extensions of our results to entropies different from the Shannon entropy.}

Equations

$$\mathtt{D}(\mathbf{p},\mathbf{q})=2W(\mathbf{p},\mathbf{q})-H(\mathbf{p})-H(\mathbf{q}),$$

$${\mathcal{P}}_n=\{(p_1,\dots ,p_n):\sum _{i=1}^n{p}_i=1,p_1\ge \dots \ge p_n\ge 0\}$$

$$\mathbf{x}⪯\mathbf{x}\vee \mathbf{y},\mathbf{y}⪯\mathbf{x}\vee \mathbf{y},$$

$$\mathbf{x}\vee \mathbf{y}⪯\mathbf{z}.$$

$$\mathbf{x}\wedge \mathbf{y}⪯\mathbf{x},\mathbf{x}\wedge \mathbf{y}⪯\mathbf{y},$$

$$\mathbf{z}⪯\mathbf{x}\wedge \mathbf{y}.$$

$$z_{i}=\min\Bigl{\{}\sum_{j=1}^{i}p_{j},\sum_{j=1}^{i}q_{j}\Bigr{\}}-\sum_{j=1}^{i-1}z_{j}.$$

$$\sum_{k=1}^{i}z_{k}=\min\Bigl{\{}\sum_{k=1}^{i}p_{k},\sum_{k=1}^{i}q_{k}\Bigr{\}}.$$

$$\sum_{k=i}^{n}z_{k}=\max\Bigl{\{}\sum_{k=i}^{n}p_{k},\sum_{k=i}^{n}q_{k}\Bigr{\}}.$$

$$H(\mathbf{x})\ge H(\mathbf{y})+D(\mathbf{y}||\mathbf{x}),$$

$${\cal C}({\bf p},{\bf q})=\Bigl{\{}\mathbf{M}=[m_{ij}]:\sum_{j}m_{ij}=p_{i},\sum_{i}m_{ij}=q_{j}\Bigr{\}}.$$

$$\mathbf{m}⪯\mathbf{p}\text{mbox}and\mathbf{m}⪯\mathbf{q}.$$

$$\mathbf{m}⪯\mathbf{p}\wedge \mathbf{q}.$$

$$H(\mathbf{m})\ge H(\mathbf{p}\wedge \mathbf{q}).$$

$$H(\mathbf{M})\ge H(\mathbf{p}\wedge \mathbf{q}).$$

$$\underset{\mathbf{N}\in \mathcal{C}(\mathbf{p},\mathbf{q})}{\min }H(\mathbf{N})\ge H(\mathbf{p}\wedge \mathbf{q}),$$

$$H(X,Y)\ge \max \{H(X),H(Y)\}.$$

$$H(\mathbf{M})\ge \max \{H(\mathbf{p}),H(\mathbf{q})\}.$$

$$H(XY)=H(\mathbf{M})\ge H(\mathbf{p}\wedge \mathbf{q})\ge \max \{H(\mathbf{p}),H(\mathbf{q})\}.$$

$$I(X;Y)\le H(\mathbf{p})+H(\mathbf{q})-H(\mathbf{p}\wedge \mathbf{q})\le \min \{H(X),H(Y)\}.$$

$$H(X|Y)$$

$$\ge H(X)-H(\mathbf{p})-H(\mathbf{q})+H(\mathbf{p}\wedge \mathbf{q})$$

$$=H(\mathbf{p}\wedge \mathbf{q})-H(\mathbf{q}).$$

$$H(X)+H(Y)$$

$$+D(X+Y||(X,Y))\ge H(X+Y),$$

$$H(Z)\ge H(f(Z))+D(f(Z)||Z),$$

$$H(\mathbf{M})\le H(\mathbf{p}\wedge \mathbf{q})+1.$$

$$H(\mathbf{M})\le \underset{\mathbf{N}\in \mathcal{C}(\mathbf{p},\mathbf{q})}{\min }H(\mathbf{N})+1.$$

$$z^{(d)}+\sum _{i\in I}A[i]=x.$$

$$\text{scalebox}0.8{\mathbf{M}}_{\mathbf{z}}=\left(\begin{array}{cccccc}0.4& & & & & 0\\ & 0.22& & & & 0\\ & & 0.18& & & 0\\ & & & 0.13& & 0\\ & & & & 0.04& 0\\ 0& 0& 0& 0& 0& 0.03\end{array}\right)$$

$$\text{scalebox}0.8{\mathbf{M}}^{(6)}=\left(\begin{array}{cccccc}0.4& & & & & 0\\ & 0.22& & & & 0\\ & & 0.18& & & 0\\ & & & 0.13& & 0\\ & & & & 0.04& 0\\ 0& 0& 0& 0& \underline{0.01}& 0.02\end{array}\right)$$

Full text

I Introduction and Motivations

Inferring an unknown joint distribution of two random variables (r.v.), when only their marginals are given, is an old problem in the area of probabilistic inference. The problem goes back at least to Hoeffding [27] and Frechet [18], who studied the question of identifying the extremal joint distribution of r.v. $X$ and $Y$ that maximizes (resp., minimizes) their correlation, given the marginal distributions of $X$ and $Y$. We refer the reader to [3, 11, 13, 37] for a (partial) account of the vast literature in the area and the many applications in the pure and applied sciences.

In this paper, we consider the following case of the problem described above. Let $X$ and $Y$ be two discrete r.v., distributed according to ${\bf p}=(p_{1},\ldots,p_{n})$ and ${\bf q}=(q_{1},\ldots,q_{m})$, respectively. We seek a minimum-entropy joint probability distribution of $X$ and $Y$, whose marginals are equal to ${\bf p}$ and ${\bf q}$.

We discuss below a few scenarios where this problem naturally arises.

I-A Entropic Causal Inference

In papers [31] and [32], the authors consider the important question of identifying the correct causal direction between two arbitrary r.v.’s $X$ and $Y$, that is, they want to discover whether it is the case that $X$ causes $Y$ or it is $Y$ that causes $X$. In general, $X$ causes $Y$ if there exists an exogenous r.v. $E$ (independent of $X$) and a deterministic function $f$ such that $Y=f(X{,}E)$. In order to identify the correct causal direction (i.e., either from $X$ to $Y$ or from $Y$ to $X$), the authors of [31] and [32] make the reasonable postulate that the entropy of the exogenous r.v. $E$ is small in the true causal direction, and empirically validate this assumption. Additionally, they prove the important fact that the problem of finding the exogenous variable $E$ with minimum entropy is equivalent to the problem of finding the minimum-entropy joint distribution of properly defined random variables, given (i.e., by fixing) their marginal distributions (see Theorem 3 of [32]). This is exactly the problem we consider in this paper. The authors of [32] observe that the latter optimization problem is NP-hard (due to results of [34] and [54]) and propose a greedy approximation algorithm to find the minimum-entropy joint distribution, given the marginals. For this greedy algorithm, the authors prove that it always finds a local minimum and that the local minimum is within an additive guaranteed gap from the unknown global optimum. The authors of [32] observe that this additive guaranteed gap can be as large as $\log n$ (here $n$ is the cardinality of the support of each involved random variable). Similar results are contained in [41], and references therein.

In this paper, we design a different greedy algorithm, and we prove that it returns a correct joint probability distribution (i.e., with the prescribed marginals) with entropy exceeding the minimum possible by at most 1 bit. Subsequently, in Section V, we extend our algorithm to the case of more than two random variables. More precisely, we consider the problem of finding a minimum–entropy joint distribution of arbitrary $k\geq 2$ discrete random variables $X_{1},\ldots,X_{k}$, consistent with the known $k$ marginal distributions of $X_{1},\ldots,X_{k}$. In this case, our algorithm has an additive guaranteed gap of at most $\log k$.

I-B On the functional representation of correlated random variables

Let $X$ and $Y$ be two arbitrary random variables with joint distribution $p(x,y)$. The functional representation lemma [16, p. 626] states that there exists a random variable $Z$ independent of $X$, and a deterministic function $f$, such that the pair of r.v.’s $(X,f(X,Z))$ is distributed like $(X,Y)$, that is, they have the same joint distribution $p(x,y)$. This lemma has been applied to establish several results in network information theory (see [16] and references therein). In several applications, it is important to find a r.v. $Z$ such that the conditional entropy $H(Y|Z)$ is close to its natural lower bound, that is, it is close to $I(X;Y)$. Recently, a very strong result to that respect was proved in [36], showing that one can indeed find a r.v. $Z$ such that $H(Y|Z)\leq I(X;Y)+\log(I(X;Y)+1)+4$ bits. Among the papers that have used (versions of) the functional representation lemma, papers [4] and [26] have considered the problem of one-shot channel simulation with unlimited common randomness.111The situation is also somewhat reminiscent of the important “reverse Shannon Theorem” of [2], where one wants to simulate an arbitrary noisy channel with a noiseless one, plus some additional source of randomness (see [2] for formal definitions). The bounds on $H(Y|Z)$ are essentially used to set a limit on the amount of bits exchanged among parties involved in the simulation. However, there is here another resource which is reasonable to bound: the amount of randomness used in the protocol, i.e.,, the amount of randomness needed to generate the auxiliary r.v. $Z$. Indeed, randomness is not free, and several clever (but expensive) methods have been devised to produce it, based on physical systems like Geiger-Muller tubes, chaotic laser, etc.; therefore it seems reasonable to require that the entropy of auxiliary r.v. $Z$ be minimum222This requirement can be made more formal by invoking the fundamental result of Knuth and Yao [30], stating that the minimum average number of unbiased random bits necessary to generate an arbitrary discrete r.v. $Z$ is sandwiched between $H(Z)$ and $H(Z)+2$.. On the basis of the (already mentioned) important result proved in [31], showing the equivalence between the problem of finding the minimum entropy auxiliary variable $Z$ such that $(X,f(X,Z))=(X,Y)$, and the problem of finding the minimum-entropy joint distribution of properly defined random variables (given their marginal distributions), it follows that the results of our paper offer a solution to the question of seeking a minimum entropy r.v. $Z$ such that the pair of r.v. $(X;Y)$ can be simulated as $(X,f(X,Z))$. The exact statement of our result is given in Corollary 2 of Section V.

I-C Metric for dimension reduction

Another work that considers the problem of finding the minimum-entropy joint distribution of two r.v. $X$ and $Y$, given the marginals of $X$ and $Y$, is the paper [54]. There, the author introduces a pseudo-metric ${\tt D}(\cdot,\cdot)$ among discrete probability distributions in the following way: given arbitrary ${\bf p}=(p_{1},\ldots,p_{n})$ and ${\bf q}=(q_{1},\ldots,q_{m})$, the distance ${\tt D}({\bf p},{\bf q})$ among ${\bf p}$ and ${\bf q}$ is defined as the quantity

[TABLE]

where $W({\bf p},{\bf q})$ is the minimum entropy of a bivariate probability distribution that has ${\bf p}$ and ${\bf q}$ as marginals, and $H$ denotes the Shannon entropy. This metric is applied in [54] to the problem of dimension-reduction of stochastic processes. The author of [54] observes that the problem of computing $W({\bf p},{\bf q})$ is NP-hard (see also [34]) and proposes another different greedy algorithm for its computation, based on some analogy with the problem of Bin Packing with overstuffing. No performance guarantee is given in [54] for the proposed algorithm. Our result directly implies that we can compute the value of the pseudometric ${\tt D}({\bf p},{\bf q})$, for arbitrary ${\bf p}$ and ${\bf q}$, with an additive gap of at most $1$ bit.333We remark that in [8, 9] we considered the different problem of computing the probability distributions ${\bf q}^{*}$ that minimizes ${\tt D}({\bf p},{\bf q})$, given ${\bf p}$.

I-D Contingency tables and transportation polytopes

In the field of Combinatorial Optimization, the set ${\cal C}({\bf p},{\bf q})$ of all couplings of given ${\bf p}=(p_{1},\ldots,p_{n})$ and ${\bf q}=(q_{1},\ldots,q_{m})$ is known as the transportation polytope $P$ defined by ${\bf p}$ and ${\bf q}$. The fact that in our case ${\bf p}$ and ${\bf q}$ are probability distributions (i.e., their components are non-negative and sum up to 1) is immaterial, since one can always normalize. A similar concepts is known in Statistics under the name of contingency tables [15]. Polytopes ${\cal C}({\bf p},{\bf q})$ are called transportation polytopes because they model the transportation of goods from $n$ supply locations (with the $i$-th location supplying a quantity of $p_{i}$) to $m$ demand locations (with the $j$-th location demanding a quantity of $q_{j}$). The feasible points $M_{i,j}$ of an element ${\bf M}=[M_{i,j}]\in{\cal C}({\bf p},{\bf q})$ model the scenario where a quantity of $M_{i,j}$ of goods is transported from the $i$-th supply location to the $j$-th demand location. Many hard combinatorial optimization problems become solvable on the transportation polytope because of its rich and well studied combinatorial structure. We refer to the survey paper [14] for an account of the vast literature on the topic. The problem we consider in this paper can be equivalently stated as the one of finding a minimum-entropy element of the polytope ${\cal C}({\bf p},{\bf q})$. To see that this is not a simple translation of a problem from one language into another, we point out a recent important trend in the area of combinatorial optimization, that is, the one that seeks sparse solutions to optimization problems. More precisely, researchers aim at finding algorithms that trade the optimality of a solution to a given problem with the sparseness of the solution (e.g., the number of variables with non-zero values in the solution, but other measure of sparseness can be employed). We address the reader to [1, 49], and references therein, for motivations and a review of this line of research. Our problem of finding a minimum entropy element in the transportation polytope ${\cal C}({\bf p},{\bf q})$ fits in the above perspective. This interpretation is possible not only because entropy can be often interpreted as a reasonable measure of sparseness (see [28]) but also because our algorithm produces an elements of ${\cal C}({\bf p},{\bf q})$ whose number of non zero elements is, in the worst case, at most twice the minimum possible. We remark that finding a matrix ${\bf M}\in{\cal C}({\bf p},{\bf q})$ with the minimum number of non-zero entries is NP-hard in general [31].

I-E Additional relations and applications

We believe that the problem of finding a minimum entropy joint distribution, with fixed marginal, is indeed a basic one. In this section we will briefly illustrate a few other scenarios where the problem matters.

Let us write the joint entropy of two r.v. $X$ and $Y$, distributed according to ${\bf p}$ and ${\bf q}$, respectively, as $H(X,Y)=H(X)+H(Y)-I(X;Y)$, where $I(X;Y)$ is the mutual information between $X$ and $Y$. Then, one sees that our original problem can be equivalently stated as the determination of a joint probability distribution of $X$ and $Y$ (having given marginals ${\bf p}$ and ${\bf q}$) that maximizes the mutual information $I(X;Y)$. In [34], this maximal mutual information is interpreted, in agreement with Renyi’s axioms for a bona fide dependence measure [44], as a measure of the largest possible correlations between two r.v. $X$ and $Y$. One can see the soundness of this interpretation also in the following way. Let $q(x,y)$ be an arbitrary joint distribution of r.v. $X$ and $Y$, with marginals equal to ${\bf p}=\{p(x)\}$ and ${\bf q}=\{q(y)\}$, respectively. Then, it is well known that the mutual information $I(X;Y)$ can be written as the relative entropy (divergence) between the joint distribution $q(x,y)$ and the joint distribution $r(x,y)=p(x)q(y)$ (i.e., a joint distribution that would make $X$ and $Y$ independent). Therefore, our problem of maximizing $I(X;Y)$ is equivalent to the one of finding a joint distribution of $X$ and $Y$ that is the farthest (in the sense of relative entropy) from $r(x,y)=p(x)q(y)$, that is, finding the joint distribution of $X$ and $Y$ that makes them “most dependent” (or correlated) as possible. Another way to see the question is to realize that we are seeking a joint distribution that minimizes the conditional entropies $H(X|Y)$ and $H(Y|X)$, that represent sensible measures of the strength of the dependence between $X$ and $Y$. Since the problem of its exact computation is NP-hard, our result implies an approximation algorithm for it. We would like to remark that there are indeed situations in which measuring the “potential” correlation between variables (as opposed to their actual correlation) can be useful. For instance, the authors of [35] introduces a measure that, in their words, “provides a score for the strength of the influence protein X has on protein Y. In many physiological conditions, only a small fraction of the cells have activated protein X in response to stimuli, and these active populations have little influence on the mutual information metric”. Since other standard measures of correlation would also fail, using a measure of potential correlation in that context could be useful.

Another situation where the need to maximize the mutual information between two r.v. (with fixed probability distributions) naturally arises, is in the area of medical imaging [43, 55].

In [46], the authors asks several algorithmic problems of this vein: given some fixed marginal distributions, find a joint probability distribution, with marginals close to the given ones, and satisfying some additional properties dictated by certain applications scenarios. Among the several problems considered in [46], the authors mention the problem of finding a minimum entropy joint distribution, with the given marginals, as an interesting open problem.

In the recent paper [57], the authors study the problem of finding a joint probability distribution of two random variables with fixed marginals, that minimizes a given function $f$ of the joint distribution. The authors study this problem in the asymptotic setting (i.e., for product marginal distributions). A strictly related problem was also studied in [52].

In [21], the authors study the problem of finding good upper and lower bounds on the mutual information $I(X;Y)$ of two r.v.’s $X$ and $Y$ when the only available knowledge consists of the marginals of $X$ and $Y$, and the pair of values $(x,y)$ for which the unknown joint distribution of $X$ and $Y$ (consistent with the given marginals) assign a non-zero probability. It is clear that our maximization problem gives an upper bound on $I(X;Y)$ when the available knowledge consists of the marginals of $X$ and $Y$, and nothing else.

Other papers considering problems somewhat related to ours are [19, 29, 38, 40, 47, 56], and [58].

I-F Structure of the paper

The rest of the paper is organized as follows. In Section II we present the mathematical tools and the auxiliary results that are needed to prove our results. In Section III we present our algorithm to find a joint probability distribution of two input r.v.’s $X$ and $Y$, with given marginal distributions, whose entropy is at most one bit away from the joint distribution of minimal entropy. We also present a worked out example to illustrate the behaviour of the algorithm in an intuitive way. The formal proofs of the correctness of the algorithm are spelled out in the successive Section III-B. In Section V we extend the algorithm presented in Section III to an arbitrary number of $k\geq 2$ input random variables. The entropy of the joint distribution produced by our algorithm is at most $\log k$ bits away from the minimum-entropy joint distribution of the $k$ r.v.’s.

Throughout this paper, the logarithms are on base 2 unless specified otherwise.

II Preliminary Results

To prove our results, we use ideas and techniques from majorization theory [39], a mathematical framework that has been proved to be very much useful in information theory (e.g., see [5, 6, 7, 8, 17, 23, 24, 48] and references therein). In this section we recall the notions and results that are relevant to our context.

Definition 1**.**

Given two probability distributions ${\bf a}=(a_{1},\ldots,a_{n})$ and ${\bf b}=(b_{1},\ldots,b_{n})$ with $a_{1}\geq\ldots\geq a_{n}\geq 0$ and $b_{1}\geq\ldots\geq b_{n}\geq 0$, $\sum_{i=1}^{n}a_{i}=\sum_{i=1}^{n}b_{i}=1$, we say that ${\bf a}$ is majorized by ${\bf b}$, and write ${\bf a}\preceq{\bf b}$, if and only if $\sum_{k=1}^{i}a_{k}\leq\sum_{k=1}^{i}b_{k},\quad\mbox{\rm for all }\ i=1,\ldots,n.$

We assume that all the probability distributions we deal with have been ordered in non-increasing order. This assumption does not affect our results, since the quantities we compute (i.e., entropies) are invariant with respect to permutations of the components of the involved probability distributions. We also use the majorization relationship between vectors of unequal lengths, by properly padding the shorter one with the appropriate number of [math]’s at the end. The majorization relation $\preceq$ is a partial ordering on the $(n-1)$-dimensional simplex

[TABLE]

of all ordered probability vectors of $n$ elements, that is, for each ${\bf x},{\bf y},{\bf z}\in{\cal P}_{n}$ it holds that

${\bf x}\preceq{\bf x}$;

• 2)

${\bf x}\preceq{\bf y}$ and ${\bf y}\preceq{\bf z}$ implies ${\bf x}\preceq{\bf z}$;

• 3)

${\bf x}\preceq{\bf y}$ and ${\bf y}\preceq{\bf x}$ implies ${\bf x}={\bf y}$.

It turns out that the partially ordered set $({\cal P}_{n},\preceq)$ is indeed a lattice [5],444The same result was independently rediscovered in [12], see also [22] for a different proof. i.e., for all ${\bf x},{\bf y}\in{\cal P}_{n}$ there exists a unique least upper bound ${\bf x}\lor{\bf y}$ and a unique greatest lower bound ${\bf x}\land{\bf y}$. We recall that the least upper bound ${\bf x}\lor{\bf y}$ is the vector in ${\cal P}_{n}$ such that

[TABLE]

and for all ${\bf z}\in{\cal P}_{n}$ for which ${\bf x}\preceq{\bf z},\ {\bf y}\preceq{\bf z}$ it holds that

[TABLE]

Analogously, the greatest lower bound ${\bf x}\land{\bf y}$ is the vector in ${\cal P}_{n}$ such that

[TABLE]

and for all ${\bf z}\in{\cal P}_{n}$ for which ${\bf z}\preceq{\bf x},\ {\bf z}\preceq{\bf y}$ it holds that

[TABLE]

In the paper [5] the authors also gave a simple and efficient algorithm to explicitly compute ${\bf x}\lor{\bf y}$ and ${\bf x}\land{\bf y}$, given arbitrary vectors ${\bf x},{\bf y}\in{\cal P}_{n}$. Due to the important role it will play in our main result, we recall how to compute the greatest lower bound.

Fact 1**.**

[5]* Let ${\bf x}=(x_{1},\dots,x_{n}),{\bf y}=(y_{1},\dots,y_{n})\in{\cal P}_{n}$ and let ${\bf z}=(z_{1},\dots,z_{n})={\bf x}\land{\bf y}$. Then, $z_{1}=\min\{p_{1},q_{1}\}$ and for each $i=2,\dots,n,$ it holds that*

[TABLE]

Equivalently, we have

[TABLE]

Moreover, using $\sum_{k}z_{k}=\sum_{k}p_{k}=\sum_{k}q_{k}=1,$ we also have that for each $i=1,\dots,n,$ it holds that

[TABLE]

We also remind the important Schur-concavity property of the entropy function [39]:

*For any ${\bf x},{\bf y}\in{\cal P}_{n}$, ${\bf x}\preceq{\bf y}$ implies that $H({\bf x})\geq H({\bf y})$, with equality if and only if ${\bf x}={\bf y}$. *

A notable strengthening of above fact has been proved in [24]. There, the authors prove that ${\bf x}\preceq{\bf y}$ implies

[TABLE]

where $D({\bf y}||{\bf x})$ is the relative entropy between ${\bf x}$ and ${\bf y}$.

We also need the concept of aggregation (see [54] and [8]), and a result from [8], whose proof is repeated here to make the paper self-contained. Given ${\bf p}=(p_{1},\ldots,p_{n})\in{\cal P}_{n}$ and an integer $2\leq m<n$, we say that ${\bf q}=(q_{1},\ldots,q_{m})\in{\cal P}_{m}$ is an aggregation of ${\bf p}$ if there is a partition of $\{1,\ldots,n\}$ into disjoint sets $I_{1},\ldots,I_{m}$ such that $q_{j}=\sum_{i\in I_{j}}p_{i}$, for $j=1,\ldots m$.

Lemma 1**.**

[8]* Let ${\bf q}\in{\cal P}_{m}$ be any aggregation of ${\bf p}\in{\cal P}_{n}$. Then it holds that ${\bf p}\preceq{\bf q}$.*

Proof:

We shall prove by induction on $i$ that $\sum_{k=1}^{i}q_{k}\geq\sum_{k=1}^{i}p_{k}$. Because ${\bf q}$ is an aggregation of ${\bf p}$, we know that there exists $I_{j}\subseteq\{1,\ldots,n\}$ such that $1\in I_{j}$. This implies that $q_{1}\geq q_{j}\geq p_{1}$. Let us suppose that $\sum_{k=1}^{i-1}q_{k}\geq\sum_{k=1}^{i-1}p_{k}$. If there exist indices $j\geq i$ and $\ell\leq i$ such that $\ell\in I_{j}$, then $q_{i}\geq q_{j}\geq p_{\ell}\geq p_{i}$, implying $\sum_{k=1}^{i}q_{k}\geq\sum_{k=1}^{i}p_{k}$. Should it be otherwise, for each $j\geq i$ and $\ell\leq i$ it holds that $\ell\not\in I_{j}$. Therefore, $\{1,\ldots,i\}\subseteq I_{1}\cup\ldots\cup I_{i-1}$. This immediately gives $\sum_{k=1}^{i-1}q_{k}\geq\sum_{k=1}^{i}p_{k}$, from which we get $\sum_{k=1}^{i}q_{k}\geq\sum_{k=1}^{i}p_{k}$. ∎

Let us now discuss some first consequences of the above framework. Given two discrete random variables $X$ and $Y$, with probability distributions ${\bf p}=(p_{1},\ldots,p_{n})$ and ${\bf q}=(q_{1},\ldots,q_{m})$, respectively, denote by ${\cal C}({\bf p},{\bf q})$ the set of all joint distributions of $X$ and $Y$ that have ${\bf p}$ and ${\bf q}$ as marginals. In the literature, elements of ${\cal C}({\bf p},{\bf q})$ are often called couplings of ${\bf p}$ and ${\bf q}$, and play an important role in many information theoretic problems, e.g, see [47]. For our purposes, each element in ${\cal C}({\bf p},{\bf q})$ can be seen as an $n\times m$ matrix $M=[m_{ij}]\in\mathbb{R}^{n\times m}$ such that its row-sums give the elements of ${\bf p}=(p_{1},\ldots,p_{n})$ and its column-sums give the elements of ${\bf q}=(q_{1},\ldots,q_{m})$, that is,

[TABLE]

Now, for any $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$, let us write its elements in a $1\times mn$ vector $\mathbf{m}\in{\cal P}_{mn}$, with its components ordered in non-increasing fashion. From (3) we obtain that both ${\bf p}$ and ${\bf q}$ are aggregations of *each * $\mathbf{m}\in{\cal P}_{mn}$ obtained from some $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$. By Lemma 1, we get that555Recall that we use the majorization relationship between vectors of unequal lengths, by properly padding the shorter one with the appropriate number of [math]’s at the end. This trick does not affect our subsequent results, since we use the customary assumption that $0\log 0=0$.

[TABLE]

Recalling the definition and properties of the greatest lower bound of two vectors in ${\cal P}_{mn}$, we also obtain

[TABLE]

From (5), and the Schur-concavity of the Shannon entropy, we also obtain that

[TABLE]

Since, obviously, the entropy of $H(\mathbf{m})$ is equal to the entropy $H(\mathbf{M})$, where $\mathbf{M}$ is the matrix in ${\cal C}({\bf p},{\bf q})$ from which the vector $\mathbf{m}$ was obtained, we get the following result ( a key one for our purposes).

Lemma 2**.**

For any* ${\bf p}$ and ${\bf q}$, and $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$, it holds that*

[TABLE]

Lemma 2 obviously implies that the minimum-entropy coupling of ${\bf p}$ and ${\bf q}$ that we are seeking satisfies the inequality

[TABLE]

and it is one of the key results towards our algorithm to find an element $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ with entropy at most $1$ bit larger than the entropy of the minimum entropy coupling, i.e., $H(\mathbf{M})\leq OPT+1$, where $OPT=\min_{\mathbf{N}\in{\cal C}({\bf p},{\bf q})}H(\mathbf{N}).$

II-A An interlude

Before describing our algorithm and its analysis, let us illustrate some consequences of Lemma 2 not directly aimed towards proving our main results, but nevertheless of some interest.

It is well known that for any joint distribution of the two r.v. $X$ and $Y$ it holds that

[TABLE]

Since $H(XY)=H(\mathbf{M})$, for some $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ (where ${\bf p}$ and ${\bf q}$ are the marginal distributions of $X$ and $Y$, respectively), we can formulate the bound (7) in the following equivalent way:

for any $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ it holds that

[TABLE]

Lemma 2 allows us to strengthen the lower bound (7). Indeed, by the definition of the greatest lower bound ${\bf p}\land{\bf q}$ of probability distributions ${\bf p}$ and ${\bf q}$, it holds that ${\bf p}\land{\bf q}\preceq{\bf p}$ and ${\bf p}\land{\bf q}\preceq{\bf q}$, and therefore, by the Schur-concavity of the entropy function and Lemma 2 we get the improved lower bound

[TABLE]

Inequality (8) also allows us to improve on the classical upper bound on the mutual information given by $I(X;Y)\leq\min\{H(X),H(Y)\},$ since (8) implies

[TABLE]

The new bounds (8) and (9) are strictly better than the usual ones, whenever ${\bf p}\not\preceq{\bf q}$ and ${\bf q}\not\preceq{\bf p}$. Technically, one could improve them even more, by using the inequality $H({\bf x})\geq H({\bf y})+D({\bf y}||{\bf x})$, whenever ${\bf x}\preceq{\bf y}$ [24]. However, in this paper we just need what we can get from the inequality $H({\bf x})\geq H({\bf y})$, if ${\bf x}\preceq{\bf y}$ holds.

Inequalities (8) and (9) could be useful also in other contexts, when one needs to bound the joint entropy (or the mutual information) of two r.v.’s $X$ and $Y$, and the only available knowledge is given by the marginal distributions of $X$ and $Y$ (and not their joint distribution). Let $X$ and $Y$ be two r.v.’s, where $X$ is distributed according to ${\bf p}$ and $Y$ according to ${\bf q}$, and let $H(X|Y)$ be the conditional entropy of $X$ given $Y$. From (9) we get

[TABLE]

The inequality $H(X|Y)\geq H({\bf p}\land{\bf q})-H({\bf q})$ gives a lower bound on $H(X|Y)$ that does not depend on the joint distribution of $X$ and $Y$. In particular, it also implies that if the probability distributions ${\bf p}$ and ${\bf q}$ of $X$ and $Y$ are such that ${\bf q}\not\preceq{\bf p}$, then the conditional entropy $H(X|Y)$ cannot be zero, no matter what the joint distribution of $X$ and $Y$ is. By the Fano inequality, one gets a lower bound of the error probability $\Pr\{X\neq Y\}$ that depends only on the “structure” of the probability distributions ${\bf p}$ and ${\bf q}$ of $X$ and $Y$ and not on the joint distribution of $X$ and $Y$. Admittedly, this lower bound is weak, but the only fact that one could derive one that is independent from the joint distribution of $X$ and $Y$ seems novel and interesting to us.

Another possible application of the framework of Section II concerns the problem of sumset estimates for Shannon entropy [33, 50]. There, one wants to find upper and lower bounds on the entropy of $H(X+Y),H(X-Y)$ (and similar expressions), in terms of the individual entropies $H(X),H(Y)$. As an example, one could somewhat improve the trivial estimate $H(X)+H(Y)\geq H(X+Y)$ in the following way. Let us consider $X+Y$ and observe that the probability mass function of $X+Y$ is an aggregation of the pmf of the joint random variable $(X,Y)$. Then, by Lemma 1 and formula (2), one immediately gets the inequality

[TABLE]

where the last inequality is strict unless the pmf of $X+Y$ is equal to that of $(X,Y)$. Similar improvements can be obtained for other expressions like $X-Y$. More in general, one has the following inequality that holds for any determinist function $f$ and discrete r.v. $Z$:

[TABLE]

where one recovers (10) when $Z=(X,Y)$ and $f(X,Y)=X+Y$.

III An Algorithm to Approximate  $OPT=\min_{\mathbf{N}\in{\cal C}({\bf p},{\bf q})}H(\mathbf{N}).$

In this section we present our main result, that is, an algorithm that from the input distributions ${\bf p}$ and ${\bf q}$, constructs a coupling $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ such that

[TABLE]

Lemma 2 will imply our desired result, that is

[TABLE]

The following lemma is technical in nature, but it turns out to be a very useful tool of our main algorithm.

Lemma 3**.**

Let $A[1\dots k]$ be an array of $k$ non-negative real numbers and $z$ a positive real number such that $z\geq A[i]$ for each $i=1,\dots,k.$ For any $x\geq 0$ such that $x\leq z+\sum_{i=1}^{k}A[i]$ there exists a subset $I\subseteq\{1,\dots k\}$ and $0\leq z^{(d)}\leq z$ such that

[TABLE]

Moreover, $I$ and $z^{(d)}$ can be computed in linear time.

Proof:

If $\sum_{i=1}^{k}A[i]<x,$ the desired result is given by setting $I=\{1,\dots,k\}$ and $z^{(d)}=x-\sum_{i=1}^{k}A[i]$ which is a positive number not larger than $z$, from the assumption that $z+\sum_{i=1}^{k}A[i]\geq x.$ Note that the condition can be checked in linear time.

Let us now assume that $\sum_{i=1}^{k}A[i]\geq x.$ Let $j$ be the minimum index such that $\sum_{i=1}^{j}A[j]\geq x.$ Then setting $I=\{1,\dots,j-1\}$ (if $j=1,$ we set $I=\emptyset$) and—using the assumption that $z\geq A[j]$—$z^{(d)}=x-\sum_{i=1}^{j}A[j]$ we have the desired result. Note that also in this case the index $j$ can be found in linear time. ∎

As said before, Lemma 3 is an important technical tool of our main algorithm. Therefore, in Algorithm 2 we give an efficient way to compute the value $z^{(d)}$ and the set of indices $I$.

By padding the probability distributions with the appropriate number of [math]’s, we can assume that both ${\bf p},{\bf q}\in{\cal P}_{n}$. We are now ready to present our main algorithm. The pseudocode is given in Algorithm 1. Since the description of Algorithm 1 might look complicated, we see fit to illustrate and comment its behavior with a worked out example. The reader is advised to go through the content of Section III-A before reading the formal proofs of Section III-B.

III-A How Algorithm 1 works: An informal description of its functioning and a numerical example

At any point during the execution of the algorithm, we say that ${\bf q}$ is $i$-satisfied if the sum of the entries on columns $i,i+1,\dots,n$ of the matrix the algorithm is constructing, is equal to $q_{i}+q_{i+1}+\dots+q_{n}$ Analogously, we say that ${\bf p}$ is $i$-satisfied if the sum of the entries on rows $i,i+1,\dots,n$ is equal to $p_{i}+p_{i+1}+\dots+p_{n}$. Clearly, a matrix $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ if and only if it holds that both ${\bf p}$ and ${\bf q}$ are $i$-satisfied for each $i=1,\dots,n$.

Let ${\bf z}$ be the vector defined in Fact 1, and $\mathbf{M}_{{\bf z}}$ be a matrix defined by setting $\mathbf{M}_{{\bf z}}[i,i]=z_{i}$ and setting all the other entries to zero. The basic observation is that for the matrix $\mathbf{M}_{{\bf z}}$ either ${\bf p}$ or ${\bf q}$ is $i$-satisfied, for each $i=1,\dots,n$, (but not necessarily both). In addition, every constraint which is not satisfied, coincides with an overflow, i.e., if for instance for some $i$ we have that for $\mathbf{M}_{{\bf z}}$ defined above ${\bf p}$ is not $i$-satisfied, it is necessarily the case that the sum of rows $i,i+1,\dots,n$ of $\mathbf{M}_{{\bf z}}$ is strictly greater than $p_{i}+p_{i+1}+\dots+p_{n}.$

We can understand our algorithm as working on how to modify $\mathbf{M}_{{\bf z}}$ in order to achieve $i$-satisfiability for both ${\bf p}$ and ${\bf q}$ for each $i=1,\dots,n,$ by splitting in at most two parts each diagonal element. The algorithm processes the vector ${\bf z}$ from the smallest component $z_{n}$ to the largest $z_{1}$. For $i=n,\dots,1$ it keeps $z_{i}$ in the diagonal entry $\mathbf{M}[i,i]$ as long as both ${\bf p}$ and ${\bf q}$ are $i$-satisfied.

When, e.g., ${\bf q}$ is not $i$-satisfied, it must be necessarily overflowed, i.e., the sum of the components on the $i$-th column is larger than $q_{i}.$ Then, the algorithm’s action is equivalent to removing the surplus from the $i$-th column and place it onto the column $i-1$ so that ${\bf q}$ becomes $i$-satisfied and ${\bf p}$ remains $i$-satisfied, as the mass moved is still on the same rows.

This operation can be accomplished using Lemma 3, i.e., by selecting a subset of the non-zero components on column $i$ together with $0<z_{i}^{\prime}<z_{i}$ so that their sum is equal to $q_{i}$. Keep this mass on column $i$ and move the remaining components and the left over of $z_{i}$ to column $i-1.$ In this process only $z_{i}$ gets split.

Analogously, when ${\bf p}$ is not $i$-satisfied, it must be necessarily overflowed, i.e., the sum of the components on the $i$-th row is larger than $p_{i}.$ Then, the algorithm’s action is equivalent to removing the surplus from the $i$-th row and place it onto the row $i-1$ so that ${\bf p}$ becomes $i$-satisfied and ${\bf q}$ remains $i$-satisfied, as the mass moved is still on the same columns.

This operation is again accomplished using Lemma 3: select a subset of the non-zero components on row $i$ together with $0<z_{i}^{\prime}<z_{i}$ so that their sum is equal to $p_{i}$. Keep this mass on row $i$ and move the remaining components and the left over of $z_{i}$ to row $i-1.$ Again in this process only $z_{i}$ gets split.

Let us consider the following example: Let $n=6$ and ${\bf p}=(0.4,0.3,0.15,0.08,0.04,0.03)$ and ${\bf q}=(0.44,0.18,0.18,0.15,0.03,0.02),$ be the two probability distributions for which we are seeking a coupling of minimum entropy. We have ${\bf z}={\bf p}\land{\bf q}=(0.4,0.22,0.18,0.13,0.04,0.03).$

In the first iteration, we process the entry $(6,6)$ containing $z_{6}$ (indicated in bold, below). In the matrix $\mathbf{M}_{z}$ (below) we have that ${\bf p}$ is $6$-satisfied but $q_{6}$ is overflowed. Therefore, we split $z_{6}$ into $0.2=q_{6}$ and $0.1$ and leave the former as entry $m_{6\,6}$ and the make the latter be entry $m_{6\,5},$ obtaining the matrix $\mathbf{M}^{(6)}$ on the right. The underlined values represent the mass that has been moved from one column to the next one on the left.

[TABLE]

[TABLE]

Then, we process entry $(5,5)$ containing $z_{5}$ (indicated in bold, below). In $\mathbf{M}^{(6)}$ we now have that ${\bf p}$ is $5$-satisfied but $q_{5}$ is overflowed. Therefore, we apply Lemma 3 to column 5, in order to find a split of $z_{5}$ and some of the other components of column 5 whose total sum is equal to $q_{5}$ and we move the remaining mass to column 4. Splitting $z_{5}$ into $0.2+0.2$ we obtain the matrix $\mathbf{M}^{(5)}$ on the right. The underlined values represent the mass that has been moved from one column to the next one on the left.

[TABLE]

[TABLE]

Then, we process entry $(4,4)$ containing $z_{4}$ (indicated in bold, below). In $\mathbf{M}^{(5)}$ we now have that ${\bf q}$ is $4$-satisfied but $p_{4}$ is overflowed. Therefore, we apply Lemma 3 to row 4, in order to find a split of $z_{4}$ such that one part is equal to $p_{4}$ and we move the remaining mass to row 3. Splitting $z_{4}$ into $0.8+0.5$ we obtain the matrix $\mathbf{M}^{(4)}$ on the right. The underlined values represent the mass that has been moved.

[TABLE]

[TABLE]

Then, we process entry $(3,3)$ containing $z_{3}$ (indicated in bold, below). In $\mathbf{M}^{(4)}$ we now have that ${\bf q}$ is $3$-satisfied but $p_{3}$ is overflowed. Therefore, we apply Lemma 3 to row 4, in order to find a split of $z_{3}$ and some of the other components of row 3 whose total sum is equal to $p_{3}$ and we move the remaining mass to row 2. If we split $z_{3}$ into $0.15+0.03$ we obtain the matrix $\mathbf{M}^{(3)}$ on the right. The underlined values represent the mass that has been moved from one row to the next one above.

[TABLE]

[TABLE]

Then, we process entry $(2,2)$ containing $z_{2}$ (indicated in bold, below). In $\mathbf{M}^{(3)}$ we now have that ${\bf p}$ is $2$-satisfied but $q_{2}$ is overflowed. Therefore, we apply Lemma 3 to column 2, in order to find a split of $z_{2}$ and some of the other components of column 2 whose total sum is equal to $q_{2}$ and we move the remaining mass to column 1. If we split $z_{2}$ into $0.18+0.04$ we obtain the matrix $\mathbf{M}^{(2)}$ on the right. The underlined values represent the mass that has been moved from one column to the next one on the left.

[TABLE]

[TABLE]

Finally, we process entry $(1,1)$ containing $z_{1}$ (indicated in bold, below). In $\mathbf{M}^{(2)}$ we have that both ${\bf p}$ and ${\bf q}$ are $1$-satisfied. Therefore we get the unmodified matrix $\mathbf{M}^{(1)}$ on the right which is our joint distribution. Notice that each component of ${\bf z}$ has been split at most into two parts. In particular only when $z_{i}$ is processed the first time it might get split, while the other components (obtained by the previous subdivision of some other components of ${\bf z}$) might be relocated but not chunked again.

[TABLE]

[TABLE]

III-B *The proof of correctness of *Algorithm 1

The following theorem shows the correctness of Algorithm 1. In particular, the equalities in (13), for the case $i=1$, imply that the matrix built by the algorithm is a coupling of ${\bf p}$ and ${\bf q}$.

Theorem 1**.**

For each $i=n,n-1,\dots,1$ at the beginning of iteration $i$ of the main while loop the following conditions hold

for each $i^{\prime}\leq i$ we have

[TABLE] 2. 2.

exactly one of the following holds

1. (a)

$\sum_{k=1}^{n}m_{k\,i}=q_{i}$* and $\sum_{k=1}^{n}m_{i\,ki}=p_{i}$* 2. (b)

$\sum_{k=1}^{n}m_{k\,i}>q_{i}$* and $\sum_{k=1}^{n}m_{i\,ki}=p_{i}$* 3. (c)

$\sum_{k=1}^{n}m_{k\,i}=q_{i}$* and $\sum_{k=1}^{n}m_{i\,ki}>p_{i}$ hence, at most one of the if conditions is true.*

Moreover at the end of the iteration $i,$ for each $i^{\prime}\geq i$ it holds that

[TABLE]

Proof:

We prove the statement by reverse induction. For $i=n$, due to the initialization in lines 1 and 4 we have that for each $i^{\prime}\leq n$ the only non-zero entry in row $i^{\prime}$ and in column $i^{\prime}$ is $m_{i^{\prime}\,i^{\prime}}=z_{i^{\prime}}$ and 1 holds.

By definition we have $m_{n\,n}=z_{n}=\max\{p_{n},q_{n}\}=p_{n},$ since, by the initialisation in line 2 we can assume $p_{n}\geq q_{n}$. Therefore, either $m_{n\,n}=p_{n}=q_{n}$ and 2a) holds; or $m_{n\,n}=p_{n}>q_{n}$ and 2b) holds. Thus, 2) holds.

Finally, if 2a) holds during the iteration $n$ then no modification of the matrix entries is performed and at the end of the iteration equation (13) holds. Otherwise, as already observed, because of the initialization in line 2, we have $z_{n}=p_{n}>q_{n}$. Then, as a result of the call to procedure Lemma3($z_{n},q_{n},A=\emptyset$) we will have $z^{(d)}_{n}=q_{n},\,z^{(r)}_{n}=p_{n}-q_{n}$ and the modifications to the matrix entries $m_{n,n}=z^{(d)}_{n}=q_{n}$ and $m_{n,n-1}=z^{(r)}_{n}=p_{n}-q_{n},$ from which (13) holds at the end of the iteration as desired. This settles the induction base.

Let us now assume that the claims hold for iteration $i+1$ and prove it for iteration $i$.

1. By induction hypothesis, at the beginning of iteration $i+1$ for each $i^{\prime}\leq i+1$, hence in particular for each $i^{\prime}\leq i$ it holds that

[TABLE]

During iteration $i+1$ the only possible changes to entries of the matrix are either in rows $i$ and $i+1$ (when the if at line 7 is satisfied) or in columns $i$ and $i+1$ (when the if at line 13 is satisfied). Moreover, such modifications do not change the total probability mass in rows $i$ and $i+1$ and the total probability mass in column $i$ and $i+1$, i.e., the sums $\sum_{k=1}^{n}(m_{k\,i}+m_{k\,i+1})$ and $\sum_{k=1}^{n}(m_{i\,k}+m_{i+1\,k})$ remain unchanged during iteration $i+1$. It follows that at the beginning of iteration $i$ equality (14) still holds for each $i^{\prime}\leq i.$ This settles the inductive steps for property 1.

2. By induction hypothesis from 1 with $i^{\prime}=i<i+1$ we have

[TABLE]

By induction hypothesis, we also have that for each $\ell=i+1,\dots,n,$

[TABLE]

From equations (15)-(16) together with (1) we get

[TABLE]

and

[TABLE]

Therefore, (a) if $\sum_{\ell=i}^{n}z_{\ell}=\sum_{\ell=i}^{n}q_{\ell}=\sum_{\ell=i}^{n}p_{\ell}$ then from (17) we have $p_{i}=\sum_{k=1}^{n}m_{i\,k}$ and from (18) we have $q_{i}=\sum_{k=1}^{n}m_{k\,i}.$

(b) If $\sum_{\ell=i}^{n}z_{\ell}=\sum_{\ell=i}^{n}p_{\ell}>\sum_{\ell=i}^{n}q_{\ell}$ then from (17) we have $p_{i}=\sum_{k=1}^{n}m_{i\,k}$ and from (18) we have $q_{i}<\sum_{k=1}^{n}m_{k\,i}.$

(c) If $\sum_{\ell=i}^{n}z_{\ell}=\sum_{\ell=i}^{n}q_{\ell}>\sum_{\ell=i}^{n}p_{\ell}$ then from (17) we have $p_{i}<\sum_{k=1}^{n}m_{i\,k}$ and from (18) we have $q_{i}=\sum_{k=1}^{n}m_{k\,i}.$

Exactly one of these three cases is possible, which proves the induction step for 2).

Let us now prove the induction step for (13).

Assume first that during the iteration $i$ case 2a) applies. From the previous point, this means that $\sum_{\ell=i}^{n}z_{\ell}=\sum_{\ell=i}^{n}q_{\ell}=\sum_{\ell=i}^{n}p_{\ell}$. Then, none of the two if compounds (lines 7-12 and lines 13-18) are executed and no matrix entry is changed in this iteration. As a result at the end of the iteration we have that, for each $i^{\prime}>i$ the formula (13) holds by induction hypothesis. Moreover it also holds for $i^{\prime}=i$ since, from 1) and 2a) we have

[TABLE]

and analogously

[TABLE]

Assume now that during the iteration $i$ case 2b) applies (the case 2c) can be dealt with symmetrically). Then, the if compound in lines 7-12 is executed. As a result, values $z^{(d)}_{i},z^{(r)}_{i}$ and set $I\subseteq[n]$ are computed such that $z^{(d)}_{i}+z^{(r)}_{i}=z_{i}$ and $z^{(d)}_{i}+\sum_{k\in I}m_{k\,i}=q_{i}$. Before the assignments in line 9 and the execution of the for loop, we had that

[TABLE]

After the assignments in line 9 and the execution of the for loop, the mass in the last two terms is moved to column $i-1$, but without changing the row. Therefore the row sums do not change and the column sums of column $i$ and $i-1$ change so that $\sum_{k=1}^{n}m_{k\,i}=z^{(d)}_{i}+\sum_{k\in I}m_{k\,i}=q_{i}$ as desired.

Finally, it is possible that during the iteration $i$ case 2c) applies. The analysis for this case is analogous to the one for the previous case, from which it can be easily obtained by symmetry swapping the roles of rows and columns, taking into account that we have to consider the result of the operations executed within the if compound in lines 13-18.

The proof is complete. ∎

III-C *The guaranteed additive gap of *Algorithm 1

We are now ready to formally prove our main result.

Theorem 2**.**

For any ${\bf p},{\bf q}\in{\cal P}_{n}$, Algorithm 1 outputs in polynomial time an $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ such that

[TABLE]

Proof:

It is not hard to see that the values of all non-zero entry of the matrix $\mathbf{M}$ are initially set in line 4 and then in lines 9 and 15—in fact, the assignments in lines 11-12 and 17-18 have the effect of shifting by one column to the left or by one row up values that had been fixed at some point earlier in lines 9 and 15. Therefore, all the final non-zero entries of $\mathbf{M}$ can be partitioned into $n$ pairs $z_{j}^{(r)},z_{j}^{(d)}$ with $z_{j}^{(r)}+z_{j}^{(d)}=z_{j}$ for $j=1,\dots,n$. By using the standard assumption $0\log\frac{1}{0}=0$ and applying Jensen inequality we have

[TABLE]

which concludes the proof of the bound on the additive gap guaranteed by Algorithm 1. Moreover, one can see that Algorithm 1 can be implemented so to run in $O(n^{2})$ time. For the time complexity of the algorithm we observe the following easily verifiable fact:

• •

the initialization in line 1 takes $O(n^{2})$;

• •

the condition in line 2 can be easily verified in $O(n)$ which is also the complexity of swapping ${\bf p}$ with ${\bf q}$, if needed;

• •

the vector $z={\bf p}\land{\bf q}$ can be computed in $O(n)$, once the suffix sums $\sum_{j=k}^{n}p_{j},\,\sum_{j=k}^{n}q_{j}$ ($k=1,\dots,n$) have been precomputed (also doable in $O(n)$);

• •

the main while loop is executed $n$ times and all the operations executed in an iteration are easily upper bounded by $O(n)$. The most expensive are the calls to the procedure Lemma3, and the for-loops. All these take $O(n)$. Therefore, the overall running time of the while loop is also $O(n^{2})$.

Therefore we can conclude that the time complexity of Algorithm 1 is $O(n^{2}),$ hence polynomial time. ∎

III-D Improving the time complexity

We note that the time complexity of Algorithm 1 can be improved if we build the coupling ${\bf M}$ in sparse form, i.e., as the set of values $\{({\bf M}[i,j],(i,j))\mid{\bf M}[i,j]\neq 0\}$ containing only the non-zero entries of ${\bf M}$ together with their coordinates.

We can keep the moved masses, i.e., the pieces $z_{i}^{(r)}$ that are iteratively moved from one column to the previous one, in line 11 (respectively, from one row to the previous one, in line 17) in a priority queue ${\cal Q}$. For each such element we store in the priority queue its row index (resp. column index) and its mass. With a standard implementation of a priority queue, we can then efficiently find the value of the minimum mass stored in constant time $O(1)$ (we refer to this operation as Min(${\cal Q}$) and extract the minimum mass in time logarithmic in the number of elements stored in the priority queue (we refer to this operation as ExtractMin(${\cal Q}$))[10]. Accordingly, procedure Lemma 3 amounts to iteratively extract from the priority queue the smallest mass as long as the queue is not empty and the sum of the masses extracted do not overcomes $\min\{p_{i},q_{i}\}$ . Whenever we split $z_{i}$, we insert $z_{i}^{(r)}$ into the priority queue (this operation can also be implemented to require time logarithmic in the size of the queue; we refer to it as Insert(${\cal Q},(z_{i}^{(r)},i)$)).

At any time, the priority queue will contain $O(n)$ elements. Therefore, each insertion (Insert) and extraction (ExtractMin) from the priority queue takes $O(\log n)$ time. Moreover, since each element enters the queue at most once, the overall time of all insertion and extraction operations is upper bounded by $O(n\log n)$. The remaining part of the algorithm takes $O(n)$, apart from the possible initial sorting of the two distribution, adding another $O(n\log n)$ term. Therefore, the resulting implementation has complexity $O(n\log n)$.

We report in appendix a pseudocode of such implementation, where, for the sake of a clearer description, we use two priority queues, ${\cal Q}^{(row)},{\cal Q}^{(col)},$ storing masses moved among rows and masses moved among columns respectively.

IV Extending the Results to other Entropy Measures

Our approach to prove entropic inequalities via majorization theory seems quite powerful. Indeed, it allows us to extend our results to different kind of entropies, with no additional effort. As an example, let us consider the order $\alpha$ Rényi entropy [45] of a probability distribution ${\bf p}=(p_{1},\ldots,p_{n})$, defined as

[TABLE]

where $\alpha\in(0,1)\cup(\,\infty)$. It is well known that the Rényi entropy is Schur-concave, for all the values of the parameter $\alpha$ [25]. Therefore, we immediately have the analogous of Lemma 2 of Section II.

Lemma 4**.**

For any* ${\bf p}$ and ${\bf q}$, for any $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ and $\alpha\in(0,1)\cup(\,\infty)$, it holds that*

[TABLE]

We now prove the analogous of Theorem 2, that is

Theorem 3**.**

For any ${\bf p},{\bf q}\in{\cal P}_{n}$, Algorithm 1 outputs in polynomial time an $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ such that

[TABLE]

Proof:

Let $\mathbf{M}$ be the matrix constructed by our Algorithm 1, and let $\alpha\in(0,1)$. Proceedings as in Theorem 2 (and with the same notations), we have:

[TABLE]

The proof for the case $\alpha\in(1,\infty)$ is the same, by noticing that for $\alpha>1$ the Jensen inequality goes into the opposite direction and that $1/(1-\alpha)<0$. ∎

One can extend our results also to other entropies, like the Tsallis entropy [51], using its Schur-concavity property proved in [20]. The mathematical details can be easily worked out by the motivated reader.

V An Extension to Multivariate Distributions

In this section we will show how the algorithm Min-Entropy-Joint-Distr can be used to guarantee an additive gap at most $\log k$ for the problem of minimizing the entropy of a joint distribution, with marginals equal to $k$ given input distributions, for any $k\geq 2.$

In what follows, for the ease of the description, we shall assume that $k=2^{\kappa}$ for some integer $\kappa\geq 1,$ i.e., $k$ is a power of $2$. A pictorial perspective on the algorithm’s behaviour is to imagine that the input distributions are in the leaves of a complete binary tree with $k=2^{\kappa}$ leaves. Each internal node $\nu$ of the tree contains the joint distribution of the distributions in the leaves of the subtree rooted at $\nu$. Such a distribution is computed by applying the algorithm Min-Entropy-Joint-Distr to the distributions in the children of $\nu$.

The algorithm builds such a tree starting from the leaves. Thus, the joint distribution of all the input distributions will be given by the distribution computed at the root of the tree.

Let us first see how in the case of four input distributions, our algorithms does indeed guarantee that the final matrix is a joint distribution with marginals equal to the input distributions. An high-level pictorial way to describe how our algorithm operates is given in Figure 1.

Given probability distributions ${\bf p}$ and ${\bf q}$, of dimension $n_{{\bf p}}$ and $n_{{\bf q}}$, respectively, the algorithm Min-Entropy-Joint-Distr described in Section III produces a matrix $M_{{\bf p}\,{\bf q}}$ such that for fixed $\overline{i}$ and $\overline{j}$ it holds that

[TABLE]

We also have that, for each $i_{1}\neq i_{2}$ the set of entries of $M_{{\bf p}\,{\bf q}}$ whose sum is $p_{i_{1}}$ is disjoint from the set of entries of $M_{{\bf p}\,{\bf q}}$ whose sum is $p_{i_{2}}$. Analogously, for each $j_{1}\neq j_{2}$ the set of entries of $M_{{\bf p}\,{\bf q}}$ whose sum is $q_{j_{1}}$ is disjoint from the set of entries of $M_{{\bf p}\,{\bf q}}$ whose sum is $q_{j_{2}}$.

Let us define the probability distribution ${\bf x}=(x_{1},x_{2},\ldots)$ whose components are all and only the non-zero entries of $M_{{\bf p}\,{\bf q}}$ sorted in non-increasing order. Let $n_{{\bf x}}$ denote the number of components of ${\bf x}$. For each $a=1,\dots n_{{\bf x}}$ let us fix a one-one mapping $a\leftrightarrow(i,j)$ recording the fact that $x_{a}=M_{{\bf p}\,{\bf q}}[i,j]$.

Consider now another pair of distributions ${\bf r},{\bf t}$, of dimension $n_{{\bf r}}$ and $n_{{\bf t}}$, respectively. Applying Algorithm Min-Entropy-Joint-Distr to ${\bf r},{\bf t}$ we obtain a matrix $M_{{\bf r}\,{\bf t}}$ such that for fixed $\overline{i^{\prime}}$ and $\overline{j^{\prime}}$ it holds that

[TABLE]

As before, for each $i^{\prime}_{1}\neq i^{\prime}_{2}$ the set of entries of $M_{{\bf r}\,{\bf t}}$ whose sum is $p_{i^{\prime}_{1}}$ is disjoint from the set of entries of $M_{{\bf r}\,{\bf t}}$ whose sum is $p_{i^{\prime}_{2}}$. Also, for each $j^{\prime}_{1}\neq j^{\prime}_{2}$ the set of entries of $M_{{\bf r}\,{\bf t}}$ whose sum is $q_{j^{\prime}_{1}}$ is disjoint from the set of entries of $M_{{\bf r}\,{\bf t}}$ whose sum is $q_{j^{\prime}_{2}}$.

Let ${\bf y}=(y_{1},y_{2},\ldots)$ be the probability distribution whose components are all and only the non-zero entries of $M_{{\bf r}\,{\bf t}}$ sorted in non-increasing order. Let $n_{{\bf y}}$ denote the number of components of ${\bf y}$. For each $b=1,\dots n_{{\bf y}}$ let us fix a one-one mapping $b\leftrightarrow(i^{\prime},j^{\prime})$ recording the fact that $y_{b}=M_{{\bf r}\,{\bf t}}[i^{\prime},j^{\prime}]$.

If we now apply algorithm Min-Entropy-Joint-Distr on the distributions ${\bf x}$ and ${\bf y}$ we get a matrix $M_{{\bf x}\,{\bf y}}$ such that for fixed $\overline{a}$ and $\overline{b}$ it holds that

[TABLE]

Let us now define a new $4$-dimensional array $\tilde{M}[i,j,i^{\prime},j^{\prime}]$ by stipulating that for each $i\in[n_{{\bf p}}],\,j\in[n_{{\bf q}}],\,i^{\prime}\in[n_{{\bf r}}],\,j^{\prime}\in[n_{{\bf t}}],$ the following equalities hold

[TABLE]

Then, applying the properties above, for each $\overline{i}\in[n_{{\bf p}}],$ we have that

[TABLE]

where the equality in (26) follows by restricting the sum over the non-zero entries of $\tilde{M}$; (27) follows by the definition of $\tilde{M}$; (28) follows by (24); the first part of (29) follows by the fact that the components of ${\bf x}$ coincide with non zero entries of $M_{{\bf p}\,{\bf q}}$; the first equality in (29) follows since we are adding to the previous term only entries $M_{{\bf p}\,{\bf q}}[\overline{i},j]=0$; finally the last equality follows from (23).

Proceeding in the same way we can show that for each $\overline{j},\overline{i^{\prime}},\overline{j^{\prime}}$ we have

[TABLE]

hence concluding that $\tilde{M}$ is a joint distribution with marginals equal to ${\bf p},{\bf q},{\bf r},{\bf t},$ as desired.

Completing the argument. We can now inductively extend the above argument to the case of more distributions. Assume that we have produced an array ${\bf M}_{{\bf p}^{(1)},\dots,{\bf p}^{(r)}}$ which is a joint distribution with marginals equal to ${\bf p}^{(1)},\dots,{\bf p}^{(r)}$, of dimension $n_{1},\dots,n_{r}$ respectively. Analogously, let us assume that we have produced ${\bf M}_{{\bf q}^{(1)},\dots,{\bf q}^{(s)}}$ which is a joint distribution with marginals equal to ${\bf q}^{(1)},\dots,{\bf q}^{(s)},$ of dimension $m_{1},\dots,m_{s}$ respectively. This means that for each $\ell=1,\dots,r$ and $1\leq i\leq n_{\ell}$ and for each $\ell^{\prime}=1,\dots,r$ and $1\leq j\leq n_{\ell^{\prime}}$ we have that

[TABLE]

Proceeding as before, let us define the probability distribution ${\bf x}$ whose components are all and only the non-zero entries of $M_{{\bf p}^{(1)},\dots,{\bf p}^{(r)}}$ sorted in non-increasing order. Let $n_{{\bf x}}$ denote the number of components of ${\bf x}$. For each $a=1,\dots n_{{\bf x}}$ let us fix a one-one mapping $a\leftrightarrow(i_{1},\dots,i_{r})$ recording the fact that $x_{a}=M_{{\bf p}^{(1)},\dots,{\bf p}^{(r)}}[i_{1},\dots,i_{r}]$.

Let ${\bf y}$ be the probability distribution whose components are all and only the non-zero entries of $M_{{\bf q}^{(1)},\dots,{\bf q}^{(s)}}$ sorted in non-increasing order. Let $n_{{\bf y}}$ denote the number of components of ${\bf y}$. For each $b=1,\dots n_{{\bf y}}$ let us fix a one-one mapping $b\leftrightarrow(j_{1},\dots,j_{s})$ recording the fact that $y_{b}=M_{{\bf q}^{(1)},\dots,{\bf q}^{(s)}}[j_{1},\dots,j_{s}]$.

Applying algorithm Min-Entropy-Joint-Distr on the distributions ${\bf x}$ and ${\bf y}$ we get a matrix $M_{{\bf x}\,{\bf y}}$ such that for fixed $\overline{k}$ and $\overline{\ell}$

[TABLE]

Therefore, we can define a new $r+s$-dimensional array $\tilde{M}[i_{1},\dots,i_{r},j_{1},\dots,j_{s}]$ by stipulating that for each $i_{1},\dots,i_{r}$ such that $i_{\ell}\in[n_{k}]$ for $\ell=1,\dots,r$ and for any $j_{1},\dots,j_{s}$ such that $j_{\ell^{\prime}}\in[m_{\ell^{\prime}}]$ for $\ell^{\prime}=1,\dots,s,$

[TABLE]

It is not hard to see that proceeding like in (26)-(29) one can show that $\tilde{M}$ is indeed a joint distribution with marginals equal to ${\bf p}^{(1)},\dots,{\bf p}^{(r)},{\bf q}{(1)},\dots,{\bf q}^{(s)}.$

V-A *The pseudocode: Algorithm 3 *

Algorithm 3 shows the pseudocode for our procedure. We denote by $m^{(i-j)}$ the non-zero components of the distribution that our algorithm builds as joint distribution of ${\bf p}^{(i)},{\bf p}^{(i+1)},\dots,{\bf p}^{(j)}.$

The vector $Ind^{(i-j)}$ is used to record for each component ${\bf m}^{(i-j)}[w]$ the indices of the component of the joint probability distribution of ${\bf p}^{(i)},\dots,{\bf p}^{(j)}$ which coincides with ${\bf m}^{(i-j)}[w]$. With respect to the description above $Ind^{(i-j)}$ is used to record the one-one mapping between the elements of ${\bf m}^{(i-j)}$ and the non-zero elements of the joint distribution of ${\bf p}^{(i)},{\bf p}^{(i+1)},\dots,{\bf p}^{(j)}.$ Therefore, in accordance to the above arguments, after the execution of line 17, for $w=1,\dots,|{\bf m}^{(i-j)}|$, if, e.g., $Ind^{(i-j)}[w]=\langle s_{i}[w],s_{i+1}[w],\dots,s_{j}[w]\rangle$ it means that setting $M^{(i-j)}[s_{i}[w],s_{i+1}[w],\dots,s_{j}[w]]\leftarrow{\bf m}^{(i-j)}[w]$ and setting the remaining components of $M^{(i-j)}$ to zero, the array $M^{(i-j)}$ is a joint distribution matrix for ${\bf p}^{(i)},\dots,{\bf p}^{(j)}$ whose non-zero components are equal to the components of ${\bf m}^{(i-j)}.$ Hence, in particular, we have that $H(M^{(i-j)})=H({\bf m}^{(i-j)}).$

Note that the algorithm explicitly uses this correspondence only for the final array $M^{(1-k)}$ representing the joint distribution of all input distributions. Based on the above discussion the correctness of the algorithm can be easily verified.

V-B The additive gap guaranteed by K-Min-Entropy-Joint-Distribution

In this section we will prove that the entropy of the joint distribution output by the algorithm guarantees an additive gap at most $\log k.$

We will prepare some definitions and lemmas which will be key tools for proving the approximation guarantee of our algorithm. The proof of these technical lemmas is deferred to the next section.

Let us define the following:

Definition 2**.**

For any ${\bf p}=(p_{1},\dots,p_{n})\in{\cal P}_{n}$ we denote by ${\sf half}({\bf p})$ the distribution $(\frac{p_{1}}{2},\frac{p_{1}}{2},\frac{p_{2}}{2},\frac{p_{2}}{2},\dots,\frac{p_{n}}{2},\frac{p_{n}}{2})$ obtained by splitting each component of ${\bf p}$ into two identical halves.

For any $i\geq 2,$ let us also define ${\sf half}^{(i)}({\bf p})={\sf half}({\sf half}^{(i-1)}({\bf p})),$ where ${\sf half}^{(1)}({\bf p})={\sf half}({\bf p})$ and ${\sf half}^{(0)}({\bf p})={\bf p}.$

We will employ the following two technical lemmas whose proofs are in the next section.

Lemma 5**.**

For any ${\bf p}\preceq{\bf q}$ we have also ${\sf half}({\bf p})\preceq{\sf half}({\bf q})$

Lemma 6**.**

For any pair of distributions ${\bf p},{\bf q}\in{\cal P}_{n}.$ and any $i\geq 0$, It holds that

[TABLE]

Theorem 4**.**

For each $\ell=0,1,\dots\kappa$ and $s=0,1,2,\dots,k/2^{\ell}-1$ let $i=i(\ell,s)=s\cdot 2^{\ell}+1$ and $j=j(\ell,s)=(s+1)\cdot 2^{\ell}=i+2^{\ell}-1.$ Then, we have

[TABLE]

Proof:

The proof is by induction on $\ell.$ The base case follows by definition of the operator ${\sf half}^{(\ell)}$ and the fact that the algorithm sets ${\bf m}^{(i-i)}={\bf p}^{(i)},$ for each $i$ hence in particular ${\bf m}^{(i-i)}={\bf p}^{(i)}={\sf half}^{(0)}({\bf p}^{(i)}),$ which proves the desired inequality.

We now prove the induction step. Let $\ell>0.$ It is enough to consider only the case $s=0,$ since the other cases are perfectly analogous.

Therefore, $i=1$ and $j=2^{\ell}.$ Using the notation employed in the pseudocode, let $j_{1}=2^{\ell-1},\,j_{2}=2^{\ell-1}+1.$ By induction hypothesis we can assume that

[TABLE]

[TABLE]

It follows that

[TABLE]

where

• •

(V-B) follows from (36) by the definition of the operator ${\sf half}$;

• •

(LABEL:eq:half3) follows from (V-B) by Lemma 6;

• •

(39) follows from (LABEL:eq:half3) by the induction hypotheses (34)-(35) ;

• •

(40) follows from (39) by observing that the components of ${\bf m}^{(i-j)}$ coincide with the components of the array $M$ output by algorithm Min-Entropy-Joint-Distribution executed on the distributions ${\bf m}^{(i-j_{1})}$ and ${\bf m}^{(j_{2}-j)}.$ Let ${\bf z}={\bf m}^{(i-j_{1})}\wedge{\bf m}^{(j_{2}-j)}$ and $|{\bf z}|$ denote the number of components of ${\bf z}.$ By the analysis presented in the previous section we have that we can partition the components of $M$ (equivalently, the components of ${\bf m}^{(i-j)}$) into subsets $M_{1},M_{2},\dots,M_{|{\bf z}|}$ such that

• –

$1\leq|M_{i}|\leq 2$

• –

for each $i=1,\dots,|{\bf z}|,$ it holds that $\sum_{x\in M_{i}}x=z_{i}$;

Therefore—assuming, w.l.o.g., that the components of ${\bf m}^{(i-j)}$ are reordered such that those in $M_{i}$ immediately precede those in $M_{i+1}$—we have ${\sf half}({\bf z})={\bf m}^{(i-j)}P$ where $P=[p_{i\,j}]$ is a doubly stochastic matrix defined by

[TABLE]

from which it follows that ${\sf half}({\bf z})\preceq{\bf m}^{(i-j)}$ yielding (40). ∎

An immediate consequence of the last theorem is the following

Corollary 1**.**

For any $k$ probability distributions ${\bf p}^{(1)},\dots,{\bf p}^{(k)}$ let $M$ be the joint distribution, with marginals equal to ${\bf p}^{(1)},\dots,{\bf p}^{(k)}$, output by algorithm K-Min-Entropy-Joint-Distribution. Then,

[TABLE]

Proof:

Let $k$ be a power of $2$. Otherwise, one can duplicate some of the probability distributions until there are $k^{\prime}=2^{\lceil\log k\rceil}$ of them. By Theorem 4 we have

[TABLE]

Therefore, by the Schur-concavity of the entropy we have

[TABLE]

where the last equality follows by the simple observation that for any probability distribution ${\bf x}$ and integer $i\geq 0$ we have $H({\sf half}^{(i)}({\bf x}))=H({\bf x})+i.$ ∎

We also have the following lower bound which, together with the previous corollary implies that our algorithm guarantees an additive gap of at most $\log k$ bits for the problem of computing the joint distribution of minimum entropy of $k$ input distributions.

Lemma 7**.**

Fix $k$ distributions ${\bf p}^{(1)},{\bf p}^{(2)},\cdots,{\bf p}^{(k)}$. For any $M$ being a joint distribution with marginals ${\bf p}^{(1)},{\bf p}^{(2)},\cdots,{\bf p}^{(k)}$, it holds that

[TABLE]

Proof:

For each $i=1,\dots,k,$ the distribution ${\bf p}^{(I)}$ is an aggregation of $M$, hence $M\preceq{\bf p}^{(i)}.$

By definition of the greatest lower bound operator $\wedge$ for any distribution ${\bf x}$ such that for each $i$ it holds that ${\bf x}\prec{\bf p}^{(i)}$ we have ${\bf x}\preceq{\bf p}^{(1)}\wedge{\bf p}^{(2)}\wedge\cdots{\bf p}^{(k)}$. Therefore, in particular we have $M\preceq{\bf p}^{(1)}\wedge{\bf p}^{(2)}\wedge\cdots{\bf p}^{(k)},$ which, by the Schur concavity of the entropy gives the desired result. ∎

Remark 1**.**

The time complexity of Algorithm 3 is dominated by the time to build the output matrix in line 19, which takes $O(n^{k}).$ However, if the output matrix is returned in sparse form and in line 9 the improved implementation of algorithm Min-Entropy-Joint-Distribution is used (see section III-D and the appendix), the time for the overall construction is upper bounded by $\sum_{\ell=1}^{\log k}O(\frac{k}{2^{\ell}}2^{\ell-1}n\log(2^{\ell-1}n))=O(nk\log(nk)).$ To see this, observe that the main for loop in line 5 is executed $O(\log k)$ times and in each iteration $\ell=1,\dots,\log k$ there are $\frac{k}{2^{\ell}}$ executions of Min-Entropy-Joint-Distribution over distributions having $O(2^{\ell-1}n)$ non-zero entries and the algorithm employs the arrays $Ind$ in order to perform the computation only considering the non-zero entries of these distributions.

Summarising we have shown the following

Theorem 5**.**

Let ${\bf p}^{(1)},\dots,{\bf p}^{(m)}\in{\cal P}_{n}$. Let $M^{*}$ be a joint distribution with marginals ${\bf p}^{(1)},\dots,{\bf p}^{(m)}$ of minimum entropy among all the joint distribution having marginals equal to ${\bf p}^{(1)},\dots,{\bf p}^{(m)}.$ Let $M$ be the joint distribution of ${\bf p}^{(1)},\dots,{\bf p}^{(m)}$ output by our algorithm. Then,

[TABLE]

Hence, our (polynomial time) algorithm guarantees an additive gap of $\log(m).$

Using Lemma 7 and Theorem 5 above, and Theorem 3 of [31], we obtain the following version of the functional representation lemma (please see the discussion in Section I-B of the present paper).

Corollary 2**.**

Let $X$ and $Y$ be two arbitrary random variables with joint distribution $p(x,y)$, where $X$ takes values $x_{1},\ldots,x_{k}$. Let ${\bf p}^{(1)},\dots,{\bf p}^{(k)}$ be the distribution of the conditioned r.v. $Y|X=x_{1},\ldots,Y|X=x_{k}$, respectively. Then, for any r.v. $Z$ independent from $X$ for which there exist a function $f$ such that $(X,Y)=(X,f(X,Z))$, it holds that

[TABLE]

Conversely, there exists a $Z$ independent from $X$ and a function $f$ for which $(X,Y)=(X,f(X,Z))$ such that

[TABLE]

V-C The proofs of the two technical lemmas about the operator ${\sf half}$

Lemma 5. For any ${\bf p}\preceq{\bf q}$ we have also ${\sf half}({\bf p})\preceq{\sf half}({\bf q})$.

Proof:

It is easy to see that assuming ${\bf p}$ and ${\bf q}$ rearranged in order to have $p_{1}\geq p_{2}\geq\dots\geq p_{n}$ and $q_{1}\geq q_{2}\geq\dots\geq q_{n}$ we also have ${\sf half}({\bf p})_{1}\geq{\sf half}({\bf p})_{2}\geq\dots\geq{\sf half}({\bf p})_{2n}$ and ${\sf half}({\bf q})_{1}\geq{\sf half}({\bf q})_{2}\geq\dots\geq{\sf half}({\bf q})_{2n}.$

By assumption we also have that for each $j=1,\dots,n$ it holds that $\sum_{i=1}^{j}p_{i}\leq\sum_{i=1}^{j}p_{i}.$

Therefore, for each $j=1,\dots 2n$ it holds that

[TABLE]

proving that ${\sf half}({\bf p})\preceq{\sf half}({\bf q}).$ ∎

Fact 2**.**

For any pair of distributions ${\bf p},{\bf q}\in{\cal P}_{n}.$ It holds that

[TABLE]

Proof:

By Lemma 5 we have that

[TABLE]

Then, by the property of the operator $\wedge$ which gives the greatest lower bound we have the desired result. ∎

On the basis of Fact 2 we can extend the result to ”powers” of the operator ${\sf half}$ and have our Lemma 6.

Lemma 6. For any pair of distributions ${\bf p},{\bf q}\in{\cal P}_{n}.$ and any $i\geq 0$, It holds that

[TABLE]

Proof:

We argue by induction on $i$. The base case $i=1$ is given by the previous Fact 2. Then, for any $i>1$

[TABLE]

from which the desired result immediately follows. The first $\preceq$-inequality follows by induction hypothesis and the second inequality by Fact 2. ∎

VI Conclusions

In this paper we have studied the problem of finding a minimum entropy joint distribution with fixed marginals. We have pointed out that this problem naturally arises in a variety of situations: causal inference, one-shot channel simulation, metric computation for dimension reduction, optimization in the transportation polytope, and several others. Our main result consists in a polynomial time algorithm to find an ${\bf M}\in{\cal C}({\bf p},{\bf q})$ such that $H({\bf M})\leq OPT+1$ bit, where $OPT=\min_{\mathbf{N}\in{\cal C}({\bf p},{\bf q})}H(\mathbf{N}).$ We ave also shown that our approach (relying on majorization among probability distributions) allows us to easily extend our results to Rényi entropies of arbitrary positive orders (thus generalizing the result for the Shannon entropy where the latter is equal to the Rényi entropy of order 1).

There are many possible extensions of our work. Firstly, although our result for the minimum entropy bivariate joint distribution with fixed two marginals seems quite tight, it is very likely that a more direct approach (i.e., that does not rely on the iterative construction of Section V) could give better results for multivariate joint distributions. Another interesting problem would be to extend our results to the case in which one seeks a minimum entropy bivariate joint distribution with marginals “close” to given ones, for appropriate measures of closeness. Finally, a natural research problem is related to the scenario considered in Section I-B: Given arbitrary correlated r.v.’s $X$ and $Y$, it would be interesting to find a r.v. $Z$, independent from $X$, such that the pair of r.v.’s $(X,f(X,Z))$ is distributed like $(X,Y)$, for appropriate deterministic function $f$, for which both $H(Z)$ and $H(Y|Z)$ are close to their lower bounds.

Acknowledgments

The authors want to thank Executive Editor Professor I. Sason, Associate Editor Professor I. Kontoyiannis, and the anonymous referees for many useful comments and suggestions.