Neural Autoregressive Distribution Estimation

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 $\boldsymbol{h}$ from some prior $p(\boldsymbol{h})$, followed by sampling the observed data $\boldsymbol{x}$ from some conditional $p(\boldsymbol{x}\,|\,\boldsymbol{h})$. Unfortunately, this approach quickly becomes intractable and requires approximations when the latent state $\boldsymbol{h}$ increases in complexity. Specifically, computing the marginal probability of the data, $p(\boldsymbol{x})=\sum_{\boldsymbol{h}}p(\boldsymbol{x}\,|\,\boldsymbol{h})\,p(\boldsymbol{h})$, is only tractable under fairly constraining assumptions on $p(\boldsymbol{x}\,|\,\boldsymbol{h})$ and $p(\boldsymbol{h})$.

Another popular approach, based on undirected graphical models, gives probabilities of the form $p(\boldsymbol{x})=\operatorname{exp}\left\{\phi(\boldsymbol{x})\right\}/ Z$, where $\phi$ is a tractable function and $Z$ 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 $p(\boldsymbol{x})$ exactly in undirected models either, due to the normalizing constant $Z$.

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 $p(\boldsymbol{x})$ under a NADE model is tractable and can be computed efficiently, given an arbitrary ordering of the dimensions of $\boldsymbol{x}$. 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.

NADE

We consider the problem of modeling the distribution $p(\boldsymbol{x})$ of input vector observations $\boldsymbol{x}$. For now, we will assume that the dimensions of $\boldsymbol{x}$ are binary, that is $x_{d}\in\{0,1\}~{}\forall d$. 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 $D$-dimensional distribution $p(\boldsymbol{x})$ can be factored into a product of one-dimensional distributions, in any order $o$ (a permutation of the integers $1,\dots,D$):

Here $o_{<d}$ contains the first $d-1$ dimensions in ordering $o$ and $\boldsymbol{x}_{o_{<d}}$ 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 $D$ conditionals $p(x_{o_{d}}\,|\,\boldsymbol{x}_{o_{<d}})$.

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 $p(x_{o_{d}}\,|\,\boldsymbol{x}_{<d})$ is parameterized as follows:

The hidden layer matrix $\boldsymbol{W}$ and bias $\boldsymbol{c}$ are shared by each hidden layer $\boldsymbol{h}_{d}$ (which are all of the same size). This parameter sharing scheme (illustrated in Figure 1) means that NADE has $O(HD)$ parameters, rather than $O(HD^{2})$ required if the neural networks were separate. Limiting the number of parameters can reduce the risk of over-fitting. Another advantage is that all $D$ hidden layers $\boldsymbol{h}_{d}$ can be computed in $O(HD)$ time instead of $O(HD^{2})$. Denoting the pre-activation of the $d^{\rm th}$ hidden layer as $\boldsymbol{a}_{d}=\boldsymbol{W}_{\cdot,o_{<d}}\boldsymbol{x}_{o_{<d}}+\boldsymbol{c}$, this complexity is achieved by using the recurrence

where Equation 5 given vector $\boldsymbol{a}_{d-1}$ can be computed in $O(H)$. Moreover, the computation of Equation 2 given $\boldsymbol{h}$ is also $O(H)$. Thus, computing $p(\boldsymbol{x})$ from $D$ conditional distributions (Equation 1) costs $O(HD)$ for NADE. This complexity is comparable to that of regular feed-forward neural network models.

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

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

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

is intractable. However, these could be approximated using mean-field variational inference. Specifically, consider the conditional over $x_{o_{d}}$, $\boldsymbol{x}_{o_{>d}}$ and $\boldsymbol{h}$ instead:

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

where $\mu_{j}(d)$ is the marginal probability of $x_{o_{j}}$ being equal to 1, given $\boldsymbol{x}_{o_{<d}}$. Similarly, $\tau_{k}(d)$ is the marginal for hidden variable $h_{k}$. The dependence on $d$ comes from conditioning on $\boldsymbol{x}_{o_{<d}}$, that is on the first $d\!-\!1$ dimensions of $\boldsymbol{x}$ in ordering $o$.

