Log In Sign Up

Bidirectional Helmholtz Machines

Efficient unsupervised training and inference in deep generative models remains a challenging problem. One basic approach, called Helmholtz machine, involves training a top-down directed generative model together with a bottom-up auxiliary model used for approximate inference. Recent results indicate that better generative models can be obtained with better approximate inference procedures. Instead of improving the inference procedure, we here propose a new model which guarantees that the top-down and bottom-up distributions can efficiently invert each other. We achieve this by interpreting both the top-down and the bottom-up directed models as approximate inference distributions and by defining the model distribution to be the geometric mean of these two. We present a lower-bound for the likelihood of this model and we show that optimizing this bound regularizes the model so that the Bhattacharyya distance between the bottom-up and top-down approximate distributions is minimized. This approach results in state of the art generative models which prefer significantly deeper architectures while it allows for orders of magnitude more efficient approximate inference.


page 4

page 8


Reweighted Expectation Maximization

Training deep generative models with maximum likelihood remains a challe...

Learning to Bound: A Generative Cramér-Rao Bound

The Cramér-Rao bound (CRB), a well-known lower bound on the performance ...

Stochastic Backpropagation and Approximate Inference in Deep Generative Models

We marry ideas from deep neural networks and approximate Bayesian infere...

Multi-objects Generation with Amortized Structural Regularization

Deep generative models (DGMs) have shown promise in image generation. Ho...

Bisimulation-based Approximate Lifted Inference

There has been a great deal of recent interest in methods for performing...

Learning and Inference Movement with Deep Generative Model

Learning and inference movement is a very challenging problem due to its...

Asymptotically exact inference in differentiable generative models

Many generative models can be expressed as a differentiable function of ...

Code Repositories


Bidirectional Helmholtz Machines

view repo

1 Introduction and background

Training good generative models and fitting them to complex and high dimensional training data with probability mass in multiple disjunct locations remains a major challenge. This is especially true for models with multiple layers of deterministic or stochastic variables, which is unfortunate because it has been argued previously 

(Hinton et al., 2006; Bengio, 2009) that deeper generative models have the potential to capture higher-level abstractions and thus generalize better. Although there has been progress in dealing with continous-valued latent variables (Kingma & Welling, 2014), building a hierarchy of representations, especially with discrete-valued latent variables, remains a challenge.

With the Helmholtz machine (Hinton et al., 1995; Dayan et al., 1995), a concept was introduced that proposed to not only fit a powerful but intractable generative model to the training data, but also to jointly train a parametric approximate inference model . The distribution is used to perform approximate inference over the latent variables of the generative model given an observed example , i.e. to approximate . This basic idea has been applied and enhanced many times; initially with the wake-sleep algorithm (WS, Hinton et al. (1995); Dayan & Hinton (1996)) and more recently with the variational autoencoder (VAE, Kingma & Welling (2014)

), stochastic backpropagation and approximate inference in deep generative models 

(Rezende et al., 2014), neural variational inference and learning (NVIL, Mnih & Gregor (2014)) and reweighted wake-sleep (RWS, Bornschein & Bengio (2015)).

Recent results indicate that significant improvements can be made when better approximate inference methods are used: Salimans et al. (2014) for example presented an iterative inference procedure that improves the samples from by employing a learned MCMC transition operator. In (Hjelm et al., 2015) the authors propose an iterative, gradient based procedure to refine the approximate posterior. Burda et al. (2015) present the importance weighted auto encoder (IWAE), an improved VAE that, similarly to RWS, uses multiple samples from to calculate gradients. And RWS already reported that autoregressive distributions lead to noticable improvements. In contrast to these previous approaches that aimed at incorporating more powerful inference methods to gain better approximations of , we here propose to regularize the top-down model such that its posterior stays close to the approximate inference distribution (and vice versa). We archive this by interpreting both and as approximate inference models for our actual generative model , which is defined to be the geometric mean over the top-down and bottom-up approximate inference models, i.e., .

