Title: AN INFORMATION THEORETIC APPROACH TO DIS-TRIBUTED REPRESENTATION LEARNING

Abstract: The problem of distributed representation learning is one in which multiple sources of information X1, . . . , XK are processed separately so as to extract useful information about some statistically correlated ground truth Y . We investigate this problem from informationtheoretic grounds. For both discrete memoryless (DM) and memoryless vector Gaussian models, we establish fundamental limits of learning in terms of optimal tradeoffs between relevance and complexity. We also develop a variational bound on the optimal tradeoff that generalizes the evidence lower bound (ELBO) to the distributed setting. Furthermore, we provide a variational inference type algorithm that allows to compute this bound and in which the mappings are parametrized by neural networks and the bound approximated by Markov sampling and optimized with stochastic gradient descent. Experimental results on synthetic and real datasets are provided to support the efficiency of the approaches and algorithms which we develop in this paper.

1. INTRODUCTION:

Let a measurable variable X ∈ X and a target variable Y ∈ Y with unknown joint distribution PX,Y be given. In the classic problem of statistical learning, one wishes to infer an accurate predictor of the target variable Y ∈ Y based on observed realizations of X ∈ X . That is, for a given class F of admissible predictors φ : X → Ŷ and an additive loss function : Y → Ŷ that measures discrepancies between true values and their estimated fits, one aims at finding the mapping φ ∈ F that minimizes the expected risk CP X,Y (φ, ) = EP X,Y [ (Y, φ(X))].

(1)

Because the joint distribution PX,Y is unknown, in practice the risk equation 1 (also called population risk) cannot be computed directly; and, in the standard approach, one usually resorts to choosing the predictor with minimal risk on a training dataset consisting of n labeled samples {(xi, yi)} n i=1 that are drawn independently from the unknown joint distribution PX,Y . Also, it is important to restrict the set F of admissible predictors to a low-complexity class to prevent overfitting. This leads to the abstract inference problem shown in Figure 1 .

In this paper, we study a generalization of this problem in which the prediction is to be performed in a distributed manner. The model is shown in Figure 2 . Here, the prediction of the target variable Y ∈ Y is to be performed on the basis of samples of statistically correlated random variables (X1, . . . , XK ) that are observed each at a distinct predictor. We investigate this problem in the case in which the loss function (•) is the logarithmic-loss fidelity measure, given by log(y, ŷ) = log 1 ŷ(y)

(2)

where ŷ(•) designates a probability distribution on Y and ŷ(y) is the value of this distribution evaluated for the outcome y ∈ Y. The choice of a 'good" loss function is often controversial in statistical learning theory, and although a complete and rigorous justification of the usage of logarithmic loss as a fidelity measure in learning theory is still awaited, partial explanations appeared in Jiao et al. (2015) and, especially in Painsky and Wornell (2018) where it is shown that, for binary classification problems, by minimizing the logarithmic-loss one actually minimizes an upper bound to any choice of loss function that is smooth, proper (i.e., unbiased and Fisher consistent) and convex. Also, we constrain the complexity of the predictors by using mutual information as a regularizer term. This is inline with recent works Xu and Raginsky (2017) ; Russo and Zou (2015) that show that the generalization error can be upper-bounded using the mutual information between the input dataset and the output of the predictor -see also Bousquet and Elisseeff (2002) ; Shalev-Shwartz et al. (2010) where the stability of an algorithm is controlled by constraining the mutual information between its input and output.

P X|Y Y ∈ Y φ ψ Y ∈ Y X ∈ X U = φ(X)

Figure 1 : An abstract inference model for learning.

P X1,...,XK |Y Y ∈ Y • • • φ 1 φ k ψ Y ∈ Y XK X1 U1

UK:

Figure 2 : A model for distributed, e.g., multi-view, learning.

1.1. AN EXAMPLE: MULTI-VIEW LEARNING:

In many data analytics problems, data is collected from various sources of information or feature extractors; and is intrinsically heterogeneous. For example, an image can be identified by its color or texture features; and a document may contain text and images. Conventional machine learning approaches concatenate all available data into one big row vector (or matrix) on which a suitable algorithm is then applied. Treating different observations as a single source might cause overfitting and is not physically meaningful because each group of data may have different statistical properties. Alternatively, one may partition the data into groups according to samples homogeneity, and each group of data be regarded as a separate view. This paradigm, termed multi-view learning Xu et al. (2013) , has received growing interest; and various algorithms exist, sometimes under references such as co-training Blum and Mitchell (1998) ; Dhillon et al. (2011) ; Kumar and Daumé (2011) ; Gönen and Alpaydın (2011) , multiple kernel learning Gönen and Alpaydın (2011) and subspace learning Jia et al. (2010) . By using distinct encoder mappings to represent distinct groups of data, and jointly optimizing over all mappings to remove redundancy, multiview learning offers a degree of flexibility that is not only desirable in practice but is likely to result in better learning capability. Actually, as shown in Vapnik (2013) , local learning algorithms produce less errors than global ones. Viewing the problem as that of function approximation, the intuition is that it is usually non-easy to find a unique function that holds good predictability properties in the entire data space.

1.2. INFORMAL SUMMARY OF RESULTS:

In this paper, first we characterize the optimal tradeoff between relevance and complexity for the distributed learning model of Figure 2 for both discrete memoryless (DM) and memoryless vector Gaussian models. While the result for the discrete data model (Theorem 1) is not difficult to establish using connections with Courtade and Weissman (2014, Appendix B) which we explicit here, the result for the multivariate Gaussian data model (Theorem 2), which provides a sharp analytic characterization of optimal tradeoffs, is new and non-trivial (the proof of the converse part is not straightforward and was missing before this work in both theoretic learning and information theoretic communities including in the scalar case). Second, we develop a variational bound on the optimal tradeoff that can be seen as a generalization of the ELBO and the β-VAE criteria Higgins et al. (2016) to the distributed setting. Furthermore, for both DM and Gaussian models, we also provide a variational inference type algorithm which is parametrized by neural networks and allows to compute the developed variational bound when the data distribution is not known. Specifically, the main contributions of this paper are:

• In Section 3.2, we find an explicit analytic characterization of optimal tradeoffs between relevance and complexity for the memoryless vector Gaussian model. The result generalizes the Gaussian Information Bottleneck method of Globerson and Tishby (2004) ; Chechik et al. (Feb. 2005) to the distributed learning scenario.

