Masked Autoregressive Flow for Density Estimation

05/19/2017 ∙ by George Papamakarios, et al. ∙ 0

Autoregressive models are among the best performing neural density estimators. We describe an approach for increasing the flexibility of an autoregressive model, based on modelling the random numbers that the model uses internally when generating data. By constructing a stack of autoregressive models, each modelling the random numbers of the next model in the stack, we obtain a type of normalizing flow suitable for density estimation, which we call Masked Autoregressive Flow. This type of flow is closely related to Inverse Autoregressive Flow and is a generalization of Real NVP. Masked Autoregressive Flow achieves state-of-the-art performance in a range of general-purpose density estimation tasks.



There are no comments yet.


page 14

page 15

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

The joint density of a set of variables

is a central object of interest in machine learning. Being able to access and manipulate

enables a wide range of tasks to be performed, such as inference, prediction, data completion and data generation. As such, the problem of estimating from a set of examples

is at the core of probabilistic unsupervised learning and generative modelling.

In recent years, using neural networks for density estimation has been particularly successful. Combining the flexibility and learning capacity of neural networks with prior knowledge about the structure of data to be modelled has led to impressive results in modelling natural images

[Theis and Bethge, 2015, van den Oord et al., 2016c, b, Dinh et al., 2017, Salimans et al., 2017] and audio data [Uria et al., 2015, van den Oord et al., 2016a]. State-of-the-art neural density estimators have also been used for likelihood-free inference from simulated data [Paige and Wood, 2016, Papamakarios and Murray, 2016], variational inference [Rezende and Mohamed, 2015, Kingma et al., 2016], and as surrogates for maximum entropy models [Loaiza-Ganem et al., 2017].

Neural density estimators differ from other approaches to generative modelling—such as variational autoencoders

[Kingma and Welling, 2014, Rezende et al., 2014] and generative adversarial networks [Goodfellow et al., 2014]

—in that they readily provide exact density evaluations. As such, they are more suitable in applications where the focus is on explicitly evaluating densities, rather than generating synthetic data. For instance, density estimators can learn suitable priors for data from large unlabelled datasets, for use in standard Bayesian inference

[Zoran and Weiss, 2011]. In simulation-based likelihood-free inference, conditional density estimators can learn models for the likelihood [Fan et al., 2013] or the posterior Papamakarios and Murray [2016] from simulated data. Density estimators can learn effective proposals for importance sampling [Papamakarios and Murray, 2015] or sequential Monte Carlo [Gu et al., 2015, Paige and Wood, 2016]; such proposals can be used in probabilistic programming environments to speed up inference [Kulkarni et al., 2015, Le et al., 2017]. Finally, conditional density estimators can be used as flexible inference networks for amortized variational inference and as part of variational autoencoders [Kingma and Welling, 2014, Rezende et al., 2014].

A challenge in neural density estimation is to construct models that are flexible enough to represent complex densities, but have tractable density functions and learning algorithms. There are mainly two families of neural density estimators that are both flexible and tractable: autoregressive models [Uria et al., 2016] and normalizing flows [Rezende and Mohamed, 2015]. Autoregressive models decompose the joint density as a product of conditionals, and model each conditional in turn. Normalizing flows transform a base density (e.g. a standard Gaussian) into the target density by an invertible transformation with tractable Jacobian.

Our starting point is the realization (as pointed out by Kingma et al. [2016]) that autoregressive models, when used to generate data, correspond to a differentiable transformation of an external source of randomness (typically obtained by random number generators). This transformation has a tractable Jacobian by design, and for certain autoregressive models it is also invertible, hence it precisely corresponds to a normalizing flow. Viewing an autoregressive model as a normalizing flow opens the possibility of increasing its flexibility by stacking multiple models of the same type, by having each model provide the source of randomness for the next model in the stack. The resulting stack of models is a normalizing flow that is more flexible than the original model, and that remains tractable.

In this paper we present Masked Autoregressive Flow (MAF), which is a particular implementation of the above normalizing flow that uses the Masked Autoencoder for Distribution Estimation (MADE) [Germain et al., 2015] as a building block. The use of MADE enables density evaluations without the sequential loop that is typical of autoregressive models, and thus makes MAF fast to evaluate and train on parallel computing architectures such as Graphics Processing Units (GPUs). We show a close theoretical connection between MAF and Inverse Autoregressive Flow (IAF) [Kingma et al., 2016], which has been designed for variational inference instead of density estimation, and show that both correspond to generalizations of the successful Real NVP [Dinh et al., 2017]. We experimentally evaluate MAF on a wide range of datasets, and we demonstrate that MAF outperforms RealNVP and achieves state-of-the-art performance on a variety of general-purpose density estimation tasks.

2 Background

2.1 Autoregressive density estimation

Using the chain rule of probability, any joint density

can be decomposed into a product of one-dimensional conditionals as . Autoregressive density estimators [Uria et al., 2016] model each conditional as a parametric density, whose parameters are a function of a hidden state . In recurrent architectures, is a function of the previous hidden state and the input variable . The Real-valued Neural Autoregressive Density Estimator (RNADE) [Uria et al., 2013]

uses mixtures of Gaussian or Laplace densities for modelling the conditionals, and a simple linear rule for updating the hidden state. More flexible approaches for updating the hidden state are based on Long Short-Term Memory recurrent neural networks

[Theis and Bethge, 2015, van den Oord et al., 2016c].

