How to Find a Joint Probability Distribution of Minimum Entropy (almost) given the Marginals

TL;DR

This paper addresses the challenge of approximating the joint distribution with minimum entropy given marginals, proposing an efficient algorithm with bounded additive error, and extends the approach to multiple variables.

Contribution

It introduces an approximation algorithm for the NP-hard problem of minimum-entropy joint distribution, with guarantees on the additive error, and generalizes to multiple variables.

Findings

Provides an efficient approximation algorithm with additive factor 1 for two variables.

Extends the approximation method to $k$ variables with additive factor $ ext{log} k$.

Discusses applications where minimum-entropy joint distributions are relevant.

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 the joint probability distribution in ${\cal C}({\bf p}, {\bf q})$ of minimum entropy (equivalently, the joint probability distribution that 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 by at most 1, thus providing an approximation algorithm with additive approximation factor of 1. 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 $X_1, \ldots , X_k$. In this case, our approximation algorithm has an additive approximation factor of $\log k$. We also discuss some related applications of our findings.

Equations

$${\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},\text{mbox}andforall\mathbf{z}\in {\mathcal{P}}_n\text{mbox}forwhich\mathbf{x}⪯\mathbf{z},\mathbf{y}⪯\mathbf{z}\text{mbox}itholdsthat\mathbf{x}\vee \mathbf{y}⪯\mathbf{z}.$$

$$\mathbf{x}\wedge \mathbf{y}⪯\mathbf{x},\mathbf{x}\wedge \mathbf{y}⪯\mathbf{y},\text{mbox}andforall\mathbf{z}\in {\mathcal{P}}_n\text{mbox}forwhich\mathbf{z}⪯\mathbf{x},\mathbf{z}⪯\mathbf{y}\text{mbox}itholdsthat\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}.$$

$${\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}).$$

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

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

$$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(\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.$$

$$\sum _{k=i}^n{z}_k=\sum _{k=i}^n{p}_k$$

$$\sum _{k=i}^n{z}_k=\sum _{k=i}^n{q}_i$$

$$z_{i_s-1}=q_{i_s-1}-\left(\sum _{k=i_s}^n{p}_k-\sum _{k=i_s}^n{q}_k\right)\ge p_{i_s-1}$$

$$z_{i_s-1}=p_{i_s-1}-\left(\sum _{k=i_s}^n{q}_k-\sum _{k=i_s}^n{p}_k\right)\ge q_{i_s-1}$$

$$\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{\}}.$$

$$z^{(d)}+\sum _{y\in Q}y=x.$$

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

$$H(\mathbf{M})$$

$$\sum _{k\ge i_s}{m}_{k{j}^{\prime }}=q_{j^{\prime }},\text{mbox}andR[j^{\prime }]+\sum _{k\ge i_s}{m}_{j^{\prime }k}=p_{j^{\prime }}.$$

$$\sum _{k\ge i_{s-1}}{m}_{i_{s-1}-1k}=\sum _{k=i_{s-1}}^n{q}_k-\sum _{k=i_{s-1}}^n{p}_k=p_{i_{s-1}-1}-z_{i_{s-1}-1}.$$

$$q_{i_{s-1}-1}\le z_{i_{s-1}-1}=p_{i_{s-1}-1}-\sum _{k\ge i_{s-1}}{m}_{i_{s-1}-1k}.$$

$$q_j=z_j^{(d)}+\sum _{\ell \in Q}R[\ell ]=m_{jj}+\sum _{k\ge j+1}{m}_{kj}=\sum _{k\ge i_s}{m}_{kj}$$

$$p_{i_{s-1}-1}$$

$$\sum _{k=i_s}^n{m}_{k{j}^{\prime }}=q_{j^{\prime }}\sum _{k=i_s}^n{m}_{j^{\prime }k}=p_{j^{\prime }}-R[j^{\prime }]$$

Full text

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 $X_{1},\ldots,X_{k}$. In this case, our approximation algorithm has an additive approximation factor of $\log k$.

We also discuss some related applications of our findings.

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 [13] and Frechet [9], 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 [1, 6, 8, 16] 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 general 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}$. This problem arises in many situations. For instance, the authors of [14] consider the important question of identifying the correct causal direction between two arbitrary r.v. $X$ and $Y$, that is, they want to discover whether it is the case that $X$ causes $Y$ or it is $Y$ that causes $Y$. In general, $X$ causes $Y$ if there exists an exogenous r.v. $E$ 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 [14] 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 interesting 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., fixed) their marginal distributions. This is exactly the problem we consider in this paper. The authors of [14] also observe that the latter optimization problem is NP-hard (due to results of [15, 20]), and evaluate experimentally a greedy approximation algorithm to find the minimum-entropy joint distribution, given the marginals. No proved performance guarantee is given in [14] for that algorithm. In this paper, we give 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 of 1. Subsequently, in Section IV 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 approximation algorithm has an additive approximation factor of $\log k$.

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 [20]. 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})$, $m\leq n$, the distance ${\tt D}({\bf p},{\bf q})$ among ${\bf p}$ and ${\bf q}$ is defined as the quantity ${\tt D}({\bf p},{\bf q})=2W({\bf p},{\bf q})-H({\bf p})-H({\bf q})$, 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 [20] to the problem of order-reduction of stochastic processes. The author of [20] observes that the problem of computing $W({\bf p},{\bf q})$ is NP-hard and proposes another different greedy algorithm for its computation, based on some analogy with the problem of Bin Packing with overstuffing. Again, no performance guarantee is given in [20] for the proposed algorithm. Our result directly implies that we can compute the value of ${\tt D}({\bf p},{\bf q})$, for arbitrary ${\bf p}$ and ${\bf q}$, with an additive error of at most $1$.111We remark that in [5] we considered the different problem of computing the probability distributions ${\bf q}^{*}$ that minimizes ${\tt D}({\bf p},{\bf q})$, given ${\bf p}$.

There are many other problems that require the computation of the minimum-entropy joint probability distribution of two random variables, whose marginals are equal to ${\bf p}$ and ${\bf q}$. We shall limit ourselves to discuss a few additional examples, postponing a more complete examination in a future version of the paper. To this purpose, let us write the joint entropy of two r.v. $X$ and $Y$, distributed according to ${\bf p}$ and ${\bf q}$, respectively, as $H(XY)=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 the paper [15] this maximal mutual information is interpreted, in agreement with Renyi axioms for a bona fide dependence measure [19], as a measure of the largest possible dependence of the two r.v. $X$ and $Y$. Since the problem of its exact computation is obviously NP-hard, our result implies an approximation algorithm for it. 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 [18, 21]. Finally, our problem could also be seen as a kind of “channel-synthesis” problem, where it is given pair of r.v. $(X,Y)$, and the goal is to construct a memoryless channel that maximizes the mutual information $I(X;Y)$ between $X$ and $Y$.

II Mathematical Preliminaries

We start by recalling a few notions of majorization theory [17] 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 probabilities 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 set

[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 that the partially ordered set $({\cal P}_{n},\preceq)$ is indeed a lattice [2],222The same result was independently rediscovered in [7], see also [10] 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]

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

[TABLE]

In the paper [2] the authors also gave a simple and efficient algorithm to 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**.**

[2]* 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]

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

*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 [12]. There, the authors prove that ${\bf x}\preceq{\bf y}$ implies $H({\bf x})\geq H({\bf y})+D({\bf y}||{\bf x})$, where $D({\bf y}||{\bf x})$ is the relative entropy between ${\bf x}$ and ${\bf y}$.

We also need the concept of aggregation [20, 5] and a result from [5], whose proof is repetead here to make the paper self-contained. Given ${\bf p}=(p_{1},\ldots,p_{n})\in{\cal P}_{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**.**

[5]* 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 consequences of 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}$). For our purposes, each element in ${\cal C}({\bf p},{\bf q})$ can be seen as $n\times m$ matrix $M=[m_{ij}]\in\mathbb{R}^{n\times m}$ such that its row-sums give the elements of ${\bf p}$ and its column-sums give the elements of ${\bf q}$, 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 (1) 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 that333Recall that we use the majorization relationship between vectors of unequal lenghts, 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\frac{1}{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 (3), 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 important result (for us).

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 is one of the key results towards our approximation algorithm to find an element $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ with entropy $H(\mathbf{M})\leq OPT+1$, where $OPT=\min_{\mathbf{N}\in{\cal C}({\bf p},{\bf q})}H(\mathbf{N}).$

Before describing our algorithm, let us illustrate some interesting consequences of Lemma 2. It is well known that for any joint distribution of the two r.v. $X$ and $Y$ it holds that

[TABLE]

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

[TABLE]

Lemma 2 strengthens the lower bound (5). Indeed, since, by definition, it holds that ${\bf p}\land{\bf q}\preceq{\bf p}$ and ${\bf p}\land{\bf q}\preceq{\bf q}$, by the Schur-concavity of the entropy function and Lemma 2 we get the (improved) lower bound

[TABLE]

Inequality (6) 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 (6) implies

[TABLE]

The new bounds 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}$ [12]. However, in this paper we just need what we can get from the inequality $H({\bf x})\geq H({\bf y})$, whenever ${\bf x}\preceq{\bf y}$ holds.

Inequalities (6) and (7) could be useful also in other contexts, when one needs to bound the joint entropy (or the mutual information) of two r.v. $X$ and $Y$, and the only available knowledge is given by the marginal distributions of $X$ and $Y$ (and not their joint distribution).

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

In this section we present an algorithm that, having in input distributions ${\bf p}$ and ${\bf q}$, constructs an $\mathbf{M}\in{\cal C}({\bf p},{\bf q})$ such that

[TABLE]

Lemma 2 will imply that

[TABLE]

We need to introduce some additional notations and state some properties which will be used in the description of our algorithm.

Definition 2**.**

Let ${\bf p}=(p_{1},\dots,p_{n})$ and ${\bf q}=(q_{1},\dots,q_{n})$ be two probability distributions in ${\cal P}_{n}$. We assume that for the maximum $i\in\{1,\dots,n\}$ such that $p_{i}\neq q_{i}$—if it exists—it holds that $p_{i}>q_{i}.$444Notice that up to swapping the role of ${\bf p}$ and ${\bf q}$, the definition applies to any pair of distinct distributions. Let $k$ be the minimum integer such that there are indices $i_{0}=n+1>i_{1}>i_{2}>\cdots>i_{k}=1$ satisfying the following conditions for each $s=1,\dots,k$:

• •

if $s$ is odd, then $i_{s}$ is the minimum index smaller than $i_{s-1}$ such that $\sum_{k=i}^{n}p_{k}\geq\sum_{k=i}^{n}q_{k}$ holds for each $i=i_{s},i_{s}+1,\dots,i_{s-1}-1$;

• •

if $s$ is even, then $i_{s}$ is the minimum index smaller than $i_{s-1}$ such that $\sum_{k=i}^{n}p_{k}{\leq}\sum_{k=i}^{n}q_{k}$ holds for each $i=i_{s},i_{s}+1,\dots,i_{s-1}-1$.

We refer to the integers $i_{0},i_{1},\dots,i_{k}$ as the inversion points of ${\bf p}$ and ${\bf q}$.555If ${\bf p}={\bf q}$ then we have $k=1$ and $i_{1}=1.$

Fact 2**.**

Let ${\bf p}$ and ${\bf q}$ be two probability distributions in ${\cal P}_{n}$ and $i_{0}=n+1>i_{1}\cdots>i_{k}=1$ be their inversion points. Let ${\bf z}={\bf p}\wedge{\bf q}.$ Then the following relationships hold:

for each odd $s\in\{1,\dots,k\}$ and $i\in\{i_{s},\dots,i_{s-1}-1\}$

[TABLE] 2. 2.

for each even $s\in\{1,\dots,k\}$ and $i\in\{i_{s},\dots,i_{s-1}-1\}$

[TABLE] 3. 3.

for each odd $s\in\{1,\dots,k\}$ and $i\in\{i_{s},\dots,i_{s-1}-2\}$ we have $z_{i}=p_{i}$ 4. 4.

for each even $s\in\{1,\dots,k\}$ and $i\in\{i_{s},\dots,i_{s-1}-2\}$ we have $z_{i}=q_{i}$ 5. 5.

for each odd $s\in\{0,\dots,k-1\}$

[TABLE] 6. 6.

for each even $s\in\{0,\dots,k-1\}$

[TABLE]

Proof:

By Fact 1 it holds that

[TABLE]

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

[TABLE]

This, together with the definition of the inversion indices $i_{0},\dots,i_{k}$, imply properties 1) and 2). The remaining properties are easily derived from 1) and 2) by simple algebraic calculations. ∎

Lemma 3**.**

Let $A$ be a multiset of non-negative real numbers and $z$ a positive real number such that $z\geq y$ for each $y\in A.$ For any $x\geq 0$ such that $x\leq z+\sum_{y\in A}y$ there exists a subset $Q\subseteq A$ and $0\leq z^{(d)}\leq z$ such that

[TABLE]

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

Proof:

If $\sum_{y\in A}y<x,$ we get $Q=A$ and the desired result directly follows from the assumption that $z+\sum_{y\in A}y\geq x.$ Note that the condition can be checked in linear time.

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

This lemma is a major technical tool of our main algorithm. We present a procedure implementing the construction of the the split of $z$ and the set $Q$ in Algorithm 2.

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 presented is given below (Algorithm 1). An informal description of it, that gives also the intuition behind its functioning, is presented in subsection III-B.

The following theorem shows the correctness of Algorithm 1. It relies on a sequence of technical results, Lemmas 4 and 5 and Corollaries 1 and 2, whose statements are deferred to the end of this section.

Theorem 1**.**

For any pair of probability distributions ${\bf p},{\bf q}\in{\cal P}_{n}$ the output of Algorithm 1 is an $n\times n$ matrix $\mathbf{M}=[m_{i\,j}]\in{\cal C}({\bf p},{\bf q})$ i.e., such that $\sum_{j}m_{i\,j}=p_{i}$ and $\sum_{i}m_{i\,j}=q_{j}.$

Proof:

Let $k$ be the value of $s$ in the last iteration of the line 6. for-loop of the algorithm. The desired result directly follows by Corollaries 1 and 2 according to whether $k$ is odd or even, respectively. ∎

We now 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 non-zero entry of the matrix $\mathbf{M}$ are all fixed in lines 12 and 21—in fact, for the assignments in lines 15 and 24 the algorithm uses values stored in $R$ or $C$ which were fixed at some point earlier in lines 12 and 21. 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 additive approximation guarantee of 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 lines 3-4 takes $O(n^{2})$;

• •

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

• •

the computation of the inversion points of ${\bf p}$ and ${\bf q}$ in line 4 can be performed in $O(n)$ following Definition 2, 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 vector $z={\bf p}\land{\bf q}$ can be computed in $O(n)$, e.g., based on the precomputed suffix sums;

• •

in the main body of the algorithm, the most expensive parts are the calls to the procedure Lemma3, and the for-loops in lines 10, 19, 14, and 23. All these take $O(n)$ and it is not hard to see that they are executed at most $O(n)$ times (once per component of ${\bf z}$). Therefore, the main body of the algorithm in lines 5-24 takes $O(n^{2})$.

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

III-A The analysis of correctness of Algorithm 1: technical lemmas

In this section we state four technical lemmas we used for the analysis of Algorithm 1 which leads to Theorem 1. In the Appendix, we give a numerical example of an execution of Algorithm 1.

Lemma 4**.**

At the end of each iteration of the for-loop of lines 7-12 in Algorithm 1 ($s=1,\dots,k,$ and $i_{s}\leq j<i_{s-1}$) we have (i) $m_{\ell\,c}=0$ for each $\ell,c$ such that $\min\{\ell,c\}<j$; and (ii) for each $j^{\prime}=j,\dots,i_{s-1}-1$

[TABLE]

Proof:

For (i) we observe that, before the first iteration ($j=i_{s-1}-1$) the condition holds (by line 3, when $s=1$, and by Corollary 2 for odd $s>1$). Then, within each iteration of the for-loop values $m_{\ell\,c}$ only change in lines 12, where (i) is clearly preserved, and in line 11 where, as the result of call to Algorithm 2, we have $\ell\in Q\subseteq\{j+1,\dots,i_{s-1}-1\},$ which again preserves (i).

We now prove (ii) by induction on the value of $j$. First we observe that at the beginning of the first iteration of the for-loop (lines 7-12), i.e., for $j=i_{s-1}-1,$ it holds that

[TABLE]

This is true when $s=1,$ since in this case we have $i_{s-1}=n+1,$ hence the two sums in the middle term are both 0; the first term is 0 since no term in $M$ has been fixed yet, and the last term is also 0, since $p_{n}=z_{n}$ by assumption. The equation is also true for each odd $s>1$ by Corollary 2. Moreover, at the beginning of the first iteration ($j=i_{s-1}-1$) it holds that $R[\ell]=0$ for each $\ell=1,\dots,n.$ This is true for $s=1$ because of the initial setting in line 2. For $s>1$ the property holds since any $R[\ell]$ is only assigned non-zero value within the for-loop (lines 7-12) and unless the algorithm stops any non-zero $R[\ell]$ is zeroed again immediately after the for-loop, in lines 14-15 unless the exit condition $i_{s}=1$ is verified which means that $s=k$ and the algorithm terminates immediately after.

When Algorithm 1 enters the for-loop at lines 7-12, $s$ is odd. Then, by point 6. of Fact 2 and (15) it holds that

[TABLE]

This implies that for $j=i_{s-1}-1$ the values $z_{j},q_{j}$ together with the values in $R[j+1\dots i_{s-1}-1]$ satisfy the hypotheses of Lemma 3. Hence, the call in line 7 to Algorithm 2 (implementing the construction in the proof of Lemma 3) correctly returns a splitting of $z_{j}$ into two parts $z^{(d)}_{j}$ and $z_{j}^{(r)}$ and a set of indices $Q\subseteq\{j+1,\dots,i_{s-1}-1\}$ s.t.

[TABLE]

where the first equality holds after the execution of lines 10-12, and the second equality holds because by (i) $m_{k\,j}=0$ for $k<j$. We have established the first equation of (14). Moreover, the second equation of (14) also holds because by the equality in (16), the result of the assignment in line 12 and (by (i), with $j=i_{s-1}-1$) $m_{i_{s-1}\,k}=0$ for $i_{s}\leq k<i_{s-1}-1$, we get

[TABLE]

We now argue for the cases $j=i_{s-1}-2,i_{s-1}-3,\dots,i_{s}.$ By induction we can assume that at the beginning of any iteration of the for-loop (lines 7-12) with $i_{s}\leq j<i_{s-1}-1,$ we have that for each $i_{s-1}-1\geq j^{\prime}>j$

[TABLE]

and (if $s>1$) by Corollary 2 for each $j^{\prime}\geq i_{s-1}$ we have

[TABLE]

Moreover, by point 1. and 3. of Fact 2 we have

[TABLE]

From these, we have

[TABLE]

where (20) follows by using (17) and (18); (21) follows from (20) by simple algebraic manipulations; finally (22) follows from (21) because, by (i), at the end of iteration $j+1$, we have $m_{\ell\,c}=0$ if $\ell<j+1$ or $c<j+1$; hence

[TABLE]

and

[TABLE]

and the equal terms and cancel out.

For each $k=j+1,\dots,i_{s-1}-1$ such that $R[k]\neq 0$ we have $R[k]=z^{(r)}_{k}\leq z_{k}\leq z_{j}$, where we are using the fact that for ${\bf z}={\bf p}\wedge{\bf q}$ it holds that $z_{1}\geq\cdots\geq z_{n}.$

Therefore $z_{j},q_{j}$ and the values in $R[j+1\dots,i_{s-1}-1]$ satisfy the hypotheses of Lemma 3. Hence, the call in line 7 to Algorithm 2 (implementing the construction of Lemma 3) correctly returns a splitting of $z_{j}$ into two parts $z^{(d)}_{j}$ and $z_{j}^{(r)}$ and a set of indices $Q\subseteq\{j+1,\dots,i_{s-1}-1\}$ s.t. $q_{j}=z^{(d)}_{j}+\sum_{\ell\in Q}R[\ell].$ Then we can use the same argument used in the first part of this proof (for the base case $j=i_{s-1}-1$) to show that the first equation of (14) holds after lines 10-12.

For $j^{\prime}=j$, the second equation in (14) is guaranteed by the assignment in line 12. Moreover, for $j^{\prime}>j$ and $j^{\prime}\not\in Q$ it holds since no entry $m_{j^{\prime}\,k}$ or $R[j^{\prime}]$ is modified in lines 10-12. Finally, for $j^{\prime}>j$ and $j^{\prime}\in Q$ before the execution of lines 10-12 we had $p_{j^{\prime}}=R[j^{\prime}]+\sum_{k\geq i_{s}}m_{j^{\prime}\,k}$ with $R[j]>0$ and $m_{j^{\prime}\,j}=0$ and after the execution of lines 10-12 the values of $R[j^{\prime}]$ and $m_{j^{\prime}\,j}$ are swapped, hence the equality still holds. The proof of the second equality of (14) is now complete. ∎

Corollary 1**.**

When algorithm 1 reaches line 16, it holds that

[TABLE]

[TABLE]

and (if $i_{s}\neq 1$, i.e., this is not the last iteration of the outermost for-loop)

[TABLE]

Proof:

We will only prove (23) and (24) for $j=i_{s},\dots,i_{s-1}-1.$ In fact, this is all we have to show if $s=1$. Furthermore, if $s>1,$ then in the previous iteration of the outermost for-loop the algorithm has reached line 24 and by Corollary 2 (23) and (24) also hold for each $j\geq i_{s-1},$ as desired.

By Lemma 4 when the algorithm reaches line 13 we have that for each $j=i_{s},\dots,i_{s-1}-1$

[TABLE]

Note that the entries of $M$ in columns $i_{s},\dots,i_{s-1}-1$ will not be changed again and within the for-loop in lines 7-12 only the values $m_{i\,j}$ with $j=i_{s},\dots i_{s-1}-1$ and $i=j,\dots,i_{s-1}-1$ may have been changed. Hence, for each $j=i_{s},\dots,i_{s-1}-1,$ it holds that $\sum_{k<i_{s}}m_{k\,j}=0$ and $\sum_{k<i_{s}-1}m_{j\,k}=0$ as desired.

Moreover, by Lemma 4 (i) with $j=i_{s}$ it holds that $m_{\ell\,c}=0,$ when $\ell<i_{s}$ or $c<i_{s}.$ Hence, for each $j=i_{s},\dots,i_{s-1}-1,$ it holds that $\sum_{k<i_{s}}m_{k\,j}=0$ and $\sum_{k<i_{s}-1}m_{j\,k}=0$ as desired.

Since the operations in lines 13-16 only change values in column $i_{s}-1$ of $M$ and in the vector $R$, the first equation in (26) directly implies that (23) holds for each $j=i_{s},\dots,i_{s-1}.$

With the aim of proving (24) let us first observe that if $i_{s}=1$ (hence $s=k$ and the algorithm is performing the last iteration of the outermost for-loop) then by (23) and (26) we have

[TABLE]

and, since $R[j]\geq 0,$ it follows that $R[k]=0,$ for each $j=i_{s},\dots,i_{s-1}-1.$

Now, first assume that $i_{s}>1$ hence $s<k$ From (26) for each $j=i_{s},\dots,i_{s-1}-1$ such that $R[j]=0$ we immediately have that (24) is also satisfied. Hence, this is the case for all $j=i_{s},\dots,i_{s-1}-1$ when $s=k$ and $i_{s}=1.$ Moreover, if there is some $j\in\{i_{s},\dots,i_{s-1}-1\}$ such that $R[j]\neq 0$ (when $s<k$ and $i_{s}>1,$), after the execution of line 15, for each $j=i_{s},\dots,i_{s-1}-1$ such that $R[j]$ was $\neq 0$ we have $m_{j\,i_{s}-1}=R[j]$, hence,

[TABLE]

completing the proof of (24).

Finally, we prove (25). By the assignments in line 15 and the fact that this is the first time that values in column $i_{s}-1$ of $M$ are set to non-zero values, from point 5 of Fact 2 we get

[TABLE]

that, together with the fact that at this point the only non-zero values in column $i_{s}-1$ of $M$ are in the rows $i_{s},\dots,i_{s-1}-1,$ completes the proof of (25). ∎

Lemma 5**.**

At the end of each iteration of the for-loop of lines lines 17-21 in Algorithm 1 ($s=1,\dots,k,$ and $i_{s}\leq j<i_{s-1}$) we have (i) $m_{\ell\,c}=0$ for each $\ell,c$ such that $\min\{\ell,c\}<j$; and (ii) for each $j^{\prime}=j,\dots,i_{s-1}-1$,

[TABLE]

Proof:

The proof can be easily obtained by proceeding like in Lemma 5 (swapping the roles of rows and columns of $\mathbf{M}$ and ${\bf p}$ and ${\bf q}$). ∎

Corollary 2**.**

When algorithm 1 reaches line 24, it holds that

[TABLE]

[TABLE]

and if $i_{s}\neq 1$ (the outermost for-loop is not in last iteration)

[TABLE]

Proof:

The proof can be easily obtained by proceeding like in Corollary 1 (swapping the roles of rows and columns of $\mathbf{M}$ and of ${\bf p}$ and ${\bf q}$). ∎

III-B How Algorithm 1 works: An informal description of its functioning

Given the inversion points $i_{0},i_{1},\dots,i_{k}$ of the two probability distributions ${\bf p}$ and ${\bf q}$, as defined in Definition 2, for each $s=1,\dots,k$ let us call the list of integers $L^{s}=\{i_{s-1}-1,i_{s-1}-2,\dots,i_{s}\}$ (listed in decreasing order) a ${\bf p}$-segment, or a ${\bf q}$-segment, according to whether $s$ is odd or even. For each $i$ belonging to a ${\bf p}$-segment we have

[TABLE]

For each $i$ belonging to a ${\bf q}$-segment we have

[TABLE]

Algorithm 1 proceeds by filling entries of the matrix $M$ with non-zero values. Other possible actions of the algorithm consist in moving probabilities from one entry of $M$ to a neighboring one. The reasons of this moving will become clear as the description of the algorithm unfolds.

At any point during the execution of the algorithm, we say that a column $i$ is satisfied if the sum of the entries on column $i$ is equal to $q_{i}$. Analogously, we say that a row $i$ is satisfied if the sum of the entries on row $i$ is equal to $p_{i}$. Obviously, the goal of the algorithm is to satisfy all rows and columns. Line 3 makes sure that the first value of $j$ in line 7 is in a ${\bf p}$-segment.

Line 3 makes sure that the first value of $j$ in line 7 is in a ${\bf p}$-segment. For each $j=n,\dots,1,$ with $j$ in a ${\bf p}$-segment, the algorithm maintains the following invariants:

• 1-p

all rows $j^{\prime}>j$ are satisfied

• 2-p

all columns $j^{\prime}>j$ are satisfied

• 3-p

the non-zero entries $M[j^{\prime},j]$ for $j^{\prime}>j$ in the same ${\bf p}$-segment, satisfy $M[j^{\prime},j]+M[j^{\prime},j^{\prime}]=z_{j^{\prime}}$

The main steps of the algorithm when $j$ is in a ${\bf p}$-segment amount to:

• • Step1-p:

Put $z_{j}$ in $M[j,j]$. By the assumption that $j$ is part of a ${\bf p}$-segment, we have that this assignment satisfies also row $j$. However, this assignment might create an excess on column $j$, (i.e., the sum of the elements on column $j$ could be greater than the value of $q_{j}$) since by the invariants 1-p and 2-p and the assigned value to the entry $M[j,j]$ we have that the sum of all the entries filled so far equals to

