Year: 2018

Abstract:
Extracting relevant information, causally inferring and predicting the future states with high accuracy is a crucial task for modeling complex systems. The endeavor to address these tasks is made even more challenging when we have to deal with high-dimensional heterogeneous data streams. Such data streams often have higher-order inter-dependencies across spatial and temporal dimensions. We propose to perform a soft-clustering of the data and learn its dynamics to produce a compact dynamical model while still ensuring the original objectives of causal inference and accurate predictions. To efficiently and rigorously process the dynamics of soft-clustering, we advocate for an information theory inspired approach that incorporates stochastic calculus and seeks to determine a trade-off between the predictive accuracy and compactness of the mathematical representation. We cast the model construction as a maximization of the compression of the state variables such that the predictive ability and causal interdependence (relatedness) constraints between the original data streams and the compact model are closely bounded. We provide theoretical guarantees concerning the convergence of the proposed learning algorithm. To further test the proposed framework, we consider a high-dimensional Gaussian case study and describe an iterative scheme for updating the new model parameters. Using numerical experiments, we demonstrate the benefits on compression and prediction accuracy for a class of dynamical systems. Finally, we apply the proposed algorithm to the real-world dataset of multimodal sentiment intensity and show improvements in prediction with reduced dimensions.

heading: 1 INTRODUCTION
text: The use of machine learning for making inference and prediction from the real-world data has shown unprecedented growth. There exist a plethora of approaches for complex system (CS) modeling (e.g., multi-input multi-output state space identification (Stoica & Jansson, 2000), expectation maximization (EM) (Martens, 2010), regularization (Chiuso, 2016), graphical models (Meinshausen & Buhlmann, 2006), combined regularization and Bayesian learning (Fox et al., 2008; 2011; Bonettini et al., 2015), kernel-based regularization (Pillonetto & Chiuso, 2015)). With the increase in the size of data, the complexity of the accurate models also increases, making inference and predictions slower. The major challenges of the upcoming era, hence, are likely to deal with the massive and diverse data sources, and still making quick decisions. Therefore, the compact modeling of time-varying complex systems1 is a challenging task and appealing for more investigation. The real-world data has complex inter-dependencies across spatial and temporal dimensions. We aim to identify such dependencies and carefully construct a compact representation of the given CS model while still ensuring accurate predictions. We do so by performing a soft-clustering of such inter-dependencies to preserve only the relevant information. For the CS model in the form of a dynamical system, we additionally argue that similar to the data, the most relevant information also gets transformed at each hop in an alternate dynamical system. From a bird’s-eye view, we track
1In this work, we interpret complex systems (networks) (Barrat et al., 2008; Gao et al., 2014) as graphs comprising of nodes interacting spatially and temporally, i.e. both inter and intra dependence, with node activities available in the form of time-series.
how the most relevant information propagates across the given dynamical system. We represent this propagation via an alternate dynamical system (compact model) and develop an unsupervised learning technique of such process.
The most relevant work in this regard is information bottleneck (IB) principle (Tishby et al., 2000). For fixed two random variables, it performs a soft-clustering to compress one variable while predicting another, given the joint probability distribution. The IB has been successfully applied to speech recognition (Hecht & Tishby, 2005), document classification (Slonim & Tishby, 2000), gene expression (Friedman et al., 2001) and deep learning (Tishby & Zaslavsky, 2015), etc., and it has shown good performance. In contrast, we aim to learn a dynamics of the soft-clustering across the given dynamical system, and propose a general optimization framework to study the trade-offs between compactness and the resulting accuracies.
The problem statement addressed in this work is: Given a dynamical system, we aim to develop a compact model by learning the dynamics of the soft-clustering in an unsupervised manner, or alternate dynamical process, through Information Bottleneck hierarchy (IBH). The main contributions of the present work are as follows: (i) By learning the dynamics of the soft-clustering, we propose an alternate compact dynamical system of the given process, with emphasis on the prediction accuracies. (ii) We formulate a novel optimization setup, compact perception problem, and characterize general solution to the information theoretic problem. (iii) We quantify how most relevant information about future gets transformed at each hop in the alternatively designed dynamical system.
A brief mention of the mathematical notations is provided in the next part.

heading: 1.1 DEFINITIONS AND NOTATIONS
text: In this manuscript, we use capital letters to denote random variables (RVs), and lowercase letters are used for the realizations. The bold letters are used for multi-variate RVs. For a RV X , with little abuse of notation, we denote the probability mass function pX(x) as p(X), unless specified otherwise. The expectation operator is denoted as E[.]. A Gaussian distributed multi-variate RV is denoted as X ∼ N (µX,ΣX), where µX and ΣX are the mean vector and covariance matrix, respectively. Next, we present a few information theoretic definitions relevant to this work. Definition 1. The Kullback-Leibler (KL) divergence (Cover & Thomas, 1991) between two probability mass functions p(·) and q(·) is written as
DKL(p||q) = ∑ x p(x) log p(x) q(x) . (1)
The KL divergence is in general not symmetric, and DKL ≥ 0 with DKL(p||q) = 0 if and only if p = q. Using (1), the mutual information between two RVs X and Y is defined as I(X;Y ) = DKL(p(X,Y )||p(X)p(Y )). Next, we state the problem statement addressed in this work.