• In Section 3.3, we study the problem of maximizing relevance under a constraint on the sum complexity for which we establish a variational bound which generalizes the ELBO and the β-VAE criteria to the distributed setting.

• Section 3.4 is algorithmic-oriented. We develop a variational inference type algorithm which enables to compute the bound. This algorithm is obtained by parametrizing the encoders, the decoder, and the prior distributions via DNNs and using Monte-Carlo sampling. Also, it makes usage of Kingma et al.'s re-parametrization trick Kingma and Welling (2013) and can be seen as a generalization of the variational information bottleneck algorithm in Alemi et al. (2017) to the distributed setting.

• Section 4 contains some experimental results on real datasets which show the efficiency of the approaches and algorithms that we develop in this paper.

Most relevant to this paper is the single-encoder Information Bottleneck (IB) method of Tishby et al. (1999) which readily and elegantly captures the above mentioned viewpoint of seeking the right balance between data fit and generalization by using the mutual information both as a cost function and as a regularizer term. Thus, the results of this paper can be seen as a generalization of those of Tishby et al. (1999) for the DM model and Globerson and Tishby (2004) ; Chechik et al. (Feb. 2005) for the Gaussian model to the distributed learning setting.

Remark: Due to space constraints, the proofs of the results of this paper are deferred to the appendices section, which also contains additional experimental results. 

1.3. NOTATION:

: k ∈ K}. If K = ∅, XK = ∅. For k ∈ K we let X K/k = (X1, . . . , X k-1 , X k+1 , .

. . , XK ), and assume that X0 = XK+1 = ∅. Also, for zero-mean random vectors X and Y, the quantities Σx, Σx,y and Σ x|y denote respectively the covariance matrix of the vector X, the covariance matric of vector (X, Y) and the conditional covariance matrix of X, conditionally on Y.

Finally, for two probability measures PX and QX on the random variable X ∈ X , the relative entropy or Kullback-Leibler divergence is denoted as DKL(PX QX ). 

2. FORMAL PROBLEM FORMULATION:

X k --Y --X K/k (3)

holds for all k ∈ K. That is, the joint pmf factorizes as PX 1 ,...,X K ,Y (x1, . . . , xK , y) = PY (y)

EQUATION

The variable Y is a target variable; and we seek to characterize how accurate it can be predicted from a measurable random vector (X1, . . . , XK ) when the components of this vector are processed separately, each by a distinct encoder. More specifically, let {(X1,i, . . . , XK,i, Yi)} n i=1 be a collection of n independent copies of (X1, . . . , XK , Y ). Encoder k ∈ K only observes the sequence X n k ; and generates a description

EQUATION

where

M (n) k

is an arbitrary set of descriptions. The range of allowable description sets will be specified below. A decoder ψ(•) collects all descriptions JK = (J1, . . . , JK ) and returns an estimate Ŷ n of Y n as

EQUATION

)

The relevance of the estimation Ŷ n is measured in terms of the relevance, defined here as the information that the descriptions φ1(X n 1 ), . . . , φK (X n K ) collectively preserve about Y n , as measured by Shannon mutual information1 

EQUATION

where Ŷ n = ψ(φ1(X n 1 ), . . . , φK (X n K )) and the subscript PX K ,Y indicates that the mutual information is computed under the joint distribution PX K ,Y .

There are various ways to control the complexity of the encoding functions {φ k } K k=1 . In this paper, we do so by restricting their ranges. This is known as minimum description length complexity measure Hinton and van Camp (1993) . Specifically, the mapping φ k (•) at Encoder k ∈ K needs to satisfy

EQUATION

)

Definition 1 A tuple (∆, R1, . . . , RK ) is said to be achievable if there exists an integer n, a family of encoding mappings {φ k } K k=1 and a decoder mapping ψ such that

EQUATION

)

The relevance-complexity region IRDIB is given by the closure of all achievable tuples (∆, R1, . . . , RK ).

In some cases, for given RK = (R1, . . . , RK ), for the ease of the exposition we will be content with the relevance-complexity function ∆(RK, PX K ,Y ) defined as

EQUATION

where the maximization is subjected to equation 8.

3. MAIN RESULTS:



3.1. DISCRETE MEMORYLESS DATA MODEL:

The following theorem (the proof of which can be found in the appendices section) provides a computable characterization of the relevance-complexity region IRDIB. The result can be seen as a generalization of Tishby et al. Tishby et al. (1999) single encoder IB to the distributed learning model with K encoders.

Theorem 1 The relevance-complexity region IRDIB of the distributed learning problem with PX K ,Y for which the Markov chain equation 3 holds is given by the union of all tuples (∆, R1, . . . , RK ) ∈ R K+1 + that satisfy for all S ⊆ K,

EQUATION

for some set of pmfs

EQUATION

)

