Automatic Differentiation Variational Inference

Introduction

We develop an automatic method that derives variational inference algorithms for complex probabilistic models. We implement our method in Stan, a probabilistic programming system that lets a user specify a model in an intuitive programming language and then compiles that model into an inference executable. Our method enables fast inference with large datasets on an expansive class of probabilistic models.

Figure 14 gives an example. Say we want to analyze how people navigate a city by car. We have a dataset of all the taxi rides taken over the course of a year: 1.7 million trajectories. To explore patterns in this data, we propose a mixture model with an unknown number of components. This is a non-conjugate model that we seek to fit to a large dataset. Previously, we would have to manually derive an inference algorithm that scales to large data. In our method, we write a Stan program and compile it; we can then fit the model in minutes and analyze the results with ease.

The context of this research is the field of probabilistic modeling, which has emerged as a powerful language for customized data analysis. Probabilistic modeling lets us express our assumptions about data in a formal mathematical way, and then derive algorithms that use those assumptions to compute about an observed dataset. It has had an impact on myriad applications in both statistics and machine learning, including natural language processing, speech recognition, computer vision, population genetics, and computational neuroscience.

Probabilistic modeling leads to a natural research cycle. A scientist first uses her domain knowledge to posit a simple model that includes latent variables; then, she uses an inference algorithm to infer those variables from her data; next, she analyzes her results and identifies where the model works and where it falls short; last, she refines the model and repeats the process. When we cycle through these steps, we find expressive, interpretable, and useful models (Gelman et al., 2013; Blei, 2014). One of the broad goals of machine learning is to make this process easy.

Looping around this cycle, however, is not easy. The data we study are often large and complex; accordingly, we want to propose rich probabilistic models and scale them up. But using such models requires complex algorithms that are difficult to derive, implement, and scale. The bottleneck is this computation that precludes the scientist from taking full advantage of the probabilistic modeling cycle.

This problem motivates the important ideas of probabilistic programming and automated inference. Probabilistic programming allows a user to write a probability model as a computer program and then compile that program into an efficient inference executable. Automated inference is the backbone of such a system—it inputs a probability model, expressed as a program, and outputs an efficient algorithm for computing with it. Previous approaches to automatic inference have mainly relied on Markov chain Monte Carlo (mcmc) algorithms. The results have been successful, but automated mcmc is too slow for many real-world applications.

We approach the problem through variational inference, a faster alternative to mcmc that has been used in many large-scale problems (Blei et al., 2016). Though it is a promising method, developing a variational inference algorithm still requires tedious model-specific derivations and implementation; it has not seen widespread use in probabilistic programming. Here we automate the process of deriving scalable variational inference algorithms. We build on recent ideas in so-called “black-box” variational inference to leverage strengths of probabilistic programming systems, namely the ability to transform the space of latent variables and to automate derivatives of the joint distribution. The result, called automatic differentiation variational inference (advi), provides an automated solution to variational inference: the inputs are a probabilistic model and a dataset; the outputs are posterior inferences about the model’s latent variables.This paper extends the method presented in (Kucukelbir et al., 2015). We implemented and deployed advi as part of Stan, a probabilistic programming system (Stan Development Team, 2015).

advi in Stan resolves the computational bottleneck of the probabilistic modeling cycle. A scientist can easily propose a probabilistic model, analyze a large dataset, and revise the model, without worrying about computation. advi enables this cycle by providing automated and scalable variational inference for an expansive class of models. Sections 3 and 4 present ten probabilistic modeling examples, including a progressive analysis of 1.7 million taxi trajectories.

Variational inference turns the task of computing a posterior into an optimization problem. We posit a parameterized family of distributions $q(\mathbf{\boldsymbol{\theta}})\in\mathcal{Q}$ and then find the member of that family that minimizes the Kullback-Leibler (kl) divergence to the exact posterior. Traditionally, using a variational inference algorithm requires the painstaking work of developing and implementing a custom optimization routine: specifying a variational family appropriate to the model, computing the corresponding objective function, taking derivatives, and running a gradient-based or coordinate-ascent optimization.

advi solves this problem automatically. The user specifies the model, expressed as a program, and advi automatically generates a corresponding variational algorithm. The idea is to first automatically transform the inference problem into a common space and then to solve the variational optimization. Solving the problem in this common space solves variational inference for all models in a large class. In more detail, advi follows these steps.

advi recasts the gradient of the variational objective function as an expectation over $q$. This involves the gradient of the log joint with respect to the latent variables $\nabla_{\mathbf{\boldsymbol{\theta}}}\log p(\mathbf{\boldsymbol{x}},\mathbf{\boldsymbol{\theta}})$. Expressing the gradient as an expectation opens the door to Monte Carlo methods for approximating it (Robert and Casella, 1999).

advi further reparameterizes the gradient in terms of a standard Gaussian. To do this, it uses another transformation, this time within the variational family. This second transformation enables advi to efficiently compute Monte Carlo approximations—it needs only to sample from a standard Gaussian (Kingma and Welling, 2014; Rezende et al., 2014).

advi uses noisy gradients to optimize the variational distribution (Robbins and Monro, 1951). An adaptively tuned step-size sequence provides good convergence in practice (Bottou, 2012).

We developed advi in the Stan system, which gives us two important types of automatic computation around probabilistic models. First, Stan provides a library of transformations—ways to convert a variety of constrained latent variables (e.g., positive reals) to be unconstrained, without changing the underlying joint distribution. Stan’s library of transformations helps us with step 1 above. Second, Stan implements automatic differentiation to calculate $\nabla_{\mathbf{\boldsymbol{\theta}}}\log p(\mathbf{\boldsymbol{x}},\mathbf{\boldsymbol{\theta}})$ (Carpenter et al., 2015; Baydin et al., 2015). These derivatives are crucial in step 2, when computing the gradient of the advi objective.

Organization of paper. Section 2 develops the recipe that makes advi. We expose the details of each of the steps above and present a concrete algorithm. Section 3 studies the properties of advi. We explore its accuracy, its stochastic nature, and its sensitivity to transformations. Section 4 applies advi to an array of probability models. We compare its speed to mcmc sampling techniques and present a case study using a dataset with millions of observations. Section 5 concludes the paper with a discussion.

Automatic Differentiation Variational Inference

Automatic differentiation variational inference (advi) offers a recipe for automating the computations involved in variational inference. The strategy is as follows: transform the latent variables of the model into a common space, choose a variational approximation in the common space, and use generic computational techniques to solve the variational problem.

Many posterior densities are not tractable to compute; their normalizing constants lack analytic (closed-form) solutions. Thus we often seek to approximate the posterior. advi approximates the posterior of differentiable probability models. Members of this class of models have continuous latent variables $\mathbf{\boldsymbol{\theta}}$ and a gradient of the log-joint with respect to them, $\nabla_{\mathbf{\boldsymbol{\theta}}}\log p(\mathbf{\boldsymbol{x}},\mathbf{\boldsymbol{\theta}})$. The gradient is valid within the support of the prior

where $K$ is the dimension of the latent variable space. This support set is important: it will play a role later in the paper. We make no assumptions about conjugacy, either full (Diaconis et al., 1979) or conditional (Hoffman et al., 2013).

Many machine learning models are differentiable. For example: linear and logistic regression, matrix factorization with continuous or discrete observations, linear dynamical systems, and Gaussian processes. (See Table 1.) At first blush, the restriction to continuous random variables may seem to leave out other common machine learning models, such as mixture models, hidden Markov models, and topic models. However, marginalizing out the discrete variables in the likelihoods of these models renders them differentiable. Marginalization is not tractable for all models, such as the Ising model, sigmoid belief network, and (untruncated) Bayesian nonparametric models, such as Dirichlet process mixtures (Blei et al., 2006). These are not differentiable probability models.

2 Variational Inference

Variational inference turns approximate posterior inference into an optimization problem (Wainwright and Jordan, 2008; Blei et al., 2016). Consider a family of approximating densities of the latent variables $q(\mathbf{\boldsymbol{\theta}}\,;\,\mathbf{\boldsymbol{\phi}})$, parameterized by a vector $\mathbf{\boldsymbol{\phi}}\in\mathbf{\boldsymbol{\Phi}}$. Variational inference (vi) finds the parameters that minimize the kl divergence to the posterior,

The optimized $q(\mathbf{\boldsymbol{\theta}}\,;\,\mathbf{\boldsymbol{\phi}}^{*})$ then serves as an approximation to the posterior.

The kl divergence lacks an analytic form because it involves the posterior. Instead we maximize the evidence lower bound (elbo)

The first term is an expectation of the joint density under the approximation, and the second is the entropy of the variational density. The elbo is equal to the negative kl divergence up to the constant $\log p(\mathbf{\boldsymbol{x}})$. Maximizing the elbo minimizes the kl divergence (Jordan et al., 1999; Bishop, 2006).

We explicitly include this constraint because we have not specified the form of the variational approximation; we must ensure that $q(\mathbf{\boldsymbol{\theta}}\,;\,\mathbf{\boldsymbol{\phi}})$ stays within the support of the posterior.

Our recipe for automating vi. The traditional way of solving Equation (3) is difficult. We begin by choosing a variational family $q(\mathbf{\boldsymbol{\theta}}\,;\,\mathbf{\boldsymbol{\phi}})$ that, by definition, satisfies the support matching constraint. We compute the expectations in the elbo, either analytically or through approximation. We then decide on a strategy to maximize the elbo. For instance, we might use coordinate ascent by iteratively updating the components of $\mathbf{\boldsymbol{\phi}}$. Or, we might follow gradients of the elbo with respect to $\mathbf{\boldsymbol{\phi}}$ while staying within $\mathbf{\boldsymbol{\Phi}}$. Finally, we implement, test, and debug software that performs the above. Each step requires expert thought and analysis in the service of a single algorithm for a single model.

In contrast, our approach allows the scientist to define any differentiable probability model, and we automate the process of developing a corresponding vi algorithm. Our recipe for automating vi has three ingredients. First, we automatically transform the support of the latent variables $\mathbf{\boldsymbol{\theta}}$ to the real coordinate space (Section 2.3); this lets us choose from a variety of variational distributions $q$ without worrying about the support matching constraint (Section 2.4). Second, we compute the elbo for any model using Monte Carlo (mc) integration, which only requires being able to sample from the variational distribution (Section 2.5). Third, we employ stochastic gradient ascent to maximize the elbo and use automatic differentiation to compute gradients without any user input (Section 2.6). With these tools, we can develop a generic method that automatically solves the variational optimization problem for a large class of models.

3 Automatic Transformation of Constrained Variables

Define a one-to-one differentiable function

and identify the transformed variables as $\mathbf{\boldsymbol{\zeta}}=T(\mathbf{\boldsymbol{\theta}})$. The transformed joint density $p(\mathbf{\boldsymbol{x}},\mathbf{\boldsymbol{\zeta}})$ is a function of $\mathbf{\boldsymbol{\zeta}}$; it has the representation

where $p(\mathbf{\boldsymbol{x}},\mathbf{\boldsymbol{\theta}}=T^{-1}(\mathbf{\boldsymbol{\zeta}}))$ is the joint density in the original latent variable space, and $J_{T^{-1}}(\mathbf{\boldsymbol{\zeta}})$ is the Jacobian of the inverse of $T$. Transformations of continuous probability densities require a Jacobian; it accounts for how the transformation warps unit volumes and ensures that the transformed density integrates to one (Olive, 2014). (See Appendix A.)

As we describe in the introduction, we implement our algorithm in Stan (Stan Development Team, 2015). Stan maintains a library of transformations and their corresponding Jacobians.Stan provides various transformations for upper and lower bounds, simplex and ordered vectors, and structured matrices such as covariance matrices and Cholesky factors. With Stan, we can automatically transforms the joint density of any differentiable probability model to one with real-valued latent variables. (See Figure 2.)

4 Variational Approximations in Real Coordinate Space

Mean-field Gaussian. One option is to posit a factorized (mean-field) Gaussian variational approximation

Full-rank Gaussian. Another option is to posit a full-rank Gaussian variational approximation

The full-rank Gaussian generalizes the mean-field Gaussian approximation. The off-diagonal terms in the covariance matrix $\mathbf{\boldsymbol{\Sigma}}$ capture posterior correlations across latent random variables.This is a form of structured mean-field variational inference (Wainwright and Jordan, 2008; Barber, 2012). This leads to a more accurate posterior approximation than the mean-field Gaussian; however, it comes at a computational cost. Various low-rank approximations to the covariance matrix reduce this cost, yet limit its ability to model complex posterior correlations (Seeger, 2010; Challis and Barber, 2013).

The choice of a Gaussian. Choosing a Gaussian distribution may call to mind the Laplace approximation technique, where a second-order Taylor expansion around the maximum-a-posteriori estimate gives a Gaussian approximation to the posterior. However, using a Gaussian variational approximation is not equivalent to the Laplace approximation (Opper and Archambeau, 2009). Our approach is distinct in another way: the posterior approximation in the original latent variable space is non-Gaussian.

The implicit variational density. The transformation $T$ from Equation (4) maps the support of the latent variables to the real coordinate space. Thus, its inverse $T^{-1}$ maps back to the support of the latent variables. This implicitly defines the variational approximation in the original latent variable space as q\left(T(\mathbf{\boldsymbol{\theta}})\,;\,\mathbf{\boldsymbol{\phi}}\right)\big{|}\det J_{T}(\mathbf{\boldsymbol{\theta}})\big{|}. The transformation ensures that the support of this approximation is always bounded by that of the posterior in the original latent variable space.

Sensitivity to $T$. There are many ways to transform the support a variable to the real coordinate space. The form of the transformation directly affects the shape of the variational approximation in the original latent variable space. We study sensitivity to the choice of transformation in Section 3.3.

5 The Variational Problem in Real Coordinate Space

Here is the story so far. We began with a differentiable probability model $p(\mathbf{\boldsymbol{x}},\mathbf{\boldsymbol{\theta}})$. We transformed the latent variables into $\mathbf{\boldsymbol{\zeta}}$, which live in the real coordinate space. We defined variational approximations in the transformed space. Now, we consider the variational optimization problem.

Write the variational objective function, the elbo, in real coordinate space as

Now, we can freely optimize the elbo in the real coordinate space without worrying about the support matching constraint. The optimization problem from Equation (3) becomes

where the parameter vector $\mathbf{\boldsymbol{\phi}}$ lives in some appropriately dimensioned real coordinate space. This is an unconstrained optimization problem that we can solve using gradient ascent. Traditionally, this would require manual computation of gradients. Instead, we develop a stochastic gradient ascent algorithm that uses automatic differentiation to compute gradients and mc integration to approximate expectations.

We cannot directly use automatic differentiation on the elbo. This is because the elbo involves an unknown expectation. However, we can automatically differentiate the functions inside the expectation. (The model $p$ and transformation $T$ are both easy to represent as computer functions (Baydin et al., 2015).) To apply automatic differentiation, we want to push the gradient operation inside the expectation. To this end, we employ one final transformation: elliptical standardizationAlso known as a “coordinate transformation” (Rezende et al., 2014), an “invertible transformation” (Titsias and Lázaro-Gredilla, 2014), and the “re-parameterization trick” (Kingma and Welling, 2014). (Härdle and Simar, 2012).

Elliptical standardization. Consider a transformation $S_{\mathbf{\boldsymbol{\phi}}}$ that absorbs the variational parameters $\mathbf{\boldsymbol{\phi}}$; this converts the Gaussian variational approximation into a standard Gaussian. In the mean-field case, the standardization is $\mathbf{\boldsymbol{\eta}}=S_{\mathbf{\boldsymbol{\phi}}}(\mathbf{\boldsymbol{\zeta}})=\text{diag}\left(\exp\left(\mathbf{\boldsymbol{\omega}}\right)\right)^{-1}(\mathbf{\boldsymbol{\zeta}}-\mathbf{\boldsymbol{\mu}})$. In the full-rank Gaussian, the standardization is $\mathbf{\boldsymbol{\eta}}=S_{\mathbf{\boldsymbol{\phi}}}(\mathbf{\boldsymbol{\zeta}})=\mathbf{\boldsymbol{L}}^{-1}(\mathbf{\boldsymbol{\zeta}}-\mathbf{\boldsymbol{\mu}})$.

In both cases, the standardization encapsulates the variational parameters; in return it gives a fixed variational density

The standardization transforms the variational problem from Equation (5) into

The expectation is now in terms of a standard Gaussian density. The Jacobian of elliptical standardization evaluates to one, because the Gaussian distribution is a member of the location-scale family: standardizing a Gaussian gives another Gaussian distribution. (See Appendix A.)

We do not need to transform the entropy term as it does not depend on the model or the transformation; we have a simple analytic form for the entropy of a Gaussian and its gradient. We implement these once and reuse for all models.

6 Stochastic Optimization

We now reach the final step: stochastic optimization of the variational objective function.

Computing gradients. Since the expectation is no longer dependent on $\mathbf{\boldsymbol{\phi}}$, we can directly calculate its gradient. Push the gradient inside the expectation and apply the chain rule to get

We obtain gradients with respect to $\mathbf{\boldsymbol{\omega}}$ (mean-field) and $\mathbf{\boldsymbol{L}}$ (full-rank) in a similar fashion

We can now compute the gradients inside the expectation with automatic differentiation. The only thing left is the expectation. mc integration provides a simple approximation: draw samples from the standard Gaussian and evaluate the empirical mean of the gradients within the expectation (Appendix D). In practice a single sample suffices. (We study this in detail in Section 3.2 and in the experiments in Section 4.)

This gives noisy unbiased gradients of the elbo for any differentiable probability model. We can now use these gradients in a stochastic optimization routine to automate variational inference.

Stochastic gradient ascent. Equipped with noisy unbiased gradients of the elbo, advi implements stochastic gradient ascent (Algorithm 1). This algorithm is guaranteed to converge to a local maximum of the elbo under certain conditions on the step-size sequence.This is also called a learning rate or schedule in the machine learning community. Stochastic gradient ascent falls under the class of stochastic approximations, where Robbins and Monro (1951) established a pair of conditions that ensure convergence: most prominently, the step-size sequence must decay sufficiently quickly. Many sequences satisfy these criteria, but their specific forms impact the success of stochastic gradient ascent in practice. We describe an adaptive step-size sequence for advi below.

Adaptive step-size sequence. Adaptive step-size sequences retain (possibly infinite) memory about past gradients and adapt to the high-dimensional curvature of the elbo optimization space (Amari, 1998; Duchi et al., 2011; Ranganath et al., 2013; Kingma and Adam, 2015). These sequences enjoy theoretical bounds on their convergence rates. However, in practice, they can be slow to converge. The empirically justified rmsprop sequence (Tieleman and Hinton, 2012), which only retains finite memory of past gradients, converges quickly in practice but lacks any convergence guarantees. We propose a new step-size sequence which effectively combines both approaches.

Consider the step-size $\mathbf{\boldsymbol{\rho}}^{(i)}$ and a gradient vector $\mathbf{\boldsymbol{g}}^{(i)}$ at iteration $i$. We define the $k$th element of $\mathbf{\boldsymbol{\rho}}^{(i)}$ as

where we apply the following recursive update

with an initialization of $s^{(1)}_{k}={g^{2}_{k}}^{(1)}$.

The middle term $i^{-1/2+\epsilon}$ decays as a function of the iteration $i$. We set $\epsilon=10^{-16}$, a small value that guarantees that the step-size sequence satisfies the Robbins and Monro (1951) conditions.

The last term adapts to the curvature of the elbo optimization space. Memory about past gradients are processed in Equation (11). The weighting factor $\alpha\in(0,1)$ defines a compromise of old and new gradient information, which we set to $0.1$. The quantity $s_{k}$ converges to a non-zero constant. Without the previous decaying term, this would lead to possibly large oscillations around a local optimum of the elbo. The additional perturbation $\tau>0$ prevents division by zero and down-weights early iterations. In practice the step-size is not very sensitive to this value (Hoffman et al., 2013), so we set $\tau=1$.

Complexity and data subsampling. advi has complexity $\mathcal{O}(NMK)$ per iteration, where $N$ is the number of data points, $M$ is the number of mc samples (typically between 1 and 10), and $K$ is the number of latent variables. Classical vi which hand-derives a coordinate ascent algorithm has complexity $\mathcal{O}(NK)$ per pass over the dataset. The added complexity of automatic differentiation over analytic gradients is roughly constant (Carpenter et al., 2015; Baydin et al., 2015).

We scale advi to large datasets using stochastic optimization with data subsampling (Hoffman et al., 2013; Titsias and Lázaro-Gredilla, 2014). The adjustment to Algorithm 1 is simple: sample a minibatch of size $B\ll N$ from the dataset and scale the likelihood of the model by $N/B$ (Hoffman et al., 2013). The stochastic extension of advi has a per-iteration complexity $\mathcal{O}(BMK)$.

In Sections 4.3 and 4.4, we apply this stochastic extension to analyze datasets with millions of observations.

7 Related Work

advi automates variational inference within the Stan probabilistic programming system. This draws on two major themes.

The first theme is probabilistic programming. One class of systems focuses on probabilistic models where the user specifies a joint probability distribution. Some examples are BUGS (Spiegelhalter et al., 1995), JAGS (Plummer et al., 2003), and Stan (Stan Development Team, 2015). Another class of systems allows the user to directly specify more general probabilistic programs. Some examples are Church (Goodman et al., 2008), Figaro (Pfeffer, 2009), Venture (Mansinghka et al., 2014), and Anglican (Wood et al., 2014). Both classes primarily rely on various forms of mcmc techniques for inference; they typically cannot scale to very large data.

The second is a body of work that generalizes variational inference. Ranganath et al. (2014) and Salimans and Knowles (2014) propose a black-box technique that only requires computing gradients of the variational approximating family. Kingma and Welling (2014) and Rezende et al. (2014) describe a reparameterization of the variational problem that simplifies optimization. Titsias and Lázaro-Gredilla (2014) leverage the gradient of the model for a class of real-valued models. Rezende and Mohamed (2015) and Tran et al. (2016) improve the accuracy of black-box variational approximations. Here we build on and extend these ideas to automate variational inference; we highlight technical connections as we study the properties of advi in Section 3.

Some notable work crosses both themes. Bishop et al. (2002) present an automated variational algorithm for graphical models with conjugate exponential relationships between all parent-child pairs. Winn and Bishop (2005); Minka et al. (2014) extend this to graphical models with non-conjugate relationships by either using custom approximations or an expensive sampling approach. advi automatically supports a more comprehensive class of nonconjugate models; see Section 2.1. Wingate and Weber (2013) study a more general setting, where the variational approximation itself is a probabilistic program.

Properties of Automatic Differentiation Variational Inference

Automatic differentiation variational inference (advi) extends classical variational inference techniques in a few directions. In this section, we use simulated data to study three aspects of advi: the accuracy of mean-field and full-rank approximations, the variance of the advi gradient estimator, and the sensitivity to the transformation $T$.

We begin by considering three models that expose how the mean-field approximation affects the accuracy of advi.

We draw $1000$ datapoints from the model and run both variants of advi, mean-field and full-rank, until convergence. Figure 4 compares the advi methods to the exact posterior. Both procedures correctly identify the mean of the analytic posterior. However, the shape of the mean-field approximation is incorrect. This is because the mean-field approximation ignores off-diagonal terms of the Gaussian covariance. advi minimizes the kl divergence from the approximation to the exact posterior; this leads to a systemic underestimation of marginal variances (Bishop, 2006).

We simulated $9$ random covariates from the prior distribution (plus a constant intercept) and drew $1000$ datapoints from the likelihood. We estimated the posterior of the coefficients with advi and Stan’s default mcmc technique, the no-U-turn sampler (nuts) (Hoffman and Gelman, 2014). Figure 5 shows the marginal posterior densities obtained from each approximation. mcmc and advi perform similarly in their estimates of the posterior mean. The mean-field approximation, as expected, underestimates marginal posterior variances on most of the coefficients. The full-rank approximation, once again, better matches the posterior.

Stochastic volatility time-series model. Finally, we study a model where the data are not exchangeable. Consider an autoregressive process to model how the latent volatility (i.e., variance) of an economic asset changes over time (Kim et al., 1998); our goal is to estimate the sequence of volatilities. We expect these posterior estimates to be correlated, especially when the volatilities trend away from their mean value.

In detail, the price data exhibit latent volatility as part of the variance of a zero-mean Gaussian

where the log volatility follows an auto-regressive process

We place the following priors on the latent variables

We set $\mu=-1.025$, $\phi=0.9$ and $\sigma=0.6$, and simulate a dataset of $500$ time-steps from the generative model above. Figure 6 plots the posterior mean of the log volatility $h_{t}$ as a function of time. Mean-field advi struggles to describe the mean of the posterior, particularly when the log volatility drifts far away from $\mu$. In contrast, full-rank advi matches the estimates obtained from sampling.

We further investigate this by studying posterior correlations of the log volatility sequence. We draw $S=1000$ samples of $500$-dimensional log volatility sequences $\{\mathbf{\boldsymbol{h}}^{(s)}\}_{1}^{S}$. Figure 7 shows the empirical posterior covariance matrix, $\nicefrac{{1}}{{S-1}}\sum_{s}(\mathbf{\boldsymbol{h}}^{(s)}-\overline{\mathbf{\boldsymbol{h}}})(\mathbf{\boldsymbol{h}}^{(s)}-\overline{\mathbf{\boldsymbol{h}}})^{\top}$ for each method. The mean-field covariance (fig. 7(a)) fails to capture the locally correlated structure of the full-rank and sampling covariance matrices (figs. 7(b) and 7(c)). All covariance matrices exhibit a blurry spread due to finite sample size.

The regions where the local correlation is strongest correspond to the regions where mean-field underestimates the log volatility. To help identify these regions, we overlay the sampling mean log volatility estimate from Figure 6 above each matrix. Both full-rank advi and sampling results exhibit correlation where the log volatility trends away from its mean value.

Recommendations. How to choose between full-rank and mean-field advi? Scientists interested in posterior variances and covariances should use the full-rank approximation. Full-rank advi captures posterior correlations, in turn producing more accurate marginal variance estimates. For large data, however, full-rank advi can be prohibitively slow.

Scientists interested in prediction should initially rely on the mean-field approximation. Mean-field advi offers a fast algorithm for approximating the posterior mean. In practice, accurate posterior mean estimates dominate predictive accuracy; underestimating marginal variances matters less.

2 Variance of the Stochastic Gradients

advi uses Monte Carlo integration to approximate gradients of the elbo, and then uses these gradients in a stochastic optimization algorithm (Section 2). The speed of advi hinges on the variance of the gradient estimates. When a stochastic optimization algorithm suffers from high-variance gradients, it must repeatedly recover from poor parameter estimates.

advi is not the only way to compute Monte Carlo approximations of the gradient of the elbo. Black box variational inference (bbvi) takes a different approach (Ranganath et al., 2014). The bbvi gradient estimator uses the gradient of the variational approximation and avoids using the gradient of the model. For example, the following bbvi estimator

and the advi gradient estimator in Equation (7) both lead to unbiased estimates of the exact gradient. While bbvi is more general—it does not require the gradient of the model and thus applies to more settings—its gradients can suffer from high variance.

Figure 8 empirically compares the variance of both estimators for two models. Figure 8a shows the variance of both gradient estimators for a simple univariate model, where the posterior is a $\text{Gamma}(10,10)$. We estimate the variance using ten thousand re-calculations of the gradient $\nabla_{\mathbf{\boldsymbol{\phi}}}\mathcal{L}$, across an increasing number of mc samples $M$. The advi gradient has lower variance; in practice, a single sample suffices. (See the experiments in Section 4.)

3 Sensitivity to Transformations

advi uses a transformation $T$ from the unconstrained space to the constrained space. We now study how the choice of this transformation affects the non-Gaussian posterior approximation in the original latent variable space.

Figure 9 show the advi approximation under both transformations. Table 2 reports the corresponding kl divergences. Both graphical and numerical results prefer $T_{2}$ over $T_{1}$. A quick analysis corroborates this. $T_{1}$ is the logarithm, which flattens out for large values. However, $T_{2}$ is almost linear for large values of $\theta$. Since both the Gamma (the posterior) and the Gaussian (the advi approximation) densities are light-tailed, $T_{2}$ is the preferable transformation.

Is there an optimal transformation? Without loss of generality, we consider fixing a standard Gaussian distribution in the real coordinate space.For two transformations $T_{1}$ and $T_{2}$ from latent variable space to real coordinate space, there always exists a transformation $T_{3}$ within the real coordinate space such that $T_{1}(\theta)=T_{3}(T_{2}(\theta))$. The optimal transformation is then

This observation motivates pairing transformations with Gaussian variational approximations; there is no need for more complex variational families. advi takes the approach of using a library and a model compiler. This is not the only option. For example, Knowles (2015) posits a factorized Gamma density for positively constrained latent variables. In theory, this is equivalent to a mean-field Gaussian density paired with the transformation $T=P_{\text{Gamma}}$, the cumulative density function of the Gamma. (In practice, $P_{\text{Gamma}}$ is difficult to compute.) Challis and Barber (2012) study Fourier transform techniques for location-scale variational approximations beyond the Gaussian. Another option is to learn the transformation during optimization. We discuss recent approaches in this direction in Section 5.

Automatic Differentiation Variational Inference in Practice

We now apply automatic differentiation variational inference (advi) to an array of nonconjugate probability models. With simulated and real data, we study linear regression with automatic relevance determination, hierarchical logistic regression, several variants of non-negative matrix factorization, mixture models, and probabilistic principal component analysis. We compare mean-field advi to two mcmc sampling algorithms: Hamiltonian Monte Carlo (hmc) (Girolami and Calderhead, 2011) and nuts, which is an adaptive extension of hmc It is the default sampler in Stan. (Hoffman and Gelman, 2014).

To place advi and mcmc on a common scale, we report predictive likelihood on held-out data as a function of time. Specifically, we estimate the predictive likelihood

using Monte Carlo estimation. With mcmc, we run the chain and plug in each sample to estimate the integral above; with advi, we draw a sample from the variational approximation at every iteration.

We conclude with a case study: an exploratory analysis of millions of taxi rides. Here we show how a scientist might use advi in practice.

We begin with two nonconjugate regression models: linear regression with automatic relevance determination (ard) (Bishop, 2006) and hierarchical logistic regression (Gelman and Hill, 2006).

Linear regression with ard. This is a linear regression model with a hierarchical prior structure that leads to sparse estimates of the coefficients. (Details in Appendix F.1.) We simulate a dataset with $250$ regressors such that half of the regressors have no predictive power. We use $10\,000$ data points for training and withhold $1000$ for evaluation.

Logistic regression with a spatial hierarchical prior. This is a hierarchical logistic regression model from political science. The prior captures dependencies, such as states and regions, in a polling dataset from the United States 1988 presidential election (Gelman and Hill, 2006). The model is nonconjugate and would require some form of approximation to derive a classical vi algorithm. (Details in Appendix F.2.)

The dataset includes $145$ regressors, with age, education, and state and region indicators. We use $10\,000$ data points for training and withhold $1536$ for evaluation.

Results. Figure 10 plots average log predictive accuracy as a function of time. For these simple models, all methods reach the same predictive accuracy. We study advi with two settings of $M$, the number of mc samples used to estimate gradients. A single sample per iteration is sufficient; it is also the fastest. (We set $M=1$ from here on.)

2 Non-negative Matrix Factorization

We continue by exploring two nonconjugate non-negative matrix factorization models (Lee and Seung, 1999): a constrained Gamma Poisson model (Canny, 2004) and a Dirichlet Exponential Poisson model. Here, we show how easy it is to explore new models using advi. In both models, we use the Frey Face dataset, which contains $1956$ frames ($28\times 20$ pixels) of facial expressions extracted from a video sequence.

Constrained Gamma Poisson. This is a Gamma Poisson matrix factorization model with an ordering constraint: each row of one of the Gamma factors goes from small to large values. (Details in Appendix F.3.)

Dirichlet Exponential Poisson. This is a nonconjugate Dirichlet Exponential factorization model with a Poisson likelihood. (Details in Appendix F.4.)

Results. Figure 11 shows average log predictive accuracy as well as ten factors recovered from both models. advi provides an order of magnitude speed improvement over nuts. nuts struggles with the Dirichlet Exponential model. In both cases, hmc does not produce any useful samples within a budget of one hour; we omit hmc from here on.

The Gamma Poisson model appears to pick significant frames out of the dataset. The Dirichlet Exponential factors are sparse and indicate components of the face that move, such as eyebrows, cheeks, and the mouth.

3 Gaussian Mixture Model

This is a nonconjugate Gaussian mixture model (gmm) applied to color image histograms. We place a Dirichlet prior on the mixture proportions, a Gaussian prior on the component means, and a lognormal prior on the standard deviations. (Details in Appendix F.5.) We explore the imageclef dataset, which has $250\,000$ images (Villegas et al., 2013). We withhold $10\,000$ images for evaluation.

In Figure 12a we randomly select $1000$ images and train a model with $10$ mixture components. advi quickly finds a good solution. nuts struggles to find an adequate solution and hmc fails altogether (not shown). This is likely due to label switching, which can affect hmc-based algorithms in mixture models (Stan Development Team, 2015).

Figure 12b shows advi results on the full dataset. We increase the number of mixture components to $30$. Here we use advi, with additional stochastic subsampling of minibatches from the data (Hoffman et al., 2013). With a minibatch size of $500$ or larger, advi reaches high predictive accuracy. Smaller minibatch sizes lead to suboptimal solutions, an effect also observed in Hoffman et al. (2013). advi converges in about two hours; nuts cannot handle such large datasets.

4 A Case Study: Exploring Millions of Taxi Trajectories

How might a scientist use advi in practice? How easy is it to develop and revise new models? To answer these questions, we apply advi to a modern exploratory data analysis task: analyzing traffic patterns. In this section, we demonstrate how advi enables a scientist to quickly develop and revise complex hierarchical models.

The city of Porto has a centralized taxi system of 442 cars. When serving customers, each taxi reports its spatial location at 15 second intervals; this sequence of $(x,y)$ coordinates describes the trajectory and duration of each trip. A dataset of trajectories is publicly available: it contains all 1.7 million taxi rides taken during the year 2014 (European Conference of Machine Learning, 2015).

The trajectories have structure; for example, major roads and highways appear frequently. This motivates an approach where we first identify a lower-dimensional representation of the data to capture aggregate features, and then we cluster the trajectories in this representation. This is easier than clustering them in the original data space.

We begin with simple dimension reduction: probabilistic principal component analysis (ppca) (Bishop, 2006). This is a Bayesian generalization of classical principal component analysis, which is easy to write in Stan. However, like its classical counterpart, ppca does not identify how many principal components to use for the subspace. To address this, we propose an extension: ppca with automatic relevance determination (ard).

ppca with ard identifies the latent dimensions that are most effective at explaining variation in the data. The strategy is similar to that in Section 4.1. We assume that there are $100$ latent dimensions (i.e., the same dimension as the data) and impose a hierarchical prior that encourages sparsity. Consequently, the model only uses a subset of the latent dimensions to describe the data. (Details in Appendix F.6.)

We randomly subsample ten thousand trajectories and use advi to infer a subspace. Figure 13 plots the progression of the elbo. advi converges in approximately an hour and finds an eleven-dimensional subspace. We omit sampling results as both hmc and nuts struggle with the model; neither produce useful samples within an hour.

Equipped with this eleven-dimensional subspace, we turn to analyzing the full dataset of 1.7 million taxi trajectories. We first project all trajectories into the subspace. We then use the gmm from Section 4.3 ($K=30$) components to cluster the trajectories. advi takes less than half an hour to converge.

Figure 14 shows a visualization of fifty thousand randomly sampled trajectories. Each color represents the set of trajectories that associate with a particular Gaussian mixture. The clustering is geographical: taxi trajectories that are close to each other are bundled together. The clusters identify frequently taken taxi trajectories.

When we processed the raw data, we interpolated each trajectory to an equal length. This discards all duration information. What if some roads are particularly prone to traffic? Do these roads lead to longer trips?

Supervised probabilistic principal component analysis (sup-ppca) is one way to model this. The idea is to regress the durations of each trip onto a subspace that also explains variation in a response variable, in this case, the duration. sup-ppca is a simple extension of ppca (Murphy, 2012). We further extend it using the same ard prior as before. (Details in Appendix F.7.)

advi enables a quick repeat of the above analysis, this time with sup-ppca. With advi, we find another set of gmm clusters in less than two hours. These clusters, however, are more informative.

Figure 15 shows two clusters that identify particularly busy roads: the bridges of Porto that cross the Duoro river. Figure 15(a) shows a group of short trajectories that use the two old bridges near the city center. Figure 15(b) show a group of longer trajectories that use the two newer bridges that connect highways that circumscribe the city.

Analyzing these taxi trajectories illustrates how exploratory data analysis is an iterative effort: we want to rapidly evaluate models and modify them based on what we learn. advi, which provides automatic and fast inference, enables effective exploration of massive datasets.

Discussion

We presented automatic differentiation variational inference (advi), a variational inference tool that works for a large class of probabilistic models. The main idea is to transform the latent variables into a common space. Solving the variational inference problem in this common space solves it for all models in the class. We studied advi using ten different probability models; this showcases how easy it is to use advi in practice. We also developed and deployed advi as part of Stan, a probabilistic programming system; this makes advi available to everyone.

Begin with accuracy. As we showed in Section 3.3, advi can be sensitive to the transformations that map the constrained parameter space to the real coordinate space. Dinh et al. (2014) and Rezende and Mohamed (2015) use a cascade of simple transformations to improve accuracy. Tran et al. (2016) place a Gaussian process to learn the optimal transformation and prove its expressiveness as a universal approximator. A class of hierarchical variational models (Ranganath et al., 2015) extend these complex distributions to discrete latent variable models.

Continue with optimization. advi uses first-order automatic differentiation to implement stochastic gradient ascent. Higher-order gradients may enable faster convergence; however computing higher-order gradients comes at a computational cost (Fan et al., 2015). Optimization using line search could also improve convergence speed and robustness (Mahsereci and Hennig, 2015), as well as natural gradient approaches for nonconjugate models (Khan et al., 2015).

Follow with practical heuristics. Two things affect advi convergence: initialization and step-size scaling. We initialize advi in the real coordinate space as a standard Gaussian. A better heuristic could adapt to the model and dataset based on moment matching. We adaptively tune the scale of the step-size sequence using a finite search. A better heuristic could avoid this additional computation.

End with probabilistic programming. We designed and deployed advi with Stan in mind. Thus, we focused on the class of differentiable probability models. How can we extend advi to discrete latent variables? One approach would be to adapt advi to use the black box gradient estimator for these variables (Ranganath et al., 2014). This requires some care as these gradients will exhibit higher variance than the gradients with respect to the differentiable latent variables. (See Section 3.2.) With support for discrete latent variables, modified versions of advi could be extended to more general probabilistic programming systems, such as Church (Goodman et al., 2008), Figaro (Pfeffer, 2009), Venture (Mansinghka et al., 2014), and Anglican (Wood et al., 2014).

Acknowledgments. We thank Bruno Jacobs, and the reviewers for their helpful comments. This work is supported by NSF IIS-0745520, IIS-1247664, IIS-1009542, SES-1424962, ONR N00014-11-1-0651, DARPA FA8750-14-2-0009, N66001-15-C-4032, Sloan G-2015-13987, IES DE R305D140059, NDSEG, Facebook, Adobe, Amazon, and the Siebel Scholar and John Templeton Foundations.

Appendix A Transformations of Continuous Probability Densities

We present a brief summary of transformations, largely based on (Olive, 2014).

Consider a scalar (univariate) random variable $X$ with probability density function $f_{X}(x)$. Let $\mathcal{X}=\text{supp}(f_{X}(x))$ be the support of $X$. Now consider another random variable $Y$ defined as $Y=T(X)$. Let $\mathcal{Y}=\text{supp}(f_{Y}(y))$ be the support of $Y$.

If $T$ is a one-to-one and differentiable function from $\mathcal{X}$ to $\mathcal{Y}$, then $Y$ has probability density function

Let us sketch a proof. Consider the cumulative density function $Y$. If the transformation $T$ is increasing, we directly apply its inverse to the cdf of $Y$. If the transformation $T$ is decreasing, we apply its inverse to one minus the cdf of $Y$. The probability density function is the derivative of the cumulative density function. These things combined give the absolute value of the derivative above.

The extension to multivariate variables $\mathbf{\boldsymbol{X}}$ and $\mathbf{\boldsymbol{Y}}$ requires a multivariate version of the absolute value of the derivative of the inverse transformation. This is the absolute determinant of the Jacobian, $|\det J_{T^{-1}}(\mathbf{\boldsymbol{Y}})|$ where the Jacobian is

Intuitively, the Jacobian describes how a transformation warps unit volumes across spaces. This matters for transformations of random variables, since probability density functions must always integrate to one. If the transformation is linear, then we can drop the Jacobian adjustment; it evaluates to one. Similarly, affine transformations, like elliptical standardizations, also have Jacobians that evaluate to one; they preserve unit volumes.

Appendix B Transformation of the Evidence Lower Bound

Recall that $\mathbf{\boldsymbol{\zeta}}=T(\mathbf{\boldsymbol{\theta}})$ and that the variational approximation in the real coordinate space is $q(\mathbf{\boldsymbol{\zeta}}\,;\,\mathbf{\boldsymbol{\phi}})$.

We begin with the elbo in the original latent variable space. We then transform the latent variable space to the real coordinate space.

Appendix C Gradients of the Evidence Lower Bound

First, consider the gradient with respect to the $\mathbf{\boldsymbol{\mu}}$ parameter. We exchange the order of the gradient and the integration through the dominated convergence theorem (Çınlar, 2011). The rest is the chain rule for differentiation.

Then, consider the gradient with respect to the mean-field $\mathbf{\boldsymbol{\omega}}$ parameter.

Finally, consider the gradient with respect to the full-rank $\mathbf{\boldsymbol{L}}$ parameter.

Appendix D Automating Expectations: Monte Carlo Integration

Expectations are integrals. We can use mc integration to approximate them (Robert and Casella, 1999). All we need are samples from $q$.

mc integration provides noisy, yet unbiased, estimates of the integral. The standard deviation of the estimates are of order $1/\sqrt{S}$.

Appendix E Running advi in Stan

Visit http://mc-stan.org/ to download the latest version of Stan. Follow instructions on how to install Stan. You are then ready to use advi.

Stan offers multiple interfaces. We describe the command line interface (cmdStan) below,

where myData.data.R is the dataset stored in the R language Rdump format. output_advi.csv contains samples from the posterior and elbo_advi.csv reports the elbo.

Appendix F Details of Studied Models

Linear regression with ard is a high-dimensional sparse regression model (Bishop, 2006; Drugowitsch, 2013). We describe the model below. Stan code is in Figure 17.

The inputs are $\mathbf{\boldsymbol{x}}=x_{1:N}$ where each $x_{n}$ is $D$-dimensional. The outputs are $\mathbf{\boldsymbol{y}}=y_{1:N}$ where each $y_{n}$ is $1$-dimensional. The weights vector $\mathbf{\boldsymbol{w}}$ is $D$-dimensional. The likelihood

describes measurements corrupted by iid Gaussian noise with unknown standard deviation $\sigma$.

The ard prior and hyper-prior structure is as follows

where $\mathbf{\boldsymbol{\alpha}}$ is a $D$-dimensional hyper-prior on the weights, where each component gets its own independent Gamma prior.

We simulate data such that only half the regressions have predictive power. The results in Figure 10 use $a_{0}=b_{0}=c_{0}=d_{0}=1$ as hyper-parameters for the Gamma priors.

F.2 Hierarchical Logistic Regression

Hierarchical logistic regression models structured datasets in an intuitive way. We study a model of voting preferences from the 1988 United States presidential election. Chapter 14.1 of (Gelman and Hill, 2006) motivates the model and explains the dataset. We also describe the model below. Stan code is in Figure 18, based on (Stan Development Team, 2015).

The standard deviation terms all have uniform hyper-priors, constrained between 0 and 100.

F.3 Non-negative Matrix Factorization: Constrained Gamma Poisson Model

The Gamma Poisson factorization model describes discrete data matrices (Canny, 2004; Cemgil, 2009).

Consider a $U\times I$ matrix of observations. We find it helpful to think of $u=\{1,\cdots,U\}$ as users and $i=\{1,\cdots,I\}$ as items, as in a recommendation system setting. The generative process for a Gamma Poisson model with $K$ factors is

For each component $k$, draw $\theta_{uk}\sim\text{Gam}(a_{0},b_{0})$.

For each component $k$, draw $\beta_{ik}\sim\text{Gam}(c_{0},d_{0})$.

Draw the observation $y_{ui}\sim\text{Poisson}(\mathbf{\boldsymbol{\theta}}_{u}^{\top}\mathbf{\boldsymbol{\beta}}_{i})$.

A potential downfall of this model is that it is not uniquely identifiable: swapping rows and columns of $\mathbf{\boldsymbol{\theta}}$ and $\mathbf{\boldsymbol{\beta}}$ give the same inner product. One way to contend with this is to constrain either vector to be an ordered vector during inference. We constrain each $\mathbf{\boldsymbol{\theta}}_{u}$ vector in our model in this fashion. Stan code is in Figure 19. We set $K=10$ and all the Gamma hyper-parameters to 1 in our experiments.

F.4 Non-negative Matrix Factorization: Dirichlet Exponential Poisson Model

Another model for discrete data is a Dirichlet Exponential model. The Dirichlet enforces uniqueness while the exponential promotes sparsity. This is a non-conjugate model that does not appear to have been studied in the literature.

The generative process for a Dirichlet Exponential model with $K$ factors is

Draw the $K$-vector $\mathbf{\boldsymbol{\theta}}_{u}\sim\text{Dir}(\mathbf{\boldsymbol{\alpha}}_{0})$.

For each component $k$, draw $\beta_{ik}\sim\text{Exponential}(\lambda_{0})$.

Draw the observation $y_{ui}\sim\text{Poisson}(\mathbf{\boldsymbol{\theta}}_{u}^{\top}\mathbf{\boldsymbol{\beta}}_{i})$.

Stan code is in Figure 20. We set $K=10$, $\alpha_{0}=1000$ for each component, and $\lambda_{0}=0.1$. With this configuration of hyper-parameters, the factors $\mathbf{\boldsymbol{\beta}}_{i}$ appear sparse.

F.5 Gaussian Mixture Model

The Gaussian mixture model (gmm) is a celebrated probability model (Bishop, 2006). We use it to group a dataset of natural images based on their color histograms. We build a high-dimensional gmm with a Gaussian prior for the mixture means, a lognormal prior for the mixture standard deviations, and a Dirichlet prior for the mixture components.

Represent the images as $\mathbf{\boldsymbol{y}}=y_{1:N}$ where each $y_{n}$ is $D$-dimensional and there are $N$ observations. The likelihood for the images is

with a Dirichlet prior for the mixture proportions

and a lognormal prior for the mixture standard deviations

The dimension of the color histograms in the imageclef dataset is $D=576$. This is a concatenation of three $192$-length histograms, one for each color channel (red, green, blue) of the images.

We scale the image histograms to have zero mean and unit variance. Setting $\alpha_{0}$ to a small value encourages the model to use fewer components to explain the data. Larger values of $\alpha_{0}$ encourage the model to use all $K$ components. We set $\alpha_{0}=1\,000$ in our experiments.

advi code is in Figure 21. The stochastic data subsampling version of the code is in Figure 22.

F.6 Probabilistic Principal Component Analysis with Automatic Relevance Determination

Probabilistic principal component analysis (ppca) is a Bayesian extension of classical principal component analysis (Bishop, 2006). The generative process is straightforward. Consider a dataset of $\mathbf{\boldsymbol{x}}=x_{1:N}$ where each $x_{n}$ is $D$-dimensional. Let $M$ be the dimension of the subspace we seek.

First define a set of latent variables $\mathbf{\boldsymbol{z}}=z_{1:N}$ where each $z_{n}$ is $M$-dimensional. Draw each $z_{n}$ from a standard normal

Then define a set of principal components $\mathbf{\boldsymbol{w}}=w_{1:D}$ where each $w_{d}$ is $M$-dimensional. Similarly, draw the principal components from a standard normal

Finally define the likelihood through an inner product as

The standard deviation $\sigma$ is also a latent variable. Place a lognormal prior on it as

We extend ppca by adding an ard hierarchical prior. The extended model introduces a $M$-dimensional vector $\mathbf{\boldsymbol{\alpha}}$ which chooses which principal components to retain. ($M<D$ now represents the maximum number of principal components to consider.) The extended extends the above by

F.7 Supervised Probabilistic Principal Component Analysis with Automatic Relevance Determination

Supervised probabilistic principal component analysis (sup-ppca) augments ppca by regressing a vector of observed random variables $y$ onto the principal component subspace. The idea is to not only find a set of principal components that describe variation in the dataset $\mathbf{\boldsymbol{x}}$, but to also predict $y$. The complete model is

References