heading: 2 PROBLEM FORMULATION
text: Complex systems consist of a large number of interacting dynamic components. In practical situations, we have high-dimensional time-series (node activities) having complex spatial and temporal correlations. The challenge lies in identifying such complicated inter-dependencies and recover the true dimensionality of the system in an unsupervised manner. We equivalently call it as identifying the compact model while preserving the relevant information across spatio-temporal states. Under stationary assumptions, (Tishby et al., 2000) proposed the IB to process information compactly. For the sake of completeness, we provide a brief overview of the IB principle.

heading: 2.1 INFORMATION BOTTLENECK APPROACH
text: The IB compresses a variableX into a new stochastic variableB via soft-clustering, while maintaining as much information as possible about another variable of interest Y . The variableB operates to minimize the information compression task and to maximize the relevant information. Provided that the relevance of one RV to another is measurable by an information-theoretic metric (i.e., mutual information), the trade-off can be written as the following variational problem
min p(B|X)
I(X;B)− βI(B;Y ), (2)
where β controls the trade-off between the tasks mentioned above. Hence, the variable B, solving the minimization problem (2), encodes the most informative part from the input X about output Y .
Inspired by the IB concept and its internal predictive coding, we propose a causal inference framework for discrete-time stochastic dynamical systems. More precisely, we label two consecutive dynamic states Xk and Xk+1 to be the input and the output of the bottleneck Bk, respectively. Hence, Bk carries the most informative part (relevant information) from the past about the future. Next, we formalize this idea to create a sequence of bottlenecks (or dynamics of soft-clustering) to compactly and accurately represent the given dynamical system.

heading: 2.2 LEARNING INFORMATION PROPAGATION: SETUP
text: We consider a stochastic dynamical system, involving large number n of random processes X = [X1, X2, ....Xn], interacting and evolving in time. A fundamental problem for decision making or prediction is to learn a compressed representation of the system dynamics from high-dimensional data (i.e., system identification). The higher-order dependencies hiding behind the high-dimensions make this task challenging, and it needs sophisticated techniques to extract meaningful information for an appropriate application. In this paper, we propose a framework that focuses only on the dynamics of the relevant information, by learning an alternate representation of the given process.
Adopting an information theoretic representation, we aim to determine Bk and Bk+1 jointly since they provide compressed predictive information about the system dynamics. Figure 1 summarizes our objectives, and we study the shaded region (3 states) of the dynamical system in isolation, i.e. without the influence of any other RVs in the complete system. An argument to generalize this study to any N -length is provided later in Section 3.1. We determine the stochastic variable Bk that not only compresses Xk−1 as much as possible while preserving the relevant information about Xk, but also delivers this information to Bk+1. The variable Bk+1 quantifies the meaningful information about the state of the system at time k + 1, building upon the compressed information received via Bk, thus forming a dynamics of the relevant information.
By construction, the new mapping Bk is designed from Xk−1 to preserve the consistent information about Xk, therefore, given Xk−1, Bk and Xk are independent. Similarly, the mapping Bk+1 is independent of Xk+1 given Bk. Since, Bk+1 carries information from Bk which is compressed representation of Xk−1, then given Bk, Xk−1 and Bk+1 are also independent. With this framework, the following Markov chains are considered: Xk–Xk−1–Bk, Xk+1–Bk–Bk+1, Xk−1–Bk–Bk+1 and Bk–Xk−1–Xk+1.
The trade-off between the compression and preservation of relevant information is defined as the minimum achievable rate I(Xk−1; Bk) subject to constraints on the information processing. We call it a compact perception problem which determines a trade-off between compression representation and predictive characteristics. Formally, this can be written as the following optimization.
min p(Bk|Xk−1),p(Bk+1|Bk) I(Xk−1; Bk)
subject to I(Xk−1; Xk)− I(Bk; Xk) ≤ 1 I(Bk; Bk+1)− I(Xk−1; Xk) ≤ 2 I(Xk−1; Xk+1)− I(Xk+1; Bk+1) ≤ 3
(3)
where the constraints characterize the bounds on the desired prediction/compression at each step. We wish to lower bound (to guarantee) the prediction accuracies via lower bounding I(Bk; Xk)
and I(Xk+1; Bk+1), and upper bound (to limit) the compression level across hop via I(Bk; Bk+1). For example, 1 bounds the accuracy of the prediction of Xk by Bk. The information flow across alternate dynamical process is controlled by 2. Lastly, 3 defines closeness of prediction of Xk+1 by Bk+1. We show in the Section 3.1, and results in Section 4, that such trade-off can be alternatively studied by choice of the Lagrange parameters.

heading: 3 LEARNING INFORMATION PROCESSING FOR DYNAMICAL SYSTEMS
text: This section provides the main results concerning solving the compact perception optimization problem in equation (3) under the most general case. Next, this general result is used to study a highdimensional continuous Gaussian distribution. Lastly, we describe an iterative method to update the corresponding parameters.

heading: 3.1 LEARNING INFORMATION PROPAGATION: SOLUTION
text: Finding the alternate dynamical representation, or Bk and Bk+1, as stated in (3) requires to solve a variational problem. To solve this problem, we introduce the Lagrange multipliers β, λ and γ for the information processing constraints. Hence, we find the alternative representation (or IBH) by minimizing the following functional written using (3): F [p(Bk|Xk−1), p(Bk), p(Xk|Bk), p(Xk+1|Bk), p(Bk+1|Bk), p(Bk+1)] =
I(Xk−1; Bk)− βI(Bk; Xk) + λ (I(Bk; Bk+1)− γI(Xk+1; Bk+1)) . (4) We can argue the following regarding the optimal solution which minimizes the equation (4). Theorem 1. The optimal solution that minimizes functional F in (4) satisfy the following selfconsistent equations:
p(Bk|Xk−1) = p(Bk)
Z1 × exp {−βDKL (p(Xk|Xk−1)||p(Xk|Bk))
−λDKL (p(Bk+1|Bk)||p(Bk+1)) −λγEBk+1|Bk [DKL (p(Xk+1|Xk−1)||p(Xk+1|Bk+1))] } , (5)
p(Bk+1|Bk) = p(Bk+1)
Z2 exp {−γDKL (p(Xk+1|Bk)||p(Xk+1|Bk+1))} , (6)
where Z1 and Z2 are normalizing partition functions.
The functional F in (4) may not be convex in the product space of the associated probability simplexes. Hence, it is difficult to obtain a global optimum, however, a stationary point (and locally optimal solution, in most of the cases) can be obtained using the following result. Corollary 1. The self-consistent equations in Theorem 1 can be used to write an iterative procedure to update the associated probabilities as equations (14)-(20). The iterative approach of Corollary 1 is detailed in the Appendix B. The idea is similar to the BlahutArimoto algorithm (Arimoto, 1972; Blahut, 1972), also observed in (Tishby et al., 2000). Finally, it remains to show the existence of a stationary point of the functional in (4). The convergence of the iterations in Corollary 1 is established through the following result. Lemma 1. The iterative procedure in Corollary 1 to minimize the functional F in (4) is convergent to a stationary point. The idea of proof (in Appendix C) is also somewhat similar to the EM McLachlan & Krishnan (1996) but with minimization of the functional. As in standard EM, in most of the cases, the stationary convergence point is local minimum (maximum in EM) of the functional. Corollary 2. The IBH solution in Theorem 1 reduces to IB in (2) upon setting Xk−1 = Xk, and β →∞. Proof. With Xk−1 = Xk, the functional of IBH in equation (4) can be written as F = (1− β)I(Xk−1; Bk) + λ (I(Bk; Bk+1)− γI(Xk+1; Bk+1)) . The functional F minimization in the limit of β → ∞ has one solution of p(Bk|Xk−1) = 1Bk=Xk−1 (such that I(Xk−1; Bk) > 0), and the problem reduces to minimize the following
F̂ = I(Xk−1; Bk+1)− γI(Xk+1; Bk+1), where λ can be dropped, as λ ≥ 0.
An advantage of using the optimization framework in (3) is that it can be generalized to any length of the given input dynamical process by properly repeating the second and third constraint. In this work, we have studied the case of length three of the input process. However, the same principles can be used to write the solution for any length N of the dynamical system (the complete Figure 1), as will be presented in the future work.

heading: 3.2 MEASURING INFORMATION FLOW FOR LINEAR DYNAMICS
text: Many applications in machine learning involve time-series with multi-dimensional observations where each dimension is very likely to be correlated with others, and the state of observations evolves in time. Assuming that these observations are corrupted by Gaussian noise, then a linear dynamical system is a promising model for analyzing the provided time-series data. Thus, the system dynamics can be modeled similarly as the evolution of a stochastic time-invariant linear system. Finally, if we have Gaussian distributed initial-state of the linear dynamical system with additive Gaussian noise at each time step, then all future states are jointly Gaussian distributed. In this case, our framework learns the IBH through Gaussian random vectors.
Recall that the states Xk of the dynamical system under study are jointly Gaussian, and without loss of generality we assume that they are centered. As explained previously, we aim to design Bk and Bk+1 to define an alternate representation which captures the dynamics of the relevant information. It is shown in the prior works of (Globerson & Tishby, 2004; Chechik et al., 2005) that for the problem setup of IB in which the input and output variables are jointly Gaussian, the optimum solution of the IB Lagrangian obtained by a stochastic transformation is also jointly Gaussian with the bottleneck’s input. Consequently, Bk is jointly Gaussian with Xk−1 and Bk+1. Since RVs in consideration are mean centered and Xk−1,Bk,Bk+1 are jointly Gaussian; the IBH variables can be very well represented as linear transformations of each other. Additionally, using the MC conditions from Section 2.2 we write the following linear relations.
Bk = ΦXk−1 + ξk,
Bk+1 = ∆Bk + ξk+1, (7)
where ξk and ξk+1 are centered Gaussian random vectors independent of Xk−1 and Bk, respectively. Given the aforementioned settings, the solution of the minimization problem in equation (3) is determined by finding the matrices Φ and ∆, and the covariance matrices Σξk and Σξk+1 . An iterative procedure using Corollary 1 to update the concerned parameters in (7) is presented as the following result. Theorem 2. Given the parameters β, λ and γ, the Gaussian bottlenecks Bk = ΦXk−1 + ξk and Bk+1 = ∆Bk + ξk+1 are obtained by performing the following iterations over the parameters
Σ (t+1) ξk
= ( βΣ
−1 (t) Bk|Xk − (β − 1) Σ −1 (t) Bk + λ∆T (t)Σ −1 (t) Bk+1
∆(t) )−1 ,
Φ(t+1) = Σ (t+1) ξk
( βΣ
−1 (t) Bk|XkΦ (t)(I −ΣXk−1|XkΣ −1 Xk−1 )
+λγ∆T (t)Σ −1 (t) Bk+1|Xk+1∆ (t)Φ(t)Σ−1Xk+1ΣXk+1Xk−1Σ −1 Xk−1
) ,
Σ (t+1) ξk+1 = ( γΣ −1 (t) Bk+1|Xk+1 − (γ − 1)Σ −1 (t) Bk+1 )−1 ,
∆(t+1) = ( I − γ − 1
γ Σ
(t) Bk+1|Xk+1Σ −1 (t) Bk+1
)−1 ∆(t) ( I −Σ(t+1)Bk|Xk+1Σ −1 (t+1) Bk ) ,
(8)
where t is the iteration index.
The detailed proof of the Theorem 2 is provided in the Appendix. In the next section, we show some numerical results generated using synthetic data.

heading: 4 SIMULATION EXPERIMENTS
text: We numerically evaluate the results of Theorem 1 and Theorem 2 using synthetically generated Gaussian distribution. Specifically, the covariance matrices for the input dynamical process
Xk−1–Xk–Xk+1 are generated numerically for a given size tuple, and we compute the parameters of (7) using (8). The quantities of interest are I(Bk+1; Xk+1) and I(Bk+1; Xk−1) which are indicators of prediction information and inverse of compression, respectively. We compare the prediction/compression behavior of the proposed approach of the alternate design of the dynamical system vs. designing local IB’s between each hop. The local IB’s are designed between two consecutive RVs in the dynamical system independently while the IBH is designed jointly. We show the distinction upon varying the dimensions of the input process (Xk−1,Xk,Xk+1) in Figure 2. It is observed that the gap between prediction (for a fixed level of compression) grows with an increase in the input dimensions. The IBH by design takes into account the entire input dynamical system, and construct an alternate representation which provides better prediction at each step, by appropriate choice of the Lagrange parameters β, λ, γ.
The Lagrange parameters (β, λ, γ) control the trade-off between compression and prediction at each step of the alternatively designed dynamical process. In the optimization problem (3), β corresponds to the first constraint, and hence plays a deciding role in prediction accuracy of Bk. The λ corresponds to the second constraint, and will control the flow of relevant information across Bk. Finally, λ and γ together tune the accuracy of prediction using Bk+1. For example, in (7), the information tapping behavior can be visualized by inspecting the ranks of Φ and ∆ matrices upon varying (β, λ, γ). In Figure 3a, the rank(∆) increases upon increasing γ for each choice of λ. Since λ appears in front of both prediction and compression expression in (4), it has little effect on the rank(∆) for a fixed γ and β. Next, in Figure 3b, we observe dynamical effects of the information flow. By fixing λ, we limit the information acceptance of Bk+1, hence the parameter β can only increase the rank(∆) up to a certain limit by allowing maximum information through Bk. We witness in Figure 3b, that with higher λ, the parameter β quickly increase the rank(∆). Now, for a fixed λ, both
β and γ intertwine with each other to decide rank(∆). By fixing β, we limit the input information through Bk to Bk+1, or availability, and hence γ can only increase the rank(∆) up to certain extent. Similarly, fixing γ limits the maximum information that Bk+1 can process, and therefore β can do best narrowly up to some extent. We witness this hyperbolic behavior in Figure 3c.

heading: 5 REAL-WORLD DATA: MULTIMODAL SENTIMENT INTENSITY
text: In this section, we apply the ideas of IBH to extract the features from the challenging multimodal datasets available in the form of time-series. Particularly, we have used the CMU Multimodal Sentiment Analysis (CMU-MOSI) dataset (Zadeh et al., 2016) which consists of a total of 2198 videos. Each video comprises one speaker expressing their opinion in front of the camera. The available modalities from the dataset are text, visual and audio. The goal of the dataset is to perform a discriminative task of predicting the speaker sentiment using the available modalities. From the three modalities of text, visual and audio, the corresponding features are extracted using GloVe word embeddings (Pennington et al., 2014), Facet (iMotions, 2017) and COVAREP (Degottex et al., 2014), respectively. The extracted feature size are 300 for text, 74 for audio, and 46 for visual component.
The extracted features in the form of time-series are aligned across modalities as shown in Figure 4. Specifically, multiple features of visual and audio modality (due to high frame/sampling rate) are time averaged with boundaries corresponding to the text component. For each speaker, we have a maximum of 20 words with three modalities and the corresponding sentiment intensity being a real number ranging from−3 to 3, with negative values representing negative sentiments and vice-versa. Some of the prior work in the multimodal representation learning include Discriminative representation learning (Zadeh et al., 2018a; 2017; 2018b; Chaplot et al., 2017) and Generative representation learning (Sohn et al., 2014; Srivastava & Salakhutdinov, 2014; Suzuki et al., 2016). Interestingly, recent work Tsai et al. (2018) proposed Multimodal Factorized Model (MFM) that exploits the fusion of both these techniques. The work factorizes the data into discriminative factors and modalityspecific generative factors. We note that, at the very core, the challenge in learning patterns from multiple modalities is to address the complex inter and well as intra dependencies across them. Since IBH is capable of compressing the given dynamical model into an alternate version stochastically, we propose to, (i) first map the multiple modalities into a time-varying linear dynamical system; and then (ii) use IBH to identify the complex inter and well as intra dependencies across them in a reduced dimensional model for better discrimination using simple machine learning classifiers.
The IBH is applied to various modalities to capture the information flow patterns across them and in the compressed fashion. We have taken three modalities to be in the following Markov Chain, text–audio–visual. The intuitive reason behind this assumption is that text is the most informative modality for sentiment and hence the first state, while audio and visual follows text in the Markov chain assuming the speaker is being honest in speaking and making visual expressions for a particular sentiment. The data is mean centered and the covariances of three modalities are estimated for each speaker in the training as well as the testing dataset. However, we have only 20 words and hence 20 samples to estimate the covariance matrix of much larger dimensions. To remedy this fewer samples problem, we have resorted to something called pooling of the covariance matrix as follows. The covariance matrix of any modality for the ith speaker is written as Σi = αΣ̂i + (1 − α)Σpool, where Σ̂i is the estimated covariance from 20 samples, and Σpool is called pooled covariance matrix
which is estimated by taking all of the training data. The parameter α is used to make a trade-off between these two matrices and is usually chosen close to 1. The covariance matrices are fed to the algorithm in Theorem 2 to estimate Φ and ∆ matrices. The Φ matrix has components for interactions across text and audio, while ∆ has entries representing interactions with the compressed version of text-audio inter-dependencies as well as with the visual modality. Therefore, the entries of ∆ matrix are good candidates for representing inter as well as intra dependencies across all three modalities. Hence, we use ∆ matrix as a feature to predict the sentiment intensity, as this matrix is central to information propagation from the text, through audio, to visual component. The low rank of ∆ will be key in reducing the number of features. The values of parameters (β, λ, γ) are chosen such that maximum information flow to first bottleneck, by setting high value for β, and λ close to 1, and low values for γ for better compression.
The results are reported for various evaluation methods, namely binary (with the positive and negative sentiment) and 7− class classification, and Mean Average Error (MAE) for regression. We compare our performance with naive early fusion of features in raw format, most recent results in the multimodal neural networks (Tsai et al., 2018) vs. features processed by IBH, and the numerical results are presented in Table 1. Support Vector Machines (SVM) is used in the cases of early fusion and IBH. The processing of features by IBH has a two-fold advantage: First, due to compression, the dimensionality of the input can be reduced from 8400 to 1150. Second, the performance with respect to various metrics is better.

heading: 6 CONCLUSION
text: In this paper, we have introduced a novel information-theoretic inspired approach to learn the compact dynamics of a time-varying complex system. The trade-off between the predictive accuracy and the compactness of the mathematical representation is formulated as a multi-hop compact perception optimization problem. A key ingredient to solve the aforementioned problem is to exploit variational calculus in order to derive the general solution expressions. Additionally, we have investigated the guaranteed convergence of the proposed iterative algorithm. Moreover, considering a specific class of distributions (Gaussian), we have provided closed-form expressions for the model parameters’ update in our algorithm. Interestingly, the proposed compact perception shows improvements in prediction with reduced dimension on challenging real-world problems.
The quantification of information flow across a dynamical system can have an enormous impact on understanding and improving the current state-of-the-art in neural networks as realized in (Tishby & Zaslavsky, 2015). Moreover, modeling with dynamical systems is a standard approach, and by using the proposed framework, we can make a better compact representation of the system. The driving force of a dynamical system can enforce different behaviors of information flow, as realized in defining dynamical entropy by (Sinai, 1959). Therefore, measuring the information flow can help in estimating/differentiating the actual driving component behind the observed activities. Such concepts are useful in predicting brain imagined tasks from observed electroencephalogram activities.

heading: A PROOF OF THEOREM 1
text: Proof. For the sake of simplicity, a sketch of the proof is given for discrete variables. The Lagrangian associated with the minimization problem is the following
L = F − ∑
Bk,Xk−1
α1(Xk−1)p(Bk|Xk−1)− ∑
Bk,Bk+1
α2(Bk)p(Bk+1|Bk), (9)
where α1(Xk−1) and α2(Bk) are Lagrange multipliers for the normalization of the distributions p(Bk|Xk−1) and p(Bk+1|Bk), respectively. Taking the derivative of each term of the Lagrangian L with respect to p(Bk|Xk−1), we have
δI(Xk−1; Bk) δp(Bk|Xk−1) = p(Xk−1) log p(Bk|Xk−1) p(Bk) ,
δI(Xk; Bk)
δp(Bk|Xk−1) = p(Xk−1) ∑ Xk p(Xk|Xk−1) log p(Xk|Bk) p(Xk) ,
δI(Bk; Bk+1) δp(Bk|Xk−1) = p(Xk−1)DKL(p(Bk+1|Bk)||p(Bk+1)),
δI(Xk+1; Bk+1)
δp(Bk|Xk−1) = −p(Xk−1) ∑ Bk+1 p(Bk+1|Bk)DKL(p(Xk+1|Xk−1)||p(Xk+1|Bk+1))
+ p(Xk−1)DKL(p(Xk+1|Xk−1)||p(Xk+1)). (10) Setting the derivative of the Lagrangian equal to zero and arranging the terms we obtain the self consistent equation (5). Note that all the constant terms in the derivative independent of Bk will be captured by the Lagrange multiplier α1(Xk−1). The derivative of the Lagrangian L with respect to p(Bk+1|Bk) involves only the two last terms from the functional F and the term that ensures the normalization condition. Then, we have
δI(Bk; Bk+1) δp(Bk+1|Bk) = p(Bk) log p(Bk+1|Bk) p(Bk+1) , (11)
δI(Xk+1; Bk+1)
δp(Bk+1|Bk) = p(Bk) ∑ Xk+1 p(Xk+1|Bk) log p(Xk+1|Bk+1) p(Xk+1) . (12)
Thus, the variational condition is written as follows
p(Bk)
( log
p(Bk+1|Bk) p(Bk+1)
+ γDKL (p(Xk+1|Bk)||p(Xk+1|Bk+1))− α2(Bk) ) = 0, (13)
where α2(Bk) is the summation of the Lagrange multiplier α2(Bk) and the terms independent of Bk+1, and hence the equation (6) follows.
B ITERATIVE SOLUTION IN COROLLARY 1
The self-consistent equations derived in Theorem 1 to minimize the functional F in (4) can be used to write the following set of iterative equations.
p(t+1)(Bk|Xk−1) = p(t)(Bk)
Z (t+1) 1 (β, λ, γ,Xk−1)
× exp { −βDKL ( p(Xk|Xk−1)||p(t)(Xk|Bk) ) −λDKL ( p(t)(Bk+1|Bk)||p(t)(Bk+1)
) −λγE(t)Bk+1|Bk [ DKL ( p(Xk+1|Xk−1)||p(t)(Xk+1|Bk+1) )]} , (14)
p(t+1)(Bk) = ∑ Xk−1 p(t+1)(Bk|Xk−1) p(Xk−1), (15)
p(t+1)(Xk|Bk) = ∑ Xk−1 p(Xk|Xk−1) p(t+1)(Xk−1|Bk), (16)
p(t+1)(Xk+1|Bk) = ∑ Xk p(Xk+1|Xk) p(t+1)(Xk|Bk), (17)
p(t+1)(Bk+1|Bk) = p(t)(Bk+1)
Z (t+1) 2 (γ,Bk)
exp { −γDKL ( p(t+1)(Xk+1|Bk)||p(t)(Xk+1|Bk+1) )} ,
(18) p(t+1)(Bk+1) = ∑ Bk p(t+1)(Bk+1|Bk) p(t+1)(Bk), (19)
p(t+1)(Bk+1|Xk+1) = ∑
Bk,Xk−1
p(t+1)(Bk+1|Bk) p(t+1)(Bk|Xk−1) p(Xk−1|Xk+1), (20)
where in equations (15)-(17), (19) and (20) we have used the Markov assumption as stated in Section 2.2.

