An Introduction to Variational Autoencoders

Diederik P. Kingma, Max Welling

Chapter 1 Introduction

One major division in machine learning is generative versus discriminative modeling. While in discriminative modeling one aims to learn a predictor given the observations, in generative modeling one aims to solve the more general problem of learning a joint distribution over all the variables. A generative model simulates how the data is generated in the real world. “Modeling” is understood in almost every science as unveiling this generating process by hypothesizing theories and testing these theories through observations. For instance, when meteorologists model the weather they use highly complex partial differential equations to express the underlying physics of the weather. Or when an astronomer models the formation of galaxies s/he encodes in his/her equations of motion the physical laws under which stellar bodies interact. The same is true for biologists, chemists, economists and so on. Modeling in the sciences is in fact almost always generative modeling.

There are many reasons why generative modeling is attractive. First, we can express physical laws and constraints into the generative process while details that we don’t know or care about, i.e. nuisance variables, are treated as noise. The resulting models are usually highly intuitive and interpretable and by testing them against observations we can confirm or reject our theories about how the world works.

Another reason for trying to understand the generative process of data is that it naturally expresses causal relations of the world. Causal relations have the great advantage that they generalize much better to new situations than mere correlations. For instance, once we understand the generative process of an earthquake, we can use that knowledge both in California and in Chile.

To turn a generative model into a discriminator, we need to use Bayes rule. For instance, we have a generative model for an earthquake of type A and another for type B, then seeing which of the two describes the data best we can compute a probability for whether earthquake A or B happened. Applying Bayes rule is however often computationally expensive.

In discriminative methods we directly learn a map in the same direction as we intend to make future predictions in. This is in the opposite direction than the generative model. For instance, one can argue that an image is generated in the world by first identifying the object, then generating the object in 3D and then projecting it onto an pixel grid. A discriminative model takes these pixel values directly as input and maps them to the labels. While generative models can learn efficiently from data, they also tend to make stronger assumptions on the data than their purely discriminative counterparts, often leading to higher asymptotic bias [banerjee2007analysis] when the model is wrong. For this reason, if the model is wrong (and it almost always is to some degree!), if one is solely interested in learning to discriminate, and one is in a regime with a sufficiently large amount of data, then purely discriminative models typically will lead to fewer errors in discriminative tasks. Nevertheless, depending on how much data is around, it may pay off to study the data generating process as a way to guide the training of the discriminator, such as a classifier. For instance, one may have few labeled examples and many more unlabeled examples. In this semi-supervised learning setting, one can use the generative model of the data to improve classification [kingma2014semi, sonderby2016train].

Generative modeling can be useful more generally. One can think of it as an auxiliary task. For instance, predicting the immediate future may help us build useful abstractions of the world that can be used for multiple prediction tasks downstream. This quest for disentangled, semantically meaningful, statistically independent and causal factors of variation in data is generally known as unsupervised representation learning, and the variational autoencoder (VAE) has been extensively employed for that purpose. Alternatively, one may view this as an implicit form of regularization: by forcing the representations to be meaningful for data generation, we bias the inverse of that process, which maps from input to representation, into a certain mould. The auxiliary task of predicting the world is used to better understand the world at an abstract level and thus to better make downstream predictions.

The VAE can be viewed as two coupled, but independently parameterized models: the encoder or recognition model, and the decoder or generative model. These two models support each other. The recognition model delivers to the generative model an approximation to its posterior over latent random variables, which it needs to update its parameters inside an iteration of “expectation maximization” learning. Reversely, the generative model is a scaffolding of sorts for the recognition model to learn meaningful representations of the data, including possibly class-labels. The recognition model is the approximate inverse of the generative model according to Bayes rule.

One advantage of the VAE framework, relative to ordinary Variational Inference (VI), is that the recognition model (also called inference model) is now a (stochastic) function of the input variables. This in contrast to VI where each data-case has a separate variational distribution, which is inefficient for large data-sets. The recognition model uses one set of parameters to model the relation between input and latent variables and as such is called “amortized inference”. This recognition model can be arbitrary complex but is still reasonably fast because by construction it can be done using a single feedforward pass from input to latent variables. However the price we pay is that this sampling induces sampling noise in the gradients required for learning. Perhaps the greatest contribution of the VAE framework is the realization that we can counteract this variance by using what is now known as the “reparameterization trick”, a simple procedure to reorganize our gradient computation that reduces variance in the gradients.

The VAE is inspired by the Helmholtz Machine [dayan1995helmholtz] which was perhaps the first model that employed a recognition model. However, its wake-sleep algorithm was inefficient and didn’t optimize a single objective. The VAE learning rules instead follow from a single approximation to the maximum likelihood objective.

VAEs marry graphical models and deep learning. The generative model is a Bayesian network of the form p(x∣z)p(z)p(\mathbf{x}|\mathbf{z})p(\mathbf{z}), or, if there are multiple stochastic latent layers, a hierarchy such as p(x∣zL)p(zL∣zL−1)p(\mathbf{x}|\mathbf{z}_{L})p(\mathbf{z}_{L}|\mathbf{z}_{L-1}) ...p(z1∣z0)...p(\mathbf{z}_{1}|\mathbf{z}_{0}). Similarly, the recognition model is also a conditional Bayesian network of the form q(z∣x)q(\mathbf{z}|\mathbf{x}) or as a hierarchy, such as q(z0∣z1)...q(zL∣X)q(\mathbf{z}_{0}|\mathbf{z}_{1})...q(\mathbf{z}_{L}|X). But inside each conditional may hide a complex (deep) neural network, e.g. z∣x∼f(x,ϵ)\mathbf{z}|\mathbf{x}\sim f(\mathbf{x},\boldsymbol{\epsilon}), with ff a neural network mapping and ϵ\boldsymbol{\epsilon} a noise random variable. Its learning algorithm is a mix of classical (amortized, variational) expectation maximization but through the reparameterization trick ends up backpropagating through the many layers of the deep neural networks embedded inside of it.

Since its inception, the VAE framework has been extended in many directions, e.g. to dynamical models [johnson2016composing], models with attention [gregor2015draw], models with multiple levels of stochastic latent variables [kingma2016improving], and many more. It has proven itself as a fertile framework to build new models in. More recently, another generative modeling paradigm has gained significant attention: the generative adversarial network (GAN) [goodfellow2014generative]. VAEs and GANs seem to have complementary properties: while GANs can generate images of high subjective perceptual quality, they tend to lack full support over the data [grover2018flow], as opposed to likelihood-based generative models. VAEs, like other likelihood-based models, generate more dispersed samples, but are better density models in terms of the likelihood criterion. As such many hybrid models have been proposed to try to represent the best of both worlds [dumoulin2016adversarially, grover2018flow, rosca2018distribution].

As a community we seem to have embraced the fact that generative models and unsupervised learning play an important role in building intelligent machines. We hope that the VAE provides a useful piece of that puzzle.

2 Aim

The framework of variational autoencoders (VAEs) [kingma2013auto, rezende2014stochastic] provides a principled method for jointly learning deep latent-variable models and corresponding inference models using stochastic gradient descent. The framework has a wide array of applications from generative modeling, semi-supervised learning to representation learning.

This work is meant as an expanded version of our earlier work [kingma2013auto], allowing us to explain the topic in finer detail and to discuss a selection of important follow-up work. This is not aimed to be a comprehensive review of all related work. We assume that the reader has basic knowledge of algebra, calculus and probability theory.

In this chapter we discuss background material: probabilistic models, directed graphical models, the marriage of directed graphical models with neural networks, learning in fully observed models and deep latent-variable models (DLVMs). In chapter 2 we explain the basics of VAEs. In chapter 3 we explain advanced inference techniques, followed by an explanation of advanced generative models in chapter 4. Please refer to section A.1 for more information on mathematical notation.

3 Probabilistic Models and Variational Inference

In the field of machine learning, we are often interested in learning probabilistic models of various natural and artificial phenomena from data. Probabilistic models are mathematical descriptions of such phenomena. They are useful for understanding such phenomena, for prediction of unknowns in the future, and for various forms of assisted or automated decision making. As such, probabilistic models formalize the notion of knowledge and skill, and are central constructs in the field of machine learning and AI.

As probabilistic models contain unknowns and the data rarely paints a complete picture of the unknowns, we typically need to assume some level of uncertainty over aspects of the model. The degree and nature of this uncertainty is specified in terms of (conditional) probability distributions. Models may consist of both continuous-valued variables and discrete-valued variables. The, in some sense, most complete forms of probabilistic models specify all correlations and higher-order dependencies between the variables in the model, in the form of a joint probability distribution over those variables.

Let’s use x\mathbf{x} as the vector representing the set of all observed variables whose joint distribution we would like to model. Note that for notational simplicity and to avoid clutter, we use lower case bold (e.g. x\mathbf{x}) to denote the underlying set of observed random variables, i.e. flattened and concatenated such that the set is represented as a single vector. See section A.1 for more on notation.

We assume the observed variable x\mathbf{x} is a random sample from an unknown underlying process, whose true (probability) distribution p∗(x)p^{*}(\mathbf{x}) is unknown. We attempt to approximate this underlying process with a chosen model pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}), with parameters θ\boldsymbol{\theta}:

Learning is, most commonly, the process of searching for a value of the parameters θ\boldsymbol{\theta} such that the probability distribution function given by the model, pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}), approximates the true distribution of the data, denoted by p∗(x)p^{*}(\mathbf{x}), such that for any observed x\mathbf{x}:

Naturally, we wish pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}) to be sufficiently flexible to be able to adapt to the data, such that we have a chance of obtaining a sufficiently accurate model. At the same time, we wish to be able to incorporate knowledge about the distribution of data into the model that is known a priori.

Often, such as in case of classification or regression problems, we are not interested in learning an unconditional model pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}), but a conditional model pθ(y∣x)p_{\boldsymbol{\theta}}(\mathbf{y}|\mathbf{x}) that approximates the underlying conditional distribution p∗(y∣x)p^{*}(\mathbf{y}|\mathbf{x}): a distribution over the values of variable y\mathbf{y}, conditioned on the value of an observed variable x\mathbf{x}. In this case, x\mathbf{x} is often called the input of the model. Like in the unconditional case, a model pθ(y∣x)p_{\boldsymbol{\theta}}(\mathbf{y}|\mathbf{x}) is chosen, and optimized to be close to the unknown underlying distribution, such that for any x\mathbf{x} and y\mathbf{y}:

A relatively common and simple example of conditional modeling is image classification, where x\mathbf{x} is an image, and y\mathbf{y} is the image’s class, as labeled by a human, which we wish to predict. In this case, pθ(y∣x)p_{\boldsymbol{\theta}}(\mathbf{y}|\mathbf{x}) is typically chosen to be a categorical distribution, whose parameters are computed from x\mathbf{x}.