Remark 1 In Theorem 1, the random variable T stands for a convexification of the region, i.e., convex combination of achievable relevance-complexity tuples is itself achievable. For given T = t, the result of Theorem1 comprises the optimization over K conditional distributions {P U K |X k ,t }. For k ∈ K, the conditional distribution P U K |X k ,t represents a stochastic encoding of the feature X k into a latent variable U k . Intuitively, the latent variable U k should capture all relevant information about Y that is contained in X k and non redundant with those carried out by {Ui} i =k . The requirement of non-redundancy is mandated by the need to operate at the minimum possible complexity at which a desired relevance level is achievable (recall that minimum complexity, as expressed by algorithm's input-output mutual information, translates directly into a better generalization capability). Collectively, however, the set of all latent variables (U1, . . . , UK ) should be expressive enough to reproduce the target variable Y to within the desired relevance level.

Remark 2 Like for the single-encoder IB problem of Tishby et al. (1999) and an increasing number of works that followed, including Courtade and Weissman (2014, Section III-F), our approach here is asymptotic. In addition to that it leads to an exact characterization, the result also readily provides a lower bound on the performance in the non-asymptotic (e.g., one shot) setting. For the latter setting known approaches (e.g., the functional representation lemma of Li and El Gamal ( 2018)) would lead to only non-matching inner and outer bounds on the region of optimal tradeoff pairs, as this is the case even for the single encoder case Li et al. (2018) .

3.2. MEMORYLESS VECTOR GAUSSIAN DATA MODEL:

We now turn to a continuous-alphabet setting. Here, (X1, . . . , XK , Y) is a zero-mean Gaussian random vector such that

EQUATION

where H k ∈ C n k ×ny models the linear model connecting the target variable Y ∈ C ny to the observation at encoder k, and N k ∈ C n k , k = 1, . . . , K, is the noise vector at encoder k, assumed to be Gaussian with zero-mean and covariance matrix Σ k , and independent from all other noises and the target variable Y. We denote by Σy the covariance matrix of of the target vector Y ∈ C ny .

For this model, we find an explicit analytic characterization of optimal tradeoffs between relevance and complexity. The proof relies on deriving an outer bound on the region described by equation 12, and showing that it is achievable with Gaussian distribution, with no time-sharing. In doing so, we use techniques that rely on the de Bruijn identity and the properties of Fisher information and minimum mean square error (MMSE).

Theorem 2 The relevance-complexity region IR G DIB for the vector Gaussian model is given by the union of all tuples (∆, R1, . . . , RL) that satisfy for all S ⊆ K

∆ ≤ R k + log I -Σ 1/2 k Ω k Σ 1/2 k + log k∈S c Σ 1/2 y H † k Ω k H k Σ 1/2 y + I , for some 0 Ω k Σ -1 k .

Proof: The proof of the direct part follows by evaluating the region of Theorem 1, which can be extended to the case of continuous alphabets using standard discretization (quantization) arguments, with the choices

T = ∅ and p(u k |x k , t) = CN (x k , Σ 1/2 k (Ω k -I)Σ 1/2 k ).

The main contribution in the proof is that of the converse part. This proof is technical and rather lengthy and, for this reason, is deferred to the appendices section.

In the special case in which K = 1, the result of Theorem 2 recovers that by Globerson and Tishby (2004) (see also Chechik et al. (Feb. 2005 )) which establishes the optimal relevance-complexity tradeoff of the single-encoder Gaussian IB problem.

3.3. A VARIATIONAL BOUND:

In this section, we consider the problem of learning encoders-and decoder mappings that maximize the relevance level for a given (fixed) complexity level, i.e., those that perform at the vicinity of the boundary of the region IRDIB. First, we derive a parametrization of the relevance-complexity region; and, then, we develop a variational bound which expresses the optimal encoders' and decoder mappings as the solution to an optimization problem -(an algorithm for solving this problem in the case of unknown distributions is given in the next section).

Let Rsum := K k=1 R k . Also, let IR sum DIB denote the region of achievable (relevance, sum-complexity) pairs,

IR sum DIB := (∆, Rsum) ∈ R 2 + : ∃(R1, . . . , RK ) ∈ R K + s.t.

(∆, R1, . . . , RK ) ∈ IRDIB and

K k=1 R k = Rsum .

Proposition 1 The relevance-complexity region under sum-complexity constraint RI sum DIB is given by the convexhull of all tuples (∆,

EQUATION

and where the maximization is over the set of pmfs

P := {P U 1 |X 1 , . . . , P U K |X K } such that the joint pmf factorizes as pY (y) K k=1 p X k |Y (x k |y) K k=1 p U k |X k (u k |x k ).

The next proposition provides a characterization of the pairs (∆, Rsum) that lie on the boundary of RI sum DIB in terms of a nonnegative parameter s ≥ 0.

Proposition 2 For every pair (∆, Rsum) ∈ R 2 + that lies on the boundary of the relevance-complexity region RI sum DIB there exist s ≥ 0 such that (∆, Rsum) = (∆s, Rs), where

EQUATION

EQUATION

and P * is the set of conditional pmfs P that maximize the cost function

EQUATION

Using Proposition 2 it is clear that the encoders {P U k |X k } k∈K that achieve the relevance-complexity pair (∆s, Rs) can be computed by maximizing the regularized cost equation 18 for the corresponding value of s ≥ 0. The corresponding optimal decoder P Y |U K for these encoders can be found as in equation ??. Different relevance-complexity pairs (∆s, Rs) on the boundary of IR sum DIB and encoders-and decoder mappings that achieve it can be found by solving equation 18 for different values of s ≥ 0 and then evaluating equation 16 and equation 17 for the obtained solution.

The optimization of equation 18 generally requires to compute marginal distributions involving the descriptions U1, . . . , UK , an aspect which can be non-easy computationally costly. To overcome this limitation, in the following we derive a tight variational bound on Ls(P) which lower bounds the DIB cost function with respect to some arbitrary distributions. Let us consider the arbitrary decoder

Q Y |U 1 ,...,U K (y|u1, . . . , uK ) for y ∈ Y, u1 ∈ U1, . . . , uK ∈ UK , the K decoders Q Y |U k (y|u k ) for k ∈ K for y ∈ Y, u k ∈ U k , and latent variable priors QU k (u k ), k ∈ K, u k ∈ U k .

For short, we denote

Q := {Q Y |U 1 ,...,U K , Q Y |U 1 , . . . , Q Y |U K , QU 1 , . . . , QU K }.

Let us define the variational DIB cost function L VB s (P, Q) as

EQUATION

The following lemma states that L VB s (P, Q) is a lower bound to Ls(P) for all distributions Q.

Lemma 1 For fixed pmfs P, we have

EQUATION

)

In addition, there exists a unique Q that achieves the maximum max Q L VB s (P, Q) = Ls(P), and is given by

EQUATION

EQUATION

where PU k , P Y |U k and P Y |U 1 ,...,U K are computed from the pmfs P.

Using the above, the optimization in equation 16 can be written in terms of the variational DIB cost function as

EQUATION

)

We close this section by noting that the cost function equation 19 can be seen as a generalization of the evidence lower bound (ELBO) as given in Rezende et al. (2014) ; Kingma and Welling (2013) for the single-encoder learning to the distributed setting. Also, in the specific case in which Y = (X1, . . . , XK ) the bound generalizes the ELBO used for VAEs to the case of an arbitrary number of encoders.

3.4. CASE OF UNKNOWN DISTRIBUTIONS: VARIATIONAL DISTRIBUTED IB ALGORITHM:

In practice only a set of training samples {(X1,i, . . . , XK,i, Yi)} n i=1 are available. In this section, we provide a method to optimize equation 23 in this case by parametrizing the encoding and decoding distributions that are to optimize using a family of distributions whose parameters are determined by Deep Neural networks (DNNs). This allows us to formulate equation 23 in terms of the DNN parameters and optimize it by using the reparametrization trick Kingma and Welling (2013) , Monte Carlo sampling, as well as stochastic gradient descent (SGD) type algorithms.

Let F e NN,k denote the parametric family of encoding probability distributions

P U k |X k over U k for each element on X k . Each member of this collection, P U k |X k ;γ e k , is described by a parameter vector γ e k ∈ Γ e k ⊆ R l e k ,

where

Γ e k ⊆ R l e

k denotes the set of allowable parameter vectors. The parameter vector γ e k is the output of a DNN

f θ k : X k → Γ e k , with network parameters θ k ∈ Θ k ⊆ R d e

k , e.g., the weights of the network at all layers. The DNN f θ k takes X k as input and outputs the parameter vector γ e k , determining one of the probability members

EQUATION

)

For example, the family of multivariate Gaussian distributions is parametrized by the mean µ θ k and covariance matrix Σ θ k , i.e., γ k := (µ θ k , Σ θ k ). Therefore, given an observation X k , γ k := (µ θ k , Σ θ k ) is determined by the output of the DNN f θ k and F e NN,k is given by

P U k |X k ;γ k (u k |x k ) = N (u k ; µ θ k , Σ θ k ).

Similarly, for decoders Q Y |U k over Y, define the family of distributions parametrized by a vector in

EQUATION

and for the distribution

EQUATION

)

Finally, for the distributions Qϕ k (u k ) we define the family of distributions with parameter

ϕ k ∈ Ψ k ⊆ R l p k F p NN,k = QU k ;ϕ k (u k ), for u k ∈ U k : ϕ k ∈ Ψ k .

In the following, for brevity we use

P θ k (u k |x k ), Q ψ k (y|u k ), Q ψ K (y|uK) and Qϕ k (u k ) to denote the distribu- tions parametrized by the DNNs f θ k , f ψ k , f ψ K and ϕ k , respectively.

By restricting the optimization of the variational DIB cost in equation 23 to the encoder, decoder and priors within the families of distributions

EQUATION

where we use the notation θ := [θ1, . . . , θK ], φ := [φ1, . . . , φK , φK] and ϕ := [ϕ1, . . . , ϕK ] to denote the DNN and prior parameters and, the cost in equation 27 is given by

EQUATION

Next, we train the DNNs to maximize a Monte Carlo approximation of equation 27 over θ, φ, ϕ using SGD.

We use the reparameterization trick Kingma and Welling (2013) , to sample from

P θ k (U k |X k ).

In particular, we consider F e NN,k to consist of a parametric family of distributions that can be sampled by first sampling a random variable Z k with distribution PZ k (z k ), z k ∈ Z k and then transforming the samples using some function

g θ k : X k × Z k → U k parameterized by θ k , such that U k = g θ k (x k , Z k ) ∼ P θ k (U k |x k ).

The reparametrization trick reduces the original optimization to estimating θ k of the deterministic function g θ k and allows to compute estimates of the gradient using backpropagation Kingma and Welling (2013) . The variational DIB cost in equation 27 can be approximated, by sampling m independent samples

EQUATION

)

4 EXPERIMENTS: RESILIENCE TO NOISE, ROTATION AND OCCLUSION

In this experiment, we test the robustness of our method against noise, rotation and random occlusion on the MNIST dataset. Specifically, we combine two types of random occlusions: the first encoder observes a digit from the MNIST that is occluded by a square which is rotated randomly (rotation angle uniformly distributed over [-45 o , 45 o ]); and the second encoder observes a noisy version of the same digit corrupted by additive noise (noise level uniform between 0 and 3). The noisy pixels are clipped between 0 and 1, with more than 60% of the pixels occluded. These occlusions make the problem significantly more involved than the standard MNIST (for which application of our algorithm leads to an relevance of about 99.9%).

We considered a CNN deterministic networks with dropout which achieves a 99.8% for test data on the clean MNIST data. Then, we have trained the same CNN architecture for each of the noisy inputs to the encoders, resulting in a relevance of 92.1% from the input to encoder 1 (randomly rotated occlusion) and 79.68% from the input to encoder 2 (noisy clipped image). We applied our D-VIB algorithm of Section 3.4 to this model with the CNN architecture of Table 1 , in which Encoder k = 1, 2 is parametrized by an nu k = 256 dimensional multivariate Gaussian distribution N (µ e k , Σ e k ) determined by the output of a DNN f θ k consisting of the concatenation of convolution, dense and maxpool layers with ReLu activations and dropout. The output of the last layer is followed by a dense layer without activation that generate µ e k and Σ e k . The prior is chosen as Qϕ k (u) = N (0, I). Each decoder takes the samples from P θ k (U k |X k ) and processes its inputs with a dense layer DNN (f φ K and f φ k ) each with 256 neurons and ReLu activation, which outputs a vector ŷi of size |Y| = 10 normalized with a softmax, corresponding to a distribution over the one-hot encoding of the digit labels {0, . . . , 9} from the K observations,

EQUATION

EQUATION

where Softmax(p) for p ∈ R d is a vector with i-th entry as [Softmax(p)]i = exp(pi)/ d j=1 exp(pj). Figure 4 shows the relevance-complexity tradeoffs obtained using our D-VIB algorithm of Section 3.4, with n = 50.000 and 15 distinct s-values randomly chosen in the range [10 -10 , 1]. For comparison, we also present the performance obtained using three methods among state-of the-art multiview learning approaches: (i) applying a deterministic CNN on the two views concatenated (deterministic CNN), (ii) applying the singleencoder variational IB method of Alemi et al. on the two views concatenated (C-VIB), and (iii) learning one function for each view via a distinct CNNs and optimize all CNNs independently (independent CNNs). The achieved relevance is reported in Table 2 . For other experimental results, see the appendices section.

We also mention that at a high level our algorithm D-VIB can be considered as performing some form of coregularization (for instance its Gaussian version is similar to the CCA of Hardoon et al. (2004) ). Comparatively, the single-view algorithm C-VIB can be viewed as belonging to the family of co-training style algorithms (such as the co-EM of Nigam and Ghani (2000) ) which, as mentioned in the recent survey Zhao et al. (2017) , override on single-view algorithms. The performance of D-VIB dominates that of C-VIB, which itself dominates co-EM. and a decoding function

EQUATION

A distortion-rate tuple (D, R1, . . . , RK ) is achievable for the DM CEO source coding problem with side information if there exist a blocklength n, encoding functions { φk } K k=1 and a decoding function ψ such that

R k ≥ 1 n log M (n) k , for k = 1, . . . , K, D ≥ E log Y n , ψ( φ1(X n 1 ), . . . , φK (X n K )) .

The distortion-rate region DRCEO of the CEO model is defined as the closure of all non-negative tuples (D, R1, . . . , RK ) that are achievable.