[TABLE]

and the entries on the columns $j^{\prime}>j$ satisfy exactly $q_{j^{\prime}}$, that is, sums up to $q_{j^{\prime}}$.

• Step2-p:

If there is an excess on column $j$, we adjust it by applying Lemma 3. Indeed, by Lemma 3 we can select entries $M[j^{\prime},j]$ that together with part of $M[j,j]=z_{j}$ sum up exactly to $q_{j}.$ The remaining part of $z_{j}=M[j,j]$ and each of the non-selected entries on column $j$ are kept on their same row but moved to column $j-1$. In the pseudocode of Algorithm 1, this operation is simulated by using the auxiliary array $R[\cdot]$. Notice that by this operation we are maintaining invariants $1$-p, $2$-p, $3$-p for $j\leftarrow j+1.$

Step1-p and Step2-p are repeated as long as $j$ is part of a ${\bf p}$-segment.

When $j$ becomes part of a ${\bf q}$-segment, the roles pf ${\bf p}$ and ${\bf q}$ are inverted, namely, we have that the following invariants hold:

• 1-q

all rows $j^{\prime}>j$ are satisfied

• 2-q

all columns $j^{\prime}>j$ are satisfied

• 3-q

the non-zero entries $M[j,j^{\prime}]$ for $j^{\prime}>j$ in the same ${\bf q}$-segment, satisfy $M[j,j^{\prime}]+M[j^{\prime},j^{\prime}]=z_{j^{\prime}}$

From now on and as long as $j$ is part of a ${\bf q}$-segment the main steps of the algorithm amount to:

• • Step1-q:

Put $z_{j}$ in $M[j,j]$. By the assumption that $j$ is part of a ${\bf q}$-segment, we have that this assignment satisfies also column $j$. Again, this assignment might create an excess on row $j$, since by the invariants 1-q and 2-q and the entry $M[j,j]$ we have that the sum of all the entries filled so far is equal to

[TABLE]

and the entries on the rows $j^{\prime}>j$ satisfy exactly $p_{j^{\prime}}.$

• Step2-q:

If there is an excess on row $j$, by Lemma 3 we can select entries $M[j,j^{\prime}]$ that together with part of $M[j,j]=z_{j}$ sum up exactly to $p_{j}.$ The remaining part of $z_{j}=M[j,j]$ and the non-selected entries are kept on the same column but are moved up to row $j-1$. In the pseudocode of Algorithm 1, this operation is simulated by using the auxiliary array $C[\cdot]$. Notice that by this operation we are maintaining invariants $1$-q, $2$-q, and $3$-q for $j\leftarrow j+1.$

Again these steps are repeated as long as $j$ is part of a ${\bf q}$-segment. When $j$ becomes part of a ${\bf p}$-segment again, we will have that once more invariants $1$-p, $2$-p, and $3$-p are satisfied. Then, the algorithm resorts to repeat steps Step1-p and Step2-p as long as $j$ is part of a ${\bf p}$-segment, and so on and so forth switching between ${\bf p}$-segments and ${\bf q}$-segments, until all the rows and columns are satisfied.

From the above description it is easy to see that all the values used to fill in entries of the matrix are created by splitting into two parts some element $z_{j}$. This key property of the algorithm implies the bound on the entropy of $M$ being at most $H({\bf z})+1.$

We note here that (for efficiency reasons) in the pseudocode of Algorithm 1 instead of moving values from one column to the next one (Step2-p) or from one row to the next one (Step2-q), the arrays $R$ and $C$ are used, where $R[j^{\prime}]$ plays the role of $M[j^{\prime},j]$, in invariant $3$-p above, and $C[j^{\prime}]$ plays the role of $M[j,j^{\prime}]$, in invariant $3$-q above.

IV Extending the result to more distributions

In this section we will show how the algorithm Min-Entropy-Joint-Distribution can be used to attain a $\log k$ additive approximation for the problem of minimising the entropy of the joint distribution of $k$ 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-Distribution 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.

We will 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)}.$ Algorithm 3 shows the pseudocode for our procedure.

The vector $Indices^{(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]$. Therefore, if after the execution of line 17, for $w=1,\dots,|{\bf m}^{(i-j)}|$, we have $Indices^{(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,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)}).$

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. In the rest of this section we will prove that the entropy of the joint distribution output by the algorithm guarantees additive $\log k$ approximation.

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 3**.**

