Year: 2020

Abstract:
Modelling highly multi-modal data is a challenging problem in machine learning. Most algorithms are based on maximizing the likelihood, which corresponds to the M(oment)-projection of the data distribution to the model distribution. The M-projection forces the model to average over modes it cannot represent. In contrast, the I(nformation)-projection ignores such modes in the data and concentrates on the modes the model can represent. Such behavior is appealing whenever we deal with highly multi-modal data where modelling single modes correctly is more important than covering all the modes. Despite this advantage, the I-projection is rarely used in practice due to the lack of algorithms that can efficiently optimize it based on data. In this work, we present a new algorithm called Expected Information Maximization (EIM) for computing the I-projection solely based on samples for general latent variable models, where we focus on Gaussian mixtures models and Gaussian mixtures of experts. Our approach applies a variational upper bound to the I-projection objective which decomposes the original objective into single objectives for each mixture component as well as for the coefficients, allowing an efficient optimization. Similar to GANs, our approach employs discriminators but uses a more stable optimization procedure, using a tight upper bound. We show that our algorithm is much more effective in computing the I-projection than recent GAN approaches and we illustrate the effectiveness of our approach for modelling multi-modal behavior on two pedestrian and traffic prediction datasets.

heading: 1 INTRODUCTION
text: Learning the density of highly multi-modal distributions is a challenging machine learning problem relevant to many fields such as modelling human behavior (Pentland & Liu, 1999). Most common methods rely on maximizing the likelihood of the data. It is well known that the maximum likelihood solution corresponds to computing the M(oment)-projection of the data distribution to the parametric model distribution (Bishop, 2006). Yet, the M-projection averages over multiple modes in case the model distribution is not rich enough to fully represent the data (Bishop, 2006). This averaging effect can result in poor models, that put most of the probability mass in areas that are not covered by the data. The counterpart of the M-projection is the I(nformation)-projection. The I-projection concentrates on the modes the model is able to represent and ignores the remaining ones. Hence, it does not suffer from the averaging effect (Bishop, 2006).
In this paper, we explore the I-projection for mixture models which are typically trained by maximizing the likelihood via expectation maximization (EM) (Dempster et al., 1977). Despite the richness of mixture models, the averaging problem remains as we typically do not know the correct number
of modes and it is hard to identify all modes of the data correctly. By the use of the I-projection, our mixture models do not suffer from this averaging effect and can generate more realistic samples that are less distinguishable from the data. In this paper we concentrate on learning Gaussian mixture models and conditional Gaussian mixtures of experts (Jacobs et al., 1991) where the mean and covariance matrix are generated by deep neural networks.
We propose Expected Information Maximization (EIM) 1, a novel approach capable of computing the I-projection between the model and the data. By exploiting the structure of the I-projection, we can derive a variational upper bound objective, which was previously used in the context of variational inference (Maaløe et al., 2016; Ranganath et al., 2016; Arenz et al., 2018). In order to work with this upper bound objective based on samples, we use a discriminator to approximate the required density ratio, relating our approach to GANs (Goodfellow et al., 2014; Nowozin et al., 2016; Uehara et al., 2016). The discriminator also allows us to use additional discriminative features to improve model quality. In our experiments, we demonstrate that EIM is much more effective in computing the I-projection than recent GAN approaches. We apply EIM to a synthetic obstacle avoidance task, an inverse kinematic task of a redundant robot arm as well as a pedestrian and car prediction task using the Stanford Drone Dataset (Robicquet et al., 2016) and a traffic dataset from the Next Generation Simulation program.

heading: 2 PRELIMINARIES
text: Our approach heavily builds on minimizing Kullback-Leibler divergences as well as the estimation of density ratios. We will therefore briefly review both concepts.
Density Ratio Estimation. Our approach relies on estimating density ratios r(x) = q(x)/p(x) based on samples of q(x) and p(x). Sugiyama et al. (2012) introduced a framework to estimate such density ratios based on the minimization of Bregman divergences (Bregman, 1967). For our work we employ one approach from this framework, namely density ratio estimation by binary logistic regression. Assume a logistic regressor C(x) = σ(φ(x)) with logits φ(x) and sigmoid activation function σ. Further, we train C(x) on predicting the probability that a given sample x was sampled from q(x). It can be shown that such a logistic regressor using a cross-entropy loss is optimal for C(x) = q(x)/ ( q(x) + p(x) ) . Using this relation, we can compute the log density ratio estimator by
log q(x)
p(x) = log
q(x)/ ( q(x) + p(x) ) p(x)/ ( q(x) + p(x) ) = log C(x) 1− C(x) = σ −1(C(x)) = φ(x). The logistic regressor is trained by minimizing the binary cross-entropy
argminφ(x)BCE(φ(x), p(x), q(x)) = −Eq(x) [log (σ(φ(x)))]− Ep(x) [log (1− σ(φ(x)))] ,
where different regularization techniques such as `2 regularization or dropout (Srivastava et al., 2014) can be used to avoid overfitting.
1Code available at https://github.com/pbecker93/ExpectedInformationMaximization
Moment and Information Projection. The Kullback-Leibler divergence (Kullback & Leibler, 1951) is a standard similarity measure for distributions. It is defined as KL (p(x)||q(x)) =∫ p(x) log p(x)/q(x)dx. Due to its asymmetry, the Kullback-Leibler Divergence provides two different optimization problems (Bishop, 2006) to fit a model distribution q(x) to a target distribution p(x), namely
argminq(x)KL (p(x)||q(x))︸ ︷︷ ︸ Moment-projection and argminq(x)KL (q(x)||p(x))︸ ︷︷ ︸ Information-projection .
Here, we will assume that p(x) is the data distribution, i.e., p(x) is unknown but we have access to samples from p(x). It can easily be seen that computing the M-projection to the data distribution is equivalent to maximizing the likelihood (ML) of the model (Bishop, 2006). ML solutions match the moments of the model with the moments of the target distribution, which results in averaging over modes that can not be represented by the model. In contrast, the I-projection forces the learned generator q(x) to have low probability whenever p(x) has low probability, which is also called zero forcing.

heading: 3 RELATED WORK
text: We will now discuss competing methods for computing the I-projection that are based on GANs. Those are, to the best of our knowledge, the only other approaches capable of computing the Iprojection solely based on samples of the target distribution. Furthermore, we will distinguish our approach from approaches based on variational inference that also use the I-projection.
Variational Inference. The I-projection is a common objective in Variational Inference (Opper & Saad, 2001; Bishop, 2006; Kingma & Welling, 2013). Those methods aim to fit tractable approximations to intractable distributions of which the unnormalized density is available. EIM, on the other hand, does not assume access to the unnormalized density of the target distributions but only to samples. Hence, it is not a variational inference approach, but a density estimation approach. However, our approach uses an upper bound that has been previously applied to variational inference (Maaløe et al., 2016; Ranganath et al., 2016; Arenz et al., 2018). EIM is especially related to the VIPS algorithm (Arenz et al., 2018), which we extend from the variational inference case to the density estimation case. Additionally, we introduce conditional latent variable models into the approach.
Generative Adversarial Networks. While the original GAN approach minimizes the JensenShannon Divergence (Goodfellow et al., 2014), GANs have since been adapted to a variety of other distance measures between distributions, such as the Wasserstein distance (Arjovsky et al., 2017), symmetric KL (Chen et al., 2018) and arbitrary f -divergences (Ali & Silvey, 1966; Nowozin et al., 2016; Uehara et al., 2016; Poole et al., 2016). Since the I-projection is a special case of an f - divergence, those approaches are of particular relevance to our work. Nowozin et al. (2016) use a variational bound for f -divergences (Nguyen et al., 2010) to derive their approach, the f -GAN. Uehara et al. (2016) use a bound that directly follows from the density ratio estimation under Bregman divergences framework introduced by Sugiyama et al. (2012) to obtain their b-GAN. While the b-GAN’s discriminator directly estimates the density ratio, the f -GAN’s discriminator estimates an invertible mapping of the density ratio. Yet, in the case of the I-projection, both the f -GAN and the b-GAN yield the same objective, as we show in Appendix C.2. For both the f -GAN and b-GAN the desired f -divergence determines the discriminator objective. Uehara et al. (2016) note that the discriminator objective, implied by the I-projection, is unstable. As both approaches are formulated in a general way to minimize any f -divergence, they do not exploit the special structure of the Iprojection. Exploiting this structure permits us to apply a tight upper bound of the I-projection for latent variable models, which results in a higher quality of the estimated models.
Li et al. (2019) introduce an adversarial approach to compute the I-projection based on density ratios, estimated by logistic regression. Yet, their approach assumes access to the unnormalized target density, i.e., they are working in a variational inference setting. The most important difference to GANs is that we do not base EIM on an adversarial formulation and no adversarial game has to be solved. This removes a major source of instability in the training process, which we discuss in more detail in Section 4.3.

heading: 4 EXPECTED INFORMATION MAXIMIZATION
text: Expected Information Maximization (EIM) is a general algorithm for minimizing the I-projection for any latent variable model. We first derive EIM for general marginal latent variable models, i.e., q(x) = ∫ q(x|z)q(z)dz and subsequently extend our derivations to conditional latent variable mod-
els, i.e., q(x|y) = ∫ q(x|z, y)q(z|y)dz. EIM uses an upper bound for the objective of the marginal distribution. Similar to Expectation-Maximization (EM), our algorithm iterates between an M-step and an E-step. In the corresponding M-step, we minimize the upper bound and in the E-step we tighten it using a variational distribution.

heading: 4.1 EIM FOR LATENT VARIABLE MODELS
text: The I-projection can be simplified using a (tight) variational upper bound (Arenz et al., 2018) which can be obtained by introducing an auxiliary distribution q̃(z|x) and using Bayes rule
KL (q(x)||p(x)) = Uq̃,p(q)︸ ︷︷ ︸ upper bound −Eq(x)[KL (q(z|x)||q̃(z|x))︸ ︷︷ ︸ ≥0 ],
where
Uq̃,p(q) =
∫∫ q(x|z)q(z) ( log
q(x|z)q(z) p(x)
− log q̃(z|x) ) dzdx. (1)
The derivation of the bound is given in Appendix B. It is easy to see that Uq̃,p(q) is an upper bound as the expected KL term is always non-negative. In the corresponding E-step, the model from the previous iteration, which we denote as qt(x), is used to tighten the bound by setting q̃(z|x) = qt(x|z)qt(z)/qt(x). In the M-step, we update the model distribution by minimizing the upper bound Uq̃,p(q). Yet, opposed to Arenz et al. (2018), we cannot work directly with the upper bound since it still depends on log p(x), which we cannot evaluate. However, we can reformulate the upper bound by setting the given relation for q̃(z|x) of the E-step into Eq. 1,
Uqt,p(q) =
∫ q(z) (∫ q(x|z) log qt(x)
p(x) dx+ KL (q(x|z)||qt(x|z))
) dz + KL (q(z)||qt(z)) . (2)
The upper bound now contains a density ratio between the old model distribution and the data. This density ratio can be estimated using samples of qt and p, for example, by using logistic regression as shown in Section 2. We can use the logits φ(x) of such a logistic regressor to estimate the log density ratio log(qt(x)/p(x)) in Equation 2. This yields an upper bound Uqt,φ(q) that depends on φ(x) instead of p(x). Optimizing this bound corresponds to the M-step of our approach. In the E-step, we set qt to the newly obtained q and retrain the density ratio estimator φ(x). Both steps formally result in the following bilevel optimization problem qt+1 ∈ argminq(x)Uqt,φ∗(q) s.t. φ∗(x) ∈ argminφ(x)BCE(φ(x), p(x), qt(x)). Using a discriminator also comes with the advantage that we can use additional discriminative features g(x) as input to our discriminator that are not directly available for the generator. For example, if x models trajectories of pedestrians, g(x) could indicate whether the trajectory reaches any positions that are not plausible such as rooftops or trees. These features simplify the discrimination task and can therefore improve our model accuracy which is not possible with M-projection based algorithms such as EM.

heading: 4.2 EIM FOR CONDITIONAL LATENT VARIABLE MODELS
text: For conditional distributions, we aim at finding the conditional I-projection argminq(x|y)Ep(y) [KL (q(x|y)||p(x|y))] .
The derivations for the conditional upper bound follow the same steps as the derivations in the marginal case, where all distributions are extended by the context variable y. We refer to the supplement for details. The log density ratio estimator φ(x, y) now discriminates between samples of the joint distribution of x and y. For training φ(x, y) we generate a new sample x for each context y, using the distribution qold(x|y). Hence, as the context distribution is the same for the true data and the generated data, the log density ratio of the conditional distributions is equal to the log density ratio of the joint distributions.

heading: 4.3 RELATION TO GANS AND EM
text: There is a close relation of EIM to GANs due to the use of a logistic discriminator for the density ratio estimation. It is therefore informative to investigate the differences in the case without latent variables. In an adversarial formulation, the density ratio estimator would directly replace the density ratio in the original I-projection equation, i.e.,
argminq(x) ∫ q(x)φ∗(x)dx s.t. φ∗(x) ∈ argminφ(x)BCE(φ(x), p(x), q∗(x)).
However, such adversarial games are often hard to optimize. In contrast, EIM offers a bilevel optimization problem where the discriminator is explicitly learned on the old data distribution qt(x),
argminq(x) ∫ q(x)φ∗(x)dx+ KL (q(x)||qt(x)) s.t. φ∗(x) ∈ argminφ(x)BCE(φ(x), p(x), qt(x)).
Thus, there is no circular dependency between the optimal generator and the optimal discriminator. Figure 2 illustrates that the proposed non-adversarial formulation does not suffer from too large model updates. Choosing the number and step-size of the updates is thus far less critical.
EIM can also be seen as the counter-part of Expectation-Maximization (EM). While EM optimizes the M-projection with latent variable models, EIM uses the I-projection. However, both approaches decompose the corresponding projections into an upper bound (or lower bound for EM) and a KLterm that depends on the conditional distribution q(z|x) to tighten this bound. The exact relationship is discussed in Appendix C.1.

heading: 4.4 EIM FOR GAUSSIAN MIXTURES MODELS
text: We consider Gaussian mixture models with d components, i.e., multivariate Gaussian distributions q(x|zi) = N (µi,Σi) and a categorical distribution q(z) = Cat(π) for the coefficients. As the latent distribution q(z) is discrete, the upper bound in EIM (Equation 2) simplifies, as the integral over z can be written as a sum. Similar to the EM-algorithm, this objective can be updated individually for the coefficients and the components. For both updates, we will use similar update rules as defined in the VIPS algorithm (Arenz et al., 2018). VIPS uses a trust region optimization for the components and the coefficients, where both updates can be solved in closed form as the components are Gaussian. The trust regions prevent the new model from going too far away from qt where the density ratio estimator is inaccurate, and hence, further stabilize the learning process. We will now sketch both updates, where we refer to Appendix B.2 for the full details.
For updating the coefficients, we assume that the components have not yet been updated, and therefore KL (q(x|zi)||qt(x|zi)) = 0 for all zi. The objective for the coefficients thus simplifies to
argminq(z) d∑ i=1 q(zi)φ(zi) + KL (q(z)||qt(z)) with φ(zi) = Eq(x|zi) [φ(x)] , (3)
where φ(zi) can be approximated using samples from the corresponding component. This objective can easily be optimized in closed form, as shown in the VIPS algorithm (Arenz et al., 2018). We also use a KL trust-region to specify the step size of the update. For updating the individual components, the objective simplifies to
argminq(x|zi)Eq(x|zi) [φ(x)] + KL (q(x|zi)||qt(x|zi)) . (4) As in VIPS, this optimization problem can be solved in closed form using the MORE algorithm (Abdolmaleki et al., 2015). The MORE algorithm uses a quadratic surrogate function that locally approximates φ(x). The resulting solution optimizes Equation 4 under a KL trust-region. The pseudo-code of EIM for GMMs can be found in Appendix A.

heading: 4.5 EIM FOR GAUSSIAN MIXTURES OF EXPERTS
text: In the conditional case, we consider mixtures of experts consisting of d multivariate Gaussians, whose parameters depend on an input y in a nonlinear fashion, i.e., q(x|zi,y) = N (ψµ,i(y), ψΣ,i(y)) and the gating is given by a neural network with softmax output. We again decompose the resulting upper bound into individual update steps for the components and the gating. Yet, closed-form solutions are no longer available and we need to resort to gradient-based updates. The objective for updating the gating is given by
argminq(z|y) d∑ i=1 ( Ep(y)q(x|zi,y) [q(zi|y)φ(x,y)] ) + Ep(y) [KL (q(z|y)||qt(z|y))] . (5)
We minimize this equation w.r.t. the parameters of the gating by gradient descent using the Adam (Kingma & Ba, 2014) algorithm. The objective for updating a single component i is given by
argminq(x|zi,y)Ep̃(y|zi) [ Eq(x|zi,y) [φ(x,y)] + KL (q(x|zi,y)||qt(x|zi,y)) ] , (6)
where p̃(y|zi) = p(y)q(zi|y)/q(zi). Note that we normalized the objective by q(zi) =∫ p(y)q(zi|y)dy to ensure that also components with a low prior q(zi) get large enough gradients for the updates. As we have access to the derivatives of the density ratio estimator w.r.t. x, we can optimize Equation 6 with gradient descent using the reparametrization trick (Kingma & Welling, 2013) and Adam.

heading: 5 EVALUATION
text: We compare our approach to GANs and perform an ablation study on a toy task, with data sampled from known mixture models. We further apply our approach to two synthetic datasets, learning the joint configurations of a planar robot as well as a non-linear obstacle avoidance task, and two real datasets, namely the Stanford Drone Dataset (Robicquet et al., 2016) and a traffic dataset from the Next Generation Simulation program. A full overview of all hyperparameters and network architectures can be found in Appendix E.

heading: 5.1 COMPARISON TO GENERATIVE ADVERSARIAL APPROACHES AND ABLATION STUDY
text: We compare to the f -GAN which is the only other method capable of minimizing the I-projection solely based on samples. We use data sampled from randomly generated GMMs with different numbers of components and dimensionalities. To study the influence of the previously mentioned differences of EIM to generative adversarial approaches, we also perform an ablation study. We compare to a version of EIM where we neglect the additional KL-term (EIM, no KL), a version were we trained all components and the coefficients jointly using gradient descent (Joint EIM), and a version where we do both (Joint EIM, no KL). The average I-projection achieved by the various approaches can be found in Figure 3.

heading: 5.2 LINE REACHING WITH PLANAR ROBOT
text: We extended the introductory example of the planar reaching task and collected expert data from a 10-link planar robot tasked with reaching a point on a line. We fitted GMMs with an increasing
number of components using EIM, EIM with additional features, where the end-effector coordinates for a given joint configuration were provided, and EM. Even for a large number of components, we see effects similar to the introductory example, i.e., the M-projection solution provided by EM fails to reach the line while EIM manages to do so. For small numbers, EIM ignores parts of the line, while more and more parts of it get covered as we increase the number of components. With the additional features, the imitation of the line reaching was even more accurate. See Figure 4, for the average distance between the end-effector and the line as well as samples from both EM and EIM.

heading: 5.3 PEDESTRIAN AND TRAFFIC PREDICTION
text: We evaluated our approach on data from the Stanford Drone Dataset (SDD) (Robicquet et al., 2016) and a traffic dataset from the Next Generation Simulation (NGS) program2. The SDD data consists of trajectories of pedestrians, bikes, and cars and we only used the data corresponding to a single video of a single scene (Video 1, deathCircle). The NGS data consists of trajectories of cars where we considered the data recorded on Lankershim Boulevard. In both cases we extracted trajectories of length 5, yielding highly multimodal data due to pedestrians, bikes, and cars moving at different speeds and in different directions. We evaluated on the achieved log-likelihood of EIM and EM, see Figure 5. EM achieves the highest likelihood as it directly optimizes this measure. However, we can already see that EM massively overfits when we increase the number of components as the test-set likelihood degrades. EIM, on the other hand, produced better models with an increasing number of components. Additionally, we generated a mask indicating whether a given point is on the road or
2 https://data.transportation.gov/Automobiles/Next-Generation-Simulation-NGSIM-Vehicle-Trajector/8ect-6jqj
not and evaluated how realistic the learned models are by measuring the amount of samples violating the mask, i.e., predicting road users outside of the road. We also evaluate a version of EIM where we provide additional features indicating if the mask is violated. EIM achieves a much better value on this mask for the NGS dataset, while we needed the additional mask features for the discriminator on the SDD dataset to outperform EM. Both experiments show that EIM can learn highly multi-modal density estimates that produce more realistic samples than EM. They further show that the models learned by EIM can be refined by additional prior knowledge provided as feature vectors.

heading: 5.4 OBSTACLE AVOIDANCE
text: We evaluate the conditional version of EIM on an artificial obstacle avoidance task. The context contains the location of three obstacles within an image. The gating, as well as the components, are given by deep neural networks. Details about the network architectures can be found in the Appendix. The data consists of trajectories going from the left to the right of the image. The trajectories are defined by setting 3 via-points such that no obstacle is hit. To generate the data we sample via-points over and under the obstacles with a probability proportional to the distance between the obstacle and the image border. Hence, for three obstacles, there are 23 = 8 different modes in the data. Note that, like in most real-world scenarios, the expert data is not perfect, and about 13% of the trajectories in the dataset hit an obstacle. We fit models with various numbers of components to this data using EIM and EM and compare their performance on generating trajectories that achieve the goal. Results are shown in Figure 6 together with a visualization of the task and samples produced by EIM and EM. EIM was able to identify most modes for the different given inputs and did not suffer from any averaging effect. In contrast, EM does not find all modes. As a consequence, some
components of the mixture model had to average over multiple modes, resulting in poor quality trajectories.

heading: 6 CONCLUSION
text: We introduced Expected Information Maximization (EIM), a novel approach for computing the Iprojection between general latent variable models and a target distribution, solely based on samples of the latter. General upper bound objectives for marginal and conditional distributions were derived, resulting in an algorithm similar to EM, but tailored for the I-projection instead of the M-projection. We introduced efficient methods to optimize these upper bound objectives for mixture models. In our experiments, we demonstrated the benefits of the I-projection for different behavior modelling tasks. The introduced approach opens various pathways for future research. While we focused on mixture models, the derived upper bounds are not exclusive to those and can be used for arbitrary latent variable models. Another possibility is an online adaptation of the number of used components. Arenz et al. (2018) propose heuristics for such an adaptation in their VIPS approach. Those could easily be adapted to our approach.

heading: A PSEUDO CODE
text: EIM-for-GMMs({x(j)p }j=1···N , q(x)); Input: Data {x(j)p }j=1···N , Initial Model q(x) = ∑d i=1 q(x|zi)q(zi) = ∑d i=1 πiN (x|µi,Σi) for i in number of iterations do E-Step: qt(z) = q(z), qt(x|zi) = q(x|zi) for all components i Update Density Ratio Estimator: sample data from model {x(j)q }j=1···N ∼ qt(x) retrain density ratio estimator φ(x) on {x(j)p }j=1···N and {x(j)q }j=1···N M-Step Coefficients: for i in number of components do
compute loss li = 1
N
∑N j=1 φ ( x (j) q ) with samples {x(j)q }j=1···N ∼ qt(x|zi)
end update q(z) using losses li and MORE equations M-Step Components: for i in number of components do
fit φ̂(x) surrogate to pairs ( x (j) q , φ ( x (j) q )) with samples {x(j)q }j=1···N ∼ qt(x|zi)
update q(x|zi) using surrogate φ̂(x) and MORE equations end
end Algorithm 1: Expected Information Maximization for Gaussian Mixture Models.
Pseudo-code for EIM for GMMs can be found in algorithm 1

heading: B DERIVATIONS
text: Derivations of the upper bound stated in Equation 1. We assume latent variable models q(x) =∫ q(x|z)q(z)dz and use the identities q(x|z)q(z) = q(z|x)q(x) and log q(x) = log q(x|z)q(z) − log q(z|x).
KL (q(x)||p(x)) = ∫ q(x) log q(x)
p(x) dx =
∫∫ q(x|z)q(z) log q(x)
p(x) dzdx
= ∫∫ q(x|z)q(z) ( log
q(x|z)q(z) p(x)
− log q(z|x) ) dzdx
= ∫∫ q(x|z)q(z) ( log
q(x|z)q(z) p(x)
− log q(z|x) + log q̃(z|x)− log q̃(z|x) ) dzdx
= ∫∫ q(x|z)q(z) ( log
q(x|z)q(z) p(x)
− log q̃(z|x) ) dzdx− ∫∫ q(x|z)q(z) (log q(z|x)− log q̃(z|x)) dzdx
= ∫∫ q(x|z)q(z) ( log
q(x|z)q(z) p(x)
− log q̃(z|x) ) dzdx− ∫ q(x) ∫ q(z|x) log q(z|x)
q̃(z|x)dzdx
=U(q, q̃, p)− Eq(x) [KL (q(z|x)||q̃(z|x))] .
After plugging the E-Step, i.e., q̃(z|x) = qt(x|z)qt(z)/qt(x), into the objective it simplifies to
U(q, q̃, p)
= ∫∫ q(x|z)q(z) ( log
q(x|z)q(z) p(x) − log qt(x|z)qt(z) qt(x)
) dzdx
= ∫∫ q(x|z)q(z) (log q(x|z) + log q(z)− log p(x)− log qt(x|z)− log qt(z) + log qt(x)) dzdx
= ∫∫ q(x|z)q(z) ( log qt(x)
p(x) + log q(x|z) qt(x|z) + log q(z) qt(z)
) dzdx
= ∫ q(z) (∫ q(x|z) ( log qt(x)
p(x) + log q(x|z) qt(x|z)
) dx+ log q(z)
qt(z)
) dz
= ∫ q(z) ∫ q(x|z) log qt(x)
p(x) dxdz +
∫ q(z) ∫ q(x|z) log q(x|z)
qt(x|z) dxdz +
∫ q(z) log q(z)
qt(z) dz
= ∫∫ q(x|z)q(z) log qt(x)
p(x) dzdx+ Eq(z) [KL (q(x|z)||qt(x|z))] + KL (q(z)||qt(z)) ,
which concludes the derivation of upper bound of latent variable models.
B.1 DERIVATIONS CONDITIONAL UPPER BOUND
By introducing an auxiliary distribution q̃(z|x, y) the upper bound to the expected KL for conditional latent variable models q(x|y) = ∫ q(x|z, y)q(z|y)dz can be derived by
Ep(y) [KL (q(x|y)||p(x|y))] = ∫∫
p(y)q(x|y) log q(x|y) p(x|y)dxdy
= ∫ p(y) ∫∫ q(x|z, y)q(z|y) ( log
q(x|z, y)q(z|y) p(x|y) − log q(z|x, y)
) dzdxdy
= ∫ p(y) ∫∫ q(x|z, y)q(z|y)
· ( log
q(x|z, y)q(z|y) p(x|y) − log q(z|x, y) + log q̃(z|x, y)− log q̃(z|x, y)
) dzdxdy
= ∫ p(y) ∫∫ q(x|z, y)q(z|y) ( log
q(x|z, y)q(z|y) p(x|y) − log q̃(z|x, y)
) dzdxdy
− ∫ p(y) ∫∫ q(x|z, y)q(z, y) (log q(z|x, y)− log q̃(z|x, y)) dzdxdy
= ∫ p(y) ∫∫ q(x|z, y)q(z|y) ( log
q(x|z, y)q(z|y) p(x|y) − log q̃(z|x, y)
) dzdxdy
− ∫∫ p(y)q(x|y) ∫ q(z|x, y) log q(z|x, y)
q̃(z|x, y)dzdxdy
=U(q, q̃, p)− Ep(y),q(x|y) [KL (q(z|x, y)||q̃(z|x, y))] .
During the E-step the bound is tightened by setting q̃(z|x, y) = qt(x|z, y)qt(z|y)/qt(x|y). U(q, q̃, p)
= ∫ p(y) ∫∫ q(x|z, y)q(z|y) ( log
q(x|z, y)q(z|y) p(x|y) − log qt(x|z, y)qt(z|y) qt(x|y)
) dzdxdy
= ∫ p(y) ∫∫ q(x|z, y)q(z|y)
· (log q(x|z, y) + log q(z|y)− log p(x|y)− log qt(x|z, y)− log qt(z|y) + log qt(x|y)) dzdxdy
= ∫ p(y) ∫∫ q(x|z, y)q(z|y) ( log
qt(x|y) p(x|y) + log q(x|z, y) qt(x|z, y) + log q(z|y) qt(z|y)
) dzdxdy
= ∫ p(y) ∫ q(z|y) (∫ q(x|z, y) ( log
qt(x|y) p(x|y) + log q(x|z, y) qt(x|z, y)
) dx+ log
q(z|y) qt(z|y)
) dzdy
= ∫ p(y) ∫ q(z|y) ∫ q(x|z, y) log qt(x|y)
p(x|y) dxdzdy
+ ∫ p(y) ∫ q(z|y) ∫ q(x|z, y) log q(x|z, y)
qt(x|z, y) dxdzdy +
∫ p(y) ∫ q(z|y) log q(z|y)
qt(z|y) dzdy
= ∫∫∫ p(y)q(z|y)q(x|z, y) log qt(x|y)
p(x|y) dxdzdy
+ Ep(y),q(z|y) [KL (q(x|z, y)||qt(x|z, y))] + Ep(y) [KL (q(z|y)||qt(z|y))] , which concludes the derivation of the upper bound for conditional latent variable models.
B.2 USING MORE FOR CLOSED FORM UPDATES FOR GMMS
The MORE algorithm, as introduced by Abdolmaleki et al. (2015), can be used to solve optimization problems of the following form
argmaxq(x)Eq(x)[f(x)] s.t. KL (q(x)||qold(x)) ≤ for an exponential family distribution q(x), some function f(x), and an upper bound on the allowed change, . Abdolmaleki et al. (2015) show that the optimal solution is given by
q(x) ∝ qold(x) exp ( f(x)
η
) = exp ( η log qold(x) + f(x)
η
) ,
where η denotes the Lagrangian multiplier corresponding to the KL constraint. In order to obtain this Lagrangian multiplier, the following, convex, dual function has to be minimized
g(η) = η + η log ∫ exp ( η log qold(x) + f(x)
η
) dx. (7)
For discrete distributions, such as the categorical distribution used to represent the coefficients of a GMM, we can directly work with those equations. For continuous distributions, Abdolmaleki et al. (2015) propose approximating f(x) with a local surrogate. The features to fit this surrogate are chosen such that they are compatible (Kakade, 2002), i.e., of the same form as the distributions sufficient statistics. For multivariate Gaussians, the sufficient statistics are squared features and thus the surrogate compatible to such a Gaussian distribution is given by
f̂(x) = −1 2 xT F̂x + f̂ T x + f0.
The parameters of this surrogate can now be used to update the natural parameters of the Gaussian, i.e, the precision matrix Q = Σ−1 and q = Σ−1µ by
Q = Qt + 1
η F̂ and q = qt +
1 η f̂ .
In order to apply the MORE algorithm to solve the optimization problems stated in Equation 3 and Equation 4 we make two trivial modifications. First, we invert the signs in Equation 3 and Equation 4, as we are now maximizing. Second, to account for the additional KL term in our objectives, we add 1 to η, everywhere except the first term of the sum in Equation 7.

heading: C ELABORATION ON RELATED WORK
text: C.1 RELATION BETWEEN EIM AND EM
Recall that the Expectation-Maximization (EM) algorithm (Dempster et al., 1977) maximizes the log-likelihood of the data by iteratively maximizing and tightening the following lower bound
Ep(x) [log q(x)] = Ep(x) [∫
q̃(z|x) log q(x, z) q̃(z|x) dz
] + Ep(x) [∫ q̃(z|x) log q̃(z|x) q(z|x)dz ]
= L(q, q̃)︸ ︷︷ ︸ lower bound +Ep(x) [KL (q̃(z|x)||q(z|x))]︸ ︷︷ ︸ ≥0 .
It is instructive to compare our upper bound (Equation 1) to this lower bound. As mentioned, maximizing the likelihood is equivalent to minimizing the M-projection, i.e., argminq(x)KL (p(x)||q(x)), where, in relation to our objective, the model and true distribution have switched places in the nonsymmetric KL objective. Like our approach, EM introduces an auxiliary distribution q̃(z|x) and bounds the objective from below by subtracting the KL between auxiliary distribution and model, i.e., KL (q̃(z|x)||q(z|x)). In contrast, we obtain our upper bound by adding KL (q(z|x)||q̃(z|x)) to the objective. Again, the distributions have exchanged places within the KL.
C.2 EQUALITY OF f -GAN AND b-GAN
As pointed out in section 3 both the f -GAN (Nowozin et al., 2016) and the b-GAN (Uehara et al., 2016) yield the same objective for the I-projection.
We start with the f -GAN. Nowozin et al. (2016) propose the following adversarial objective, based on a variatonal bound for f -divergences (Nguyen et al., 2010)
argminq(x)argmaxV (x)F (q(x), V (x)) = Ep(x) [g(V (x))]− Eq(x) [f∗(g(V (x)))] . Here V (x) denotes a neural network with linear output, g(v) the output activation and f∗(t) the Fenchel conjugate (Hiriart-Urruty & Lemaréchal, 2012) of f(u), i.e., the generator function of the f -divergence. For the I-projection f(u) = − log u and f∗(t) = −1 − log(t). In theory, the only restriction posed on the choice of g(v) is that it outputs only values within the domain of f∗(t) Nowozin et al. (2016) suggest g’s for various f− divergences and chose exclusively monotony increasing functions which output large values for samples that are believed to be from the data distribution. For the I-projection they suggest g(v) = −exp(−v). Thus the f -GAN objective for the I-projection is given by
argminq(x)argmaxV (x)F (q(x), V (x)) = −Ep(x) [exp(−V (x)] + Eq(x) [1− V (x)] .
The b-GAN objective follows from the density ratio estimation framework given by Sugiyama et al. (2012) and is given by
argminq(x)argmaxr(x)Ep(x) [f ′(r(x))]− Eq(x) [f ′(r(x))r(x)− f(r(x))]
Here f ′(u) denotes the derivative of f(u) and r(x) denotes an density ratio estimator. We need to enforce r(x) > 0 for all x to obtain a valid density ratio estimate. In practice this is usually done by learning rl(x) = log r(x) instead. Plugging rl(x), f(u) and f ′(u) = 1/u into the general b-GAN objective yields
argminq(x)argmaxrl(x)F (q(x), rl(x)) = −Ep(x) [exp(−rl(x))] + Eq(x) [1− rl(x)] . Which is the same objective as the f -GAN uses. Yet, f -GAN and b-GAN objectives are not identical for arbitrary f -divergences.
D VISUALIZATION OF SAMPLES
D.1 PEDESTRIAN AND TRAFFIC PREDICTION
Samples from the Stanford Drone Dataset can be found in Figure 7

heading: E HYPERPARAMETERS
text: In all experiments, we realize the density ratio estimator as fully connected neural networks which we train using Adam (Kingma & Ba, 2014) and early stopping using a validation set.
Comparison to Generative Adversarial Approaches and Ablation Study
• Data: 10, 000 Train Samples, 5, 000 Test Samples, 5, 000 Validation samples (for early stopping the density ratio estimator)
• Density Ratio Estimator (EIM) / Variational function V (x) (f -GAN): 3 fully connected layers, 50 neurons each, trained with L2 regularization with factor 0.001, early stopping and batch size 1, 000
• Updates EIM: MORE-like updates with = 0.05 for components and coefficients, 1, 000 samples per component and update
• Updates FGAN: Iterate single update steps for generator and discriminator using learning rates of 1e− 3 and batch size of 1, 000.
Line Reaching with Planar Robot
• Data: 10, 000 train samples, 5, 000 test samples, 5, 000 validation samples (for early stopping the density ratio estimator)
• Density Ratio Estimation: 2 fully connected layers of width 100, early stopping and batch size 1, 000
• Updates: MORE-like updates with = 0.005 for components and coefficients, 1, 000 samples per component and update
Pedestrian and Traffic Prediction
• Data SDD: 7, 500 train samples, 3, 801 test samples, 3, 801 validation samples (for early stopping the density ratio estimator)
• Data NGS: 10, 000 train samples, 5, 000 test samples, 5, 000 validation samples (for early stopping the density ratio estimator)
• Density Ratio Estimation: 3 fully connected layers of width 256, trained with L2 regularization with factor 0.0005 early stopping and batch size 1, 000.
• Updates: MORE-like updates with = 0.01 for components and coefficients, 1, 000 samples per component and update
Obstacle Avoidance
• Data: 1, 000 train contexts with 10 samples each, 500 test contexts with 10 samples each, 500 validation contexts with 10 samples each (for early stopping the density ratio estimator).
• Density Ratio Estimation: 3 fully connected layers of width 256, trained with L2 regularization with factor 0.0005 early stopping and batch size 1, 000.
• Component and Gating networks: 2 fully connected layer of width 64 for each component and the gating. Trained with Adam (α = 1e− 3, β0 = 0.5) for 10 epochs in each iteration.

Review summaries:
The paper presents an algorithm to match two distributions with latent variables, named expected information maximization (EIM). Specifically, EIM is based on the I-Projection, which basically is equivalent to minimizing the reverse KL divergence (i.e. min KL[p_model || p_data]); to handle latent variables, an upper-bound is derived, which is the corresponding reverse KL divergence in the joint space. To minimize that joint reverse KL, a specific procedure is developed, leading to the presented EIM. EIM variants for different applications are discussed. Fancy robot-related experiments are used to evaluate the presented algorithm.

This paper propose EIM an analog to EM but to perform the I-projection (i.e. reverse-KL) instead of the usual M-projection for EM. The motivation is that the reverse-KL is mode-seeking in contrast to the forward-KL which is mode-covering. The authors argue that in the case that the model is mis-specified, I-projection is sometimes desired as to avoid putting mass on very unlikely regions of the space under the target p.