A drawback of autoregressive models is that they are sensitive to the order of the variables. For example, the order of the variables matters when learning the density of Figure (a)a if we assume a model with Gaussian conditionals. As Figure (b)b shows, a model with order cannot learn this density, even though the same model with order can represent it perfectly. In practice is it hard to know which of the factorially many orders is the most suitable for the task at hand. Autoregressive models that are trained to work with an order chosen at random have been developed, and the predictions from different orders can then be combined in an ensemble [Uria et al., 2014, Germain et al., 2015]. Our approach (Section 3) can use a different order in each layer, and using random orders would also be possible.

Straightforward recurrent autoregressive models would update a hidden state sequentially for every variable, requiring sequential computations to compute the probability of a

-dimensional vector, which is not well-suited for computation on parallel architectures such as GPUs. One way to enable parallel computation is to start with a fully-connected model with

inputs and outputs, and drop out connections in order to ensure that output will only be connected to inputs . Output can then be interpreted as computing the parameters of the conditional . By construction, the resulting model will satisfy the autoregressive property, and at the same time it will be able to calculate efficiently on a GPU. An example of this approach is the Masked Autoencoder for Distribution Estimation (MADE) [Germain et al., 2015], which drops out connections by multiplying the weight matrices of a fully-connected autoencoder with binary masks. Other mechanisms for dropping out connections include masked convolutions [van den Oord et al., 2016c] and causal convolutions [van den Oord et al., 2016a].

2.2 Normalizing flows

A normalizing flow [Rezende and Mohamed, 2015] represents as an invertible differentiable transformation of a base density . That is, where . The base density is chosen such that it can be easily evaluated for any input (a common choice for is a standard Gaussian). Under the invertibility assumption for , the density can be calculated as


In order for Equation (1) to be tractable, the transformation must be constructed such that (a) it is easy to invert, and (b) the determinant of its Jacobian is easy to compute. An important point is that if transformations and have the above properties, then their composition also has these properties. In other words, the transformation can be made deeper by composing multiple instances of it, and the result will still be a valid normalizing flow.

There have been various approaches in developing normalizing flows. An early example is Gaussianization [Chen and Gopinath, 2001]

, which is based on successive application of independent component analysis. Enforcing invertibility with nonsingular weight matrices has been proposed

[Rippel and Adams, 2013, Ballé et al., 2016], however in such approaches calculating the determinant of the Jacobian scales cubicly with data dimensionality in general. Planar/radial flows [Rezende and Mohamed, 2015] and Inverse Autoregressive Flow (IAF) [Kingma et al., 2016] are models whose Jacobian is tractable by design. However, they were developed primarily for variational inference and are not well-suited for density estimation, as they can only efficiently calculate the density of their own samples and not of externally provided datapoints. The Non-linear Independent Components Estimator (NICE) [Dinh et al., 2014] and its successor Real NVP [Dinh et al., 2017] have a tractable Jacobian and are also suitable for density estimation. IAF, NICE and Real NVP are discussed in more detail in Section 3.

3 Masked Autoregressive Flow

(a) Target density
(b) MADE with Gaussian conditionals
(c) MAF with layers
Figure 1: (a) The density to be learnt, defined as . (b) The density learnt by a MADE with order and Gaussian conditionals. Scatter plot shows the train data transformed into random numbers

; the non-Gaussian distribution indicates that the model is a poor fit.

(c) Learnt density and transformed train data of a layer MAF with the same order .

3.1 Autoregressive models as normalizing flows

Consider an autoregressive model whose conditionals are parameterized as single Gaussians. That is, the conditional is given by


In the above, and

are unconstrained scalar functions that compute the mean and log standard deviation of the

conditional given all previous variables. We can generate data from the above model using the following recursion:


In the above, is the vector of random numbers the model uses internally to generate data, typically by making calls to a random number generator often called randn().

Equation (3) provides an alternative characterization of the autoregressive model as a transformation from the space of random numbers to the space of data . That is, we can express the model as where . By construction, is easily invertible. Given a datapoint , the random numbers that were used to generate it are obtained by the following recursion:


Due to the autoregressive structure, the Jacobian of is triangular by design, hence its absolute determinant can be easily obtained as follows:


It follows that the autoregressive model can be equivalently interpreted as a normalizing flow, whose density can be obtained by substituting Equations (4) and (5) into Equation (1). This observation was first pointed out by Kingma et al. [2016].

A useful diagnostic for assessing whether an autoregressive model of the above type fits the target density well is to transform the train data into corresponding random numbers using Equation (4), and assess whether the ’s come from independent standard normals. If the ’s do not seem to come from independent standard normals, this is evidence that the model is a bad fit. For instance, Figure (b)b shows that the scatter plot of the random numbers associated with the train data can look significantly non-Gaussian if the model fits the target density poorly.

Here we interpret autoregressive models as a flow, and improve the model fit by stacking multiple instances of the model into a deeper flow. Given autoregressive models , we model the density of the random numbers of with , model the random numbers of with and so on, finally modelling the random numbers of with a standard Gaussian. This stacking adds flexibility: for example, Figure (c)c demonstrates that a flow of autoregressive models is able to learn multimodal conditionals, even though each model has unimodal conditionals. Stacking has previously been used in a similar way to improve model fit of deep belief nets Hinton et al. [2006] and deep mixtures of factor analyzers Tang et al. [2012].