In Section 2 we show that this model definition leads to an objective that can be interpreted as using a regularization term that encurages solutions where and are close to each other in terms of the Bhattacharyya distance. In Section 3 we will explain how to perform importance sampling based training and inference. The ability to model complex distributions and the computational efficiency of this approach are demonstrated empirically in Section 4.

2 Model definition and properties

We introduce the bidirectional Helmholtz Machine (BiHM) by defining a joint probability distribution over three variable vectors, an observed vector

and two latent variable vectors and

. Analogous to a Deep Boltzmann Machine (DBM, 

Salakhutdinov & Hinton (2009)

), we think of these as layers in a neural network with links between

and on the one side, and and on the other side. We will present our approach for the specific case of an architecture with two hidden layers, but it can be applied to arbitrary graphs of variables without loops. It can especially be used to train architectures with more than two stacked layers of latent variables.

Let be a joint probability distribution constructed in a specific way from two constituent distributions and ,

where is a normalization constant and and are directed graphical models from to and vice versa,

We assume that the prior distribution and all conditional distributions belong to parametrized families of distributions which can be evaluated and sampled from efficiently. For we do not assume an explicit form but define it to be the marginal

The normalization constant guarantees that . Using the Cauchy-Schwarz inequality and identifying with and with , it becomes clear that for arbitrary and . Furthermore, we see that if and only if . We can therefore obtain a lower bound on the marginal probability by defining


This suggests that the model distribution can be fitted to some training data by maximizing the bound of the log-likelihood (LL) instead of , as we elaborate in the following section. Since can reach the maximum only when , the model is implicitly pressured to find a maximum likelihood solution that yields .

2.1 Alternative view based on the Bhattacharyya distance

Recalling the Bhattacharyya distance (for which holds for arbitrary distributions and only if ) the model LL can be decomposed into


where we clearly see that the proposed training objective corresponds to the LL minus 2 times the Bhattacharyya distance , i.e., it is maximzing the true LL and minimizing the distance between and . We can compare this to the variational approach, where the marginal probability of some model containing latent variables is rewritten in terms of the KL-divergence to obtain a lower bound


Analogous to variational methods that maximize the lower bound (3), we can thus maximize , and it will tighten the bound as approaches zero. While this seems very similar to the variational lower bound, we should highlight that there are some important conceptual differences: 1) The KL-divergence in variational methods measures the distance between distributions given some training data. The Bhattacharyya distance here in contrast quantifies a property of the model independently of any training data. In fact, we saw that . 2) The variational lower bound is typically used to construct approximate inference algorithms. We here use our bound just to remove the normalization constant from our target distribution . Even after applying the lower-bound, we still have to tackle the inference problem which manifests itself in form of the full combinatorial sum over and in equation (1). Although it seems intuitively reasonable to use a variational approximation on top of the bound we will here not follow this direction but rather use importance sampling to perform approximate inference and learning (see section 3). Combining a variational method with the bound is therefore subject to future work.

We can also argue that optimizing instead of is beneficial in the light of the original goal we formulated in section 1: To learn a generative model that is regularized to be close to the model which we use to perform approximate inference for . Let us assume we have two equally well trained models and , i.e., in expectation over the empirical distribution , but the expected bound for the first model is closer to the LL than the expected bound for the second model: . Using equation (2) we see that which indicates that is closer to than is to (when we measure their distance using the Bhattacharyya distance). According to our original goal, we thus prefer solution , where the bound is maximized and the distance minimized.

Note that the decomposition (2) also emphasizes why our recursive definition is a consistent and reasonable one: minimizing

during learning means that the joint distributions

and approach each other. This implies that the marginals and for all layers become more similar. This also implies in the limit of ; a requirement that most simple parametrized distributions could never fulfill.

Figure 1: MNIST experiments: A) Random samples from top-down model . B) Generally improved samples after running iterations of Gibbs sampling to obtain approximate samples from the joint model . In A and B

we show expected samples instead of sampling from the bottom Bernoulli distribution, as usual.