Conditional models become more difficult to learn when the predicted variables are very high-dimensional, such as images, video or sound. One example is the reverse of the image classification problem: prediction of a distribution over images, conditioned on the class label. Another example with both high-dimensional input, and high-dimensional output, is time series prediction, such as text or video prediction.

To avoid notational clutter we will often assume unconditional modeling, but one should always keep in mind that the methods introduced in this work are, in almost all cases, equally applicable to conditional models. The data on which the model is conditioned, can be treated as inputs to the model, similar to the parameters of the model, with the obvious difference that one doesn’t optimize over their value.

4 Parameterizing Conditional Distributions with Neural Networks

Differentiable feed-forward neural networks, from here just called neural networks, are a particularly flexible and computationally scalable type of function approximator. Learning of models based on neural networks with multiple ’hidden’ layers of artificial neurons is often referred to as deep learning [goodfellow2016deeplearning, lecun2015deep]. A particularly interesting application is probabilistic models, i.e. the use of neural networks for probability density functions (PDFs) or probability mass functions (PMFs) in probabilistic models. Probabilistic models based on neural networks are computationally scalable since they allow for stochastic gradient-based optimization which, as we will explain, allows scaling to large models and large datasets. We will denote a deep neural network as a vector function: NeuralNet(.)\text{NeuralNet}(.).

At the time of writing, deep learning has been shown to work well for a large variety of classification and regression problems, as summarized in [lecun2015deep, goodfellow2016deeplearning]. In case of neural-network based image classification [lecun1998gradient], for example, neural networks parameterize a categorical distribution pθ(y∣x)p_{\boldsymbol{\theta}}(y|\mathbf{x}) over a class label yy, conditioned on an image x\mathbf{x}.

where the last operation of NeuralNet(.)\text{NeuralNet}(.) is typically a softmax() function such that ∑ipi=1\sum_{i}p_{i}=1.

5 Directed Graphical Models and Neural Networks

We work with directed probabilistic models, also called directed probabilistic graphical models (PGMs), or Bayesian networks. Directed graphical models are a type of probabilistic models where all the variables are topologically organized into a directed acyclic graph. The joint distribution over the variables of such models factorizes as a product of prior and conditional distributions:

where Pa(xj)Pa(\mathbf{x}_{j}) is the set of parent variables of node jj in the directed graph. For non-root-nodes, we condition on the parents. For root nodes, the set of parents is the empty set, such that the distribution is unconditional.

Traditionally, each conditional probability distribution pθ(xj∣Pa(xj))p_{\boldsymbol{\theta}}(\mathbf{x}_{j}|Pa(\mathbf{x}_{j})) is parameterized as a lookup table or a linear model [koller2009probabilistic]. As we explained above, a more flexible way to parameterize such conditional distributions is with neural networks. In this case, neural networks take as input the parents of a variable in a directed graph, and produce the distributional parameters η\boldsymbol{\eta} over that variable:

We will now discuss how to learn the parameters of such models, if all the variables are observed in the data.

6 Learning in Fully Observed Models with Neural Nets

If all variables in the directed graphical model are observed in the data, then we can compute and differentiate the log-probability of the data under the model, leading to relatively straightforward optimization.

We often collect a dataset D\mathcal{D} consisting of N≥1N\geq 1 datapoints:

The datapoints are assumed to be independent samples from an unchanging underlying distribution. In other words, the dataset is assumed to consist of distinct, independent measurements from the same (unchanging) system. In this case, the observations D={x(i)}i=1N\mathcal{D}=\{\mathbf{x}^{(i)}\}_{i=1}^{N} are said to be i.i.d., for independently and identically distributed. Under the i.i.d. assumption, the probability of the datapoints given the parameters factorizes as a product of individual datapoint probabilities. The log-probability assigned to the data by the model is therefore given by:

6.2 Maximum Likelihood and Minibatch SGD

The most common criterion for probabilistic models is maximum log-likelihood (ML). As we will explain, maximization of the log-likelihood criterion is equivalent to minimization of a Kullback Leibler divergence between the data and model distributions.

Under the ML criterion, we attempt to find the parameters θ\boldsymbol{\theta} that maximize the sum, or equivalently the average, of the log-probabilities assigned to the data by the model. With i.i.d. dataset D\mathcal{D} of size NDN_{\mathcal{D}}, the maximum likelihood objective is to maximize the log-probability given by equation (1.10).

Using calculus’ chain rule and automatic differentiation tools, we can efficiently compute gradients of this objective, i.e. the first derivatives of the objective w.r.t. its parameters θ\boldsymbol{\theta}. We can use such gradients to iteratively hill-climb to a local optimum of the ML objective. If we compute such gradients using all datapoints, ∇θlog⁡pθ(D)\nabla_{\boldsymbol{\theta}}\log p_{\boldsymbol{\theta}}(\mathcal{D}), then this is known as batch gradient descent. Computation of this derivative is, however, an expensive operation for large dataset size NDN_{\mathcal{D}}, since it scales linearly with NDN_{\mathcal{D}}.

A more efficient method for optimization is stochastic gradient descent (SGD) (section A.3), which uses randomly drawn minibatches of data M⊂D\mathcal{M}\subset\mathcal{D} of size NMN_{\mathcal{M}}. With such minibatches we can form an unbiased estimator of the ML criterion:

The ≃\simeq symbol means that one of the two sides is an unbiased estimator of the other side. So one side (in this case the right-hand side) is a random variable due to some noise source, and the two sides are equal when averaged over the noise distribution. The noise source, in this case, is the randomly drawn minibatch of data M\mathcal{M}. The unbiased estimator log⁡pθ(M)\log p_{\boldsymbol{\theta}}(\mathcal{M}) is differentiable, yielding the unbiased stochastic gradients:

These gradients can be plugged into stochastic gradient-based optimizers; see section A.3 for further discussion. In a nutshell, we can optimize the objective function by repeatedly taking small steps in the direction of the stochastic gradient.

6.3 Bayesian inference

From a Bayesian perspective, we can improve upon ML through maximum a posteriori (MAP) estimation (section section A.2.1), or, going even further, inference of a full approximate posterior distribution over the parameters (see section A.1.4).

7 Learning and Inference in Deep Latent Variable Models

We can extend fully-observed directed models, discussed in the previous section, into directed models with latent variables. Latent variables are variables that are part of the model, but which we don’t observe, and are therefore not part of the dataset. We typically use z\mathbf{z} to denote such latent variables. In case of unconditional modeling of observed variable x\mathbf{x}, the directed graphical model would then represent a joint distribution pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}) over both the observed variables x\mathbf{x} and the latent variables z\mathbf{z}. The marginal distribution over the observed variables pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}), is given by:

This is also called the (single datapoint) marginal likelihood or the model evidence, when taken as a function of θ\boldsymbol{\theta}.

Such an implicit distribution over x\mathbf{x} can be quite flexible. If z\mathbf{z} is discrete and pθ(x∣z)p_{\boldsymbol{\theta}}(\mathbf{x}|\mathbf{z}) is a Gaussian distribution, then pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}) is a mixture-of-Gaussians distribution. For continuous z\mathbf{z}, pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}) can be seen as an infinite mixture, which are potentially more powerful than discrete mixtures. Such marginal distributions are also called compound probability distributions.

7.2 Deep Latent Variable Models

We use the term deep latent variable model (DLVM) to denote a latent variable model pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}) whose distributions are parameterized by neural networks. Such a model can be conditioned on some context, like pθ(x,z∣y)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}|\mathbf{y}). One important advantage of DLVMs, is that even when each factor (prior or conditional distribution) in the directed model is relatively simple (such as conditional Gaussian), the marginal distribution pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}) can be very complex, i.e. contain almost arbitrary dependencies. This expressivity makes deep latent-variable models attractive for approximating complicated underlying distributions p∗(x)p^{*}(\mathbf{x}).

Perhaps the simplest, and most common, DLVM is one that is specified as factorization with the following structure:

where pθ(z)p_{\boldsymbol{\theta}}(\mathbf{z}) and/or pθ(x∣z)p_{\boldsymbol{\theta}}(\mathbf{x}|\mathbf{z}) are specified. The distribution p(z)p(\mathbf{z}) is often called the prior distribution over z\mathbf{z}, since it is not conditioned on any observations.

7.3 Example DLVM for multivariate Bernoulli data

A simple example DLVM, used in [kingma2013auto] for binary data x\mathbf{x}, is with a spherical Gaussian latent space, and a factorized Bernoulli observation model:

where ∀pj∈p:0≤pj≤1\forall p_{j}\in\mathbf{p}:0\leq p_{j}\leq 1 (e.g. implemented through a sigmoid nonlinearity as the last layer of the DecoderNeuralNetθ(.)\text{DecoderNeuralNet}_{\boldsymbol{\theta}}(.)), where DD is the dimensionality of x\mathbf{x}, and Bernoulli(.;p)\text{Bernoulli}(.;p) is the probability mass function (PMF) of the Bernoulli distribution.

8 Intractabilities

The main difficulty of maximum likelihood learning in DLVMs is that the marginal probability of data under the model is typically intractable. This is due to the integral in equation (1.13) for computing the marginal likelihood (or model evidence), pθ(x)=∫pθ(x,z) dzp_{\boldsymbol{\theta}}(\mathbf{x})=\int p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z})\,d\mathbf{z}, not having an analytic solution or efficient estimator. Due to this intractability, we cannot differentiate it w.r.t. its parameters and optimize it, as we can with fully observed models.

The intractability of pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}), is related to the intractability of the posterior distribution pθ(z∣x)p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x}). Note that the joint distribution pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}) is efficient to compute, and that the densities are related through the basic identity:

Since pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}) is tractable to compute, a tractable marginal likelihood pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}) leads to a tractable posterior pθ(z∣x)p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x}), and vice versa. Both are intractable in DLVMs.

Approximate inference techniques (see also section A.2) allow us to approximate the posterior pθ(z∣x)p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x}) and the marginal likelihood pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}) in DLVMs. Traditional inference methods are relatively expensive. Such methods, for example, often require a per-datapoint optimization loop, or yield bad posterior approximations. We would like to avoid such expensive procedures.

Likewise, the posterior over the parameters of (directed models parameterized with) neural networks, p(θ∣D)p(\boldsymbol{\theta}|\mathcal{D}), is generally intractable to compute exactly, and requires approximate inference techniques.

Chapter 2 Variational Autoencoders

In this chapter we explain the basics of variational autoencoders (VAEs).

In the previous chapter, we introduced deep latent-variable models (DLVMs), and the problem of estimating the log-likelihood and posterior distributions in such models. The framework of variational autoencoders (VAEs) provides a computationally efficient way for optimizing DLVMs jointly with a corresponding inference model using SGD.