We choose to implement the set of functions with masking, following the approach used by MADE [Germain et al., 2015]. MADE is a feedforward network that takes as input and outputs and for all with a single forward pass. The autoregressive property is enforced by multiplying the weight matrices of MADE with suitably constructed binary masks. In other words, we use MADE with Gaussian conditionals as the building layer of our flow. The benefit of using masking is that it enables transforming from data to random numbers and thus calculating in one forward pass through the flow, thus eliminating the need for sequential recursion as in Equation (4). We call this implementation of stacking MADEs into a flow Masked Autoregressive Flow (MAF).

3.2 Relationship with Inverse Autoregressive Flow

Like MAF, Inverse Autoregressive Flow (IAF) [Kingma et al., 2016] is a normalizing flow which uses MADE as its component layer. Each layer of IAF is defined by the following recursion:


Similarly to MAF, functions are computed using a MADE with Gaussian conditionals. The difference is architectural: in MAF and are directly computed from previous data variables , whereas in IAF and are directly computed from previous random numbers .

The consequence of the above is that MAF and IAF are different models with different computational trade-offs. MAF is capable of calculating the density of any datapoint in one pass through the model, however sampling from it requires performing sequential passes (where is the dimensionality of ). In contrast, IAF can generate samples and calculate their density with one pass, however calculating the density of an externally provided datapoint requires passes to find the random numbers associated with . Hence, the design choice of whether to connect and directly to (obtaining MAF) or to (obtaining IAF) depends on the intended usage. IAF is suitable as a recognition model for stochastic variational inference [Kingma and Welling, 2014, Rezende et al., 2014], where it only ever needs to calculate the density of its own samples. In contrast, MAF is more suitable for density estimation, because each example requires only one pass through the model whereas IAF requires .

A theoretical equivalence between MAF and IAF is that training a MAF with maximum likelihood corresponds to fitting an implicit IAF to the base density with stochastic variational inference. Let be the data density we wish to learn, be the base density, and be the transformation from to as implemented by MAF. The density defined by MAF (with added subscript for disambiguation) is


The inverse transformation from to can be seen as describing an implicit IAF with base density , which defines the following implicit density over the space:


Training MAF by maximizing the total log likelihood on train data corresponds to fitting to by stochastically minimizing . In Section A of the appendix, we show that


Hence, stochastically minimizing is equivalent to fitting to by minimizing

. Since the latter is the loss function used in variational inference, and

can be seen as an IAF with base density and transformation , it follows that training MAF as a density estimator of is equivalent to performing stochastic variational inference with an implicit IAF, where the posterior is taken to be the base density and the transformation implements the reparameterization trick [Kingma and Welling, 2014, Rezende et al., 2014]. This argument is presented in more detail in Section A of the appendix.

3.3 Relationship with Real NVP

Real NVP [Dinh et al., 2017] (NVP stands for Non Volume Preserving) is a normalizing flow obtained by stacking coupling layers. A coupling layer is an invertible transformation from random numbers to data with a tractable Jacobian, defined by


In the above, denotes elementwise multiplication, and the is applied to each element of . The transformation copies the first elements, and scales and shifts the remaining elements, with the amount of scaling and shifting being a function of the first elements. When stacking coupling layers into a flow, the elements are permuted across layers so that a different set of elements is copied each time. A special case of the coupling layer where is used by NICE [Dinh et al., 2014].

We can see that the coupling layer is a special case of both the autoregressive transformation used by MAF in Equation (3), and the autoregressive transformation used by IAF in Equation (6). Indeed, we can recover the coupling layer from the autoregressive transformation of MAF by setting for and making and functions of only for (for IAF we need to make and functions of instead for ). In other words, both MAF and IAF can be seen as more flexible (but different) generalizations of Real NVP, where each element is individually scaled and shifted as a function of all previous elements. The advantage of Real NVP compared to MAF and IAF is that it can both generate data and estimate densities with one forward pass only, whereas MAF would need passes to generate data and IAF would need passes to estimate densities.

3.4 Conditional MAF

Given a set of example pairs , conditional density estimation is the task of estimating the conditional density . Autoregressive modelling extends naturally to conditional density estimation. Each term in the chain rule of probability can be conditioned on side-information , decomposing any conditional density as . Therefore, we can turn any unconditional autoregressive model into a conditional one by augmenting its set of input variables with and only modelling the conditionals that correspond to . Any order of the variables can be chosen, as long as comes before . In masked autoregressive models, no connections need to be dropped from the inputs to the rest of the network.

We can implement a conditional version of MAF by stacking MADEs that were made conditional using the above strategy. That is, in a conditional MAF, the vector becomes an additional input for every layer. As a special case of MAF, Real NVP can be made conditional in the same way.

4 Experiments

4.1 Implementation and setup

We systematically evaluate three types of density estimator (MADE, Real NVP and MAF) in terms of density estimation performance on a variety of datasets. Code for reproducing our experiments (which uses Theano

[Theano Development Team, 2016]) can be found at

MADE. We consider two versions: (a) a MADE with Gaussian conditionals, denoted simply by MADE, and (b) a MADE whose conditionals are each parameterized as a mixture of Gaussians, denoted by MADE MoG. We used in all our experiments. MADE can be seen either as a MADE MoG with , or as a MAF with only one autoregressive layer. Adding more Gaussian components per conditional or stacking MADEs to form a MAF are two alternative ways of increasing the flexibility of MADE, which we are interested in comparing.