heading: C PROOF OF LEMMA 1
text: Proof. The proof of the Lemma 1 can be divided into two parts. First, we show that the F in (4) is lower-bounded. Next, we show that each iteration using the equations (14)-(20) monotonically decrease the functional.
Lower bound: Let us consider the following alternate functional.
F̃ = I(Xk−1; Bk) + βEXk−1,Bk DKL(p(Xk|Xk−1)||p(Xk|Bk)) + λ I(Bk; Bk+1) λγ EXk−1,Bk EBk+1|Bk DKL(p(Xk+1|Xk−1)||p(Xk+1|Bk+1)). (21)
It can be readily verified that the functional F̃ ≥ 0 for given non-negative constants β, λ and γ. Also, the F̃ in (21) can be expanded as
F̃ = I(Xk−1; Bk) + βEXk−1,Bk ∑ Xk p(Xk|Xk−1)log p(Xk|Xk−1) p(Xk|Bk)
+λ I(Bk; Bk+1) + λγ EXk−1,Bk ∑
Bk+1,Xk+1
p(Xk+1|Xk−1)log p(Xk+1|Xk−1) p(Xk+1|Bk+1)
= I(Xk−1; Bk) + β(I(Xk−1; Xk)− I(Xk; Bk)) + λ I(Bk; Bk+1) +λγ(I(Xk−1; Xk+1)− I(Xk+1; Bk+1)). (22)
Since the functional F̃ in the equation (22) differs from the functional F in (4) only in constants, therefore, F is lower bounded as well. Monotonicity: For proving monotonic decrement of the functional F , we will use the formulation similar to (Arimoto, 1972). First, let us consider the following observation made in (Tishby et al., 2000).
For a given joint distribution p(X,Y ), we can write the following.
p(Y ) = arg min φ(Y )
EX DKL(p(Y |X)||φ(Y )), (23)
where the minimization is performed over the probability simplex of φ(Y ) such that the joint distribution is p(X,Y ). The functional F̃ in (22) can be written in its most general form as
F̃(p1, p2, φ1, φ2) = EXk−1,Bk|Xk−1∼p1 log p1
p(Bk)
+ βEXk−1,Bk|Xk−1∼p1 DKL(p(Xk|Xk−1)||φ1)
+ λEXk−1,Bk|Xk−1∼p1EBk+1|Bk∼p2 log p2
p(Bk)
+ λγ EXk−1,Bk|Xk−1∼p1 EBk+1|Bk∼p2 DKL(p(Xk+1|Xk−1)||φ2), (24)
where p1 = p(Bk|Xk−1), p2 = p(Bk+1|Bk), φ1 = φ1(Xk|Bk) and φ2 = φ2(Xk+1|Bk+1), therefore, F̃ reduces to the form in (22) upon setting φ1 = p(Xk|Bk) and φ2 = p(Xk+1|Bk+1). With the objective of minimizing the functional F̃ , we can write its value at iteration t as F̃ (t) = F̂(p(t)1 , p (t) 2 , φ (t) 1 , φ (t) 2 ). The iterations to minimize F̃ will involve the successive choice of tuple (p1, p2, φ1, φ2). At iteration t, let us assume that we have chosen p (t) 1 and p (t) 2 , and we define the following
G(t, t) = min φ1,φ2 F̃(p(t)1 , p (t) 2 , φ1, φ2)
= F̃(p(t)1 , p (t) 2 , φ (t) 1 , φ (t) 2 ). (25)
Using equation (23), it easily follows that φ(t)1 = p (t)(Xk|Bk) and φ(t)2 = p(t)(Xk+1|Bk+1), and hence G(t, t) = F̃(t). Now, for fixed p(t)2 , φ (t) 1 , φ (t) 2 , it can be easily realized that the F̃ is convex in p1, therefore, minimizing F̃ with respect to (w.r.t.) p1 will involve writing Lagrangian, and then differentiation, and setting to zero. This step is similar to the Theorem 1, and we can write that
p (t+1) 1 = arg min
p1
F̃(p1, p(t)2 , φ (t) 1 , φ (t) 2 ), (26)
where the resulting solution is (5), and hence p(t+1)1 will have the expression as in (14). Similarly, for fixed p(t+1)1 , φ (t) 1 , φ (t) 2 , the F̃ is convex in p2, and the same steps follow to obtain p (t+1) 2 as written in (18). It should be noted that the choice to perform minimization w.r.t. p1 before p2 is arbitrary, and the reverse can also be performed. This will change the order of iteration index in equations (14)-(20) accordingly. Using equations (25) and (26), we can conclude that
F̃ (t+1) = G(t+ 1, t+ 1) ≤ G(t+ 1, t) ≤ G(t, t) = F̃ (t),
and therefore, iterating equations (14)-(20) written using the self-consistent equations of Theorem 1 minimizes F̃ (t) monotonically. Since F̃ and F differs only in constant, it reduces F as well monotonically from above.