To turn the DLVM’s intractable posterior inference and learning problems into tractable problems, we introduce a parametric inference model qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}). This model is also called an encoder or recognition model. With ϕ\boldsymbol{\phi} we indicate the parameters of this inference model, also called the variational parameters. We optimize the variational parameters ϕ\boldsymbol{\phi} such that:

As we will explain, this approximation to the posterior help us optimize the marginal likelihood.

Like a DLVM, the inference model can be (almost) any directed graphical model:

where Pa(zj)Pa(\mathbf{z}_{j}) is the set of parent variables of variable zj\mathbf{z}_{j} in the directed graph. And also similar to a DLVM, the distribution qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) can be parameterized using deep neural networks. In this case, the variational parameters ϕ\boldsymbol{\phi} include the weights and biases of the neural network. For example:

Typically, we use a single encoder neural network to perform posterior inference over all of the datapoints in our dataset. This can be contrasted to more traditional variational inference methods where the variational parameters are not shared, but instead separately and iteratively optimized per datapoint. The strategy used in VAEs of sharing variational parameters across datapoints is also called amortized variational inference [gershman2014amortized]. With amortized inference we can avoid a per-datapoint optimization loop, and leverage the efficiency of SGD.

2 Evidence Lower Bound (ELBO)

The optimization objective of the variational autoencoder, like in other variational methods, is the evidence lower bound, abbreviated as ELBO. An alternative term for this objective is variational lower bound. Typically, the ELBO is derived through Jensen’s inequality. Here we will use an alternative derivation that avoids Jensen’s inequality, providing greater insight about its tightness.

For any choice of inference model qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}), including the choice of variational parameters ϕ\boldsymbol{\phi}, we have:

The second term in eq. (2.8) is the Kullback-Leibler (KL) divergence between qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) and pθ(z∣x)p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x}), which is non-negative:

and zero if, and only if, qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) equals the true posterior distribution.

The first term in eq. (2.8) is the variational lower bound, also called the evidence lower bound (ELBO):

Due to the non-negativity of the KL divergence, the ELBO is a lower bound on the log-likelihood of the data.

So, interestingly, the KL divergence DKL(qϕ(z∣x)∣∣pθ(z∣x))D_{KL}(q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x})||p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x})) determines two ’distances’:

By definition, the KL divergence of the approximate posterior from the true posterior;

The gap between the ELBO Lθ,ϕ(x)\mathcal{L}_{\boldsymbol{\theta},\boldsymbol{\phi}}(\mathbf{x}) and the marginal likelihood log⁡pθ(x)\log p_{\boldsymbol{\theta}}(\mathbf{x}); this is also called the tightness of the bound. The better qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) approximates the true (posterior) distribution pθ(z∣x)p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x}), in terms of the KL divergence, the smaller the gap.

By looking at equation 2.11, it can be understood that maximization of the ELBO Lθ,ϕ(x)\mathcal{L}_{\boldsymbol{\theta},\boldsymbol{\phi}}(\mathbf{x}) w.r.t. the parameters θ\boldsymbol{\theta} and ϕ\boldsymbol{\phi}, will concurrently optimize the two things we care about:

It will approximately maximize the marginal likelihood pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}). This means that our generative model will become better.

It will minimize the KL divergence of the approximation qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) from the true posterior pθ(z∣x)p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x}), so qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) becomes better.

3 Stochastic Gradient-Based Optimization of the ELBO

An important property of the ELBO, is that it allows joint optimization w.r.t. all parameters (ϕ\boldsymbol{\phi} and θ\boldsymbol{\theta}) using stochastic gradient descent (SGD). We can start out with random initial values of ϕ\boldsymbol{\phi} and θ\boldsymbol{\theta}, and stochastically optimize their values until convergence.

Given a dataset with i.i.d. data, the ELBO objective is the sum (or average) of individual-datapoint ELBO’s:

Unbiased gradients of the ELBO w.r.t. the generative model parameters θ\boldsymbol{\theta} are simple to obtain:

The last line (eq. (2.17)) is a simple Monte Carlo estimator of the second line (eq. (2.15)), where z\mathbf{z} in the last two lines (eq. (2.16) and eq. (2.17)) is a random sample from qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}).

Unbiased gradients w.r.t. the variational parameters ϕ\boldsymbol{\phi} are more difficult to obtain, since the ELBO’s expectation is taken w.r.t. the distribution qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}), which is a function of ϕ\boldsymbol{\phi}. I.e., in general:

In the case of continuous latent variables, we can use a reparameterization trick for computing unbiased estimates of ∇θ,ϕLθ,ϕ(x)\nabla_{\boldsymbol{\theta},\boldsymbol{\phi}}\mathcal{L}_{\boldsymbol{\theta},\boldsymbol{\phi}}(\mathbf{x}), as we will now discuss. This stochastic estimate allows us to optimize the ELBO using SGD; see algorithm 1. See section 2.9.1 for a discussion of variational methods for discrete latent variables.

4 Reparameterization Trick

For continuous latent variables and a differentiable encoder and generative model, the ELBO can be straightforwardly differentiated w.r.t. both ϕ\boldsymbol{\phi} and θ\boldsymbol{\theta} through a change of variables, also called the reparameterization trick ([kingma2013auto] and [rezende2014stochastic]).

First, we express the random variable z∼qϕ(z∣x)\mathbf{z}\sim q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) as some differentiable (and invertible) transformation of another random variable ϵ\boldsymbol{\epsilon}, given z\mathbf{z} and ϕ\boldsymbol{\phi}:

where the distribution of random variable ϵ\boldsymbol{\epsilon} is independent of x\mathbf{x} or ϕ\boldsymbol{\phi}.

4.2 Gradient of expectation under change of variable

Given such a change of variable, expectations can be rewritten in terms of ϵ\boldsymbol{\epsilon},

where z=g(ϵ,ϕ,x)\mathbf{z}=\mathbf{g}(\boldsymbol{\epsilon},\boldsymbol{\phi},\mathbf{x}). and the expectation and gradient operators become commutative, and we can form a simple Monte Carlo estimator:

where in the last line, z=g(ϕ,x,ϵ)\mathbf{z}=\mathbf{g}(\boldsymbol{\phi},\mathbf{x},\boldsymbol{\epsilon}) with random noise sample ϵ∼p(ϵ)\boldsymbol{\epsilon}\sim p(\boldsymbol{\epsilon}). See figure 2.3 for an illustration and further clarification, and figure 3.2 for an illustration of the resulting posteriors for a 2D toy problem.

4.3 Gradient of ELBO

Under the reparameterization, we can replace an expectation w.r.t. qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) with one w.r.t. p(ϵ)p(\boldsymbol{\epsilon}). The ELBO can be rewritten as:

where z=g(ϵ,ϕ,x)\mathbf{z}=g(\boldsymbol{\epsilon},\boldsymbol{\phi},\mathbf{x}).

This gradient is an unbiased estimator of the exact single-datapoint ELBO gradient; when averaged over noise ϵ∼p(ϵ)\boldsymbol{\epsilon}\sim p(\boldsymbol{\epsilon}), this gradient equals the single-datapoint ELBO gradient:

Computation of the (estimator of) the ELBO requires computation of the density log⁡qϕ(z∣x)\log q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}), given a value of x\mathbf{x}, and given a value of z\mathbf{z} or equivalently ϵ\boldsymbol{\epsilon}. This log-density is a simple computation, as long as we choose the right transformation g()\mathbf{g}().

Note that we typically know the density p(ϵ)p(\boldsymbol{\epsilon}), since this is the density of the chosen noise distribution. As long as g(.)\mathbf{g}(.) is an invertible function, the densities of ϵ\boldsymbol{\epsilon} and z\mathbf{z} are related as:

where the second term is the log of the absolute value of the determinant of the Jacobian matrix (∂z/∂ϵ)(\partial\mathbf{z}/\partial\boldsymbol{\epsilon}):

We call this the log-determinant of the transformation from ϵ\boldsymbol{\epsilon} to z\mathbf{z}. We use the notation log⁡dϕ(x,ϵ)\log d_{\boldsymbol{\phi}}(\mathbf{x},\boldsymbol{\epsilon}) to make explicit that this log-determinant, similar to g()\mathbf{g}(), is a function of x\mathbf{x}, ϵ\boldsymbol{\epsilon} and ϕ\boldsymbol{\phi}. The Jacobian matrix contains all first derivatives of the transformation from ϵ\boldsymbol{\epsilon} to z\mathbf{z}:

As we will show, we can build very flexible transformations g()\mathbf{g}() for which log⁡dϕ(x,ϵ)\log d_{\boldsymbol{\phi}}(\mathbf{x},\boldsymbol{\epsilon}) is simple to compute, resulting in highly flexible inference models qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}).

5 Factorized Gaussian posteriors

A common choice is a simple factorized Gaussian encoder qϕ(z∣x)=N(z;μ,diag(σ2))q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x})=\mathcal{N}(\mathbf{z};\boldsymbol{\mu},\text{diag}(\boldsymbol{\sigma}^{2})):

where N(zi;μi,σi2)\mathcal{N}(z_{i};\mu_{i},\sigma_{i}^{2}) is the PDF of the univariate Gaussian distribution. After reparameterization, we can write:

where ⊙\odot is the element-wise product. The Jacobian of the transformation from ϵ\boldsymbol{\epsilon} to z\mathbf{z} is:

i.e. a diagonal matrix with the elements of σ\boldsymbol{\sigma} on the diagonal. The determinant of a diagonal (or more generally, triangular) matrix is the product of its diagonal terms. The log determinant of the Jacobian is therefore:

when z=g(ϵ,ϕ,x)z=g(\boldsymbol{\epsilon},\phi,\mathbf{x}).

The factorized Gaussian posterior can be extended to a Gaussian with full covariance:

A reparameterization of this distribution is given by:

where L\mathbf{L} is a lower (or upper) triangular matrix, with non-zero entries on the diagonal. The off-diagonal elements define the correlations (covariances) of the elements in z\mathbf{z}.

The reason for this parameterization of the full-covariance Gaussian, is that the Jacobian determinant is remarkably simple. The Jacobian in this case is trivial: ∂z∂ϵ=L\frac{\partial\mathbf{z}}{\partial\boldsymbol{\epsilon}}=\mathbf{L}. Note that the determinant of a triangular matrix is the product of its diagonal elements. Therefore, in this parameterization:

This parameterization corresponds to the Cholesky decomposition Σ=LLT\boldsymbol{\Sigma}=\mathbf{L}\mathbf{L}^{T} of the covariance of z\mathbf{z}:

One way to build a matrix L\mathbf{L} with the desired properties, namely triangularity and non-zero diagonal entries, is by constructing it as follows:

and then proceeding with z=μ+Lϵ\mathbf{z}=\boldsymbol{\mu}+\mathbf{L}\boldsymbol{\epsilon} as described above. Lmask\mathbf{L}_{mask} is a masking matrix with zeros on and above the diagonal, and ones below the diagonal. Note that due to the masking L\mathbf{L}, the Jacobian matrix (∂z/∂ϵ)(\partial\mathbf{z}/\partial\boldsymbol{\epsilon}) is triangular with the values of σ\boldsymbol{\sigma} on the diagonal. The log-determinant is therefore identical to the factorized Gaussian case:

