Neural Autoregressive Distribution Estimation

05/07/2016 ∙ by Benigno Uria, et al. ∙ Université de Sherbrooke Twitter 0

We present Neural Autoregressive Distribution Estimation (NADE) models, which are neural network architectures applied to the problem of unsupervised distribution and density estimation. They leverage the probability product rule and a weight sharing scheme inspired from restricted Boltzmann machines, to yield an estimator that is both tractable and has good generalization performance. We discuss how they achieve competitive performance in modeling both binary and real-valued observations. We also present how deep NADE models can be trained to be agnostic to the ordering of input dimensions used by the autoregressive product rule decomposition. Finally, we also show how to exploit the topological structure of pixels in images using a deep convolutional architecture for NADE.



There are no comments yet.


page 19

page 23

page 29

page 32

Code Repositories


Neural Autoregressive Distribution Estimation

view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Distribution estimation is one of the most general problems addressed by machine learning. From a good and flexible distribution estimator, in principle it is possible to solve a variety of types of inference problem, such as classification, regression, missing value imputation, and many other predictive tasks.

Currently, one of the most common forms of distribution estimation is based on directed graphical models. In general these models describe the data generation process as sampling a latent state from some prior , followed by sampling the observed data from some conditional . Unfortunately, this approach quickly becomes intractable and requires approximations when the latent state increases in complexity. Specifically, computing the marginal probability of the data, , is only tractable under fairly constraining assumptions on and .

Another popular approach, based on undirected graphical models, gives probabilities of the form , where is a tractable function and is a normalizing constant. A popular choice for such a model is the restricted Boltzmann machine (RBM), which substantially out-performs mixture models on a variety of binary datasets (Salakhutdinov and Murray, 2008). Unfortunately, we often cannot compute probabilities exactly in undirected models either, due to the normalizing constant .

In this paper, we advocate a third approach to distribution estimation, based on autoregressive models and feed-forward neural networks. We refer to our particular approach as Neural Autoregressive Distribution Estimation (NADE). Its main distinguishing property is that computing

under a NADE model is tractable and can be computed efficiently, given an arbitrary ordering of the dimensions of

. We show that the framework is flexible and can model both binary and real-valued observations, can be made order-agnostic, and can be adapted to the case of 2D images using convolutional neural networks. In each case, we’re able to reach competitive results, compared to popular directed and undirected graphical model alternatives.

2 Nade

We consider the problem of modeling the distribution

of input vector observations

. For now, we will assume that the dimensions of are binary, that is . The model generalizes to other data types, which is explored later (Section 3) and in other work (Section 8).

NADE begins with the observation that any -dimensional distribution can be factored into a product of one-dimensional distributions, in any order (a permutation of the integers ):


Here contains the first dimensions in ordering and is the corresponding subvector for these dimensions. Thus, one can define an ‘autoregressive’ generative model of the data simply by specifying a parameterization of all conditionals .

Frey et al. (1996)

followed this approach and proposed using simple (log-)linear logistic regression models for these conditionals. This choice yields surprisingly competitive results, but are not competitive with non-linear models such as an RBM.

Bengio and Bengio (2000)

proposed a more flexible approach, with a single-layer feed-forward neural network for each conditional. Moreover, they allowed connections between the output of each network and the hidden layer of networks for the conditionals appearing earlier in the autoregressive ordering. Using neural networks led to some improvements in modeling performance, though at the cost of a really large model for very high-dimensional data.

In NADE, we also model each conditional using a feed-forward neural network. Specifically, each conditional is parameterized as follows:


where is the logistic sigmoid, and with as the number of hidden units, , , , are the parameters of the NADE model.

The hidden layer matrix and bias are shared by each hidden layer (which are all of the same size). This parameter sharing scheme (illustrated in Figure 1) means that NADE has parameters, rather than required if the neural networks were separate. Limiting the number of parameters can reduce the risk of over-fitting. Another advantage is that all hidden layers can be computed in time instead of . Denoting the pre-activation of the hidden layer as , this complexity is achieved by using the recurrence


where Equation 5 given vector can be computed in . Moreover, the computation of Equation 2 given is also . Thus, computing from conditional distributions (Equation 1) costs for NADE. This complexity is comparable to that of regular feed-forward neural network models.

Figure 1: Illustration of a NADE model. In this example, in the input layer, units with value 0 are shown in black while units with value 1 are shown in white. The dashed border represents a layer pre-activation.The outputs give predictive probabilities for each dimension of a vector , given elements earlier in some ordering. There is no path of connections between an output and the value being predicted, or elements of later in the ordering. Arrows connected together correspond to connections with shared (tied) parameters.

NADE can be trained by maximum likelihood, or equivalently by minimizing the average negative log-likelihood,


usually by stochastic (minibatch) gradient descent. As probabilities cost , gradients of the negative log-probability of training examples can also be computed in . Algorithm 1 describes the computation of both and the gradients of with respect to NADE’s parameters.

Input: training observation vector and ordering of the input dimensions.
Output: and gradients of on parameters.
# Computing
for  from 1 to  do
end for
# Computing gradients of
for  from to 1 do
end for
Algorithm 1 Computation of and learning gradients for NADE.

2.1 Relationship with the RBM

The proposed weight-tying for NADE isn’t simply motivated by computational reasons. It also reflects the computations of approximation inference in the RBM.

Denoting the energy function and distribution under an RBM as


computing all conditionals


is intractable. However, these could be approximated using mean-field variational inference. Specifically, consider the conditional over , and instead:


A mean-field approach could first approximate this conditional with a factorized distribution


where is the marginal probability of being equal to 1, given . Similarly, is the marginal for hidden variable . The dependence on comes from conditioning on , that is on the first dimensions of in ordering .

For some , a mean-field approximation is obtained by finding the parameters for and for which minimize the KL divergence between and . This is usually done by finding message passing updates that each set the derivatives of the KL divergence to 0 for some of the parameters of given others.

For some , let us fix for , leaving only for to be found. The KL-divergence develops as follows:

Then, we can take the derivative with respect to and set it to 0, to obtain:


where in the last step we have used the fact that for . Equation 14 would correspond to the message passing updates of the hidden unit marginals given the marginals of input .

Similarly, we can set the derivative with respect to for to 0 and obtain:


Equation 15 would correspond to the message passing updates of the input marginals given the hidden layer marginals . The complete mean-field algorithm would thus alternate between applying the updates of Equations 14 and 15, right to left.

We now notice that Equation 14 corresponds to NADE’s hidden layer computation (Equation 3) where . Also, Equation 15 corresponds to NADE’s output layer computation (Equation 2) where , and . Thus, in short, NADE’s forward pass is equivalent to applying a single pass of mean-field inference to approximate all the conditionals of an RBM, where initially and where a separate matrix is used for the hidden-to-input messages. A generalization of NADE based on this connection to mean field inference has been further explored by Raiko et al. (2014).

3 NADE for non-binary observations

So far we have only considered the case of binary observations . However, the framework of NADE naturally extends to distributions over other types of observations.

In the next section, we discuss the case of real-valued observations, which is one of the most general cases of non-binary observations and provides an illustrative example of the technical considerations one faces when extending NADE to new observations.

3.1 RNADE: Real-valued NADE

A NADE model for real-valued data could be obtained by applying the derivations shown in Section 2.1 to the Gaussian-RBM (Welling et al., 2005)

. The resulting neural network would output the mean of a Gaussian with fixed variance for each of the conditionals in Equation 

1. Such a model is not competitive with mixture models, for example on perceptual datasets (Uria, 2015). However, we can explore alternative models by making the neural network for each conditional distribution output the parameters of a distribution that’s not a fixed-variance Gaussian.

In particular, a mixture of one-dimensional Gaussians for each autoregressive conditional provides a flexible model. Given enough components, a mixture of Gaussians can model any continuous distribution to arbitrary precision. The resulting model can be interpreted as a sequence of mixture density networks (Bishop, 1994) with shared parameters. We call this model RNADE-MoG. In RNADE-MoG, each of the conditionals is modeled by a mixture of Gaussians:


where the parameters are set by the outputs of a neural network:


Parameter sharing conveys the same computational and statistical advantages as it does in the binary NADE.

Different one dimensional conditional forms may be preferred, for example due to limited dataset size or domain knowledge about the form of the conditional distributions. Other choices, like single variable-variance Gaussians, sinh-arcsinh distributions, and mixtures of Laplace distributions, have been examined by Uria (2015).

Training an RNADE can still be done by stochastic gradient descent on the parameters of the model with respect to the negative log-density of the training set. It was found empirically 

(Uria et al., 2013) that stochastic gradient descent leads to better parameter configurations when the gradient of the mean

