Variational Dropout and\\ the Local Reparameterization Trick
Diederik P. Kingma$^*$
Tim Salimans$^\times$ and Max Welling$^{*\dagger
Diederik P. Kingma$^*$
Tim Salimans$^\times$ and Max Welling$^{*\dagger
We investigate a local reparameterizaton technique for greatly reducing the variance of stochastic gradients for variational Bayesian inference (SGVB) of a posterior over model parameters, while retaining parallelizability. This local reparameterization translates uncertainty about global parameters into local noise that is independent across datapoints in the minibatch. Such parameterizations can be trivially parallelized and have variance that is inversely proportional to the minibatch size, generally leading to much faster convergence. Additionally, we explore a connection with dropout: Gaussian dropout objectives correspond to SGVB with local reparameterization, a scale-invariant prior and proportionally fixed posterior variance. Our method allows inference of more flexibly parameterized posteriors; specifically, we propose variational dropout, a generalization of Gaussian dropout where the dropout rates are learned, often leading to better models. The method is demonstrated through several experiments.
fsgvb.tex
\documentclass{article} % For LaTeX2e
\usepackage{nips15submit_e,times}
\usepackage{hyperref}
\usepackage{url}

\title{Variational Dropout and\\ the Local Reparameterization Trick}

\author{Diederik P. Kingma$^*$, Tim Salimans$^\times$ and Max Welling$^{*\dagger}$ \\
$^*$ Machine Learning Group, University of Amsterdam\\
$^\times$ Algoritmica\\
$^\dagger$ University of California, Irvine, and the Canadian Institute for Advanced Research (CIFAR)\\
\texttt{D.P.Kingma@uva.nl, salimans.tim@gmail.com, M.Welling@uva.nl}
}

\newcommand{\fix}{\marginpar{FIX}}
\newcommand{\new}{\marginpar{NEW}}

\usepackage[sort&compress, square, numbers]{natbib}
\usepackage{amsmath, amsthm}
\usepackage{amsfonts}
\usepackage{amssymb}

\usepackage[pdftex]{graphicx} % more modern
%\usepackage{subfigure} 
\usepackage{subcaption}
\graphicspath{{figures/}} 

%\usepackage{algpseudocode}
\usepackage{algorithm} 
\usepackage{algorithmic}
\usepackage{hyperref}
\usepackage{balance}
\usepackage{url}
%\newcommand{\theHalgorithm}{\arabic{algorithm}}

\usepackage{float}
\usepackage{placeins}


% Package for highlighting text
\usepackage{soul}

\newcommand{\ba}{\mathbf{a}}
\newcommand{\bb}{\mathbf{b}}
\newcommand{\bd}{\mathbf{d}}
\newcommand{\bw}{\mathbf{w}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\by}{\mathbf{y}}
\newcommand{\Exp}[2]{\mathbb{E}_{#1}\left[#2\right]}
\newcommand{\Var}[2]{\text{Var}_{#1}\left[#2\right]}
\newcommand{\Cov}[2]{\text{Cov}_{#1}\left[#2\right]}
\newcommand{\Corr}[2]{\text{Corr}_{#1}\left[#2\right]}

\newcommand{\bA}{\mathbf{A}}
\newcommand{\bB}{\mathbf{B}}
\newcommand{\bV}{\mathbf{V}}
\newcommand{\bW}{\mathbf{W}}
\newcommand{\bv}{\mathbf{v}}

\newcommand{\bphi}{\mathbf{\phi}}
\newcommand{\bt}{\mathbf{\theta}}
\newcommand{\bdelta}{\mathbf{\delta}}
\newcommand{\beps}{\mathbf{\epsilon}}
\newcommand{\q}{q_\mathbf{\phi}}
\newcommand{\bmu}{\mathbf{\mu}}
\newcommand{\bsigma}{\mathbf{\sigma}}
\newcommand{\bSigma}{\mathbf{\Sigma}}

\newcommand{\D}{\mathcal{D}}

\newcommand{\eqnr}{\addtocounter{equation}{1}\tag{\theequation}}

\nipsfinalcopy % Uncomment for camera-ready version

\begin{document}


\maketitle

\begin{abstract}
We investigate a local reparameterizaton technique for greatly reducing the variance of stochastic gradients for variational Bayesian inference (SGVB) of a posterior over model parameters, while retaining parallelizability. This local reparameterization translates uncertainty about global parameters into local noise that is independent across datapoints in the minibatch. Such parameterizations can be trivially parallelized and have variance that is inversely proportional to the minibatch size, generally leading to much faster convergence.
Additionally, we explore a connection with dropout: Gaussian dropout objectives correspond to SGVB with local reparameterization, a scale-invariant prior and proportionally fixed posterior variance. Our method allows inference of more flexibly parameterized posteriors; specifically, we propose \emph{variational dropout}, a generalization of Gaussian dropout where the dropout rates are learned, often leading to better models. The method is demonstrated through several experiments.
\end{abstract}

\section{Introduction}
Deep neural networks are a flexible family of models that easily scale to millions of parameters and datapoints, but are still tractable to optimize using minibatch-based stochastic gradient ascent. Due to their high flexibility, neural networks have the capacity to fit a wide diversity of nonlinear patterns in the data. This flexbility often leads to \emph{overfitting} when left unchecked: spurious patterns are found that happen to fit well to the training data, but are not predictive for new data. Various regularization techniques for controlling this overfitting are used in practice; a currently popular and empirically effective technique being \emph{dropout}~\cite{hinton2012improving}. In \cite{wang2013fast} it was shown that regular (binary) dropout has a Gaussian approximation called \emph{Gaussian dropout} with virtually identical regularization performance but much faster convergence. In section 5 of \cite{wang2013fast} it is shown that Gaussian dropout optimizes a lower bound on the marginal likelihood of the data. In this paper we show that a relationship between dropout and Bayesian inference can be extended and exploited to greatly improve the efficiency of variational Bayesian inference on the model parameters. This work has a direct interpretation as a generalization of Gaussian dropout, with the same fast convergence but now with the freedom to specify more flexibly parameterized posterior distributions.

Bayesian posterior inference over the neural network parameters is a theoretically attractive method for controlling overfitting; exact inference is computationally intractable, but efficient approximate schemes can be designed. Markov Chain Monte Carlo (MCMC) is a class of approximate inference methods with asymptotic guarantees, pioneered by~\cite{neal1995bayesian} for the application of regularizing neural networks. Later useful refinements include ~\cite{welling2011bayesian} and \cite{ahn2012bayesian}.

An alternative to MCMC is variational inference~\cite{hinton1993keeping} or the equivalent \emph{minimum description length} (MDL) framework. Modern variants of stochastic variational inference have been applied to neural networks with some succes~\cite{graves2011practical}, but have been limited by high variance in the gradients. Despite their theoretical attractiveness, Bayesian methods for inferring a posterior distribution over neural network weights have not yet been shown to outperform simpler methods such as dropout. Even a new crop of efficient variational inference algorithms based on stochastic gradients with minibatches of data~\cite{salimans2013fixedform,kingma2013auto,rezende2014stochastic} have not yet been shown to significantly improve upon simpler dropout-based regularization.

In section~\ref{sec:variational_dropout} we explore an as yet unexploited trick for improving the efficiency of stochastic gradient-based variational inference with minibatches of data, by translating uncertainty about global parameters into local noise that is independent across datapoints in the minibatch. The resulting method has an optimization speed on the same level as fast dropout~\cite{wang2013fast}, and indeed has the original Gaussian dropout method as a special case. An advantage of our method is that it allows for full Bayesian analysis of the model, and that it's significantly more flexible than standard dropout. The approach presented here is closely related to several popular methods in the literature that regularize by adding random noise; these relationships are discussed in section~\ref{sec:relatedwork}. 

\section{Efficient and Practical Bayesian Inference}
\label{sec:variational_dropout}
We consider Bayesian analysis of a dataset $\D$, containing a set of $N$ i.i.d. observations of tuples $(\bx,\by)$, where the goal is to learn a model with \emph{parameters} or \emph{weights} $\bw$ of the conditional probability $p(\by|\bx,\bw)$ (standard classification or regression)\footnote{Note that the described method is not limited to classification or regression and is straightforward to apply to other modeling settings like unsupervised models and temporal models.}. Bayesian inference in such a model consists of updating some initial belief over parameters $\bw$ in the form of a prior distribution $p(\bw)$, after observing data $\D$, into an updated belief over these parameters in the form of (an approximation to) the posterior distribution $p(\bw|\D)$. Computing the true posterior distribution through Bayes' rule $p(\bw|\D) = p(\bw)p(\D|\bw)/p(\D)$ involves computationally intractable integrals, so good approximations are necessary. In \emph{variational inference}, inference is cast as an optimization problem where we optimize the parameters $\bphi$ of some parameterized model $\q(\bw)$ such that $\q(\bw)$ is a close approximation to $p(\bw|\D)$ as measured by the Kullback-Leibler divergence $D_{KL}(\q(\bw)||p(\bw|\D))$. This divergence of our posterior $\q(\bw)$ to the true posterior is minimized in practice by maximizing the so-called variational lower bound $\mathcal{L}(\bphi)$ of the marginal likelihood of the data:
\begin{align}
\mathcal{L}(\bphi) &= - D_{KL}(\q(\bw)||p(\bw)) + L_\D(\bphi)\label{eq:bound}\\
\text{where\quad} L_\D(\bphi) &=  \sum_{(\bx, \by) \in \D} \Exp{\q(\bw)}{ \log p(\by|\bx,\bw)}
\label{eq:bound2}
\end{align}
We'll call $L_\D(\bphi)$ the \emph{expected log-likelihood}. The bound $\mathcal{L}(\bphi)$ plus $D_{KL}(\q(\bw)||p(\bw|\D))$ equals the (conditional) marginal log-likelihood $\sum_{(\bx, \by) \in \D} \log p(\by|\bx)$. Since this marginal log-likelihood is constant w.r.t. $\bphi$, maximizing the bound w.r.t. $\bphi$ will minimize $D_{KL}(\q(\bw)||p(\bw|\D))$. %, which is our goal.

\subsection{Stochastic Gradient Variational Bayes (SGVB)}\label{sec:sgvb}
Various algorithms for gradient-based optimization of the variational bound (eq.~\eqref{eq:bound}) with differentiable $q$ and $p$ exist. See section~\ref{sec:relatedwork} for an overview. A recently proposed efficient method for minibatch-based optimization with differentiable models is the \emph{stochastic gradient variational Bayes} (SGVB) method introduced in~\cite{kingma2013auto} (especially appendix F) and ~\cite{rezende2014stochastic}. The basic trick in SGVB is to parameterize the random parameters $\bw \sim \q(\bw)$ as: $\bw = f(\beps, \bphi)$ where $f(.)$ is a differentiable function and $\beps \sim p(\beps)$ is a random noise variable. In this new parameterisation, an unbiased differentiable minibatch-based Monte Carlo estimator of the expected log-likelihood can be formed:
\begin{align}
L_\D(\bphi) \simeq L_\D^\text{SGVB}(\bphi) = \frac{N}{M} \sum_{i=1}^M \log p(\by^i|\bx^i,\bw = f(\beps, \bphi)),
\label{eq:objective_sgvb}
\end{align}
where $(\bx^i,\by^i)_{i=1}^M$ is a minibatch of data with $M$ random datapoints $(\bx^i,\by^i) \sim \D$, and $\beps$ is a noise vector drawn from the noise distribution $p(\beps)$. We'll assume that the remaining term in the variational lower bound, $D_{KL}(\q(\bw)||p(\bw))$, can be computed deterministically, but otherwise it may be approximated similarly. The estimator \eqref{eq:objective_sgvb} is differentiable w.r.t. $\bphi$ and unbiased, so its gradient is also unbiased: $\nabla_\phi L_\D(\bphi) \simeq \nabla_\bphi L_\D^\text{SGVB}(\bphi)$. We can proceed with variational Bayesian inference by randomly initializing $\bphi$ and performing stochastic gradient ascent on $\mathcal{L}(\bphi)$~\eqref{eq:bound}.

\subsection{Variance of the SGVB estimator}\label{sec:variance}
The theory of stochastic approximation tells us that stochastic gradient ascent using \eqref{eq:objective_sgvb} will asymptotically converge to a local optimum for an appropriately declining step size and sufficient weight updates \cite{robbinsmonro}, but in practice the performance of stochastic gradient ascent crucially depends on the variance of the gradients. If this variance is too large, stochastic gradient descent will fail to make much progress in any reasonable amount of time. Our objective function consists of an expected log likelihood term that we approximate using Monte Carlo, and a KL divergence term $D_{KL}(\q(\bw)||p(\bw))$ that we assume can be calculated analytically and otherwise be approximated with Monte Carlo with similar reparameterization. 

Assume that we draw minibatches of datapoints with replacement; see appendix~\ref{sec:variance_without_replacement} for a similar analysis for minibatches without replacement.
Using $L_i$ as shorthand for $\log p(\by^i|\bx^i,\bw = f(\beps^i, \bphi))$, the contribution to the likelihood for the $i$-th datapoint in the minibatch, the Monte Carlo estimator \eqref{eq:objective_sgvb} may be rewritten as
$L_\D^\text{SGVB}(\bphi) = \frac{N}{M} \sum_{i=1}^M L_{i}$, whose variance is given by
\begin{align}
\Var{}{L_\D^\text{SGVB}(\bphi)}
=& \frac{N^2}{M^2} \Big( \sum_{i=1}^M \Var{}{L_{i}} + 2 \sum_{i=1}^M \sum_{j=i+1}^M \Cov{}{L_{i},L_{j}} \Big)\\
=& N^2 \Big( \frac{1}{M} \Var{}{L_{i}} + \frac{M-1}{M} \Cov{}{L_{i},L_{j}} \Big),
\label{eq:objective_sgvb_var}
\end{align}
where the variances and covariances are w.r.t.\ both the data distribution and $\beps$ distribution, i.e. $\Var{}{L_{i}} = \Var{\beps,\bx^i,\by^i}{\log p(\by^i|\bx^i,\bw = f(\beps, \bphi))}$, with $\bx^i,\by^i$ drawn from the empirical distribution defined by the training set. As can be seen from \eqref{eq:objective_sgvb_var}, the total contribution to the variance by $\Var{}{L_i}$ is inversely proportional to the minibatch size $M$. However, the total contribution by the covariances does not decrease with $M$. In practice, this means that the variance of $L_\D^\text{SGVB}(\bphi)$ can be dominated by the covariances for even moderately large $M$.

\subsection{Local Reparameterization Trick}\label{sec:fastsgvb}

We therefore propose an alternative estimator for which we have $\Cov{}{L_{i},L_{j}} = 0$, so that the variance of our stochastic gradients scales as $1/M$. We then make this new estimator computationally efficient by not sampling $\beps$ directly, but only sampling the intermediate variables $f(\beps)$ through which $\beps$ influences $L_\D^\text{SGVB}(\bphi)$. By doing so, the global uncertainty in the weights is translated into a form of local uncertainty that is independent across examples and easier to sample. We refer to such a reparameterization from global noise to local noise as the \emph{local reparameterization trick}. Whenever a source of global noise can be translated to local noise in the intermediate states of computation ($\beps \rightarrow f(\beps)$), a local reparameterization can be applied to yield a computationally and statistically efficient gradient estimator.

Such local reparameterization applies to a fairly large family of models, but is best explained through a simple example: Consider a standard fully connected neural network containing a hidden layer consisting of 1000 neurons. This layer receives an $M \times 1000$ input feature matrix $\bA$ from the layer below, which is multiplied by a $1000 \times 1000$ weight matrix $\bW$, before a nonlinearity is applied, i.e. $\bB = \bA\bW$. We then specify the posterior approximation on the weights to be a fully factorized Gaussian, i.e.\ $\q(w_{i,j}) = N(\mu_{i,j},\sigma^{2}_{i,j}) \hspace{0.1cm} \forall w_{i,j} \in \bW$, which means the weights are sampled as $w_{i,j} = \mu_{i,j} + \sigma_{i,j}\epsilon_{i,j}$, with $\epsilon_{i,j} \sim N(0,1)$. In this case we could make sure that $\Cov{}{L_{i},L_{j}} = 0$ by sampling a separate weight matrix $\bW$ for each example in the minibatch, but this is not computationally efficient: we would need to sample $M$ million random numbers for just a single layer of the neural network. Even if this could be done efficiently, the computation following this step would become much harder: Where we originally performed a simple matrix-matrix product of the form $\bB = \bA\bW$, this now turns into $M$ separate local vector-matrix products. The theoretical complexity of this computation is higher, but, more importantly, such a computation can usually not be performed in parallel using fast device-optimized BLAS (Basic Linear Algebra Subprograms). This also happens with other neural network architectures such as convolutional neural networks, where optimized libraries for convolution cannot deal with separate filter matrices per example.

Fortunately, the weights (and therefore $\beps$) only influence the expected log likelihood through the neuron activations $\bB$, which are of much lower dimension. If we can therefore sample the random activations $\bB$ directly, without sampling $\bW$ or $\beps$, we may obtain an efficient Monte Carlo estimator at a much lower cost. For a factorized Gaussian posterior on the weights, the posterior for the activations (conditional on the input $\bA$) is also factorized Gaussian:
\begin{align}
\q(w_{i,j}) &= N(\mu_{i,j},\sigma^{2}_{i,j}) \hspace{0.1cm} \forall w_{i,j} \in \bW \hspace{0.2cm} \Longrightarrow \hspace{0.2cm} \q(b_{m,j}|\bA) = N(\gamma_{m,j},\delta_{m,j}), \text{ with } \nonumber\\
\gamma_{m,j} &= \sum_{i=1}^{1000} a_{m,i}\mu_{i,j}, \hspace{0.2cm} \text{ and } \hspace{0.2cm} \delta_{m,j} = \sum_{i=1}^{1000} a^{2}_{m,i}\sigma^{2}_{i,j}.
\end{align}
Rather than sampling the Gaussian weights and then computing the resulting activations, we may thus sample the activations from their implied Gaussian distribution directly, using $b_{m,j} = \gamma_{m,j} + \sqrt{\delta_{m,j}}\zeta_{m,j}$, with $\zeta_{m,j} \sim N(0,1)$. Here, $\zeta$ is an $M \times 1000$ matrix, so we only need to sample $M$ thousand random variables instead of $M$ million: a thousand fold savings.

In addition to yielding a gradient estimator that is more \textit{computationally efficient} than drawing separate weight matrices for each training example, the local reparameterization trick also leads to an estimator that has \textit{lower variance}. To see why, consider the stochastic gradient estimate with respect to the posterior parameter $\sigma^{2}_{i,j}$ for a minibatch of size $M=1$. Drawing random weights $\bW$, we get
\begin{align}
\label{eq:stochest_weights}
\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}} = \frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\frac{\epsilon_{i,j}a_{m,i}}{2\sigma_{i,j}}.
\end{align}
If, on the other hand, we form the same gradient using the local reparameterization trick, we get
\begin{align}
\label{eq:stochest_local}
\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}} = \frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\frac{\zeta_{m,j}a_{m,i}^{2}}{2\sqrt{\delta_{m,j}}}.
\end{align}
Here, there are two stochastic terms: The first is the backpropagated gradient $\partial L_{\D}^\text{SGVB} / \partial b_{m,j}$, and the second is the sampled random noise ($\epsilon_{i,j}$ or $\zeta_{m,j}$). Estimating the gradient with respect to $\sigma^{2}_{i,j}$ then basically comes down to estimating the covariance between these two terms. This is much easier to do for $\zeta_{m,j}$ as there are much fewer of these: individually they have higher correlation with the backpropagated gradient $\partial L_\D^\text{SGVB}/\partial b_{m,j}$, so the covariance is easier to estimate. In other words, measuring the effect of $\zeta_{m,j}$ on $\partial L_\D^\text{SGVB}/\partial b_{m,j}$ is easy as $\zeta_{m,j}$ is the only random variable directly influencing this gradient via $b_{m,j}$. On the other hand, when sampling random weights, there are a thousand $\epsilon_{i,j}$ influencing each gradient term, so their individual effects get lost in the noise. In appendix \ref{sec:low_variance_reparam} we make this argument more rigorous, and in section~\ref{sec:experiments} we show that it holds experimentally.


%The problem of approximating the gradient thus boils down to estimating one of the correlation coefficients $\rho_{1} = \Corr{q,\D}{\partial L_\D^\text{SGVB}/\partial b_{m,j}, \epsilon_{i,j}}$ or $\rho_{2} = \Corr{q,\D}{\partial L_\D^\text{SGVB}/\partial b_{m,j}, \zeta_{m,j}}$. As is well known, accurately estimating a correlation coefficient is easier if the absolute correlation is large: in general the standard error for the sample correlation coefficient is given by
%\[
%\text{SE}{\rho} = \frac{1 - \rho^{2}}{\sqrt{M - 1}}.
%\]
%Since there are a thousand times fewer $\zeta_{m,j}$ than $\epsilon_{i,j}$, individually they have much higher correlation with the gradient $\partial L_\D^\text{SGVB}/\partial b_{m,j}}$, so they are much easier to estimate. In other words, measuring the effect of $\zeta_{m,j}$ on $\partial L_\D^\text{SGVB}/\partial b_{m,j}}$ is easy as $\zeta_{m,j}$ is the only random variable influencing $b_{m,j}$. On the other hand, when sampling random weights, there are a thousand $\epsilon_{i,j}$ influencing each $b_{m,j}$, so their individual effects get lost in the noise.



%Trivial parallelizability can be retained by translating weight uncertainty in $\q(\bw)$ into local uncertainty in the intermediate states of computation. We call this the \emph{local reparameterisation trick}. Let $\bb^i = f(\ba^i, \bw^i)$ be an intermediate variable of computation for datapoint $i$, with $\ba^i$ a previous state of computation (typically a function of the model input $\bx^i$) and $\bw^i$ the model parameters with distribution $\bw^i \sim \q(\bw)$ that is equal for all $i$. Let $L(\bb^i)$ be (part of) an objective, such as a log-likelihood function, whose expectation w.r.t. $\q(\bw)$ we would like to estimate (equations ~\eqref{eq:bound} and \eqref{eq:objective_sgvb}). Given the distribution $\q(\bw)$, we equivalently have some distribution over $\bb^i$, namely $\q(\bb^i|\ba^i)$. The trick is to parameterize $\bb^i \sim \q(\bb^i|\ba^i)$ as a differentiable function $\bb^i = g^\text{VD}(\beps^i, \ba^i, \bphi)$ with random noise $\beps^i \sim p^\text{VD}(\beps)$, such that:
%\begin{align}
%\Exp{\q(\bw)}{L(f(\ba^i, \bw^i))}
%&= \Exp{\q(\bb^i|\ba^i)}{L(\bb^i)}
%= \Exp{p^\text{VD}(\beps)}{L(g^\text{VD}(\beps^i,\ba^i,\bphi))}
%\end{align}
%The distribution of the original SGVB parameterisation $\bb^i = f(\ba^i, f(\beps^i, \bphi))$ with $\beps^i \sim p(\beps)$, equals the distribution of the new parameterisation $\bb^i = g^\text{VD}(\beps^i, \ba^i, \bphi)$ with $\beps^i \sim p^\text{VD}(\beps)$. The key difference is that computation of $g^\text{VD}(\beps^i, \ba^i, \bphi)$ can often be trivially parallelized, whereas $f(\ba^i, f(\beps^i, \bphi))$ cannot. The resulting operations can be combined into a \emph{variational dropout} (VD) estimator
%\begin{align}
%L_\D(\bphi) \simeq L_\D^\text{VD}(\bphi)
%\label{eq:objective_fsgvb}
%\end{align}
%that is generally much more efficient than $L_\D^\text{SGVB}(\bphi)$ with minibatches.


%\subsection{Fast SGVB through Reparameterization with Local Noise[todo: improve clarity]}\label{sec:fastsgvb}
%
%There is a \emph{de facto} standard strategy for efficiently optimizing objective functions with minibatches of data: vectorize the steps of operations such that the bulk of computation is performed in a small number of large operations, such as large matrix-matrix operations or convolutional operations. Using device-optimized BLAS (Basic Linear Algebra Subprograms) libraries, these few large operations can be automatically and efficiently parallelized across multiple CPU or GPU cores. This strategy is at the heart of efficient implementations in software like Matlab, Theano~\cite{bastien2012theano} and Torch~\cite{collobert2011torch7}, and typically leads to orders of magnitude faster computation than non-vectorized implementations.
%
%A property of $L_\D^\text{SGVB}$ is that the stochasticity is global: weights are sampled from $\q$ only once, then used across all datapoints in the minibatch. This is sub-optimal in terms of variance. A better estimator is one where the noise is independent per datapoint: the variance of the estimator is then inversely proportional to the minibatch size $M$. A first na\"ive attempt is:
%$L_\D^{\text{Na\"ive}}(\bphi) = \frac{N}{M} \sum_{i=1}^M \log p(\by^i|\bx^i,\bw^i)$
%with an explicit separate sample $\bw^i \sim \q(\bw)$ for each individual datapoint $i$ in the minibatch, using a parameterization $\bw^i = f(\beps^i, \bphi)$.
%however, we then lose the ability to easily parallelize this computation, generally leading to drastically slower computation as explained in the previous paragraph. 
%
%[Todo: equation that shows variance reduction effect of separate $\bw^i$ per element in the minibatch.]
%
%Trivial parallelizability can be retained by translating weight uncertainty in $\q(\bw)$ into local uncertainty in the intermediate states of computation. We call this the \emph{local reparameterisation trick}. Let $\bb^i = f(\ba^i, \bw^i)$ be an intermediate variable of computation for datapoint $i$, with $\ba^i$ a previous state of computation (typically a function of the model input $\bx^i$) and $\bw^i$ the model parameters with distribution $\bw^i \sim \q(\bw)$ that is equal for all $i$. Let $L(\bb^i)$ be (part of) an objective, such as a log-likelihood function, whose expectation w.r.t. $\q(\bw)$ we would like to estimate (equations ~\eqref{eq:bound} and \eqref{eq:objective_sgvb}). Given the distribution $\q(\bw)$, we equivalently have some distribution over $\bb^i$, namely $\q(\bb^i|\ba^i)$. The trick is to parameterize $\bb^i \sim \q(\bb^i|\ba^i)$ as a differentiable function $\bb^i = g^\text{local}(\beps^i, \ba^i, \bphi)$ with random noise $\beps^i \sim p^\text{local}(\beps)$, such that:
%\begin{align}
%\Exp{\q(\bw)}{L(f(\ba^i, \bw^i))}
%&= \Exp{\q(\bb^i|\ba^i)}{L(\bb^i)}
%= \Exp{p^\text{local}(\beps)}{L(g^\text{local}(\beps^i,\ba^i,\bphi))}
%\end{align}
%The distribution of the original SGVB parameterisation $\bb^i = f(\ba^i, f(\beps^i, \bphi))$ with $\beps^i \sim p(\beps)$, equals the distribution of the new parameterisation $\bb^i = g^\text{local}(\beps^i, \ba^i, \bphi)$ with $\beps^i \sim p^\text{local}(\beps)$. The key difference is that computation of $g^\text{local}(\beps^i, \ba^i, \bphi)$ can often be trivially parallelized, whereas $f(\ba^i, f(\beps^i, \bphi))$ cannot. The resulting operations can be combined into a \emph{fast SGVB} estimator
%\begin{align}
%L_\D(\bphi) \simeq L_\D^\text{F-SGVB}(\bphi)
%\label{eq:objective_fsgvb}
%\end{align}
%that is generally much more efficient than $L_\D^\text{SGVB}(\bphi)$ with minibatches.

%\subsection{Example: Dot Product, Gaussian Weight Uncertainty}
%Let us give an example of the local reparameterization trick, with the common linear operation $\bb^i = \ba^i \bW$, where the row vector $\ba^i$ is of dimension $1 \times N_A$ and is the result of earlier operations, and $\bb^i$ is of dimension $1 \times N_B$. Let $\bW$ be a parameter matrix, with approximate posterior $\q(\bW) = \prod_{jk}\q(W_{jk})$, with $\q(W_{jk}) = \mathcal{N}(\mu_{jk},\bSigma_{jk})$ where $\bmu$ and $\bSigma$ are variational parameter matrices. With $(\odot, (.)^2, \sqrt{(.)})$ we denote element-wise product, square and square-root.
%
%A non-parallelizable parameterisation of $\bb^i$, is: $\bb^i = \ba^i (\bmu+\beps^i \odot \sqrt{\bSigma})$.
%
%A trivially parallelizable reparameterisation of $\bb^i$ with local noise $\epsilon^i \sim \mathcal{N}(0,I)$ is: $\bb^i = \ba^i \bmu + \epsilon^i \odot \sqrt{(\ba^i)^2 \bSigma}$. Let $\bB = \bA \bW$, where $\bA$ is a matrix of size $M \times N_A$ with $\ba^i$ as $i$-th row, and $\bB$ is the corresponding result with size $M \times N_B$. The vectorized form is $\bB = \bA \bmu + \epsilon \odot \sqrt{\bA^2 \bSigma}$. In this form, most of the computation is contained in a few large matrix-matrix operations, whose computation is trivially parallelized.

%\section{Gaussian Dropout as a Special Case}
%\label{sec:gauss_dropout}
%In \cite{wang2013fast} it was shown that regularization by regular (binary) dropout \cite{srivastava2014dropout} has a fast approximation called \emph{Gaussian dropout}. In Gaussian dropout, linear operations $\bb^i = \ba^i \bW$ are replaced by a stochastic relation $p(\bb^i|\ba^i) = \mathcal{N}(\ba^i\bW,\alpha (\ba^i)^2 \bW^2)$ where $(.)^2$ is again a component-wise square, and $\alpha = p/(1-p)$ where $p$ is the dropout rate in original dropout. Gaussian dropout was shown to empirically perform as good as regular dropout, but generally with much faster convergence. 
%Gaussian uncertainty over $\bb^i$ can be straightforwardly reparameterized as Gaussian uncertainty over the weights $\bW$. See appendix~\ref{sec:dropout} for a derivation. 
%Gaussian dropout's implicit Gaussian weight uncertainty has an interpretation as a Gaussian variational posterior $\q(w_j) = \mathcal{N}(\theta_j, \alpha \theta_j^2)$, with fixed $\alpha=p/(1-p)$.

%As we will explain, Gaussian dropout corresponds to variational dropout with this particular Gaussian posterior and the \emph{log-uniform prior}. With the \emph{log-uniform prior}, we denote a prior for each individual parameter $w_j$: $p(\text{sign}(w_j)) = .5$ and $p(\log |w_j|) = \mathcal{U}(\text{interval})$, where the interval covers the complete log-domain of the magnitude of floating point numbers. Appendix~\ref{sec:floatingpoint} explains the relationship of the log-uniform distribution with floating-point numbers.

%Of key importance is that $D_{KL}(\q(\bw)||p(\bw))$, for this choice of prior and posterior, is just a function of $\alpha$. See appendix appendix~\ref{sec:kldiv} for proof. Since $\alpha=p/(1-p)$ is constant, $D_{KL}(\q(\bw)||p(\bw))$ is constant. As a result, all that is required for maximization of the bound (equation ~\eqref{eq:bound}), with this Gaussian posterior and a log-uniform prior is the maximization of $L_\D(\bphi)$ (equation ~\eqref{eq:bound2}); the rest of the terms are constant. The term $L_\D(\bphi)$ is exactly the term that is maximized by Gaussian dropout. 

%Thus, we have shown that optimizing a model with Gaussian dropout is a special case of variational dropout; i.e., SGVB with the local reparameterization trick, a log-uniform prior and a posterior on weights $w_j$: $\q(w_j) = \mathcal{N}(\theta_j, \alpha \theta_j^2)$, with fixed $\alpha = p/(1-p)$ where $p$ is the dropout rate. 

%\section{Examples of Variational Dropout}

%Despite its rigid implicit posterior variance (as we have shown), standard Gaussian dropout is known as an effective regularizer. For this reason, we will use standard Gaussian dropout as a starting point for exploration of different types of variational dropout.

%\subsection{Diagonal Gaussian Posterior}
%One promising direction is to keep the log-uniform prior and Gaussian posterior, and infer both a posterior mean and variance for each individual weight $w_j$, as opposed to just a mean:
%\begin{align}
%q_{\phi_j}(w_j) &= \mathcal{N}(\theta_j,\sigma_j^2)
%\text{\quad where\quad} \sigma_j^2 = \exp(\beta_j) \theta_j^2
%\end{align}
%where we parameterize $\alpha_j = \exp(\beta_j)$ such that $\phi_j = (\theta_j,\beta_j)$ are the variational parameters for weight $w_j$. As shown in appendix~\ref{sec:kldiv}, the KL-divergence between this Gaussian $\q(w_j)$ and log-uniform $p(w_j)$ is intractable, but can be extremely closely approximated by the 3rd degree polynomial:
%\begin{align}
%-D_{KL}[\q(w_j)|p(w_j)] \approx \text{constant} + 0.5\log(\alpha_j) + c_{1}\alpha_j + c_{2}\alpha_j^{2} + c_{3}\alpha_j^{3}
%\end{align}
%The coefficients $c1$, $c2$ and $c3$ are given in the appendix.

%Say something about relatation to Jeffrey's prior? Also scale-invariant, but our log-uniform prior is not improper since it has finite domain and thus normalizable.

%\subsection{Old: log-uniform posterior}
 %We parameterize the approximate posterior for each $w_j$ as $\q(\log |w_j|) = \mathcal{U}(\log |\theta_j|,\log |\theta_j|+ \exp(\delta_j))$, with variational parameters $\bphi=\{\bt,\bdelta\}$: a location parameter $\theta_j$ in regular space (including a sign), coupled with a scale of $\exp(\delta_j)$ in the logarithm of its magnitude. We constrain each $\delta_j$ to be non-positive for reasons that we will explain.
%We parameterize the scale in log-space of each $\q(w_j)$ as $\exp(\delta_j)$. 
%This results in an extraordinarily simple KL-divergence with the log-uniform prior: $D_{KL}(\q(|w_j|)||p(|w_j|)) = - \delta_j$. Note the scale-invariant property of this divergence: it is constant w.r.t. $\bt$.  When plugging this choice of $\q(\bw)$ and $p(\bw)$ in the variational bound of eq.~\eqref{eq:bound}, we arrive at our objective function:
%\begin{align}
%\mathcal{L}(\bphi) = \sum_j \delta_j + \sum_{(\bx, \by) \in \D} \Exp{\q(\bw)}{ \log p(\by|\bx,\bw)}
%\label{eq:objective}
%\end{align}
%where the scale parameter $\delta_j$ is constrained to be non-positive. We can use the parameterization $\bw_j = \theta_j \cdot \exp(\epsilon_j \cdot \delta_j)$ using uniform random draws $\epsilon_j \sim \mathcal{U}(0,1)$ for a Monte Carlo estimator of this bound.


\section{Variational Dropout}
\label{sec:gauss_dropout}
\textit{Dropout} is a technique for regularization of neural network parameters, which works by adding multiplicative noise to the input of each layer of the neural network during optimization. Using the notation of section~\ref{sec:fastsgvb}, for a fully connected neural network dropout corresponds to:
\begin{align}
\label{eq:dropout}
\bB = (\bA \circ \xi)\theta, \hspace{0.1cm} \text{ with }\hspace{0.1cm} \xi_{i,j} \sim p(\xi_{i,j})
\end{align}
where $\bA$ is the $M \times K$ matrix of input features for the current minibatch, $\theta$ is a $K \times L$ weight matrix, and $\bB$ is the $M \times L$ output matrix for the current layer (before a nonlinearity is applied). The $\circ$ symbol denotes the elementwise (Hadamard) product of the input matrix with a $M \times K$ matrix of independent noise variables $\xi$. By adding noise to the input during training, the weight parameters $\theta$ are less likely to overfit to the training data, as shown empirically by previous publications. Originally, \cite{hinton2012improving} proposed drawing the elements of $\xi$ from a Bernoulli distribution with probability $1-p$, with $p$ the \textit{dropout rate}. Later it was shown that using a continuous distribution with the same relative mean and variance, such as a Gaussian $N(1,\alpha)$ with $\alpha = p/(1-p)$, works as well or better \cite{srivastava2014dropout}.

Here, we re-interpret dropout with continuous noise as a variational method, and propose a generalization that we call \emph{variational dropout}. In developing variational dropout we provide a firm Bayesian justification for dropout training by deriving its implicit prior distribution and variational objective. This new interpretation allows us to propose several useful extensions to dropout, such as a principled way of making the normally fixed dropout rates $p$ \emph{adaptive} to the data.

\subsection{Variational dropout with independent weight noise} \label{sec:dropout_indep}
If the elements of the noise matrix $\xi$ are drawn independently from a Gaussian $N(1,\alpha)$, the marginal distributions of the activations $b_{m,j} \in \bB$ are Gaussian as well:
\begin{align}
\label{eq:marg_act}
\q(b_{m,j}|\bA) = N(\gamma_{m,j},\delta_{m,j}), \text{ with } \gamma_{m,j} = \sum_{i=1}^{K} a_{m,i}\theta_{i,j}, \text{ and } \delta_{m,j} = \alpha \sum_{i=1}^{K} a^{2}_{m,i}\theta_{i,j}^{2}.
\end{align}
Making use of this fact, \cite{wang2013fast} proposed \emph{Gaussian dropout}, a regularization method where, instead of applying \eqref{eq:dropout}, the activations are directly drawn from their (approximate or exact) marginal distributions as given by \eqref{eq:marg_act}. \cite{wang2013fast} argued that these marginal distributions are exact for Gaussian noise $\xi$, and for Bernoulli noise still approximately Gaussian because of the central limit theorem. This ignores the dependencies between the different elements of $\bB$, as present using \eqref{eq:dropout}, but \cite{wang2013fast} report good results nonetheless.  %Gaussian dropout was shown empirically to perform as well as regular dropout, but generally with faster convergence.

As noted by \cite{wang2013fast}, and explained in appendix~\ref{sec:dropout}, this Gaussian dropout noise can also be interpreted as arising from a Bayesian treatment of a neural network with weights $\bW$ that multiply the input to give $\bB = \bA\bW$, where the posterior distribution of the weights is given by a factorized Gaussian with $\q(w_{i,j}) = \mathcal{N}(\theta_{i,j}, \alpha \theta_{i,j}^2)$. From this perspective, the marginal distributions \eqref{eq:marg_act} then arise through the application of the local reparameterization trick, as introduced in section~\ref{sec:fastsgvb}. The variational objective corresponding to this interpretation is discussed in section~\ref{sec:dropout_prior}.

\subsection{Variational dropout with correlated weight noise}\label{sec:dropout_corr}
Instead of ignoring the dependencies of the activation noise, as in section~\ref{sec:dropout_indep}, we may retain the dependencies by interpreting dropout \eqref{eq:dropout} as a form of \emph{correlated weight noise}:
\begin{align}
\bB = (\bA \circ \xi)\theta, \hspace{0.1cm} \xi_{i,j} \sim N(1,\alpha) \hspace{0.2cm} \Longleftrightarrow \hspace{0.2cm} \bb^{m} = \ba^{m}\bW, \text{ with }\nonumber\\
\bW = (\bw_{1}',\bw_{2}',\ldots,\bw_{K}')', \text{ and } \bw_{i} = s_{i}\theta_{i}, \text{ with } \q(s_{i}) = N(1,\alpha),
\end{align}
where $\ba^{m}$ is a row of the input matrix and $\bb^{m}$ a row of the output. The $\bw_{i}$ are the rows of the weight matrix, each of which is constructed by multiplying a non-stochastic parameter vector $\theta_{i}$ by a stochastic scale variable $s_{i}$. The distribution on these scale variables we interpret as a Bayesian posterior distribution. The weight parameters $\theta_{i}$ (and the biases) are estimated using maximum likelihood. The original Gaussian dropout sampling procedure \eqref{eq:dropout} can then be interpreted as arising from a local reparameterization of our posterior on the weights $\bW$.

\subsection{Dropout's scale-invariant prior and variational objective} \label{sec:dropout_prior}
The posterior distributions $\q(\bW)$ proposed in sections~\ref{sec:dropout_indep} and \ref{sec:dropout_corr} have in common that they can be decomposed into a parameter vector $\theta$ that captures the mean, and a multiplicative noise term determined by parameters $\alpha$. Any posterior distribution on $\bW$ for which the noise enters this multiplicative way, we will call a \emph{dropout posterior}. Note that many common distributions, such as univariate Gaussians (with nonzero mean), can be reparameterized to meet this requirement.

During dropout training, $\theta$ is adapted to maximize the expected log likelihood $\Exp{q_\alpha}{L_\D(\theta)}$. For this to be consistent with the optimization of a variational lower bound of the form in \eqref{eq:bound2}, the prior on the weights $p(\bw)$ has to be such that $D_{KL}(\q(\bw)||p(\bw))$ does not depend on $\theta$. In appendix~\ref{sec:kldiv} we show that the only prior that meets this requirement is the scale invariant log-uniform prior:
\[
p(\log(|w_{i,j}|)) \propto c,
\]
i.e. a prior that is uniform on the log-scale of the weights (or the weight-scales $s_{i}$ for section~\ref{sec:dropout_corr}). As explained in appendix~\ref{sec:floatingpoint}, this prior has an interesting connection with the floating point format for storing numbers: From an MDL perspective, the floating point format is optimal for communicating numbers drawn from this prior. Conversely, the KL divergence $D_{KL}(\q(\bw)||p(\bw))$ with this prior has a natural interpretation as regularizing the number of significant digits our posterior $\q$ stores for the weights $w_{i,j}$ in the floating-point format.% In section~\ref{sec:experiments} we examine the empirical performance of this prior.

Putting the expected log likelihood and KL-divergence penalty together, we see that dropout training maximizes the following variatonal lower bound w.r.t.\ $\theta$:
\begin{align}
\label{eq:dropout_bound}
\Exp{q_{\alpha}}{L_\D(\theta)} - D_{KL}(q_{\alpha}(\bw)||p(\bw)),  
\end{align}
where we have made the dependence on the $\theta$ and $\alpha$ parameters explicit. The noise parameters $\alpha$ (e.g.\ the dropout rates) are commonly treated as hyperparameters that are kept fixed during training. For the log-uniform prior this then corresponds to a fixed limit on the number of significant digits we can learn for each of the weights $w_{i,j}$. In section~\ref{sec:adaptive_dropout} we discuss the possibility of making this limit adaptive by also maximizing the lower bound with respect to $\alpha$.

For the choice of a factorized Gaussian approximate posterior with $\q(w_{i,j}) = \mathcal{N}(\theta_{i,j}, \alpha \theta_{i,j}^2)$, as discussed in section~\ref{sec:dropout_indep}, the lower bound~\eqref{eq:dropout_bound} is analyzed in detail in appendix~\ref{sec:kldiv}. There, it is shown that for this particular choice of posterior the negative KL-divergence $-D_{KL}(q_{\alpha}(\bw)||p(\bw))$ is not analytically tractable, but can be approximated extremely accurately using
\[
-D_{KL}[\q(w_i)|p(w_i)] \approx \text{constant} + 0.5\log(\alpha) + c_{1}\alpha + c_{2}\alpha^{2} + c_{3}\alpha^{3},
\]
with
\[
c_{1} = 1.16145124, \hspace{0.5cm} c_{2} = -1.50204118, \hspace{0.5cm} c_{3} = 0.58629921.
\]
The same expression may be used to calculate the corresponding term $-D_{KL}(q_{\alpha}(\mathbf{s})||p(\mathbf{s}))$ for the posterior approximation of section~\ref{sec:dropout_corr}.

\subsection{Adaptive regularization through optimizing the dropout rate}\label{sec:adaptive_dropout}
The noise parameters $\alpha$ used in dropout training (e.g.\ the dropout rates) are usually treated as fixed hyperparameters, but now that we have derived dropout's variational objective \eqref{eq:dropout_bound}, making these parameters adaptive is trivial: simply maximize the variational lower bound with respect to $\alpha$. We can use this to learn a separate dropout rate per layer, per neuron, of even per separate weight. In section~\ref{sec:experiments} we look at the predictive performance obtained by making $\alpha$ adaptive.

We found that very large values of $\alpha$ correspond to local optima from which it is hard to escape due to large-variance gradients. To avoid such local optima, we found it beneficial to set a constraint $\alpha \leq 1$ during training, i.e.\ we maximize the posterior variance at the square of the posterior mean, which corresponds to a dropout rate of $0.5$.


%\subsection{Gaussian Dropout}
%In \cite{wang2013fast} it was shown that regularization by regular (binary) dropout \cite{srivastava2014dropout} has a fast approximation called \emph{Gaussian dropout}. In Gaussian dropout, linear operations $\bb^i = \ba^i \bW$ are replaced by a stochastic relation $p(\bb^i|\ba^i) = \mathcal{N}(\ba^i\bW,\alpha (\ba^i)^2 \bW^2)$ where $(.)^2$ is again a component-wise square, and $\alpha = p/(1-p)$ where $p$ is the dropout rate in original dropout. Gaussian dropout was shown to empirically perform as good as regular dropout, but generally with much faster convergence. 
%Gaussian uncertainty over $\bb^i$ can be straightforwardly reparameterized as Gaussian uncertainty over the weights $\bW$. See appendix~\ref{sec:dropout} for a derivation. 
%Gaussian dropout's implicit Gaussian weight uncertainty has an interpretation as a Gaussian variational posterior $\q(w_j) = \mathcal{N}(\theta_j, \alpha \theta_j^2)$, with fixed $\alpha=p/(1-p)$.

%As we will explain, Gaussian dropout corresponds to variational dropout with this particular Gaussian posterior and the \emph{log-uniform prior}. With the \emph{log-uniform prior}, we denote a prior for each individual parameter $w_j$: $p(\text{sign}(w_j)) = .5$ and $p(\log |w_j|) = \mathcal{U}(\text{interval})$, where the interval covers the complete log-domain of the magnitude of floating point numbers. Appendix~\ref{sec:floatingpoint} explains the relationship of the log-uniform distribution with floating-point numbers.

%Of key importance is that $D_{KL}(\q(\bw)||p(\bw))$, for this choice of prior and posterior, is just a function of $\alpha$. See appendix appendix~\ref{sec:kldiv} for proof. Since $\alpha=p/(1-p)$ is constant, $D_{KL}(\q(\bw)||p(\bw))$ is constant. As a result, all that is required for maximization of the bound (equation ~\eqref{eq:bound}), with this Gaussian posterior and a log-uniform prior is the maximization of $L_\D(\bphi)$ (equation ~\eqref{eq:bound2}); the rest of the terms are constant. The term $L_\D(\bphi)$ is exactly the term that is maximized by Gaussian dropout. 

%Thus, we have shown that optimizing a model with Gaussian dropout is a special case of variational dropout; i.e., SGVB with the local reparameterization trick, a log-uniform prior and a posterior on weights $w_j$: $\q(w_j) = \mathcal{N}(\theta_j, \alpha \theta_j^2)$, with fixed $\alpha = p/(1-p)$ where $p$ is the dropout rate. 

%\subsection{Examples of Variational Dropout}

%Despite its rigid implicit posterior variance (as we have shown), standard Gaussian dropout is known as an effective regularizer. For this reason, we will use standard Gaussian dropout as a starting point for exploration of various types of variational dropout.

%\subsection{Diagonal Gaussian Posterior}
%One promising direction is to keep the log-uniform prior and Gaussian posterior, and infer both a posterior mean and variance for each individual weight $w_j$, as opposed to just a mean:
%\begin{align}
%q_{\phi_j}(w_j) &= \mathcal{N}(\theta_j,\sigma_j^2)
%\text{\quad where\quad} \sigma_j^2 = \exp(\beta_j) \theta_j^2
%\end{align}
%where we parameterize $\alpha_j = \exp(\beta_j)$ such that $\phi_j = (\theta_j,\beta_j)$ are the variational parameters for weight $w_j$. As shown in appendix~\ref{sec:kldiv}, the KL-divergence between this Gaussian $\q(w_j)$ and log-uniform $p(w_j)$ is intractable, but can be extremely closely approximated by the 3rd degree polynomial:
%\begin{align}
%-D_{KL}[\q(w_j)|p(w_j)] \approx \text{constant} + 0.5\log(\alpha_j) + c_{1}\alpha_j + c_{2}\alpha_j^{2} + c_{3}\alpha_j^{3}
%\end{align}
%The coefficients $c_1$, $c_2$ and $c_3$ are given in the appendix.

%Say something about relatation to Jeffrey's prior? Also scale-invariant, but our log-uniform prior is not improper since it has finite domain and thus normalizable.

%\subsection{Old: log-uniform posterior}
 %We parameterize the approximate posterior for each $w_j$ as $\q(\log |w_j|) = \mathcal{U}(\log |\theta_j|,\log |\theta_j|+ \exp(\delta_j))$, with variational parameters $\bphi=\{\bt,\bdelta\}$: a location parameter $\theta_j$ in regular space (including a sign), coupled with a scale of $\exp(\delta_j)$ in the logarithm of its magnitude. We constrain each $\delta_j$ to be non-positive for reasons that we will explain.
%We parameterize the scale in log-space of each $\q(w_j)$ as $\exp(\delta_j)$. 
%This results in an extraordinarily simple KL-divergence with the log-uniform prior: $D_{KL}(\q(|w_j|)||p(|w_j|)) = - \delta_j$. Note the scale-invariant property of this divergence: it is constant w.r.t. $\bt$.  When plugging this choice of $\q(\bw)$ and $p(\bw)$ in the variational bound of eq.~\eqref{eq:bound}, we arrive at our objective function:
%\begin{align}
%\mathcal{L}(\bphi) = \sum_j \delta_j + \sum_{(\bx, \by) \in \D} \Exp{\q(\bw)}{ \log p(\by|\bx,\bw)}
%\label{eq:objective}
%\end{align}
%where the scale parameter $\delta_j$ is constrained to be non-positive. We can use the parameterization $\bw_j = \theta_j \cdot \exp(\epsilon_j \cdot \delta_j)$ using uniform random draws $\epsilon_j \sim \mathcal{U}(0,1)$ for a Monte Carlo estimator of this bound.


\section{Related Work}\label{sec:relatedwork}
Pioneering work in practical variational inference for neural networks was done in~\cite{graves2011practical}, where a (biased) variational lower bound estimator was introduced with good results on recurrent neural network models. In later work~\cite{kingma2013auto,rezende2014stochastic} it was shown that even more practical estimators can be formed for most types of continuous latent variables or parameters using a (non-local) reparameterization trick, leading to efficient and unbiased stochastic gradient-based variational inference. These works focused on an application to latent-variable inference; extensive empirical results on inference of global model parameters were reported in ~\cite{blundell2015weight}, including succesful application to reinforcement learning. These earlier works used the relatively high-variance estimator~\eqref{eq:objective_sgvb}, upon which we improve. Variable reparameterizations have a long history in the statistics literature, but have only recently found use for efficient \emph{gradient-based} machine learning and inference~\cite{bengio2013estimating,kingma2013fast,salimans2013fixedform}. Related is also \emph{probabilistic backpropagation}~\cite{hernandez2015probabilistic}, an algorithm for inferring marginal posterior probabilities; however, it requires certain tractabilities in the network making it insuitable for the type of models under consideration in this paper.

As we show here, regularization by dropout~\cite{srivastava2014dropout, wang2013fast} can be interpreted as variational inference. DropConnect~\cite{wan2013regularization} is similar to dropout, but with binary noise on the weights rather than hidden units. DropConnect thus has a similar interpretation as variational inference, with a uniform prior over the weights, and a mixture of two Dirac peaks as posterior. In ~\cite{ba2013adaptive}, \emph{standout} was introduced, a variation of dropout where a binary belief network is learned for producing dropout rates. Recently, \cite{maeda2014bayesian} proposed another Bayesian perspective on dropout. In recent work~\cite{bayer2015fast}, a similar reparameterization is described and used for variational inference; their focus is on closed-form approximations of the variational bound, rather than unbiased Monte Carlo estimators. ~\cite{maeda2014bayesian} and \cite{gal2015dropout} also investigate a Bayesian perspective on dropout, but focus on the binary variant. ~\cite{gal2015dropout} reports various encouraging results on the utility of dropout's implied prediction uncertainty.

%We show how SGVB can be used to infer low-precision floating-point parameters; a potential application is low-precision hardware\cite{courbariaux2014low}.

%"Building Blocks for Variational Bayesian Learning of Latent Variable Models " (Tapani et al)”

\section{Experiments}\label{sec:experiments}
We compare our method to standard binary dropout and two popular versions of Gaussian dropout, which we'll denote with type A and type B. With Gaussian dropout type A we denote the pre-linear Gaussian dropout from ~\cite{srivastava2014dropout}; type B denotes the post-linear Gaussian dropout from ~\cite{wang2013fast}. This way, the method names correspond to the matrix names in section 2 ($\bA$ or $\bB$) where noise is injected. Models were implemented in Theano~\cite{bergstra2010theano}, and optimization was performed using Adam~\cite{kingma2015adam} with default hyper-parameters and temporal averaging.

Two types of variational dropout were included. Type A is correlated weight noise as introduced in section~\ref{sec:dropout_corr}: an adaptive version of Gaussian dropout type A. Variational dropout type B has independent weight uncertainty as introduced in section~\ref{sec:dropout_indep}, and corresponds to Gaussian dropout type B.

A \emph{de facto} standard benchmark for regularization methods is the task of MNIST hand-written digit classification. We choose the same architecture as~\cite{srivastava2014dropout}: a fully connected neural network with 3 hidden layers and rectified linear units (ReLUs). We follow the dropout hyper-parameter recommendations from these earlier publications, which is a dropout rate of $p=0.5$ for the hidden layers and $p=0.2$ for the input layer. We used early stopping with all methods, where the amount of epochs to run was determined based on performance on a validation set. 

\paragraph{Variance.} We start out by empirically comparing the variance of the different available stochastic estimators of the gradient of our variational objective. To do this we train the neural network described above for either 10 epochs (test error 3\%) or 100 epochs (test error 1.3\%), using variational dropout with independent weight noise. After training, we calculate the gradients for the weights of the top and bottom level of our network on the full training set, and compare against the gradient estimates per batch of $M=1000$ training examples. Appendix~\ref{sec:variance_corr} contains the same analysis for the case of variational dropout with correlated weight noise.

Table \ref{tab:mnist_var_indep} shows that the local reparameterization trick yields the lowest variance among all variational dropout estimators for all conditions, although it is still substantially higher compared to not having any dropout regularization. The $1/M$ variance scaling achieved by our estimator is especially important early on in the optimization when it makes the largest difference (compare \textit{weight sample per minibatch} and \textit{weight sample per data point}). The additional variance reduction obtained by our estimator through drawing fewer random numbers (section~\ref{sec:fastsgvb}) is about a factor of 2, and this remains relatively stable as training progresses (compare \textit{local reparameterization} and \textit{weight sample per data point}).
\begin{table}[htb]
\begin{tabular}{lllll}
 & top layer & top layer & bottom layer & bottom layer \\
stochastic gradient estimator & 10 epochs & 100 epochs & 10 epochs & 100 epochs \\
\hline
local reparameterization (ours) & $7.8\times 10^{3}$ & $1.2\times 10^{3}$ & $1.9\times 10^{2}$ & $1.1\times 10^{2}$ \\
weight sample per data point (slow) & $1.4\times 10^{4}$ & $2.6\times 10^{3}$ & $4.3\times 10^{2}$ & $2.5\times 10^{2}$ \\
weight sample per minibatch (standard) & $4.9\times 10^{4}$ & $4.3\times 10^{3}$ & $8.5\times 10^{2}$ & $3.3\times 10^{2}$ \\
no dropout noise (minimal var.) & $2.8\times 10^{3}$ & $5.9\times 10^{1}$ & $1.3\times 10^{2}$ & $9.0\times 10^{0}$
\end{tabular}
\caption{Average empirical variance of minibatch stochastic gradient estimates (1000 examples) for a fully connected neural network, regularized by variational dropout with independent weight noise.}
\label{tab:mnist_var_indep}
\end{table}
\FloatBarrier
%\begin{table}[htb]
%\begin{tabular}{lllll}
% & top layer & top layer & bottom layer & bottom layer \\
%stochastic gradient estimator & 10 epochs & 100 epochs & 10 epochs & 100 epochs \\
%\hline
%local reparameterization (ours) & $3.2\times 10^{4}$ & $5.3\times 10^{3}$ & $8.2\times 10^{2}$ & $5.2\times 10^{2}$ \\
%separate weight samples (slow) & $6.6\times 10^{4}$ & $1.2\times 10^{4}$ & $2.1\times 10^{3}$ & $1.2\times 10^{3}$ \\
%single weight sample (standard) & $1.1\times 10^{5}$ & $1.3\times 10^{4}$ & $2.8\times 10^{3}$ & $1.3\times 10^{3}$\\
%no dropout noise (minimal var.) & $1.3\times 10^{4}$ & $3.0\times 10^{2}$ & $6.4\times 10^{2}$ & $4.5\times 10^{1}$
%\end{tabular}
%\caption{Average empirical variance of minibatch stochastic gradient estimates (200 examples) for a fully connected neural network, regularized by variational dropout with independent weight noise.}
%\label{tab:mnist_var_indep}
%\end{table}
%\begin{table}[htb]
%\begin{tabular}{lllll}
% & top layer & top layer & bottom layer & bottom layer \\
%stochastic gradient estimator & 10 epochs & 100 epochs & 10 epochs & 100 epochs \\
%\hline
%local reparameterization (ours) & $1.3\times 10^{5}$ & $2.9\times 10^{4}$ & $4.0\times 10^{3}$ & $2.7\times 10^{3}$ \\
%single weight sample (standard) & $4.4\times 10^{5}$ & $4.8\times 10^{4}$ & $9.3\times 10^{3}$ & $3.4\times 10^{3}$ \\
%no dropout noise (minimal var.) & $3.5\times 10^{5}$ & $4.7\times 10^{3}$ & $2.1\times 10^{3}$ & $5.9\times 10^{2}$
%\end{tabular}
%\caption{Average empirical variance of minibatch stochastic gradient estimates (200 examples) for a fully connected neural network, regularized by variational dropout with correlated weight noise. For this specification, the variance of the stochastic gradient estimator with separate weight samples is identical to that using the local reparameterization trick: Ignoring optimization of dropout rates, both methods correspond to the standard dropout training procedure.}
%\label{tab:mnist_var_corr}
%\end{table}

\paragraph{Speed.} We compared the regular SGVB estimator, with separate weight samples per datapoint with the efficient estimator based on local reparameterization, in terms of wall-clock time efficiency. With our implementation on a modern GPU, optimization with the na\"ive estimator took \textbf{1635 seconds} per epoch, while the efficient estimator took \textbf{7.4 seconds}: an over 200 fold speedup.

\paragraph{Classification error.} Figure~\ref{fig:mnist_classification} shows test-set classification error for the tested regularization methods, for various choices of number of hidden units. Our adaptive variational versions of Gaussian dropout perform equal or better than their non-adaptive counterparts and standard dropout under all tested conditions. The difference is especially noticable for the smaller networks. In these smaller networks, we observe that variational dropout infers dropout rates that are on average far lower than the dropout rates for larger networks. This adaptivity comes at negligable computational cost.

We found that slightly downscaling the KL divergence part of the variational objective can be beneficial. \emph{Variational (A2)} in figure~\ref{fig:mnist_classification} denotes performance of type A variational dropout but with a KL-divergence downscaled with a factor of 3; this small modification seems to prevent underfitting, and beats all other dropout methods in the tested models.



\begin{figure}[t]
  \vspace{-1cm}
\centering
\begin{subfigure}[t]{0.49\textwidth}\centering
\includegraphics[width=7cm]{graphs_classification/mnist}
\caption{Classification error on the MNIST dataset}
\end{subfigure}
\begin{subfigure}[t]{0.49\textwidth}\centering
\includegraphics[width=7cm]{graphs_classification/cifar10}
\caption{Classification error on the CIFAR-10 dataset}
\end{subfigure}
\caption{Best viewed in color.  \textbf{(a)} Comparison of various dropout methods, when applied to fully-connected neural networks for classification on the MNIST dataset. Shown is the classification error of networks with 3 hidden layers, averaged over 5 runs. he variational versions of Gaussian dropout perform equal or better than their non-adaptive counterparts; the difference is especially large with smaller models, where regular dropout often results in severe underfitting. \textbf{(b)} Comparison of dropout methods when applied to convolutional net a trained on the CIFAR-10 dataset, for different settings of network size $k$. The network has two convolutional layers with each $32 k$ and $64 k$ feature maps, respectively, each with stride 2 and followed by a softplus nonlinearity. This is followed by two fully connected layers with each $128 k$ hidden units.} \label{fig:mnist_classification}
\end{figure}


\FloatBarrier

\section{Conclusion}
Efficiency of posterior inference using stochastic gradient-based variational Bayes (SGVB) can often be significantly improved through a \emph{local reparameterization} where global parameter uncertainty is translated into local uncertainty per datapoint. By injecting noise locally, instead of globally at the model parameters, we obtain an efficient estimator that has low computational complexity, can be trivially parallelized and has low variance. We show how dropout is a special case of SGVB with local reparameterization, and suggest \emph{variational dropout}, a straightforward extension of regular dropout where optimal dropout rates are inferred from the data, rather than fixed in advance. We report encouraging empirical results.

\subsubsection*{Acknowledgments}
We thank the reviewers and Yarin Gal for valuable feedback. Diederik Kingma is supported by the Google European Fellowship in Deep Learning, Max Welling is supported by research grants from Google and Facebook, and the NWO project in Natural AI (NAI.14.108).

\small
\bibliographystyle{apalike}
\bibliography{bib}

\appendix

\newpage
\section{Floating-point numbers and compression}\label{sec:floatingpoint}
The \emph{floating-point} format is by far the most common number format used for model parameters in neural network type models. The absolute value of floating-point numbers in base 2 equals $\text{\emph{sign}} \cdot (s / 2^{p-1}) \cdot 2^e$ where $s$ is the mantissa with $p$ bits and $e$ is the exponent, both non-negative integer-valued. For random bit patterns for $s$ and $e$, floating point numbers closely follow a log-uniform distribution with finite range. The floating-point format is thus close to optimal, from a compression point of view, for communicating parameters for which the receiver has a log-uniform prior distribution, with finite range, over the magnitude of parameters. Specifically, this prior for each individual floating-point parameter $w_j$ is $p(\log |w_j|) = \mathcal{U}(\text{interval})$, where the interval is approximately $(\log(2^{-8}),\log(2^8))$ for single-point precision and $(\log(2^{-11}),\log(2^{11}))$ for double-precision floating-point numbers in the most common IEEE specification. 

The KL-divergence between a variational posterior $\q$ (eq.~\eqref{eq:bound}) and this log-uniform prior is relatively simple: it equals the entropy under $\q$ of the log-transformed magnitude of $\bw$ plus a constant:
\begin{align}
- D_{KL}(\q(\bw)||p(\bw)) = \mathcal{H}(\q(\log |\bw|)) + \text{constant}
\end{align}
This divergence has a natural interpretation as controlling the number of significant digits under $\q$ of the corresponding floating-point number.

\section{Derivation of dropout's implicit variational posterior}\label{sec:dropout}
In the common version of dropout~\cite{srivastava2014dropout}, linear operations $\bb^i = \ba^i \bW$ in a model, where $\ba^i$ is a column vector that is the result of previous operations on the input for datapoint $i$, are replaced by a stochastic operation:
\begin{align}
\bb^i &= (\ba^i \odot (\bd^i/(1-p)) \bW \\
\text{i.e.:\quad} b_k^i &= \sum_j W_{jk} a_j^i d_j^i /(1-p)\\
d_j^i &\sim \text{Bernoulli}(1-p)\\
p &= \text{dropout rate (hyper-parameter, e.g. 0.5)}
\end{align}
where column vector $\ba^i$ is the result of previous operations on the input for datapoint $i$. Expected values and variance w.r.t. the dropout rate are:
\begin{align}
\Exp{}{d_j^i} &= (1-p)\\
\Exp{}{d_j^i/(1-p)} &= 1\\
\Var{}{d_j^i/(1-p)} &= \Var{}{d_j^i}/(1-p)^2 = p/(1-p)\\
\end{align}
The expected value and variance for the elements of the resulting vector $\bb^i$ are:
\begin{align}
\Exp{}{b_k^i} &= \Exp{}{\sum_j W_{jk} a_j^i d_j^i / (1-p)}\\
&= \sum_j W_{jk} a_j^i \Exp{}{d_j^i / (1-p)}\\
&= \sum_j W_{jk} a_j^i\\
\Var{}{b_k^i} &= \Var{}{\sum_j W_{jk} a_j^i d_j^i / (1-p)}\\
&= \sum_j \Var{}{W_{jk} a_j^i d_j^i / (1-p)}\\
&= \sum_j W_{jk}^2 (a_j^i)^2 \Var{}{d_j^i / (1-p)}\\
&= p/(1-p) \sum_j W_{jk}^2 (a_j^i)^2
\end{align}
\subsection{Gaussian dropout}
Thus, ignoring dependencies between elements of $\bb^i$, dropout can be approximated by a Gaussian ~\cite{wang2013fast}:
\begin{align}
p(\bb^i|\ba^i,W) &= \mathcal{N}(\ba^i W,  p(1-p)^{-1} (\ba^i)^2 (\bW)^2)\\
\text{i.e.:\quad} \bb^i &= \ba^i W + \beps^i \odot \sqrt{p(1-p)^{-1} (\ba^i)^2 (\bW)^2}\\
\text{where:\quad} \beps &\sim \mathcal{N}(0,I)
\label{eq:gaussiandropout}
\end{align}
where $(.)^2$, $\sqrt(.)$ and $\odot$ are again component-wise operations.
\subsection{Weight uncertainty}
Equivalently, the Gaussian dropout relationship $p(\bb^i|\ba^i,W)$ can be parameterized as weight uncertainty:
\begin{align}
\bb^i &= \ba^i \bV\\
p(V_{jk}|W_{jk},\alpha_{jk}) &= \mathcal{N}(W_{jk}, \alpha_{jk} W_{jk}^2)\\
\alpha_{jk} &= p/(1-p)\\
\Rightarrow\text{\quad} p(\bb^i|\ba^i,W) &= \mathcal{N}(\ba^i W,  p(1-p)^{-1} (\ba^i)^2 (\bW)^2)
\end{align}
In the main text, we show how we can treat the distribution over weights $p(V_{jk}|W_{jk},\alpha_{jk})$ as a variational posterior $\q$, which we can optimize w.r.t. each $\alpha_{jk}$ (the dropout rate for individual parameters) to optimize a bound on the marginal likelihood of the data.

\section{Negative KL-divergence for the log-uniform prior}\label{sec:kldiv}
The variational lower bound \eqref{eq:bound} that we wish to optimize is given by
\[
\mathcal{L}(\bphi) = - D_{KL}(\q(\bw)||p(\bw)) + L_\D(\bphi),
\]
where $L_\D(\bphi)$ is the expected log likelihood term which we approximate using Monte Carlo, as discussed in section~\ref{sec:fastsgvb}. Here we discuss how to deterministically compute the negative KL divergence term $-D_{KL}(\q(\bw)||p(\bw))$ for the \textit{log-uniform} prior proposed in section~\ref{sec:dropout}. For simplicity we only discuss the case of factorizing approximate posterior $\q(\bw) = \prod_{i=1}^{k} \q(w_{i})$, although the extension to more involved posterior approximations is straight forward.

Our log-uniform prior $p(w_{i})$ (appendix \ref{sec:floatingpoint}) consists of a Bernoulli distribution on the \textit{sign bit} $s_{i} \in \{-1,+1\}$ which captures the sign of $w_{i}$,  combined with a uniform prior on the log-scale of $w_{i}$:
\begin{eqnarray}
w_{i} & = & s_{i}|w_{i}| \nonumber\\
p(s_{i}) & = & \text{Bernoulli}(0.5) \text{ on } \{-1,+1\}\nonumber\\
p(\log(|w_{i}|)) & \propto & c
\end{eqnarray}
For an approximate posterior $\q(w_{i})$ with support on the real line, the KL divergence will thus consist of two parts: the KL-divergence between the posterior and prior on the sign bit $s_{i}$ and the KL-divergence between the conditional posterior and prior on the log scale $\log(|w_{i}|)$.

The negative KL divergence between the posterior and prior on the sign bit is given by:
\begin{align}
-D_{KL}[\q(s_{i})|p(s_{i})] = \log(0.5) - Q_{\phi}(0)\log[Q_{\phi}(0)] - [1-Q_{\phi}(0)]\log[1-Q_{\phi}(0)],
\end{align}
with $Q_{\phi}(0)$ the posterior CDF of $w_{i}$ evaluated at 0, i.e. the probability of drawing a negative weight from our approximate posterior.

The negative KL divergence of the second part (the log-scale) is then
\begin{align}
-D_{KL}[\q(\log(|w_{i}|) | s_{i})|p(\log(|w_{i}|))] = \log(c) - Q_{\phi}(0)\mathbb{E}_{\q(w_{i}|w_{i}<0)}\left[\log(\q(w_{i})/Q_{\phi}(0)) + \log(|w_{i}|)\right] \nonumber \\
 - \left(1-Q_{\phi}(0)\right)\mathbb{E}_{\q(w_{i}|w_{i}>0)}\left[\log(\q(w_{i})/(1-Q_{\phi}(0))) + \log(|w_{i}|)\right], \nonumber
\end{align}
where $\q(w_{i})/Q_{\phi}(0)$ is the truncation of our approximate posterior to the negative part of the real-line, and the $\log(|w_{i}|)$ term enters through transforming the uniform prior from the log-space to the normal space.

Putting both parts together the terms involving the CDF $Q_{\phi}(0)$ cancel, and we get:
\begin{align}
\label{eq:combined_kl}
-D_{KL}[\q(w_i)|p(w_i)] = \log(c) + \log(0.5) + \mathcal{H}(\q(w_{i})) - \mathbb{E}_{q} \log(|w_{i}|),
\end{align}
where $\mathcal{H}(\q(w_{i}))$ is the entropy of our approximate posterior.

We now consider the special case of a \textit{dropout posterior}, as proposed in section~\ref{sec:gauss_dropout}, which we defined as any approximate posterior distribution $\q(w_{i})$ for which the weight noise is \textit{multiplicative}:
\begin{align}
w_{i} = \theta_{i}\epsilon_{i}, \text{ with } \epsilon_{i} \sim q_{\alpha}(\epsilon_{i}),
\end{align}
where we have divided the parameters $\phi_{i} = (\theta_{i}, \alpha)$ into parameters governing the mean of $\q(w_i)$ (assuming $\mathbb{E}_{q} \epsilon_{i} = 1$) and the multiplicative noise. Note that certain families of distributions, such as Gaussians, can always be reparameterized in this way.

Assuming $q_{\alpha}(\epsilon_{i})$ is a continuous distribution, this choice of posterior gives us
\begin{align}
\label{eq:h_eps}
\mathcal{H}(\q(w_{i})) = \mathcal{H}(q_{\alpha}(\epsilon_{i})) + \log(|\theta_{i}|),
\end{align}
where the last term can be understood as the log jacobian determinant of a linear transformation of the random variable $\epsilon_{i}$.

Furthermore, we can rewrite the $-\mathbb{E}_{q} \log(|w_{i}|)$ term of \eqref{eq:combined_kl} as
\begin{align}
\label{eq:e_w}
-\mathbb{E}_{q} \log(|w_{i}|) =  -\mathbb{E}_{q} \log(|\theta_{i}\epsilon_{i}|) = -\log(|\theta_{i}|) - \mathbb{E}_{q} \log(|\epsilon_{i}|).
\end{align}
Taking \eqref{eq:h_eps} and \eqref{eq:e_w} together, the terms depending on $\theta$ thus cancel exactly, proving that the KL divergence between posterior and prior is indeed independent of these parameters:
\begin{align}
-D_{KL}[\q(w_i)|p(w_i)] & = \log(c) + \log(0.5) + \mathcal{H}(\q(w_{i})) - \mathbb{E}_{q} \log(|w_{i}|) \nonumber\\
& = \log(c) + \log(0.5) + \mathcal{H}(q_{\alpha}(\epsilon_{i})) - \mathbb{E}_{q_{\alpha}} \log(|\epsilon_{i}|).
\end{align}
In fact, the log-uniform prior is the only prior consistent with dropout. If the prior had any additional terms involving $w_{i}$, their expectation w.r.t.\ $\q$ would not cancel with the entropy term in $-D_{KL}[\q(w_i)|p(w_i)]$.

For the specific case of \textit{Gaussian dropout}, as proposed in section~\ref{sec:gauss_dropout}, we have $q_{\alpha}(\epsilon_{i}) = N(1,\alpha)$. This then gives us
\[
-D_{KL}[\q(w_i)|p(w_i)] = \log(c) + \log(0.5) + 0.5(1 + \log(2\pi) + \log(\alpha)) - \mathbb{E}_{q_{\alpha}} \log(|\epsilon_{i}|).
\]
%For this choice of posterior we have $\mathcal{H}(\q(w_{i})) = 0.5(1 + \log(2\pi) + \log(\alpha) + \log(\theta_{i}^{2}))$, and we can rewrite the last term in the KL divergence $-\mathbb{E}_{q} \log(|w_{i}|)$ using the \textit{reparameterization trick} as
%\begin{align}
%-\mathbb{E}_{q} \log(|w_{i}|) & = & -0.5\mathbb{E}_{q} \log(w_{i}^{2}) = -0.5\mathbb{E}_{\epsilon \sim \mathcal{N}(0,1)}\log[(\theta_{i}+\alpha\theta_{i}\epsilon)^{2}] \\
%& = & -0.5\log(\theta_{i}^{2}) - 0.5\mathbb{E}_{\epsilon \sim \mathcal{N}(0,1)}\log[(1 + \sqrt{\alpha}\epsilon)^{2}].
%\end{align}
%The terms depending on $\theta$ thus cancel exactly in $\mathcal{H}(\q(w_{i}))-\mathbb{E}_{q} \log(|w_{i}|)$, proving that the KL divergence between posterior and prior is indeed independent of these parameters.
%The final term $-0.5\mathbb{E}_{\epsilon \sim \mathcal{N}(0,1)}\log[(1 + \sqrt{\alpha}\epsilon)^{2}]$ is not analytically tractable, but it can be pre-computed numerically for all values of $\alpha$ in the relevant range and then approximated using a 3rd degree polynomial in $\alpha$. As shown in figure~\ref{fig:gausskl}, this approximation is extremely close.
The final term $\mathbb{E}_{q_{\alpha}} \log(|\epsilon_{i}|)$ in this expression is not analytically tractable, but it can be pre-computed numerically for all values of $\alpha$ in the relevant range and then approximated using a 3rd degree polynomial in $\alpha$. As shown in figure~\ref{fig:gausskl}, this approximation is extremely close.

For Gaussian dropout this then gives us:
\[
-D_{KL}[\q(w_i)|p(w_i)] \approx \text{constant} + 0.5\log(\alpha) + c_{1}\alpha + c_{2}\alpha^{2} + c_{3}\alpha^{3},
\]
with
%\[
%c_{1} = 1.161451241083230, \hspace{0.5cm} c_{2} = -1.502041176441722, \hspace{0.5cm} c_{3} = 0.586299206427007.
%\]
\[
c_{1} = 1.16145124, \hspace{0.5cm} c_{2} = -1.50204118, \hspace{0.5cm} c_{3} = 0.58629921.
\]
Alternatively, we can make use of the fact that $-\mathbb{E}_{q_{\alpha}} \log(|\epsilon_{i}|) \geq 0$ for all values of $\alpha$, and use the following lower bound:
\[
-D_{KL}[\q(w_i)|p(w_i)] \geq \text{constant} + 0.5\log(\alpha).
\]
\begin{figure}[thb]
  \begin{minipage}[c]{0.5\linewidth}
    \includegraphics[width=\linewidth]{approx_gauss_kl.pdf}
  \end{minipage}\hfill
  \begin{minipage}[c]{0.5\linewidth}
    \caption{Exact and approximate negative KL divergence $-D_{KL}[\q(w_i)|p(w_i)]$ for the \textit{log-uniform} prior and \textit{Gaussian dropout} approximate posterior. Since we use an improper prior, the constant is arbitrary; we choose it so that the exact KL divergence is 0 at $\alpha=1$, which is the maximum value we allow during optimization of the lower bound. The exact and approximate results are indistinguishable on the domain shown, which corresponds to dropout rates between $p=0.05$ and $p=0.5$. For $\log(\alpha)\rightarrow-\infty$ the approximated term vanishes, so that both our polynomial approximation and the lower bound become exact.}
\label{fig:gausskl}
  \end{minipage}
\end{figure}
\FloatBarrier

\section{Variance reduction of local parameterization compared to separate weight samples}
\label{sec:low_variance_reparam}

In section~\ref{sec:fastsgvb} we claim that our local reparameterization not only is more computationally efficient than drawing separate weight matrices for each training example, but that it also leads to an estimator that has lower variance. To see why, consider the stochastic gradient estimate with respect to the posterior parameter $\sigma^{2}_{i,j}$ for a minibatch of size $M=1$. Drawing random weights $\bW$ to form the estimate, we get
\begin{align}
\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}} = \frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\frac{\epsilon_{i,j}a_{m,i}}{2\sigma_{i,j}}.
\end{align}
If, on the other hand, we form the same gradient using the local reparameterization trick, we get
\begin{align}
\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}} = \frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\frac{\zeta_{m,j}a_{m,i}^{2}}{2\sqrt{\delta_{m,j}}}.
\end{align}
We can decompose the variance of these estimators as
\[
\Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}} = \Var{b_{m,j}}{\Exp{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}|b_{m,j}}} + \Exp{b_{m,j}}{\Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}|b_{m,j}}}.
\]
Here, the first term is identical for both gradient estimators, but the second term is not. When drawing separate weight matrices per training example, it is
\begin{align}
\Exp{b_{m,j}}{\Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}|b_{m,j}}} = \Exp{b_{m,j}}{ \Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}|b_{m,j}}q} \nonumber\\
+ \Exp{b_{m,j}}{ \Exp{\q,\D}{\left(\frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\right)^{2}|b_{m,j}}\frac{\Var{\q,\D}{\epsilon_{i,j}^{2}|b_{m,j}}a_{m,i}^{2}}{4\sigma_{i,j}^{2}} },
\end{align}
with $q = (b_{m,j}-\gamma_{m,j})^{2}a_{m,i}^{4}/4$.

When using the local reparameterization trick, we simply have
\begin{align}
\Exp{b_{m,j}}{\Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}|b_{m,j}}} = \Exp{b_{m,j}}{ \Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}|b_{m,j}}q}.
\end{align}
In this case, the second term vanishes because the random variable $\zeta_{m,j}$ is uniquely determined given $b_{m,j}$, as there is only a single $\zeta_{m,j}$ for each $b_{m,j}$. Because there are many $\epsilon_{i,j}$ for each $b_{m,j}$, this variable is not uniquely determined by conditioning on $b_{m,j}$, so there is an additional (positive) contribution to the variance.

\section{Variance of stochastic gradients for variational dropout with correlated weight noise}
\label{sec:variance_corr}

In section~\ref{sec:experiments} we empirically compare the variance of the different available stochastic estimators of the gradient of our variational objective for a model regularized using variational dropout with independent weight noise (section\ref{sec:dropout_indep}). Here we present the same analysis for variational dropout with correlated weight noise (section\ref{sec:dropout_corr}). As table~\ref{tab:mnist_var_corr} shows, the relative performance of the different estimators is similar to that reported in section~\ref{sec:experiments}, although the noise level is generally much higher using this regularization method.
\begin{table}[htb]
\begin{tabular}{lllll}
 & top layer & top layer & bottom layer & bottom layer \\
stochastic gradient estimator & 10 epochs & 100 epochs & 10 epochs & 100 epochs \\
\hline
local reparameterization (ours) & $2.9\times 10^{4}$ & $5.8\times 10^{3}$ & $8.8\times 10^{2}$ & $5.3\times 10^{2}$ \\
weight sample per minibatch (standard) & $3.1\times 10^{5}$ & $2.0\times 10^{4}$ & $5.0\times 10^{3}$ & $1.1\times 10^{3}$\\
no dropout noise (minimal var.) & $7.1\times 10^{3}$ & $9.8\times 10^{2}$ & $4.3\times 10^{2}$ & $1.2\times 10^{2}$
\end{tabular}
\caption{Average empirical variance of minibatch stochastic gradient estimates (1000 examples) for a fully connected neural network, regularized by variational dropout with correlated weight noise. For this specification, the variance of the stochastic gradient estimator with separate weight samples for each data point is identical to that using the local reparameterization trick: Ignoring optimization of dropout rates, both methods correspond to the standard dropout training procedure.}
\label{tab:mnist_var_corr}
\end{table}

\section{Variance of SGVB estimator with minibatches of datapoints without replacement}
\label{sec:variance_without_replacement}
Using the indicator variables $s_i \in S$, to denote whether the $i$-th training example is included in our minibatch, and $L_{i}(\beps,\phi)$ as shorthand for $\log p(\by^i|\bx^i,\bw = f(\beps, \bphi))$, the Monte Carlo estimator \eqref{eq:objective_sgvb} may be rewritten as
\begin{align}
L_\D^\text{SGVB}(\bphi) = \frac{N}{M} \sum_{i=1}^N s_{i}L_{i}(\beps,\phi).
\end{align}
The variance of which is given by
\begin{align}
\Var{\beps,S}{L_\D^\text{SGVB}(\bphi)} \leq & \frac{N^{2}}{M^{2}} \sum_{i=1}^N \Var{S}{s_{i}}\frac{L_\D(\bphi)^{2}}{N^{2}} + \Var{\beps}{L_{i}(\beps,\phi)}\Exp{S}{s_{i}^{2}} \nonumber\\
& + \frac{2N^{2}}{M^{2}} \sum_{i>j} \Exp{S}{s_{i}}\Exp{S}{s_{j}}\Cov{\beps}{L_{i}(\beps,\phi),L_{j}(\beps,\phi)} \nonumber\\
\leq & \frac{N}{M}\left( \sum_{i=1}^N \frac{L_\D(\bphi)^{2}}{N^{2}} + \Var{\beps}{L_{i}(\beps,\phi)} \right) + 2 \sum_{i>j} \Cov{\beps}{L_{i}(\beps,\phi),L_{j}(\beps,\phi)},
\end{align}
where the inequalities arise through ignoring the negative correlations between the $s_{i}$, and where we use that $\Exp{S}{s_{i}} = \Exp{S}{s_{i}^{2}} = M/N$, and $\Var{S}{s_{i}} \leq M/N$, for $N>2M$.

As can be seen, the variance contributed by the random selection of data points from $\D$ is inversely proportional to the minibatch size $M$. However, the variance contributed by our random sampling of $\beps$ does not decrease with $M$, as we are using a single random sample for the entire minibatch and $\Cov{\beps}{L_{i}(\beps,\phi),L_{j}(\beps,\phi)}$ is positive on average\footnote{We have $\Exp{\bx^i,\by^i}{L_{i}(\beps,\phi)} = \Exp{\bx^j,\by^j}{L_{j}(\beps,\phi)}$ if the examples in the training data are i.i.d., which means that $\Exp{\bx^i,\by^i,\bx^j,\by^j}{ \Cov{\beps}{L_{i}(\beps,\phi),L_{j}(\beps,\phi)} } = \Var{\beps}{ \Exp{\bx^i,\by^i}{L_{i}(\beps,\phi)} } \geq 0$.}. In practice, this means that the variance of $L_\D^\text{SGVB}(\bphi)$ can be dominated by $\beps$ for even moderately large $M$.
\end{document}
fsgvb.tex
\documentclass{article} % For LaTeX2e
\usepackage{nips15submit_e,times}
\usepackage{hyperref}
\usepackage{url}

\title{Variational Dropout and\\ the Local Reparameterization Trick}

\author{Diederik P. Kingma$^*$, Tim Salimans$^\times$ and Max Welling$^{*\dagger}$ \\
$^*$ Machine Learning Group, University of Amsterdam\\
$^\times$ Algoritmica\\
$^\dagger$ University of California, Irvine, and the Canadian Institute for Advanced Research (CIFAR)\\
\texttt{D.P.Kingma@uva.nl, salimans.tim@gmail.com, M.Welling@uva.nl}
}

\newcommand{\fix}{\marginpar{FIX}}
\newcommand{\new}{\marginpar{NEW}}

\usepackage[sort&compress, square, numbers]{natbib}
\usepackage{amsmath, amsthm}
\usepackage{amsfonts}
\usepackage{amssymb}

\usepackage[pdftex]{graphicx} % more modern
%\usepackage{subfigure} 
\usepackage{subcaption}
\graphicspath{{figures/}} 

%\usepackage{algpseudocode}
\usepackage{algorithm} 
\usepackage{algorithmic}
\usepackage{hyperref}
\usepackage{balance}
\usepackage{url}
%\newcommand{\theHalgorithm}{\arabic{algorithm}}

\usepackage{float}
\usepackage{placeins}


% Package for highlighting text
\usepackage{soul}

\newcommand{\ba}{\mathbf{a}}
\newcommand{\bb}{\mathbf{b}}
\newcommand{\bd}{\mathbf{d}}
\newcommand{\bw}{\mathbf{w}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\by}{\mathbf{y}}
\newcommand{\Exp}[2]{\mathbb{E}_{#1}\left[#2\right]}
\newcommand{\Var}[2]{\text{Var}_{#1}\left[#2\right]}
\newcommand{\Cov}[2]{\text{Cov}_{#1}\left[#2\right]}
\newcommand{\Corr}[2]{\text{Corr}_{#1}\left[#2\right]}

\newcommand{\bA}{\mathbf{A}}
\newcommand{\bB}{\mathbf{B}}
\newcommand{\bV}{\mathbf{V}}
\newcommand{\bW}{\mathbf{W}}
\newcommand{\bv}{\mathbf{v}}

\newcommand{\bphi}{\mathbf{\phi}}
\newcommand{\bt}{\mathbf{\theta}}
\newcommand{\bdelta}{\mathbf{\delta}}
\newcommand{\beps}{\mathbf{\epsilon}}
\newcommand{\q}{q_\mathbf{\phi}}
\newcommand{\bmu}{\mathbf{\mu}}
\newcommand{\bsigma}{\mathbf{\sigma}}
\newcommand{\bSigma}{\mathbf{\Sigma}}

\newcommand{\D}{\mathcal{D}}

\newcommand{\eqnr}{\addtocounter{equation}{1}\tag{\theequation}}

\nipsfinalcopy % Uncomment for camera-ready version

\begin{document}


\maketitle

\begin{abstract}
We investigate a local reparameterizaton technique for greatly reducing the variance of stochastic gradients for variational Bayesian inference (SGVB) of a posterior over model parameters, while retaining parallelizability. This local reparameterization translates uncertainty about global parameters into local noise that is independent across datapoints in the minibatch. Such parameterizations can be trivially parallelized and have variance that is inversely proportional to the minibatch size, generally leading to much faster convergence.
Additionally, we explore a connection with dropout: Gaussian dropout objectives correspond to SGVB with local reparameterization, a scale-invariant prior and proportionally fixed posterior variance. Our method allows inference of more flexibly parameterized posteriors; specifically, we propose \emph{variational dropout}, a generalization of Gaussian dropout where the dropout rates are learned, often leading to better models. The method is demonstrated through several experiments.
\end{abstract}

\section{Introduction}
Deep neural networks are a flexible family of models that easily scale to millions of parameters and datapoints, but are still tractable to optimize using minibatch-based stochastic gradient ascent. Due to their high flexibility, neural networks have the capacity to fit a wide diversity of nonlinear patterns in the data. This flexbility often leads to \emph{overfitting} when left unchecked: spurious patterns are found that happen to fit well to the training data, but are not predictive for new data. Various regularization techniques for controlling this overfitting are used in practice; a currently popular and empirically effective technique being \emph{dropout}~\cite{hinton2012improving}. In \cite{wang2013fast} it was shown that regular (binary) dropout has a Gaussian approximation called \emph{Gaussian dropout} with virtually identical regularization performance but much faster convergence. In section 5 of \cite{wang2013fast} it is shown that Gaussian dropout optimizes a lower bound on the marginal likelihood of the data. In this paper we show that a relationship between dropout and Bayesian inference can be extended and exploited to greatly improve the efficiency of variational Bayesian inference on the model parameters. This work has a direct interpretation as a generalization of Gaussian dropout, with the same fast convergence but now with the freedom to specify more flexibly parameterized posterior distributions.

Bayesian posterior inference over the neural network parameters is a theoretically attractive method for controlling overfitting; exact inference is computationally intractable, but efficient approximate schemes can be designed. Markov Chain Monte Carlo (MCMC) is a class of approximate inference methods with asymptotic guarantees, pioneered by~\cite{neal1995bayesian} for the application of regularizing neural networks. Later useful refinements include ~\cite{welling2011bayesian} and \cite{ahn2012bayesian}.

An alternative to MCMC is variational inference~\cite{hinton1993keeping} or the equivalent \emph{minimum description length} (MDL) framework. Modern variants of stochastic variational inference have been applied to neural networks with some succes~\cite{graves2011practical}, but have been limited by high variance in the gradients. Despite their theoretical attractiveness, Bayesian methods for inferring a posterior distribution over neural network weights have not yet been shown to outperform simpler methods such as dropout. Even a new crop of efficient variational inference algorithms based on stochastic gradients with minibatches of data~\cite{salimans2013fixedform,kingma2013auto,rezende2014stochastic} have not yet been shown to significantly improve upon simpler dropout-based regularization.

In section~\ref{sec:variational_dropout} we explore an as yet unexploited trick for improving the efficiency of stochastic gradient-based variational inference with minibatches of data, by translating uncertainty about global parameters into local noise that is independent across datapoints in the minibatch. The resulting method has an optimization speed on the same level as fast dropout~\cite{wang2013fast}, and indeed has the original Gaussian dropout method as a special case. An advantage of our method is that it allows for full Bayesian analysis of the model, and that it's significantly more flexible than standard dropout. The approach presented here is closely related to several popular methods in the literature that regularize by adding random noise; these relationships are discussed in section~\ref{sec:relatedwork}. 

\section{Efficient and Practical Bayesian Inference}
\label{sec:variational_dropout}
We consider Bayesian analysis of a dataset $\D$, containing a set of $N$ i.i.d. observations of tuples $(\bx,\by)$, where the goal is to learn a model with \emph{parameters} or \emph{weights} $\bw$ of the conditional probability $p(\by|\bx,\bw)$ (standard classification or regression)\footnote{Note that the described method is not limited to classification or regression and is straightforward to apply to other modeling settings like unsupervised models and temporal models.}. Bayesian inference in such a model consists of updating some initial belief over parameters $\bw$ in the form of a prior distribution $p(\bw)$, after observing data $\D$, into an updated belief over these parameters in the form of (an approximation to) the posterior distribution $p(\bw|\D)$. Computing the true posterior distribution through Bayes' rule $p(\bw|\D) = p(\bw)p(\D|\bw)/p(\D)$ involves computationally intractable integrals, so good approximations are necessary. In \emph{variational inference}, inference is cast as an optimization problem where we optimize the parameters $\bphi$ of some parameterized model $\q(\bw)$ such that $\q(\bw)$ is a close approximation to $p(\bw|\D)$ as measured by the Kullback-Leibler divergence $D_{KL}(\q(\bw)||p(\bw|\D))$. This divergence of our posterior $\q(\bw)$ to the true posterior is minimized in practice by maximizing the so-called variational lower bound $\mathcal{L}(\bphi)$ of the marginal likelihood of the data:
\begin{align}
\mathcal{L}(\bphi) &= - D_{KL}(\q(\bw)||p(\bw)) + L_\D(\bphi)\label{eq:bound}\\
\text{where\quad} L_\D(\bphi) &=  \sum_{(\bx, \by) \in \D} \Exp{\q(\bw)}{ \log p(\by|\bx,\bw)}
\label{eq:bound2}
\end{align}
We'll call $L_\D(\bphi)$ the \emph{expected log-likelihood}. The bound $\mathcal{L}(\bphi)$ plus $D_{KL}(\q(\bw)||p(\bw|\D))$ equals the (conditional) marginal log-likelihood $\sum_{(\bx, \by) \in \D} \log p(\by|\bx)$. Since this marginal log-likelihood is constant w.r.t. $\bphi$, maximizing the bound w.r.t. $\bphi$ will minimize $D_{KL}(\q(\bw)||p(\bw|\D))$. %, which is our goal.

\subsection{Stochastic Gradient Variational Bayes (SGVB)}\label{sec:sgvb}
Various algorithms for gradient-based optimization of the variational bound (eq.~\eqref{eq:bound}) with differentiable $q$ and $p$ exist. See section~\ref{sec:relatedwork} for an overview. A recently proposed efficient method for minibatch-based optimization with differentiable models is the \emph{stochastic gradient variational Bayes} (SGVB) method introduced in~\cite{kingma2013auto} (especially appendix F) and ~\cite{rezende2014stochastic}. The basic trick in SGVB is to parameterize the random parameters $\bw \sim \q(\bw)$ as: $\bw = f(\beps, \bphi)$ where $f(.)$ is a differentiable function and $\beps \sim p(\beps)$ is a random noise variable. In this new parameterisation, an unbiased differentiable minibatch-based Monte Carlo estimator of the expected log-likelihood can be formed:
\begin{align}
L_\D(\bphi) \simeq L_\D^\text{SGVB}(\bphi) = \frac{N}{M} \sum_{i=1}^M \log p(\by^i|\bx^i,\bw = f(\beps, \bphi)),
\label{eq:objective_sgvb}
\end{align}
where $(\bx^i,\by^i)_{i=1}^M$ is a minibatch of data with $M$ random datapoints $(\bx^i,\by^i) \sim \D$, and $\beps$ is a noise vector drawn from the noise distribution $p(\beps)$. We'll assume that the remaining term in the variational lower bound, $D_{KL}(\q(\bw)||p(\bw))$, can be computed deterministically, but otherwise it may be approximated similarly. The estimator \eqref{eq:objective_sgvb} is differentiable w.r.t. $\bphi$ and unbiased, so its gradient is also unbiased: $\nabla_\phi L_\D(\bphi) \simeq \nabla_\bphi L_\D^\text{SGVB}(\bphi)$. We can proceed with variational Bayesian inference by randomly initializing $\bphi$ and performing stochastic gradient ascent on $\mathcal{L}(\bphi)$~\eqref{eq:bound}.

\subsection{Variance of the SGVB estimator}\label{sec:variance}
The theory of stochastic approximation tells us that stochastic gradient ascent using \eqref{eq:objective_sgvb} will asymptotically converge to a local optimum for an appropriately declining step size and sufficient weight updates \cite{robbinsmonro}, but in practice the performance of stochastic gradient ascent crucially depends on the variance of the gradients. If this variance is too large, stochastic gradient descent will fail to make much progress in any reasonable amount of time. Our objective function consists of an expected log likelihood term that we approximate using Monte Carlo, and a KL divergence term $D_{KL}(\q(\bw)||p(\bw))$ that we assume can be calculated analytically and otherwise be approximated with Monte Carlo with similar reparameterization. 

Assume that we draw minibatches of datapoints with replacement; see appendix~\ref{sec:variance_without_replacement} for a similar analysis for minibatches without replacement.
Using $L_i$ as shorthand for $\log p(\by^i|\bx^i,\bw = f(\beps^i, \bphi))$, the contribution to the likelihood for the $i$-th datapoint in the minibatch, the Monte Carlo estimator \eqref{eq:objective_sgvb} may be rewritten as
$L_\D^\text{SGVB}(\bphi) = \frac{N}{M} \sum_{i=1}^M L_{i}$, whose variance is given by
\begin{align}
\Var{}{L_\D^\text{SGVB}(\bphi)}
=& \frac{N^2}{M^2} \Big( \sum_{i=1}^M \Var{}{L_{i}} + 2 \sum_{i=1}^M \sum_{j=i+1}^M \Cov{}{L_{i},L_{j}} \Big)\\
=& N^2 \Big( \frac{1}{M} \Var{}{L_{i}} + \frac{M-1}{M} \Cov{}{L_{i},L_{j}} \Big),
\label{eq:objective_sgvb_var}
\end{align}
where the variances and covariances are w.r.t.\ both the data distribution and $\beps$ distribution, i.e. $\Var{}{L_{i}} = \Var{\beps,\bx^i,\by^i}{\log p(\by^i|\bx^i,\bw = f(\beps, \bphi))}$, with $\bx^i,\by^i$ drawn from the empirical distribution defined by the training set. As can be seen from \eqref{eq:objective_sgvb_var}, the total contribution to the variance by $\Var{}{L_i}$ is inversely proportional to the minibatch size $M$. However, the total contribution by the covariances does not decrease with $M$. In practice, this means that the variance of $L_\D^\text{SGVB}(\bphi)$ can be dominated by the covariances for even moderately large $M$.

\subsection{Local Reparameterization Trick}\label{sec:fastsgvb}

We therefore propose an alternative estimator for which we have $\Cov{}{L_{i},L_{j}} = 0$, so that the variance of our stochastic gradients scales as $1/M$. We then make this new estimator computationally efficient by not sampling $\beps$ directly, but only sampling the intermediate variables $f(\beps)$ through which $\beps$ influences $L_\D^\text{SGVB}(\bphi)$. By doing so, the global uncertainty in the weights is translated into a form of local uncertainty that is independent across examples and easier to sample. We refer to such a reparameterization from global noise to local noise as the \emph{local reparameterization trick}. Whenever a source of global noise can be translated to local noise in the intermediate states of computation ($\beps \rightarrow f(\beps)$), a local reparameterization can be applied to yield a computationally and statistically efficient gradient estimator.

Such local reparameterization applies to a fairly large family of models, but is best explained through a simple example: Consider a standard fully connected neural network containing a hidden layer consisting of 1000 neurons. This layer receives an $M \times 1000$ input feature matrix $\bA$ from the layer below, which is multiplied by a $1000 \times 1000$ weight matrix $\bW$, before a nonlinearity is applied, i.e. $\bB = \bA\bW$. We then specify the posterior approximation on the weights to be a fully factorized Gaussian, i.e.\ $\q(w_{i,j}) = N(\mu_{i,j},\sigma^{2}_{i,j}) \hspace{0.1cm} \forall w_{i,j} \in \bW$, which means the weights are sampled as $w_{i,j} = \mu_{i,j} + \sigma_{i,j}\epsilon_{i,j}$, with $\epsilon_{i,j} \sim N(0,1)$. In this case we could make sure that $\Cov{}{L_{i},L_{j}} = 0$ by sampling a separate weight matrix $\bW$ for each example in the minibatch, but this is not computationally efficient: we would need to sample $M$ million random numbers for just a single layer of the neural network. Even if this could be done efficiently, the computation following this step would become much harder: Where we originally performed a simple matrix-matrix product of the form $\bB = \bA\bW$, this now turns into $M$ separate local vector-matrix products. The theoretical complexity of this computation is higher, but, more importantly, such a computation can usually not be performed in parallel using fast device-optimized BLAS (Basic Linear Algebra Subprograms). This also happens with other neural network architectures such as convolutional neural networks, where optimized libraries for convolution cannot deal with separate filter matrices per example.

Fortunately, the weights (and therefore $\beps$) only influence the expected log likelihood through the neuron activations $\bB$, which are of much lower dimension. If we can therefore sample the random activations $\bB$ directly, without sampling $\bW$ or $\beps$, we may obtain an efficient Monte Carlo estimator at a much lower cost. For a factorized Gaussian posterior on the weights, the posterior for the activations (conditional on the input $\bA$) is also factorized Gaussian:
\begin{align}
\q(w_{i,j}) &= N(\mu_{i,j},\sigma^{2}_{i,j}) \hspace{0.1cm} \forall w_{i,j} \in \bW \hspace{0.2cm} \Longrightarrow \hspace{0.2cm} \q(b_{m,j}|\bA) = N(\gamma_{m,j},\delta_{m,j}), \text{ with } \nonumber\\
\gamma_{m,j} &= \sum_{i=1}^{1000} a_{m,i}\mu_{i,j}, \hspace{0.2cm} \text{ and } \hspace{0.2cm} \delta_{m,j} = \sum_{i=1}^{1000} a^{2}_{m,i}\sigma^{2}_{i,j}.
\end{align}
Rather than sampling the Gaussian weights and then computing the resulting activations, we may thus sample the activations from their implied Gaussian distribution directly, using $b_{m,j} = \gamma_{m,j} + \sqrt{\delta_{m,j}}\zeta_{m,j}$, with $\zeta_{m,j} \sim N(0,1)$. Here, $\zeta$ is an $M \times 1000$ matrix, so we only need to sample $M$ thousand random variables instead of $M$ million: a thousand fold savings.

In addition to yielding a gradient estimator that is more \textit{computationally efficient} than drawing separate weight matrices for each training example, the local reparameterization trick also leads to an estimator that has \textit{lower variance}. To see why, consider the stochastic gradient estimate with respect to the posterior parameter $\sigma^{2}_{i,j}$ for a minibatch of size $M=1$. Drawing random weights $\bW$, we get
\begin{align}
\label{eq:stochest_weights}
\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}} = \frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\frac{\epsilon_{i,j}a_{m,i}}{2\sigma_{i,j}}.
\end{align}
If, on the other hand, we form the same gradient using the local reparameterization trick, we get
\begin{align}
\label{eq:stochest_local}
\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}} = \frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\frac{\zeta_{m,j}a_{m,i}^{2}}{2\sqrt{\delta_{m,j}}}.
\end{align}
Here, there are two stochastic terms: The first is the backpropagated gradient $\partial L_{\D}^\text{SGVB} / \partial b_{m,j}$, and the second is the sampled random noise ($\epsilon_{i,j}$ or $\zeta_{m,j}$). Estimating the gradient with respect to $\sigma^{2}_{i,j}$ then basically comes down to estimating the covariance between these two terms. This is much easier to do for $\zeta_{m,j}$ as there are much fewer of these: individually they have higher correlation with the backpropagated gradient $\partial L_\D^\text{SGVB}/\partial b_{m,j}$, so the covariance is easier to estimate. In other words, measuring the effect of $\zeta_{m,j}$ on $\partial L_\D^\text{SGVB}/\partial b_{m,j}$ is easy as $\zeta_{m,j}$ is the only random variable directly influencing this gradient via $b_{m,j}$. On the other hand, when sampling random weights, there are a thousand $\epsilon_{i,j}$ influencing each gradient term, so their individual effects get lost in the noise. In appendix \ref{sec:low_variance_reparam} we make this argument more rigorous, and in section~\ref{sec:experiments} we show that it holds experimentally.


%The problem of approximating the gradient thus boils down to estimating one of the correlation coefficients $\rho_{1} = \Corr{q,\D}{\partial L_\D^\text{SGVB}/\partial b_{m,j}, \epsilon_{i,j}}$ or $\rho_{2} = \Corr{q,\D}{\partial L_\D^\text{SGVB}/\partial b_{m,j}, \zeta_{m,j}}$. As is well known, accurately estimating a correlation coefficient is easier if the absolute correlation is large: in general the standard error for the sample correlation coefficient is given by
%\[
%\text{SE}{\rho} = \frac{1 - \rho^{2}}{\sqrt{M - 1}}.
%\]
%Since there are a thousand times fewer $\zeta_{m,j}$ than $\epsilon_{i,j}$, individually they have much higher correlation with the gradient $\partial L_\D^\text{SGVB}/\partial b_{m,j}}$, so they are much easier to estimate. In other words, measuring the effect of $\zeta_{m,j}$ on $\partial L_\D^\text{SGVB}/\partial b_{m,j}}$ is easy as $\zeta_{m,j}$ is the only random variable influencing $b_{m,j}$. On the other hand, when sampling random weights, there are a thousand $\epsilon_{i,j}$ influencing each $b_{m,j}$, so their individual effects get lost in the noise.



%Trivial parallelizability can be retained by translating weight uncertainty in $\q(\bw)$ into local uncertainty in the intermediate states of computation. We call this the \emph{local reparameterisation trick}. Let $\bb^i = f(\ba^i, \bw^i)$ be an intermediate variable of computation for datapoint $i$, with $\ba^i$ a previous state of computation (typically a function of the model input $\bx^i$) and $\bw^i$ the model parameters with distribution $\bw^i \sim \q(\bw)$ that is equal for all $i$. Let $L(\bb^i)$ be (part of) an objective, such as a log-likelihood function, whose expectation w.r.t. $\q(\bw)$ we would like to estimate (equations ~\eqref{eq:bound} and \eqref{eq:objective_sgvb}). Given the distribution $\q(\bw)$, we equivalently have some distribution over $\bb^i$, namely $\q(\bb^i|\ba^i)$. The trick is to parameterize $\bb^i \sim \q(\bb^i|\ba^i)$ as a differentiable function $\bb^i = g^\text{VD}(\beps^i, \ba^i, \bphi)$ with random noise $\beps^i \sim p^\text{VD}(\beps)$, such that:
%\begin{align}
%\Exp{\q(\bw)}{L(f(\ba^i, \bw^i))}
%&= \Exp{\q(\bb^i|\ba^i)}{L(\bb^i)}
%= \Exp{p^\text{VD}(\beps)}{L(g^\text{VD}(\beps^i,\ba^i,\bphi))}
%\end{align}
%The distribution of the original SGVB parameterisation $\bb^i = f(\ba^i, f(\beps^i, \bphi))$ with $\beps^i \sim p(\beps)$, equals the distribution of the new parameterisation $\bb^i = g^\text{VD}(\beps^i, \ba^i, \bphi)$ with $\beps^i \sim p^\text{VD}(\beps)$. The key difference is that computation of $g^\text{VD}(\beps^i, \ba^i, \bphi)$ can often be trivially parallelized, whereas $f(\ba^i, f(\beps^i, \bphi))$ cannot. The resulting operations can be combined into a \emph{variational dropout} (VD) estimator
%\begin{align}
%L_\D(\bphi) \simeq L_\D^\text{VD}(\bphi)
%\label{eq:objective_fsgvb}
%\end{align}
%that is generally much more efficient than $L_\D^\text{SGVB}(\bphi)$ with minibatches.


%\subsection{Fast SGVB through Reparameterization with Local Noise[todo: improve clarity]}\label{sec:fastsgvb}
%
%There is a \emph{de facto} standard strategy for efficiently optimizing objective functions with minibatches of data: vectorize the steps of operations such that the bulk of computation is performed in a small number of large operations, such as large matrix-matrix operations or convolutional operations. Using device-optimized BLAS (Basic Linear Algebra Subprograms) libraries, these few large operations can be automatically and efficiently parallelized across multiple CPU or GPU cores. This strategy is at the heart of efficient implementations in software like Matlab, Theano~\cite{bastien2012theano} and Torch~\cite{collobert2011torch7}, and typically leads to orders of magnitude faster computation than non-vectorized implementations.
%
%A property of $L_\D^\text{SGVB}$ is that the stochasticity is global: weights are sampled from $\q$ only once, then used across all datapoints in the minibatch. This is sub-optimal in terms of variance. A better estimator is one where the noise is independent per datapoint: the variance of the estimator is then inversely proportional to the minibatch size $M$. A first na\"ive attempt is:
%$L_\D^{\text{Na\"ive}}(\bphi) = \frac{N}{M} \sum_{i=1}^M \log p(\by^i|\bx^i,\bw^i)$
%with an explicit separate sample $\bw^i \sim \q(\bw)$ for each individual datapoint $i$ in the minibatch, using a parameterization $\bw^i = f(\beps^i, \bphi)$.
%however, we then lose the ability to easily parallelize this computation, generally leading to drastically slower computation as explained in the previous paragraph. 
%
%[Todo: equation that shows variance reduction effect of separate $\bw^i$ per element in the minibatch.]
%
%Trivial parallelizability can be retained by translating weight uncertainty in $\q(\bw)$ into local uncertainty in the intermediate states of computation. We call this the \emph{local reparameterisation trick}. Let $\bb^i = f(\ba^i, \bw^i)$ be an intermediate variable of computation for datapoint $i$, with $\ba^i$ a previous state of computation (typically a function of the model input $\bx^i$) and $\bw^i$ the model parameters with distribution $\bw^i \sim \q(\bw)$ that is equal for all $i$. Let $L(\bb^i)$ be (part of) an objective, such as a log-likelihood function, whose expectation w.r.t. $\q(\bw)$ we would like to estimate (equations ~\eqref{eq:bound} and \eqref{eq:objective_sgvb}). Given the distribution $\q(\bw)$, we equivalently have some distribution over $\bb^i$, namely $\q(\bb^i|\ba^i)$. The trick is to parameterize $\bb^i \sim \q(\bb^i|\ba^i)$ as a differentiable function $\bb^i = g^\text{local}(\beps^i, \ba^i, \bphi)$ with random noise $\beps^i \sim p^\text{local}(\beps)$, such that:
%\begin{align}
%\Exp{\q(\bw)}{L(f(\ba^i, \bw^i))}
%&= \Exp{\q(\bb^i|\ba^i)}{L(\bb^i)}
%= \Exp{p^\text{local}(\beps)}{L(g^\text{local}(\beps^i,\ba^i,\bphi))}
%\end{align}
%The distribution of the original SGVB parameterisation $\bb^i = f(\ba^i, f(\beps^i, \bphi))$ with $\beps^i \sim p(\beps)$, equals the distribution of the new parameterisation $\bb^i = g^\text{local}(\beps^i, \ba^i, \bphi)$ with $\beps^i \sim p^\text{local}(\beps)$. The key difference is that computation of $g^\text{local}(\beps^i, \ba^i, \bphi)$ can often be trivially parallelized, whereas $f(\ba^i, f(\beps^i, \bphi))$ cannot. The resulting operations can be combined into a \emph{fast SGVB} estimator
%\begin{align}
%L_\D(\bphi) \simeq L_\D^\text{F-SGVB}(\bphi)
%\label{eq:objective_fsgvb}
%\end{align}
%that is generally much more efficient than $L_\D^\text{SGVB}(\bphi)$ with minibatches.

%\subsection{Example: Dot Product, Gaussian Weight Uncertainty}
%Let us give an example of the local reparameterization trick, with the common linear operation $\bb^i = \ba^i \bW$, where the row vector $\ba^i$ is of dimension $1 \times N_A$ and is the result of earlier operations, and $\bb^i$ is of dimension $1 \times N_B$. Let $\bW$ be a parameter matrix, with approximate posterior $\q(\bW) = \prod_{jk}\q(W_{jk})$, with $\q(W_{jk}) = \mathcal{N}(\mu_{jk},\bSigma_{jk})$ where $\bmu$ and $\bSigma$ are variational parameter matrices. With $(\odot, (.)^2, \sqrt{(.)})$ we denote element-wise product, square and square-root.
%
%A non-parallelizable parameterisation of $\bb^i$, is: $\bb^i = \ba^i (\bmu+\beps^i \odot \sqrt{\bSigma})$.
%
%A trivially parallelizable reparameterisation of $\bb^i$ with local noise $\epsilon^i \sim \mathcal{N}(0,I)$ is: $\bb^i = \ba^i \bmu + \epsilon^i \odot \sqrt{(\ba^i)^2 \bSigma}$. Let $\bB = \bA \bW$, where $\bA$ is a matrix of size $M \times N_A$ with $\ba^i$ as $i$-th row, and $\bB$ is the corresponding result with size $M \times N_B$. The vectorized form is $\bB = \bA \bmu + \epsilon \odot \sqrt{\bA^2 \bSigma}$. In this form, most of the computation is contained in a few large matrix-matrix operations, whose computation is trivially parallelized.

%\section{Gaussian Dropout as a Special Case}
%\label{sec:gauss_dropout}
%In \cite{wang2013fast} it was shown that regularization by regular (binary) dropout \cite{srivastava2014dropout} has a fast approximation called \emph{Gaussian dropout}. In Gaussian dropout, linear operations $\bb^i = \ba^i \bW$ are replaced by a stochastic relation $p(\bb^i|\ba^i) = \mathcal{N}(\ba^i\bW,\alpha (\ba^i)^2 \bW^2)$ where $(.)^2$ is again a component-wise square, and $\alpha = p/(1-p)$ where $p$ is the dropout rate in original dropout. Gaussian dropout was shown to empirically perform as good as regular dropout, but generally with much faster convergence. 
%Gaussian uncertainty over $\bb^i$ can be straightforwardly reparameterized as Gaussian uncertainty over the weights $\bW$. See appendix~\ref{sec:dropout} for a derivation. 
%Gaussian dropout's implicit Gaussian weight uncertainty has an interpretation as a Gaussian variational posterior $\q(w_j) = \mathcal{N}(\theta_j, \alpha \theta_j^2)$, with fixed $\alpha=p/(1-p)$.

%As we will explain, Gaussian dropout corresponds to variational dropout with this particular Gaussian posterior and the \emph{log-uniform prior}. With the \emph{log-uniform prior}, we denote a prior for each individual parameter $w_j$: $p(\text{sign}(w_j)) = .5$ and $p(\log |w_j|) = \mathcal{U}(\text{interval})$, where the interval covers the complete log-domain of the magnitude of floating point numbers. Appendix~\ref{sec:floatingpoint} explains the relationship of the log-uniform distribution with floating-point numbers.

%Of key importance is that $D_{KL}(\q(\bw)||p(\bw))$, for this choice of prior and posterior, is just a function of $\alpha$. See appendix appendix~\ref{sec:kldiv} for proof. Since $\alpha=p/(1-p)$ is constant, $D_{KL}(\q(\bw)||p(\bw))$ is constant. As a result, all that is required for maximization of the bound (equation ~\eqref{eq:bound}), with this Gaussian posterior and a log-uniform prior is the maximization of $L_\D(\bphi)$ (equation ~\eqref{eq:bound2}); the rest of the terms are constant. The term $L_\D(\bphi)$ is exactly the term that is maximized by Gaussian dropout. 

%Thus, we have shown that optimizing a model with Gaussian dropout is a special case of variational dropout; i.e., SGVB with the local reparameterization trick, a log-uniform prior and a posterior on weights $w_j$: $\q(w_j) = \mathcal{N}(\theta_j, \alpha \theta_j^2)$, with fixed $\alpha = p/(1-p)$ where $p$ is the dropout rate. 

%\section{Examples of Variational Dropout}

%Despite its rigid implicit posterior variance (as we have shown), standard Gaussian dropout is known as an effective regularizer. For this reason, we will use standard Gaussian dropout as a starting point for exploration of different types of variational dropout.

%\subsection{Diagonal Gaussian Posterior}
%One promising direction is to keep the log-uniform prior and Gaussian posterior, and infer both a posterior mean and variance for each individual weight $w_j$, as opposed to just a mean:
%\begin{align}
%q_{\phi_j}(w_j) &= \mathcal{N}(\theta_j,\sigma_j^2)
%\text{\quad where\quad} \sigma_j^2 = \exp(\beta_j) \theta_j^2
%\end{align}
%where we parameterize $\alpha_j = \exp(\beta_j)$ such that $\phi_j = (\theta_j,\beta_j)$ are the variational parameters for weight $w_j$. As shown in appendix~\ref{sec:kldiv}, the KL-divergence between this Gaussian $\q(w_j)$ and log-uniform $p(w_j)$ is intractable, but can be extremely closely approximated by the 3rd degree polynomial:
%\begin{align}
%-D_{KL}[\q(w_j)|p(w_j)] \approx \text{constant} + 0.5\log(\alpha_j) + c_{1}\alpha_j + c_{2}\alpha_j^{2} + c_{3}\alpha_j^{3}
%\end{align}
%The coefficients $c1$, $c2$ and $c3$ are given in the appendix.

%Say something about relatation to Jeffrey's prior? Also scale-invariant, but our log-uniform prior is not improper since it has finite domain and thus normalizable.

%\subsection{Old: log-uniform posterior}
 %We parameterize the approximate posterior for each $w_j$ as $\q(\log |w_j|) = \mathcal{U}(\log |\theta_j|,\log |\theta_j|+ \exp(\delta_j))$, with variational parameters $\bphi=\{\bt,\bdelta\}$: a location parameter $\theta_j$ in regular space (including a sign), coupled with a scale of $\exp(\delta_j)$ in the logarithm of its magnitude. We constrain each $\delta_j$ to be non-positive for reasons that we will explain.
%We parameterize the scale in log-space of each $\q(w_j)$ as $\exp(\delta_j)$. 
%This results in an extraordinarily simple KL-divergence with the log-uniform prior: $D_{KL}(\q(|w_j|)||p(|w_j|)) = - \delta_j$. Note the scale-invariant property of this divergence: it is constant w.r.t. $\bt$.  When plugging this choice of $\q(\bw)$ and $p(\bw)$ in the variational bound of eq.~\eqref{eq:bound}, we arrive at our objective function:
%\begin{align}
%\mathcal{L}(\bphi) = \sum_j \delta_j + \sum_{(\bx, \by) \in \D} \Exp{\q(\bw)}{ \log p(\by|\bx,\bw)}
%\label{eq:objective}
%\end{align}
%where the scale parameter $\delta_j$ is constrained to be non-positive. We can use the parameterization $\bw_j = \theta_j \cdot \exp(\epsilon_j \cdot \delta_j)$ using uniform random draws $\epsilon_j \sim \mathcal{U}(0,1)$ for a Monte Carlo estimator of this bound.


\section{Variational Dropout}
\label{sec:gauss_dropout}
\textit{Dropout} is a technique for regularization of neural network parameters, which works by adding multiplicative noise to the input of each layer of the neural network during optimization. Using the notation of section~\ref{sec:fastsgvb}, for a fully connected neural network dropout corresponds to:
\begin{align}
\label{eq:dropout}
\bB = (\bA \circ \xi)\theta, \hspace{0.1cm} \text{ with }\hspace{0.1cm} \xi_{i,j} \sim p(\xi_{i,j})
\end{align}
where $\bA$ is the $M \times K$ matrix of input features for the current minibatch, $\theta$ is a $K \times L$ weight matrix, and $\bB$ is the $M \times L$ output matrix for the current layer (before a nonlinearity is applied). The $\circ$ symbol denotes the elementwise (Hadamard) product of the input matrix with a $M \times K$ matrix of independent noise variables $\xi$. By adding noise to the input during training, the weight parameters $\theta$ are less likely to overfit to the training data, as shown empirically by previous publications. Originally, \cite{hinton2012improving} proposed drawing the elements of $\xi$ from a Bernoulli distribution with probability $1-p$, with $p$ the \textit{dropout rate}. Later it was shown that using a continuous distribution with the same relative mean and variance, such as a Gaussian $N(1,\alpha)$ with $\alpha = p/(1-p)$, works as well or better \cite{srivastava2014dropout}.

Here, we re-interpret dropout with continuous noise as a variational method, and propose a generalization that we call \emph{variational dropout}. In developing variational dropout we provide a firm Bayesian justification for dropout training by deriving its implicit prior distribution and variational objective. This new interpretation allows us to propose several useful extensions to dropout, such as a principled way of making the normally fixed dropout rates $p$ \emph{adaptive} to the data.

\subsection{Variational dropout with independent weight noise} \label{sec:dropout_indep}
If the elements of the noise matrix $\xi$ are drawn independently from a Gaussian $N(1,\alpha)$, the marginal distributions of the activations $b_{m,j} \in \bB$ are Gaussian as well:
\begin{align}
\label{eq:marg_act}
\q(b_{m,j}|\bA) = N(\gamma_{m,j},\delta_{m,j}), \text{ with } \gamma_{m,j} = \sum_{i=1}^{K} a_{m,i}\theta_{i,j}, \text{ and } \delta_{m,j} = \alpha \sum_{i=1}^{K} a^{2}_{m,i}\theta_{i,j}^{2}.
\end{align}
Making use of this fact, \cite{wang2013fast} proposed \emph{Gaussian dropout}, a regularization method where, instead of applying \eqref{eq:dropout}, the activations are directly drawn from their (approximate or exact) marginal distributions as given by \eqref{eq:marg_act}. \cite{wang2013fast} argued that these marginal distributions are exact for Gaussian noise $\xi$, and for Bernoulli noise still approximately Gaussian because of the central limit theorem. This ignores the dependencies between the different elements of $\bB$, as present using \eqref{eq:dropout}, but \cite{wang2013fast} report good results nonetheless.  %Gaussian dropout was shown empirically to perform as well as regular dropout, but generally with faster convergence.

As noted by \cite{wang2013fast}, and explained in appendix~\ref{sec:dropout}, this Gaussian dropout noise can also be interpreted as arising from a Bayesian treatment of a neural network with weights $\bW$ that multiply the input to give $\bB = \bA\bW$, where the posterior distribution of the weights is given by a factorized Gaussian with $\q(w_{i,j}) = \mathcal{N}(\theta_{i,j}, \alpha \theta_{i,j}^2)$. From this perspective, the marginal distributions \eqref{eq:marg_act} then arise through the application of the local reparameterization trick, as introduced in section~\ref{sec:fastsgvb}. The variational objective corresponding to this interpretation is discussed in section~\ref{sec:dropout_prior}.

\subsection{Variational dropout with correlated weight noise}\label{sec:dropout_corr}
Instead of ignoring the dependencies of the activation noise, as in section~\ref{sec:dropout_indep}, we may retain the dependencies by interpreting dropout \eqref{eq:dropout} as a form of \emph{correlated weight noise}:
\begin{align}
\bB = (\bA \circ \xi)\theta, \hspace{0.1cm} \xi_{i,j} \sim N(1,\alpha) \hspace{0.2cm} \Longleftrightarrow \hspace{0.2cm} \bb^{m} = \ba^{m}\bW, \text{ with }\nonumber\\
\bW = (\bw_{1}',\bw_{2}',\ldots,\bw_{K}')', \text{ and } \bw_{i} = s_{i}\theta_{i}, \text{ with } \q(s_{i}) = N(1,\alpha),
\end{align}
where $\ba^{m}$ is a row of the input matrix and $\bb^{m}$ a row of the output. The $\bw_{i}$ are the rows of the weight matrix, each of which is constructed by multiplying a non-stochastic parameter vector $\theta_{i}$ by a stochastic scale variable $s_{i}$. The distribution on these scale variables we interpret as a Bayesian posterior distribution. The weight parameters $\theta_{i}$ (and the biases) are estimated using maximum likelihood. The original Gaussian dropout sampling procedure \eqref{eq:dropout} can then be interpreted as arising from a local reparameterization of our posterior on the weights $\bW$.

\subsection{Dropout's scale-invariant prior and variational objective} \label{sec:dropout_prior}
The posterior distributions $\q(\bW)$ proposed in sections~\ref{sec:dropout_indep} and \ref{sec:dropout_corr} have in common that they can be decomposed into a parameter vector $\theta$ that captures the mean, and a multiplicative noise term determined by parameters $\alpha$. Any posterior distribution on $\bW$ for which the noise enters this multiplicative way, we will call a \emph{dropout posterior}. Note that many common distributions, such as univariate Gaussians (with nonzero mean), can be reparameterized to meet this requirement.

During dropout training, $\theta$ is adapted to maximize the expected log likelihood $\Exp{q_\alpha}{L_\D(\theta)}$. For this to be consistent with the optimization of a variational lower bound of the form in \eqref{eq:bound2}, the prior on the weights $p(\bw)$ has to be such that $D_{KL}(\q(\bw)||p(\bw))$ does not depend on $\theta$. In appendix~\ref{sec:kldiv} we show that the only prior that meets this requirement is the scale invariant log-uniform prior:
\[
p(\log(|w_{i,j}|)) \propto c,
\]
i.e. a prior that is uniform on the log-scale of the weights (or the weight-scales $s_{i}$ for section~\ref{sec:dropout_corr}). As explained in appendix~\ref{sec:floatingpoint}, this prior has an interesting connection with the floating point format for storing numbers: From an MDL perspective, the floating point format is optimal for communicating numbers drawn from this prior. Conversely, the KL divergence $D_{KL}(\q(\bw)||p(\bw))$ with this prior has a natural interpretation as regularizing the number of significant digits our posterior $\q$ stores for the weights $w_{i,j}$ in the floating-point format.% In section~\ref{sec:experiments} we examine the empirical performance of this prior.

Putting the expected log likelihood and KL-divergence penalty together, we see that dropout training maximizes the following variatonal lower bound w.r.t.\ $\theta$:
\begin{align}
\label{eq:dropout_bound}
\Exp{q_{\alpha}}{L_\D(\theta)} - D_{KL}(q_{\alpha}(\bw)||p(\bw)),  
\end{align}
where we have made the dependence on the $\theta$ and $\alpha$ parameters explicit. The noise parameters $\alpha$ (e.g.\ the dropout rates) are commonly treated as hyperparameters that are kept fixed during training. For the log-uniform prior this then corresponds to a fixed limit on the number of significant digits we can learn for each of the weights $w_{i,j}$. In section~\ref{sec:adaptive_dropout} we discuss the possibility of making this limit adaptive by also maximizing the lower bound with respect to $\alpha$.

For the choice of a factorized Gaussian approximate posterior with $\q(w_{i,j}) = \mathcal{N}(\theta_{i,j}, \alpha \theta_{i,j}^2)$, as discussed in section~\ref{sec:dropout_indep}, the lower bound~\eqref{eq:dropout_bound} is analyzed in detail in appendix~\ref{sec:kldiv}. There, it is shown that for this particular choice of posterior the negative KL-divergence $-D_{KL}(q_{\alpha}(\bw)||p(\bw))$ is not analytically tractable, but can be approximated extremely accurately using
\[
-D_{KL}[\q(w_i)|p(w_i)] \approx \text{constant} + 0.5\log(\alpha) + c_{1}\alpha + c_{2}\alpha^{2} + c_{3}\alpha^{3},
\]
with
\[
c_{1} = 1.16145124, \hspace{0.5cm} c_{2} = -1.50204118, \hspace{0.5cm} c_{3} = 0.58629921.
\]
The same expression may be used to calculate the corresponding term $-D_{KL}(q_{\alpha}(\mathbf{s})||p(\mathbf{s}))$ for the posterior approximation of section~\ref{sec:dropout_corr}.

\subsection{Adaptive regularization through optimizing the dropout rate}\label{sec:adaptive_dropout}
The noise parameters $\alpha$ used in dropout training (e.g.\ the dropout rates) are usually treated as fixed hyperparameters, but now that we have derived dropout's variational objective \eqref{eq:dropout_bound}, making these parameters adaptive is trivial: simply maximize the variational lower bound with respect to $\alpha$. We can use this to learn a separate dropout rate per layer, per neuron, of even per separate weight. In section~\ref{sec:experiments} we look at the predictive performance obtained by making $\alpha$ adaptive.

We found that very large values of $\alpha$ correspond to local optima from which it is hard to escape due to large-variance gradients. To avoid such local optima, we found it beneficial to set a constraint $\alpha \leq 1$ during training, i.e.\ we maximize the posterior variance at the square of the posterior mean, which corresponds to a dropout rate of $0.5$.


%\subsection{Gaussian Dropout}
%In \cite{wang2013fast} it was shown that regularization by regular (binary) dropout \cite{srivastava2014dropout} has a fast approximation called \emph{Gaussian dropout}. In Gaussian dropout, linear operations $\bb^i = \ba^i \bW$ are replaced by a stochastic relation $p(\bb^i|\ba^i) = \mathcal{N}(\ba^i\bW,\alpha (\ba^i)^2 \bW^2)$ where $(.)^2$ is again a component-wise square, and $\alpha = p/(1-p)$ where $p$ is the dropout rate in original dropout. Gaussian dropout was shown to empirically perform as good as regular dropout, but generally with much faster convergence. 
%Gaussian uncertainty over $\bb^i$ can be straightforwardly reparameterized as Gaussian uncertainty over the weights $\bW$. See appendix~\ref{sec:dropout} for a derivation. 
%Gaussian dropout's implicit Gaussian weight uncertainty has an interpretation as a Gaussian variational posterior $\q(w_j) = \mathcal{N}(\theta_j, \alpha \theta_j^2)$, with fixed $\alpha=p/(1-p)$.

%As we will explain, Gaussian dropout corresponds to variational dropout with this particular Gaussian posterior and the \emph{log-uniform prior}. With the \emph{log-uniform prior}, we denote a prior for each individual parameter $w_j$: $p(\text{sign}(w_j)) = .5$ and $p(\log |w_j|) = \mathcal{U}(\text{interval})$, where the interval covers the complete log-domain of the magnitude of floating point numbers. Appendix~\ref{sec:floatingpoint} explains the relationship of the log-uniform distribution with floating-point numbers.

%Of key importance is that $D_{KL}(\q(\bw)||p(\bw))$, for this choice of prior and posterior, is just a function of $\alpha$. See appendix appendix~\ref{sec:kldiv} for proof. Since $\alpha=p/(1-p)$ is constant, $D_{KL}(\q(\bw)||p(\bw))$ is constant. As a result, all that is required for maximization of the bound (equation ~\eqref{eq:bound}), with this Gaussian posterior and a log-uniform prior is the maximization of $L_\D(\bphi)$ (equation ~\eqref{eq:bound2}); the rest of the terms are constant. The term $L_\D(\bphi)$ is exactly the term that is maximized by Gaussian dropout. 

%Thus, we have shown that optimizing a model with Gaussian dropout is a special case of variational dropout; i.e., SGVB with the local reparameterization trick, a log-uniform prior and a posterior on weights $w_j$: $\q(w_j) = \mathcal{N}(\theta_j, \alpha \theta_j^2)$, with fixed $\alpha = p/(1-p)$ where $p$ is the dropout rate. 

%\subsection{Examples of Variational Dropout}

%Despite its rigid implicit posterior variance (as we have shown), standard Gaussian dropout is known as an effective regularizer. For this reason, we will use standard Gaussian dropout as a starting point for exploration of various types of variational dropout.

%\subsection{Diagonal Gaussian Posterior}
%One promising direction is to keep the log-uniform prior and Gaussian posterior, and infer both a posterior mean and variance for each individual weight $w_j$, as opposed to just a mean:
%\begin{align}
%q_{\phi_j}(w_j) &= \mathcal{N}(\theta_j,\sigma_j^2)
%\text{\quad where\quad} \sigma_j^2 = \exp(\beta_j) \theta_j^2
%\end{align}
%where we parameterize $\alpha_j = \exp(\beta_j)$ such that $\phi_j = (\theta_j,\beta_j)$ are the variational parameters for weight $w_j$. As shown in appendix~\ref{sec:kldiv}, the KL-divergence between this Gaussian $\q(w_j)$ and log-uniform $p(w_j)$ is intractable, but can be extremely closely approximated by the 3rd degree polynomial:
%\begin{align}
%-D_{KL}[\q(w_j)|p(w_j)] \approx \text{constant} + 0.5\log(\alpha_j) + c_{1}\alpha_j + c_{2}\alpha_j^{2} + c_{3}\alpha_j^{3}
%\end{align}
%The coefficients $c_1$, $c_2$ and $c_3$ are given in the appendix.

%Say something about relatation to Jeffrey's prior? Also scale-invariant, but our log-uniform prior is not improper since it has finite domain and thus normalizable.

%\subsection{Old: log-uniform posterior}
 %We parameterize the approximate posterior for each $w_j$ as $\q(\log |w_j|) = \mathcal{U}(\log |\theta_j|,\log |\theta_j|+ \exp(\delta_j))$, with variational parameters $\bphi=\{\bt,\bdelta\}$: a location parameter $\theta_j$ in regular space (including a sign), coupled with a scale of $\exp(\delta_j)$ in the logarithm of its magnitude. We constrain each $\delta_j$ to be non-positive for reasons that we will explain.
%We parameterize the scale in log-space of each $\q(w_j)$ as $\exp(\delta_j)$. 
%This results in an extraordinarily simple KL-divergence with the log-uniform prior: $D_{KL}(\q(|w_j|)||p(|w_j|)) = - \delta_j$. Note the scale-invariant property of this divergence: it is constant w.r.t. $\bt$.  When plugging this choice of $\q(\bw)$ and $p(\bw)$ in the variational bound of eq.~\eqref{eq:bound}, we arrive at our objective function:
%\begin{align}
%\mathcal{L}(\bphi) = \sum_j \delta_j + \sum_{(\bx, \by) \in \D} \Exp{\q(\bw)}{ \log p(\by|\bx,\bw)}
%\label{eq:objective}
%\end{align}
%where the scale parameter $\delta_j$ is constrained to be non-positive. We can use the parameterization $\bw_j = \theta_j \cdot \exp(\epsilon_j \cdot \delta_j)$ using uniform random draws $\epsilon_j \sim \mathcal{U}(0,1)$ for a Monte Carlo estimator of this bound.


\section{Related Work}\label{sec:relatedwork}
Pioneering work in practical variational inference for neural networks was done in~\cite{graves2011practical}, where a (biased) variational lower bound estimator was introduced with good results on recurrent neural network models. In later work~\cite{kingma2013auto,rezende2014stochastic} it was shown that even more practical estimators can be formed for most types of continuous latent variables or parameters using a (non-local) reparameterization trick, leading to efficient and unbiased stochastic gradient-based variational inference. These works focused on an application to latent-variable inference; extensive empirical results on inference of global model parameters were reported in ~\cite{blundell2015weight}, including succesful application to reinforcement learning. These earlier works used the relatively high-variance estimator~\eqref{eq:objective_sgvb}, upon which we improve. Variable reparameterizations have a long history in the statistics literature, but have only recently found use for efficient \emph{gradient-based} machine learning and inference~\cite{bengio2013estimating,kingma2013fast,salimans2013fixedform}. Related is also \emph{probabilistic backpropagation}~\cite{hernandez2015probabilistic}, an algorithm for inferring marginal posterior probabilities; however, it requires certain tractabilities in the network making it insuitable for the type of models under consideration in this paper.

As we show here, regularization by dropout~\cite{srivastava2014dropout, wang2013fast} can be interpreted as variational inference. DropConnect~\cite{wan2013regularization} is similar to dropout, but with binary noise on the weights rather than hidden units. DropConnect thus has a similar interpretation as variational inference, with a uniform prior over the weights, and a mixture of two Dirac peaks as posterior. In ~\cite{ba2013adaptive}, \emph{standout} was introduced, a variation of dropout where a binary belief network is learned for producing dropout rates. Recently, \cite{maeda2014bayesian} proposed another Bayesian perspective on dropout. In recent work~\cite{bayer2015fast}, a similar reparameterization is described and used for variational inference; their focus is on closed-form approximations of the variational bound, rather than unbiased Monte Carlo estimators. ~\cite{maeda2014bayesian} and \cite{gal2015dropout} also investigate a Bayesian perspective on dropout, but focus on the binary variant. ~\cite{gal2015dropout} reports various encouraging results on the utility of dropout's implied prediction uncertainty.

%We show how SGVB can be used to infer low-precision floating-point parameters; a potential application is low-precision hardware\cite{courbariaux2014low}.

%"Building Blocks for Variational Bayesian Learning of Latent Variable Models " (Tapani et al)”

\section{Experiments}\label{sec:experiments}
We compare our method to standard binary dropout and two popular versions of Gaussian dropout, which we'll denote with type A and type B. With Gaussian dropout type A we denote the pre-linear Gaussian dropout from ~\cite{srivastava2014dropout}; type B denotes the post-linear Gaussian dropout from ~\cite{wang2013fast}. This way, the method names correspond to the matrix names in section 2 ($\bA$ or $\bB$) where noise is injected. Models were implemented in Theano~\cite{bergstra2010theano}, and optimization was performed using Adam~\cite{kingma2015adam} with default hyper-parameters and temporal averaging.

Two types of variational dropout were included. Type A is correlated weight noise as introduced in section~\ref{sec:dropout_corr}: an adaptive version of Gaussian dropout type A. Variational dropout type B has independent weight uncertainty as introduced in section~\ref{sec:dropout_indep}, and corresponds to Gaussian dropout type B.

A \emph{de facto} standard benchmark for regularization methods is the task of MNIST hand-written digit classification. We choose the same architecture as~\cite{srivastava2014dropout}: a fully connected neural network with 3 hidden layers and rectified linear units (ReLUs). We follow the dropout hyper-parameter recommendations from these earlier publications, which is a dropout rate of $p=0.5$ for the hidden layers and $p=0.2$ for the input layer. We used early stopping with all methods, where the amount of epochs to run was determined based on performance on a validation set. 

\paragraph{Variance.} We start out by empirically comparing the variance of the different available stochastic estimators of the gradient of our variational objective. To do this we train the neural network described above for either 10 epochs (test error 3\%) or 100 epochs (test error 1.3\%), using variational dropout with independent weight noise. After training, we calculate the gradients for the weights of the top and bottom level of our network on the full training set, and compare against the gradient estimates per batch of $M=1000$ training examples. Appendix~\ref{sec:variance_corr} contains the same analysis for the case of variational dropout with correlated weight noise.

Table \ref{tab:mnist_var_indep} shows that the local reparameterization trick yields the lowest variance among all variational dropout estimators for all conditions, although it is still substantially higher compared to not having any dropout regularization. The $1/M$ variance scaling achieved by our estimator is especially important early on in the optimization when it makes the largest difference (compare \textit{weight sample per minibatch} and \textit{weight sample per data point}). The additional variance reduction obtained by our estimator through drawing fewer random numbers (section~\ref{sec:fastsgvb}) is about a factor of 2, and this remains relatively stable as training progresses (compare \textit{local reparameterization} and \textit{weight sample per data point}).
\begin{table}[htb]
\begin{tabular}{lllll}
 & top layer & top layer & bottom layer & bottom layer \\
stochastic gradient estimator & 10 epochs & 100 epochs & 10 epochs & 100 epochs \\
\hline
local reparameterization (ours) & $7.8\times 10^{3}$ & $1.2\times 10^{3}$ & $1.9\times 10^{2}$ & $1.1\times 10^{2}$ \\
weight sample per data point (slow) & $1.4\times 10^{4}$ & $2.6\times 10^{3}$ & $4.3\times 10^{2}$ & $2.5\times 10^{2}$ \\
weight sample per minibatch (standard) & $4.9\times 10^{4}$ & $4.3\times 10^{3}$ & $8.5\times 10^{2}$ & $3.3\times 10^{2}$ \\
no dropout noise (minimal var.) & $2.8\times 10^{3}$ & $5.9\times 10^{1}$ & $1.3\times 10^{2}$ & $9.0\times 10^{0}$
\end{tabular}
\caption{Average empirical variance of minibatch stochastic gradient estimates (1000 examples) for a fully connected neural network, regularized by variational dropout with independent weight noise.}
\label{tab:mnist_var_indep}
\end{table}
\FloatBarrier
%\begin{table}[htb]
%\begin{tabular}{lllll}
% & top layer & top layer & bottom layer & bottom layer \\
%stochastic gradient estimator & 10 epochs & 100 epochs & 10 epochs & 100 epochs \\
%\hline
%local reparameterization (ours) & $3.2\times 10^{4}$ & $5.3\times 10^{3}$ & $8.2\times 10^{2}$ & $5.2\times 10^{2}$ \\
%separate weight samples (slow) & $6.6\times 10^{4}$ & $1.2\times 10^{4}$ & $2.1\times 10^{3}$ & $1.2\times 10^{3}$ \\
%single weight sample (standard) & $1.1\times 10^{5}$ & $1.3\times 10^{4}$ & $2.8\times 10^{3}$ & $1.3\times 10^{3}$\\
%no dropout noise (minimal var.) & $1.3\times 10^{4}$ & $3.0\times 10^{2}$ & $6.4\times 10^{2}$ & $4.5\times 10^{1}$
%\end{tabular}
%\caption{Average empirical variance of minibatch stochastic gradient estimates (200 examples) for a fully connected neural network, regularized by variational dropout with independent weight noise.}
%\label{tab:mnist_var_indep}
%\end{table}
%\begin{table}[htb]
%\begin{tabular}{lllll}
% & top layer & top layer & bottom layer & bottom layer \\
%stochastic gradient estimator & 10 epochs & 100 epochs & 10 epochs & 100 epochs \\
%\hline
%local reparameterization (ours) & $1.3\times 10^{5}$ & $2.9\times 10^{4}$ & $4.0\times 10^{3}$ & $2.7\times 10^{3}$ \\
%single weight sample (standard) & $4.4\times 10^{5}$ & $4.8\times 10^{4}$ & $9.3\times 10^{3}$ & $3.4\times 10^{3}$ \\
%no dropout noise (minimal var.) & $3.5\times 10^{5}$ & $4.7\times 10^{3}$ & $2.1\times 10^{3}$ & $5.9\times 10^{2}$
%\end{tabular}
%\caption{Average empirical variance of minibatch stochastic gradient estimates (200 examples) for a fully connected neural network, regularized by variational dropout with correlated weight noise. For this specification, the variance of the stochastic gradient estimator with separate weight samples is identical to that using the local reparameterization trick: Ignoring optimization of dropout rates, both methods correspond to the standard dropout training procedure.}
%\label{tab:mnist_var_corr}
%\end{table}

\paragraph{Speed.} We compared the regular SGVB estimator, with separate weight samples per datapoint with the efficient estimator based on local reparameterization, in terms of wall-clock time efficiency. With our implementation on a modern GPU, optimization with the na\"ive estimator took \textbf{1635 seconds} per epoch, while the efficient estimator took \textbf{7.4 seconds}: an over 200 fold speedup.

\paragraph{Classification error.} Figure~\ref{fig:mnist_classification} shows test-set classification error for the tested regularization methods, for various choices of number of hidden units. Our adaptive variational versions of Gaussian dropout perform equal or better than their non-adaptive counterparts and standard dropout under all tested conditions. The difference is especially noticable for the smaller networks. In these smaller networks, we observe that variational dropout infers dropout rates that are on average far lower than the dropout rates for larger networks. This adaptivity comes at negligable computational cost.

We found that slightly downscaling the KL divergence part of the variational objective can be beneficial. \emph{Variational (A2)} in figure~\ref{fig:mnist_classification} denotes performance of type A variational dropout but with a KL-divergence downscaled with a factor of 3; this small modification seems to prevent underfitting, and beats all other dropout methods in the tested models.



\begin{figure}[t]
  \vspace{-1cm}
\centering
\begin{subfigure}[t]{0.49\textwidth}\centering
\includegraphics[width=7cm]{graphs_classification/mnist}
\caption{Classification error on the MNIST dataset}
\end{subfigure}
\begin{subfigure}[t]{0.49\textwidth}\centering
\includegraphics[width=7cm]{graphs_classification/cifar10}
\caption{Classification error on the CIFAR-10 dataset}
\end{subfigure}
\caption{Best viewed in color.  \textbf{(a)} Comparison of various dropout methods, when applied to fully-connected neural networks for classification on the MNIST dataset. Shown is the classification error of networks with 3 hidden layers, averaged over 5 runs. he variational versions of Gaussian dropout perform equal or better than their non-adaptive counterparts; the difference is especially large with smaller models, where regular dropout often results in severe underfitting. \textbf{(b)} Comparison of dropout methods when applied to convolutional net a trained on the CIFAR-10 dataset, for different settings of network size $k$. The network has two convolutional layers with each $32 k$ and $64 k$ feature maps, respectively, each with stride 2 and followed by a softplus nonlinearity. This is followed by two fully connected layers with each $128 k$ hidden units.} \label{fig:mnist_classification}
\end{figure}


\FloatBarrier

\section{Conclusion}
Efficiency of posterior inference using stochastic gradient-based variational Bayes (SGVB) can often be significantly improved through a \emph{local reparameterization} where global parameter uncertainty is translated into local uncertainty per datapoint. By injecting noise locally, instead of globally at the model parameters, we obtain an efficient estimator that has low computational complexity, can be trivially parallelized and has low variance. We show how dropout is a special case of SGVB with local reparameterization, and suggest \emph{variational dropout}, a straightforward extension of regular dropout where optimal dropout rates are inferred from the data, rather than fixed in advance. We report encouraging empirical results.

\subsubsection*{Acknowledgments}
We thank the reviewers and Yarin Gal for valuable feedback. Diederik Kingma is supported by the Google European Fellowship in Deep Learning, Max Welling is supported by research grants from Google and Facebook, and the NWO project in Natural AI (NAI.14.108).

\small
\bibliographystyle{apalike}
\bibliography{bib}

\appendix

\newpage
\section{Floating-point numbers and compression}\label{sec:floatingpoint}
The \emph{floating-point} format is by far the most common number format used for model parameters in neural network type models. The absolute value of floating-point numbers in base 2 equals $\text{\emph{sign}} \cdot (s / 2^{p-1}) \cdot 2^e$ where $s$ is the mantissa with $p$ bits and $e$ is the exponent, both non-negative integer-valued. For random bit patterns for $s$ and $e$, floating point numbers closely follow a log-uniform distribution with finite range. The floating-point format is thus close to optimal, from a compression point of view, for communicating parameters for which the receiver has a log-uniform prior distribution, with finite range, over the magnitude of parameters. Specifically, this prior for each individual floating-point parameter $w_j$ is $p(\log |w_j|) = \mathcal{U}(\text{interval})$, where the interval is approximately $(\log(2^{-8}),\log(2^8))$ for single-point precision and $(\log(2^{-11}),\log(2^{11}))$ for double-precision floating-point numbers in the most common IEEE specification. 

The KL-divergence between a variational posterior $\q$ (eq.~\eqref{eq:bound}) and this log-uniform prior is relatively simple: it equals the entropy under $\q$ of the log-transformed magnitude of $\bw$ plus a constant:
\begin{align}
- D_{KL}(\q(\bw)||p(\bw)) = \mathcal{H}(\q(\log |\bw|)) + \text{constant}
\end{align}
This divergence has a natural interpretation as controlling the number of significant digits under $\q$ of the corresponding floating-point number.

\section{Derivation of dropout's implicit variational posterior}\label{sec:dropout}
In the common version of dropout~\cite{srivastava2014dropout}, linear operations $\bb^i = \ba^i \bW$ in a model, where $\ba^i$ is a column vector that is the result of previous operations on the input for datapoint $i$, are replaced by a stochastic operation:
\begin{align}
\bb^i &= (\ba^i \odot (\bd^i/(1-p)) \bW \\
\text{i.e.:\quad} b_k^i &= \sum_j W_{jk} a_j^i d_j^i /(1-p)\\
d_j^i &\sim \text{Bernoulli}(1-p)\\
p &= \text{dropout rate (hyper-parameter, e.g. 0.5)}
\end{align}
where column vector $\ba^i$ is the result of previous operations on the input for datapoint $i$. Expected values and variance w.r.t. the dropout rate are:
\begin{align}
\Exp{}{d_j^i} &= (1-p)\\
\Exp{}{d_j^i/(1-p)} &= 1\\
\Var{}{d_j^i/(1-p)} &= \Var{}{d_j^i}/(1-p)^2 = p/(1-p)\\
\end{align}
The expected value and variance for the elements of the resulting vector $\bb^i$ are:
\begin{align}
\Exp{}{b_k^i} &= \Exp{}{\sum_j W_{jk} a_j^i d_j^i / (1-p)}\\
&= \sum_j W_{jk} a_j^i \Exp{}{d_j^i / (1-p)}\\
&= \sum_j W_{jk} a_j^i\\
\Var{}{b_k^i} &= \Var{}{\sum_j W_{jk} a_j^i d_j^i / (1-p)}\\
&= \sum_j \Var{}{W_{jk} a_j^i d_j^i / (1-p)}\\
&= \sum_j W_{jk}^2 (a_j^i)^2 \Var{}{d_j^i / (1-p)}\\
&= p/(1-p) \sum_j W_{jk}^2 (a_j^i)^2
\end{align}
\subsection{Gaussian dropout}
Thus, ignoring dependencies between elements of $\bb^i$, dropout can be approximated by a Gaussian ~\cite{wang2013fast}:
\begin{align}
p(\bb^i|\ba^i,W) &= \mathcal{N}(\ba^i W,  p(1-p)^{-1} (\ba^i)^2 (\bW)^2)\\
\text{i.e.:\quad} \bb^i &= \ba^i W + \beps^i \odot \sqrt{p(1-p)^{-1} (\ba^i)^2 (\bW)^2}\\
\text{where:\quad} \beps &\sim \mathcal{N}(0,I)
\label{eq:gaussiandropout}
\end{align}
where $(.)^2$, $\sqrt(.)$ and $\odot$ are again component-wise operations.
\subsection{Weight uncertainty}
Equivalently, the Gaussian dropout relationship $p(\bb^i|\ba^i,W)$ can be parameterized as weight uncertainty:
\begin{align}
\bb^i &= \ba^i \bV\\
p(V_{jk}|W_{jk},\alpha_{jk}) &= \mathcal{N}(W_{jk}, \alpha_{jk} W_{jk}^2)\\
\alpha_{jk} &= p/(1-p)\\
\Rightarrow\text{\quad} p(\bb^i|\ba^i,W) &= \mathcal{N}(\ba^i W,  p(1-p)^{-1} (\ba^i)^2 (\bW)^2)
\end{align}
In the main text, we show how we can treat the distribution over weights $p(V_{jk}|W_{jk},\alpha_{jk})$ as a variational posterior $\q$, which we can optimize w.r.t. each $\alpha_{jk}$ (the dropout rate for individual parameters) to optimize a bound on the marginal likelihood of the data.

\section{Negative KL-divergence for the log-uniform prior}\label{sec:kldiv}
The variational lower bound \eqref{eq:bound} that we wish to optimize is given by
\[
\mathcal{L}(\bphi) = - D_{KL}(\q(\bw)||p(\bw)) + L_\D(\bphi),
\]
where $L_\D(\bphi)$ is the expected log likelihood term which we approximate using Monte Carlo, as discussed in section~\ref{sec:fastsgvb}. Here we discuss how to deterministically compute the negative KL divergence term $-D_{KL}(\q(\bw)||p(\bw))$ for the \textit{log-uniform} prior proposed in section~\ref{sec:dropout}. For simplicity we only discuss the case of factorizing approximate posterior $\q(\bw) = \prod_{i=1}^{k} \q(w_{i})$, although the extension to more involved posterior approximations is straight forward.

Our log-uniform prior $p(w_{i})$ (appendix \ref{sec:floatingpoint}) consists of a Bernoulli distribution on the \textit{sign bit} $s_{i} \in \{-1,+1\}$ which captures the sign of $w_{i}$,  combined with a uniform prior on the log-scale of $w_{i}$:
\begin{eqnarray}
w_{i} & = & s_{i}|w_{i}| \nonumber\\
p(s_{i}) & = & \text{Bernoulli}(0.5) \text{ on } \{-1,+1\}\nonumber\\
p(\log(|w_{i}|)) & \propto & c
\end{eqnarray}
For an approximate posterior $\q(w_{i})$ with support on the real line, the KL divergence will thus consist of two parts: the KL-divergence between the posterior and prior on the sign bit $s_{i}$ and the KL-divergence between the conditional posterior and prior on the log scale $\log(|w_{i}|)$.

The negative KL divergence between the posterior and prior on the sign bit is given by:
\begin{align}
-D_{KL}[\q(s_{i})|p(s_{i})] = \log(0.5) - Q_{\phi}(0)\log[Q_{\phi}(0)] - [1-Q_{\phi}(0)]\log[1-Q_{\phi}(0)],
\end{align}
with $Q_{\phi}(0)$ the posterior CDF of $w_{i}$ evaluated at 0, i.e. the probability of drawing a negative weight from our approximate posterior.

The negative KL divergence of the second part (the log-scale) is then
\begin{align}
-D_{KL}[\q(\log(|w_{i}|) | s_{i})|p(\log(|w_{i}|))] = \log(c) - Q_{\phi}(0)\mathbb{E}_{\q(w_{i}|w_{i}<0)}\left[\log(\q(w_{i})/Q_{\phi}(0)) + \log(|w_{i}|)\right] \nonumber \\
 - \left(1-Q_{\phi}(0)\right)\mathbb{E}_{\q(w_{i}|w_{i}>0)}\left[\log(\q(w_{i})/(1-Q_{\phi}(0))) + \log(|w_{i}|)\right], \nonumber
\end{align}
where $\q(w_{i})/Q_{\phi}(0)$ is the truncation of our approximate posterior to the negative part of the real-line, and the $\log(|w_{i}|)$ term enters through transforming the uniform prior from the log-space to the normal space.

Putting both parts together the terms involving the CDF $Q_{\phi}(0)$ cancel, and we get:
\begin{align}
\label{eq:combined_kl}
-D_{KL}[\q(w_i)|p(w_i)] = \log(c) + \log(0.5) + \mathcal{H}(\q(w_{i})) - \mathbb{E}_{q} \log(|w_{i}|),
\end{align}
where $\mathcal{H}(\q(w_{i}))$ is the entropy of our approximate posterior.

We now consider the special case of a \textit{dropout posterior}, as proposed in section~\ref{sec:gauss_dropout}, which we defined as any approximate posterior distribution $\q(w_{i})$ for which the weight noise is \textit{multiplicative}:
\begin{align}
w_{i} = \theta_{i}\epsilon_{i}, \text{ with } \epsilon_{i} \sim q_{\alpha}(\epsilon_{i}),
\end{align}
where we have divided the parameters $\phi_{i} = (\theta_{i}, \alpha)$ into parameters governing the mean of $\q(w_i)$ (assuming $\mathbb{E}_{q} \epsilon_{i} = 1$) and the multiplicative noise. Note that certain families of distributions, such as Gaussians, can always be reparameterized in this way.

Assuming $q_{\alpha}(\epsilon_{i})$ is a continuous distribution, this choice of posterior gives us
\begin{align}
\label{eq:h_eps}
\mathcal{H}(\q(w_{i})) = \mathcal{H}(q_{\alpha}(\epsilon_{i})) + \log(|\theta_{i}|),
\end{align}
where the last term can be understood as the log jacobian determinant of a linear transformation of the random variable $\epsilon_{i}$.

Furthermore, we can rewrite the $-\mathbb{E}_{q} \log(|w_{i}|)$ term of \eqref{eq:combined_kl} as
\begin{align}
\label{eq:e_w}
-\mathbb{E}_{q} \log(|w_{i}|) =  -\mathbb{E}_{q} \log(|\theta_{i}\epsilon_{i}|) = -\log(|\theta_{i}|) - \mathbb{E}_{q} \log(|\epsilon_{i}|).
\end{align}
Taking \eqref{eq:h_eps} and \eqref{eq:e_w} together, the terms depending on $\theta$ thus cancel exactly, proving that the KL divergence between posterior and prior is indeed independent of these parameters:
\begin{align}
-D_{KL}[\q(w_i)|p(w_i)] & = \log(c) + \log(0.5) + \mathcal{H}(\q(w_{i})) - \mathbb{E}_{q} \log(|w_{i}|) \nonumber\\
& = \log(c) + \log(0.5) + \mathcal{H}(q_{\alpha}(\epsilon_{i})) - \mathbb{E}_{q_{\alpha}} \log(|\epsilon_{i}|).
\end{align}
In fact, the log-uniform prior is the only prior consistent with dropout. If the prior had any additional terms involving $w_{i}$, their expectation w.r.t.\ $\q$ would not cancel with the entropy term in $-D_{KL}[\q(w_i)|p(w_i)]$.

For the specific case of \textit{Gaussian dropout}, as proposed in section~\ref{sec:gauss_dropout}, we have $q_{\alpha}(\epsilon_{i}) = N(1,\alpha)$. This then gives us
\[
-D_{KL}[\q(w_i)|p(w_i)] = \log(c) + \log(0.5) + 0.5(1 + \log(2\pi) + \log(\alpha)) - \mathbb{E}_{q_{\alpha}} \log(|\epsilon_{i}|).
\]
%For this choice of posterior we have $\mathcal{H}(\q(w_{i})) = 0.5(1 + \log(2\pi) + \log(\alpha) + \log(\theta_{i}^{2}))$, and we can rewrite the last term in the KL divergence $-\mathbb{E}_{q} \log(|w_{i}|)$ using the \textit{reparameterization trick} as
%\begin{align}
%-\mathbb{E}_{q} \log(|w_{i}|) & = & -0.5\mathbb{E}_{q} \log(w_{i}^{2}) = -0.5\mathbb{E}_{\epsilon \sim \mathcal{N}(0,1)}\log[(\theta_{i}+\alpha\theta_{i}\epsilon)^{2}] \\
%& = & -0.5\log(\theta_{i}^{2}) - 0.5\mathbb{E}_{\epsilon \sim \mathcal{N}(0,1)}\log[(1 + \sqrt{\alpha}\epsilon)^{2}].
%\end{align}
%The terms depending on $\theta$ thus cancel exactly in $\mathcal{H}(\q(w_{i}))-\mathbb{E}_{q} \log(|w_{i}|)$, proving that the KL divergence between posterior and prior is indeed independent of these parameters.
%The final term $-0.5\mathbb{E}_{\epsilon \sim \mathcal{N}(0,1)}\log[(1 + \sqrt{\alpha}\epsilon)^{2}]$ is not analytically tractable, but it can be pre-computed numerically for all values of $\alpha$ in the relevant range and then approximated using a 3rd degree polynomial in $\alpha$. As shown in figure~\ref{fig:gausskl}, this approximation is extremely close.
The final term $\mathbb{E}_{q_{\alpha}} \log(|\epsilon_{i}|)$ in this expression is not analytically tractable, but it can be pre-computed numerically for all values of $\alpha$ in the relevant range and then approximated using a 3rd degree polynomial in $\alpha$. As shown in figure~\ref{fig:gausskl}, this approximation is extremely close.

For Gaussian dropout this then gives us:
\[
-D_{KL}[\q(w_i)|p(w_i)] \approx \text{constant} + 0.5\log(\alpha) + c_{1}\alpha + c_{2}\alpha^{2} + c_{3}\alpha^{3},
\]
with
%\[
%c_{1} = 1.161451241083230, \hspace{0.5cm} c_{2} = -1.502041176441722, \hspace{0.5cm} c_{3} = 0.586299206427007.
%\]
\[
c_{1} = 1.16145124, \hspace{0.5cm} c_{2} = -1.50204118, \hspace{0.5cm} c_{3} = 0.58629921.
\]
Alternatively, we can make use of the fact that $-\mathbb{E}_{q_{\alpha}} \log(|\epsilon_{i}|) \geq 0$ for all values of $\alpha$, and use the following lower bound:
\[
-D_{KL}[\q(w_i)|p(w_i)] \geq \text{constant} + 0.5\log(\alpha).
\]
\begin{figure}[thb]
  \begin{minipage}[c]{0.5\linewidth}
    \includegraphics[width=\linewidth]{approx_gauss_kl.pdf}
  \end{minipage}\hfill
  \begin{minipage}[c]{0.5\linewidth}
    \caption{Exact and approximate negative KL divergence $-D_{KL}[\q(w_i)|p(w_i)]$ for the \textit{log-uniform} prior and \textit{Gaussian dropout} approximate posterior. Since we use an improper prior, the constant is arbitrary; we choose it so that the exact KL divergence is 0 at $\alpha=1$, which is the maximum value we allow during optimization of the lower bound. The exact and approximate results are indistinguishable on the domain shown, which corresponds to dropout rates between $p=0.05$ and $p=0.5$. For $\log(\alpha)\rightarrow-\infty$ the approximated term vanishes, so that both our polynomial approximation and the lower bound become exact.}
\label{fig:gausskl}
  \end{minipage}
\end{figure}
\FloatBarrier

\section{Variance reduction of local parameterization compared to separate weight samples}
\label{sec:low_variance_reparam}

In section~\ref{sec:fastsgvb} we claim that our local reparameterization not only is more computationally efficient than drawing separate weight matrices for each training example, but that it also leads to an estimator that has lower variance. To see why, consider the stochastic gradient estimate with respect to the posterior parameter $\sigma^{2}_{i,j}$ for a minibatch of size $M=1$. Drawing random weights $\bW$ to form the estimate, we get
\begin{align}
\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}} = \frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\frac{\epsilon_{i,j}a_{m,i}}{2\sigma_{i,j}}.
\end{align}
If, on the other hand, we form the same gradient using the local reparameterization trick, we get
\begin{align}
\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}} = \frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\frac{\zeta_{m,j}a_{m,i}^{2}}{2\sqrt{\delta_{m,j}}}.
\end{align}
We can decompose the variance of these estimators as
\[
\Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}} = \Var{b_{m,j}}{\Exp{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}|b_{m,j}}} + \Exp{b_{m,j}}{\Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}|b_{m,j}}}.
\]
Here, the first term is identical for both gradient estimators, but the second term is not. When drawing separate weight matrices per training example, it is
\begin{align}
\Exp{b_{m,j}}{\Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}|b_{m,j}}} = \Exp{b_{m,j}}{ \Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}|b_{m,j}}q} \nonumber\\
+ \Exp{b_{m,j}}{ \Exp{\q,\D}{\left(\frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}\right)^{2}|b_{m,j}}\frac{\Var{\q,\D}{\epsilon_{i,j}^{2}|b_{m,j}}a_{m,i}^{2}}{4\sigma_{i,j}^{2}} },
\end{align}
with $q = (b_{m,j}-\gamma_{m,j})^{2}a_{m,i}^{4}/4$.

When using the local reparameterization trick, we simply have
\begin{align}
\Exp{b_{m,j}}{\Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial \sigma^{2}_{i,j}}|b_{m,j}}} = \Exp{b_{m,j}}{ \Var{\q,\D}{\frac{\partial L_\D^\text{SGVB}}{\partial b_{m,j}}|b_{m,j}}q}.
\end{align}
In this case, the second term vanishes because the random variable $\zeta_{m,j}$ is uniquely determined given $b_{m,j}$, as there is only a single $\zeta_{m,j}$ for each $b_{m,j}$. Because there are many $\epsilon_{i,j}$ for each $b_{m,j}$, this variable is not uniquely determined by conditioning on $b_{m,j}$, so there is an additional (positive) contribution to the variance.

\section{Variance of stochastic gradients for variational dropout with correlated weight noise}
\label{sec:variance_corr}

In section~\ref{sec:experiments} we empirically compare the variance of the different available stochastic estimators of the gradient of our variational objective for a model regularized using variational dropout with independent weight noise (section\ref{sec:dropout_indep}). Here we present the same analysis for variational dropout with correlated weight noise (section\ref{sec:dropout_corr}). As table~\ref{tab:mnist_var_corr} shows, the relative performance of the different estimators is similar to that reported in section~\ref{sec:experiments}, although the noise level is generally much higher using this regularization method.
\begin{table}[htb]
\begin{tabular}{lllll}
 & top layer & top layer & bottom layer & bottom layer \\
stochastic gradient estimator & 10 epochs & 100 epochs & 10 epochs & 100 epochs \\
\hline
local reparameterization (ours) & $2.9\times 10^{4}$ & $5.8\times 10^{3}$ & $8.8\times 10^{2}$ & $5.3\times 10^{2}$ \\
weight sample per minibatch (standard) & $3.1\times 10^{5}$ & $2.0\times 10^{4}$ & $5.0\times 10^{3}$ & $1.1\times 10^{3}$\\
no dropout noise (minimal var.) & $7.1\times 10^{3}$ & $9.8\times 10^{2}$ & $4.3\times 10^{2}$ & $1.2\times 10^{2}$
\end{tabular}
\caption{Average empirical variance of minibatch stochastic gradient estimates (1000 examples) for a fully connected neural network, regularized by variational dropout with correlated weight noise. For this specification, the variance of the stochastic gradient estimator with separate weight samples for each data point is identical to that using the local reparameterization trick: Ignoring optimization of dropout rates, both methods correspond to the standard dropout training procedure.}
\label{tab:mnist_var_corr}
\end{table}

\section{Variance of SGVB estimator with minibatches of datapoints without replacement}
\label{sec:variance_without_replacement}
Using the indicator variables $s_i \in S$, to denote whether the $i$-th training example is included in our minibatch, and $L_{i}(\beps,\phi)$ as shorthand for $\log p(\by^i|\bx^i,\bw = f(\beps, \bphi))$, the Monte Carlo estimator \eqref{eq:objective_sgvb} may be rewritten as
\begin{align}
L_\D^\text{SGVB}(\bphi) = \frac{N}{M} \sum_{i=1}^N s_{i}L_{i}(\beps,\phi).
\end{align}
The variance of which is given by
\begin{align}
\Var{\beps,S}{L_\D^\text{SGVB}(\bphi)} \leq & \frac{N^{2}}{M^{2}} \sum_{i=1}^N \Var{S}{s_{i}}\frac{L_\D(\bphi)^{2}}{N^{2}} + \Var{\beps}{L_{i}(\beps,\phi)}\Exp{S}{s_{i}^{2}} \nonumber\\
& + \frac{2N^{2}}{M^{2}} \sum_{i>j} \Exp{S}{s_{i}}\Exp{S}{s_{j}}\Cov{\beps}{L_{i}(\beps,\phi),L_{j}(\beps,\phi)} \nonumber\\
\leq & \frac{N}{M}\left( \sum_{i=1}^N \frac{L_\D(\bphi)^{2}}{N^{2}} + \Var{\beps}{L_{i}(\beps,\phi)} \right) + 2 \sum_{i>j} \Cov{\beps}{L_{i}(\beps,\phi),L_{j}(\beps,\phi)},
\end{align}
where the inequalities arise through ignoring the negative correlations between the $s_{i}$, and where we use that $\Exp{S}{s_{i}} = \Exp{S}{s_{i}^{2}} = M/N$, and $\Var{S}{s_{i}} \leq M/N$, for $N>2M$.

As can be seen, the variance contributed by the random selection of data points from $\D$ is inversely proportional to the minibatch size $M$. However, the variance contributed by our random sampling of $\beps$ does not decrease with $M$, as we are using a single random sample for the entire minibatch and $\Cov{\beps}{L_{i}(\beps,\phi),L_{j}(\beps,\phi)}$ is positive on average\footnote{We have $\Exp{\bx^i,\by^i}{L_{i}(\beps,\phi)} = \Exp{\bx^j,\by^j}{L_{j}(\beps,\phi)}$ if the examples in the training data are i.i.d., which means that $\Exp{\bx^i,\by^i,\bx^j,\by^j}{ \Cov{\beps}{L_{i}(\beps,\phi),L_{j}(\beps,\phi)} } = \Var{\beps}{ \Exp{\bx^i,\by^i}{L_{i}(\beps,\phi)} } \geq 0$.}. In practice, this means that the variance of $L_\D^\text{SGVB}(\bphi)$ can be dominated by $\beps$ for even moderately large $M$.
\end{document}
fsgvb.bbl
\begin{thebibliography}{}

\bibitem[Ahn et~al., 2012]{ahn2012bayesian}
Ahn, S., Korattikara, A., and Welling, M. (2012).
\newblock {B}ayesian posterior sampling via stochastic gradient {F}isher
  scoring.
\newblock {\em arXiv preprint arXiv:1206.6380}.

\bibitem[Ba and Frey, 2013]{ba2013adaptive}
Ba, J. and Frey, B. (2013).
\newblock Adaptive dropout for training deep neural networks.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
  3084--3092.

\bibitem[Bayer et~al., 2015]{bayer2015fast}
Bayer, J., Karol, M., Korhammer, D., and Van~der Smagt, P. (2015).
\newblock Fast adaptive weight noise.
\newblock {\em arXiv preprint arXiv:1507.05331}.

\bibitem[Bengio, 2013]{bengio2013estimating}
Bengio, Y. (2013).
\newblock Estimating or propagating gradients through stochastic neurons.
\newblock {\em arXiv preprint arXiv:1305.2982}.

\bibitem[Bergstra et~al., 2010]{bergstra2010theano}
Bergstra, J., Breuleux, O., Bastien, F., Lamblin, P., Pascanu, R., Desjardins,
  G., Turian, J., Warde-Farley, D., and Bengio, Y. (2010).
\newblock Theano: a {CPU} and {GPU} math expression compiler.
\newblock In {\em Proceedings of the Python for Scientific Computing Conference
  (SciPy)}, volume~4.

\bibitem[Blundell et~al., 2015]{blundell2015weight}
Blundell, C., Cornebise, J., Kavukcuoglu, K., and Wierstra, D. (2015).
\newblock Weight uncertainty in neural networks.
\newblock {\em arXiv preprint arXiv:1505.05424}.

\bibitem[Gal and Ghahramani, 2015]{gal2015dropout}
Gal, Y. and Ghahramani, Z. (2015).
\newblock Dropout as a {B}ayesian approximation: Representing model uncertainty
  in deep learning.
\newblock {\em arXiv preprint arXiv:1506.02142}.

\bibitem[Graves, 2011]{graves2011practical}
Graves, A. (2011).
\newblock Practical variational inference for neural networks.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
  2348--2356.

\bibitem[Hern{\'a}ndez-Lobato and Adams, 2015]{hernandez2015probabilistic}
Hern{\'a}ndez-Lobato, J.~M. and Adams, R.~P. (2015).
\newblock Probabilistic backpropagation for scalable learning of {B}ayesian
  neural networks.
\newblock {\em arXiv preprint arXiv:1502.05336}.

\bibitem[Hinton et~al., 2012]{hinton2012improving}
Hinton, G.~E., Srivastava, N., Krizhevsky, A., Sutskever, I., and
  Salakhutdinov, R.~R. (2012).
\newblock Improving neural networks by preventing co-adaptation of feature
  detectors.
\newblock {\em arXiv preprint arXiv:1207.0580}.

\bibitem[Hinton and Van~Camp, 1993]{hinton1993keeping}
Hinton, G.~E. and Van~Camp, D. (1993).
\newblock Keeping the neural networks simple by minimizing the description
  length of the weights.
\newblock In {\em Proceedings of the sixth annual conference on Computational
  learning theory}, pages 5--13. ACM.

\bibitem[Kingma and Ba, 2015]{kingma2015adam}
Kingma, D. and Ba, J. (2015).
\newblock Adam: A method for stochastic optimization.
\newblock {\em Proceedings of the International Conference on Learning
  Representations 2015}.

\bibitem[Kingma, 2013]{kingma2013fast}
Kingma, D.~P. (2013).
\newblock Fast gradient-based inference with continuous latent variable models
  in auxiliary form.
\newblock {\em arXiv preprint arXiv:1306.0733}.

\bibitem[Kingma and Welling, 2014]{kingma2013auto}
Kingma, D.~P. and Welling, M. (2014).
\newblock Auto-encoding variational {B}ayes.
\newblock {\em Proceedings of the 2nd International Conference on Learning
  Representations}.

\bibitem[Maeda, 2014]{maeda2014bayesian}
Maeda, S.-i. (2014).
\newblock A {B}ayesian encourages dropout.
\newblock {\em arXiv preprint arXiv:1412.7003}.

\bibitem[Neal, 1995]{neal1995bayesian}
Neal, R.~M. (1995).
\newblock {\em {B}ayesian learning for neural networks}.
\newblock PhD thesis, University of Toronto.

\bibitem[Rezende et~al., 2014]{rezende2014stochastic}
Rezende, D.~J., Mohamed, S., and Wierstra, D. (2014).
\newblock Stochastic backpropagation and approximate inference in deep
  generative models.
\newblock In {\em Proceedings of the 31st International Conference on Machine
  Learning (ICML-14)}, pages 1278--1286.

\bibitem[Robbins and Monro, 1951]{robbinsmonro}
Robbins, H. and Monro, S. (1951).
\newblock A stochastic approximation method.
\newblock {\em The Annals of Mathematical Statistics}, 22(3):400--407.

\bibitem[Salimans and Knowles, 2013]{salimans2013fixedform}
Salimans, T. and Knowles, D.~A. (2013).
\newblock Fixed-form variational posterior approximation through stochastic
  linear regression.
\newblock {\em Bayesian Analysis}, 8(4).

\bibitem[Srivastava et~al., 2014]{srivastava2014dropout}
Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov,
  R. (2014).
\newblock Dropout: A simple way to prevent neural networks from overfitting.
\newblock {\em The Journal of Machine Learning Research}, 15(1):1929--1958.

\bibitem[Wan et~al., 2013]{wan2013regularization}
Wan, L., Zeiler, M., Zhang, S., Cun, Y.~L., and Fergus, R. (2013).
\newblock Regularization of neural networks using dropconnect.
\newblock In {\em Proceedings of the 30th International Conference on Machine
  Learning (ICML-13)}, pages 1058--1066.

\bibitem[Wang and Manning, 2013]{wang2013fast}
Wang, S. and Manning, C. (2013).
\newblock Fast dropout training.
\newblock In {\em Proceedings of the 30th International Conference on Machine
  Learning (ICML-13)}, pages 118--126.

\bibitem[Welling and Teh, 2011]{welling2011bayesian}
Welling, M. and Teh, Y.~W. (2011).
\newblock {B}ayesian learning via stochastic gradient {L}angevin dynamics.
\newblock In {\em Proceedings of the 28th International Conference on Machine
  Learning (ICML-11)}, pages 681--688.

\end{thebibliography}

fsgvb.bbl
\begin{thebibliography}{}

\bibitem[Ahn et~al., 2012]{ahn2012bayesian}
Ahn, S., Korattikara, A., and Welling, M. (2012).
\newblock {B}ayesian posterior sampling via stochastic gradient {F}isher
  scoring.
\newblock {\em arXiv preprint arXiv:1206.6380}.

\bibitem[Ba and Frey, 2013]{ba2013adaptive}
Ba, J. and Frey, B. (2013).
\newblock Adaptive dropout for training deep neural networks.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
  3084--3092.

\bibitem[Bayer et~al., 2015]{bayer2015fast}
Bayer, J., Karol, M., Korhammer, D., and Van~der Smagt, P. (2015).
\newblock Fast adaptive weight noise.
\newblock {\em arXiv preprint arXiv:1507.05331}.

\bibitem[Bengio, 2013]{bengio2013estimating}
Bengio, Y. (2013).
\newblock Estimating or propagating gradients through stochastic neurons.
\newblock {\em arXiv preprint arXiv:1305.2982}.

\bibitem[Bergstra et~al., 2010]{bergstra2010theano}
Bergstra, J., Breuleux, O., Bastien, F., Lamblin, P., Pascanu, R., Desjardins,
  G., Turian, J., Warde-Farley, D., and Bengio, Y. (2010).
\newblock Theano: a {CPU} and {GPU} math expression compiler.
\newblock In {\em Proceedings of the Python for Scientific Computing Conference
  (SciPy)}, volume~4.

\bibitem[Blundell et~al., 2015]{blundell2015weight}
Blundell, C., Cornebise, J., Kavukcuoglu, K., and Wierstra, D. (2015).
\newblock Weight uncertainty in neural networks.
\newblock {\em arXiv preprint arXiv:1505.05424}.

\bibitem[Gal and Ghahramani, 2015]{gal2015dropout}
Gal, Y. and Ghahramani, Z. (2015).
\newblock Dropout as a {B}ayesian approximation: Representing model uncertainty
  in deep learning.
\newblock {\em arXiv preprint arXiv:1506.02142}.

\bibitem[Graves, 2011]{graves2011practical}
Graves, A. (2011).
\newblock Practical variational inference for neural networks.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
  2348--2356.

\bibitem[Hern{\'a}ndez-Lobato and Adams, 2015]{hernandez2015probabilistic}
Hern{\'a}ndez-Lobato, J.~M. and Adams, R.~P. (2015).
\newblock Probabilistic backpropagation for scalable learning of {B}ayesian
  neural networks.
\newblock {\em arXiv preprint arXiv:1502.05336}.

\bibitem[Hinton et~al., 2012]{hinton2012improving}
Hinton, G.~E., Srivastava, N., Krizhevsky, A., Sutskever, I., and
  Salakhutdinov, R.~R. (2012).
\newblock Improving neural networks by preventing co-adaptation of feature
  detectors.
\newblock {\em arXiv preprint arXiv:1207.0580}.

\bibitem[Hinton and Van~Camp, 1993]{hinton1993keeping}
Hinton, G.~E. and Van~Camp, D. (1993).
\newblock Keeping the neural networks simple by minimizing the description
  length of the weights.
\newblock In {\em Proceedings of the sixth annual conference on Computational
  learning theory}, pages 5--13. ACM.

\bibitem[Kingma and Ba, 2015]{kingma2015adam}
Kingma, D. and Ba, J. (2015).
\newblock Adam: A method for stochastic optimization.
\newblock {\em Proceedings of the International Conference on Learning
  Representations 2015}.

\bibitem[Kingma, 2013]{kingma2013fast}
Kingma, D.~P. (2013).
\newblock Fast gradient-based inference with continuous latent variable models
  in auxiliary form.
\newblock {\em arXiv preprint arXiv:1306.0733}.

\bibitem[Kingma and Welling, 2014]{kingma2013auto}
Kingma, D.~P. and Welling, M. (2014).
\newblock Auto-encoding variational {B}ayes.
\newblock {\em Proceedings of the 2nd International Conference on Learning
  Representations}.

\bibitem[Maeda, 2014]{maeda2014bayesian}
Maeda, S.-i. (2014).
\newblock A {B}ayesian encourages dropout.
\newblock {\em arXiv preprint arXiv:1412.7003}.

\bibitem[Neal, 1995]{neal1995bayesian}
Neal, R.~M. (1995).
\newblock {\em {B}ayesian learning for neural networks}.
\newblock PhD thesis, University of Toronto.

\bibitem[Rezende et~al., 2014]{rezende2014stochastic}
Rezende, D.~J., Mohamed, S., and Wierstra, D. (2014).
\newblock Stochastic backpropagation and approximate inference in deep
  generative models.
\newblock In {\em Proceedings of the 31st International Conference on Machine
  Learning (ICML-14)}, pages 1278--1286.

\bibitem[Robbins and Monro, 1951]{robbinsmonro}
Robbins, H. and Monro, S. (1951).
\newblock A stochastic approximation method.
\newblock {\em The Annals of Mathematical Statistics}, 22(3):400--407.

\bibitem[Salimans and Knowles, 2013]{salimans2013fixedform}
Salimans, T. and Knowles, D.~A. (2013).
\newblock Fixed-form variational posterior approximation through stochastic
  linear regression.
\newblock {\em Bayesian Analysis}, 8(4).

\bibitem[Srivastava et~al., 2014]{srivastava2014dropout}
Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov,
  R. (2014).
\newblock Dropout: A simple way to prevent neural networks from overfitting.
\newblock {\em The Journal of Machine Learning Research}, 15(1):1929--1958.

\bibitem[Wan et~al., 2013]{wan2013regularization}
Wan, L., Zeiler, M., Zhang, S., Cun, Y.~L., and Fergus, R. (2013).
\newblock Regularization of neural networks using dropconnect.
\newblock In {\em Proceedings of the 30th International Conference on Machine
  Learning (ICML-13)}, pages 1058--1066.

\bibitem[Wang and Manning, 2013]{wang2013fast}
Wang, S. and Manning, C. (2013).
\newblock Fast dropout training.
\newblock In {\em Proceedings of the 30th International Conference on Machine
  Learning (ICML-13)}, pages 118--126.

\bibitem[Welling and Teh, 2011]{welling2011bayesian}
Welling, M. and Teh, Y.~W. (2011).
\newblock {B}ayesian learning via stochastic gradient {L}angevin dynamics.
\newblock In {\em Proceedings of the 28th International Conference on Machine
  Learning (ICML-11)}, pages 681--688.

\end{thebibliography}
