The Variational Gaussian Process

Introduction

Variational inference is a powerful tool for approximate posterior inference. The idea is to posit a family of distributions over the latent variables and then find the member of that family closest to the posterior. Originally developed in the 1990s (Hinton & Van Camp, 1993; Waterhouse et al., 1996; Jordan et al., 1999), variational inference has enjoyed renewed interest around developing scalable optimization for large datasets (Hoffman et al., 2013), deriving generic strategies for easily fitting many models (Ranganath et al., 2014), and applying neural networks as a flexible parametric family of approximations (Kingma & Welling, 2014; Rezende et al., 2014). This research has been particularly successful for computing with deep Bayesian models (Neal, 1990; Ranganath et al., 2015a), which require inference of a complex posterior distribution (Hinton et al., 2006).

Classical variational inference typically uses the mean-field family, where each latent variable is independent and governed by its own variational distribution. While convenient, the strong independence limits learning deep representations of data. Newer research aims toward richer families that allow dependencies among the latent variables. One way to introduce dependence is to consider the variational family itself as a model of the latent variables (Lawrence, 2000; Ranganath et al., 2015b). These variational models naturally extend to Bayesian hierarchies, which retain the mean-field “likelihood” but introduce dependence through variational latent variables.

In this paper we develop a powerful new variational model—the variational Gaussian process (vgp). The vgp is a Bayesian nonparametric variational model; its complexity grows efficiently and towards any distribution, adapting to the inference problem at hand. We highlight three main contributions of this work:

We prove a universal approximation theorem: under certain conditions, the vgp can capture any continuous posterior distribution—it is a variational family that can be specified to be as expressive as needed.

We derive an efficient stochastic optimization algorithm for variational inference with the vgp. Our algorithm can be used in a wide class of models. Inference with the vgp is a black box variational method (Ranganath et al., 2014).

We study the vgp on standard benchmarks for unsupervised learning, applying it to perform inference in deep latent Gaussian models (Rezende et al., 2014) and DRAW (Gregor et al., 2015), a latent attention model. For both models, we report the best results to date.

Technical summary. Generative models hypothesize a distribution of observations $\boldsymbol{\mathbf{x}}$ and latent variables $\boldsymbol{\mathbf{z}}$, $p(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{z}})$. Variational inference posits a family of the latent variables $q(\boldsymbol{\mathbf{z}};\boldsymbol{\mathbf{\lambda}})$ and tries to find the variational parameters $\boldsymbol{\mathbf{\lambda}}$ that are closest in KL divergence to the posterior. When we use a variational model, $q(\boldsymbol{\mathbf{z}};\boldsymbol{\mathbf{\lambda}})$ itself might contain variational latent variables; these are implicitly marginalized out in the variational family (Ranganath et al., 2015b).

The vgp is a flexible variational model. It draw inputs from a simple distribution, warps those inputs through a non-linear mapping, and then uses the output of the mapping to govern the distribution of the latent variables $\boldsymbol{\mathbf{z}}$. The non-linear mapping is itself a random variable, constructed from a Gaussian process. The vgp is inspired by ideas from both the Gaussian process latent variable model (Lawrence, 2005) and Gaussian process regression (Rasmussen & Williams, 2006).

The variational parameters of the vgp are the kernel parameters for the Gaussian process and a set of variational data, which are input-output pairs. The variational data is crucial: it anchors the non-linear mappings at given inputs and outputs. It is through these parameters that the vgp learns complex representations. Finally, given data $\boldsymbol{\mathbf{x}}$, we use stochastic optimization to find the variational parameters that minimize the KL divergence to the model posterior.

Variational Gaussian Process

Variational models introduce latent variables to the variational family, providing a rich construction for posterior approximation (Ranganath et al., 2015b). Here we introduce the variational Gaussian process (vgp), a Bayesian nonparametric variational model that is based on the Gaussian process. The Gaussian process (gp) provides a class of latent variables that lets us capture downstream distributions with varying complexity.

We first review variational models and Gaussian processes. We then outline the mechanics of the vgp and prove that it is a universal approximator.

Let $p(\boldsymbol{\mathbf{z}}\,|\,\boldsymbol{\mathbf{x}})$ denote a posterior distribution over $d$ latent variables $\boldsymbol{\mathbf{z}}=(z_{1},\ldots,z_{d})$ conditioned on a data set $\boldsymbol{\mathbf{x}}$. For a family of distributions $q(\boldsymbol{\mathbf{z}};\boldsymbol{\mathbf{\lambda}})$ parameterized by $\boldsymbol{\mathbf{\lambda}}$, variational inference seeks to minimize the divergence $\operatorname{KL}(q(\boldsymbol{\mathbf{z}};\boldsymbol{\mathbf{\lambda}})\,\|\,p(\boldsymbol{\mathbf{z}}\,|\,\boldsymbol{\mathbf{x}}))$. This is equivalent to maximizing the evidence lower bound (elbo) (Wainwright & Jordan, 2008). The elbo can be written as a sum of the expected log likelihood of the data and the KL divergence between the variational distribution and the prior,

Traditionally, variational inference considers a tractable family of distributions with analytic forms for its density. A common specification is a fully factorized distribution $\prod_{i}q(z_{i};\lambda_{i})$, also known as the mean-field family. While mean-field families lead to efficient computation, they limit the expressiveness of the approximation.

The variational family of distributions can be interpreted as a model of the latent variables $\boldsymbol{\mathbf{z}}$, and it can be made richer by introducing new latent variables. Hierarchical variational models consider distributions specified by a variational prior of the mean-field parameters $q(\boldsymbol{\mathbf{\lambda}};\boldsymbol{\mathbf{\theta}})$ and a factorized “likelihood” $\prod_{i}q(z_{i}\mid\lambda_{i})$. This specifies the variational model,

which is governed by prior hyperparameters $\boldsymbol{\mathbf{\theta}}$. Hierarchical variational models are richer than classical variational families—their expressiveness is determined by the complexity of the prior $q(\boldsymbol{\mathbf{\lambda}})$. Many expressive variational approximations can be viewed under this construct (Saul & Jordan, 1996; Jaakkola & Jordan, 1998; Rezende & Mohamed, 2015; Tran et al., 2015).

2 Gaussian Processes

with parameters $\boldsymbol{\mathbf{\theta}}=(\sigma^{2}_{\textsc{ard}},\omega_{1},\ldots,\omega_{c})$. The weights $\omega_{j}$ tune the importance of each dimension. They can be driven to zero during inference, leading to automatic dimensionality reduction.

Given data $\mathcal{D}$, the conditional distribution of the gp forms a distribution over mappings which interpolate between input-output pairs,

Here, $\mathbf{K}_{\boldsymbol{\mathbf{\xi}}s}$ denotes the covariance function $k(\boldsymbol{\mathbf{\xi}},\mathbf{s})$ for an input $\boldsymbol{\mathbf{\xi}}$ and over all data inputs $\mathbf{s}_{n}$, and $\mathbf{t}_{i}$ represents the $i^{th}$ output dimension.

3 Variational Gaussian Processes

We describe the variational Gaussian process (vgp), a Bayesian nonparametric variational model that admits arbitrary structures to match posterior distributions. The vgp generates $\boldsymbol{\mathbf{z}}$ by generating latent inputs, warping them with random non-linear mappings, and using the warped inputs as parameters to a mean-field distribution. The random mappings are drawn conditional on “variational data,” which are variational parameters. We will show that the vgp enables samples from the mean-field to follow arbitrarily complex posteriors.

The vgp specifies the following generative process for posterior latent variables $\boldsymbol{\mathbf{z}}$:

Draw approximate posterior samples $\boldsymbol{\mathbf{z}}\in\operatorname{supp}(p)$: $\boldsymbol{\mathbf{z}}=(z_{1},\ldots,z_{d})\sim\prod_{i=1}^{d}q(f_{i}(\boldsymbol{\mathbf{\xi}})).$

Figure 1 displays a graphical model for the vgp. Here, $\mathcal{D}=\{(\mathbf{s}_{n},\mathbf{t}_{n})\}_{n=1}^{m}$ represents variational data, comprising input-output pairs that are parameters to the variational distribution. Marginalizing over all latent inputs and non-linear mappings, the vgp is

The vgp is parameterized by kernel hyperparameters $\boldsymbol{\mathbf{\theta}}$ and variational data.

As a variational model, the vgp forms an infinite ensemble of mean-field distributions. A mean-field distribution is given in the first term of the integrand above. It is conditional on a fixed function $f(\cdot)$ and input $\boldsymbol{\mathbf{\xi}}$; the $d$ outputs $f_{i}(\boldsymbol{\mathbf{\xi}})=\lambda_{i}$ are the mean-field’s parameters. The vgp is a form of a hierarchical variational model (Eq.2) (Ranganath et al., 2015b). It places a continuous Bayesian nonparametric prior over mean-field parameters.

Unlike the mean-field, the vgp can capture correlation between the latent variables. The reason is that it evaluates the $d$ independent gp draws at the same latent input $\boldsymbol{\mathbf{\xi}}$. This induces correlation between their outputs, the mean-field parameters, and thus also correlation between the latent variables. Further, the vgp is flexible. The complex non-linear mappings drawn from the gp allow it to capture complex discrete and continuous posteriors.

We emphasize that the vgp needs variational data. Unlike typical gp regression, there are no observed data available to learn a distribution over non-linear mappings of the latent variables $\boldsymbol{\mathbf{z}}$. Thus the "data" are variational parameters that appear in the conditional distribution of $f$ in Eq.4. They anchor the random non-linear mappings at certain input-ouput pairs. When optimizing the vgp, the learned variational data enables finds a distribution of the latent variables that closely follows the posterior.

4 Universal approximation theorem

To understand the capacity of the vgp for representing complex posterior distributions, we analyze the role of the Gaussian process. For simplicity, suppose the latent variables $\boldsymbol{\mathbf{z}}$ are real-valued, and the vgp treats the output of the function draws from the gp as posterior samples. Consider the optimal function $f^{*}$, which is the transformation such that when we draw $\boldsymbol{\mathbf{\xi}}\sim\mathcal{N}(\boldsymbol{\mathbf{0}},\boldsymbol{\mathbf{I}})$ and calculate $\boldsymbol{\mathbf{z}}=f^{*}(\boldsymbol{\mathbf{\xi}})$, the resulting distribution of $\boldsymbol{\mathbf{z}}$ is the posterior distribution.

An explicit construction of $f^{*}$ exists if the dimension of the latent input $\boldsymbol{\mathbf{\xi}}$ is equal to the number of latent variables. Let $P^{-1}$ denote the inverse posterior CDF and $\Phi$ the standard normal CDF. Using techniques common in copula literature (Nelsen, 2006), the optimal function is

Imagine generating samples $\boldsymbol{\mathbf{z}}$ using this function. For latent input $\boldsymbol{\mathbf{\xi}}\sim\mathcal{N}(\boldsymbol{\mathbf{0}},\boldsymbol{\mathbf{I}})$, the standard normal CDF $\Phi$ applies the probability integral transform: it squashes $\xi_{i}$ such that its output $u_{i}=\Phi(\xi_{i})$ is uniformly distributed on $$. The inverse posterior CDF then transforms the uniform random variables$P^{-1}(u_{1},\ldots,u_{d})=\boldsymbol{\mathbf{z}}$ to follow the posterior. The function produces exact posterior samples.

In the vgp, the random function interpolates the values in the variational data, which are optimized to minimize the KL divergence. Thus, during inference, the distribution of the gp learns to concentrate around this optimal function. This perspective provides intuition behind the following result.

Let $q(\boldsymbol{\mathbf{z}};\boldsymbol{\mathbf{\theta}},\mathcal{D})$ denote the variational Gaussian process. Consider a posterior distribution $p(\boldsymbol{\mathbf{z}}\,|\,\boldsymbol{\mathbf{x}})$ with a finite number of latent variables and continuous quantile function (inverse CDF). There exists a sequence of parameters $(\boldsymbol{\mathbf{\theta}}_{k},\mathcal{D}_{k})$ such that

See Appendix B for a proof. Theorem 1 states that any posterior distribution with strictly positive density can be represented by a vgp. Thus the vgp is a flexible model for learning posterior distributions.

Black box inference

We derive an algorithm for black box inference over a wide class of generative models.

The original elbo (Eq.1) is analytically intractable due to the log density, $\log q_{\textsc{vgp}}(\boldsymbol{\mathbf{z}})$ (Eq.5). To address this, we present a tractable variational objective inspired by auto-encoders (Kingma & Welling, 2014).

A tractable lower bound to the model evidence $\log p(\boldsymbol{\mathbf{x}})$ can be derived by subtracting an expected KL divergence term from the elbo,

where $r(\boldsymbol{\mathbf{\xi}},f\,|\,\boldsymbol{\mathbf{z}})$ is an auxiliary model (we describe $r$ in the next subsection). Various versions of this objective have been considered in the literature (Jaakkola & Jordan, 1998; Agakov & Barber, 2004), and it has been recently revisited by Salimans et al. (2015) and Ranganath et al. (2015b). We perform variational inference in the posterior latent variable space, minimizing $\operatorname{KL}(q\|p)$ to learn the variational model; for this to occur we perform auxiliary inference in the variational latent variable space, minimizing $\operatorname{KL}(q\|r)$ to learn an auxiliary model. See Figure 2.

Unlike previous approaches, we rewrite this variational objective to connect to auto-encoders:

where the KL divergences are now taken over tractable distributions (see Appendix C). In auto-encoder parlance, we maximize the expected negative reconstruction error, regularized by two terms: an expected divergence between the variational model and the original model’s prior, and an expected divergence between the auxiliary model and the variational model’s prior. This is simply a nested instantiation of the variational auto-encoder bound (Kingma & Welling, 2014): a divergence between the inference model and a prior is taken as regularizers on both the posterior and variational spaces. This interpretation justifies the previously proposed bound for variational models; as we shall see, it also enables lower variance gradients during stochastic optimization.

2 Auto-encoding variational models

An inference network provide a flexible parameterization of approximating distributions as used in Helmholtz machines (Hinton & Zemel, 1994), deep Boltzmann machines (Salakhutdinov & Larochelle, 2010), and variational auto-encoders (Kingma & Welling, 2014; Rezende et al., 2014). It replaces local variational parameters with global parameters coming from a neural network. For latent variables $\boldsymbol{\mathbf{z}}_{n}$ (which correspond to a data point $\boldsymbol{\mathbf{x}}_{n}$), an inference network specifies a neural network which takes $\boldsymbol{\mathbf{x}}_{n}$ as input and its local variational parameters $\boldsymbol{\mathbf{\lambda}}_{n}$ as output. This amortizes inference by only defining a set of global parameters.

To auto-encode the vgp we specify inference networks to parameterize both the variational and auxiliary models:

3 Stochastic optimization

We maximize the variational objective $\widetilde{\mathcal{L}}(\boldsymbol{\mathbf{\theta}},\boldsymbol{\mathbf{\phi}})$ over both $\boldsymbol{\mathbf{\theta}}$ and $\boldsymbol{\mathbf{\phi}}$, where $\boldsymbol{\mathbf{\theta}}$ newly denotes both the kernel hyperparameters and the inference network’s parameters for the vgp, and $\boldsymbol{\mathbf{\phi}}$ denotes the inference network’s parameters for the auxiliary model. Following the black box methods, we write the gradient as an expectation and apply stochastic approximations (Robbins & Monro, 1951), sampling from the variational model and evaluating noisy gradients.

First, we reduce variance of the stochastic gradients by analytically deriving any tractable expectations. The KL divergence between $q(\boldsymbol{\mathbf{z}}\,|\,f(\boldsymbol{\mathbf{\xi}}))$ and $p(\boldsymbol{\mathbf{z}})$ is commonly used to reduce variance in traditional variational auto-encoders: it is analytic for deep generative models such as the deep latent Gaussian model (Rezende et al., 2014) and deep recurrent attentive writer (Gregor et al., 2015). The KL divergence between $r(f\,|\,\boldsymbol{\mathbf{\xi}},\boldsymbol{\mathbf{z}})$ and $q(f\,|\,\boldsymbol{\mathbf{\xi}})$ is analytic as the distributions are both Gaussian. The difference $\log q(\boldsymbol{\mathbf{\xi}})-\log r(\boldsymbol{\mathbf{\xi}}\,|\,\boldsymbol{\mathbf{z}})$ is simply a difference of Gaussian log densities. See Appendix C for more details.

To derive black box gradients, we can first reparameterize the vgp, separating noise generation of samples from the parameters in its generative process (Kingma & Welling, 2014; Rezende et al., 2014). The gp easily enables reparameterization: for latent inputs $\boldsymbol{\mathbf{\xi}}\sim\mathcal{N}(\boldsymbol{\mathbf{0}},\boldsymbol{\mathbf{I}})$, the transformation $\mathbf{f}(\boldsymbol{\mathbf{\xi}};\boldsymbol{\mathbf{\theta}})=\boldsymbol{\mathbf{L}}\boldsymbol{\mathbf{\xi}}+\mathbf{K}_{\boldsymbol{\mathbf{\xi}}s}\mathbf{K}_{ss}^{-1}\mathbf{t}_{i}$ is a location-scale transform, where $\boldsymbol{\mathbf{L}}\boldsymbol{\mathbf{L}}^{\top}=\mathbf{K}_{\boldsymbol{\mathbf{\xi}}\boldsymbol{\mathbf{\xi}}}-\mathbf{K}_{\boldsymbol{\mathbf{\xi}}s}\mathbf{K}_{ss}^{-1}\mathbf{K}_{\boldsymbol{\mathbf{\xi}}s}^{\top}$. This is equivalent to evaluating $\boldsymbol{\mathbf{\xi}}$ with a random mapping from the gp. Suppose the mean-field $q(\boldsymbol{\mathbf{z}}\,|\,f(\boldsymbol{\mathbf{\xi}}))$ is also reparameterizable, and let $\boldsymbol{\mathbf{\epsilon}}\sim w$ such that $\boldsymbol{\mathbf{z}}(\boldsymbol{\mathbf{\epsilon}};\mathbf{f})$ is a function of $\boldsymbol{\mathbf{\xi}}$ whose output $\boldsymbol{\mathbf{z}}\sim q(\boldsymbol{\mathbf{z}}\,|\,f(\boldsymbol{\mathbf{\xi}}))$. This two-level reparameterization is equivalent to the generative process for $\boldsymbol{\mathbf{z}}$ outlined in Section 2.3.

We now rewrite the variational objective as

Eq.7 enables gradients to move inside the expectations and backpropagate over the nested reparameterization. Thus we can take unbiased stochastic gradients, which exhibit low variance due to both the analytic KL terms and reparameterization. The gradients are derived in Appendix D, including the case when the first KL is analytically intractable.

We outline the method in Algorithm 1. For massive data, we apply subsampling on $\boldsymbol{\mathbf{x}}$ (Hoffman et al., 2013). For gradients of the model log-likelihood, we employ convenient differentiation tools such as those in Stan and Theano (Carpenter et al., 2015; Bergstra et al., 2010). For non-differentiable latent variables $\boldsymbol{\mathbf{z}}$, or mean-field distributions without efficient reparameterizations, we apply the black box gradient estimator from Ranganath et al. (2014) to take gradients of the inner expectation.

4 Computational and storage complexity

The algorithm has $\mathcal{O}(d+m^{3}+LH^{2})$ complexity, where $d$ is the number of latent variables, $m$ is the size of the variational data, and $L$ is the number of layers of the neural networks with $H$ the average hidden layer size. In particular, the algorithm is linear in the number of latent variables, which is competitive with other variational inference methods. The number of variational and auxiliary parameters has $\mathcal{O}(c+LH)$ complexity; this complexity comes from storing the kernel hyperparameters and the neural network parameters.

Unlike most gp literature, we require no low rank constraints, such as the use of inducing variables for scalable computation (Quiñonero-Candela & Rasmussen, 2005). The variational data serve a similar purpose, but inducing variables reduce the rank of a (fixed) kernel matrix; the variational data directly determine the kernel matrix and thus the kernel matrix is not fixed. Although we haven’t found it necessary in practice, see Appendix E for scaling the size of variational data.

Related work

Recently, there has been interest in applying parametric transformations for approximate inference. Parametric transformations of random variables induce a density in the transformed space, with a Jacobian determinant that accounts for how the transformation warps unit volumes. Kucukelbir et al. (2016) consider this viewpoint for automating inference, in which they posit a transformation from the standard normal to a possibly constrained latent variable space. In general, however, calculating the Jacobian determinant incurs a costly $\mathcal{O}(d^{3})$ complexity, cubic in the number of latent variables. Dinh et al. (2015) consider volume-preserving transformations which avoid calculating Jacobian determinants. Salimans et al. (2015) consider volume-preserving transformations defined by Markov transition operators. Rezende & Mohamed (2015) consider a slightly broader class of parametric transformations, with Jacobian determinants having at most $\mathcal{O}(d)$ complexity.

Instead of specifying a parametric class of mappings, the vgp posits a Bayesian nonparametric prior over all continuous mappings. The vgp can recover a certain class of parametric transformations by using kernels which induce a prior over that class. In the context of the vgp, the gp is an infinitely wide feedforward network which warps latent inputs to mean-field parameters. Thus, the vgp offers complete flexibility on the space of mappings—there are no restrictions such as invertibility or linear complexity—and is fully Bayesian. Further, it is a hierarchical variational model, using the gp as a variational prior over mean-field parameters (Ranganath et al., 2015b). This enables inference over both discrete and continuous latent variable models.

In addition to its flexibility over parametric methods, the vgp is more computationally efficient. Parametric methods must consider transformations with Jacobian determinants of at most $\mathcal{O}(d)$ complexity. This restricts the flexibility of the mapping and therefore the flexibility of the variational model (Rezende & Mohamed, 2015). In comparison, the distribution of outputs using a gp prior does not require any Jacobian determinants (following Eq.4); instead it requires auxiliary inference for inferring variational latent variables (which is fast). Further, unlike discrete Bayesian nonparametric priors such as an infinite mixture of mean-field distributions, the gp enables black box inference with lower variance gradients—it applies a location-scale transform for reparameterization and has analytically tractable KL terms.

Transformations, which convert samples from a tractable distribution to the posterior, is a classic technique in Bayesian inference. It was first studied in Monte Carlo methods, where it is core to the development of methods such as path sampling, annealed importance sampling, and sequential Monte Carlo (Gelman & Meng, 1998; Neal, 1998; Chopin, 2002). These methods can be recast as specifying a discretized mapping $f_{t}$ for times $t_{0}<\ldots<t_{k}$, such that for draws $\boldsymbol{\mathbf{\xi}}$ from the tractable distribution, $f_{t_{0}}(\boldsymbol{\mathbf{\xi}})$ outputs the same samples and $f_{t_{k}}(\boldsymbol{\mathbf{\xi}})$ outputs exact samples following the posterior. By applying the sequence in various forms, the transformation bridges the tractable distribution to the posterior. Specifying a good transformation—termed “schedule” in the literature—is crucial to the efficiency of these methods. Rather than specify it explicitly, the vgp adaptively learns this transformation and avoids discretization.

Limiting the vgp in various ways recovers well-known probability models as variational approximations. Specifically, we recover the discrete mixture of mean-field distributions (Bishop et al., 1998; Jaakkola & Jordan, 1998). We also recover a form of factor analysis (Tipping & Bishop, 1999) in the variational space. Mathematical details are in Appendix A.

Experiments

Following standard benchmarks for variational inference in deep learning, we learn generative models of images. In particular, we learn the deep latent Gaussian model (dlgm) (Rezende et al., 2014), a layered hierarchy of Gaussian random variables following neural network architecures, and the recently proposed Deep Recurrent Attentive Writer (draw) (Gregor et al., 2015), a latent attention model that iteratively constructs complex images using a recurrent architecture and a sequence of variational auto-encoders (Kingma & Welling, 2014).

For the learning rate we apply a version of RMSProp (Tieleman & Hinton, 2012), in which we scale the value with a decaying schedule $1/t^{1/2+\epsilon}$ for $\epsilon>0$. We fix the size of variational data to be $500$ across all experiments and set the latent input dimension equal to the number of latent variables.

The binarized MNIST data set (Salakhutdinov & Murray, 2008) consists of 28x28 pixel images with binary-valued outcomes. Training a dlgm, we apply two stochastic layers of 100 random variables and 50 random variables respectively, and in-between each stochastic layer is a deterministic layer with 100 units using tanh nonlinearities. We apply mean-field Gaussian distributions for the stochastic layers and a Bernoulli likelihood. We train the vgp to learn the dlgm for the cases of one stochastic layer and two stochastic layers.

For draw (Gregor et al., 2015), we augment the mean-field Gaussian distribution originally used to generate the latent samples at each time step with the vgp, as it places a complex variational prior over its parameters. The encoding recurrent neural network now outputs variational data (used for the variational model) as well as mean-field Gaussian parameters (used for the auxiliary model). We use the same architecture hyperparameters as in Gregor et al. (2015).

After training we evaluate test set log likelihood, which are lower bounds on the true value. See Table 1 which reports both approximations and lower bounds of $\log p(\boldsymbol{\mathbf{x}})$ for various methods. The vgp achieves the highest known results on log-likelihood using draw, reporting a value of -79.88 compared to the original highest of -80.97. The vgp also achieves the highest known results among the class of non-structure exploiting models using the dlgm, with a value of -81.32 compared to the previous best of -82.90 reported by Burda et al. (2016).

2 Sketch

As a demonstration of the vgp’s complexity for learning representations, we also examine the Sketch data set (Eitz et al., 2012). It consists of 20,000 human sketches equally distributed over 250 object categories. We partition it into 18,000 training examples and 2,000 test examples. We fix the architecture of draw to have a 2x2 read window, 5x5 write attention window, and 64 glimpses—these values were selected using a coarse grid search and choosing the set which lead to the best training log likelihood. For inference we use the original auto-encoder version as well as the augmented version with the vgp.

See Table 2. draw with the vgp achieves a significantly better lower bound, performing better than the original version which has seen state-of-the-art success in many computer vision tasks. (Until the results presented here, the results from the original draw were the best reported performance for this data set.). Moreover, the model inferred using the vgp is able to generate more complex images than the original version—it not only performs better but maintains higher visual fidelity.

Discussion

We present the variational Gaussian process (vgp), a variational model which adapts its shape to match complex posterior distributions. The vgp draws samples from a tractable distribution, and posits a Bayesian nonparametric prior over transformations from the tractable distribution to mean-field parameters. The vgp learns the transformations from the space of all continuous mappings—it is a universal approximator and finds good posterior approximations via optimization.

In future work the vgp will be explored for application in Monte Carlo methods, where it may be an efficient proposal distribution for importance sampling and sequential Monte Carlo. An important avenue of research is also to characterize local optima inherent to the objective function. Such analysis will improve our understanding of the limits of the optimization procedure and thus the limits of variational inference.

We thank David Duvenaud, Alp Kucukelbir, Ryan Giordano, and the anonymous reviewers for their helpful comments. This work is supported by NSF IIS-0745520, IIS-1247664, IIS-1009542, ONR N00014-11-1-0651, DARPA FA8750-14-2-0009, N66001-15-C-4032, Facebook, Adobe, Amazon, and the Seibel and John Templeton Foundations.

References

Appendix A Special cases of the variational Gaussian process

We now analyze two special cases of the vgp: by limiting its generative process in various ways, we recover well-known models. This provides intuition behind the vgp’s complexity. In Section 4 we show many recently proposed models can also be viewed as special cases of the vgp.

A mixture of mean-field distributions is a vgp without a kernel.

A discrete mixture of mean-field distributions (Bishop et al., 1998; Jaakkola & Jordan, 1998; Lawrence, 2000) is a classically studied variational model with dependencies between latent variables. Instead of a mapping which interpolates between inputs of the variational data, suppose the vgp simply performs nearest-neighbors for a latent input $\boldsymbol{\mathbf{\xi}}$—selecting the output $t_{n}$ tied to the nearest variational input $s_{n}$. This induces a multinomial distribution of outputs, which samples one of the variational outputs’ mean-field parameters.Formally, given variational input-output pairs $\{(\mathbf{s}_{n},\mathbf{t}_{n})\}$, the nearest-neighbor function is defined as $f(\boldsymbol{\mathbf{\xi}})=\mathbf{t}_{j}$, such that $\|\boldsymbol{\mathbf{\xi}}-\mathbf{s}_{j}\|<\|\boldsymbol{\mathbf{\xi}}-\mathbf{s}_{k}\|$ for all $k$. Then the output’s distribution is multinomial with probabilities $P(f(\boldsymbol{\mathbf{\xi}})=\mathbf{t}_{j})$, proportional to areas of the partitioned nearest-neighbor space. Thus, with a gp prior that interpolates between inputs, the vgp can be seen as a kernel density smoothing of the nearest-neighbor function.

Variational factor analysis is a vgp with linear kernel and no variational data.

Consider factor analysis (Tipping & Bishop, 1999) in the variational space: For simplicity, we avoid discussion of the vgp’s underlying mean-field distribution, i.e., we specify each mean-field factor to be a degenerate point mass at its parameter value.

Marginalizing over the latent inputs induces linear dependence in $\boldsymbol{\mathbf{z}}$, $q(\boldsymbol{\mathbf{z}};\boldsymbol{\mathbf{w}})=\mathcal{N}(\boldsymbol{\mathbf{z}};\boldsymbol{\mathbf{0}},\boldsymbol{\mathbf{w}}\boldsymbol{\mathbf{w}}^{\top})$. Consider the dual interpretation

with $q(\boldsymbol{\mathbf{z}}\,|\,\boldsymbol{\mathbf{\xi}})=\mathcal{N}(\boldsymbol{\mathbf{z}};\boldsymbol{\mathbf{0}},\boldsymbol{\mathbf{\xi}}\boldsymbol{\mathbf{\xi}}^{\top})$. The maximum likelihood estimate of $\boldsymbol{\mathbf{w}}$ in factor analysis is the maximum a posteriori estimate of $\boldsymbol{\mathbf{\xi}}$ in the gp formulation. More generally, use of a non-linear kernel induces non-linear dependence in $\boldsymbol{\mathbf{z}}$. Learning the set of kernel hyperparameters $\boldsymbol{\mathbf{\theta}}$ thus learns the set capturing the most variation in its latent embedding of $\boldsymbol{\mathbf{z}}$ (Lawrence, 2005).

Appendix B Proof of Theorem 1

Let $q(\boldsymbol{\mathbf{z}};\boldsymbol{\mathbf{\theta}},\mathcal{D})$ denote the variational Gaussian process. Consider a posterior distribution $p(\boldsymbol{\mathbf{z}}\,|\,\boldsymbol{\mathbf{x}})$ with a finite number of latent variables and continuous quantile function (inverse CDF). There exists a sequence of parameters $(\boldsymbol{\mathbf{\theta}}_{k},\mathcal{D}_{k})$ such that

Let the mean-field distribution be given by degenerate delta distributions

Let the size of the latent input be equivalent to the number of latent variables $c=d$ and fix $\sigma^{2}_{\textsc{ard}}=1$ and $\boldsymbol{\mathbf{\omega}}_{j}=1$. Furthermore for simplicity, we assume that $\boldsymbol{\mathbf{\xi}}$ is drawn uniformly on the $d$-dimensional hypercube. Then as explained in Section 2.4, if we let $P^{-1}$ denote the inverse posterior cumulative distribution function, the optimal $f$ denoted $f^{*}$ such that