Key to the proof of Theorem 1 is the following proposition which states that IRDIB and DRCEO can be inferred from each other.

Proposition 3 (∆, R1, . . . , RK ) ∈ IRDIB if and only if H(Y ) -∆, R1, . . . , RK ∈ DRCEO.

Proof: Let, for k = 1, . . . , K, J k = φ k (X n k ) and J = (J1, . . . , JK ). Then,

EQUATION

= DKL(P (y n |j) P ( The result of Theorem 1 follows easily by combining (Courtade and Weissman, 2014, Theorem 10) , which provides a single-letter characterization of the rate distortion region DR CEO of the CEO problem, and Proposition 3.

EQUATION

5.3. PROOF OF THEOREM 2:

The proof of the direct part of Theorem 2 follows by evaluating the region of Theorem 1 with the choice T = ∅ and p(u

k |x k , t) = CN (x k , Σ 1/2 k (Ω k -I)Σ 1/2 k ).

The proof of the converse part is as follows. Fix t ∈ T , S ⊆ K and a family of distributions {p(u k |x k , t)} K k=1 such that the joint distribution factorizes as equation 13. Also, let 0

EQUATION

)

Such Ω k,t always exists since

EQUATION

)

Then, we have

EQUATION

where the inequality is due to Lemma 2; and equation 41 is due to equation 39.

Also, we have

EQUATION

where equation 42 follows by using Lemma 2; and equation 43 holds by using the following equality

EQUATION

)

the proof of which uses a connection between MMSE and Fisher information as shown next.

For the proof of equation 44, first note that from the MMSE estimation of Gaussian random vectors El Gamal and Kim (2011) , we have

EQUATION

)

where

EQUATION

)

Note that ZSc is independent of YSc due to the orthogonality principle of the MMSE and its Gaussian distribution. Hence, it is also independent of USc,q. We have mmse

EQUATION

where equation 47 follows since the cross terms are zero due to the Markov chain (U k,t , X k ) --Y --(U , X K/k ); and equation 48 follows due to equation 39 and G k . Finally,

EQUATION

where equation 49 is due to Lemma 3; equation 50 is due to equation 48; and equation 51 follows due to equation 46.

Now, let Ωk := t∈T p(t)Ω k,t . The rest of the converse proof follows by averaging over the time sharing random variable to get

EQUATION

where equation 52 follows from the concavity of the log-det function and Jensen's inequality. Similarly to equation 52, from equation 43 and Jensen's Inequality we have

EQUATION

)

Finally, using equation 52 and equation 53 in equation ??, noting that

Ω k = t∈T p(t)Ω k,t Σ -1 k since 0 Ω k,t Σ -1

k , and taking the union over Ω k satisfying 0 Ω k Σ -1 k , completes the proof of the converse part; and, hence, that of Theorem 2.

Finally, equation 70 is valid for any Rsum ≥ 0 and s ≥ 0. Given s, and hence (∆s, Rs), choosing R = Rs yields ∆(Rs, PX K ,Y ) ≤ ∆s. Together with equation 66, this completes the proof of Proposition 2.

5.6. PROOF OF LEMMA 1:

The proof follows by deriving the following bounds. For any conditional pmf Q Y |Z (y|z), y ∈ Y and z ∈ Z, e.g., Z = UK or Z = U k , proceeding similarly to equation 38 and averaging over Z, we have

H(Y |Z) = E[-log Q Y |Z (Y |Z)] -DKL(P Y |Z Q Y |Z ).

(71)

Similarly, we have

I(X k ; U k ) = H(U k ) -H(U k |X k ) (72) = E[-log QU k (U k )] -DKL(PU k QU k ) -H(X k |UK ) (73) = DKL(P Y |U k QU k ) -DKL(PU k QU k ) (74)

Thus, we get

EQUATION

)

where equation 75 holds by the non-negativity of relative entropy: and the equality is met if and only if Q * is as given by equation 21 and equation 22.

6. OTHER EXPERIMENTAL RESULTS (REGRESSION FOR UNKNOWN GAUSSIAN MODEL):

6.1 D-VIB ALGORITHM FOR VECTOR GAUSSIAN MODEL

For the vector Gaussian data model equation 14 the optimal distributions P and Q in equation 23 lie within the family of multivariate Gaussian distributions. Motivated by this observation, we consider the following parameterization for k ∈ K:

EQUATION

EQUATION

EQUATION

where µ e k , Σ e k are the output of a DNN f θ k with input X k that encodes the observations in a nu k -dimensional Gaussian distribution, µ d K , Σ d K are the outputs of a DNN f φ K with inputs U1, . . . , UK , sampled from

P θ k (u k |x k ), and µ d k , Σ e k are the output of a DNN f φ k with input U k , k = 1, . . . , K.

With the above choice of parametric encoders and decoders, and using a single sample m = 1, the empirical DIB cost in equation 29 is given for the sample (x1,i, . . . , xK,i, yi) by

L emp s,i (θ, φ, ϕ) := - 1 2 (yi -µ d 12,i ) T Σ d,-1 12,i (yi -µ d 12,i ) + log det(Σ d 12,i ) -s K k=1 1 2 (yi -µ d k,i ) T Σ d-1 k,i (yi -µ d k,i ) + log det(Σ d k,i ) -s K k=1 1 2 (µ e k,i -I) T (µ e k,i -I) + log |Σ e,-1 k,i | -nu k + tr{Σ e k,i } - ny 2 (1 + sK) log(2π),

where (µ d 12,i , Σ d 12,i ) denote the output of the DNN f φ K for the i-th sample (x1,i, . . . , xK,i, yi), and similarly for the other mean and covariance terms; and where we have used that each term in the empirical DIB cost equation 29 can be computed noting that for d-dimensional Gaussian pmfs N (y; µ, Σ) we have

log N (y; µ, Σ) = - 1 2 (y -µ) T Σ -1 (y -µ) + d log(2π) + log det(Σ) ,

and the KL divergence between two multivariate Gaussian pmfs P1 ∼ N (µ1, Σ1) and

EQUATION

The multivariate Gaussian parametrization of the encoders, decoders and prior distribution as given by equation 76-equation 79 can be used for other data models that are not necessary Gaussian. For example, it is particularly suitable for regression problems in which Y lies on a continuous space. Also, it is very often used in conjunction with VAE generative problems Rezende et al. (2014) ; Kingma and Welling (2013) .

6.2. REGRESSION FOR VECTOR GAUSSIAN DATA MODEL:

Consider a distributed learning model with K = 2 encoders, each observing a noisy version of an ny-dimensional Gaussian vector Y ∼ N (y; 0, I), as

X k = H k Y + N k , where H k ∈ R n k ×ny and the noises are distributed as N k ∼ N (0, I) for k = 1, 2.

For this model, the optimal accuracy-complexity region can be computed using Theorem 2. In what follows, we evaluate the performance of our D-VIB of the previous section for regression. The algorithm is trained using a dataset of n i.i.d. samples {(X1,i, X2,i, Yi)} n i=1 form the described vector Gaussian data model. We train the DNNs for various values of the parameter s. We use the multivariate Gaussian parameterization in equation 76-equation 79 for the DNNs architecture shown in Table 6 .2. Specifically, Encoder k, k = 1, 2, consists of three dense layers of 512 neurons each followed by rectified linear unit (ReLu) activations. The output of encoder k is processed by a dense layer without nonlinear activation to generate µ e k and Σ e k of size 512 and 512 × 512, respectively. Each decoder consists of two dense layers of 512 neurons with ReLu activations. The output of decoder 1, 2 and 12 is processed, each, by a fully connected layer without activation to generate Figure 5 shows the optimal relevance-complexity region of tuples (∆, Rsum) obtained from Theorem 2 for a vector Gaussian model with K = 2 encoders, target variable dimension ny = 1, and observations dimension n1 = n2 = 3. A set of 40.000 samples split among training (30.000 samples) and test (10.000 samples). The figure depicts all accuracy-complexity pairs obtained by application of our algorithm D-VIB to this setting. The results are compared to the case of inference with known joint distribution (referred to as D-IB, see next section) as well as the case of centralized inference (C-IB). For the D-VIB algorithm, the the DNN architecture for the coders is shown in Table 6 .2. Figure 6 shows the evolution of the associated mean squared error (MSE) in the estimation of the label Y using our D-VIB algorithm. As it can bee seen from both figures the performance of our D-VIB algorithm (which does not require knowledge of the joint label-feature distribution) is very close to that predicted by the theory, i.e., our Theorem 2.

Figure 7 shows similar curves for ny = 2, n1 = n2 = 3 dimensions, for various sizes of the training datset.

As expected large training sets allow a more accurate prediction. Noteworthy, that the performance during the training phase might be better than that of the centralized learning scenario is an indicator can be caused by overfitting. Related to this aspect, recall that although the D-VIB algorithm does not estimate the underlying distribution explicitly, intuitively it does for the computation of the cost function. This is related to that universal compressors also learn the actual distribution of the data that is being compressed. Recall that since the plug-in estimator of entropy is biased downward, estimations of the mutual information terms that are involved in the cost function are then biased upward, which is an alternate explanation to the observed overfitting during the training phase. In this section, we derive an iterative method to optimize the variational DIB cost function in equation 23 when the data model is discrete and the joint distribution PX K ,Y is either known, or a good estimation of it can be obtained from the training samples. In these cases, the maximizing distributions P, Q of the variational DIB cost in equation 23 can be efficiently found by an alternating optimization procedure over P and Q similar to the expectation-maximization (EM) algorithm Dempster et al. (1977) and the standard Blahut-Arimoto (BA) methodBlahut (1972) . An extension to the vector Gaussian data model, which involves random variable with continuous alphabets, is also provided. The main idea of the algorithm is that at iteration t, the optimal distributions P (t) that maximize the variational D-IB bound L VB s (P, Q (t) ) for fixed Q (t) can be optimized in closed form and, next, the maximizing pmfs Q (t) for given P (t) can be also found analytically. So, starting from an initialization P (0) and Q (0) the algorithms performs the following computations successively and in this order, until convergence,

EQUATION

We refer to such algorithm as "Blahut-Arimoto Distributed Information Bottleneck Algorithm (BA-DIB)". Algorithm 1 describes the steps taken by BA-DIB to successively maximize L VB s (P, Q) by solving a concave optimization problem over P and over Q at each iteration. We have the following lemma whose proof follows essentially by using the log-sum inequality Cover and Thomas (1991) and the convexity of the mapping x → x log x.

Lemma 4 The function L VB s (P, Q) is concave in P and in Q.

For fixed P (t) , the optimal Q (t) maximizing the variational D-IB bound in equation 19 follows from Lemma 1 as given by equation 21-equation 22. For fixed Q (t) , the optimal P (t) can be found using the following lemma.

Lemma 5 For fixed Q, there exists a P that achieves the maximum max P L VB s (P, Q), where

EQUATION

for u k ∈ U k and x k ∈ X k , k ∈ K, and where we define

EQUATION

Proof: Due to its concavity, to maximize L VB s (P, Q) with respect to P for given Q, we add the Lagrange multipliers λx k ≥ 0 for each constraint u k ∈U k p(u k |x k ) = 1 with x k ∈ X k . For each s, λx k ≥ 0 and p(u k |x k ) can be explicitly found by solving the KKT conditions, e.g.,

∂ ∂p(u k |x k )   L VB s (P, Q) + x k ∈X k λx k   u k ∈U k p(u k |x k ) -1     = 0.

This completes the proof. 

|x k ) = 1 |U k | for u k ∈ U k , x k ∈ X k , k = 1, .

. . , K. 4: repeat 5: Compute Q (t+1) using equation 21 and equation 22. 6: Compute P (t+1) using equation 82. 7: t ← t + 1 8: until convergence.

7.1.1. CONVERGENCE:

Algorithm 1 essentially falls into the class of the Successive Upper-Bound Minimization (SUM) algorithms Razaviyayn et al. (2013) in which L VB s (P, Q) acts as a globally tight lower bound on Ls(P). Algorithm 1 provides a sequence P (t) for each iteration t, which converges to a stationary point of the optimization problem equation 23.

Proposition 4 Every limit point of the sequence P (t) generated by Algorithm 1 converges to a stationary point of equation 23.

Proof: Let Q * (P) = arg max Q L VB s (P, Q).

Using Lemma 1, for every P = P, it holds that

EQUATION

Since Ls(P) and L VB s (P, Q * (P )) satisfy the assumptions of (Razaviyayn et al., 2013 , Proposition 1), then L VB s (P, Q * (P )) satisfies A1-A4 in Razaviyayn et al. (2013) . Convergence to a stationary point of equation 23 follows from (Razaviyayn et al., 2013 , Theorem 1).