More generally, we can replace z=Lϵ+μ\mathbf{z}=\mathbf{L}\boldsymbol{\epsilon}+\boldsymbol{\mu} with a chain of (differentiable and nonlinear) transformations; as long as the Jacobian of each step in the chain is triangular with non-zero diagonal entries, the log determinant remains simple. This principle is used by inverse autoregressive flow (IAF) as explored by [kingma2016improving] and discussed in chapter 3.

6 Estimation of the Marginal Likelihood

After training a VAE, we can estimate the probability of data under the model using an importance sampling technique, as originally proposed by [rezende2014stochastic]. The marginal likelhood of a datapoint can be written as:

Taking random samples from qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}), a Monte Carlo estimator of this is:

where each z(l)∼qϕ(z∣x)\mathbf{z}^{(l)}\sim q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) is a random sample from the inference model. By making LL large, the approximation becomes a better estimate of the marginal likelihood, and in fact since this is a Monte Carlo estimator, for L→∞L\to\infty this converges to the actual marginal likelihood.

Notice that when setting L=1L=1, this equals the ELBO estimator of the VAE. We can also use the estimator of eq. (2.57) as our objective function; this is the objective used in importance weighted autoencoders [burda2015importance] (IWAE). In that paper, it was also shown that the objective has increasing tightness for increasing value of LL. It was later shown by [cremer2017reinterpreting] that the IWAE objective can be re-interpreted as an ELBO objective with a particular inference model. The downside of these approaches for optimizing a tighter bound, is that importance weighted estimates have notoriously bad scaling properties to high-dimensional latent spaces.

7 Marginal Likelihood and ELBO as KL Divergences

One way to improve the potential tightness of the ELBO, is increasing the flexibility of the generative model. This can be understood through a connection between the ELBO and the KL divergence.

With i.i.d. dataset D\mathcal{D} of size NDN_{\mathcal{D}}, the maximum likelihood criterion is:

where qD(x)q_{\mathcal{D}}(\mathbf{x}) is the empirical (data) distribution, which is a mixture distribution:

where each component qD(i)(x)q_{\mathcal{D}}^{(i)}(\mathbf{x}) typically corresponds to a Dirac delta distribution centered at value x(i)\mathbf{x}^{(i)} in case of continuous data, or a discrete distribution with all probability mass concentrated at value x(i)\mathbf{x}^{(i)} in case of discrete data. The Kullback Leibler (KL) divergence between the data and model distributions, can be rewritten as the negative log-likelihood, plus a constant:

where constant=−H(qD(x))\text{constant}=-\mathcal{H}(q_{\mathcal{D}}(\mathbf{x})). So minimization of the KL divergence above is equivalent to maximization of the data log-likelihood log⁡pθ(D)\log p_{\boldsymbol{\theta}}(\mathcal{D}).

Taking the combination of the empirical data distribution qD(x)q_{\mathcal{D}}(\mathbf{x}) and the inference model, we get a joint distribution over data x\mathbf{x} and latent variables z\mathbf{z}: qD,ϕ(x,z)=qD(x)q(z∣x)q_{\mathcal{D},\boldsymbol{\phi}}(\mathbf{x},\mathbf{z})=q_{\mathcal{D}}(\mathbf{x})q(\mathbf{z}|\mathbf{x}).

The KL divergence of qD,ϕ(x,z)q_{\mathcal{D},\boldsymbol{\phi}}(\mathbf{x},\mathbf{z}) from pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}) can be written as the negative ELBO, plus a constant:

where constant=−H(qD(x))\text{constant}=-\mathcal{H}(q_{\mathcal{D}}(\mathbf{x})). So maximization of the ELBO, is equivalent to the minimization of this KL divergence DKL(qD,ϕ(x,z)∣∣pθ(x,z))D_{KL}(q_{\mathcal{D},\boldsymbol{\phi}}(\mathbf{x},\mathbf{z})||p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z})). The relationship between the ML and ELBO objectives can be summarized in the following simple equation:

One additional perspective is that the ELBO can be viewed as a maximum likelihood objective in an augmented space. For some fixed choice of encoder qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}), we can view the joint distribution pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}) as an augmented empirical distribution over the original data x\mathbf{x} and (stochastic) auxiliary features z\mathbf{z} associated with each datapoint. The model pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}) then defines a joint model over the original data, and the auxiliary features. See figure 2.4.

8 Challenges

In our work, consistent with findings in [bowman2015generating] and [sonderby2016train], we found that stochastic optimization with the unmodified lower bound objective can gets stuck in an undesirable stable equilibrium. At the start of training, the likelihood term log⁡p(x∣z)\log p(\mathbf{x}|\mathbf{z}) is relatively weak, such that an initially attractive state is where q(z∣x)≈p(z)q(\mathbf{z}|\mathbf{x})\approx p(\mathbf{z}), resulting in a stable equilibrium from which it is difficult to escape. The solution proposed in [bowman2015generating] and [sonderby2016train] is to use an optimization schedule where the weights of the latent cost DKL(q(z∣x)∣∣p(z))D_{KL}(q(\mathbf{z}|\mathbf{x})||p(\mathbf{z})) is slowly annealed from to 11 over many epochs.

An alternative proposed in [kingma2016improving] is the method of free bits: a modification of the ELBO objective, that ensures that on average, a certain minimum number of bits of information are encoded per latent variable, or per group of latent variables.

The latent dimensions are divided into the KK groups. We then use the following minibatch objective, which ensures that using less than λ\lambda nats of information per subset jj (on average per minibatch M\mathcal{M}) is not advantageous:

8.2 Blurriness of generative model

Issues with ’blurriness’ can thus can be countered by choosing a sufficiently flexible inference model, and/or a sufficiently flexible generative model. In the next two chapters we will discuss techniques for constructing flexible inference models and flexible generative models.

9 Related prior and concurrent work

Here we briefly discuss relevant literature prior to and concurrent with the work in [kingma2013auto].

The wake-sleep algorithm [hinton1995wake] is another on-line learning method, applicable to the same general class of continuous latent variable models. Like our method, the wake-sleep algorithm employs a recognition model that approximates the true posterior. A drawback of the wake-sleep algorithm is that it requires a concurrent optimization of two objective functions, which together do not correspond to optimization of (a bound of) the marginal likelihood. An advantage of wake-sleep is that it also applies to models with discrete latent variables. Wake-Sleep has the same computational complexity as AEVB per datapoint.

Variational inference has a long history in the field of machine learning. We refer to [wainwright2008graphical] for a comprehensive overview and synthesis of ideas around variational inference for exponential family graphical models. Among other connections, [wainwright2008graphical] shows how various inference algorithms (such as expectation propagation, sum-product, max-product and many others) can be understood as exact or approximate forms of variational inference.

Stochastic variational inference [hoffman2013stochastic] has received increasing interest. [blei2012variational] introduced a control variate schemes to reduce the variance of the score function gradient estimator, and applied the estimator to exponential family approximations of the posterior. In [ranganath2013black] some general methods, e.g. a control variate scheme, were introduced for reducing the variance of the original gradient estimator. In [salimans2013fixedform], a similar reparameterization as in this work was used in an efficient version of a stochastic variational inference algorithm for learning the natural parameters of exponential-family approximating distributions.

In [graves2011practical] a similar estimator of the gradient is introduced; however the estimator of the variance is not an unbiased estimator w.r.t. the ELBO gradient.

The VAE training algorithm exposes a connection between directed probabilistic models (trained with a variational objective) and autoencoders. A connection between linear autoencoders and a certain class of generative linear-Gaussian models has long been known. In [roweis1998algorithms] it was shown that PCA corresponds to the maximum-likelihood (ML) solution of a special case of the linear-Gaussian model with a prior p(z)=N(0,I)p(\mathbf{z})=\mathcal{N}(0,\mathbf{I}) and a conditional distribution p(x∣z)=N(x;Wz,ϵI)p(\mathbf{x}|\mathbf{z})=\mathcal{N}(\mathbf{x};\mathbf{W}\mathbf{z},\epsilon\mathbf{I}), specifically the case with infinitesimally small ϵ\epsilon. In this limiting case, the posterior over the latent variables p(z∣x)p(\mathbf{z}|\mathbf{x}) is a Dirac delta distribution: p(z∣x)=δ(z−W′x)p(\mathbf{z}|\mathbf{x})=\delta(\mathbf{z}-\mathbf{W}^{\prime}\mathbf{x}) where W′=(WTW)−1WT\mathbf{W}^{\prime}=(\mathbf{W}^{T}\mathbf{W})^{-1}\mathbf{W}^{T}, i.e., given W\mathbf{W} and x\mathbf{x} there is no uncertainty about latent variable z\mathbf{z}. [roweis1998algorithms] then introduces an EM-type approach to learning W\mathbf{W}. Much earlier work [bourlard1988auto] showed that optimization of linear autoencoders retrieves the principal components of data, from which it follows that learning linear autoencoders correspond to a specific method for learning the above case of linear-Gaussian probabilistic model of the data. However, this approach using linear autoencoders is limited to linear-Gaussian models, while our approach applies to a much broader class of continuous latent variable models.

When using neural networks for both the inference model and the generative model, the combination forms a type of autoencoder [goodfellow2016deeplearning] with a specific regularization term:

In an analysis of plain autoencoders [vincent2010stacked] it was shown that the training criterion of unregularized autoencoders corresponds to maximization of a lower bound (see the infomax principle [linsker1989application]) of the mutual information between input XX and latent representation ZZ. Maximizing (w.r.t. parameters) of the mutual information is equivalent to maximizing the conditional entropy, which is lower bounded by the expected log-likelihood of the data under the autoencoding model [vincent2010stacked], i.e. the negative reconstruction error. However, it is well known that this reconstruction criterion is in itself not sufficient for learning useful representations [bengio2013representation]. Regularization techniques have been proposed to make autoencoders learn useful representations, such as denoising, contractive and sparse autoencoder variants [bengio2013representation]. The VAE objective contains a regularization term dictated by the variational bound, lacking the usual nuisance regularization hyper-parameter required to learn useful representations. Related are also encoder-decoder architectures such as the predictive sparse decomposition (PSD) [koray-psd-08], from which we drew some inspiration. Also relevant are the recently introduced Generative Stochastic Networks [bengio2014deep] where noisy autoencoders learn the transition operator of a Markov chain that samples from the data distribution. In [salakhutdinov2010efficient] a recognition model was employed for efficient learning with Deep Boltzmann Machines. These methods are targeted at either unnormalized models (i.e. undirected models like Boltzmann machines) or limited to sparse coding models, in contrast to our proposed algorithm for learning a general class of directed probabilistic models.