For some $d$, a mean-field approximation is obtained by finding the parameters $\mu_{j}(d)$ for $j\in\{d,\dots,D\}$ and $\tau_{k}(d)$ for $k\in\{1,\dots,H\}$ which minimize the KL divergence between $q(x_{o_{d}},\boldsymbol{x}_{o_{>d}},\boldsymbol{h}\,|\,\boldsymbol{x}_{o_{<d}})$ and $p(x_{o_{d}},\boldsymbol{x}_{o_{>d}},\boldsymbol{h}\,|\,\boldsymbol{x}_{o_{<d}})$. 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 $q(x_{o_{d}},\boldsymbol{x}_{o_{>d}},\boldsymbol{h}\,|\,\boldsymbol{x}_{o_{<d}})$ given others.

For some $d$, let us fix $\mu_{j}(d)=x_{o_{d}}$ for $j<d$, leaving only $\mu_{j}(d)$ for $j>d$ to be found. The KL-divergence develops as follows:

Then, we can take the derivative with respect to $\tau_{k}(d)$ and set it to 0, to obtain:

where in the last step we have used the fact that $\mu_{j}(d)=x_{o_{j}}$ for $j<d$. Equation 14 would correspond to the message passing updates of the hidden unit marginals $\tau_{k}(d)$ given the marginals of input $\mu_{j}(d)$.

Similarly, we can set the derivative with respect to $\mu_{j}(d)$ for $j\geq d$ to 0 and obtain:

Equation 15 would correspond to the message passing updates of the input marginals $\mu_{j}(d)$ given the hidden layer marginals $\tau_{k}(d)$. 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 $\mu_{j}(d)=0~{}\;\forall j\geq d$. Also, Equation 15 corresponds to NADE’s output layer computation (Equation 2) where $j=d$, $\tau_{k}(d)=h_{d,k}$ and $\boldsymbol{W}^{\top}=\boldsymbol{V}$. Thus, in short, NADE’s forward pass is equivalent to applying a single pass of mean-field inference to approximate all the conditionals $p(\boldsymbol{x}_{o_{d}}\,|\,\boldsymbol{x}_{o_{<d}})$ of an RBM, where initially $\mu_{j}(d)=0$ and where a separate matrix $\boldsymbol{V}$ 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).

NADE for non-binary observations

So far we have only considered the case of binary observations $x_{i}$. 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.

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 $\left(\frac{\partial J}{\partial\mu_{o_{d},c}}\right)$ was multiplied by the standard deviation ($\sigma_{o_{d},c}$).

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 $O(DH^{2}L)$, where $L$ stands for the number of hidden layers and can be assumed to be a small constant, $D$ is the number of variables modeled, and $H$ is the number of hidden units, which we assumed to be of the same order as $D$. 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 $D!$ 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 $O(H^{2}L)$, while remaining quadratic in the size of the network’s layers.

We first describe the model for an $L$-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 $\boldsymbol{x}_{o_{<d}}$ 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:

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 $o$ and the network’s parameters $\boldsymbol{\theta}$ explicit, $D!$ stands for the set of all orderings (the permutations of $D$ elements) and $\boldsymbol{x}$ is a uniformly sampled datapoint from the training set $\cal{X}$. Using NADE’s expression for the density of a datapoint in Equation 1 we have