C) Sensitivity of the test set LL estimates to the number of samples . We plot the test set and estimates of our best BiHM model together with the estimates of our best RWS and VAE models.

Figure 2:

Inpainting of binarized MNIST digits. The left column in each block shows the original digit randomly sampled from the MNIST test set; the first column shows the masked version presented to the algorithm, and the next three columns show different, independent reconstructions after 100 Gibbs iterations. (see section

3). All images were selected randomly.

3 Inference and training

  for number of training iterations do
      Sample from the training distribution (i.e. )
     for  do
         Sample       (for layers to )
         Compute and       (for )
     end for
      Compute unnormalized importance weights      
      Normalize the weights
      Update parameters of and : gradient descent     with gradient estimator      
  end for
Algorithm 1 Training using importance samples

Based on the construction of outlined in the previous section we can define a wide range of possible models. Furthermore, we have a wide range of potential training and appropriate inference methods we could employ to maximize .

In this text we concentrate on binary latent and observed variables and model all our conditional distributions by simple sigmoid belief network layers, e.g., where refers to the Bernoulli distribution with , are the connection weights between the latent variables and the visible variable ; is the bias of , and

is the sigmoid function. For our top-level prior

, we use a factorized Bernoulli distribution: .

We form an estimate of by using importance sampling instead of the exhaustive sum over and in equation (1). We use as the proposal distribution which is by construction easy to evaluate and to sample from:


with and . Using the same approach, we can also derive the well known estimator for the marginal probability of a datapoint under the top-down generative model :


Comparing (4) and (5) and making use of Jensen’s inequality it becomes clear that .

Analogous to the parameter updates in RWS (Bornschein & Bengio, 2015), we can derive an importance sampling based estimate for the LL gradient with respect to the parameters of and (jointly denoted by ) and use it to optimize our objective (we use to jointly denote the latent variables of all layers)


with and importance weights

for .

In contrast to VAEs and IWAEs, the updates do not require any form of backpropagation through layers because, as far as the gradient computation is concerned, these samples are considered fully observed. The gradient approximation (6) computes the weighted average over the individual gradients. These properties are basically inherited from the RWS training algorithm. But in contrast to RWS, and in contrast to most other algorithms which employ a generative model and an approximate inference model , we here automatically obtain parameter updates for both and because we optimize which contains both. The resulting training method is summarized in algorithm 1.

Estimating the partition function

To compute and to monitor the training progress it is desirable to estimate the normalization constant . In stark contrast to undirected models like RBMs or DBMs, we can here derive an unbiased importance sampling based estimator for :


We denote the number of samples used to approximate the outer expectation and the inner expectation with and respectively. In the experimental section, we show that we obtain high quality estimates for with =1 and a relatively small number of samples . By taking the logarithm, we obtain a biased estimator for , which will, unfortunately, underestimate

on average due to the concavity of the logarithm and the variance of the

estimate. This can lead to overestimates for (see equation (4)) if we are not careful. Fortunately, the bias on the estimated is induced only by the concavity of the logarithm; the underlying estimator for is unbiased. We can thus effectively minimize the bias by minimizing the variance of the estimate (e.g. by taking more samples). This is a much better situation than for

-estimating methods that rely on Markov chains in high dimensional spaces, which might miss entire modes because of mixing issues.

Figure 3: estimates for different values of as a function of A) the number of samples , B) the total number of samples

for the BiHM trained on MNIST; the gray region shows the mean and the standard deviation for 10 runs with

=1. This shows that, from the point of view of total computation, convergence is fastest with =1; and that we obtain a high quality estimate of the partition function with only a few million samples. C) Evolution of the estimates of , , and during training on MNIST.
auto regressive
 NADE 13.19 11.99 84.81 9.81 273.08 27.22 46.66 28.39
 EoNADE 13.19 12.58 82.31 9.68 272.38 27.31 46.12 27.87
 DARN 13.19 11.91 81.04 9.55 274.68 28.17 46.10 28.83
 RWS - NADE 13.16 11.68 84.26 9.71 271.11 26.43 46.09 27.92