The proposed DARN method [gregor2014deep], also learns a directed probabilistic model using an autoencoding structure, however their method applies to binary latent variables. In concurrent work, [rezende2014stochastic] also make the connection between autoencoders, directed probabilistic models and stochastic variational inference using the reparameterization trick we describe in [kingma2013auto]. Their work was developed independently of ours and provides an additional perspective on the VAE.

An alternative unbiased stochastic gradient estimator of the ELBO is the score function estimator [kleijnen1996optimization]:

where z∼qϕ(z∣x)\mathbf{z}\sim q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}).

This is also known as the likelihood ratio estimator [glynn1990likelihood, fu2006gradient] and the REINFORCE gradient estimator [williams1992simple]. The method has been successfully used in various methods like neural variational inference [mnih2014neural], black-box variational inference [ranganath2013black], automated variational inference [wingate2013automated], and variational stochastic search [paisley2012variational], often in combination with various novel control variate techniques [glasserman2013monte] for variance reduction. An advantage of the likelihood ratio estimator is its applicability to discrete latent variables.

We do not directly compare to these techniques, since we concern ourselves with continuous latent variables, in which case we have (computationally cheap) access to gradient information ∇zlog⁡pθ(x,z)\nabla_{\mathbf{z}}\log p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}), courtesy of the backpropagation algorithm. The score function estimator solely uses the scalar-valued log⁡pθ(x,z)\log p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}), ignoring the gradient information about the function log⁡pθ(x,z)\log p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}), generally leading to much higher variance. This has been experimentally confirmed by e.g. [kucukelbir2016automatic], which finds that a sophisticated score function estimator requires two orders of magnitude more samples to arrive at the same variance as a reparameterization based estimator.

The difference in efficiency of our proposed reparameterization-based gradient estimator, compared to score function estimators, can intuitively be understood as removing an information bottleneck during the computation of gradients of the ELBO w.r.t. ϕ\boldsymbol{\phi} from current parameters θ\boldsymbol{\theta}: in the latter case, this computation is bottlenecked by the scalar value log⁡pθ(x,z)\log p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}), while in the former case it is bottlenecked by the much wider vector ∇zlog⁡pθ(x,z)\nabla_{\mathbf{z}}\log p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}).

Chapter 3 Beyond Gaussian Posteriors

In this chapter we discuss techniques for improving the flexibility of the inference model qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}). Increasing the flexibility and accuracy of the inference model wel generally improve the tightness of the variational bound (ELBO), bringing it closer the true marginal likelihood objective.

Requirements for the inference model, in order to be able to efficiently optimize the ELBO, are that it is (1) computationally efficient to compute and differentiate its probability density qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}), and (2) computationally efficient to sample from, since both these operations need to be performed for each datapoint in a minibatch at every iteration of optimization. If z\mathbf{z} is high-dimensional and we want to make efficient use of parallel computational resources like GPUs, then parallelizability of these operations across dimensions of z\mathbf{z} is a large factor towards efficiency. This requirement restricts the class of approximate posteriors q(z∣x)q(\mathbf{z}|\mathbf{x}) that are practical to use. In practice this often leads to the use of simple Gaussian posteriors. However, as explained, we also need the density q(z∣x)q(\mathbf{z}|\mathbf{x}) to be sufficiently flexible to match the true posterior p(z∣x)p(\mathbf{z}|\mathbf{x}), in order to arrive at a tight bound.

2 Improving the Flexibility of Inference Models

Here we will review two general techniques for improving the flexibility of approximate posteriors in the context of gradient-based variational inference: auxiliary latent variables, and normalizing flows.

One method for improving the flexibility of inference models, is through the introduction of auxiliary latent variables, as explored by [salimans2015markov], [ranganath2016hierarchical] and [maaloe2016auxiliary].

The methods work by augmenting both the inference model and the generative model with a continuous auxiliary variable, here denoted with u\mathbf{u}.

The inference model defines a distribution over both u\mathbf{u} and and z\mathbf{z}, which can, for example, factorize as:

This inference model augmented with u\mathbf{u}, implicitly defines a potentially powerful implicit marginal distribution:

Likewise, we introduce an additional distribution in the generative model: such that our generative model is now over the joint distribution pθ(x,z,u)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z},\mathbf{u}). This can, for example, factorize as:

The ELBO objective with auxiliary variables, given empirical distribution qD(x)q_{\mathcal{D}}(\mathbf{x}), is then (again) equivalent to minimization of a KL divergence:

Recall that maximization of the original ELBO objective, without auxiliary variables, is equivalent to minimization of DKL(qD,ϕ(x,z)∣∣pθ(x,z))D_{KL}(q_{\mathcal{D},\boldsymbol{\phi}}(\mathbf{x},\mathbf{z})||\allowbreak p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z})), and that maximization of the expected marginal likelihood is equivalent to minimization of DKL(qD,ϕ(x)∣∣pθ(x))D_{KL}(q_{\mathcal{D},\boldsymbol{\phi}}(\mathbf{x})||p_{\boldsymbol{\theta}}(\mathbf{x})).

We can gain additional understanding into the relationship between the objectives, through the following equation:

From this equation it can be seen that in principle, the ELBO gets worse by augmenting the VAE with an auxiliary variable u\mathbf{u}:

The introduction of auxiliary latent variables in the graph, are a special case of VAEs with multiple layers of latent variables, which are discussed in chapter 4. In our experiment with CIFAR-10, we make use of multiple layers of stochastic variables.

2.2 Normalizing Flows

An alternative approach towards flexible approximate posteriors is Normalizing Flow (NF), introduced by [rezende2015variational] in the context of stochastic gradient variational inference. In normalizing flows, we build flexible posterior distributions through an iterative procedure. The general idea is to start off with an initial random variable with a relatively simple distribution with a known (and computationally cheap) probability density function, and then apply a chain of invertible parameterized transformations ft\mathbf{f}_{t}, such that the last iterate zT\mathbf{z}_{T} has a more flexible distributionwhere x\mathbf{x} is the context, such as the value of the datapoint. In case of models with multiple levels of latent variables, the context also includes the value of the previously sampled latent variables.:

The Jacobian of the transformation factorizes:

So its determinant also factorizes as well:

As long as the Jacobian determinant of each of the transformations ft\mathbf{f}_{t} can be computed, we can still compute the probability density function of the last iterate:

[rezende2015variational] experimented with a transformation of the form:

where u\mathbf{u} and w\mathbf{w} are vectors, wT\mathbf{w}^{T} is w\mathbf{w} transposed, bb is a scalar and h(.)h(.) is a nonlinearity, such that uh(wTzt−1+b)\mathbf{u}h(\mathbf{w}^{T}\mathbf{z}_{t-1}+b) can be interpreted as a MLP with a bottleneck hidden layer with a single unit. This flow does not scale well to a high-dimensional latent space: since information goes through the single bottleneck, a long chain of transformations is required to capture high-dimensional dependencies.

3 Inverse Autoregressive Transformations

In order to find a type of normalizing flow that scales well to a high-dimensional space, [kingma2016improving] consider Gaussian versions of autoregressive autoencoders such as MADE [germain2015made] and the PixelCNN [pixelrnn]. Let y\mathbf{y} be a variable modeled by such a model, with some chosen ordering on its elements y={yi}i=1D\mathbf{y}=\{y_{i}\}_{i=1}^{D}. We will use [μ(y),σ(y)][\boldsymbol{\mu}(\mathbf{y}),\boldsymbol{\sigma}(\mathbf{y})] to denote the function of the vector y\mathbf{y}, to the vectors μ\boldsymbol{\mu} and σ\boldsymbol{\sigma}. Due to the autoregressive structure, the Jacobian matrix is triangular with zeros on the diagonal: ∂[μi,σi]/∂yj=\partial[\boldsymbol{\mu}_{i},\boldsymbol{\sigma}_{i}]/\partial\mathbf{y}_{j}= for j≥ij\geq i. The elements [μi(y1:i−1),σi(y1:i−1)][\mu_{i}(\mathbf{y}_{1:i-1}),\sigma_{i}(\mathbf{y}_{1:i-1})] are the predicted mean and standard deviation of the ii-th element of y\mathbf{y}, which are functions of only the previous elements in y\mathbf{y}.

Sampling from such a model is a sequential transformation from a noise vector ϵ∼N(0,I)\boldsymbol{\epsilon}\sim\mathcal{N}(0,\mathbf{I}) to the corresponding vector y\mathbf{y}: y0=μ0+σ0⊙ϵ0y_{0}=\mu_{0}+\sigma_{0}\odot\epsilon_{0}, and for i>0i>0, yi=μi(y1:i−1)+σi(y1:i−1)⋅ϵiy_{i}=\mu_{i}(\mathbf{y}_{1:i-1})+\sigma_{i}(\mathbf{y}_{1:i-1})\cdot\epsilon_{i}. The computation involved in this transformation is clearly proportional to the dimensionality DD. Since variational inference requires sampling from the posterior, such models are not interesting for direct use in such applications. However, the inverse transformation is interesting for normalizing flows. As long as we have σi>0\sigma_{i}>0 for all ii, the sampling transformation above is a one-to-one transformation, and can be inverted:

[kingma2016improving] make two key observations, important for normalizing flows. The first is that this inverse transformation can be parallelized, since (in case of autoregressive autoencoders) computations of the individual elements ϵi\epsilon_{i} do not depend on each other. The vectorized transformation is:

where the subtraction and division are element-wise.

The second key observation, is that this inverse autoregressive operation has a simple Jacobian determinant. Note that due to the autoregressive structure, ∂[μi,σi]/∂yj=\partial[\mu_{i},\sigma_{i}]/\partial y_{j}= for j≥ij\geq i. As a result, the transformation has a lower triangular Jacobian (∂ϵi/∂yj=0\partial\epsilon_{i}/\partial y_{j}=0 for j>ij>i), with a simple diagonal: ∂ϵi/∂yi=1σi\partial\epsilon_{i}/\partial y_{i}=\frac{1}{\sigma_{i}}. The determinant of a lower triangular matrix equals the product of the diagonal terms. As a result, the log-determinant of the Jacobian of the transformation is remarkably simple and straightforward to compute:

The combination of model flexibility, parallelizability across dimensions, and simple log-determinant, makes this transformation interesting for use as a normalizing flow over high-dimensional latent space.

For the following section we will use a slightly different, but equivalently flexible, transformation of the type:

4 Inverse Autoregressive Flow (IAF)