where $d$ indexes the elements in the ordering, $o$, 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: $o_{<d}$ (the indices of the $d-1$ first dimensions in the ordering), $o_{d}$ (the index of the $d$-th variable) and $o_{>d}$ (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 $o_{>d}$. 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 $\boldsymbol{x}$, $d$ and $o_{<d}$. The innermost expectation over values of $o_{d}$ can be calculated cheaply, because all of the neural network computations depend only on the masked input $\boldsymbol{x}_{o_{<d}}$, and can be reused for each possible $o_{d}$. Assuming all orderings are equally probable, we will estimate $\mathcal{J}(\boldsymbol{\theta})$ by:

which is an unbiased estimator of Equation 29. Therefore, training can be done by descent on the gradient of $\widehat{\mathcal{J}}(\boldsymbol{\theta})$.

For binary observations, we use the cross-entropy scaled by a factor of $\frac{D}{D-d+1}$ as the training loss which corresponds to minimizing $\widehat{\mathcal{J}}$:

Differentiating this cost involves backpropagating the gradients of the cross-entropy only from the outputs in $o_{\geq d}$ and rescaling them by $\frac{D}{D-d+1}$.

The resulting training procedure resembles that of a denoising autoencoder (Vincent et al., 2008). Like the autoencoder, $D$ outputs are used to predict $D$ inputs corrupted by a random masking process ($\boldsymbol{m}_{o_{<d}}$ in Equation 25). A single forward pass can compute $\boldsymbol{h}^{(L)}_{o_{\geq d}}$, which provides a prediction $p(x_{o_{d}}\!=\!1\,|\,\boldsymbol{x}_{o_{<d}},\boldsymbol{\theta},o_{<d},o_{d})$ for every masked variable, which could be used next in an ordering starting with $o_{<d}$. 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 $\boldsymbol{x}\odot\boldsymbol{m}_{o_{<d}}$) with the mask $\boldsymbol{m}_{o_{<d}}$. 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 LABEL:tab:deepnade-mnist-results). The resulting neural network is illustrated in Figure 2.

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 $\{o^{(k)}\}_{k=1}^{K}$ over the input variables and use the average probability $\frac{1}{K}\sum_{k=1}^{K}p(\boldsymbol{x}\,|\,\theta,o^{(k)})$ 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 $K$ 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.

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).

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

Moreover, in order for ConvNADE to output conditional probabilities as shown in Equation 36, the output layer must have only one feature map $\boldsymbol{H}^{(L)}_{1}$, whose dimension matches the dimension of the input $\boldsymbol{X}$. This can be achieved by carefully combining layers that use either “valid” or “full” convolutions.

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, $\boldsymbol{x}$ corresponds to the input image having been flattened.

In the end, the output layer $\boldsymbol{g}$ 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 $\boldsymbol{h}^{(L)}$ being substituted for $\boldsymbol{g}$ as defined in Equation 42.

As for the DeepNADE model, we found that providing the mask $\boldsymbol{M}_{o_{<d}}$ 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.

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.

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 $d^{\rm th}$ conditional to all of the hidden layers of the $d-1$ 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).

Results

In this section, we evaluate the performance of our different NADE models on a variety of 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 sitehttp://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/, except for OCR-lettershttp://ai.stanford.edu/~btaskar/ocr/ and NIPS-0-12 http://www.cs.nyu.edu/~roweis/data.html. Code to download these datasets is available here: http://info.usherbrooke.ca/hlarochelle/code/nade.tar.gz. Table LABEL:tab:uci_datasets_info summarizes the main statistics for these datasets.

For these experiments, we only consider tractable distribution estimators, where we can evaluate $p(\boldsymbol{x})$ 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 $\{32,64,128,256,512,1024\}$ 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 $\{0.005,0.0005,0.00005\}$, and the number of iterations over the training set from $\{100,500,1000\}$.

FVSBN: Fully visible sigmoid belief network, that models each conditional $p(x_{o_{d}}\,|\,\boldsymbol{x}_{o_{<d}})$ 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 $\eta\in\{0.05,0.005,0.0005\}$, and a decreasing schedule constant $\gamma$ from $\{0,0.001,0.000001\}$ for the learning rate schedule $\eta/(1+\gamma t)$ for the $t^{\rm th}$ 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 $O(D^{2})$ 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-$\alpha$” smoothing:

selecting $\alpha$ for each dataset from $\{10^{-20},0.001,0.01,0.1\}$ 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 $\gamma$ 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 $100,000$ 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.

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.

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:

How does NADE compare to other intractable generative models?

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 $n_{h}=500$). 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 ($2\times 2$ pixels) are provided as inputs to the encoder instead of the whole image and the decoder also produces patches ($5\times 5$ 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.

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 $10^{-4}$. 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.

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.

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).

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 $2,4,\dots D$, where choosing $D$ results in a mixture of Gaussians, and the number of components was chosen among $2,4,\dots 50$. Cross-validation selected fewer than $50$ 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), $sinh$-$arcsinh$ 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 $50$, 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 $\{2,5,10,20\}$, the weight-decay (used only to regularize the input to hidden weights) from $\left\{2.0,1.0,0.1,0.01,0.001,0\right\}$, and the learning rate from $\left\{0.1,0.05,0.025,0.0125\right\}$. Learning rates were decreased linearly to reach 0 after the last epoch.

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.

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 $255$. We added uniformly distributed noise between and $1$ to the brightness of each pixel. We then divided by $256$, 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 outputsEmpirically, 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: $\boldsymbol{\theta_{t}}=(1-\eta)\boldsymbol{\theta_{t-1}}+\eta\boldsymbol{\theta_{EM}}$. The step size, $\eta$, 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.

The RNADE-FV model with fixed conditional variances obtained very low statistical performance. Adding an output parameter per dimension to have variable standard deviations made our models competitive with MoG with 100 full-covariance components. However, in order to obtain results superior to the mixture of Gaussians model trained by Zoran and Weiss, we had to use richer conditional distributions: one-dimensional mixtures of Gaussians (RNADE-MoG). On average, the best RNADE model obtained 3.3 nats per patch higher log-density than a MoG fitted with the same training data.

In Figure 9, we show one hundred examples from the test set, one hundred examples from Zoran and Weiss’ mixture of Gaussians, and a hundred samples from our best RNADE-MoG model. Similar patterns can be observed in the three cases: uniform patches, edges, and locally smooth noisy patches.

3.3 Speech acoustics

We also measured the ability of RNADE to model small patches of speech spectrograms, extracted from the TIMIT dataset (Garofolo et al., 1993). The patches contained 11 frames of 20 filter-banks plus energy; totalling 231 dimensions per datapoint. A good generative model of speech acoustics could be used, for example, in denoising, or speech detection tasks.

We fitted the models using the standard TIMIT training subset, which includes recordings from 605 speakers of American English. We compare RNADE with a mixture of Gaussians by measuring their log-likelihood on the complete TIMIT core-test dataset: a held-out set of 25 speakers.

The RNADE models have 512 hidden units, ReLU activations, and a mixture of 20 one-dimensional Gaussian components per output. Given the large scale of this dataset, hyperparameter choices were again made manually using validation data. The same minibatch training procedures for RNADE and mixture of Gaussians were used as for natural image patches.

The RNADE models were trained by stochastic gradient descent, with 25 datapoints per minibatch, for a total of 200 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 momentum 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. Again, we found that multiplying the gradient of the mean output parameters by the standard deviation improved results. RNADE training was early stopped but didn’t show signs of overfitting.

As for the MoG model, it was trained exactly as in Section 7.3.2.

The results are shown in Table 8. The best RNADE (which would have been selected based on validation results) has 15 nats higher likelihood per test example than the best mixture of Gaussians. Examples from the test set, and samples from the MoG and RNADE-MoG models are shown in Figure 10. In contrast with the log-likelihood measure, there are no marked differences between the samples from each model. Both sets of samples look like blurred spectrograms, but RNADE seems to capture sharper formant structures (peaks of energy at the lower frequency bands characteristic of vowel sounds).

Conclusion

We’ve described the Neural Autoregressive Distribution Estimator, a tractable, flexible and competitive alternative to directed and undirected graphical models for unsupervised distribution estimation.

Since the publication of the first formulation of NADE (Larochelle and Murray, 2011), it has been extended to many more settings, other than those described in this paper. Larochelle and Lauly (2012); Zheng et al. (2015b) adapted NADE for topic modeling of documents and images, while Boulanger-Lewandowski et al. (2012) used NADE for modeling music sequential data. Theis and Bethge (2015) and Oord et al. (2016) proposed different NADE models for images than the one we presented, applied to natural images and based on convolutional and LSTM hidden units. Zheng et al. (2015a) used a NADE model to integrate an attention mechanism into an image classifier. Bornschein and Bengio (2015) showed that NADE could serve as a powerful prior over the latent state of directed graphical model. These are just a few examples of many possible ways one can leverage the flexibility and effectiveness of NADE models.

References