non AR
 RBM 16.26 22.66 96.74 15.15 277.37 43.05 48.88 29.38
 RWS - SBN 13.65 12.68 90.63 9.90 272.54 29.99 46.16 28.18
 - 13.78 12.43 86.49 9.40 272.66 27.10 46.12 28.14
 - 13.82 12.31 86.92 9.40 272.71 27.30 46.98 28.22
 - 0.20 0.27 0.56 0.09 1.97 1.87 0.41 0.54
  ess 81.5% 89.1% 11.2% 92.5% 16.8% 22.5% 70.6% 55.9%
Table 1:

Negative log-likelihood (NLL) on various binary datasets from the UCI repository: The top rows quote results from shallow models with autoregressive weights between their units within one layer. The second block shows results from non-autoregressive models (quoted from

Bornschein & Bengio (2015)). In the third block we show the results obtained by training a BiHMs. We report the estimated test set NLL when evaluating just the top-down model, , and when evaluating . Our BiHM models consistently obtain similar or better results than RWS while they prefer deeper architectures. All effective sample size (see 4.1) estimates are with error bounds smaller than .

Sampling and inpainting

We now discuss two general approaches for approximate sampling from a BiHM. One can either easily and efficiently sample from the directed model , or one can use Gibbs sampling to draw higher-quality samples from the undirected model . For the latter, importance resampling is used to approximately draw samples from the conditional distributions, e.g. from:

Here we choose to draw the proposal samples from the mixture distribution , which ensures that we have a symmetric chance of covering the high probability configurations of induced by and . We resample a final sample from propotionally to their importance weights which are thus given by

where is randomly drawn from or and the constant collects all terms not containing and can be ignored. For we choose to sample by drawing the proposal samples from

. We iteratively update all odd layers followed by all even layers until we consider the chain to be in equilibrium.

Equipped with approximate sampling procedures for the conditional distributions, it is straightforward to construct an algorithm for inpainting: Given a corrupted input datapoint , we first initialize a Markov chain by drawing and then run the Gibbs sampling procedure. Whenever we sample the bottom layer , we keep the non-corrupted elements of fixed.

4 Experimental results

In this section we present experimental results obtained when applying the algorithm to various binary datasets. Our main goal is to ensure that the theoretical properties discussed in section 2 translate into a robust algorithm that yields competitive results even when used with simple sigmoid belief network layers as conditional distributions. We train all models using Adam (Kingma & Ba, 2014) with a mini-batch size of 100. We initialize the weights according to Glorot & Bengio (2010), set the biases to -1, and use regularization = on all the weights. Our implementation is available at

UCI binary datasets

To ascertain that importance sampling based training of BiHMs works in general, we applied it to the 8 binary datasets from the UCI dataset repository that were evaluated e.g. in (Larochelle & Murray, 2011). We use a learning rate of or for all the experiments. The architectures, layer sizes and final LL estimates can be found in tables 1 and 2.

Dataset BiHM layer sizes RWS layer sizes
 ADULT 100, 70, 50, 25 100, 20, 5
 CONNECT4 300, 110, 50, 20 150, 50, 10
 DNA 200,150,130,100,70,50,30,20,10 150, 10
 MUSHROOM 150, 100, 90, 60, 40, 20 150, 50, 10
 NIPS-0-12 200, 100, 50, 25 150, 50, 10
 OCR 600, 500, 100, 50, 30, 10 300, 100, 10
 RCV1 500, 120, 70, 30 200, 50, 10
 WEB 650, 580, 70, 30, 10 300, 50, 10
Table 2: Architectures for our best UCI BiHM models compared to our best RWS models. We observe that BiHMs prefer significantly deeper architectures than RWS.

Binarized MNIST

We use the MNIST dataset that was binarized according to (Murray & Salakhutdinov, 2009) and which we downloaded in binarized form  (Larochelle, 2011). Compared to RWS, we again observe that BiHMs prefer significantly deeper and narrower models. Our best model consists of 1 visible and 12 latent layers with 300,200,100,75,50,35,30,25,20,15,10,10 latent variables. We follow the same experimental procedure as in the RWS paper: First train the model with =10 samples and a learning rate of until convergence and then fine-tune the parameters with =100 samples and a learning rate of . All layers are actually used to model the empirical distribution; we confirmed that training shallower models (obtained by leaving out individual layers) decreases the performance. We obtain test set log-loglikelihoods of -84.6 0.23 and -84.3 0.22 111Recently, a version of MNIST where binary observed variables are resampled during training has been used (e.g., Burda et al. (2015)). On this dataset we obtain -82.9. The next section presents a more detailed analysis of these estimates and their dependency on the number of samples from the proposal distribution . Note that even though this model is relatively deep, it is not particularly large, with about parameters in total. The DBMs in (Salakhutdinov & Hinton, 2009) contain about and 1.1 million parameters; a variational autoencoder with two deterministic, 500 units wide encoder and decoder layers, and with 100 top level latent units contains more than 1.4 million parameters.

To highlight the models ability to generate crisp (non-blurry) digits we draw approximate samples from which are visualized in Fig. 1 B. Fig. 1 A shows samples obtained when drawing from the top-down generative model before running any Gibbs iterations. Fig. 2 visualizes the results when running the inpainting algorihm to reconstruct partially occluded images. Our sourcecode package contains additional results and animations.

Model - -
autoregressive models
  NADE 88.9
  DARN (1 latent layer) 88.3 84.1
continuous latent variables
  VAE 96.2 88.7
  VAE + HMC (8 iterations) 88.3 85.5
  IWAE (2 latent layers) 96.1 85.3
binary latent variables
  NVIL (2 latent layers) 99.6
  RWS (5 latent layers) 96.3 85.4
BiHM (12 latent layers) 89.2 84.3
Table 3: Comparison of BiHMs to other recent methods in the literature. We report the lower bounds and estimates for the marginal log probability on the binarized MNIST test set.

Toronto Face Database

We also trained models on the 98,058 examples from the unlabeled section of the Toronto face database (TFD, Susskind et al. (2010)). Each training example is of size

pixels and we interpret the gray-level as Bernoulli probability for the bottom layer. We observe that training proceeds rapidly during the first few epochs but mostly only learns the mean-face. During the next few hundred epochs training proceeds much slower but the estimated log-likelihood

increases steadily. Fig. 5 A shows random samples from a model with respectively 1000,700,700,300 latent variables in the 4 hidden layers. It was trained with a learning rate of

; all other hyperparameters were set to the same values as before. Fig.

5 B shows the results from inpainting experiments with this model.

Figure 4: Learning curves for MNIST experiments: For BiHM, RWS and VAE we chose the learning rate such that we get optimal results after update steps; for IWAE we use the original learning rate schedule published in (Burda et al., 2015). BiHM and RWS use =10 samples per datapoint; IWAE uses =5 and a batch size of 20 (according to the original publication). We generally observe that BiHM show very good convergence in terms of progress per update step and in terms of total training time.
Figure 5: Results after training on TFD: A) Random selection of 12 samples drawn from (10 iterations of Gibbs sampling). B) The left column in each block shows the input; the right column shows a random output sample generated by the inpaiting algorithm (see section 3).

4.1 Analysis of importance sampling-based estimates

Estimating the partition function

In Fig. 3 A we plot estimates (equation (7)) over the number of outer samples for our best MNIST model and for 3 different choices of , i.e.,. In Fig. 3 B we plot the estimates over the total number of samples . We observe that choosing =1 and using only about 10 million samples results in high quality estimates for

with an standard error far below 0.1 nats. Estimating based on 10 million samples takes less than 2 minutes on a GTX980 GPU. Fig.

3 C shows the development of the estimate during learning and in relation to the LL estimates.

Importance sampling efficiency

A widely used metric to estimate the quality of an importance sampling estimator is the effective sampling size (ESS), given by (see, e.g., Robert & Casella (2009)). Larger values indicate more efficient sampling (more information extracted per sample). We compute the ESS over the MNIST test set for =100,000 proposal samples from . For our best RWS model, a model with 5 stochastic layers (400,300,200,100,10), we obtain ; for the BiHM model we obtain . When we estimate the ESS for using from the BiHM as a proposal distribution for , we obtain =. The estimated ESS values indicate that training BiHM models indeed results in distributions whose intractable posterior as well as top-down model are much better modeled by the learned

. We also estimated the ESS for a VAE with two determninistic, 500 units wide ReLU layers in the encoder and decoder. This model has a single stochastic layer with 100 continuous variables at the top; it reaches a final estimated test set LL of

. The final variational lower bound, which corresponds exactly to the importance sampling estimate of with =1 sample, is . For this model we obtain an ESS of . These results indicate that we need thousands of samples to obtain reliable LL estimates with low approximation error. In Fig. 1 C we plot the estimated test set LL over the number of samples used to estimate and . For all the models and for small a number of samples we significantly underestimate the LL; but, in comparison to RWS, the estimates for the BiHM model are much higher and less sensitive to . E.g, using =10 samples to evaluate the BiHM model results in a higher LL estimate than using =10,000 samples to evaluate the RWS model.

Computational cost

To demonstrate the computational efficiency of our approach we show typical MNIST learning curves in Fig. 4. For BiHM, RWS and VAE the learning rate was chosen within a factor of 2 to obtain optimal results after update steps ( for BiHM and RWS, for VAE; =10 for BiHM and RWS). For the IWAE experiment we use the original code, hyperparameters and learning rate schedule from (Burda et al., 2015): This experiment thus uses a mini-batch size of 20 instead of 100, =5 training samples and and 8 different learning rates over the course of 3300 epochs. In all cases we used =1000 samples to evaluate the test set log-likelihoods. We generally observe that BiHM show very good convergence in terms of progress per update step and competitive performance in terms of total training time. Note that BiHMs and RWS allow for an efficient distributed implementation in the future: per sample, only the binary activations and a single floating point number (the importance weight) need to be communicated between layers. VAEs and IWAEs need to communicate continuous activations during the forward pass and continuous partial gradients during the backward pass. At test time BiHMs are typically much more effective than the other methods: BiHMs obtain good LL estimates with =10 or 100 samples per datapoint while VAE, RWS and IWAE models require 10,000 samples to obtain competitive results (compare Fig. 1 C).

5 Conclusion and future work

We introduced a new scheme to construct probabilistic generative models which are automatically regularized to be close to approximate inference distributions we have at our disposal. Using the Bhattacharyya distance we derived a lower-bound on the log-likelihood, and we demonstrated that the bound can be used to fit deep generative models with many layers of latent variables to complex training distributions.

Compared to RWS, BiHM models typically prefer many more latent layers. After training a BiHM, the directed top-down model shows better performance than a RWS trained model; both in terms of log-likelihood and sample quality. Sample quality can be further improved by approximately sampling from the full undirected BiHM model . The high similarity between and , enforced by the training objective, allows BiHMs to be evaluated orders of magnitude more efficiently than RWS, VAE and IWAE models.

Possible directions for future research could involve semisupervised learning: the symmetric nature of the generative model (it is always close to the bottom-up and top-down directed models and ) might make it particularly interesting for learning tasks that require inference given changing sets of observed and hidden variables. We also have a wide range of potential choices for our parametrized conditional distributions. Assuming continuous latent variables for example and eventually choosing an alternative inference method might make a better suited model for some training distributions.


We thank the developers of Theano 

(Bergstra et al., 2010) and Blocks (Van Merriënboer et al., 2015) for their great work. We thank NSERC, Compute Canada, CIFAR and Samsung for their support.