Define ${\cal O}_{k}$ to be the set of points $j/2^{k}$ for $j=0$ to $2^{k}$, and define ${\cal S}_{k}$ to be the $d$-dimensional product of ${\cal O}_{k}$. Let $\mathcal{D}_{k}$ be the set containing the pairs $(s_{i},f^{*}(s_{i}))$, for each element $s_{i}$ in ${\cal S}_{k}$. Denote $f^{k}$ as the gp mapping conditioned on the dataset $\mathcal{D}_{k}$, this random mapping satisfies $f^{k}(s_{i})=f^{*}(s_{i})$ for all $s_{i}\in{\cal S}_{k}$ by the noise free prediction property of Gaussian processes (Rasmussen & Williams, 2006). Then by continuity, as $k\to\infty$, $f^{k}$ converges to $f^{*}$. ∎

A broad condition under which the quantile function of a distribution is continuous is if that distribution has positive density with respect to the Lebesgue measure.

The rate of convergence for finite sizes of the variational data can be studied via posterior contraction rates for gps under random covariates (Van Der Vaart & Van Zanten, 2011). Only an additional assumption using stronger continuity conditions for the posterior quantile and the use of Matern covariance functions is required for the theory to be applicable in the variational setting.

Appendix C Variational objective

We derive the tractable lower bound to the model evidence $\log p(\boldsymbol{\mathbf{x}})$ presented in Eq.6. To do this, we first penalize the elbo with an expected KL term,

We can combine all terms into the expectations as follows:

where we apply the product rule $q(\boldsymbol{\mathbf{z}})q(\boldsymbol{\mathbf{\xi}},f\,|\,\boldsymbol{\mathbf{z}})=q(\boldsymbol{\mathbf{z}}\,|\,f(\boldsymbol{\mathbf{\xi}}))q(\boldsymbol{\mathbf{\xi}},f)$. Recombining terms as KL divergences, and written with parameters $(\boldsymbol{\mathbf{\theta}},\boldsymbol{\mathbf{\phi}})$, this recovers the auto-encoded variational objective in Section 3:

The KL divergence between the mean-field $q(\boldsymbol{\mathbf{z}}\,|\,f(\boldsymbol{\mathbf{\xi}}))$ and the model prior $p(\boldsymbol{\mathbf{z}})$ is analytically tractable for certain popular models. For example, in the deep latent Gaussian model (Rezende et al., 2014) and draw (Gregor et al., 2015), both the mean-field distribution and model prior are Gaussian, leading to an analytic KL term: for Gaussian random variables of dimension $d$,

In general, when the KL is intractable, we combine the KL term with the reconstruction term, and maximize the variational objective

We expect that this experiences slightly higher variance in the stochastic gradients during optimization.

Appendix D Gradients of the variational objective

We derive gradients for the variational objective (Eq.7). This follows trivially by backpropagation:

where we assume the KL terms are analytically written from Appendix C and gradients are propagated similarly through their computational graph. In practice, we need only be careful about the expectations, and the gradients of the functions written above are taken care of with automatic differentiation tools.

We also derive gradients for the general variational bound of Eq.8—it assumes that the first KL term, measuring the divergence between $q$ and the prior for $p$, is not necessarily tractable. Following the reparameterizations described in Section 3.3, this variational objective can be rewritten as

We calculate gradients by backpropagating over the nested reparameterizations:

Appendix E Scaling the size of variational data

If massive sizes of variational data are required, e.g., when its cubic complexity due to inversion of a $m\times m$ matrix becomes the bottleneck during computation, we can scale it further. Consider fixing the variational inputs to lie on a grid. For stationary kernels, this allows us to exploit Toeplitz structure for fast $m\times m$ matrix inversion. In particular, one can embed the Toeplitz matrix into a circulant matrix and apply conjugate gradient combined with fast Fourier transforms in order to compute inverse-matrix vector products in $\mathcal{O}(m\log m)$ computation and $\mathcal{O}(m)$ storage (Cunningham et al., 2008). For product kernels, we can further exploit Kronecker structure to allow fast $m\times m$ matrix inversion in $\mathcal{O}(Pm^{1+1/P})$ operations and $\mathcal{O}(Pm^{2/P})$ storage, where $P>1$ is the number of kernel products (Osborne, 2010). The ard kernel specifically leads to $\mathcal{O}(cm^{1+1/c})$ complexity, which is linear in $m$.