heading: D PROOF OF THEOREM 2
text: Proof. For a multivariate random variable X, X ∈ Rn with Gaussian distribution, i.e. X ∼ N (µ,Σ), the entropy can be written as
H(X) = 1
2 log det(Σ) + c, (27)
where c is constant for a given dimension n. Using equation (27), the KL-divergence between two Gaussian distributed random variables, X1 ∼ N (µ1,Σ1) and X2 ∼ N (µ2,Σ2) of the same dimensions, is written as
DKL(p(X1)||p(X2)) = 1
2 log det(Σ2) det(Σ1) + 1 2 tr(Σ−12 Σ1) + 1 2 (µ2 − µ1)TΣ−12 (µ2 − µ1). (28)
We have assumed that the given data is centered, hence all considered random variables will have zero mean, i.e., for characterizing each random variable, we only need the corresponding covariance matrix. Let us revisit the linear transformation model for the IBH.
Bk = ΦXk−1 + ξk,
Bk+1 = ∆Bk + ξk+1.
Now, to completely specifying the model, we have to determine the constant matrices Φ,∆ and the covariances of ξk, ξk+1. Since entropy is well defined for Gaussian distribution, as in (27), due to their nice tail distribution, we can use Theorem 1 and equation (28) to write the self-consistent
transition probabilities. First, expressions of the necessary KL-divergence terms are expanded as follows.
DKL(p(Xk|Xk−1)||p(Xk|Bk)) = c+ 1
2 (EXk|Xk−1 − EXk|Bk)TΣ−1Xk|Bk
(EXk|Xk−1 − EXk|Bk), (29)
DKL(p(Bk+1|Bk)||p(Bk+1)) = c+ 1
2 (EBk+1|Bk)TΣ−1Bk+1EBk+1|Bk, (30)
DKL(p(Xk+1|Xk−1)||p(Xk+1|Bk+1)) = c+ 1
2 (EXk+1|Xk−1 − EXk+1|Bk+1)TΣ−1Xk+1|Bk+1
(EXk+1|Xk−1 − EXk+1|Bk+1), (31)
where c are different constants at each step and,
EXk|Bk = ΣXkBkΣ −1 Bk Bk = ΞkBk,
ΣXk|Bk = ΣXk −ΣXkBkΣ −1 Bk ΣBkXk
= ΣXk −ΣXkXk−1ΦTΣ −1 Bk ΦΣXk−1Xk ,
EBk+1|Bk = ∆Bk, ΣBk+1|Bk = Σξk+1 ,
EXk+1|Bk+1 = ΣXk+1Xk−1ΦT∆TΣ −1 Bk+1 Bk+1 = Θk+1Bk+1,
ΣXk+1|Bk+1 = ΣXk+1 −ΣXk+1Xk−1Φ T∆TΣ−1Bk+1∆ΦΣXk−1Xk+1 .
Also, using (31), we can write that
EBk+1|Bk DKL (p(Xk+1|Xk−1)||p(Xk+1|Bk+1)) = c − EBk+1|Bk B T k+1Θ T k+1Σ −1 Xk+1|Bk+1EXk+1|Xk−1
+ 1
2 EBk+1|Bk B T k+1Θ T k+1Σ −1 Xk+1|Bk+1Θk+1Bk+1
= c′ −BTk∆TΘTk+1Σ−1Xk+1|Bk+1EXk+1|Xk−1. (32)
Upon substituting the equations (29)-(31) and (32) in Theorem 1, we obtain the following selfconsistent equations.
log p(Bk|Xk−1) = c− 1
2 BTkΣ −1 Bk
Bk − β
2 BTkΞ T kΣ −1 Xk|BkΞkBk
− λ 2 BTk∆ TΣ−1Bk+1∆Bk + βB T kΞ T kΣ −1 Xk|BkEXk|Xk−1 + λγBTk∆ TΘTk+1Σ −1 Xk+1|Bk+1EXk+1|Xk−1. (33)
Now, Bk|Xk−1 ∼ N (ΦXk−1,Σξk), therefore upon comparing terms, we obtain the following self-consistent equations.
Σ−1ξk = Σ −1 Bk + βΞTkΣ −1 Xk|BkΞk + λ∆ TΣ−1Bk+1∆,
Φ = Σξk(βΞ T kΣ −1 Xk|BkΣXkXk−1Σ −1 Xk−1
+λγ∆TΘTk+1Σ −1 Xk+1|Bk+1ΣXk+1Xk−1Σ −1 Xk−1 ). (34)
Setting up the iterations, similar to Blahut-Arimoto equations, we obtain the following iterative procedure with t as iteration index.
Σ (t+1) ξk = ( Σ −1 (t) Bk + βΞ T (t) k Σ −1 (t) Xk|BkΞ (t) k + λ∆ T (t)Σ −1 (t) Bk+1 ∆(t) )−1
(a) = ( βΣ
−1 (t) Bk|Xk − (β − 1) Σ −1 (t) Bk + λ∆T (t)Σ −1 (t) Bk+1
∆(t) )−1 ,
Φ(t+1) = Σ (t+1) ξk (βΞ T (t) k Σ −1 (t) Xk|BkΣXkXk−1Σ −1 Xk−1
+λγ∆T (t)Θ T (t) k+1 Σ −1 (t) Xk+1|Bk+1ΣXk+1Xk−1Σ −1 Xk−1 )
(b) = Σ (t+1) ξk (βΣ −1 (t) Bk|XkΦ (t)(I −ΣXk−1|XkΣ −1 Xk−1 )
+λγ∆T (t)Σ −1 (t) Bk+1|Xk+1∆ (t)Φ(t)Σ−1Xk+1ΣXk+1Xk−1Σ −1 Xk−1 ), (35)
where in (a) and (b) we have used the expansion of Ξk,Θk+1, matrix inversion lemma (Golub & Van Loan, 1996), and the identity ΣX|Y ΣXY Σ −1 Y = Σ −1 X ΣXY Σ −1 Y |X , where X,Y and Z are Normal distributed random variables. We can also write the transition probability of the bottleneck variables, Bk+1|Bk using Theorem 1 and the following useful expansion of the required KL-divergence term.
DKL(p(Xk+1|Bk)||p(Xk+1|Bk+1)) = c+ 1
2 (EXk+1|Bk − EXk+1|Bk+1)TΣ−1Xk+1|Bk+1
(EXk+1|Bk − EXk+1|Bk+1), (36)
where,
EXk+1|Bk = ΣXk+1Xk−1ΦTΣ −1 Bk Bk = ΥkBk,
ΣXk+1|Bk = ΣXk+1 −ΣXk+1Xk−1Φ TΣ−1BkΦΣXk−1Xk+1 .
After substituting equation (36) in (6) of Theorem 1, we obtain
log p(Bk+1|Bk) = c− 1
2 BTk+1Σ −1 Bk+1
Bk+1 − γ
2 BTk+1Θ T k+1Σ −1 Xk+1|Bk+1Θk+1Bk+1
+ γBTk+1Θ T k+1Σ −1 Xk+1|Bk+1ΥkBk. (37)
Again, since Bk+1|Bk ∼ N (∆Bk,Σξk+1), after comparing the terms on both sides of the equation (37), we obtain
Σ−1ξk+1 = Σ −1 Bk+1 + γΘTk+1Σ −1 Xk+1|Bk+1Θk+1,
∆ = γΣξk+1Θ T k+1Σ −1 Xk+1|Bk+1Υk. (38)
After setting up the iterations like we did for equation (34), and using the (18)-(20) we first write that
Σ (t) Bk+1 = ∆(t)Σ (t+1) Bk ∆(t) + Σ (t) ξk+1 ,
Σ (t) Bk+1|Xk+1 = ∆ (t)Σ (t+1) Bk|Xk+1∆ (t) + Σ (t) ξk+1 ,
Σ (t) Bk+1Xk+1 = ∆(t)Σ (t+1) BkXk+1 . (39)
Using (38) and (39), we obtain the following.
Σ (t+1) ξk+1 = ( γΣ −1 (t) Bk+1|Xk+1 − (γ − 1)Σ −1 (t) Bk+1 )−1 ,
∆(t+1) = γΣ (t+1) ξk+1 Σ −1 (t) Bk+1|Xk+1Σ (t) Bk+1Xk+1 Σ−1Xk+1ΣXk+1Xk−1Φ T (t)Σ −1 (t+1) Bk
(a) = ( I − γ − 1
γ Σ
(t) Bk+1|Xk+1Σ −1 (t) Bk+1
)−1 ∆(t) ( I −Σ(t+1)Bk|Xk+1Σ −1 (t+1) Bk ) , (40)
where (a) is written using the previously used identity ΣX|Y ΣXY Σ −1 Y = Σ −1 X ΣXY Σ −1 Y |X , and matrix inversion lemma is used at each step.

Review summaries:
This paper studies the problem of compactly represent the model of a complex dynamic system while preserving information. The method is based on the information bottleneck method. Basically, for a dynamic system whose states changing from X_{k-1}, X_k to X_{k+1}, the "information bottleneck hierarchy" method learns a variable B_k and B_{k+1} such that B_k predicts B_{k+1} well, B_k predicts X_k well, and B_{k+1} predicts X_{k+1} well, while minimizing the information of X_{k-1} contained in B_k. 

This paper studied an extension of the Information Bottleneck Principle called Information Bottleneck Hierarchy (IBH).  The goal of IBH is to extract meaningful information from a Markov Chain. Then the authors studied case of the Gaussian linear dynamic and proposed an algorithm for computing the IBH. Then an experiment was conducted to show the usage of IBH to practical problems.