Real NVP. We consider a general-purpose implementation of the coupling layer, which uses two feedforward neural networks, implementing the scaling function and the shifting function respectively. Both networks have the same architecture, except that has hyperbolic tangent hidden units, whereas has rectified linear hidden units (we found this combination to perform best). Both networks have a linear output. We consider Real NVPs with either or coupling layers, denoted by Real NVP () and Real NVP (

) respectively, and in both cases the base density is a standard Gaussian. Successive coupling layers alternate between (a) copying the odd-indexed variables and transforming the even-indexed variables, and (b) copying the even-indexed variables and transforming the odd-indexed variables.

It is important to clarify that this is a general-purpose implementation of Real NVP which is different and thus not comparable to its original version [Dinh et al., 2017], which was designed specifically for image data. Here we are interested in comparing coupling layers with autoregressive layers as building blocks of normalizing flows for general-purpose density estimation tasks, and our design of Real NVP is such that a fair comparison between the two can be made.

MAF. We consider three versions: (a) a MAF with autoregressive layers and a standard Gaussian as a base density , denoted by MAF (), (b) a MAF with autoregressive layers and a standard Gaussian as a base density, denoted by MAF (), and (c) a MAF with autoregressive layers and a MADE MoG with Gaussian components in each conditional as a base density, denoted by MAF MoG (). MAF MoG () can be thought of as a MAF () stacked on top of a MADE MoG and trained jointly with it.

In all experiments, MADE and MADE MoG order the inputs using the order that comes with the dataset by default; no alternative orders were considered. MAF uses the default order for the first autoregressive layer (i.e. the layer that directly models the data) and reverses the order for each successive layer (the same was done for IAF by Kingma et al. [2016]).

MADE, MADE MoG and each layer in MAF is a feedforward neural network with masked weight matrices, such that the autoregressive property holds. The procedure for designing the masks (due to Germain et al. [2015]) is as follows. Each input or hidden unit is assigned a degree, which is an integer ranging from to , where is the data dimensionality. The degree of an input is taken to be its index in the order. The outputs have degrees that sequentially range from to . A unit is allowed to receive input only from units with lower or equal degree, which enforces the autoregressive property. In order for output to be connected to all inputs with degree less than , and thus make sure that no conditional independences are introduced, it is both necessary and sufficient that every hidden layer contains every degree. In all experiments except for CIFAR-10, we sequentially assign degrees within each hidden layer and use enough hidden units to make sure that all degrees appear. Because CIFAR-10 is high-dimensional, we used fewer hidden units than inputs and assigned degrees to hidden units uniformly at random (as was done by Germain et al. [2015]).

We added batch normalization

[Ioffe and Szegedy, 2015] after each coupling layer in Real NVP and after each autoregressive layer in MAF. Batch normalization is an elementwise scaling and shifting operation, which is easily invertible and has a tractable Jacobian, and thus it is suitable for use in a normalizing flow. We found that batch normalization in Real NVP and MAF reduces training time, increases stability during training and improves performance (a fact that was also observed by Dinh et al. [2017] for Real NVP). Section B of the appendix discusses our implementation of batch normalization and its use in normalizing flows.

All models were trained with the Adam optimizer [Kingma and Ba, 2015], using a minibatch size of 100, and a step size of for MADE and MADE MoG, and of for Real NVP and MAF. A small amount of regularization was added, with coefficient . Each model was trained with early stopping until no improvement occurred for

consecutive epochs on the validation set. For each model, we selected the number of hidden layers and number of hidden units based on validation performance (we gave the same options to all models), as described in Section 

D of the appendix.

4.2 Unconditional density estimation

Real NVP ()
Real NVP ()
MAF ()
MAF ()
MAF MoG ()
Table 1: Average test log likelihood (in nats) for unconditional density estimation. The best performing model for each dataset is shown in bold (multiple models are highlighted if the difference is not statistically significant according to a paired -test). Error bars correspond to standard deviations.

Following Uria et al. [2013], we perform unconditional density estimation on four UCI datasets (POWER, GAS, HEPMASS, MINIBOONE) and on a dataset of natural image patches (BSDS300).

UCI datasets. These datasets were taken from the UCI machine learning repository [Lichman, 2013]. We selected different datasets than Uria et al. [2013]

, because the ones they used were much smaller, resulting in an expensive cross-validation procedure involving a separate hyperparameter search for each fold. However, our data preprocessing follows

Uria et al. [2013]. The sample mean was subtracted from the data and each feature was divided by its sample standard deviation. Discrete-valued attributes were eliminated, as well as every attribute with a Pearson correlation coefficient greater than . These procedures are meant to avoid trivial high densities, which would make the comparison between approaches hard to interpret. Section D of the appendix gives more details about the UCI datasets and the individual preprocessing done on each of them.

Image patches. This dataset was obtained by extracting random monochrome patches from the BSDS300 dataset of natural images [Martin et al., 2001]. We used the same preprocessing as by Uria et al. [2013]. Uniform noise was added to dequantize pixel values, which was then rescaled to be in the range . The mean pixel value was subtracted from each patch, and the bottom-right pixel was discarded.