For any ${\bf p}=(p_{1},\dots,p_{n})\in{\cal P}_{n}$ we denote by ${\rm 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 ${\rm half}^{(i)}({\bf p})={\rm half}({\rm half}^{(i-1)}({\bf p})),$ where ${\rm half}^{(1)}({\bf p})={\rm half}({\bf p})$ and ${\rm half}^{(0)}({\bf p})={\bf p}.$

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

Lemma 6**.**

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

Lemma 7**.**

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

[TABLE]

Theorem 3**.**

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 ${\rm 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)}={\rm 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

• •

(35) follows from (34) by the definition of the operator ${\rm half}$;

• •

(36) follows from (35) by Lemma 7;

• •

(37) follows from (36) by the induction hypotheses (32)-(33) ;

• •

(38) follows from (37) 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 ${\rm 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 ${\rm half}({\bf z})\preceq{\bf m}^{(i-j)}$ yielding 38. ∎

An immediate consequence of the last theorem is the following

Corollary 3**.**

For any $k$ probability distributions ${\bf p}^{(1)},\dots,{\bf p}^{(k)}$ let $M$ be the joint distribution of ${\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 repeat some of the probability distribution until there are $k^{\prime}=2^{\lceil\log k\rceil}$ of them. By Theorem 3 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({\rm 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 $\log k$ approximation for the problem of computing the joint distribution of minimum entropy of $k$ input distributions.

Lemma 8**.**

Fix $k$ distributions ${\bf p}^{(1)},{\bf p}^{(2)},\cdots,{\bf p}^{(k)}$. For any $M$ being a joint distribution of ${\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. ∎

Summarising we have shown the following

Theorem 4**.**

Let ${\bf p}^{(1)},\dots,{\bf p}^{(m)}\in{\cal P}_{n}$. Let $M^{*}$ be a joint distribution of ${\bf p}^{(1)},\dots,{\bf p}^{(m)}$ of minimum entropy among all the joint distribution of ${\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) algorithm provides an additive $\log(m)$ approximation.

IV-A The proofs of the two technical lemmas about the operator ${\rm half}$

Lemma 6. For any ${\bf p}\preceq{\bf q}$ we have also ${\rm half}({\bf p})\preceq{\rm 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 ${\rm half}({\bf p})_{1}\geq{\rm half}({\bf p})_{2}\geq\dots\geq{\rm half}({\bf p})_{2n}$ and ${\rm half}({\bf q})_{1}\geq{\rm half}({\bf q})_{2}\geq\dots\geq{\rm 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 ${\rm half}({\bf p})\preceq{\rm half}({\bf q}).$ ∎

Lemma 9**.**

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

[TABLE]

Proof.

By the previous lemma 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 this Lemma we can extend the result to ”powers” of the operator ${\rm half}$ and have our Lemma 7.

Lemma 7. 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 lemma. 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 Lemma 9. ∎

-B A numerical example

Fix $n=13$ and let

[TABLE]

and

[TABLE]

be the two probability distribution for which we are seeking a joint probability of minimum entropy. By Fact 1 we have

[TABLE]

By Definition 2 we have that the inversion points are $i_{0}=14,i_{1}=11,i_{2}=9,i_{3}=6,i_{4}=1.$

The resulting joint distribution produced by Algorithm 1 is given by the following matrix $\mathbf{M}=[m_{i\,j}]$, satisfying the property that $\sum_{j}m_{i\,j}=p_{i}$ and $\sum_{i}m_{i\,j}=q_{j}.$

[TABLE]

Notice that by construction

• •

for the submatrix $\mathbf{M}^{(i_{1})}=[m_{i\,j}]_{i_{1}\leq i\leq i_{0}-1,\,i_{1}-1\leq j\leq i_{0}-1}$ we have that each row $i$ contains at most two elements and the sum of the elements on the row equals $z_{i}$

• •

for the submatrix $\mathbf{M}^{(i_{2})}=[m_{i\,j}]_{i_{2}-1\leq i\leq i_{1}-1,\,i_{2}\leq j\leq i_{1}-1}$ we have that each column $i$ contains at most two elements and the sum of the elements on the column equals $z_{i}$

• •

for the submatrix $\mathbf{M}^{(i_{3})}=[m_{i\,j}]_{i_{3}\leq i\leq i_{2}-1,\,i_{3}-1\leq j\leq i_{2}-1}$ we have that each row $i$ contains at most two elements and the sum of the elements on the row equals $z_{i}$

• •

for the submatrix $\mathbf{M}^{(i_{4})}=[m_{i\,j}]_{i_{4}\leq i\leq i_{3}-1,\,i_{4}\leq j\leq i_{1}-1}$ we have that each column $i$ c contains at most two elements and the sum of the elements on the column $z_{i}$

Notice that these four sub-matrices cover all the non-zero entries of $\mathbf{M}$. This easily shows that the entries of the matrix $\mathbf{M}$ are obtained by splitting into at most two pieces the components of ${\bf z}$, implying the desired bound $H(\mathbf{M})\leq H({\bf z})+1.$