[kingma2016improving] propose inverse autoregressive flow (IAF) based on a chain of transformations that are each equivalent to an inverse autoregressive transformation of eq. (3.19) and eq. (3.21). See algorithm 3 for pseudo-code of an approximate posterior with the proposed flow. We let an initial encoder neural network output μ0\boldsymbol{\mu}_{0} and σ0\boldsymbol{\sigma}_{0}, in addition to an extra output h\mathbf{h}, which serves as an additional input to each subsequent step in the flow. The chain is initialized with a factorized Gaussian qϕ(z0∣x)=N(μ0,diag(σ0)2)q_{\boldsymbol{\phi}}(\mathbf{z}_{0}|\mathbf{x})=\mathcal{N}(\boldsymbol{\mu}_{0},\text{diag}(\boldsymbol{\sigma}_{0})^{2}):

IAF then consists of a chain of TT of the following transformations:

Each step of this flow is an inverse autoregressive transformation of the type of eq. (3.19) and eq. (3.21), and each step uses a separate autoregressive neural network. Following eq. (3.16), the density under the final iterate is:

The flexibility of the distribution of the final iterate ϵT\boldsymbol{\epsilon}_{T}, and its ability to closely fit to the true posterior, increases with the expressivity of the autoregressive models and the depth of the chain. See figure 3.1 for an illustration of the computation.

A numerically stable version, inspired by the LSTM-type update, is where we let the autoregressive network output (mt,st)(\mathbf{m}_{t},\mathbf{s}_{t}), two unconstrained real-valued vectors, and compute ϵt\boldsymbol{\epsilon}_{t} as:

This version is shown in algorithm 3. Note that this is just a particular version of the update of eq. (3.27), so the simple computation of the final log-density of eq. (3.29) still applies.

It was found beneficial for results to parameterize or initialize the parameters of each AutoregressiveNeuralNett\text{AutoregressiveNeuralNet}_{t} such that its outputs st\mathbf{s}_{t} are, before optimization, sufficiently positive, such as close to +1 or +2. This leads to an initial behavior that updates ϵ\boldsymbol{\epsilon} only slightly with each step of IAF. Such a parameterization is known as a ’forget gate bias’ in LSTMs, as investigated by [jozefowicz2015empirical].

It is straightforward to see that a special case of IAF with one step, and a linear autoregressive model, is the fully Gaussian posterior discussed earlier. This transforms a Gaussian variable with diagonal covariance, to one with linear dependencies, i.e. a Gaussian distribution with full covariance.

Autoregressive neural networks form a rich family of nonlinear transformations for IAF. For non-convolutional models, the family of masked autoregressive network introduced in [germain2015made] was used as the autoregressive neural networks. For CIFAR-10 experiments, which benefits more from scaling to high dimensional latent space, the family of convolutional autoregressive autoencoders introduced by [pixelrnn, van2016conditional] was used.

It was found that results improved when reversing the ordering of the variables after each step in the IAF chain. This is a volume-preserving transformation, so the simple form of eq. (3.29) remains unchanged.

5 Related work

As we explained, inverse autoregressive flow (IAF) is a member of the family of normalizing flows, first discussed in [rezende2015variational] in the context of stochastic variational inference. In [rezende2015variational] two specific types of flows are introduced: planar flow (eq. (3.17)) and radial flow. These flows are shown to be effective to problems with a relatively low-dimensional latent space. It is not clear, however, how to scale such flows to much higher-dimensional latent spaces, such as latent spaces of generative models of larger images, and how planar and radial flows can leverage the topology of latent space, as is possible with IAF. Volume-conserving neural architectures were first presented in in [deco1995higher], as a form of nonlinear independent component analysis.

Another type of normalizing flow, introduced by [dinh2014nice] (NICE), uses similar transformations as IAF. In contrast with IAF, NICE was directly applied to the observed variables in a generative model. NICE is type of transformations that updates only half of the variables z1:D/2\mathbf{z}_{1:D/2} per step, adding a vector f(zD/2+1:D)f(\mathbf{z}_{D/2+1:D}) which is a neural network based function of the remaining latent variables zD/2+1:D\mathbf{z}_{D/2+1:D}. Such large blocks have the advantage of computationally cheap inverse transformation, and the disadvantage of typically requiring longer chains. In experiments, [rezende2015variational] found that this type of transformation is generally less powerful than other types of normalizing flow, in experiments with a low-dimensional latent space. Concurrently to our work, NICE was extended to high-dimensional spaces in [dinh2016density] (Real NVP).

A potentially powerful transformation is the Hamiltonian flow used in Hamiltonian Variational Inference [salimans2015markov]. Here, a transformation is generated by simulating the flow of a Hamiltonian system consisting of the latent variables z\mathbf{z}, and a set of auxiliary momentum variables. This type of transformation has the additional benefit that it is guided by the exact posterior distribution, and that it leaves this distribution invariant for small step sizes. Such a transformation could thus take us arbitrarily close to the exact posterior distribution if we can apply it a sufficient number of times. In practice, however, Hamiltonian Variational Inference is very demanding computationally. Also, it requires an auxiliary variational bound to account for the auxiliary variables, which can impede progress if the bound is not sufficiently tight.

An alternative method for increasing the flexibility of variational inference is the introduction of auxiliary latent variables [salimans2015markov, ranganath2016hierarchical, tran2015variational], discussed in 3.2.1, and corresponding auxiliary inference models. Latent variable models with multiple layers of stochastic variables, such as the one used in our experiments, are often equivalent to such auxiliary-variable methods. We combine deep latent variable models with IAF in our experiments, benefiting from both techniques.

Chapter 4 Deeper Generative Models

In the previous chapter we explain advanced strategies for improving inference models. In this chapter, we review strategies for learning deeper generative models, such as inference and learning with multiple latent variables or observed variables, and techniques for improving the flexibility of the generative models pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}).

The generative model pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}), and corresponding inference model qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) can be parameterized as any directed graph. Both x\mathbf{x} and z\mathbf{z} can be composed of multiple variables with some topological ordering. It may not be immediately obvious how to optimize such models in the VAE framework; it is, however, quite straightforward, as we will now explain.

Sampling z∼qϕ(z∣x)\mathbf{z}\sim q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}). In case of multiple latent variables, this means ancestral sampling the latent variables one by one, in topological ordering defined by the inference model’s directed graph. In pseudo-code, the ancestral sampling step looks like:

where Pa(zi)Pa(\mathbf{z}_{i}) are the parents of variable zi\mathbf{z}_{i} in the inference model, which may include x\mathbf{x}. In reparameterized (and differentiable) form, this is:

Evaluating the scalar value (log⁡pθ(x,z)−log⁡qϕ(z∣x))(\log p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z})-\log q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x})) at the resulting sample z\mathbf{z} and datapoint x\mathbf{x}. This scalar is the unbiased stochastic estimate lower bound on log⁡pθ(x)\log p_{\boldsymbol{\theta}}(\mathbf{x}). It is also differentiable and optimizable with SGD.

It should be noted that the choice of latent variables’ topological ordering for the inference model can be different from the choice of ordering for the generative model.

Since the inference model has the data as root node, while the generative model has the data as leaf node, one (in some sense) logical choice would be to let the topological ordering of the latent variables in the inference model be the reverse of the ordering in the generative model.

In multiple works [salimans2016structured, sonderby2016train, kingma2016improving] it has been shown that it can be advantageous to let the generative model and inference model share the topological ordering of latent variables. The two choices of ordering are illustrated in figure 4.1. One advantage of shared ordering, as explained in these works, is that this allows us to easily share parameters between the inference and generative models, leading to faster learning and better solutions.

To see why this might be a good idea, note that the true posterior over the latent variables, is a function of the prior:

Likewise, the posterior of a latent variable given its parents (in the generative model), is:

Optimization of the generative model changes both pθ(zi∣Pa(zi))p_{\boldsymbol{\theta}}(\mathbf{z}_{i}|Pa(\mathbf{z}_{i})) and pθ(x∣zi,Pa(zi))p_{\boldsymbol{\theta}}(\mathbf{x}|\mathbf{z}_{i},Pa(\mathbf{z}_{i})). By coupling the inference model qϕ(zi∣x,Pa(zi))q_{\boldsymbol{\phi}}(\mathbf{z}_{i}|\mathbf{x},Pa(\mathbf{z}_{i})) and prior pθ(zi∣Pa(zi))p_{\boldsymbol{\theta}}(\mathbf{z}_{i}|Pa(\mathbf{z}_{i})), changes in pθ(zi∣Pa(zi))p_{\boldsymbol{\theta}}(\mathbf{z}_{i}|Pa(\mathbf{z}_{i})) can be directly reflected in changes in qϕ(zi∣Pa(zi))q_{\boldsymbol{\phi}}(\mathbf{z}_{i}|Pa(\mathbf{z}_{i})).

This coupling is especially straightforward when pθ(zi∣Pa(zi))p_{\boldsymbol{\theta}}(\mathbf{z}_{i}|Pa(\mathbf{z}_{i})) is Gaussian distributed. The inference model can be directly specified as the product of this Gaussian distribution, with a learned quadratic pseudo-likelihood term:

where ZZ is tractable to compute. This idea is explored by [salimans2016structured] and [sonderby2016train]. In principle this idea could be extended to a more general class of conjugate priors, but no work on this is known at the time of writing.

A less constraining variant, explored by [kingma2016improving], is to simply let the neural network that parameterizes qϕ(zi∣Pa(zi),x)q_{\boldsymbol{\phi}}(\mathbf{z}_{i}|Pa(\mathbf{z}_{i}),\mathbf{x}) be partially specified by a part of the neural network that parameterizes pθ(zi∣Pa(zi))p_{\boldsymbol{\theta}}(\mathbf{z}_{i}|Pa(\mathbf{z}_{i})). In general, we can let the two distributions share parameters. This allows for more complicated posteriors, like normalizing flows or IAF.

2 Alternative methods for increasing expressivity

Typically, especially with large data sets, we wish to choose an expressive class of directed models, such that it can feasibly approximate the true distribution. Popular strategies for specifying expressive models are:

Introduction of latent variables into the directed models, and optimization through (amortized) variational inference, as explained in this work.

Full autoregression: factorization of distributions into univariate (one-dimensional) conditionals, or at least very low-dimensional conditionals (section 4.3).

Specification of distributions through invertible transformations with tractable Jacobian determinant (section 4.4).

Synthesis from fully autoregressive models models is relatively slow, since the length of computation for synthesis from such models is linear in the dimensionality of the data. The length of computation of the log-likelihood of fully autoregressive models does not necesarilly scale with the dimensionality of the data. In this respect, introduction of latent variables for improving expressivity is especially interesting when x\mathbf{x} is very high-dimensional. It is relatively straightforward and computationally attractive, due to parallelizability, to specify directed models over high-dimensional variables where each conditional factorizes into independent distributions. For example, if we let pθ(xj∣Pa(xj))=∏kpθ(xj,k∣Pa(xj))p_{\boldsymbol{\theta}}(\mathbf{x}_{j}|Pa(\mathbf{x}_{j}))=\prod_{k}p_{\boldsymbol{\theta}}(x_{j,k}|Pa(\mathbf{x}_{j})), where each factor is a univariate Gaussian whose means and variance are nonlinear functions (specified by a neural network) of the parents Pa(xj)Pa(\mathbf{x}_{j}), then computations for both synthesis and evaluation of log-likelihood can be fully parallelized across dimensions kk. See [kingma2016improving] for experiments demonstrating a 100x improvement in speed of synthesis.