Table 1 shows the performance of each model on each dataset. A Gaussian fitted to the train data is reported as a baseline. We can see that on out of datasets MAF is the best performing model, with MADE MoG being the best performing model on the other . On all datasets, MAF outperforms Real NVP. For the MINIBOONE dataset, due to overlapping error bars, a pairwise comparison was done to determine which model performs the best, the results of which are reported in Section E of the appendix. MAF MoG () achieves the best reported result on BSDS300 for a single model with nats, followed by Deep RNADE [Uria et al., 2014] with . An ensemble of Deep RNADEs was reported to achieve nats [Uria et al., 2014]. The UCI datasets were used for the first time in the literature for density estimation, so no comparison with existing work can be made yet.

4.3 Conditional density estimation

For conditional density estimation, we used the MNIST dataset of handwritten digits [LeCun et al., ] and the CIFAR-10 dataset of natural images [Krizhevsky and Hinton, 2009]. In both datasets, each datapoint comes from one of distinct classes. We represent the class label as a

-dimensional, one-hot encoded vector

, and we model the density , where represents an image. At test time, we evaluate the probability of a test image by , where is a uniform prior over the labels. For comparison, we also train every model as an unconditional density estimator and report both results.

For both MNIST and CIFAR-10, we use the same preprocessing as by Dinh et al. [2017]. We dequantize pixel values by adding uniform noise, and then rescale them to

. We transform the rescaled pixel values into logit space by

, where for MNIST and for CIFAR-10, and perform density estimation in that space. In the case of CIFAR-10, we also augment the train set with horizontal flips of all train examples (as also done by Dinh et al. [2017]).

Table 2 shows the results on MNIST and CIFAR-10. The performance of a class-conditional Gaussian is reported as a baseline for the conditional case. Log likelihoods are calculated in logit space. MADE MoG is the best performing model on MNIST, whereas MAF is the best performing model on CIFAR-10. On CIFAR-10, both MADE and MADE MoG performed significantly worse than the Gaussian baseline. MAF outperforms Real NVP in all cases. To facilitate comparison with the literature, Section E of the appendix reports results in bits/pixel.111In earlier versions of the paper, results marked with * were reported incorrectly, due to an error in the calculation of the test log likelihood of conditional MAF. Thus error has been corrected in the current version.

unconditional conditional unconditional conditional
Real NVP ()
Real NVP ()
MAF () * *
MAF () * *
MAF MoG ()
Table 2: Average test log likelihood (in nats) for conditional density estimation. The best performing model for each dataset is shown in bold. Error bars correspond to standard deviations.

5 Discussion

We showed that we can improve MADE by modelling the density of its internal random numbers. Alternatively, MADE can be improved by increasing the flexibility of its conditionals. The comparison between MAF and MADE MoG showed that the best approach is dataset specific; in our experiments MAF outperformed MADE MoG in out of cases, which is strong evidence of its competitiveness. MADE MoG is a universal density approximator; with sufficiently many hidden units and Gaussian components, it can approximate any continuous density arbitrarily well. It is an open question whether MAF with a Gaussian base density has a similar property (MAF MoG clearly does).

We also showed that the coupling layer used in Real NVP is a special case of the autoregressive layer used in MAF. In fact, MAF outperformed Real NVP in all our experiments. Real NVP has achieved impressive performance in image modelling by incorporating knowledge about image structure. Our results suggest that replacing coupling layers with autoregressive layers in the original version of Real NVP is a promising direction for further improving its performance. Real NVP maintains however the advantage over MAF (and autoregressive models in general) that samples from the model can be generated efficiently in parallel.

Density estimation is one of several types of generative modelling, with the focus on obtaining accurate densities. However, we know that accurate densities do not necessarily imply good performance in other tasks, such as in data generation [Theis et al., 2016]. Alternative approaches to generative modelling include variational autoencoders [Kingma and Welling, 2014, Rezende et al., 2014], which are capable of efficient inference of their (potentially interpretable) latent space, and generative adversarial networks [Goodfellow et al., 2014], which are capable of high quality data generation. Choice of method should be informed by whether the application at hand calls for accurate densities, latent space inference or high quality samples. Masked Autoregressive Flow is a contribution towards the first of these goals.


We thank Maria Gorinova for useful comments, and Johann Brehmer for discovering the error in the calculation of the test log likelihood for conditional MAF. George Papamakarios and Theo Pavlakou were supported by the Centre for Doctoral Training in Data Science, funded by EPSRC (grant EP/L016427/1) and the University of Edinburgh. George Papamakarios was also supported by Microsoft Research through its PhD Scholarship Programme.

Appendix A Equivalence between MAF and IAF

In this section, we present the equivalence between MAF and IAF in full mathematical detail. Let be the true density the train data is sampled from. Suppose we have a MAF whose base density is , and whose transformation from to is . The MAF defines the following density over the space:


Using the definition of in Equation (11

), we can write the Kullback–Leibler divergence from

to as follows:


The inverse transformation from to can be seen as describing an implicit IAF with base density , which would define the following density over the space:


By making the change of variables in Equation (13) and using the definition of in Equation (14) we obtain


Equation (16) is the definition of the KL divergence from to , hence


Suppose now that we wish to fit the implicit density to the base density by minimizing the above KL. This corresponds exactly to the objective minimized when employing IAF as a recognition network in stochastic variational inference [Kingma et al., 2016], where would be the (typically intractable) posterior. The first step in stochastic variational inference would be to rewrite the expectation in Equation (16) with respect to the base distribution used by IAF, which corresponds exactly to Equation (13). This is often referred to as the reparameterization trick [Kingma and Welling, 2014, Rezende et al., 2014]. The second step would be to approximate Equation (13) with Monte Carlo, using samples drawn from , as follows:


Using the definition of in Equation (11), we can rewrite Equation (19) as


Since samples drawn from correspond precisely to the train data for MAF, we can recognize in Equation (20) the training objective for MAF. In conclusion, training a MAF by maximizing its total log likelihood on train data is equivalent to variationally training an implicit IAF with MAF’s base distribution as its target.

Appendix B Batch normalization

In our implementation of MAF, we inserted a batch normalization layer [Ioffe and Szegedy, 2015] between every two autoregressive layers, and between the last autoregressive layer and the base distribution. We did the same for Real NVP (the original implementation of Real NVP also uses batch normalization layers between coupling layers [Dinh et al., 2017]). The purpose of a batch normalization layer is to normalize its inputs

to have approximately zero mean and unit variance. In this section, we describe in full detail our implementation of batch normalization and its use as a layer in normalizing flows.

A batch normalization layer can be thought of as a transformation between two vectors of the same dimensionality. For consistency with our notation for autoregressive and coupling layers, let be the vector closer to the data, and be the vector closer to the base distribution. Batch normalization implements the transformation defined by


In the above, denotes elementwise multiplication. All other operations are to be understood elementwise. The inverse transformation is given by


and the absolute determinant of its Jacobian is


Vectors and are parameters of the transformation that are learnt during training. In typical implementations of batch normalization, parameter is not exponentiated. In our implementation, we chose to exponentiate in order to ensure its positivity and simplify the expression of the log absolute determinant. Parameters and correspond to the mean and variance of respectively. During training, we set and equal to the sample mean and variance of the current minibatch (we used minibatches of 100 examples). At validation and test time, we set them equal to the sample mean and variance of the entire train set. Other implementations use averages over minibatches [Ioffe and Szegedy, 2015] or maintain running averages during training [Dinh et al., 2017]. Finally, is a hyperparameter that ensures numerical stability if any of the elements of is near zero. In our experiments, we used .

Appendix C Number of parameters

To get a better idea of the computational trade-offs between different model choices versus the performance gains they achieve, we compare the number of parameters for each model. We only count connection weights, as they contribute the most, and ignore biases and batch normalization parameters. We assume that masking reduces the number of connections by approximately half.

For all models, let be the number of inputs, be the number of units in a hidden layer and be the number of hidden layers. We assume that all hidden layers have the same number of units (as we did in our experiments). For MAF MoG, let be the number of components per conditional. For Real NVP and MAF, let be the number of coupling layers/autoregressive layers respectively. Table 3 lists the number of parameters for each model.

For each extra component we add to MADE MoG, we increase the number of parameters by . For each extra autoregressive layer we add to MAF, we increase the number of parameters by . If we have one or two hidden layers (as we did in our experiments) and assume that is comparable to , the number of extra parameters in both cases is about the same. In other words, increasing flexibility by stacking has a parameter cost that is similar to adding more components to the conditionals, as long as the number of hidden layers is small.

Comparing Real NVP with MAF, we can see that Real NVP has about to times more parameters than a MAF of comparable size. Given that our experiments show that Real NVP is less flexible than a MAF of comparable size, we can conclude that MAF makes better use of its available capacity. The number of parameters of Real NVP could be reduced by tying weights between the scaling and shifting networks.

# of parameters
Real NVP
Table 3: Approximate number of parameters for each model, as measured by number of connection weights. Biases and batch normalization parameters are ignored.

Appendix D Additional experimental details

d.1 Models

MADE, MADE MoG and each autoregressive layer in MAF is a feedforward neural network (with masked weight matrices), with hidden layers of hidden units each. Similarly, each coupling layer in Real NVP contains two feedforward neural networks, one for scaling and one for shifting, each of which also has hidden layers of hidden units each. For each dataset, we gave a number of options for and (the same options where given to all models) and for each model we selected the option that performed best on the validation set. Table 4 lists the combinations of and that were given as options for each dataset.

In terms of nonlinearity for the hidden units, MADE, MADE MoG and MAF used rectified linear units, except for the GAS datasets where we used hyperbolic tangent units. In the coupling layer of Real NVP, we used hyberbolic tangent hidden units for the scaling network and rectified linear hidden units for the shifting network.

Table 4: Number of hidden layers and number of hidden units given as options for each dataset. Each combination is reported in the format .

d.2 Datasets

In the following paragraphs, we give a brief description of the four UCI datasets (POWER, GAS, HEPMASS, MINIBOONE) and of the way they were preprocessed.

POWER. The POWER dataset [POW, ] contains measurements of electric power consumption in a household over a period of months. It is actually a time series but was treated as if each example were an i.i.d. sample from the marginal distribution. The time feature was turned into an integer for the number of minutes in the day and then uniform random noise was added to it. The date was discarded, along with the global reactive power parameter, which seemed to have many values at exactly zero, which could have caused arbitrarily large spikes in the learnt distribution. Uniform random noise was added to each feature in the interval , where is large enough to ensure that with high probability there are no identical values for the feature but small enough to not change the data values significantly.

GAS. Created by Fonollosa et al. [2015], this dataset represents the readings of an array of chemical sensors exposed to gas mixtures over a hour period. Similarly to POWER, it is a time series but was treated as if each example were an i.i.d. sample from the marginal distribution. Only the data from the file ethylene_CO.txt was used, which corresponds to a mixture of ethylene and carbon monoxide. After removing strongly correlated attributes, the dimensionality was reduced to .