was multiplied by the standard deviation (


4 Orderless and Deep NADE

The fixed ordering of the variables in a NADE model makes the exact calculation of arbitrary conditional probabilities computationally intractable. Only a small subset of conditional distributions, those where the conditioned variables are at the beginning of the ordering and marginalized variables at the end, are computationally tractable.

Another limitation of NADE is that a naive extension to a deep version, with multiple layers of hidden units, is computationally expensive. Deep neural networks (Bengio, 2009; LeCun et al., 2015) are at the core of state-of-the-art models for supervised tasks like image recognition (Krizhevsky et al., 2012) and speech recognition (Dahl et al., 2013). The same inductive bias should also provide better unsupervised models. However, extending the NADE framework to network architectures with several hidden layers, by introducing extra non-linear calculations between Equations (3) and (2), increases its complexity to cubic in the number of units per layer. Specifically, the cost becomes , where stands for the number of hidden layers and can be assumed to be a small constant, is the number of variables modeled, and is the number of hidden units, which we assumed to be of the same order as . This increase in complexity is caused by no longer being able to share hidden layer computations across the conditionals in Equation 1, after the non-linearity in the first layer.

In this section we introduce an order-agnostic training procedure, DeepNADE, which will address both of the issues above. This procedure trains a single deep neural network that can assign a conditional distribution to any variable given any subset of the others. This network can then provide the conditionals in Equation 1 for any ordering of the input observations. Therefore, the network defines a factorial number of different models with shared parameters, one for each of the orderings of the inputs. At test time, given an inference task, the most convenient ordering of variables can be used. The models for different orderings will not be consistent with each other: they will assign different probabilities to a given test vector. However, we can use the models’ differences to our advantage by creating ensembles of NADE models (Section 4.1), which results in better estimators than any single NADE. Moreover, the training complexity of our procedure increases linearly with the number of hidden layers , while remaining quadratic in the size of the network’s layers.

We first describe the model for an

-layer neural network modeling binary variables. A conditional distribution is obtained directly from a hidden unit in the final layer:


This hidden unit is computed from previous layers, all of which can only depend on the variables that are currently being conditioned on. We remove the other variables from the computation using a binary mask,


which is element-wise multiplied with the inputs before computing the remaining layers as in a standard neural network:


The network is specified by a free choice of the activation function

, and learnable parameters and , where is the number of units in the -th layer. As layer zero is the masked input, . The final -th layer needs to be able to provide predictions for any element (Equation 23) and so also has units.

To train a DeepNADE, the ordering of the variables is treated as a stochastic variable with a uniform distribution. Moreover, since we wish DeepNADE to provide good predictions for any ordering, we optimize the expected likelihood over the ordering of variables:


where we’ve made the dependence on the ordering and the network’s parameters explicit, stands for the set of all orderings (the permutations of elements) and is a uniformly sampled datapoint from the training set . Using NADE’s expression for the density of a datapoint in Equation 1 we have


where indexes the elements in the ordering, , of the variables. By moving the expectation over orderings inside the sum over the elements of the ordering, the ordering can be split in three parts: (the indices of the first dimensions in the ordering), (the index of the -th variable) and

(the indices of the remaining dimensions). Therefore, the loss function can be rewritten as:


The value of each of these terms does not depend on . Therefore, it can be simplified as:


In practice, this loss function will have a very high number of terms and will have to be approximated by sampling , and . The innermost expectation over values of can be calculated cheaply, because all of the neural network computations depend only on the masked input , and can be reused for each possible . Assuming all orderings are equally probable, we will estimate by:


which is an unbiased estimator of Equation 

29. Therefore, training can be done by descent on the gradient of .

For binary observations, we use the cross-entropy scaled by a factor of as the training loss which corresponds to minimizing :


Differentiating this cost involves backpropagating the gradients of the cross-entropy only from the outputs in

and rescaling them by .

The resulting training procedure resembles that of a denoising autoencoder

(Vincent et al., 2008)

. Like the autoencoder,

outputs are used to predict inputs corrupted by a random masking process ( in Equation 25). A single forward pass can compute , which provides a prediction for every masked variable, which could be used next in an ordering starting with . Unlike the autoencoder, the outputs for variables corresponding to those provided in the input (not masked out) are ignored.

In this order-agnostic framework, missing variables and zero-valued observations are indistinguishable by the network. This shortcoming can be alleviated by concatenating the inputs to the network (masked variables ) with the mask . Therefore we advise substituting the input described in Equation 25 with


We found this modification to be important in order to obtain competitive statistical performance (see Table 


). The resulting neural network is illustrated in Figure 2.

Figure 2: Illustration of a DeepNADE model with two hidden layers. The dashed border represents a layer pre-activation. A mask specifies a subset of variables to condition on. A conditional or predictive probability of the remaining variables is given in the final layer. Note that the output units with a corresponding input mask of value 1 (shown with dotted contour) are effectively not involved in DeepNADE’s training loss (Equation 34).

4.1 Ensembles of NADE models

As mentioned, the DeepNADE parameter fitting procedure effectively produces a factorial number of different NADE models, one for each ordering of the variables. These models will not, in general, assign the same probability to any particular datapoint. This disagreement is undesirable if we require consistent inferences for different inference problems, as it will preclude the use of the most convenient ordering of variables for each inference task.

However, it is possible to use this variability across the different orderings to our advantage by combining several models. A usual approach to improve on a particular estimator is to construct an ensemble of multiple, strong but different estimators, e.g. using bagging (Ormoneit and Tresp, 1995) or stacking (Smyth and Wolpert, 1999). The DeepNADE training procedure suggests a way of generating ensembles of NADE models: take a set of uniformly distributed orderings over the input variables and use the average probability as an estimator.

The use of an ensemble increases the test-time cost of density estimation linearly with the number of orderings used. The complexity of sampling does not change however: after one of the orderings is chosen at random, the single corresponding NADE is sampled. Importantly, the cost of training also remains the same, unlike other ensemble methods such as bagging. Furthermore, the number of components can be chosen after training and even adapted to a computational budget on the fly.

5 ConvNADE: Convolutional NADE

One drawback of NADE (and its variants so far) is the lack of a mechanism for truly exploiting the high-dimensional structure of the data. For example, when using NADE on binarized MNIST, we first need to flatten the 2D images before providing them to the model as a vector. As the spatial topology is not provided to the network, it can’t use this information to share parameters and may learn less quickly.

Recently, convolutional neural networks (CNN) have achieved state-of-the-art performance on many supervised tasks related to images Krizhevsky et al. (2012). Briefly, CNNs are composed of convolutional layers, each one having multiple learnable filters. The outputs of a convolutional layer are feature maps and are obtained by the convolution on the input image (or previous feature maps) of a linear filter, followed by the addition of a bias and the application of a non-linear activation function. Thanks to the convolution, spatial structure in the input is preserved and can be exploited. Moreover, as per the definition of a convolution the same filter is reused across all sub-regions of the entire image (or previous feature maps), yielding a parameter sharing that is natural and sensible for images.

The success of CNNs raises the question: can we exploit the spatial topology of the inputs while keeping NADE’s autoregressive property? It turns out we can, simply by replacing the fully connected hidden layers of a DeepNADE model with convolutional layers. We thus refer to this variant as Convolutional NADE (ConvNADE).

First we establish some notation that we will use throughout this section. Without loss of generality, let the input be a square binary image of size and every convolution filter connecting two feature maps and also be square with their size varying for each layer . We also define the following mask , which is 1 for the locations of the first pixels in the ordering .

Formally, Equation 26 is modified to use convolutions instead of dot products. Specifically for an -layer convolutional neural network that preserves the input shape (explained below) we have




where is the number of feature maps output by the -th layer and , , with denoting the element-wise multiplication, being any activation function and is the concatenation of every row in . Note that corresponds to the number of channels the input images have.

For notational convenience, we use to denote both “valid” convolutions and “full” convolutions, instead of introducing bulky notations to differentiate these cases. The “valid” convolutions only apply a filter to complete patches of the image, resulting in a smaller image (its shape is decreased to

). Alternatively, “full” convolutions zero-pad the contour of the image before applying the convolution, thus expanding the image (its shape is increased to

). Which one is used should be self-explanatory depending on the context. Note that we only use convolutions with a stride of 1.

Moreover, in order for ConvNADE to output conditional probabilities as shown in Equation 36, the output layer must have only one feature map , whose dimension matches the dimension of the input . This can be achieved by carefully combining layers that use either “valid” or “full” convolutions.

To explore different model architectures respecting that constraint, we opted for the following strategy. Given a network, we ensured the first half of its layers was using “valid” convolutions while the other half would use “full” convolutions. In addition to that, we made sure the network was symmetric with respect to its filter shapes (the filter shape used in layer matched the one used in layer ).

For completeness, we wish to mention that ConvNADE can also include pooling and upsampling layers, but we did not see much improvement when using them. In fact, recent research suggests that these types of layers are not essential to obtain state-of-the-art results (Springenberg et al., 2015).

The flexibility of DeepNADE allows us to easily combine both convolutional and fully connected layers. To create such hybrid models, we used the simple strategy of having two separate networks, with their last layer fused together at the end. The ‘convnet’ part is only composed of convolutional layers whereas the ‘fullnet’ part is only composed of fully connected layers. The forward pass of both networks follows respectively Equations (37)–(39) and Equations (25)–(27). Note that in the ‘fullnet’ network case, corresponds to the input image having been flattened.

In the end, the output layer of the hybrid model corresponds to the aggregation of the last layer pre-activation of both ‘convnet’ and ‘fullnet’ networks. The conditionals are slightly modified as follows:


The same training procedure as for DeepNADE model can also be used for ConvNADE. For binary observations, the training loss is similar to Equation 34, with being substituted for as defined in Equation 42.

As for the DeepNADE model, we found that providing the mask as an input to the model improves performance (see Table LABEL:tab:mnist-2d-results). For the ‘convnet’ part, the mask was provided as an additional channel to the input layer. For the ‘fullnet’ part, the inputs were concatenated with the mask as shown in Equation 35.

Figure 3: Illustration of a ConvNADE that combines a convolutional neural network with three hidden layers and a fully connected feed-forward neural network with two hidden layers. The dashed border represents a layer pre-activation. Units with a dotted contour are not valid conditionals since they depend on themselves they were given in the input.

The final architecture is shown in Figure 3. In our experiments, we found that this type of hybrid model works better than only using convolutional layers (see Table LABEL:tab:mnist-2d-results). Certainly, more complex architectures could be employed but this is a topic left for future work.

6 Related Work

As we mentioned earlier, the development of NADE and its extensions was motivated by the question of whether a tractable distribution estimator could be designed to match a powerful but intractable model such as the restricted Boltzmann machine.

The original inspiration came from the autoregressive approach taken by fully visible sigmoid belief networks (FVSBN), which were shown by Frey et al. (1996) to be surprisingly competitive, despite the simplicity of the distribution family for its conditionals. Bengio and Bengio (2000) later proposed using more powerful conditionals, modeled as single layer neural networks. Moreover, they proposed connecting the output of each conditional to all of the hidden layers of the neural networks for the preceding conditionals. More recently, Germain et al. (2015) generalized this model by deriving a simple procedure for making it deep and orderless (akin to DeepNADE, in Section 4). We compare with all of these approaches in Section 7.1.

There exists, of course, more classical and non-autoregressive approaches to tractable distribution estimation, such as mixture models and Chow–Liu trees (Chow and Liu, 1968). We compare with these as well in Section 7.1.

This work also relates directly to the recently growing literature on generative neural networks. In addition to the autoregressive approach described in this paper, there exists three other types of such models: directed generative networks, undirected generative networks and hybrid networks.

Work on directed generative networks dates back to the original work on sigmoid belief networks (Neal, 1992) and the Helmholtz machine (Hinton et al., 1995; Dayan et al., 1995). Helmholtz machines are equivalent to a multilayer sigmoid belief network, with each using binary stochastic units. Originally they were trained using Gibbs sampling and gradient descent (Neal, 1992), or with the so-called wake sleep algorithm (Hinton et al., 1995). More recently, many alternative directed models and training procedures have been proposed. Kingma and Welling (2014); Rezende et al. (2014) proposed the variational autoencoder (VAE), where the model is the same as the Helmholtz machine, but with real-valued (usually Gaussian) stochastic units. Importantly, Kingma and Welling (2014) identified a reparameterization trick making it possible to train the VAE in a way that resembles the training of an autoencoder. This approach falls in the family of stochastic variational inference methods, where the encoder network corresponds to the approximate variational posterior. The VAE optimizes a bound on the likelihood which is estimated using a single sample from the variational posterior, though recent work has shown that a better bound can be obtained using an importance sampling approach (Burda et al., 2016). Gregor et al. (2015) later exploited the VAE approach to develop DRAW, a directed generative model for images based on a read-write attentional mechanism. Goodfellow et al. (2014)

also proposed an adversarial approach to training directed generative networks, that relies on a discriminator network simultaneously trained to distinguish between data and model samples. Generative networks trained this way are referred to as Generative Adversarial Networks (GAN). While the VAE optimizes a bound of the likelihood (which is the KL divergence between the empirical and model distributions), it can be shown that GAN optimizes the Jensen–Shannon (JS) divergence between the empirical and model distributions.

Li et al. (2015) instead propose a training objective derived from Maximum Mean Discrepancy (MMD; Gretton et al., 2007). Recently, the directed generative model approach has been very successfully applied to model images (Denton et al., 2015; Sohl-Dickstein et al., 2011).

The undirected paradigm has also been explored extensively for developing powerful generative networks. These include the restricted Boltzmann machine (Smolensky, 1986; Freund and Haussler, 1992) and its multilayer extension, the deep Boltzmann machine (Salakhutdinov and Hinton, 2009), which dominate the literature on undirected neural networks. Salakhutdinov and Murray (2008) provided one of the first quantitative evidence of the generative modeling power of RBMs, which motivated the original parameterization for NADE (Larochelle and Murray, 2011)

. Efforts to train better undirected models can vary in nature. One has been to develop alternative objectives to maximum likelihood. The proposal of Contrastive Divergence

(CD; Hinton, 2002) was instrumental in the popularization of the RBM. Other proposals include pseudo-likelihood (Besag, 1975; Marlin et al., 2010), score matching (Hyvärinen, 2005, 2007a, 2007b)

, noise contrastive estimation 

(Gutmann and Hyvärinen, 2010) and probability flow minimization (Sohl-Dickstein et al., 2011). Another line of development has been to optimize likelihood using Robbins–Monro stochastic approximation (Younes, 1989), also known as Persistent CD (Tieleman, 2008), and develop good MCMC samplers for deep undirected models (Salakhutdinov, 2009, 2010; Desjardins et al., 2010; Cho et al., 2010). Work has also been directed towards proposing improved update rules or parameterization of the model’s energy function (Tieleman and Hinton, 2009; Cho et al., 2013; Montavon and Müller, 2012) as well as improved approximate inference of the hidden layers (Salakhutdinov and Larochelle, 2010). The work of Ngiam et al. (2011) also proposed an undirected model that distinguishes itself from deep Boltzmann machines by having deterministic hidden units, instead of stochastic.

Finally, hybrids of directed and undirected networks are also possible, though much less common. The most notable case is the Deep Belief Network

(DBN; Hinton et al., 2006), which corresponds to a sigmoid belief network for which the prior over its top hidden layer is an RBM (whose hidden layer counts as an additional hidden layer). The DBN revived interest in RBMs, as they were required to successfully initialize the DBN.

NADE thus substantially differs from this literature focusing on directed and undirected models, benefiting from a few properties that these approaches lack. Mainly, NADE does not rely on latent stochastic hidden units, making it possible to tractably compute its associated data likelihood for some given ordering. This in turn makes it possible to efficiently produce exact samples from the model (unlike in undirected models) and get an unbiased gradient for maximum likelihood training (unlike in directed graphical models).

7 Results

In this section, we evaluate the performance of our different NADE models on a variety of datasets.

7.1 Binary vectors datasets

We start by evaluating the performance of NADE models on a set of benchmark datasets where the observations correspond to binary vectors. These datasets were mostly taken from the LIBSVM datasets web site111, except for OCR-letters222 and NIPS-0-12 333 Code to download these datasets is available here: Table LABEL:tab:uci_datasets_info summarizes the main statistics for these datasets.

Name # Inputs Train Valid. Test
Adult 123 5000 1414 26147
Connect4 126 16000 4000 47557
DNA 180 1400 600 1186
Mushrooms 112 2000 500 5624
NIPS-0-12 500 400 100 1240
OCR-letters 128 32152 10000 10000
RCV1 150 40000 10000 150000
Web 300 14000 3188 32561
Table 1: Statistics on the binary vector datasets of Section 7.1.

For these experiments, we only consider tractable distribution estimators, where we can evaluate on test items exactly. We consider the following baselines:

  • MoB: A mixture of multivariate Bernoullis, trained using the EM algorithm. The number of mixture components was chosen from based on validation set performance, and early stopping was used to determine the number of EM iterations.

  • RBM: A restricted Boltzmann machine made tractable by using only 23 hidden units, trained by contrastive divergence with up to 25 steps of Gibbs sampling. The validation set performance was used to select the learning rate from , and the number of iterations over the training set from .

  • FVSBN: Fully visible sigmoid belief network, that models each conditional with logistic regression. The ordering of inputs was selected randomly. Training was by stochastic gradient descent. The validation set was used for early stopping, as well as for choosing the base learning rate , and a decreasing schedule constant from for the learning rate schedule for the update.

  • Chow–Liu: A Chow–Liu tree is a graph over the observed variables, where the distribution of each variable, except the root, depends on a single parent node. There is an fitting algorithm to find the maximum likelihood tree and conditional distributions (Chow and Liu, 1968). We adapted an implementation provided by Harmeling and Williams (2011), who found Chow–Liu to be a strong baseline.

    The maximum likelihood parameters are not defined when conditioning on events that haven’t occurred in the training set. Moreover, conditional probabilities of zero are possible, which could give infinitely bad test set performance. We re-estimated the conditional probabilities on the Chow–Liu tree using Lidstone or “add-” smoothing:


    selecting for each dataset from based on performance on the validation set.

  • MADE (Germain et al., 2015): Generalization of the neural network approach of Bengio and Bengio (2000), to multiple layers. We consider a version using a single (fixed) input ordering and another trained on multiple orderings from which an ensemble was constructed (which was inspired from the order-agnostic approach of Section 4) that we refer to as MADE-E. See Germain et al. (2015) for more details.

We compare these baselines with the two following NADE variants:

  • NADE (fixed order): Single layer NADE model, trained on a single (fixed) randomly generated order, as described in Section 2. The sigmoid activation function was used for the hidden layer, of size 500. Much like for FVSBN, training relied on stochastic gradient descent and the validation set was used for early stopping, as well as for choosing the learning rate from {0.05, 0.005, 0.0005}, and the decreasing schedule constant from {0,0.001,0.000001}.

  • NADE-E: Single layer NADE trained according to the order-agnostic procedure described in Section 4. The rectified linear activation function was used for the hidden layer, also of size 500. Minibatch gradient descent was used for training, with minibatches of size 100. The initial learning rate, chosen among {0.016, 0.004, 0.001, 0.00025, 0.0000675}, was linearly decayed to zero over the course of parameter updates. Early stopping was used, using Equation 34 to get a stochastic estimate of the validation set average log-likelihood. An ensemble using 16 orderings was used to compute the test-time log-likelihood.

Model Adult Connect4 DNA Mushrooms NIPS-0-12 OCR-letters RCV1 Web
MoB -20.44 -23.41 -98.19 -14.46 -290.02 -40.56 -47.59 -30.16
RBM -16.26 -22.66 -96.74 -15.15 -277.37 -43.05 -48.88 -29.38
FVSBN -13.17 -12.39 -83.64 -10.27 -276.88 -39.30 -49.84 -29.35
Chow–Liu -18.51 -20.57 -87.72 -20.99 -281.01 -48.87 -55.60 -33.92
MADE -13.12 -11.90 -83.63 -9.68 -280.25 -28.34 -47.10 -28.53
MADE-E -13.13 -11.90 -79.66 -9.69 -277.28 -30.04 -46.74 -28.25
NADE -13.19 -11.99 -84.81 -9.81 -273.08 -27.22 -46.66 -28.39
NADE-E -13.19 -12.58 -82.31 -9.69 -272.39 -27.32 -46.12 -27.87
Table 2:

Average log-likelihood performance of tractable distribution baselines and NADE models, on binary vector datasets. The best result is shown in bold, along with any other result with an overlapping confidence interval.

Table LABEL:tab:bin-results presents the results. We observe that NADE restricted to a fixed ordering of the inputs achieves very competitive performance compared to the baselines. However, the order-agnostic version of NADE is overall the best method, being among the top performing model for 5 datasets out of 8.

The performance of fixed-order NADE is surprisingly robust to variations of the chosen input ordering. The standard deviation on the average log-likelihood when varying the ordering was small: on Mushrooms, DNA and NIPS-0-12, we observed standard deviations of 0.045, 0.05 and 0.15, respectively. However, models with different orders can do well on different test examples, which explains why ensembling can still help.

7.2 Binary image dataset

We now consider the case of an image dataset, constructed by binarizing the MNIST digit dataset, as generated by Salakhutdinov and Murray (2008). This benchmark has been a popular choice for the evaluation of generative neural network models. Here, we investigate two questions:

  1. How does NADE compare to other intractable generative models?

  2. Does the use of a convolutional architecture improve the performance of NADE?

For these experiments, in addition to the baselines already described in Section 7.1, we consider the following:

  • DARN (Gregor et al., 2014): This deep generative autoencoder has two hidden layers, one deterministic and one with binary stochastic units. Both layers have 500 units (denoted as ). Adaptive weight noise (adaNoise) was either used or not to avoid the need for early stopping (Graves, 2011). Evaluation of exact test probabilities is intractable for large latent representations. Hence, Monte Carlo was used to approximate the expected description length, which corresponds to an upper bound on the negative log-likelihood.

  • DRAW (Gregor et al., 2015): Similar to a variational autoencoder where both the encoder and the decoder are LSTMs, guided (or not) by an attention mechanism. In this model, both LSTMs (encoder and decoder) are composed of 256 recurrent hidden units and always perform 64 timesteps. When the attention mechanism is enabled, patches ( pixels) are provided as inputs to the encoder instead of the whole image and the decoder also produces patches ( pixels) instead of a whole image.

  • Pixel RNN (Oord et al., 2016): NADE-like model for natural images that is based on convolutional and LSTM hidden units. This model has 7 hidden layers, each composed of 16 units. Oord et al. (2016) proposed a novel two-dimensional LSTM, named Diagonal BiLSTM, which is used in this model. Unlike our ConvNADE, the ordering is fixed before training and at test time, and corresponds to a scan of the image in a diagonal fashion starting from a corner at the top and reaching the opposite corner at the bottom.

Figure 4: (Left): samples from NADE trained on binarized MNIST. (Right): probabilities from which each pixel was sampled. Ancestral sampling was used with the same fixed ordering used during training.
procedure PRETRAIN()
     if  then
          TRAIN-ALL() Train for 20 iterations.
     end if
end procedure
Algorithm 2 Pretraining of a NADE with hidden layers on dataset .
(a) ConvNADE
(b) ConvNADE + DeepNADE
Figure 5: Network architectures for binarized MNIST. (a) ConvNADE with 8 convolutional layers (depicted in blue). The number of feature maps for a given layer is given by the number before the “@” symbol followed by the filter size and the type of convolution is specified in parentheses. (b) The same ConvNADE combined with a DeepNADE consisting of three fully-connected layers of respectively 500, 500 and 784 units.

We compare these baselines with some NADE variants. The performance of a basic (fixed-order, single hidden layer) NADE model is provided in Table LABEL:tab:deepnade-mnist-results and samples are illustrated in Figure 4. More importantly, we will focus on whether the following variants achieve better test set performance:

  • DeepNADE: Multiple layers (1hl, 2hl, 3hl or 4hl) trained according to the order-agnostic procedure described in Section 4. Information about which inputs are masked was either provided or not (no input masks) to the model. The rectified linear activation function was used for all hidden layers. Minibatch gradient descent was used for training, with minibatches of size 1000. Training consisted of 200 iterations of 1000 parameter updates. Each hidden layer was pretrained according to Algorithm 2. We report an average of the average test log-likelihoods over ten different random orderings.

  • EoNADE: This variant is similar to DeepNADE except for the log-likelihood on the test set, which is instead computed from an ensemble that averages predictive probabilities over 2 or 128 orderings. To clarify, the DeepNADE results report the typical performance of one ordering, by averaging results after taking the log, and so do not combine the predictions of the models like EoNADE does.

  • ConvNADE: Multiple convolutional layers trained according to the order-agnostic procedure described in Section 4. The exact architecture is shown in Figure 5(a). Information about which inputs are masked was either provided or not (no input masks). The rectified linear activation function was used for all hidden layers. The Adam optimizer (Kingma and Ba, 2015) was used with a learning rate of

    . Early stopping was used with a look ahead of 10 epochs, using Equation 

    34 to get a stochastic estimate of the validation set average log-likelihood. An ensemble using 128 orderings was used to compute the log-likelihood on the test set.

  • ConvNADE + DeepNADE: This variant is similar to ConvNADE except for the aggregation of a separate DeepNADE model at the end of the network. The exact architecture is shown in Figure 5(b). The training procedure is the same as with ConvNADE.

Table LABEL:tab:deepnade-mnist-results presents the results obtained by models ignorant of the 2D topology, such as the basic NADE model. Addressing the first question, we observe that the order-agnostic version of NADE with two hidden layers is competitive with intractable generative models. Moreover, examples of the ability of DeepNADE to solve inference tasks by marginalization and conditional sampling are shown in Figure 6.

MoBernoullis K=10 168.95
MoBernoullis K=500 137.64
Chow–Liu tree 134.99
MADE 2hl (32 masks) 86.64
RBM (500 h, 25 CD steps) 86.34
DBN 2hl 84.55
DARN 84.71
DARN (adaNoise) 84.13
NADE (fixed order) 88.33
DeepNADE 1hl (no input masks) 99.37
DeepNADE 2hl (no input masks) 95.33
DeepNADE 1hl 92.17
DeepNADE 2hl 89.17
DeepNADE 3hl 89.38
DeepNADE 4hl 89.60
EoNADE 1hl (2 orderings) 90.69
EoNADE 1hl (128 orderings) 87.71
EoNADE 2hl (2 orderings) 87.96
EoNADE 2hl (128 orderings) 85.10
Table 3: Negative log-likelihood test results of models ignorant of the 2D topology on the binarized MNIST dataset.

Figure 6: Example of marginalization and sampling. The first column shows five examples from the test set of the MNIST dataset. The second column shows the density of these examples when a random 1010 pixel region is marginalized. The right-most five columns show samples for the hollowed region. Both tasks can be done easily with a NADE where the pixels to marginalize are at the end of the ordering.

Now, addressing the second question, we can see from Table LABEL:tab:mnist-2d-results that convolutions do improve the performance of NADE. Moreover, we observe that providing information about which inputs are masked is essential to obtaining good results. We can also see that combining convolutional and fully-connected layers helps. Even though ConvNADE+DeepNADE performs slightly worst than Pixel RNN, we note that our proposed approach is order-agnostic, whereas Pixel RNN requires a fixed ordering. Figure 7 shows samples obtained from the ConvNADE+DeepNADE model using ancestral sampling on a random ordering.

DRAW (without attention) 87.40
DRAW 80.97
Pixel RNN 79.20
ConvNADE+DeepNADE (no input masks) 85.25
ConvNADE 81.30
ConvNADE+DeepNADE 80.82
Table 4: Negative log-likelihood test results of models exploiting 2D topology on the binarized MNIST dataset.
Figure 7: (Left): samples from ConvNADE+DeepNADE trained on binarized MNIST. (Right): probabilities from which each pixel was sampled. Ancestral sampling was used with a different random ordering for each sample.

7.3 Real-valued observations datasets

In this section, we compare the statistical performance of RNADE to mixtures of Gaussians (MoG) and factor analyzers (MFA), which are surprisingly strong baselines in some tasks (Tang et al., 2012; Zoran and Weiss, 2012).

7.3.1 Low-dimensional data

We start by considering three UCI datasets (Bache and Lichman, 2013), previously used to study the performance of other density estimators (Silva et al., 2011; Tang et al., 2012), namely: red wine, white wine and parkinsons. These are low dimensional datasets (see Table 5) with hard thresholds and non-linear dependencies that make it difficult to fit mixtures of Gaussians or factor analyzers.

Following Tang et al. (2012), we eliminated discrete-valued attributes and an attribute from every pair with a Pearson correlation coefficient greater than 0.98. We normalized each dimension of the data by subtracting its training-subset sample mean and dividing by its standard deviation. All results are reported on the normalized data.

We use full-covariance Gaussians and mixtures of factor analysers as baselines. Models were compared on their log-likelihood on held-out test data. Due to the small size of the datasets (see Table 5), we used 10-folds, using 90% of the data for training, and 10% for testing.

We chose the hyperparameter values for each model by doing per-fold cross-validation, using a ninth of the training data as validation data. Once the hyperparameter values have been chosen, we train each model using all the training data (including the validation data) and measure its performance on the 10% of held-out testing data. In order to avoid overfitting, we stopped the training after reaching a training likelihood higher than the one obtained on the best validation-wise iteration of the best validation run. Early stopping was important to avoid overfitting the RNADE models. It also improved the results of the MFAs, but to a lesser degree.

The MFA models were trained using the EM algorithm (Ghahramani and Hinton, 1996; Verbeek, 2005). We cross-validated the number of components and factors. We also selected the number of factors from , where choosing results in a mixture of Gaussians, and the number of components was chosen among . Cross-validation selected fewer than components in every case.

We report the performance of several RNADE models using different parametric forms for the one-dimensional conditionals: Gaussian with fixed variance (RNADE-FV), Gaussian with variable variance (RNADE-Gaussian), - distribution (RNADE-SAS), mixture of Gaussians (RNADE-MoG), and mixture of Laplace distributions (RNADE-MoL). All RNADE models were trained by stochastic gradient descent, using minibatches of size 100, for 500 epochs, each epoch comprising 10 minibatches. We fixed the number of hidden units to

, and the non-linear activation function of the hidden units to ReLU. Three hyperparameters were cross-validated using grid-search: the number of components on each one-dimensional conditional (only applicable to the RNADE-MoG and RNADE-MoL models) was chosen from

, the weight-decay (used only to regularize the input to hidden weights) from , and the learning rate from . Learning rates were decreased linearly to reach 0 after the last epoch.

Red wine White wine Parkinsons
Total number of datapoints
Table 5: Dimensionality and size of the UCI datasets used in Section 7.3.1
Model Red wine White wine Parkinsons
Table 6: Average test set log-likelihoods per datapoint for seven models on three UCI datasets. Performances not in bold can be shown to be significantly worse than at least one of the results in bold as per a paired -test on the ten mean-likelihoods (obtained from each data fold), with significance level 0.05.

The results are shown in Table 6. RNADE with mixture of Gaussian conditionals was among the statistically significant group of best models on all datasets. As shown in Figure 8

, RNADE-SAS and RNADE-MoG models are able to capture hard thresholds and heteroscedasticity.

Figure 8: Scatter plot of dimensions vs of the red wine dataset. A thousand datapoints from the dataset are shown in black in all subfigures. As can be observed, this conditional distribution

is heteroscedastic, skewed and has hard thresholds. In red, a thousand samples from four RNADE models with different one-dimensional conditional forms are shown.

Top-left: In red, one thousand samples from a RNADE-FV model. Top-right: In red, one thousand samples from a RNADE-Gaussian model. Bottom-left: In red, one thousand samples from a RNADE-SAS (sinh-arcsinh distribution) model. Bottom-right: In red, one thousand samples from a RNADE-MoG model with 20 components per one-dimensional conditional. The RNADE-SAS and RNADE-MoG models successfully capture all the characteristics of the data.

7.3.2 Natural image patches

We also measured the ability of RNADE to model small patches of natural images. Following the work of Zoran and Weiss (2011), we use 8-by-8-pixel patches of monochrome natural images, obtained from the BSDS300 dataset Martin et al., 2001; Figure 9 gives examples.

Pixels in this dataset can take a finite number of brightness values ranging from to . We added uniformly distributed noise between and  to the brightness of each pixel. We then divided by , making the pixels take continuous values in the range . Adding noise prevents deceivingly high-likelihood solutions that assign narrow high-density spikes around some of the possible discrete values.

We subtracted the mean pixel value from each patch. Effectively reducing the dimensionality of the data. Therefore we discarded the 64th (bottom-right) pixel, which would be perfectly predictable and models could fit arbitrarily high densities to it. All of the results in this section were obtained by fitting the pixels in a raster-scan order.

Experimental details follow. We trained our models by using patches randomly drawn from 180 images in the training subset of BSDS300. We used the remaining 20 images in the training subset as validation data. We used 1000 random patches from the validation subset to early-stop training of RNADE. We measured the performance of each model by their log-likelihood on one million patches drawn randomly from the test subset of 100 images not present in the training data. Given the larger scale of this dataset, hyperparameters of the RNADE and MoG models were chosen manually using the performance of preliminary runs on the validation data, rather than by grid search.

All RNADE models reported use ReLU activations for the hidden units. The RNADE models were trained by stochastic gradient descent, using 25 datapoints per minibatch, for a total of 1,000 epochs, each comprising 1,000 minibatches. The learning rate was initialized to 0.001, and linearly decreased to reach 0 after the last epoch. Gradient momentum with factor 0.9 was used, but initiated after the first epoch. A weight decay rate of 0.001 was applied to the input-to-hidden weight matrix only. We found that multiplying the gradient of the mean output parameters by the standard deviation improves results of the models with mixture outputs444Empirically, we found this to work better than regular gradients and also better than multiplying by the variances, which would provide a step with the right units.. RNADE training was early stopped but didn’t show signs of overfitting. Even larger models might perform better.

The MoG models were trained using 1,000 iterations of minibatch EM. At each iteration 20,000 randomly sampled datapoints were used in an EM update. A step was taken from the previous parameters’ value towards the parameters resulting from the M-step: . The step size, , was scheduled to start at 0.1 and linearly decreased to reach 0 after the last update. The training of the MoG was early-stopped and also showed no signs of overfitting.

The results are shown in Table 7. We report the average log-likelihood of each model for a million image patches from the test set. The ranking of RNADE models is maintained when ordered by validation likelihood: the model with best test-likelihood would have been chosen using crossvalidation across all the RNADE models shown in the table. We also compared RNADE with a MoG trained by Zoran and Weiss (downloaded from Daniel Zoran’s website) from which we removed the 64th row and column of each covariance matrix. There are two differences in the set-up of our experiments and those of Zoran and Weiss. First, we learned the means of the MoG components, while Zoran and Weiss (2011) fixed them to zero. Second, we held-out 20 images from the training set to do early-stopping and hyperparameter optimisation, while they used the 200 images for training.

Model Test log-likelihood
MoG (Zoran and Weiss, 2012)555This model was trained using the full 200 images in the BSDS training dataset, the rest of the models were trained using 180, reserving 20 for hyperparameter crossvalidation and early-stopping. 152.8
MoG 144.7