Efficient Gradient-Based Inference through Transformations between Bayes Nets and Neural Nets

Introduction

Bayesian networks (also called belief networks) are probabilistic graphical models where the conditional dependencies within a set of random variables are described by a directed acyclic graph (DAG). Many supervised and unsupervised models can be considered as special cases of Bayesian networks.

In this paper we focus on the problem of efficient inference in Bayesian networks with multiple layers of continuous latent variables, where exact posterior inference is intractable (e.g. the conditional dependencies between variables are nonlinear) but the joint distribution is differentiable. Algorithms for approximate inference in Bayesian networks can be roughly divided into two categories: sampling approaches and parametric approaches. Parametric approaches include Belief Propagation (Pearl, 1982) or the more recent Expectation Propagation (EP) (Minka, 2001). When it is not reasonable or possible to make assumptions about the posterior (which is often the case), one needs to resort to sampling approaches such as Markov Chain Monte Carlo (MCMC) (Neal, 1993). In high-dimensional spaces, gradient-based samplers such as Hybrid Monte Carlo (Duane et al., 1987) and the recently proposed no-U-turn sampler (Hoffman & Gelman, 2011) are known for their relatively fast mixing properties. When just interested in finding a mode of the posterior, vanilla gradient-based optimization methods can be used. The alternative parameterizations suggested in this paper can dramatically improve the efficiency of any of these algorithms.

After reviewing background material in 2, we introduce a generally applicable differentiable reparameterization of continuous latent variables into a differentiable non-centered form in section 3. In section 4 we analyze the posterior dependencies in this reparameterized form. Experimental results are shown in section 6.

Background

We use bold lower case (e.g. $\mathbf{x}$ or $\mathbf{y}$) notation for random variables and instantiations (values) of random variables. We write $p_{\boldsymbol{\theta}}(\mathbf{x}|\mathbf{y})$ and $p_{\boldsymbol{\theta}}(\mathbf{x})$ to denote (conditional) probability density (PDF) or mass (PMF) functions of variables. With $\boldsymbol{\theta}$ we denote the vector containing all parameters; each distribution in the network uses a subset of $\boldsymbol{\theta}$’s elements. Sets of variables are capitalized and bold, matrices are capitalized and bold, and vectors are written in bold and lower case.

1 Bayesian networks

A Bayesian network models a set of random variables $\mathbf{V}$ and their conditional dependencies as a directed acyclic graph, where each variable corresponds to a vertex and each edge to a conditional dependency. Let the distribution of each variable $\mathbf{v}_{j}$ be $p_{\boldsymbol{\theta}}(\mathbf{v}_{j}|\mathbf{pa}_{j})$, where we condition on $\mathbf{v}_{j}$’s (possibly empty) set of parents $\mathbf{pa}_{j}$. Given the factorization property of Bayesian networks, the joint distribution over all variables is simply:

Let the graph consist of one or more (discrete or continuous) observed variables $\mathbf{x}_{j}$ and continuous latent variables $\mathbf{z}_{j}$, with corresponding conditional distributions $p_{\boldsymbol{\theta}}(\mathbf{x}_{j}|\mathbf{pa}_{j})$ and $p_{\boldsymbol{\theta}}(\mathbf{z}_{j}|\mathbf{pa}_{j})$. We focus on the case where both the marginal likelihood $p_{\boldsymbol{\theta}}(\mathbf{x})=\int_{\mathbf{z}}p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z})\,d\mathbf{z}$ and the posterior $p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x})$ are intractable to compute or differentiate directly w.r.t. $\boldsymbol{\theta}$ (which is true in general), and where the joint distribution $p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z})$ is at least once differentiable, so it is still possible to efficiently compute first-order partial derivatives $\nabla_{\boldsymbol{\theta}}\log p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z})$ and $\nabla_{\mathbf{z}}\log p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z})$.

2 Conditionally deterministic variables

