Learning and inference are useful concepts rooted in the graphical model literature. We are generally trying to create probabilistic models for various systems, including e.g. many of the generative models that we see today.These are notes of An Introduction to Variational Methods for Graphical Models from Jordan et al., 1999 and Auto-Encoding Variational Bayes from Kingma and Welling, 2013.IntroFor a given graphical model, we have known evidence nodes 𝐸 and hidden nodes 𝐻, and parameters modeling their relationships. Treat nodes as just random variables that can influence each other. We will discover how to set up graphs to model the dependencies.We are interested in two problems.1.In inference, we are interested in the posterior distribution𝑝(𝐻|𝐸)=𝑝(𝐻,𝐸)𝑝(𝐸)Here, we fix the model parameters and observe evidence 𝐸 to infer hidden 𝐻.2.In learning, we are interested in the probability of the observed evidence, alternatively known as the likelihood when seen as a function of the model parameters.𝑝(𝐸)=∑{𝐻}𝑝(𝐻,𝐸)We will first discuss exact inference methods, where they struggle, and motivate variational methods.Exact inferenceGraphical models have two main variants—directed and undirected.Directed graphical models (also known as Bayesian networks) are set up as directed acyclic graphs (DAGs).𝑝(𝑠1)𝑝(𝑠2)𝑝(𝑠3|𝑠1,𝑠2)𝑝(𝑠4|𝑠3)𝑝(𝑠5|𝑠3)𝑝(𝑠6|𝑠4,𝑠5)Figure 1: A directed graphical model with local conditional probabilities.Each node is labeled 𝑠𝑖 and has parent nodes 𝑠𝜋(𝑖), where 𝜋 returns the set of parent indices for 𝑖. We usually write the joint distributions for these graphs as 𝑝(𝑠) where 𝑠 contains all nodes; for the directed graphical model, we have:𝑝(𝑠)=∏𝑁𝑖=1𝑝(𝑠𝑖|𝑠𝜋(𝑖))chain ruleNote that each node 𝑠𝑖 is conditionally independent of its non-descendants given its parents 𝑠𝜋(𝑖). In other words, we should think of arrows as dependencies. If 𝑠𝑗 points to 𝑠𝑖, then 𝑠𝑖 depends on 𝑠𝑗.Undirected graphical models (also known as Markov random fields) associate what we call “potentials” with cliques of the graph 1 . A potential is a function on clique configurations. In other words, for a given value assignment to all the clique’s nodes 𝑐𝑖, the potential will give that clique some value 𝜙𝑖(𝑐𝑖).𝑠1𝑠2𝑠3𝑠4𝑠5𝑠6𝜙1(𝑐1)𝜙2(𝑐2)𝜙3(𝑐3)Figure 2: An undirected graphical model with a potential on each clique.In this example, we have:𝑐1={𝑠1,𝑠2,𝑠3}𝑐2={𝑠3,𝑠4,𝑠5}𝑐3={𝑠4,𝑠5,𝑠6}For undirected graphs, our joint distribution will be the normalized product of all our clique potentials:𝑝(𝑠)=1𝑍∏𝑀𝑖=1𝜙𝑖(𝑐𝑖)where 𝑀 is the number of cliques and 𝑍 is our normalizing factor:𝑍=∑{𝑠}{∏𝑀𝑖=1𝜙𝑖(𝑐𝑖)}𝑍 is also known as a partition function, taken from statistical mechanics. {𝑠} refers to all possible configurations of our nodes (i.e. value assignments to the node random variables). For example, if each node 𝑠∈𝒳︀, then {𝑠}=𝒳︀×⋯×𝒳︀.The primary algorithm we are interested in for exact inference is the junction tree algorithm. We won’t cover it in detail here, but here is an overview.1.Moralization: We first compile a directed graphical model into a corresponding undirected graphical model.We observe that in both directed and undirected graphs, we are using some kind of product of local functions. For directed, we use 𝑝(𝑠𝑖|𝑠𝜋(𝑖)), and in undirected, we use potentials. This conditional probability is almost a potential function, but note that the parents in 𝑠𝜋(𝑖) may not share the same clique.So, we add new undirected edges between parents and drop all other arrows. This gives us a “moral graph”.2.Triangulation: for a naive moral graph, we run into issues writing joint probabilities for 4-cycles.𝐴𝐵𝐶𝐷Figure 3: An undirected 4-cycle.We try to write the joint probability in a sequential manner:𝑝(𝐴,𝐵,𝐶,𝐷)=𝑝(𝐴)𝑝(𝐵|𝐴)𝑝(𝐶|𝐵)𝑝(𝐷|𝐶)Here, 𝑝(𝐷|𝐶) will not work since 𝐷 depends on 𝐴. Instead, we can insert a chord from 𝐵 to 𝐷:𝐴𝐵𝐶𝐷Figure 4: The same 4-cycle, triangulated with a chord from 𝐵 to 𝐷.Now, going sequentially from 𝐴 to 𝐶, our joint probability is𝑝(𝐴,𝐵,𝐶,𝐷)=𝑝(𝐴)𝑝(𝐵,𝐷|𝐴)𝑝(𝐶|𝐵,𝐷)We can handle 𝐵 and 𝐷 together as (𝐵,𝐷) since 𝐴 and 𝐶 are conditionally independent given 𝐵,𝐷. Intuitively, if we fix 𝐵 and 𝐷 to some values, then wiggling 𝐴 can’t reach 𝐶 through any intermediate edges.Triangulation is the process of making the graph chordal, which means the graph does not admit any 4-cycles without a chord. The purpose of this triangulation is to allow us to arrange the cliques into a junction tree.𝐴𝐵𝐹𝐶𝐺𝐷𝐸𝐻𝐵𝐶𝐶𝐸𝐵𝐸𝐵𝐺𝐸𝐺𝐴𝐵𝐶𝐶𝐷𝐸𝐵𝐶𝐸𝐵𝐸𝐺𝐵𝐹𝐺𝐸𝐺𝐻Figure 5: A triangulated graph (left) and a junction tree of its cliques (right), with each edge labeled by shared nodes.The core property of a junction tree is that if cliques 𝑐𝑖 and 𝑐𝑗 both contain some node 𝐴, then every clique on the unique path between them also has 𝐴. We can visually see this in the above example—since 𝐵 is in the top left and top right cliques, it exists in the middle two as well.There are standard algorithms to build such a tree (Maximum Cardinality Search and Kruskal’s algorithm), we will omit these. Similarly, we can use message passing algorithms to perform exact inference on a built junction tree by enforcing local consistency bubbling up into global consistency. See sum-product message passing algorithms for a more detailed explanation.Limits of exact inferenceThe original Jordan et al. 1999 paper describes several examples, including QMR-DT diagnosis models, neural networks, Boltzmann machines, and variations of Hidden Markov Models (HMMs).Here, we will consider a basic modified neural network.𝐻𝐸𝑠1𝑠2𝑠3Figure 6: A neural network with two hidden nodes and one evidence node.In this directed graphical model, we describe conditional probabilities as𝑝(𝑠𝑖=1|𝑠𝜋(𝑖))=11+exp{−∑𝑗∈𝜋(𝑖)𝜃𝑖𝑗𝑠𝑗}=𝜎(∑𝑗∈𝜋(𝑖)𝜃𝑖𝑗𝑠𝑗)where 𝜃𝑖𝑗 are model parameters and we exclude bias for simplicity. We can call this form a sigmoid belief network where 𝑠𝑖∈{0,1} with logistic non-linearities.If we wish to perform exact inference, we need to moralize this network, which creates a dependency between the parents.𝐻𝐸𝑠1𝑠2𝑠3Figure 7: Moralized sigmoid belief network.We can try to run exact inference for the joint and marginal distributions. First, our joint distribution is modeled as:𝑝(𝑠)=∏𝑖𝑝(𝑠𝑖|𝑠𝜋(𝑖))Going backward from the evidence nodes, our exact joint is:𝑝(𝑠)=𝑝(𝑠3|𝑠1,𝑠2)𝑝(𝑠1)𝑝(𝑠2)Using our expression for conditional probabilities, we have:𝑝(𝑠1)=12𝑝(𝑠2)=12𝑝(𝑠3|𝑠1,𝑠2)=(𝜎(𝜃13𝑠1+𝜃23𝑠2))𝑠3(1−𝜎(𝜃13𝑠1+𝜃23𝑠2))1−𝑠3The final expression is the Bernoulli mass. Our joint distribution is just𝑝(𝑠)=14(𝜎(𝜃13𝑠1+𝜃23𝑠2))𝑠3(1−𝜎(𝜃13𝑠1+𝜃23𝑠2))1−𝑠3Now the marginal distribution:𝑝(𝐸)=∑{𝐻}∏𝑖𝑝(𝑠𝑖|𝑠𝜋(𝑖))=∑{𝐻}14(𝜎(𝜃13𝑠1+𝜃23𝑠2))𝑠3(1−𝜎(𝜃13𝑠1+𝜃23𝑠2))1−𝑠3Here, we need to enumerate over all configurations of 𝐻 (𝑠1 and 𝑠2). Since both are binary, we have 2|𝐻|=22=4 configurations:𝑝(𝑠3=1)=14𝜎(0)⏟{0,0}+14𝜎(𝜃13)⏟{1,0}+14𝜎(𝜃23)⏟{0,1}+14𝜎(𝜃13+𝜃23)⏟{1,1}𝑝(𝑠3=0)=1−𝑝(𝑠3=1)Note the issue—we had to enumerate over 2|𝐻| configurations. Since the logistic function doesn’t distribute over addition/multiplication, we cannot swap the sum and product to make this computation tractable:𝜎(𝑎+𝑏)≠𝜎(𝑎)+𝜎(𝑏)𝜎(𝑎+𝑏)≠𝜎(𝑎)𝜎(𝑏)In other graphical models like HMMs and Boltzmann machines, we encounter similar issues when simplifying. From here, we will motivate variational methods.A brief intro to variational methodsVariational methods are approximation methods that transform a complex problem into a simpler one. We expand the original problem to include what we call “variational parameters” that we will also fit to our simpler problem.The terminology originates from calculus of variations, but similar to Jordan et al. 1999, we will refrain from starting from those fundamentals.We’ll start with simple functions to illustrate the method.Consider log𝑥. Notice that we can write this concave function as a pointwise minimum of linear functions by introducing a variational parameter 𝜆:log𝑥=min𝜆{𝜆𝑥−log𝜆−1}Given any 𝑥, we can solve for 𝜆 by optimizing our linear function:dd𝜆(𝜆𝑥−log𝜆−1)=0𝑥−1𝜆∗=0𝜆∗=1𝑥Then, substituting 𝜆∗ yields log𝑥. But notice! We are effectively tuning the slope 𝜆 of linear function 𝜆𝑥 (with varying intercepts).00.511.522.53−3−2−1012𝑥Figure 8: Solid log𝑥 with dotted tangent lines 𝜆𝑥−log𝜆−1.This family of linear functions forms the upper envelope of log𝑥, giving an upper bound to the function:log𝑥≤𝜆𝑥−log𝜆−1for all 𝜆, with 𝜆∗ giving equality for a fixed 𝑥.Intuitively, we have replaced a nonlinear function with a linear function with the cost of an extra parameter 𝜆 that must be fit per 𝑥.Convex dualityHow generally can we do this? Pulling from convex analysis, a concave function 𝑓(𝑥) can be represented by the pointwise minimum of a family of conjugate (also called dual) function 𝑓∗(𝜆):𝑓(𝑥)=min𝜆{𝜆⊤𝑥−𝑓∗(𝜆)}where 𝑥 and 𝜆 are now vectors. This new conjugate function has the similar form𝑓∗(𝜆)=min𝑥{𝜆⊤𝑥−𝑓(𝑥)}Intuitively, what are we minimizing here for the conjugate? We can break our objective function into two pieces: linear 𝜆𝑥 and concave 𝑓(𝑥) (for the univariate case).𝑓(𝑥)=log𝑥𝜆𝑥00.511.522.53−2−10123𝑥Figure 9: Our conjugate function 𝑓∗(𝜆) minimizes the gap between 𝜆𝑥 and log𝑥.In this visualization, we see that for a given slope 𝜆 and point 𝑥, there is a gap between 𝜆𝑥 and 𝑓(𝑥). The conjugate function minimizes this gap for each 𝑥 (pointwise) by wiggling the variational parameter 𝜆. We can see that this gap will be minimized when the slopes match.Neatly, it now makes sense why 𝑓∗(𝜆) is the intercept in 𝜆⊤𝑥−𝑓∗(𝜆). 𝑓∗(𝜆) represents the gap between the two when the slopes match, which is exactly the amount we need to shift 𝜆𝑥 down so it just touches 𝑓(𝑥).Similarly, for convex 𝑓(𝑥), we have:𝑓(𝑥)=max𝜆{𝜆⊤𝑥−𝑓∗(𝜆)}𝑓∗(𝜆)=max𝑥{𝜆⊤𝑥−𝑓(𝑥)}Here are the general steps you can take to find upper/lower bounds on a function of interest:1.If the function is already convex/concave, find its conjugate function.2.If not, find an invertible transform to make it convex/concave.3.Then, find the conjugate function.4.Transform the conjugate back with the inverse.We will typically use a transform function like the logarithm that has useful properties.As an example of this process, consider logistic function 𝑓(𝑥)=11+𝑒−𝑥. While the logistic function is neither convex nor concave, it is log concave2 .𝑔(𝑥)=log𝜎(𝑥)Now, we find the conjugate function to 𝑔(𝑥):𝑔∗(𝜆)=min𝑥{𝜆𝑥−𝑔(𝑥)}=min𝑥{𝜆𝑥−log𝜎(𝑥)}Optimizing, we find3 :𝑥∗=−log𝜆+log(1−𝜆)𝑔∗(𝜆)=−𝜆log𝜆−(1−𝜆)log(1−𝜆)=𝐻(𝜆)where 𝐻(𝜆) is just the entropy of a Bernoulli variable with parameter 𝜆.Substituting in 𝑔(𝑥), we have:𝑔(𝑥)=min𝜆{𝜆𝑥−𝐻(𝜆)}𝜎(𝑥)≤𝑒𝜆𝑥−𝐻(𝜆)where we used the exponential (inverse of log) to obtain a upper bound on the original logistic function.𝑔(𝑥)−4−2024−5−4−3−2−101𝑥𝜎(𝑥)−4−202400.20.40.60.811.21.4𝑥Figure 10: Left: solid log logistic function with dashed family of linear bounding functions. Right: solid logistic function with dashed family of exponential upper bounds 𝑒𝜆𝑥−𝐻(𝜆); comes from exponentiating the left graph.Variational approximations of joint and conditional probabilitiesRecall that 𝑝(𝐻|𝐸) for inference and 𝑝(𝐸) for learning are of interest to us. In our neural network case, we noticed that conditional probabilities in 𝑝(𝐸) became intractable. Now, we will use our understanding of convex duality to rewrite them with tractable variational bounds4 .𝑝(𝑠𝑖|𝑠𝜋(𝑖))≤𝑝(𝑠𝑖|𝑠𝜋(𝑖),𝜆𝑖)𝑝(𝑠𝑖|𝑠𝜋(𝑖))=𝜎(∑𝑗∈𝜋(𝑖)𝜃𝑖𝑗𝑠𝑗)≤exp{𝜆𝑖(∑𝑗∈𝜋(𝑖)𝜃𝑖𝑗𝑠𝑗)−𝐻(𝜆𝑖)}𝑝(𝐸)=∑{𝐻}∏𝑖𝑝(𝑠𝑖|𝑠𝜋(𝑖))≤∑{𝐻}∏𝑖exp{𝜆𝑖(∑𝑗∈𝜋(𝑖)𝜃𝑖𝑗𝑠𝑗)−𝐻(𝜆𝑖)}Since probabilities are within (0,1), a product of upper bounds is also a upper bound for 𝑝(𝐸). Notice here: instead of a sum of product of logistic functions, we now have a sum of product of exponentials. We can simplify this5 :𝑝(𝐸)≤𝐾1𝐾2∏𝑗∈𝐻(1+exp{∑𝑖∈𝑐(𝑗)𝜆𝑖𝜃𝑖𝑗})where 𝐾1 and 𝐾2 are constants computed in 𝑂(𝑁+|𝐸|) time, where 𝑁 is the number of nodes and |𝐸| is the number of edges. Similarly, the final time complexity is 𝑂(𝑁+|𝐸|)—the double loop iterates over all hidden nodes and their children, so we must at least go through all the edges.Based on the variational approach, we reduced the time complexity from 𝑂(2|𝐻|) to be tractable in the number of edges. We can use a similar approach for the joint distribution. Our final question here is: what is the error in the approximation?In the joint distribution, we are modeling𝑝(𝑠)=∏𝑖𝑝(𝑠𝑖|𝑠𝜋(𝑖))Since all of 𝑠 is constant here, we know the exact value of each 𝑠𝑗∈𝑠𝜋(𝑖).𝜆∗=argmin𝜆{𝜆𝑥−𝐻(𝜆)}𝑥𝑖=∑𝑗∈𝜋(𝑖)𝜃𝑖𝑗𝑠𝑗So, we can fit the variational parameters 𝜆 for each argument 𝑥𝑖 and recover the exact joint probability with no error.In the marginal distribution, we are varying over {𝐻}, so the parents 𝑠𝑗 vary in every term of the sum (i.e. 𝑥𝑖 is a function of configuration 𝑠).𝑥𝑖(𝑠)=∑𝑗∈𝜋(𝑖)𝜃𝑖𝑗𝑠𝑗This complication has two implications:•In the algebraic simplification, we cannot factor out 𝐾1=∏𝑖exp{−𝐻(𝜆𝑖)} if 𝜆𝑖 depends on 𝑠.•We need to optimize for 𝜆𝑖 using the evidence 𝐸 to remove 𝜆𝑖 dependence on 𝐻.Therefore, our upper bound may not achieve equality in this case. Instead, we can reason about how tight the bound is based on how sensitive the model is to our fit of 𝜆𝑖.Variational algorithmsWe now turn to algorithms we can use to apply variational methods to inference problems. Jordan et al. 1999 describes sequential and block approaches. We are interested in the block approach as it will aid us in neural methods for variational inference.The block approach identifies some subgraph or substructure within our broader graphical model where we can easily deploy exact inference. We will then use this “family” of subgraphs as our variational bound.Concretely, recall that 𝑝(𝑠) is our joint distribution. Since we desire 𝑝(𝐻|𝐸), we will introduce some 𝑞(𝐻|𝐸,𝜆) where 𝜆 are variational parameters. Note that the graph for 𝑞 may not match the graph for 𝑝. To select the right distribution from the 𝑞 family, we will minimize the relative entropy between the two with respect to 𝜆:𝜆∗=argmin𝜆{𝐷(𝑞(𝐻|𝐸,𝜆)‖𝑝(𝐻|𝐸))}where relative entropy (KL divergence) is defined as𝐷(𝑞‖𝑝)=∑{𝑆}𝑞(𝑆)log𝑞(𝑆)𝑝(𝑆)discrete=∫𝑞(𝑆)log𝑞(𝑆)𝑝(𝑆)d𝑆continuousUsing 𝑞, we can establish a lower bound on the log marginal:log𝑝(𝐸)=log∑{𝐻}𝑝(𝐻,𝐸)=log∑{𝐻}𝑞(𝐻|𝐸,𝜆)𝑝(𝐻,𝐸)𝑞(𝐻|𝐸,𝜆)=log𝔼𝑞(𝐻|𝐸,𝜆)[𝑝(𝐻,𝐸)𝑞(𝐻|𝐸,𝜆)]≥𝔼𝑞(𝐻|𝐸,𝜆)[log𝑝(𝐻,𝐸)𝑞(𝐻|𝐸,𝜆)]Jensen's inequality, concaveThe error of the bound can be found as follows:log𝑝(𝐸)−𝔼𝑞(𝐻|𝐸,𝜆)[log𝑝(𝐻,𝐸)𝑞(𝐻|𝐸,𝜆)]=𝔼𝑞(𝐻|𝐸,𝜆)[log𝑝(𝐸)]−𝔼𝑞(𝐻|𝐸,𝜆)[log𝑝(𝐻,𝐸)𝑞(𝐻|𝐸,𝜆)]=𝔼𝑞(𝐻|𝐸,𝜆)[log𝑝(𝐸)𝑞(𝐻|𝐸,𝜆)𝑝(𝐸)𝑝(𝐻|𝐸)]=𝔼𝑞(𝐻|𝐸,𝜆)[log𝑞(𝐻|𝐸,𝜆)𝑝(𝐻|𝐸)]=𝐷(𝑞(𝐻|𝐸,𝜆)‖𝑝(𝐻|𝐸))Therefore, by choosing 𝜆 that minimizes this error, we obtain the tightest lower bound.Auto-encoding variational BayesWe finally turn to modern variational techniques. We will focus on variational autoencoders, a general framework for learning and inference for intractable posteriors with large datasets (Kingma and Welling, 2013).The paper assumes a simple directed graphical model with many independent hidden 𝑧 and evidence 𝑥.𝑧𝑥×𝑁Figure 11: 𝑁 independent copies of hidden 𝑧 generating observed 𝑥.The probabilistic graph of interest is modeled by 𝑝 parameterized by 𝜃:𝑝𝜃(𝑥)intractable evidence𝑝𝜃(𝑧)prior𝑝𝜃(𝑥|𝑧)conditional likelihood𝑝𝜃(𝑧|𝑥)intractable posterior𝑝𝜃(𝑥,𝑧)jointWe wish to solve the following issues:•The true data-generating 𝑝𝜃 is usually intractable for our desired inference and learning problems. For example, if we choose a Gaussian prior for 𝑧 and a nonlinear neural network for 𝑝𝜃(𝑥|𝑧), we can try computing𝑝𝜃(𝑥)=∫𝑝𝜃(𝑥|𝑧)𝑝𝜃(𝑧)d𝑧but we cannot find a closed-form analytical solution due to the non-linearities in the neural network for 𝑝𝜃(𝑥|𝑧). Typical integral estimation approaches require an intractable number of forward passes per gradient step. Inference requires this intractable marginal distribution, and is analytically intractable as well.•Dataset sizes are usually massive. We usually would like a minibatch-compatible approach to fit our parameters, typically using gradient descent.Instead, we will use a variational approach.Assume we have a dataset 𝑋={𝑥(𝑖)}𝑁𝑖=1 with 𝑁 i.i.d. samples of our variable 𝑥 generated by:1.𝑧(𝑖)∼𝑝𝜃(𝑧)2.𝑥(𝑖)∼𝑝𝜃(𝑥|𝑧)where both densities are differentiable with respect to 𝜃 and 𝑧. The true 𝜃∗ and 𝑧(𝑖) are hidden.We employ a modified block approach by introducing 𝑞𝜙(𝑧|𝑥) as our variational model to solve our inference problem and approximate the true posterior 𝑝𝜃(𝑧|𝑥).The variational boundWe follow the steps from Jordan et al., 1999 exactly to find our log marginal (we wish to maximize the probability of seeing our data 𝑝𝜃(𝑥)):log𝑝𝜃(𝑥)=𝐷(𝑞𝜙(𝑧|𝑥)‖𝑝𝜃(𝑧|𝑥))⏟distance from true posterior+𝔼𝑞𝜙(𝑧|𝑥)[log𝑝𝜃(𝑥,𝑧)𝑞𝜙(𝑧|𝑥)]⏟variational boundWe can further break down the variational lower bound:𝔼𝑞𝜙(𝑧|𝑥)[log𝑝𝜃(𝑥,𝑧)𝑞𝜙(𝑧|𝑥)]=𝔼𝑞𝜙(𝑧|𝑥)[log𝑝𝜃(𝑥|𝑧)+log𝑝𝜃(𝑧)−log𝑞𝜙(𝑧|𝑥)]=𝔼𝑞𝜙(𝑧|𝑥)[log𝑝𝜃(𝑧)𝑞𝜙(𝑧|𝑥)]+𝔼𝑞𝜙(𝑧|𝑥)[log𝑝𝜃(𝑥|𝑧)]=−𝔼𝑞𝜙(𝑧|𝑥)[log𝑞𝜙(𝑧|𝑥)𝑝𝜃(𝑧)]+𝔼𝑞𝜙(𝑧|𝑥)[log𝑝𝜃(𝑥|𝑧)]ℒ︀(𝜃,𝜙)=−𝐷(𝑞𝜙(𝑧|𝑥)‖𝑝𝜃(𝑧))⏟distance from prior+𝔼𝑞𝜙(𝑧|𝑥)[log𝑝𝜃(𝑥|𝑧)]⏟reconstruction errorWe wish to maximize our lower bound, so we maximize ℒ︀(𝜃,𝜙). This maximization will:1.Minimize the relative entropy between our variational model and the prior on 𝑧2.Maximize our generative model’s ability to reconstruct 𝑥 from hidden 𝑧, which itself was constructed by 𝑞𝜙 from 𝑥.max𝜃,𝜙ℒ︀(𝜃,𝜙)Gradient estimationUsing the variational lower bound, we now have an optimization problem that we can optimize via stochastic gradient descent.ℒ︀(𝜃,𝜙)=−𝐷(𝑞𝜙(𝑧|𝑥)‖𝑝𝜃(𝑧))+𝔼𝑞𝜙(𝑧|𝑥)[log𝑝𝜃(𝑥|𝑧)]Usually, 𝑝𝜃(𝑧) is fixed to be Gaussian, so we can use standard autodifferentiation libraries to compute the gradient of ℒ︀ with respect to 𝜃.To compute the gradient with respect to 𝜙, notice that both terms are expectations over the 𝑞𝜙(𝑧|𝑥) distribution with 𝑧 being used inside the expectation. Since 𝑧 is stochastically sampled, we need a different way to estimate the gradient. The standard Monte Carlo gradient estimator using the REINFORCE trick is6 :∇𝜙𝔼𝑞𝜙(𝑧)[𝑓(𝑧)]≈1𝐿∑𝐿ℓ=1𝑓(𝑧(ℓ))∇𝜙log𝑞𝜙(𝑧(ℓ))where 𝑓(𝑧) is written by inspecting ℒ︀. However, this estimator has high variance, and so we wish to find a better estimator.The reparameterization trick used in the paper can be generally described by my previous post on Reparameterization tricks. We will briefly describe the method here, and give an example for the Gaussian case.We attempt to “standardize” 𝑧 by establishing a relationship between 𝑧 and 𝜀 that moves the stochasticity of 𝑧 into 𝜀. Formally, we choose a standardization function 𝒮︀𝜙𝜀=𝒮︀𝜙(𝑧)𝑧=𝒮︀−1𝜙(𝜀)𝜀∼𝑝(𝜀)Then, we can write the gradient as∇𝜙𝔼𝑞𝜙(𝑧)[𝑓(𝑧)]=𝔼𝑞𝜙(𝑧)[−(∇𝑧𝑓(𝑧))(∇𝑧𝒮︀𝜙(𝑧))−1(∇𝜙𝒮︀𝜙(𝑧))]In the paper, they require the standardization function to be differentiable and invertible, but the general requirement only requires differentiability.As an example, consider variational autoencoders (VAEs). We choose Gaussian distributions for our desired quantities:𝑝𝜃(𝑧)=𝒩︀(𝑧;0,𝐼)𝑝𝜃(𝑥|𝑧)=𝒩︀(𝑥;𝜇𝜃(𝑧),Σ𝜃(𝑧))MLP𝑞𝜙(𝑧|𝑥)=𝒩︀(𝑧;𝜇𝜙(𝑥),diag(𝜎2𝜙(𝑥)))MLPThen, we can write𝒮︀𝜙(𝑧)=𝑧−𝜇𝜎∇𝑧𝒮︀𝜙(𝑧)=1𝜎∇𝜇𝒮︀𝜙(𝑧)=−1𝜎∇𝜎𝒮︀𝜙(𝑧)=−𝑧−𝜇𝜎2∇𝜇𝑧=1∇𝜎𝑧=𝑧−𝜇𝜎=𝜀We can substitute these quantities to find the gradient with respect to 𝜙. Alternatively, we can use 𝑧=𝜇𝜙+𝜎𝜙⊙𝜀 where 𝜀∼𝒩︀(0,𝐼). This gives the differentiable lossℒ︀(𝜃,𝜙)=12∑𝑗(1+log(𝜎2𝜙𝑗)−𝜇2𝜙𝑗−𝜎2𝜙𝑗)+1𝐿∑ℓlog𝑝𝜃(𝑥|𝑧(ℓ))DiscussionVariational methods allow us to solve intractable inference and learning problems on probabilistic graphical models. These have been extended to deep learning techniques in recent years, as seen with variational autoencoders. Hopefully, the ideas explored here will be useful when understanding generative models broadly.