The self consistent equations equation 21, equation 22 and equation 83 satisfied by any stationary point of the D-IB problem extend those of the standard point-to-point IB problem Globerson and Tishby (2004) to the distributed IB problem with K ≥ 2 encoders. In particular, note the additional divergence term in equation 83.

7.2. GAUSSIAN SETTING:

Recall Algorithm 1. For finite alphabet sources the updating rules of Q (t+1) and P (t+1) in Algorithm 1 are relatively easy, but they become unfeasible for continuous alphabet sources. We leverage on the optimality of Gaussian test channels, shown in Theorem 2, to restrict the optimization of P to Gaussian distributions, which are easily represented by a finite set of parameters, namely mean and covariance. We show that if P (t) are Gaussian distributions, then P (t+1) are also Gaussian distributions, which can be computed with an efficient update algorithm of its representing parameters. In particular, if at time t the k-th distributions P

EQUATION

where

Z t k ∼ CN (0, Σ z t k

), we show that at t + 1, for P (t+1) updated as in equation 82, the encoder P

EQUATION

A t+1 k = Σ z t+1 k 1 + 1 s Σ -1 u t k |y A t k (I -Σ x k |y Σ -1 x k ) - 1 s Σ -1 u t k |u t K\k A t k (I -Σ x k |u t K\k Σ -1 x k ) . (87)

The detailed update procedure is given in Algorithm 2 (see the following section for the details of the derivations).

Algorithm 2 BA-DIB algorithm for the Gaussin Vector D-IB 

EQUATION

6: Compute Σ z t+1 k as in equation 86 for k ∈ K.

7: Compute A t+1 k as equation 87, k ∈ K. 8: t ← t + 1. 9: until convergence.

7.2.1. DERIVATION OF ALGORITHM 2:

We derive the update rules of Algorithm 2 and show that the Gaussian distribution is invariant to the update rules in Algorithm 1, in line with Theorem 2. First, we recall that if (X1, X2) are jointly Gaussian, then

EQUATION

where

µ x 2 |x 1 := K x 2 |x 1 x1, with K x 2 |x 1 := Σx 2 ,x 1 Σ -1 x 1 .

Then, for Q (t+1) computed as in equation 21 and equation 22 from P (t) , which is a set of Gaussian distributions, we have

EQUATION

Y|u K = CN (µ y|u t K , Σ y|u t

K:

).

Next, we look at the update P (t+1) as in equation 82 from given Q (t+1) . First, we have that p(u 

and that for two generic multivariate Gaussian distributions P1 ∼ CN (µ1, Σ1) and P2 ∼ CN (µ2, Σ2) in C N , the KL divergence is computed as in equation 80 below.

Applying equation 91 and equation 80 in equation 83 and noting that all involved distributions are Gaussian, it follows that ψs(u t k , x k ) is a quadratic form. Then, since p(u t k ) is Gaussian, the product log(p(u t k ) exp(-ψs(u t k , x k ))) is also a quadratic form, and identifying constant, first and second order terms, we can write

EQUATION

where Z(x k ) is a normalization term independent of u k , ).

Σ -1 z t+1 k = Σ -1 u t k + K H y|u t k Σ -

Next, we simplify equation 93 and equation 94 to obtain the update rules equation 86 and equation 87. From the matrix inversion lemma, similarly to Chechik et al. (Feb. 2005) , for (X1, X2) jointly Gaussian we have

EQUATION

)

Applying equation 95, in equation 93 we have

EQUATION

EQUATION

where equation 97 is due to the Markov chain U k --Y --U K\k .

Then, also from the matrix inversion lemma, we have for jointly Gaussian (X1, X2),

EQUATION

)

Applying equation 98 to equation 94, for the first term in equation 94, we have

EQUATION

where Σ y,u 

:

Alternatively, the relevance could be defined in a more operational manner by the average logarithmic loss distortion or errorEP X K ,Y [ log (Y n , Ŷ n )] = H(Y n | Ŷ n ).

Figures and Tables:

Figure 3: View 1: occluded. View 2: noisy.

Figure 4: relevance v.s. sum-complexity for n = 50.000 and s ∈ [10 -10 , 1].

where equation 38 is due to the non-negativity of the Kullback-Leibler divergence and the equality holds if and only if for P (y n |j) = P (y n |j) where P (y n |j) = Pr{Y n = y n |J = j} for all j and y n ∈ Y n . Let an achievable tuple (∆, R1, . . . , RK ) ∈ IRDIB be given. Then, there must exist functions {φ k } K k=1 such that equation 9 and equation 10 hold. Using equation 38 that by letting the decoding function ψ(JK) = {P Y n |J K (y n |JK)}, we have E[ log (Y n , Ŷ n )|JK] = H(Y n |JK), which implies (H(Y ) -∆, R1, . . . , RK ) ∈ DRCEO.

Figure 5: Accuracy vs. sum-complexity tradeoffs for the Gaussian data model. Parameters: K = 2 encoders, ny = 1, n1 = n2 = 3 and n = 30.000.

Figure 7: Effect of varying training set size n = {5.000, 10.000, 50.000} on Gaussian inference for D-VIB.

Used CNN architecture.

Table content: <table/>

Achieved relevance levels.

Table content: <table><tr><td/><td>relevance (%)</td><td/></tr><tr><td/><td>1 shot</td><td>avg.</td></tr><tr><td>D-VIB</td><td>96.16</td><td>97.24</td></tr><tr><td>D-VIB-noReg</td><td>96.04</td><td>96.72</td></tr><tr><td>C-VIB</td><td>96.01</td><td>96.68</td></tr><tr><td>Deterministic CNN</td><td>93.18</td><td>93.18</td></tr><tr><td>Independent CNNs</td><td>92.1 / 79.68</td><td>93.1 / 82.01</td></tr></table>

Used DNN architecture.

Table content: <table><tr><td/><td>DNN Layers</td></tr><tr><td>Encoder k</td><td>dense [512]-ReLu</td></tr><tr><td/><td>dense [512]-ReLu</td></tr><tr><td/><td>dense [512]-ReLu</td></tr><tr><td colspan="2">Lat. space k dense [256]-ReLu</td></tr><tr><td>Decoder 12</td><td>dense [256]-ReLu</td></tr><tr><td>Decoder k</td><td>dense [256]-ReLu</td></tr></table>

Algorithm 1 BA-DIB training algorithm for discrete data 1: inputs: discrete pmf P X1,...,X k ,Y , parameter s ≥ 0. 2: output: optimal P * U k |X k , pair (∆ s , R s ).

Table content: <table><tr><td>3: initialization</td></tr><tr><td>Set t = 0 and set P (0) with p(u k</td></tr></table>

1: inputs: covariance Σ y,x1,...,x k , parameter s ≥ 0. 2: output: optimal pairs (A * k , Σ z * k ), k = 1, . . . , K.

Table content: <table><tr><td>3: initialization</td><td/></tr><tr><td>Set randomly A 0 k and Σ z 0 k 4: repeat</td><td>0, k ∈ K.</td></tr><tr><td colspan="2">5: Compute Σ x k |u t K\k and update for k ∈ K</td></tr></table>

K\k |x k [DKL(P Y |U K\k ,x k ||Q Y |U K\k ,u k )] = DKL(P Y,U K\k |x k ||Q Y,U K\k |u k )-DKL(P U K\k |x k ||Q U K\k |u k )

Table content: <table><tr><td>marginal of U t k , given by U t k ∼ CN (0, Σ u t k</td><td>) where Σ u t k</td><td>= A t k Σx k A t,H k</td><td>+ Σ z t</td><td>t k ) is the</td></tr></table>

This shows that p (t+1) (u k |x k ) is a multivariate Gaussian distribution and that U t+1 k |{X k = x k } is also a multivariate Gaussian distributed as CN (µ u t+1

Table content: <table><tr><td/><td/><td/><td/><td/><td/><td colspan="2">+</td><td>1 s</td><td colspan="3">K y,u t K\k |u t k</td><td>Σ -1 y,u t K\k |u t k</td><td>µ y,u t K\k |x k -</td><td>1 s</td><td>K u t K\k |u t k</td><td>Σ -1 u t K\k |u t k</td><td>µ u t K\k |x k</td><td>.</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td>(94)</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td>k</td><td colspan="2">|x k</td><td>, Σ z t+1 k</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td colspan="3">1 y|u t k</td><td>K y|u t k</td><td/></tr><tr><td/><td/><td/><td>+</td><td colspan="2">1 s</td><td colspan="4">K H yu t K\k |u t k</td><td colspan="3">Σ -1 yu t K\k |u t k</td><td>K yu t K\k |u t k</td><td>-</td><td>1 s</td><td>K H u t K\k |u t k</td><td>Σ -1 u t K\k |u t k</td><td>K u t K\k |u t k</td><td>,</td><td>(93)</td></tr><tr><td>and</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>µ u t+1 k</td><td>|x k</td><td>= Σ z t+1 k</td><td colspan="2">K H y|u t k</td><td colspan="2">Σ -1 y|u t k</td><td colspan="4">µ y|x k</td><td/></tr></table>

Σx k ,y ; and equation 100 is due to the definition of Σ x k |y . Similarly, for the second term in equation 94, we have K yu t K\k ; and equation 102 is due to the Markov chain U k --Y --U K\k . Finally, we note that due to equation 85, Σ u t k |y and Σ u t k |u t K\k are given as in equation 88 and equation 89, where Σ x k |y = Σ k and Σ x k |u t K\k can be computed from its definition. This completes the proof.

Table content: <table><tr><td>t k</td><td colspan="2">= A t k K\k |u t k</td><td colspan="2">Σ -1 yu t K\k |u t k</td><td colspan="3">µ y,u t K\k |x k = Σ -1 u t k |yu t K\k</td><td>A t k (I -Σ x k |yu t K\k</td><td>Σ -1 x k )x k ,</td><td>(101)</td></tr><tr><td/><td/><td/><td/><td/><td/><td colspan="2">= Σ -1 u t k |y A t k (I -Σ x k |y Σ -1 x k )x k ,</td><td>(102)</td></tr><tr><td colspan="2">where we use Σ u t k ,yu t K\k</td><td colspan="5">= A t k Σ x k ,yu t</td></tr><tr><td colspan="5">For the third term in equation 94,</td><td/><td/></tr><tr><td/><td colspan="3">K u t K\k |u t k</td><td colspan="2">Σ -1 u t K\k |u t k</td><td colspan="2">µ u t K\k |x k = Σ -1 u t k |u t K\k</td><td>A t k (I -Σ x k |u t K\k</td><td>Σ -1 x k )x k .</td><td>(103)</td></tr><tr><td colspan="7">Equation equation 87 follows by noting that µ u t+1 k</td><td>|x k</td><td>= A t+1 k x k , and that from equation 94 A t+1 k</td><td>can be</td></tr><tr><td colspan="4">identified as in equation 87.</td><td/><td/><td/></tr></table>

Peer review (iclr2020):
I am not an expert in this area and the paper involves a lot of derivations and proofs, but I did not check the correctness of those derivations. In summary, this paper proposed a framework for integrating multiple data sources for representing data. In the framework, each data source was mapped to a latent data variable by using a nonlinear function which is called an encoder; then the mapped latent variables were jointly mapped to the target data by using another nonlinear function which is called a decoder. To make this idea to work, the paper used mutual information as the objective function to control the accuracy if the model, and at the same time to avoid overfitting the paper proposed to use MDL as a measure to control the complexity of the model. If I was right, this was the whole picture of the proposed model. My questions are the following: 1) I am not very clear how the model complexity was automatically incorporated with the objective function. It seems to me that the objective function was finally the equation (29) and then the neural networks for encoder and decoder were optimized. If this was the case, how the model complexity was incorporated, that is, how the R_k was used in the model? Was the values R_k constant in the model - I mean they are fixed constant values? How these values,i.e.,R_k, were chosen? 2) I am a mathematician, but to be honest, I feel that the Maths in the paper is huge and heavy and I thought it could not be that complex for the model. The consequence is that it make the paper to be hard to read. This is a personal feeling, you could just ignore this point. 3) Experiments: there are a lot of papers describing to integrate data sources for at least the MNIST example. It would be interesting to compare the proposed method to the literature. The experiment in 4.1 obviously is a toy data problem - I mean although the data is real, but the data generated was using noisy and rotations. It would be more interesting to apply the method to a real-world problem. 4) I think it would be more friendly to explicitly define the concepts of Discrete Memoryless and Memoryless Vector Gaussian Models. 5) The Markov chain represented in equation (3) is not well defined. I do not understand these notations. 6) Before the equation (4), is it equivalent X_k^n and X_{k,n}? I am confused by these notations 7) In equation (6), it is more readable to explicitly define the Shannon Mutual Information. 8) The second paragraph on Page 5: you use Gaussian pmfs here, but pmf denotes discrete variable. But Gaussian I assume is continuous.