HEPMASS. Used by Baldi et al. [2016], this dataset describes particle collisions in high energy physics. Half of the data are examples of particle-producing collisions (positive), whereas the rest come from a background source (negative). Here we used the positive examples from the “” dataset, where the particle mass is . Five features were removed because they had too many reoccurring values; values that repeat too often can result in spikes in the density and misleading results.

MINIBOONE. Used by Roe et al. [2005]

, this dataset comes from the MiniBooNE experiment at Fermilab. Similarly to HEPMASS, it contains a number of positive examples (electron neutrinos) and a number of negative examples (muon neutrinos). Here we use the positive examples. These had some obvious outliers (

) which had values at exactly for every column and were removed. Also, seven of the features had far too high a count for a particular value, e.g. , so these were removed as well.

Table 5 lists the dimensionality and the number of train, validation and test examples for all seven datasets. The first three datasets in Table 5 were subsampled so that the product of the dimensionality and number of examples would be approximately . For the four UCI datasets, of the data was held out and used as test data and of the remaining data was used as validation data. From the BSDS300 dataset we randomly extracted patches for training, patches for validation and patches for testing. For MNIST and CIFAR-10 we held out of the train data for validation. We augmented the CIFAR-10 train set with the horizontal flips of all remaining train examples.

train validation test
Table 5: Dimensionality and number of examples for each dataset.

Appendix E Additional results

e.1 Pairwise comparison

On the MINIBOONE dataset, the model with highest average test log likelihood is MAF MoG (). However, due to the relatively small size of this dataset, the average test log likelihoods of some other models have overlapping error bars with that of MAF MoG (). To assess whether the differences are statistically significant, we performed a pairwise comparison, which is a more powerful statistical test. In particular, we calculated the difference in test log probability between every other model and MAF MoG () on each test example, and assessed whether this difference is significantly positive, which would indicate that MAF MoG () performs significantly better. The results of this comparison are shown in Table 6. We can see that MAF MoG () is significantly better than all other models except for MAF ().

Real NVP ()
Real NVP ()
MAF ()
MAF ()
MAF MoG ()
Table 6: Pairwise comparison results for MINIBOONE. Values correspond to average difference in log probability (in nats) from the best performing model, i.e. MAF MoG (). Error bars correspond to standard deviations. Significantly positive values indicate that MAF MoG () performs better.

e.2 Bits per pixel

In the main text, the results for MNIST and CIFAR-10 were reported in log likelihoods in logit space, since this is the objective that the models were trained to optimize. For comparison with other results in the literature, in Table 7 we report the same results in bits per pixel. For CIFAR-10, different colour components count as different pixels (i.e. an image is thought of as having pixels).

In order to calculate bits per pixel, we need to transform the densities returned by a model (which refer to logit space) back to image space in the range . Let be an image of pixels in logit space and be the corresponding image in image space. The transformation from to is


where for MNIST and for CIFAR-10. If is the density in logit space as returned by the model, using the above transformation the density of can be calculated as



is the logistic sigmoid function. From that, we can calculate the bits per pixel

of image as follows:


The above equation was used to convert between the average log likelihoods reported in the main text and the results of Table 7.

unconditional conditional unconditional conditional
Real NVP ()
Real NVP ()
MAF () * *
MAF () * *
MAF MoG ()
Table 7: Bits per pixel for conditional density estimation (lower is better). The best performing model for each dataset is shown in bold. Error bars correspond to standard deviations.

e.3 Generated images

Figures 2, 3 and 4 show generated images and real examples for BSDS300, MNIST and CIFAR-10 respectively. Images were generated by MAF MoG () for BSDS300, conditional MAF () for MNIST, and conditional MAF () for CIFAR-10.

The BSDS300 generated images are visually indistinguishable from the real ones. For MNIST and CIFAR-10, generated images lack the fidelity produced by modern image-based generative approaches, such as RealNVP [Dinh et al., 2017] or PixelCNN++ [Salimans et al., 2017]. This is because our version of MAF has no knowledge about image structure, as it was designed for general-purpose density estimation and not for realistic-looking image synthesis. However, if the latter is desired, it would be possible to incorporate image modelling techniques in the design of MAF (such as convolutions or a multi-scale architecture as used by Real NVP [Dinh et al., 2017]) in order to improve quality of generated images.

(a) Generated images
(b) Real images
Figure 2: Generated and real images from BSDS300.
(a) Generated images
(b) Real images
Figure 3: Class-conditional generated and real images from MNIST. Rows are different classes. Generated images are sorted by decreasing log likelihood from left to right.
(a) Generated images
(b) Real images
Figure 4: Class-conditional generated and real images from CIFAR-10. Rows are different classes. Generated images are sorted by decreasing log likelihood from left to right.


  • [1] Individual household electric power consumption data set. Accessed on 15 May 2017.
  • Baldi et al. [2016] P. Baldi, K. Cranmer, T. Faucett, P. Sadowski, and D. Whiteson. Parameterized machine learning for high-energy physics. arXiv:1601.07913, 2016.
  • Ballé et al. [2016] J. Ballé, V. Laparra, and E. P. Simoncelli. Density modeling of images using a generalized normalization transformation. Proceedings of the 4nd International Conference on Learning Representations, 2016.
  • Chen and Gopinath [2001] S. S. Chen and R. A. Gopinath. Gaussianization. Advances in Neural Information Processing Systems 13, pages 423–429, 2001.
  • Dinh et al. [2014] L. Dinh, D. Krueger, and Y. Bengio. NICE: Non-linear Independent Components Estimation. arXiv:1410.8516, 2014.
  • Dinh et al. [2017] L. Dinh, J. Sohl-Dickstein, and S. Bengio. Density estimation using Real NVP. Proceedings of the 5th International Conference on Learning Representations, 2017.
  • Fan et al. [2013] Y. Fan, D. J. Nott, and S. A. Sisson. Approximate Bayesian computation via regression density estimation. Stat, 2(1):34–48, 2013.
  • Fonollosa et al. [2015] J. Fonollosa, S. Sheik, R. Huerta, and S. Marco. Reservoir computing compensates slow response of chemosensor arrays exposed to fast varying gas concentrations in continuous monitoring. Sensors and Actuators B: Chemical, 215:618–629, 2015.
  • Germain et al. [2015] M. Germain, K. Gregor, I. Murray, and H. Larochelle. MADE: Masked Autoencoder for Distribution Estimation. Proceedings of the 32nd International Conference on Machine Learning, pages 881–889, 2015.
  • Goodfellow et al. [2014] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio. Generative adversarial nets. Advances in Neural Information Processing Systems 27, pages 2672–2680, 2014.
  • Gu et al. [2015] S. Gu, Z. Ghahramani, and R. E. Turner. Neural adaptive sequential Monte Carlo. Advances in Neural Information Processing Systems 28, pages 2629–2637, 2015.
  • Hinton et al. [2006] G. Hinton, S. Osindero, and Y.-W. Teh. A fast learning algorithm for deep belief nets. Neural Computation, 18(7):1527–1554, 2006.
  • Ioffe and Szegedy [2015] S. Ioffe and C. Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. Proceedings of the 32nd International Conference on Machine Learning, pages 448–456, 2015.
  • Kingma and Ba [2015] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. Proceedings of the 3rd International Conference on Learning Representations, 2015.
  • Kingma and Welling [2014] D. P. Kingma and M. Welling. Auto-encoding variational Bayes. Proceedings of the 2nd International Conference on Learning Representations, 2014.
  • Kingma et al. [2016] D. P. Kingma, T. Salimans, R. Jozefowicz, X. Chen, I. Sutskever, and M. Welling. Improved variational inference with Inverse Autoregressive Flow. Advances in Neural Information Processing Systems 29, pages 4743–4751, 2016.
  • Krizhevsky and Hinton [2009] A. Krizhevsky and G. Hinton. Learning multiple layers of features from tiny images. Technical report, University of Toronto, 2009.
  • Kulkarni et al. [2015] T. D. Kulkarni, P. Kohli, J. B. Tenenbaum, and V. Mansinghka. Picture: A probabilistic programming language for scene perception.

    IEEE Conference on Computer Vision and Pattern Recognition

    , pages 4390–4399, 2015.
  • Le et al. [2017] T. A. Le, A. G. Baydin, and F. Wood. Inference compilation and universal probabilistic programming.

    Proceedings of the 20th International Conference on Artificial Intelligence and Statistics

    , 2017.
  • [20] Y. LeCun, C. Cortes, and C. J. C. Burges.

    The MNIST database of handwritten digits.

  • Lichman [2013] M. Lichman. UCI machine learning repository, 2013. URL
  • Loaiza-Ganem et al. [2017] G. Loaiza-Ganem, Y. Gao, and J. P. Cunningham. Maximum entropy flow networks. Proceedings of the 5th International Conference on Learning Representations, 2017.
  • Martin et al. [2001] D. Martin, C. Fowlkes, D. Tal, and J. Malik. A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. pages 416–423, 2001.
  • Paige and Wood [2016] B. Paige and F. Wood. Inference networks for sequential Monte Carlo in graphical models. Proceedings of the 33rd International Conference on Machine Learning, 2016.
  • Papamakarios and Murray [2015] G. Papamakarios and I. Murray. Distilling intractable generative models, 2015. Probabilistic Integration Workshop at Neural Information Processing Systems 28.
  • Papamakarios and Murray [2016] G. Papamakarios and I. Murray. Fast -free inference of simulation models with Bayesian conditional density estimation. Advances in Neural Information Processing Systems 29, 2016.
  • Rezende and Mohamed [2015] D. J. Rezende and S. Mohamed. Variational inference with normalizing flows. Proceedings of the 32nd International Conference on Machine Learning, pages 1530–1538, 2015.
  • Rezende et al. [2014] D. J. Rezende, S. Mohamed, and D. Wierstra.

    Stochastic backpropagation and approximate inference in deep generative models.

    Proceedings of the 31st International Conference on Machine Learning, pages 1278–1286, 2014.
  • Rippel and Adams [2013] O. Rippel and R. P. Adams. High-dimensional probability estimation with deep density models. arXiv:1302.5125, 2013.
  • Roe et al. [2005] B. P. Roe, H.-J. Yang, J. Zhu, Y. Liu, I. Stancu, and G. McGregor.

    Boosted decision trees as an alternative to artificial neural networks for particle identification.

    Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 543(2–3):577–584, 2005.
  • Salimans et al. [2017] T. Salimans, A. Karpathy, X. Chen, and D. P. Kingma.