The best models to date, in terms of attained log-likelihood on test data, employ a combination of the three approaches listed above.

3 Autoregressive Models

A powerful strategy for modeling high-dimensional data is to divide up the high-dimensional observed variables into small constituents (often single dimensional parts, or otherwise just parts with a small number of dimensions), impose a certain ordering, and to model their dependencies as a directed graphical model. The resulting directed graphical model breaks up the joint distribution into a product of a factors:

where DD is the dimensionality of the data. This is known as an autoregressive (AR) model. In case of neural network based autoregressive models, we let the conditional distributions be parameterized with a neural network:

In case of continuous data, autoregressive models can be interpreted as a special case of a more general approach: learning an invertible transformation from the data to a simpler, known distribution such as a Gaussian or Uniform distribution; this approach with invertible transformations is discussed in section 4.4. The techniques of autoregressive models and invertible transformations can be naturally combined with variational autoencoders, and at the time of writing, the best systems use a combination [rezende2015variational, kingma2016improving, gulrajani2016pixelvae].

A disadvantage of autoregressive models, compared to latent-variable models, is that ancestral sampling from autoregressive models is a sequential operation computation of O(D)\mathcal{O}(D) length, i.e. proportional to the dimensionality of the data. Autoregressive models also require choosing a specific ordering of input elements (equation (4.8)). When no single natural one-dimensional ordering exists, like in two-dimensional images, this leads to a model with a somewhat awkward inductive bias.

4 Invertible transformations with tractable Jacobian determinant

In case of continuous data, autoregressive models can be interpreted as a special case of a more general approach: learning an invertible transformation with tractable Jacobian determinant (also called normalizing flow) from the data to a simpler, known distribution such as a Gaussian or Uniform distribution. If we use neural networks for such invertible mappings, this is a powerful and flexible approach towards probabilistic modeling of continuous data and nonlinear independent component analysis [deco1995higher].

Such normalizing flows iteratively update a variable, which is constrained to be of the same dimensionality as the data, to a target distribution. This constraint on the dimensionality of intermediate states of the mapping can make such transformations more challenging to optimize than methods without such constraint. An obvious advantage, on the other hand, is that the likelihood and its gradient are tractable. In [dinh2014nice, dinh2016density], particularly interesting flows (NICE and Real NVP) were introduced, with equal computational cost and depth in both directions, making it both relatively cheap to optimize and to sample from such models. At the time of writing, no such model has yet been demonstrated to lead to the similar performance as purely autoregressive or VAE-based models in terms of data log-likelihood, but this remains an active area of research.

5 Follow-Up Work

Some important applications and motivations for deep generative models and variational autoencoders are:

Representation learning: learning better representations of the data. Some uses of this are:

Data-efficient learning, such as semi-supervised learning

Visualisation of data as low-dimensional manifolds

Artificial creativity: plausible interpolation between data and extrapolation from data.

Here we will now highlight some concrete applications to representation learning and artificial creativity.

In the case of supervised learning, we typically aim to learn a conditional distribution: to predict the distribution over the possible values of a variable, given the value of some another variable. One such problem is that of image classification: given an image, the prediction of a distribution over the possible class labels. Through the yearly ImageNet competion [russakovsky2015imagenet], it has become clear that deep convolutional neural networks [lecun1998gradient, goodfellow2016deeplearning] (CNNs), given a large amount of labeled images, are extraordinarily good at solving the image classification task. Modern versions of CNNs based on residual networks, which is a variant of LSTM-type neural networks [hochreiter1997long], now arguably achieves human-level classification accuracy on this task [he2015delving, he2015deep].

When the number of labeled examples is low, solutions found with purely supervised approaches tend to exhibit poor generalization to new data. In such cases, generative models can be employed as an effective type of regularization. One particular strategy, presented in [kingma2014semi], is to optimize the classification model jointly with a variational autoencoder over the input variables, sharing parameters between the two. The variational autoencoder, in this case, provides an auxiliary objective, improving the data efficiency of the classification solution. Through sharing of statistical strength between modeling problems, this can greatly improve upon the supervised classification error. Techniques based on VAEs are now among state of the art for semi-supervised classification [maaloe2016auxiliary], with on average under 1% classification error in the MNIST classification problem, when trained with only 10 labeled images per class, i.e. when more than 99.8% of the labels in the training set were removed. In concurrent work [rezende2016one], it was shown that VAE-based semi-supervised learning can even do well when only a single sample per class is presented.

A standard supervised approach, GoogLeNet [szegedy2015going], which normally achieves near state-of-the-art performance on the ImageNet validation set, achieves only around 5% top-1 classification accuracy when trained with only 1% of the labeled images, as shown by [pu2016variational]. In contrast, they show that a semi-supervised approach with VAEs achieves around 45% classification accuracy on the same task, when modeling the labels jointly with the labeled and unlabeled input images.

5.2 Understanding of data, and artificial creativity

Generative models with latent spaces allow us to transform the data into a simpler latent space, explore it in that space, and understand it better. A related branch of applications of deep generative models is the synthesis of plausible pseudo-data with certain desirable properties, sometimes coined as artificial creativity.

One example of a recent scientific application of artificial creativity, is shown in [gomez2016automatic]. In this paper, a fairly straightforward VAE is trained on hundreds of thousands of existing chemical structures. The resulting continuous representation (latent space) is subsequently used to perform gradient-based optimization towards certain properties; the method is demonstrated on the design of drug-like molecules and organic light-emitting diodes. See figure 4.2.

A similar approach was used to generating natural-language sentences from a continuous space by [bowman2015generating]. In this paper, it is shown how a VAE can be successfully trained on text. The model is shown to succesfully interpolate between sentences, and for imputation of missing words. See figure 4.3.

In [ravanbakhsh2016enabling], VAEs are applied to simulate observations of distant galaxies. This helps with the calibration of systems that need to indirectly detect the shearing of observations of distant galaxies, caused by weak gravitational lensing in the presence of dark matter between earth and those galaxies. Since the lensing effects are so weak, such systems need to be calibrated with ground-truth images with a known amount of shearing. Since real data is still limited, the proposed solution is to use deep generative models for synthesis of pseudo-data.

A popular application is image (re)synthesis. One can optimize a VAE to form a generative model over images. One can synthesize images from the generative model, but the inference model (or encoder) also allows one to encode real images into a latent space. One can modify the encoding in this latent space, then decode the image back into the observed space. Relatively simple transformations in the observed space, such as linear transformations, often translate into semantically meaningful modifications of the original image. One example, as demonstrated by [white2016sampling], is the modification of images in latent space along a "smile vector" in order to make them more happy, or more sad looking. See figure 4.4 for an example.

5.3 Other relevant follow-up work

We unfortunately do not have space to discuss all follow-up work in depth, but will here highlight a selection of relevant recent work.

In addition to our original publication [kingma2013auto], two later papers have proposed equivalent algorithms [rezende2014stochastic, lazaro2014doubly], where the latter work applies the same reparameterization gradient method to the estimation of parameter posteriors, rather than amortized latent-variable inference.

In the appendix of [kingma2013auto] we proposed to apply the reparameterization gradients to estimation of parameter posteriors. In [blundell2015weight] this method, with a mixture-of-Gaussians prior and named Bayes by Backprop, was used in experiments with some promising early results. In [kingma2015variational] we describe a refined method, the local reparameterization trick, for further decreasing the variance of the gradient estimator, and applied it to estimation of Gaussian parameter posteriors. Further results were presented in [louizos2017bayesian, louizos2017multiplicative, louizos2016structured] with increasingly sophisticated choices of priors and approximate posteriors. In [kingma2015variational, gal2016theoretically], a similar reparameterization was used to analyze Dropout as a Bayesian method, coined Variational Dropout. In [molchanov2017variational] this method was further analyzed and refined. Various papers have applied reparameterization gradients for estimating parameter posteriors, including [fortunato2017bayesian] in the context of recurrent neural networks and [kucukelbir2016automatic] more generally for Bayesian models and in [tran2017deep] for deep probabilistic programming. A Bayesian nonparametric variational family based in the Gaussian Process using reparameterization gradients was proposed in [tran2015variational].

Normalizing flows [rezende2015variational] were proposed as a framework for improving the flexibility of inference models. In [kingma2016improving], the first normalizing flow was proposed that scales well to high-dimensional latent spaces. The same principle was later applied in [papamakarios2017masked] for density estimation, and further refined in [huang2018neural]. Various other flows were proposed in [tomczak2016improving, tomczak2017improving] and [berg2018sylvester].

As an alternative to (or in conjunction with) normalizing flows, one can use auxiliary variables to improve posterior flexibility. This principle was, to the best of our knowledge, first proposed in [salimans2015markov]. In this paper, the principle was used in a combination of variational inference with Hamiltonian Monte Carlo (HMC), with the momentum variables of HMC as auxiliary variables. Auxiliary variables were more elaborately discussed in in [maaloe2016auxiliary] as Auxiliary Deep Generative Models. Similarly, one can use deep models with multiple stochastic layers to improve the variational bound, as demonstrated in [sonderby2016train] and [sonderby2016ladder] as Ladder VAEs.

There has been plenty of follow-up work on gradient variance reduction for the variational parameters of discrete latent variables, as opposed to continuous latent variables for which reparameterization gradients apply. These proposals include NVIL [mnih2014neural], MuProp [gu2015muprop], Variational inference for Monte Carlo objectives [mnih2016variational], the Concrete distribution [maddison2016concrete] and Categorical Reparameterization with Gumbel-Softmax [jang2016categorical].

The ELBO objective can be generalized into an importance-weighted objective, as proposed in [burda2015importance] (Importance-Weighted Autoencoders). This potentially reduces the variance in the gradient, but has not been discussed in-depth here since (as often the case with importance-weighted estimators) it can be difficult to scale to high-dimensional latent spaces. Other objectives have been proposed such as Rényi divergence variational inference [li2016renyi], Generative Moment Matching Networks [li2015generative], objectives based on normalizing such as NICE and RealNVP flows [sohl2015deep, dinh2014nice], black-box α\alpha-divergence minimization [hernandez2016black] and Bi-directional Helmholtz Machines [bornschein2016bidirectional].

Various combinations with adversarial objectives have been proposed. In [makhzani2015adversarial], the "adversarial autoencoder" (AAE) was proposed, a probabilistic autoencoder that uses a generative adversarial network (GAN) [goodfellow2014generative] to perform variational inference. In [dumoulin2016adversarially] Adversarially Learned Inference (ALI) was proposed, which aims to minimize a GAN objective between the joint distributions qϕ(x,z)q_{\boldsymbol{\phi}}(\mathbf{x},\mathbf{z}) and pθ(x,z)p_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{z}). Other hybrids have been proposed as well [larsen2015autoencoding, brock2016neural, hsu2017voice].