A conditionally deterministic variable $\mathbf{v}_{j}$ with parents $\mathbf{pa}_{j}$ is a variable whose value is a (possibly nonlinear) deterministic function $g_{j}(.)$ of the parents and the parameters: $\mathbf{v}_{j}=g_{j}(\mathbf{pa}_{j},\boldsymbol{\theta})$. The PDF of a conditionally deterministic variable is a Dirac delta function, which we define as a Gaussian PDF $\mathcal{N}(.;\mu,\sigma)$ with infinitesimal $\sigma$:

which equals $+\infty$ when $\mathbf{v}_{j}=g_{j}(\mathbf{pa}_{j},\boldsymbol{\theta})$ and equals 0 everywhere else such that $\int_{\mathbf{v}_{j}}p_{\boldsymbol{\theta}}(\mathbf{v}_{j}|\mathbf{pa}_{j})\,d\mathbf{v}_{j}=1$.

3 Inference problem under consideration

We are often interested in performing posterior inference, which most frequently consists of either optimization (finding a mode $\operatorname*{argmax}_{\mathbf{z}}p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x})$) or sampling from the posterior $p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x})$. Gradients of the log-posterior w.r.t. the latent variables can be easily acquired using the equality:

In words, the gradient of the log-posterior w.r.t. the latent variables is simply the sum of gradients of individual factors w.r.t. the latent variables. These gradients can then be followed to a mode if one is interested in finding a MAP solution. If one is interested in sampling from the posterior then the gradients can be plugged into a gradient-based sampler such as Hybrid Monte Carlo (Duane et al., 1987); if also interested in learning parameters, the resulting samples can be used for the E-step in Monte Carlo EM (Wei & Tanner, 1990) (MCEM).

Problems arise when strong posterior dependencies exist between latent variables. From eq. (3) we can see that the Hessian $\mathbf{H}$ of the posterior is:

Suppose a factor $\log p_{\boldsymbol{\theta}}(z_{i}|z_{j})$ connecting two scalar latent variables $\mathbf{z}_{i}$ and $\mathbf{z}_{j}$ exists, and $z_{i}$ is strongly dependent on $z_{j}$, then the Hessian’s corresponding element $\frac{\partial^{2}\log p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x})}{\partial z_{i}\partial z_{j}}$ will have a large (positive or negative) value. This is bad for gradient-based inference since it means that changes in $z_{j}$ have a large effect on the gradient $\frac{\partial\log p_{\boldsymbol{\theta}}(z_{i}|z_{j})}{\partial z_{i}}$ and changes in $z_{i}$ have a large effect on the gradient $\frac{\partial\log p_{\boldsymbol{\theta}}(z_{i}|z_{j})}{\partial z_{j}}$. In general, strong conditional dependencies lead to ill-conditioning of the posterior, resulting in smaller optimal stepsizes for first-order gradient-based optimization or sampling methods, making inference less efficient.

The differentiable non-centered parameterization (DNCP)

In this section we introduce a generally applicable transformation between continuous latent random variables and deterministic units with auxiliary parent variables. In rest of the paper we analyze its ramifications for gradient-based inference.

Let $\mathbf{z}_{j}$ be some continuous latent variable with parents $\mathbf{pa}_{j}$, and corresponding conditional PDF:

This is also known in the statistics literature as the centered parameterization (CP) of the latent variable $\mathbf{z}_{j}$. Let the differentiable non-centered parameterization (DNCP) of the latent variable $\mathbf{z}_{j}$ be:

where $g_{j}(.)$ is some differentiable function. Note that in the DNCP, the value of $\mathbf{z}_{j}$ is deterministic given both $\mathbf{pa}_{j}$ and the newly introduced auxiliary variable $\boldsymbol{\epsilon}_{j}$ which is distributed as $p(\boldsymbol{\epsilon}_{j})$. See figure 1 for an illustration of the two parameterizations.

By the change of variables, the relationship between the original PDF $p_{\boldsymbol{\theta}}(\mathbf{z}_{j}|\mathbf{pa}_{j})$, the function $g_{j}(\mathbf{pa}_{j},\boldsymbol{\epsilon}_{j})$ and the PDF $p(\boldsymbol{\epsilon}_{j})$ is:

where $det(\mathbf{J})$ is the determinant of Jacobian of $g_{j}(.)$ w.r.t. $\boldsymbol{\epsilon}_{j}$. If $z_{j}$ is a scalar variable, then $\epsilon_{j}$ is also scalar and $|det(\mathbf{J})|=|\frac{\partial z_{j}}{\partial\epsilon_{j}}|$.

In the DNCP, the original latent variable $\mathbf{z}_{j}$ has become deterministic, and its PDF $p_{\boldsymbol{\theta}}(\mathbf{z}_{j}|\mathbf{pa}_{j},\boldsymbol{\epsilon}_{j})$ can be described as a Dirac delta function (see section 2.2).

The joint PDF over the random and deterministic variables can be integrated w.r.t. the determinstic variables. If for simplicity we assume that observed variables are always leaf nodes of the network, and that all latent variables are reparameterized such that the only random variables left are the observed and auxiliary variables $\mathbf{x}$ and $\boldsymbol{\epsilon}$, then the marginal joint $p_{\boldsymbol{\theta}}(\mathbf{x},\boldsymbol{\epsilon})$ is obtained as follows:

In the last step of eq. (LABEL:eq:integrating_out_z), the inputs $\mathbf{pa}_{j}$ to the factors of observed variables $p_{\boldsymbol{\theta}}(\mathbf{x}_{j}|\mathbf{pa}_{j})$ are defined in terms of functions $\mathbf{z}_{k}=g_{k}(.)$, whose values are all recursively computed from auxiliary variables $\boldsymbol{\epsilon}$.

2 Approaches to DNCPs

There are a few basic approaches to transforming CP of a latent variable $\mathbf{z}_{j}$ to a DNCP:

Tractable and differentiable inverse CDF. In this case, let $\epsilon_{j}\sim\mathcal{U}(0,1)$, and let $g_{j}(\mathbf{z}_{j},\mathbf{pa}_{j},\boldsymbol{\theta})=F^{-1}(\mathbf{z}_{j}|\mathbf{pa}_{j};\boldsymbol{\theta})$ be the inverse CDF of the conditional distribution. Examples: Exponential, Cauchy, Logistic, Rayleigh, Pareto, Weibull, Reciprocal, Gompertz, Gumbel and Erlang distributions.

For any ”location-scale” family of distributions (with differentiable log-PDF) we can choose the standard distribution (with $\text{location}=0$, $\text{scale}=1$) as the auxiliary variable $\boldsymbol{\epsilon}_{j}$, and let $g_{j}(.)=\text{location}+\text{scale}\cdot\boldsymbol{\epsilon}_{j}$. Examples: Gaussian, Uniform, Laplace, Elliptical, Student’s t, Logistic and Triangular distributions.

Composition: It is often possible to express variables as functions of component variables with different distributions. Examples: Log-Normal (exponentiation of normally distributed variable), Gamma (a sum over exponentially distributed variables), Beta distribution, Chi-Squared, F distribution and Dirichlet distributions.

When the distribution is not in the families above, accurate differentiable approximations to the inverse CDF can be constructed, e.g. based on polynomials, with time complexity comparable to the CP (see e.g. (Devroye, 1986) for some methods).

For the exact approaches above, the CP and DNCP forms have equal time complexities. In practice, the difference in CPU time depends on the relative complexity of computing derivatives of $\log p_{\boldsymbol{\theta}}(\mathbf{z}_{j}|\mathbf{pa}_{j})$ versus computing $g_{j}(.)$ and derivatives of $\log p(\epsilon_{j})$, which can be easily verified to be similar in most cases below. Iterations with the DNCP form were slightly faster in our experiments.

3 DNCP and neural networks

It is instructive to interpret the DNCP form of latent variables as ”hidden units” of a neural network. The network of hidden units together form a neural network with inserted noise $\boldsymbol{\epsilon}$, which we can differentiate efficiently using the backpropagation algorithm (Rumelhart et al., 1986).

There has been recent increase in popularity of deep neural networks with stochastic hidden units (e.g. (Krizhevsky et al., 2012; Goodfellow et al., 2013; Bengio, 2013)). Often, the parameters $\boldsymbol{\theta}$ of such neural networks are optimized towards maximum-likelihood objectives. In that case, the neural network can be interpreted as a probabilistic model $\log p_{\boldsymbol{\theta}}(\mathbf{t}|\mathbf{x},\boldsymbol{\epsilon})$ computing a conditional distribution over some target variable $\mathbf{t}$ (e.g. classes) given some input $\mathbf{x}$. In (Bengio & Thibodeau-Laufer, 2013), stochastic hidden units are used for learning the parameters of a Markov chain transition operator that samples from the data distribution.

While ’dropout’ is designed as a regularization method, other work on stochastic neural networks exploit the power of stochastic hidden units for generative modeling, e.g. (Frey & Hinton, 1999; Rezende et al., 2014; Tang & Salakhutdinov, 2013) applying (partially) MCMC or (partically) factorized variational approaches to modelling the posterior. As we will see in sections 4 and 6, the choice of parameterization has a large impact on the posterior dependencies and the efficiency of posterior inference. However, current publications lack a good justification for their choice of parameterization. The analysis in section 4 offers some important insight in where the centered or non-centered parameterizations of such networks are more appropriate.

4 A differentiable MC likelihood estimator

We showed that many hierarchical continuous latent-variable models can be transformed into a DNCP $p_{\boldsymbol{\theta}}(\mathbf{x},\boldsymbol{\epsilon})$, where all latent variables (the introduced auxiliary variables $\boldsymbol{\epsilon}$) are root nodes (see eq. (LABEL:eq:integrating_out_z)). This has an important implication for learning since (contrary to a CP) the DNCP can be used to form a differentiable Monte Carlo estimator of the marginal likelihood:

where the parents $\mathbf{pa}_{j}^{(l)}$ of the observed variables are either root nodes or functions of root nodes whose values are sampled from their marginal: $\boldsymbol{\epsilon}^{(l)}\sim p(\boldsymbol{\epsilon})$. This MC estimator can be differentiated w.r.t. $\boldsymbol{\theta}$ to obtain an MC estimate of the log-likelihood gradient $\nabla_{\boldsymbol{\theta}}\log p_{\boldsymbol{\theta}}(\mathbf{x})$, which can be plugged into stochastic optimization methods such as Adagrad for approximate ML or MAP. When performed one datapoint at a time, we arrive at our on-line Maximum Monte Carlo Likelihood (MMCL) algorithm.

Effects of parameterizations on posterior dependencies

What is the effect of the proposed reparameterization on the efficiency of inference? If the latent variables have linear-Gaussian conditional distributions, we can use the metric of squared correlation between the latent variable and any of its children in their posterior distribution. If after reparameterization the squared correlation is decreased, then in general this will also result in more efficient inference.

For non-linear Gaussian conditional distributions, the log-PDF can be locally approximated as a linear-Gaussian using a second-order Taylor expansion. Results derived for the linear case can therefore also be applied to the non-linear case; the correlation computed using this approximation is a local dependency between the two variables.

Denote by $z$ a scalar latent variable we are going to reparameterize, and by $\mathbf{y}$ its parents, where $y_{i}$ is one of the parents. The log-PDF of the corresponding conditional distribution is

A reparameterization of $z$ using an auxiliary variable $\epsilon$ is $z=g(.)=(\mathbf{w}^{T}\mathbf{y}+b)+\sigma\epsilon$ where $\epsilon\sim\mathcal{N}(0,1)$. With (7) it can be confirmed that this change of variables is correct:

First we will derive expressions for the squared correlations between $z$ and its parents, for the CP and DNCP case, and subsequently show how they relate.

The covariance $C$ between two jointly Gaussian distributed variables $A$ and $B$ equals the negative inverse of the Hessian matrix of the log-joint PDF:

The correlation $\rho$ between two jointly Gaussian distributed variables $A$ and $B$ is given by: $\rho=\sigma_{AB}^{2}/(\sigma_{A}\sigma_{B})$. Using the equation above, the squared correlation can be computed from the elements of the Hessian matrix:

Important to note is that derivatives of the log-posterior w.r.t. the latent variables are equal to the derivatives of log-joint w.r.t. the latent variables, therefore,

The following shorthand notation is used in this section:

In the CP case, the relevant Hessian elements are as follows:

Therefore, using eq. (10), the squared correlation between $y_{i}$ and $z$ is:

1.2 Non-centered case

In the DNCP case, the Hessian elements are:

The squared correlation between $y_{i}$ and $\epsilon$ is therefore:

2 Correlation inequality

We can now compare the squared correlation, between $z$ and some parent $y_{i}$, before and after the reparameterization. Assuming $\alpha<0$ and $\beta<0$ (i.e. $L^{(\setminus z)}$ and $L^{(z\rightarrow)}$ are concave, e.g. exponential families):

Thus we have shown the surprising fact that the correlation inequality takes on an extremely simple form where the parent-dependent values $\alpha$ and $w_{i}$ play no role; the inequality only depends on two properties of $z$: the relative strenghts of $\sigma$ (its noisiness) and $\beta$ (its influence on children’s factors). Informally speaking, if the noisiness of $z$’s conditional distribution is large enough compared to other factors’ dependencies on $z$, then the reparameterized form is beneficial for inference.

3 A beauty-and-beast pair

Additional insight into the properties of the CP and DNCP can be gained by taking the limits of the squared correlations (12) and (14). Limiting behaviour of these correlations is shown in table 1. As becomes clear in these limits, the CP and DNCP often form a beauty-and-beast pair: when posterior correlations are high in one parameterization, they are low in the other. This is especially true in the limits of $\sigma\to 0$ and $\beta\to-\infty$, where squared correlations converge to either or $1$, such that posterior inference will be extremely inefficient in either CP or DNCP, but efficient in the other. This difference in shapes of the log-posterior is illustrated in figure 3.

4 Example: Simple Linear Dynamical System

Take a simple model with scalar latent variables $z_{1}$ and $z_{2}$, and scalar observed variables $x_{1}$ and $x_{2}$. The joint PDF is defined as $p(x_{1},x_{2},z_{1},z_{2})=p(z_{1})p(x_{1}|z_{1})p(z_{2}|z_{1})p(x_{2}|z_{2})$, where $p(z_{1})=\mathcal{N}(0,1)$, $p(x_{1}|z_{1})=\mathcal{N}(z_{1},\sigma_{x}^{2})$, $p(z_{2}|z_{1})=\mathcal{N}(z_{1},\sigma_{z}^{2})$ and $p(x_{2}|z_{2})=\mathcal{N}(z_{2},\sigma_{x}^{2})$. Note that the parameter $\sigma_{z}$ determines the dependency between the latent variables, and $\sigma_{x}$ determines the dependency between latent and observed variables.

We reparameterize $z_{2}$ such that it is conditionally deterministic given a new auxiliary variable $\epsilon_{2}$. Let $p(\epsilon_{2})=\mathcal{N}(0,1)$. let $z_{2}=g_{2}(z_{1},\epsilon_{2},\sigma_{z})=z_{1}+\sigma_{z}\cdot\epsilon_{2}$ and let $\epsilon_{1}=z_{1}$. See figure 3 for plots of the original and auxiliary posterior log-PDFs, for different choices of $\sigma_{z}$, along with the resulting posterior correlation $\rho$.

For what choice of parameters does the reparameterization yield smaller posterior correlation? We use equation (LABEL:eq:corr_inequality) and plug in $\sigma\leftarrow\sigma_{z}$ and $-\beta\leftarrow\sigma^{-2}_{x}$, which results in:

i.e. the posterior correlation in DNCP form $\rho_{\epsilon_{1},\epsilon_{2}}^{2}$ is smaller when the latent-variable noise parameter $\sigma^{2}_{z}$ is smaller than the oberved-variable noise parameter $\sigma^{2}_{x}$. Less formally, this means that the DNCP is preferred when the latent variable is more strongly coupled to the data (likelihood) then to its parents.

Related work

This is, to the best of our knowledge, the first work to investigate the implications of the different differentiable non-centered parameterizations on the efficiency of gradient-based inference. However, the topic of centered vs non-centered parameterizations has been investigated for efficient (non-gradient based) Gibbs Sampling in work by Papaspiliopoulos et al. (2003; 2007), which also discusses some strategies for constructing parameterization for those cases. There have been some publications for parameterizations of specific models; (Gelfand et al., 1995), for example, discusses parameterizations of mixed models, and (Meng & Van Dyk, 1998) investigate several rules for choosing an appropriate parameterization for mixed-effects models for faster EM. In the special case where Gibbs sampling is tractable, efficient sampling is possible by interleaving between centered and non-centered parameterizations, as was shown in (Yu & Meng, 2011).

Auxiliary variables are used for data augmentation (see (Van Dyk & Meng, 2001) or slice sampling (Neal, 2003)) where, in contrast with our method, sampling is performed in a higher-dimensional augmented space. Auxiliary variables are used in a similar form under the name exogenous variables in Structural Causal Models (SCMs) (Pearl, 2000). In SCMs the functional form of exogenous variables is more restricted than our auxiliary variables. The concept of conditionally deterministic variables has been used earlier in e.g. (Cobb & Shenoy, 2005), although not as a tool for efficient inference in general Bayesian networks with continuous latent variables. Recently, (Raiko et al., 2012) analyzed the elements of the Hessian w.r.t. the parameters in neural network context.

The differentiable reparameterization of latent variables in this paper was introduced earlier in (Kingma, 2013) and independently in (Bengio, 2013), but these publications lack a theoretic analysis of the impact on the efficiency of inference. In (Kingma & Welling, 2013), the reparameterization trick was used in an efficient algorithm for stochastic variational inference and learning.

Experiments

From the derived posterior correlations in the previous sections we can conclude that depending on the parameters of the model, posterior sampling can be extremely inefficient in one parameterization while it is efficient in the other. When the parameters are known, one can choose the best parameterization (w.r.t. posterior correlations) based on the correlation inequality (LABEL:eq:corr_inequality).

In practice, model parameters are often subject to change, e.g. when optimizing the parameters with Monte Carlo EM; in these situations where there is uncertainty over the value of the model parameters, it is impossible to choose the best parameterization in advance. The ”beauty-beast” duality from section 4.3 suggests a solution in the form of a very simple sampling strategy: mix the two parameterizations. Let $Q_{CP}(\mathbf{z}^{\prime}|\mathbf{z})$ be the MCMC/HMC proposal distribution based on $p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x})$ (the CP), and let $Q_{DNCP}(\mathbf{z}^{\prime}|\mathbf{z})$ be the proposal distribution based on $p_{\boldsymbol{\theta}}(\boldsymbol{\epsilon}|\mathbf{x})$ (the DNCP). Then the new MCMC proposal distribution based on the mixture is:

where we use $\rho=0.5$ in experiments. The mixing efficiency might be half that of the oracle solution (where the optimal parameterization is known), nonetheless when taking into account the uncertainty over the parameters, the expected efficiency of the mixture proposal can be better than a single parameterization chosen ad hoc.

We applied a Hybrid Monte Carlo (HMC) sampler to a Dynamic Bayesian Network (DBN) with nonlinear transition probabilities with the same structure as the illustrative model in figure 2. The prior and conditional probabilities are: $\mathbf{z}_{1}\sim\mathcal{N}(0,\mathbf{I})$, $\mathbf{z}_{t}|\mathbf{z}_{t-1}\sim\mathcal{N}(tanh(\mathbf{W}_{z}\mathbf{z}_{t-1}+\mathbf{b}_{z}),\sigma_{z}^{2}\mathbf{I})$ and $\mathbf{x}_{t}|\mathbf{z}_{t}\sim\text{Bernoulli}(sigmoid(\mathbf{W}_{x}\mathbf{z}_{t-1}))$. The parameters were intialized randomly by sampling from $\mathcal{N}(0,\mathbf{I})$. Based on the derived limiting behaviour (see table 1, we can expect that such a network in CP can have very large posterior correlations if the variance of the latent variables $\sigma^{2}_{z}$ is very small, resulting in slow sampling.

To validate this result, we performed HMC inference with different values of $\sigma^{2}_{z}$, sampling the latent variables while holding the parameters fixed. For HMC we used 10 leapfrog steps per sample, and the stepsize was automatically adjusted while sampling to obtain a HMC acceptance rate of around 0.9. At each sampling run, the first 1000 HMC samples were thrown away (burn-in); the subsequent 4000 HMC samples were kept. To estimate the efficiency of sampling, we computed the effective sample size (ESS); see e.g. (Kass et al., 1998) for a discussion on ESS.

Results. See table 2 and figure 4 for results. It is clear that the choice of parameterization has a large effect on posterior dependencies and the efficiency of inference. Sampling was very inefficient for small values of $\sigma_{z}$ in the CP, which can be understood from the limiting behaviour in table 1.

2 Generative multilayer neural net

As explained in section 3.4, a hierarchical model in DNCP form can be learned using a MC likelihood estimator which can be differentiated and optimized w.r.t. the parameters $\boldsymbol{\theta}$. We compare this Maximum Monte Carlo Likelihood (MMCL) method with the MCEM method for learning the parameters of a 4-layer hierarchical model of the MNIST dataset, where $\mathbf{x}|\mathbf{z}_{3}\sim\text{Bernoulli}(sigmoid(\mathbf{W}_{x}\mathbf{z}_{3}+\mathbf{b}_{x}))$ and $\mathbf{z}_{t}|\mathbf{z}_{t-1}\sim\mathcal{N}(\tanh(\mathbf{W}_{i}\mathbf{z}_{t-1}+\mathbf{b}_{i}),\sigma^{2}_{z_{t}}\mathbf{I})$. For MCEM, we used HMC with 10 leapfrog steps followed by a weight update using Adagrad (Duchi et al., 2010). For MMCL, we used $L\in\{10,100,500\}$. We observed that DNCP was a better parameterization than CP in this case, in terms of fast mixing. However, even in the DNCP, HMC mixed very slowly when the dimensionality of latent space become too high. For this reason, $\mathbf{z}_{1}$ and $\mathbf{z}_{2}$ were given a dimensionality of 3, while $\mathbf{z}_{3}$ was 100-dimensional but noiseless ($\sigma^{2}_{z_{1}}=0$) such that only $\mathbf{z}_{3}$ and $\mathbf{z}_{2}$ are random variables that require posterior inference by sampling. The model was trained on a small (1000 datapoints) and large (50000 datapoints) version of the MNIST dataset.

Results. We compared train- and testset marginal likelihood. See figure 5 for experimental results. As was expected, MCEM attains asymptotically better results. However, despite its simplicity, the on-line nature of MMCL means it scales better to large datasets, and (contrary to MCEM) is trivial to implement.

Conclusion

We have shown how Bayesian networks with continuous latent variables and generative neural networks are related through two different parameterizations of the latent variables: CP and DNCP. A key result is that the differentiable non-centered parameterization (DNCP) of a latent variable is preferred, in terms of its effect on decreased posterior correlations, when the variable is more strongly linked to its parents than its children. Through theoretical analysis we have also shown that the two parameterizations are complementary to each other: when posterior correlations are large in one form, they are small in the other. We have also illustrated that this theoretical result can be exploited in practice by designing a MCMC strategy that mixes between both parameterizations, making it robust to situations where MCMC can otherwise be inefficient.

Acknowledgments

The authors thank the reviewers for their excellent feedback and Joris Mooij, Ted Meeds and Taco Cohen for invaluable discussions and input.

References