One of the most prominent, and most difficult, applications of generative models is image modeling. In [kulkarni2015deep] (Deep convolutional inverse graphics network), a convolutional VAE was applied to modeling images with some success, building on work by [dosovitskiy2015learning] proposing convolutional networks for image synthesis. In [gregor2015draw] (DRAW), an attention mechanism was combined with a recurrent inference model and recurrent generative model for image synthesis. This approach was further extended in [gregor2016towards] (Towards Conceptual Compression) with convolutional networks, scalable to larger images, and applied to image compression. In [kingma2016improving], deep convolutional inference models and generative models were also applied to images. Furthermore, [gulrajani2016pixelvae] (PixelVAE) and [chen2016variational] (Variational Lossy Autoencoder) combined convolutional VAEs with the PixelCNN model [pixelrnn, van2016conditional]. Methods and VAE architectures for controlled image generation from attributes or text were studied in [kingma2014semi, yan2016attribute2image, mansimov2015generating, brock2016neural, white2016sampling]. Predicting the color of pixels based on a grayscale image is another promising application [deshpande2016learning]. The application to semi-supervised learning has been studied in [kingma2014semi, pu2016variational, xu2017variational] among other work.

Another prominent application of VAEs is modeling of text and or sequential data [bayer2014learning, bowman2015generating, serban2016hierarchical, johnson2016composing, karl2016deep, fraccaro2016sequential, miao2016neural, semeniuta2017hybrid, zhao2017learning, yang2017improved, hu2017controllable]. VAEs have also been applied to speech and handwriting [chung2015recurrent]. Sequential models typically use recurrent neural networks, such as LSTMs [hochreiter1997long], as encoder and/or decoder. When modeling sequences, the validity of a sequence can sometimes be constrained by a context-free grammar. In this case, incorporation of the grammar in VAEs can lead to better models, as shown in [kusner2017grammar] (Grammar VAEs), and applied to modeling molecules in textual representations.

Since VAEs can transform discrete observation spaces to continuous latent-variable spaces with approximately known marginals, they are interesting for use in model-based control [watter2015embed, pritzel2017neural]. In [heess2015learning] (Stochastic Value Gradients) it was shown that the re-parameterization of the observed variables, together with an observation model, can be used to compute novel forms of policy gradients. Variational inference and reparameterization gradients have also been used for variational information maximisation for intrinsically motivated reinforcement learning [mohamed2015variational] and VIME [houthooft2016vime] for improved exploration. Variational autoencoders have also been used as components in models that perform iterative reasoning about objects in a scene [eslami2016attend].

In [higgins2016beta] (β\beta-VAE) it was proposed to strengthen the contribution of DKL(qϕ(z∣x)∣∣pθ(z))D_{KL}(q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x})||p_{\boldsymbol{\theta}}(\mathbf{z})), thus restricting the information flow through the latent space, which was shown to improve disentanglement of latent factors, further studied in [chen2018isolating].

Other applications include modeling of graphs [kipf2016variational] (Variational Graph Autoencoders), learning of 3D structure from images [rezende2016unsupervised], one-shot learning [rezende2016one], learning nonlinear state space models [krishnan2017structured], voice conversion from non-parallel corpora [hsu2016voice], discrimination-aware (fair) representations [louizos2015variational] and transfer learning [edwards2016towards].

The reparameterization gradient estimator discussed in this work has been extended in various directions [ruiz2016generalized], including acceptance-rejection sampling algorithms [naesseth2017reparameterization]. The gradient variance can in some cases be reduced by ’carving up the ELBO’ [hoffman2016elbo, roeder2017sticking] and using a modified gradient estimator. A second-order gradient estimator has also been proposed in [fan2015fast].

All in all, this remains an actively researched area with frequently exciting developments.

Chapter 5 Conclusion

Directed probabilistic models form an important aspect of modern artificial intelligence. Such models can be made incredibly flexible by parameterizing the conditional distributions with differentiable deep neural networks.

Optimization of such models towards the maximum likelihood objective is straightforward in the fully-observed case. However, one is often more interested in flexible models with latent variables, such as deep latent-variable models, or Bayesian models with random parameters. In both cases one needs to perform approximate posterior estimation for which variational inference (VI) methods are suitable. In VI, inference is cast as an optimization problem over newly introduced variational parameters, typically optimized towards the ELBO, a lower bound on the model evidence, or marginal likelihood of the data. Existing methods for such posterior inference were either relatively inefficient, or not applicable to models with neural networks as components. Our main contribution is a framework for efficient and scalable gradient-based variational posterior inference and approximate maximum likelihood learning.

In this work we describe the variational autoencoder (VAE) and some of its extensions. A VAE is a combination of a deep latent-variable model (DLVM) with continuous latent variables, and an associated inference model. The DLVM is a type of generative model over the data. The inference model, also called encoder or recognition model, approximates the posterior distribution of the latent variables of the generative model. Both the generative model and the inference model are directed graphical models that are wholly or partially parameterized by deep neural networks. The parameters of the models, including the parameters of the neural networks such as the weights and biases, are jointly optimized by performing stochastic gradient ascent on the so-called evidence lower bound (ELBO). The ELBO is a lower bound on the marginal likelihood of the data, also called the variational lower bound. Stochastic gradients, necessary for performing SGD, are obtained through a basic reparameterization trick. The VAE framework is now a commonly used tool for various applications of probabilistic modeling and artificial creativity, and basic implementations are available in most major deep learning software libraries.

For learning flexible inference models, we proposed inverse autoregressive flows (IAF), a type of normalizing flow that allows scaling to high-dimensional latent spaces. An interesting direction for further exploration is comparison with transformations with computationally cheap inverses, such as NICE [dinh2014nice] and Real NVP [dinh2016density]. Application of such transformations in the VAE framework can potentially lead to relatively simple VAEs with a combination of powerful posteriors, priors and decoders. Such architectures can potentially rival or surpass purely autoregressive architectures [van2016conditional], while allowing much faster synthesis.

The proposed VAE framework remains the only framework in the literature that allows for both discrete and continuous observed variables, allows for efficient amortized latent-variable inference and fast synthesis, and which can produce close to state-of-the-art performance in terms of the log-likelihood of data.

Appendix A Appendix

A.1.2 Definitions

A.1.3 Distributions

We overload the notation of distributions (e.g. p(x)=N(x;μ,Σ)p(\mathbf{x})=\mathcal{N}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma})) with two meanings: (1) a distribution from which we can sample, and (2) the probability density function (PDF) of that distribution.

A.1.4 Bayesian Inference

Let p(θ)p(\theta) be a chosen marginal distribution over its parameters θ\theta, called a prior distribution. Let D\mathcal{D} be observed data, p(D∣θ)≡pθ(D)p(\mathcal{D}|\theta)\equiv p_{\boldsymbol{\theta}}(\mathcal{D}) be the probability assigned to the data under the model with parameters θ\theta. Recall the chain rule in probability:

Simply re-arranging terms above, the posterior distribution over the parameters θ\theta, taking into account the data D\mathcal{D}, is:

where the proportionality (∝\propto) holds since p(D)p(\mathcal{D}) is a constant that is not dependent on parameters θ\theta. The formula above is known as Bayes’ rule, a fundamental formula in machine learning and statistics, and is of special importance to this work.

A principal application of Bayes’ rule is that it allows us to make predictions about future data x′\mathbf{x}^{\prime}, that are optimal as long as the prior p(θ)p(\theta) and model class pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}) are correct:

A.2 Alternative methods for learning in DLVMs

From a Bayesian perspective, we can improve upon the maximum likelihood objective through maximum a posteriori (MAP) estimation, which maximizes the log-posterior w.r.t. θ\theta. With i.i.d. data D\mathcal{D}, this is:

The prior p(θ)p(\theta) in equation (A.5) has diminishing effect for increasingly large NN. For this reason, in case of optimization with large datasets, we often choose to simply use the maximum likelihood criterion by omitting the prior from the objective, which is numerically equivalent to setting p(θ)=constantp(\theta)=\text{constant}.

A.2.2 Variational EM with local variational parameters

Expectation Maximization (EM) is a general strategy for learning parameters in partially observed models [dempster1977em]. See section A.2.3 for a discussion of EM using MCMC. The method can be explained as coordinate ascent on the ELBO [neal1998em]. In case of of i.i.d. data, traditional variational EM methods estimate local variational parameters ϕ(i)\boldsymbol{\phi}^{(i)}, i.e. a separate set of variational parameters per datapoint ii in the dataset. In contrast, VAEs employ a strategy with global variational parameters.

EM starts out with some (random) initial choice of θ\boldsymbol{\theta} and ϕ(1:N)\boldsymbol{\phi}^{(1:N)}. It then iteratively applies updates:

until convergence. Why does this work? Note that at the E-step:

so the EE-step, sensibly, minimizes the KL divergence of qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) from the true posterior.

Secondly, note that if qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}) equals pθ(z∣x)p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x}), the ELBO equals the marginal likelihood, but that for any choice of qϕ(z∣x)q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x}), the MM-step optimizes a bound on the marginal likelihood. The tightness of this bound is defined by DKL(qϕ(z∣x)∣∣pθ(z∣x))D_{KL}(q_{\boldsymbol{\phi}}(\mathbf{z}|\mathbf{x})||p_{\boldsymbol{\theta}}(\mathbf{z}|\mathbf{x})).

A.2.3 MCMC-EM

Another Bayesian approach towards optimizing the likelihood pθ(x)p_{\boldsymbol{\theta}}(\mathbf{x}) with DLVMs is Expectation Maximization (EM) with Markov Chain Monte Carlo (MCMC). In case of MCMC, the posterior is approximated by a mixture of a set of approximately i.i.d. samples from the posterior, acquired by running a Markov chain. Note that posterior gradients in DLVMs are relatively affordable to compute by differentiating the log-joint distribution w.r.t. z\mathbf{z}:

One version of MCMC which uses such posterior for relatively fast convergence, is Hamiltonian MCMC [neal2011mcmc]. A disadvantage of this approach is the requirement for running an independent MCMC chain per datapoint.

A.3 Stochastic Gradient Descent

We work with directed models where the objective per datapoint is scalar, and due to the differentiability of neural networks that compose them, the objective is differentiable w.r.t. its parameters θ\theta. Due to the remarkable efficiency of reverse-mode automatic differentiation (also known as the backpropagation algorithm [rumelhart1988learning]), the value and gradient (i.e. the vector of partial derivatives) of differentiable scalar objectives can be computed with equal time complexity. In SGD, we iteratively update parameters θ\theta: