summaries: This manuscript proposes a new approach to fitting Markov Random Fields (MRFs).  The general structure of the algorithm is amenable to many MRF structures and can be fairly straightforwardly applied to learning on a wide variety of problems.  The theoretical analysis supports that the algorithm is reasonable.  Experimental results show strong results on several different MRF models, albeit on relatively small problems. The work proposes using variational distributions to model the model the inference of latent variables and model the partition function building on NVIL, thereby providing an algorithm that would work on general MRFs for both inference and learning. Since the two terms in the NLL are opposite in sign, it is a minimax operation and GAN like adversial training can be used. The paper shows providing tighter results to estimate the log partition function and comparisons on the digits dataset and Anneal importance sampling. This paper presents a black-box style learning algorithm for Markov Random Fields (MRF). The approach doubles down on the variational approach with variational approximations for both the positive phase and negative phase of the log likelihood objective function. For the negative phase, the authors use two separate variational approximations, one of which involves the modeling of the latent variable prior under the approximating distribution,
abstractText: We propose a black-box algorithm called Adversarial Variational Inference and Learning (AdVIL) to perform inference and learning in a general Markov random field (MRF). AdVIL employs two variational distributions to approximately infer the latent variables and estimate the partition function of an MRF, respectively. The two variational distributions provide an estimate of the negative log-likelihood of the MRF as a minimax optimization problem, which is solved by stochastic gradient descent. AdVIL is proven convergent under certain conditions. On one hand, compared to the contrastive divergence, AdVIL requires minimal assumptions about the model structure and can deal with a broader family of MRFs. On the other hand, compared to existing black-box methods, AdVIL provides a tighter estimate of the log partition function and achieves much better empirical results.
authors.name: RELIEVE YOUR HEADACHE
authors.name: TAKE ADVIL
authors.name: Chongxuan Li
authors.name: Chao Du
authors.name: Kun Xu
authors.name: Max Welling
authors.name: Jun Zhu
authors.name: Bo Zhang
references.authors: Martı́n Abadi Paul Barham Jianmin Chen Zhifeng Chen Andy Davis Jeffrey Dean Matthieu Devin Sanjay Ghemawat Geoffrey Irving Michael Isard
references.title: TensorFlow: A system for large-scale machine
references.authors: David H Ackley Geoffrey E Hinton Terrence J Sejnowski
references.title: A learning algorithm for boltzmann machines
references.venue: Cognitive science,
references.authors: Léon Bottou Frank E Curtis Jorge Nocedal
references.title: Optimization methods for large-scale machine learning
references.venue: Siam Review,
references.authors: Bob Carpenter Andrew Gelman Matthew D Hoffman Daniel Lee Ben Goodrich Michael Betancourt Marcus Brubaker Jiqiang Guo Peter Li Allen Riddell
references.title: Stan: A probabilistic programming language
references.venue: Journal of statistical software,
references.authors: Kyung Hyun Cho Tapani Raiko Alexander Ilin
references.title: Gaussian-bernoulli deep boltzmann machine
references.venue: In Neural Networks (IJCNN), The 2013 International Joint Conference on,
references.authors: KyungHyun Cho Alexander Ilin Tapani Raiko
references.title: Improved learning of gaussian-bernoulli restricted boltzmann machines
references.venue: In International conference on artificial neural networks,
references.authors: Zihang Dai Amjad Almahairi Philip Bachman Eduard Hovy Aaron Courville
references.title: Calibrating energy-based generative adversarial networks
references.venue: arXiv preprint arXiv:1702.01691,
references.authors: Yann N Dauphin Razvan Pascanu Caglar Gulcehre Kyunghyun Cho Surya Ganguli Yoshua Bengio
references.title: Identifying and attacking the saddle point problem in high-dimensional non-convex optimization
references.venue: In Advances in neural information processing systems,
references.authors: Jeff Donahue Philipp Krähenbühl Trevor Darrell
references.title: Adversarial feature learning
references.venue: arXiv preprint arXiv:1605.09782,
references.authors: Vincent Dumoulin Ishmael Belghazi Ben Poole Olivier Mastropietro Alex Lamb Martin Arjovsky Aaron Courville
references.title: Adversarially learned inference
references.venue: arXiv preprint arXiv:1606.00704,
references.authors: Chelsea Finn Paul Christiano Pieter Abbeel Sergey Levine
references.title: A connection between generative adversarial networks, inverse reinforcement learning, and energy-based models
references.venue: arXiv preprint arXiv:1611.03852,
references.authors: Ian Goodfellow Jean Pouget-Abadie Mehdi Mirza Bing Xu David Warde-Farley Sherjil Ozair Aaron Courville Yoshua Bengio
references.title: Generative adversarial nets
references.venue: In Advances in neural information processing systems,
references.authors: Martin Heusel Hubert Ramsauer Thomas Unterthiner Bernhard Nessler Sepp Hochreiter
references.title: Gans trained by a two time-scale update rule converge to a local nash equilibrium
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Geoffrey E Hinton
references.title: Training products of experts by minimizing contrastive divergence
references.venue: Neural computation,
references.authors: Geoffrey E Hinton Ruslan R Salakhutdinov
references.title: Reducing the dimensionality of data with neural networks
references.authors: Eric Jang Shixiang Gu Ben Poole
references.title: Categorical reparameterization with gumbel-softmax
references.venue: arXiv preprint arXiv:1611.01144,
references.authors: Michael I Jordan Zoubin Ghahramani Tommi S Jaakkola Lawrence K Saul
references.title: An introduction to variational methods for graphical models
references.venue: Machine learning,
references.authors: Taesup Kim Yoshua Bengio
references.title: Deep directed generative models with energy-based probability estimation
references.venue: arXiv preprint arXiv:1606.03439,
references.authors: Diederik P Kingma Jimmy Ba
references.title: Adam: A method for stochastic optimization
references.venue: arXiv preprint arXiv:1412.6980,
references.authors: Diederik P Kingma Max Welling
references.title: Auto-encoding variational bayes
references.venue: arXiv preprint arXiv:1312.6114,
references.authors: Philipp Krähenbühl Vladlen Koltun
references.title: Efficient inference in fully connected crfs with gaussian edge potentials
references.venue: In Advances in neural information processing systems,
references.authors: Volodymyr Kuleshov Stefano Ermon
references.title: Neural variational inference and learning in undirected graphical models
references.venue: In Advances in Neural Information Processing Systems,
references.authors: John Lafferty Andrew McCallum Fernando CN Pereira
references.title: Conditional random fields: Probabilistic models for segmenting and labeling sequence data
references.authors: Hugo Larochelle Iain Murray
references.title: The neural autoregressive distribution estimator
references.venue: In The Proceedings of the 14th International Conference on Artificial Intelligence and Statistics, volume 15 of JMLR: W&CP,
references.authors: Qiang Liu Dilin Wang
references.title: Learning deep energy models: Contrastive divergence vs. amortized mle
references.venue: arXiv preprint arXiv:1707.00797,
references.authors: Chris J Maddison Andriy Mnih Yee Whye Teh
references.title: The concrete distribution: A continuous relaxation of discrete random variables
references.venue: arXiv preprint arXiv:1611.00712,
references.authors: Xiao-Li Meng Wing Hung Wong
references.title: Simulating ratios of normalizing constants via a simple identity: a theoretical exploration
references.venue: Statistica Sinica,
references.authors: Andriy Mnih Karol Gregor
references.title: Neural variational inference and learning in belief networks
references.venue: arXiv preprint arXiv:1402.0030,
references.authors: Radford M Neal
references.title: Probabilistic inference using markov chain monte carlo methods
references.authors: Radford M Neal
references.title: Annealed importance sampling
references.venue: Statistics and computing,
references.authors: Jiquan Ngiam Zhenghao Chen Pang W Koh Andrew Y Ng
references.title: Learning deep energy models
references.venue: In Proceedings of the 28th international conference on machine learning
references.authors: Rajesh Ranganath Sean Gerrish David Blei
references.title: Black box variational inference
references.venue: In Artificial Intelligence and Statistics,
references.authors: Sashank J Reddi Ahmed Hefny Suvrit Sra Barnabas Poczos Alex Smola
references.title: Stochastic variance reduction for nonconvex optimization
references.venue: In International conference on machine learning,
references.authors: Danilo Jimenez Rezende Shakir Mohamed Daan Wierstra
references.title: Stochastic backpropagation and approximate inference in deep generative models
references.venue: arXiv preprint arXiv:1401.4082,
references.authors: Carsten Rother Vladimir Kolmogorov Victor Lempitsky Martin Szummer
references.title: Optimizing binary mrfs via extended roof duality
references.venue: IEEE Conference on Computer Vision and Pattern Recognition,
references.authors: Ruslan Salakhutdinov Geoffrey Hinton
references.title: Deep Boltzmann machines
references.venue: In Proceedings of the twelfth international conference on artificial intelligence and statistics,
references.authors: Ruslan Salakhutdinov Hugo Larochelle
references.title: Efficient learning of deep Boltzmann machines
references.venue: In Proceedings of the thirteenth international conference on artificial intelligence and statistics,
references.authors: Ruslan Salakhutdinov Iain Murray
references.title: On the quantitative analysis of deep belief networks
references.venue: In Proceedings of the 25th international conference on Machine learning,
references.authors: John Schulman Nicolas Heess Theophane Weber Pieter Abbeel
references.title: Gradient estimation using stochastic computation graphs
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Tijmen Tieleman
references.title: Training restricted boltzmann machines using approximations to the likelihood gradient
references.venue: In Proceedings of the 25th international conference on Machine learning,
references.authors: Martin J Wainwright Michael I Jordan
references.title: Log-determinant relaxation for approximate inference in discrete markov random fields
references.venue: IEEE transactions on signal processing,
references.authors: Martin J Wainwright Tommi S Jaakkola Alan S Willsky
references.title: A new class of upper bounds on the log partition function
references.venue: IEEE Transactions on Information Theory,
references.authors: Xiao Wang Shiqian Ma Donald Goldfarb Wei Liu
references.title: Stochastic quasi-newton methods for nonconvex stochastic optimization
references.venue: SIAM Journal on Optimization,
references.authors: Max Welling Geoffrey E Hinton
references.title: A new learning algorithm for mean field boltzmann machines
references.venue: In International Conference on Artificial Neural Networks,
references.authors: Max Welling Charles A Sutton
references.title: Learning in markov random fields with contrastive free energies
references.venue: In AISTATS. Citeseer,
references.authors: John Winn Christopher M Bishop
references.title: Variational message passing
references.venue: Journal of Machine Learning Research,
references.authors: Shuangfei Zhai Yu Cheng Rogerio Feris Zhongfei Zhang
references.title: Generative adversarial networks as variational training of energy based models
references.venue: arXiv preprint arXiv:1611.01799,
references.authors: Heusel
references.title: 2017) propose a two-time scale update rule to train minimax optimization problems with a convergence guarantee even using K1 = 1. AdVIL converges if using the same training method as in (Heusel et al., 2017), which is summarized in Proposition 2. Proposition 2. AdVIL trained with a two-time scale update rule (Heusel et al., 2017) converges to a stationary local Nash equilibrium almost surely under the following assumptions
references.authors: supk ||θk supk ||ψk supk ||φk
references.title: Remark Compared to Theorem 1, Proposition 2 ensures the convergence of AdVIL without assuming q(v, h) is sufficiently close to P (v, h) in each step. However, a two time-sclae update rule (Heusel et al., 2017) is required to satisfy Assumption 2 and extra weight decay terms are needed to satisfy
references.venue: (Heusel et al.,
references.authors: Q P i.e
references.title: We use the ADAM (Kingma & Ba, 2014) optimizer with the learning rate α = 0.0003, the moving average ratios β1 = 0.5 and β2 = 0.999, and the batch size of 500. We use a continuous z and the sigmoid activation
sections.heading: 1 INTRODUCTION
sections.text: Markov random fields (MRFs) find applications in a variety of machine learning areas (Krähenbühl & Koltun, 2011; Salakhutdinov & Larochelle, 2010; Lafferty et al., 2001). In particular, one famous example is conditional random fields (Lafferty et al., 2001), a conditional version of MRFs that was developed to address the limitations (e.g., local dependency and label bias) of directed models for sequential data (e.g., hidden Markov models and other discriminative Markov models based on directed graphical models). However, the inference and learning of general MRFs are challenging due to the presence of a global normalizing factor, i.e. partition function, especially when latent variables are present. Extensive efforts have been devoted to developing approximate methods. On one hand, sample-based methods (Neal, 1993) and variational approaches (Jordan et al., 1999; Welling & Sutton, 2005; Salakhutdinov & Larochelle, 2010) are proposed to infer the latent variables. On the other hand, extensive work (Meng & Wong, 1996; Neal, 2001; Hinton, 2002; Tieleman, 2008; Wainwright et al., 2005; Wainwright & Jordan, 2006) has been done to estimate the partition function. Among these methods, contrastive divergence (Hinton, 2002) is proven effective in certain types of models.
  Most of the existing methods highly depend on the model structure and require model-specific analysis in new applications, which makes it important to develop black-box inference and learning methods. Previous work (Ranganath et al., 2014; Schulman et al., 2015) shows the ability to automatically infer the latent variables and obtain gradient estimate in directed models. However, there is no black-box learning method for undirected models except the recent work of NVIL (Kuleshov & Ermon, 2017).
  NVIL introduces a variational distribution and derives an upper bound of the partition function in a general MRF, in the same spirit as amortized inference (Kingma & Welling, 2013; Rezende et al., 2014; Mnih & Gregor, 2014) for directed models. NVIL has several advantages over existing methods, including the ability of black-box learning, tracking the partition function during training and getting approximate samples efficiently during testing. However, NVIL also comes with two disadvantages: (1) it leaves the inference problem of MRFs unsolved1 and only trains simple MRFs with tractable
  ∗Dept. of Comp. Sci. & Tech., BNRist Center, Institute for AI, THBI Lab, Tsinghua University, Beijing, 100084, China †University of Amsterdam, and the Canadian Institute for Advanced Research (CIFAR). 1NVIL (Kuleshov & Ermon, 2017) presents a hybrid model. The “inference” in the title refers to directed part but not for an MRF.
  posteriors, and (2) the upper bound of the partition function can be underestimated (Kuleshov & Ermon, 2017), resulting in sub-optimal solutions on high-dimensional data.
  We propose Adversarial Variational Inference and Learning (AdVIL) to relieve some headache of learning an MRF model. AdVIL is a black-box inference and learning method that partly solves the two problems of NVIL and retains the advantages of NVIL at the same time. First, AdVIL introduces a variational encoder to infer the latent variables, which provides an upper bound of the free energy. Second, AdVIL introduces a variational decoder for the MRF, which provides a lower bound of the log partition function. The two variational distributions provide an estimate of the negative log-likelihood of the MRF. On one hand, the estimate is in an intuitive form of an approximate contrastive free energy, which is expressed in terms of the expected energy and the (conditional) entropy of the corresponding variational distribution. On the other hand, similar to GAN (Goodfellow et al., 2014), the estimate is a minimax optimization problem, which is solved by stochastic gradient descent (SGD) in an alternating manner. Theoretically, our algorithm is convergent if the variational decoder approximates the model well. This motivates us to introduce an auxiliary variable to enhance the flexibility of the variational decoder, whose entropy is approximated by the third variational trick.
  We evaluate AdVIL in various undirected generative models, including restricted Boltzmann machines (RBM) (Ackley et al., 1985), deep Boltzmann machines (DBM) (Salakhutdinov & Hinton, 2009), and Gaussian restricted Boltzmann machines (GRBM) (Hinton & Salakhutdinov, 2006), on several real datasets. We empirically demonstrate that (1) compared to the black-box NVIL (Kuleshov & Ermon, 2017) method, AdVIL provides a tighter estimate of the log partition function and achieves much better log-likelihood results; and (2) compared to contrastive divergence based methods (Hinton, 2002; Welling & Sutton, 2005), AdVIL can deal with a broader family of MRFs without model-specific analysis and obtain better results when the model structure gets complex as in DBM.
sections.heading: 2 BACKGROUND
sections.text: We consider a general case where the model consists of both visible variables v and latent variables h. An MRF defines the joint distribution over v and h as P (v, h) = e −E(v,h)
  Z , where E denotes the associated energy function that assigns a scalar value for a given configuration of (v, h) and Z is the partition function such that Z = ∫ v,h e−E(v,h)dvdh.
  Let PD(v) denote the empirical distribution of the training data. Minimizing the negative loglikelihood (NLL) of an MRF is a commonly chosen learning criterion and it is given by:
  L(θ) :=− EPD(v) [ log ∫ h e−E(v,h) Z dh ] , (1)
  where θ denotes the trainable parameters in E . Further, the gradient of θ is:
  ∇θL(θ) = EPD(v) [∇θF(v)]− EP (v) [∇θF(v)] , (2) where F(v) = − log ∫ h e−E(v,h)dh denotes the free energy and the gradient in Eqn. (2) is the difference of the free energy in two phases. In the first positive phase, the expectation of the free energy under the data distribution is decreased. In the second negative phase, the expectation of the free energy under the model distribution is increased.
  Unfortunately, both the NLL in Eqn. (1) and its gradient in Eqn. (2) are intractable in general for two reasons. First, the integral of the latent variables in Eqn. (1) or equivalently the computation of the free energy in Eqn. (2) is intractable. Second, the computation of the partition function in Eqn. (1) or equivalently the negative phase in Eqn. (2) is intractable.
  Variational inference. Extensive work introduces deterministic approximations for the intractability of inference, including the mean-field approximation (Welling & Hinton, 2002; Salakhutdinov & Hinton, 2009), the Kikuchi and Bethe approximations (Welling & Sutton, 2005) and the recognition model approach (Salakhutdinov & Larochelle, 2010). In this line of work, the intractability of the partition function is addressed using Monte Carlo based methods.
  Contrastive free energy. Contrastive divergence (CD) (Hinton, 2002) addresses the intractability of the partition function by approximating the negative phase in Eqn. (2) as follows:
  ∇θL(θ) = EPD(v) [∇θF(v)]− EPCD(v) [∇θF(v)] , (3)
  where PCD(v) denotes the empirical distribution obtained by starting from a data point and running several steps of Gibbs sampling according to the model distribution and the free energy F(v) is assumed to be tractable. Existing methods (Welling & Hinton, 2002; Welling & Sutton, 2005) approximate F(v) using certain function G(v) and the gradient of θ is:
  ∇θL(θ) ≈ EPD(v) [∇θG(v)]− EPCD(v) [∇θG(v)] . (4)
  Although these generalized methods exist, it is nontrivial to extend CD-based methods to general MRFs because the Gibbs sampling procedure is highly dependent on the model structure.
  Black-box learning. The recent work of NVIL (Kuleshov & Ermon, 2017) addresses the intractability of the partition function in a black-box manner via a variational upper bound of the partition function:
  Eq(v)
  [ P̃ (v)2
  q(v)2
  ] ≥ Z2, (5)
  where P̃ (v) = e−F(v) is the unnormalized marginal distribution on v and q(v) is a neural variational distribution. As a black-box learning method, NVIL potentially allows application to broader model families and improves the capabilities of probabilistic programming systems (Carpenter et al., 2017). Though promising, NVIL leaves the intractability of inference in an MRF unsolved, and the bound in Eqn. (5) is of high variance and is easily underestimated (Kuleshov & Ermon, 2017).
sections.heading: 3 METHOD
sections.text: As stated above, the black-box inference and learning of MRFs are still largely open. In this paper, we make a step towards solving the problems by a new variational approach. For simplicity, we focus on the resulting objective function in this section. See Appendix A for detailed derivation.
sections.heading: 3.1 ADVERSARIAL VARIATIONAL INFERENCE AND LEARNING
sections.text: First, we rewrite the NLL of the MRF (See an illustration in Fig. 1) as follows:
  L(θ) = −EPD(v) [−F(v)] + logZ, (6)
  where the negative free energy and the log partition function are in the form of a logarithm of an integral. Naturally, we can apply the variational trick (Jordan et al., 1999) twice and approximate the two terms individually. Due to the presence of the minus before the first term in Eqn. (6), the two variational tricks bound the two parts of the NLL in the opposite directions, detailed as below.
  Formally, on one hand, we introduce an approximate posterior for the latent variables Q(h|v), which is parameterized as a neural variational encoder (See an illustration in Fig. 1), to address the intractability of inference as follows:
  L(θ) ≤ EPD(v)Q(h|v) [E(v, h) + logQ(h|v)] + logZ := L1(θ, φ), (7)
  where φ denotes the trainable parameters in Q(h|v). The upper bound is derived via applying the Jensen inequality and the equality holds if and only if Q(h|v) = P (h|v) for all v. In the bound, the first term is the expected energy, which encourages Q(h|v) to infer latent variables that have low values of the energy function E(v, h), or equivalently high probabilities of P (v, h). The second term corresponds to the negative conditional entropy of Q(h|v), which increases the uncertainty of Q(h|v). In the paper, we denote the conditional entropy of Q(h|v) asH(Q) := −EPD(v)Q(h|v)[logQ(h|v)].
  On the other hand, we introduce an approximate sampler q(v, h), which is parameterized by a neural variational decoder (See Fig. 1), to address the intractability of the partition function as follows:
  L1(θ, φ)≥EPD(v)Q(h|v) energy term︷ ︸︸ ︷E(v, h) + entropy term︷ ︸︸ ︷logQ(h|v)︸ ︷︷ ︸ Positive Phase −Eq(v,h) energy term︷ ︸︸ ︷E(v, h) + entropy term︷ ︸︸ ︷log q(v, h)︸ ︷︷ ︸ Negative Phase := L2(θ, φ, ψ), (8) where ψ denotes the trainable parameters in q(v, h). The lower bound is derived via applying the Jensen inequality as well, and the equality holds if and only if q(v, h) = P (v, h). It can be seen that the lower bound given by q(v, h) consists of the entropy (denoted asH(q)) and energy terms, which is similar to the upper bound in Eqn. (7), and the overall objective is in the form of approximate contrastive free energy (Hinton, 2002; Welling & Sutton, 2005). Because the double variational trick bounds the NLL in opposite directions as above, we have a minimax optimization problem:
  min θ min φ max ψ L2(θ, φ, ψ). (9)
  The minimax formulation has been investigated in GAN (Goodfellow et al., 2014) and it is interpreted as an adversarial game between two networks. We name our framework adversarial variational inference and learning (AdVIL) following the well-established literature.
  Note that L2(θ, φ, ψ) is neither an upper bound, nor a lower bound of L(θ) due to the double variational trick. However, we argue that solving the optimization problem in Eqn. (9) is reasonable because (1) it is equivalent to optimizing L(θ) under the nonparametric assumption, which is similar to GAN (Goodfellow et al., 2014); and (2) it converges to a stationary point of L1(θ, φ), which is an upper bound of L(θ), under a weaker assumption, as stated in the following theoretical analysis.
sections.heading: 3.2 THEORETICAL ANALYSIS OF ADVIL
sections.text: In this section, we present our main theoretical results and the proofs can be found in Appendix C. Firstly, similarly to GAN (Goodfellow et al., 2014), we can prove that L2 is a tight estimate of L under the nonparametric assumption, which is summarized in Proposition 1 in Appendix C.1. However, the nonparametric assumption does not tolerate any approximation error between P (v, h) and q(v, h) during training and no guarantee can be obtained in finite steps. To this end, we establish a convergence theorem based on a weaker assumption that allows non-zero approximation error before convergence. A key insight is that the angle between ∂L2(θ,φ,ψ)∂θ and ∂L1(θ,φ) ∂θ is positive if q(v, h) approximates P (v, h) well, as stated in the following Lemma 1.
  Lemma 1. For any (θ, φ), there exists a symmetric positive definite matrix H such that ∂L2(θ,φ,ψ)∂θ = H ∂L1(θ,φ)∂θ under the assumption: || ∑ v,h δ(v, h) ∂E(v,h) ∂θ ||2 < || ∂L1(θ,φ) ∂θ ||2 if || ∂L1(θ,φ) ∂θ ||2 > 0 and
  || ∑ v,h δ(v, h) ∂E(v,h) ∂θ ||2 = 0 if || ∂L1(θ,φ) ∂θ ||2 = 0, where δ(v, h) = q(v, h)− P (v, h).
  Based on Lemma 1 and other commonly used assumptions in the analysis of stochastic optimization (Bottou et al., 2018), AdVIL converges to a stationary point of L1(θ, φ), as stated in Theorem 1. Theorem 1. Solving the optimization problem in Eqn. (9) using stochastic gradient descent, then (θ, φ) converges to a stationary point of L1(θ, φ) under the assumptions in the general stochastic optimization (Bottou et al., 2018) and that the condition of Lemma 1 holds in each step.
  Please see Appendix C.2 for a detailed and formal version of Theorem 1. Compared to Proposition 1 and the analysis in GAN (Goodfellow et al., 2014), Theorem 1 has a weaker statement that AdVIL converges to a stationary point of the negative evidence lower bound (i.e., L1) instead of L. Nevertheless, we argue that converging to L1 is sufficiently good for variational approaches in general. Besides, Theorem 1 states that AdVIL can at least decrease L1 in expectation if the assumption holds in finite steps. Indeed, we empirically justify Theorem 1, as detailed in Appendix E.1. Theorem 1 also provides insights for the implementation of AdVIL. Indeed, its assumption motivates us to use a sufficiently powerful q(v, h) with neural networks and auxiliary variables, and update q(v, h) multiple times per update of P (v, h), as detailed in Sec. 3.3 and Sec. 5.1 respectively.
sections.heading: 3.3 SPECIFYING THE VARIATIONAL DISTRIBUTIONS
sections.text: To efficiently get samples, both variational distributions are directed models. We use a directed neural network that maps v to h as the variational encoder Q(h|v) (Kingma & Welling, 2013).
  As for the variational decoder, we first factorize it as the product of a prior over h and a conditional distribution, namely q(v, h) = q(v|h)q(h). It is nontrivial to specify the prior q(h) because the marginal distribution of h in the MRF, i.e. P (h), can be correlated across different dimensions. Consequently, a simple q(h) is not flexible enough to track P (h) and can violate the condition of Lemma 1. To this end, we introduce an auxiliary variable z, which can be discrete or continuous, on top of h and define q(v, h) = ∫ z q(z)q(h|z)q(v|h)dz.2 (See an illustration in Fig. 1.) However, the entropy term of q(v, h) is intractable because we need to integrate out the auxiliary variable z. Therefore, we introduce the third variational distribution r(z|h) to approximate the entropy of q(v, h). As in Eqn. (7), applying the standard variational trick gives an upper bound:
  −Eq(v,h) log q(v, h) ≤ −Eq(v,h) log q(v|h)− Eq(h)r(z|h) log [ q(h, z)
  r(z|h)
  ] , (10)
  which is unsatisfactory because the estimate is minimized w.r.t r(z|h) while maximized w.r.t q(v, h). Instead, after some transformations (See details in Appendix A) we get a lower bound as follows:
  −Eq(v,h) log q(v, h) ≥ −Eq(v,h) log q(v|h)− Eq(h,z) log [ q(h, z)
  r(z|h)
  ] . (11)
  The equality holds if and only if r(z|h) = q(z|h) for all h. The difference between the two bounds is subtle: the last expectation in Eqn. (10) is over q(h)r(z|h) but that in Eqn. (11) is over q(h, z). Here, a lower bound is preferable because the estimate is maximized with respect to both r(z|h) and q(v, h) and we can train them simultaneously. For simplicity, we absorb the trainable parameters of r(z|h) into ψ. Note that after introducing z and r(z|h), we can still obtain a convergence theorem of AdVIL under the conditions that r(z|h) approximates q(z|h) well and q(v, h) = ∫ q(v, h, z)dz is sufficiently close to P (v, h) in every step, together the assumptions in general stochastic optimization.
  Following GAN (Goodfellow et al., 2014), we optimize θ, φ and ψ jointly using stochastic gradient descent (SGD) in an alternating manner. The partial derivatives of φ and ψ are estimated via the reparameterization trick (Kingma & Welling, 2013) for the continuous variables and the GumbelSoftmax trick (Jang et al., 2016; Maddison et al., 2016) for the discrete variables. See Algorithm 1 in Appendix B for the whole training procedure. Note that ψ is updated K1 > 1 times per update of θ.
sections.heading: 4 RELATED WORK
sections.text: Existing traditional methods (Neal, 2001; Hinton, 2002; Winn & Bishop, 2005; Wainwright & Jordan, 2006; Rother et al., 2007) can be used to estimate the log partition function but are nontrivial to be extended to learn general MRFs. Some methods (Winn & Bishop, 2005; Neal, 2001) require an expensive inference procedure for each update of the model and others (Hinton, 2002; Rother et al., 2007) cannot be directly applied to general cases (e.g., DBM). Among these methods, contrastive divergence (CD) (Hinton, 2002) is proven effective in certain types of models and it is closely related to AdVIL. Indeed, the partial derivative of θ in AdVIL is:
  ∂L2(θ, φ, ψ) ∂θ = EPD(v)Q(h|v) [ ∂ ∂θ E(v, h) ] − Eq(v,h) [ ∂ ∂θ E(v, h) ] , (12)
  which also involves a positive phase and a negative phase naturally and is quite similar to Eqn. (3). However, notably, the two phases average over the (v, h) pairs and only require the knowledge of the energy function without any further assumption of the model in AdVIL. Therefore, AdVIL is more suitable to general MRFs than CD (See empirical evidence in Sec. 5.3).
  In the context of black-box learning in MRFs, AdVIL competes directly with NVIL (Kuleshov & Ermon, 2017). It seems that the upper bound in Eqn. (5) is suitable for optimization because P and q share the same training direction. However, the bound holds only if the support of P̃ is a subset of the support of q. Further, the Monte Carlo estimate of the upper bound is of high variance. Therefore, the bound of NVIL can be easily underestimated, which results in sub-optimal solutions (Kuleshov & Ermon, 2017). In contrast, though AdVIL arrives at a minimax optimization problem, the estimate of Eqn. (8) is tighter and of lower variance. We empirically verify this argument (See Fig. 3) and systematically compare the two methods (See Tab. 1) in Sec.5.4.
  2An alternative way is to use an autoregressive model as q(h). See results and analysis in Appendix E.3.
  Apart from the work on approximate inference and learning in MRFs as mentioned above, AdVIL is also related to some directed models. Kim & Bengio (2016) jointly trains a deep energy model (Ngiam et al., 2011) and a directed generative model by minimizing the KL-divergence between them. Similar ideas have been highlighted in (Finn et al., 2016; Zhai et al., 2016; Dai et al., 2017; Liu & Wang, 2017). In comparison, firstly, AdVIL obtains the objective function in a unified perspective on the black-box inference and learning in general MRFs. Note that dealing with latent variables in MRFs is nontrivial (Kim & Bengio, 2016) and therefore existing work focuses on fully observable models. Secondly, AdVIL uses a sophisticated decoder with auxiliary variables to handle the latent variables and derives a principled variational approximation of the entropy term instead of the heuristics (Kim & Bengio, 2016; Zhai et al., 2016). Lastly, the convergence of AdVIL is formally characterized by Theorem 1 while the effect of the approximation error in inference is not well understood in existing methods. Adversarially learned inference (ALI) (Donahue et al., 2016; Dumoulin et al., 2016) is also formulated as a minimax optimization problem but focuses on directed models.
sections.heading: 5 EXPERIMENTS
sections.text: In this section, we evaluate AdVIL in restricted Boltzmann machines (RBM) (Ackley et al., 1985), deep Boltzmann machines (DBM) (Salakhutdinov & Hinton, 2009) and Gaussian restricted Boltzmann machines (GRBM) (Hinton & Salakhutdinov, 2006) on the Digits dataset, the UCI binary databases (Dheeru & Karra, 2017) and the Frey faces datasets (See detailed settings in Appendix D and the source code3). We compare AdVIL with strong baseline methods systematically and show the promise of AdVIL to learn a broad family of models effectively as a black-box method.
sections.heading: 5.1 EMPIRICAL ANALYSIS OF ADVIL
sections.text: We present a detailed analysis of AdVIL in RBM, whose energy function is defined as E(v, h) = −b>v−v>Wh−c>h. The conditional distributions of an RBM are tractable, but we still treat P (h|v) as unknown and train AdVIL in a fully black-box manner. The analysis is performed on the Digits dataset and we augment the data of five times by shifting the digits following the protocol in (Kuleshov & Ermon, 2017). The dimensions of v, h and z are 64, 15 and 10, respectively. Therefore, the log partition function of the RBM and the entropy of the decoder can be computed by brute force.
  Firstly, we empirically validate AdVIL in Fig. 2. Specifically, Panel (a) shows that the variational encoder Q(h|v) provides a tight upper bound of the free energy after 2,000 iterations. Panel (b) demonstrates that the variational distribution r(z|h) estimate the entropy of q(v, h) accurately. Panel (c) shows that q(v, h) can successfully track the log partition function after 5,000 iterations. Panel (d) presents that the RBM loss balances well between the negative phase and positive phase, and the model converges gradually. See Appendix E.1 for an empirical test of the condition in Lemma 1.
  Secondly, we empirically show that both P and q can generate data samples in Appendix E.2.
  Lastly, we analyze the sensitivity of K1. Theoretically, enlarging K1 will make q(v, h) and P (v, h) to be close and then help the convergence according to Theorem 1. As shown in Fig. 3 (a), a larger K1 at least won’t hurt the convergence, which agrees with Theorem 1. Though K1 = 15 is sufficient on the Digits dataset, we use K1 = 100 as a default setting for AdVIL on larger datasets.
  3See the source code in https://anonymous.4open.science/r/8c779fbc-6394-40c7-8273-e52504814703/.
sections.heading: 5.2 RBM RESULTS
sections.text: To the best of our knowledge, NVIL (Kuleshov & Ermon, 2017) is the only existing black-box learning method for MRFs and hence it is the most direct competitor of AdVIL. In this section, we provide a systematic comparison and analysis of these two methods in terms of the log-likelihood results on the UCI databases (Dheeru & Karra, 2017).
  For a fair comparison, we use the widely-adopted anneal importance sampling (AIS) (Salakhutdinov & Murray, 2008) metric for quantitative evaluation. Besides, we carefully perform grid search over the default settings of NVIL (Kuleshov & Ermon, 2017) and our settings based on their code, and choose the best configuration including K1 = 100 (See details in Appendix D). We directly compare with the best version of NVIL in Tab. 1. It can be seen that AdVIL consistently outperforms NVIL on all datasets, which demonstrate the effectiveness of AdVIL. Besides, the time complexity of AdVIL is comparable to that of NVIL with the same hyperparameters.
  We compare the learning curves of NVIL and AdVIL on the Mushroom dataset. As shown in Fig. 3 (b), the upper bound of NVIL is underestimated after 4,000 iterations and then the model can get worse or even diverge. In contrast, as shown in Fig. 3 (c), the lower bound of AdVIL is consistently valid. Besides, the estimate of NVIL is looser and of higher variance than that of AdVIL. The results agree with our analysis in Sec. 4 and explain why AdVIL significantly outperforms NVIL. Further, as shown in Fig. 3 (d), AdVIL is comparable to CD-10 and persistent contrastive divergence (PCD) (Tieleman, 2008), which leverage the tractability of the conditional distributions in an RBM.
sections.heading: 5.3 DBM RESULTS
sections.text: We would like to demonstrate that AdVIL has the ability to deal with highly intractable models such as a DBM conveniently and effectively, compared to standard CD-based methods (Hinton, 2002; Welling & Hinton, 2002; Welling & Sutton, 2005) and NVIL (Kuleshov & Ermon, 2017).
  DBM (Salakhutdinov & Hinton, 2009) is a powerful family of deep models that stack multiple RBMs together. The energy function of a two-layer DBM is defined as E(v, h1, h2) = −b>v − v>W1h1 − c>1 h1 − h>1 W2h2 − c>2 h2. Learning a DBM is challenging because P (h1, h2|v) is not tractable and CD (Hinton, 2002) is not applicable. Inspired by (Welling & Hinton, 2002; Welling & Sutton, 2005), we construct a variational CD (VCD) baseline by employing the same variational encoder Q(h1, h2|v) as in AdVIL. The free energy is approximated by the same upper bound as in Eqn. (7), which is minimized with respect to the parameters in Q(h1, h2|v). The gradient of the parameters in the DBM is given by Eqn. (4), where the Gibbs sampling procedure is approximated by h1 ∼ Q(h1|v) and v ∼ P (v|h1). Note that AdVIL can be directly applied to this case. As for
  the time complexity, the training speed of AdVIL is around ten times slower than that of VCD in our implementation. However, the approximate inference and sampling procedure of AdVIL is very efficient thanks to the directed variational distributions.
  The log-likelihood results on the UCI databases are shown in Tab. 2. It can be seen that AdVIL has a better overall performance even trained in a black-box manner, which shows the promise of AdVIL. See Appendix E.4 for learning curves and a detailed analysis of the results.
  We also extend NVIL by using the same Q(h1, h2|v) and q(v, h1, h2) as AdVIL. However, NVIL diverges after 300 iterations and gets bad AIS results (e.g., less than −40 on Digits) in our implementation. A potential reason is that the upper bound given by q in NVIL can be underestimated if q is high-dimensional, as analyzed in Sec. 4 and Fig. 3. Note that q(v, h1, h2) in DBM involves latent variables and has a higher dimension (e.g. 164 on the Digits dataset) than q(v) in RBM (e.g. 64 on the Digits dataset). The results again demonstrate the advantages of AdVIL over NVIL.
sections.heading: 5.4 GRBM RESULTS
sections.text: We now show the ability of AdVIL to learn a GRBM on the continuous Frey faces dataset. The energy function of a GRBM is E(v, h) = 12σ2 ||v− b|| 2− c>h− 1σv >Wh, where σ is the standard deviation of the Gaussian likelihood and is set as 1 manually. We standardize the data by subtracting the mean and dividing by the standard deviation. The dimensions of h and z are 200 and 50, respectively.
  Though a GRBM is more sensitive to the hyperparameters and hence harder to train than an RBM (Cho et al., 2011; 2013), AdVIL can successfully capture the underlying data distribution using the default hyperparameters (See Appendix D). As shown in Fig. 4, the samples from both the GRBM (via Gibbs sampling after 100,000 burn-in steps) and the decoder are meaningful faces. Besides, the filters of the GRBM outline diverse prototypes of faces, which accords with our expectation.
  In summary, the results of the three models together demonstrate that AdVIL can learn a broad family of models conveniently and effectively in a fully black-box manner.
sections.heading: 6 CONCLUSION AND DISCUSSION
sections.text: A novel black-box learning and inference method for undirected graphical models, called adversarial variational inference and learning (AdVIL), is proposed. The key to AdVIL is a double variational trick that approximates the negative free energy and the log partition function separately. A formal convergence theorem, which provides insights for implementation, is established for AdVIL. Empirical results show that AdVIL can deal with a broad family of MRFs in a fully black-box manner and outperforms both the standard contrastive divergence method and the black-box NVIL algorithm.
  Though AdVIL shows promising results, we emphasize that the black-box learning and inference of the MRFs are far from completely solved, especially on high-dimensional data. The two intractability
  problems of MRFs are distinct since the posterior of the latent variables is local in terms of v but the partition function is global by integrating out v. The additional integral makes estimating the partition function much more challenging. In AdVIL, simply increasing the number of updates of the decoder to obtain a tighter estimate of the partition function on high-dimensional data can be expensive. A potential future work to avoid the problem is adopting recent advances on non-convex optimization (Dauphin et al., 2014; Reddi et al., 2016; Wang et al., 2017) to accelerate the inner loop optimization. We conjecture that AdVIL is comparable to CD in RBM and superior to VCD in DBM on larger datasets if AdVIL can be trained to nearly converge based on our current results.
sections.heading: ACKNOWLEDGEMENTS
sections.text: This work was supported by the National Key Research and Development Program of China (No. 2017YFA0700904), NSFC Projects (Nos. 61620106010, U19B2034, U1811461), Beijing NSF Project (No. L172037), Beijing Academy of Artificial Intelligence (BAAI), Tsinghua-Huawei Joint Research Program, a grant from Tsinghua Institute for Guo Qiang, Tiangong Institute for Intelligent Computing, the JP Morgan Faculty Research Program and the NVIDIA NVAIL Program with GPU/DGX Acceleration. C. Li was supported by the Chinese postdoctoral innovative talent support program and Shuimu Tsinghua Scholar.
sections.heading: A DERIVATION OF THE OBJECTIVE FUNCTION
sections.text: Here we derive the objective function of AdVIL in detail. Let θ, φ and ψ denote the trainable parameters in the MRF, the variational encoder and the variational decoder, respectively. The first variational trick bounds the free energy as follows:
  L(θ) = −EPD(v) [ log( ∫ h e−E(v,h)dh) ] + logZ
  = −EPD(v) log[ ∫ h Q(h|v)e −E(v,h) Q(h|v) dh] + logZ ≤ EPD(v)Q(h|v) [E(v, h) + logQ(h|v)] + logZ := L1(θ, φ),
  where the bound is derived via applying the Jensen inequality and the equality holds if and only if Q(h|v) = P (h|v) for all v. The second variational trick bounds the log partition function as follows:
  L1(θ, φ) = EPD(v)Q(h|v) [E(v, h) + logQ(h|v)] + log( ∫ v ∫ h e−E(v,h)dvdh)
  = EPD(v)Q(h|v) [E(v, h) + logQ(h|v)] + log( ∫ v ∫ h q(v, h) e−E(v,h) q(v, h) dvdh)
  ≥ EPD(v)Q(h|v) [E(v, h) + logQ(h|v)] + Eq(v,h) [ log( e−E(v,h)
  q(v, h) )
  ]
  = EPD(v)Q(h|v) energy term︷ ︸︸ ︷E(v, h) + entropy term︷ ︸︸ ︷logQ(h|v)︸ ︷︷ ︸ Positive Phase − Eq(v,h) energy term︷ ︸︸ ︷E(v, h) + entropy term︷ ︸︸ ︷log q(v, h)︸ ︷︷ ︸ Negative Phase  := L2(θ, φ, ψ),
  where the bound is also derived via applying the Jensen inequality and the equality holds if and only if q(v, h) = P (v, h).
  To enhance the expressive power of the variational decoder, we introduce an auxiliary variable z and define q(v, h) = ∫ z q(z)q(h|z)q(v|h)dz, which makes the entropy term in the negative phase intractable. To address the problem, we propose the third variational approximation. First, we can decompose the entropy of q(v, h) as −Eq(v,h) log q(v, h) = −Eq(v,h) log q(v|h) − Eq(h) log q(h) and we only need to approximate−Eq(h) log q(h). However, simply applying the standard variational trick as above, we get an upper bound as follows:
  −Eq(h) log q(h) = −Eq(h) log ∫ z q(h, z)dz
  = −Eq(h) log ∫ z r(z|h)q(h, z) r(z|h) dz
  ≤ −Eq(h)r(z|h) log [ q(h, z)
  r(z|h)
  ] ,
  which is not satisfactory because the optimization problem will be minP minQmaxqminr. Instead, we derive a lower bound as follows:
  −Eq(h) log q(h) = −Eq(h) log q(h)− Eq(h,z) log q(z|h) + Eq(h,z) log q(z|h) = −Eq(h,z) log q(h, z) + Eq(h,z) log q(z|h)
  = −Eq(h,z) log [ q(h, z)
  r(z|h)
  ] + DKL(q(z|h)||r(z|h))
  ≥ −Eq(h,z) log [ q(h, z)
  r(z|h)
  ] ,
  where DKL(·||·) denotes the KL-divergence and the equality holds if and only if r(z|h) = q(z|h) for all h. The difference between the two bounds is that the expectation is taken over q(h)r(z|h) in the upper bound while over q(h, z) in the lower bound. Using the lower bound, the optimization problem will be minP minQmaxpmaxr.
sections.heading: B FORMAL TRAINING PROCEDURE
sections.text: The formal training procedure of AdVIL is presented in Algorithm 1.
  Algorithm 1 Adversarial variational inference and learning by stochastic gradient descent 1: Input: Constants K1 and K2, learning rate schemes α and γ, randomly initialized θ, φ and ψ 2: repeat 3: for i = 1, ...,K1 do 4: Sample a batch of (v, h, z) ∼ q(v, h, z) 5: Estimate the objective of q and r according to Eqn. (11) and the negative phase in Eqn. (8) 6: Update ψ to maximize the objective according to α 7: end for 8: for i = 1, ...,K2 do 9: Sample a batch of (v, h) ∼ PD(v)Q(h|v) 10: Estimate the objective of Q according to the positive phase in Eqn. (8) 11: Update φ to minimize the objective according to γ 12: end for 13: Sample a batch of (v, h) ∼ PD(v)Q(h|v) and another batch of (v, h) ∼ q(v, h) 14: Estimate the objective of P according to Eqn. (8) 15: Update θ to minimize the objective according to γ 16: until Convergence or reaching certain threshold
sections.heading: C DETAILED THEORETICAL ANALYSIS
sections.text: For simplicity, we consider discrete v and h (e.g., in an RBM) and the analysis can be extended to the continuous cases. We assume v ∈ {0, 1}dv and h ∈ {0, 1}dh , where dv and dh are the dimensions of the visible and latent variables respectively.
  C.1 ANALYSIS IN THE NONPARAMETRIC CASE
  We first analyze the nonparametric case in Proposition 1 as follows. Proposition 1. For any P (v, h) = exp(−E(v, h))/Z , L2(θ, φ, ψ) is a tight estimate of the negative log-likelihood of P (v), under the following assumptions
  1. Q(h|v) and q(v, h) are nonparametric.
  2. The inner optimization over Q(h|v) and q(v, h) can get their optima.
  Proof. Given P (v, h), i.e., E(v, h), to find q∗(v, h), we optimize L2 over {q(v, h)|v ∈ {0, 1}dv , h ∈ {0, 1}dh} (we will use a shortcut {q(v, h)} for simplicity). The optimization problem is equivalent to:
  min {q(v,h)} ∑ v,h q(v, h) [E(v, h) + log q(v, h)]
  subject to: ∑ v,h q(v, h) = 1,
  q(v, h) ≥ 0,∀v, h. Note that the objective function is convex since its Hessian matrix is positive semi-definite. Besides, the constraints are linear. Therefore, it is a convex optimization problem. Further, we can verify that the Slater’s condition (Boyd & Vandenberghe, 2004) holds when q is uniform and then the strong duality holds. Then, we can use the KKT conditions to solve the optimization problem.
  The Lagrangian G({q(v, h)}, λ, {µ(v, h)}) is:∑ v,h q(v, h) [E(v, h) + log q(v, h)] + λ( ∑ v,h q(v, h)− 1) + ∑ v,h µ(v, h)q(v, h),
  where λ and {µ(v, h)} are the associated Lagrange multipliers. To satisfy the stationarity, we take gradients with respect to q(v, h) for all (v, h) and get:
  [E(v, h) + log q∗(v, h) + 1] + λ+ µ(v, h) = 0,
  which implies q∗(v, h) = exp(−E(v, h)− (1 + λ+ µ(v, h))).
  According to the complementary slackness, we have µ(v, h)q∗(v, h) = 0,∀v, h,
  which implies µ(v, h) = 0,∀v, h, since q∗(v, h) > 0,∀v, h. To satisfy the primal equality constraint, we have∑
  v,h q∗(v, h) = ∑ v,h exp(−E(v, h)− (1 + λ)) = 1,
  which implies
  q∗(v, h) = exp(−E(v, h))∑
  v′,h′ exp(−E(v′, h′)) = P (v, h),∀v, h.
  To find Q∗(h|v), we optimize L2 over {Q(h|v)|v ∈ {0, 1}dv , h ∈ {0, 1}dh} (we will use a shortcut {Q(h|v)} for simplicity). The optimization problem is equivalent to:
  min {Q(h|v)} ∑ v PD(v) ∑ h Q(h|v) [E(v, h) + logQ(h|v)]
  subject to: ∑ h Q(h|v) = 1,∀v,
  Q(h|v) ≥ 0,∀v, h. Similar to the above procedure, we can get
  Q∗(h|v) = exp(−E(v, h))∑ h′ exp(−E(v, h′)) = P (h|v),∀v, h.
  Under the assumptions that (1) Q(h|v) and q(v, h) are nonparametric, and (2) the inner optimization over ψ and φ can get the optimum, the optimal variational distributions P (v, h) and P (h|v) can be obtained. Plugging them back into L2, we get
  L2 = EPD(v)P (h|v) [E(v, h) + logP (h|v)]− EP (v,h) [E(v, h) + logP (v, h)]
  = EPD(v)P (h|v)
  [ − log
  ∑ h e−E(v,h)
  ] + EP (v,h) [logZ]
  = EPD(v) [F(v)] + logZ = L.
  Remark Similar to Theorem 1 in (Goodfellow et al., 2014), Proposition 1 is under the nonparametric assumption, which is relaxed in our following analysis. Namely, we will consider more practical cases where q(v, h) may not be exactly the same as P (v, h) during training.
  C.2 MAIN CONVERGENCE THEOREM
  For convenience, we summarize the training dynamics of Algorithm 1 with K1 = 1,K2 = 1 and the exact gradients (not the stochastic ones), as follows:
  ψk+1 = ψk + αk ∂L2(θk, φk, ψk)
  ∂ψ ,
  φk+1 = φk − γk ∂L2(θk, φk, ψk+1)
  ∂φ ,
  θk+1 = θk − γk ∂L2(θk, φk, ψk+1)
  ∂θ , (13)
  where k = 1, 2, ... . We will prove that even though we are optimizingL2(θ, φ, ψ), (θk, φk) converges to a stationary point of L1(θ, φ) under certain conditions. To establish this, we first prove that the angle between ∂L2(θ,φ,ψ)∂θ and ∂L1(θ,φ) ∂θ are sufficiently positive if q(v, h) and P (v, h) satisfy certain conditions, as summarized in Lemma 1.
  Lemma 1. For any (θ, φ), there exists a symmetric positive definite matrix H such that ∂L2(θ,φ,ψ)∂θ = H ∂L1(θ,φ)∂θ under the assumption: || ∑ v,h δ(v, h) ∂E(v,h) ∂θ ||2 < || ∂L1(θ,φ) ∂θ ||2 if || ∂L1(θ,φ) ∂θ ||2 > 0 and
  || ∑ v,h δ(v, h) ∂E(v,h) ∂θ ||2 = 0 if || ∂L1(θ,φ) ∂θ ||2 = 0, where δ(v, h) = q(v, h)− P (v, h).
  Proof. According to the Cauchy-Schwarz inequality, we have
  〈∂L1(θ, φ) ∂θ , ∑ v,h δ(v, h) ∂E(v, h) ∂θ 〉 ≤ ||∂L1(θ, φ) ∂θ ||2|| ∑ v,h δ(v, h) ∂E(v, h) ∂θ ||2.
  If ||∂L1(θ,φ)∂θ ||2 > 0, according to the assumption || ∑ v,h δ(v, h) ∂E(v,h) ∂θ ||2 < || ∂L1(θ,φ) ∂θ ||2, we have
  〈∂L1(θ, φ) ∂θ , ∑ v,h δ(v, h) ∂E(v, h) ∂θ 〉 < ||∂L1(θ, φ) ∂θ ||22 = 〈 ∂L1(θ, φ) ∂θ , ∂L1(θ, φ) ∂θ 〉,
  which implies that
  〈∂L1(θ, φ) ∂θ , ∂L1(θ, φ) ∂θ − ∑ v,h δ(v, h) ∂E(v, h) ∂θ 〉 > 0.
  According to the definitions of L1(θ, φ) and L2(θ, φ, ψ), we have ∂L1(θ, φ) ∂θ − ∑ v,h δ(v, h) ∂E(v, h) ∂θ
  = ∑ v,h [PD(v)Q(h|v)− P (v, h)] ∂E(v, h) ∂θ − ∑ v,h [q(v, h)− P (v, h)]∂E(v, h) ∂θ
  = ∑ v,h [PD(v)Q(h|v)− q(v, h)] ∂E(v, h) ∂θ = ∂L2(θ, φ, ψ) ∂θ ,
  which implies that
  〈∂L1(θ, φ) ∂θ , ∂L2(θ, φ, ψ) ∂θ 〉 > 0.
  Equivalently, there exists a symmetric positive definite matrix H such that ∂L2(θ,φ,ψ)∂θ = H ∂L1(θ,φ) ∂θ . Note that this also holds when ||∂L1(θ,φ)∂θ ||2 = 0 (i.e., ∂L1(θ,φ) ∂θ = ~0) because ||∂L2(θ,φ,ψ)∂θ ||2 ≤
  ||∂L1(θ,φ)∂θ ||2+ || ∑ v,h δ(v, h) ∂E(v,h) ∂θ ||2 = 0 (i.e., ∂L2(θ,φ,ψ) ∂θ = ~0), according to the assumption.
  Remark Lemma 1 assumes that q(v, h) and P (v, h) are sufficiently close, which is encouraged by choosing a sufficiently powerful family of q(v, h) and updating ψ multiple times per update of θ, i.e. K1 > 1. If Lemma 1 holds, optimizing L2(θ, φ, ψ) with respect to θ will decrease L1(θ, φ) in expectation with a sufficiently small stepsize. Also note that for any (θ, ψ), ∂L2(θ,φ,ψ)∂φ = ∂L1(θ,φ) ∂φ and therefore, optimizing L2(θ, φ, ψ) with respect to φ will decrease L1(θ, φ) in expectation with a sufficiently small stepsize.
  Based on Lemma 1 and other commonly used assumptions in the analysis of stochastic gradient descent (Bottou et al., 2018), Algorithm 1 converges to a stationary point of L1(θ, φ), as stated in Theorem 1. Theorem 1. Solving the optimization problem in Eqn. (9) using stochastic gradient descent according to Algorithm 1, then
  lim k→∞ E[||∂L1(θk, φk) ∂θ ||22] = 0,
  under the following assumptions.
  1. The condition of Corollary 4.12 in (Bottou et al., 2018): L2(θ, φ, ψ) is twice differentiable with respect to θ, φ and ψ.
  2. Assumption 4.1 in (Bottou et al., 2018): the gradients of L2(θ, φ, ψ) with respect to θ, φ and ψ are Lipschitz.
  3. Assumption 4.3 in (Bottou et al., 2018): the first and second moments of the stochastic gradients are bounded by the expected gradients. 4. The stepsize satisfies the diminishing condition (Bottou et al., 2018), i.e., αk = γk,∑∞ k=1 γk =∞, ∑∞ k=1 γ 2 k <∞.
  5. The condition of Lemma 1 holds in each step k. Therefore,
  ∀k,∃Hk, ∂L2(θk, φk, ψk+1)
  ∂θ = Hk ∂L1(θk, φk) ∂θ .
  Proof. See Corollary 4.12 in (Bottou et al., 2018).
  Remark Assumption 1 and Assumption 2 in Theorem 1 are ensured because we use the sigmoid and tanh activation functions. Assumption 3 and Assumption 4 in Theorem 1 are ensured by the sampling and learning rate schemes of the stochastic gradient descent. Assumption 5 in Theorem 1 is weaker than the nonparametric assumption of Proposition 1 but still requires a large K1. Also note that the statement of converging to L1 in Theorem 1 is weaker than that in Proposition 1.
  C.3 COMPLEMENTARY CONVERGENCE THEOREM
  Heusel et al. (2017) propose a two-time scale update rule to train minimax optimization problems with a convergence guarantee even using K1 = 1. AdVIL converges if using the same training method as in (Heusel et al., 2017), which is summarized in Proposition 2. Proposition 2. AdVIL trained with a two-time scale update rule (Heusel et al., 2017) converges to a stationary local Nash equilibrium almost surely under the following assumptions.
  1. The gradients with respect to θ, φ and ψ are Lipschitz. 2. ∑ k αk =∞, ∑ k α 2 k <∞, ∑ k γk =∞, ∑ k γ 2 k <∞,γk = o(αk).
  3. The stochastic gradient errors are bounded in expectation.
  4. For each θ, the ordinary differentiable equation corresponding to Equation 13 has a local asymptotically stable attractor within a domain of attraction such that the attractor is Lipschitz. Similar assumptions are required for φ and ψ.
  5. supk ||θk|| <∞, supk ||ψk|| <∞, supk ||φk|| <∞.
  Proof. See Theorem 1 in (Heusel et al., 2017).
  Remark Compared to Theorem 1, Proposition 2 ensures the convergence of AdVIL without assuming q(v, h) is sufficiently close to P (v, h) in each step. However, a two time-sclae update rule (Heusel et al., 2017) is required to satisfy Assumption 2 and extra weight decay terms are needed to satisfy Assumption 4. Further, the convergence point is not necessarily a stationary point of L1(θ, φ).
sections.heading: D DATASETS AND EXPERIMENTAL SETTINGS
sections.text: We evaluate our method on the Digits dataset4, the UCI binary databases and the Frey faces datasets5. The information of the datasets is summarized in Tab. 3. We implement our model using the TensorFlow (Abadi et al., 2016) library. In all experiments, q and r are updated 100 times per update of P and Q, i.e. K1 = 100 and K2 = 1. We use the ADAM (Kingma & Ba, 2014) optimizer with the learning rate α = 0.0003, the moving average ratios β1 = 0.5 and β2 = 0.999, and the batch size of 500. We use a continuous z and the sigmoid activation function . All these hyperparameters are
  4https://scikit-learn.org/stable/modules/generated/ sklearn.datasets.load digits.html#sklearn.datasets.load digits 5http://www.cs.nyu.edu/∼roweis/data.html
  set according to the validation performance of an RBM on the Digits dataset and fixed throughout the paper unless otherwise stated. The sizes of the variational distributions depend on the structure of the MRF and are chosen according to the validation performance. The model structures in RBM and DBM experiments are summarized in Tab. 4 and Tab. 5, respectively.
  In a two-layer DBM with variables v, h1 and h2, we use an encoderQ(h1, h2|v) = Q(h1|v)Q(h2|h1) for both AdVIL and VCD. The decoder for AdVIL is the inverse of the encoder with one extra layer on the top, namely q(v, h1, h2) = ∫ q(v|h1)q(h1|h2)q(h2|z)q(z)dz. In our implementation, both AdVIL and VCD exploits that v and h2 are conditionally independent given h1. The layerwise structure potentially benefits the training of both methods. Nevertheless, in principle, any differentiable variational distributions be used in AdVIL and a systematical study is left for future work.
  The authors of NVIL (Kuleshov & Ermon, 2017) propose two variants. The first one employs a mixture of Bernoulli as q. The second one involves auxiliary variables and employs a neural network as q. Both variants scale up to an RBM of at most 64 visible units as reported in their paper (Kuleshov & Ermon, 2017). For a fair comparison, we carefully perform grid search over the default settings of NVIL and our settings based on their code and choose the best configuration. In this setting, the first variant of NVIL still fails to scale up to larger datasets and the best version of the second variant shares the same key hyperparameters as AdVIL, including K1 = 100 and a batch size of 500.
sections.heading: E MORE RESULTS
sections.text: E.1 EMPIRICAL VERIFICATION OF THEOREM 1 We empirically test the assumption || ∑ v,h δ(v, h) ∂E(v,h) ∂θk
  ||2 < ||∂L1(θk,φk)∂θk ||2 for k < 10000, where δ(v, h) = q(v, h)− P (v, h). Note that computing ||∂L1(θk,φk)∂θk ||2 exactly requires summing over v and h and therefore we train a small RBM ,where the dimensions of v, h and z are 4 on a synthetic dataset. The data distribution is a categorical distribution over {0, 1}4 and it is sampled from a Dirichlet distribution with all concentration parameters to be one. We get 10,000, 1,000 and 1,000 i.i.d samples from the categorical distribution for training, validation and test respectively. We find K1 = 10 is sufficient in this case.
  The results are shown in Fig. 5. It can be seen that a decoder with neural networks and auxiliary variables are sufficiently powerful to track the model distribution and therefore Assumption 5 in Theorem 1 holds in 10,000 steps. Besides, the model converges gradually, which agrees with our convergence analysis. The gap between the red curve in Fig. 5 (a) and the horizontal axis can be further reduced by using a more powerful q(v, h) and advanced optimization techniques for q(v, h).
  Theoretically, how Lemma 1 holds as the number of variables increases is not clear. Intuitively, we agree that it may get harder to satisfy this condition in a high-dimensional space. However, it is still possible with advanced optimization methods because the MRF is randomly initialized and learned gradually, and the variational decoder is trained to track the MRF (based on an old version of the decoder) after every update of the MRF. Empirically, computing both sides of the condition requires the value of the partition function, which is as hard as training the MRF. Therefore, verifying this condition during training in a high-dimensional case is highly nontrivial.
  E.2 SAMPLES IN RBM
  We present the samples from the RBM P and the decoder q in Fig. 6. In this case, we set the number of the hidden units to 50 and other settings remain the same as in Sec. 5.1. The first column demonstrates that the decoder is a good approximate sampler for the RBM. Note that the samples from the decoder are obtained from efficient ancestral sampling but those from the RBM is obtained by Gibbs sampling after 100,000 burn-in steps. The second column shows that ifH(Q) is removed, both models collapse to a certain mode of the data. The third column shows that ifH(q) is removed, both models fail to generate meaningful digits. These results demonstrate the importance of the entropy terms and the necessity of approximatingH(q) in a principled way.
  E.3 ADVIL WITH AN AUTOREGRESSIVE PRIOR
  Here we present the results of AdVIL with a neural autoregressive distribution estimator (NADE) (Larochelle & Murray, 2011) as the prior on the Digits dataset. We use the same RBM as in Sec. 5.1.
  The dimension of the latent units in NADE is 15, which is the same as the dimension of the auxiliary variables in the hierarchical decoder presented in Sec. 3.3.
  Compared to the hierarchical decoder, the NADE decoder has a tractable entropy and hence does not require r(z|h). However, getting samples from NADE is slow while AdVIL requires samples during training. Therefore, K1 = 5 for the NADE decoder has a similar training cost as the hierarchical decoder.
  Fig. 7 compares the two decoders. AdVIL with the NADE decoder achieves a slightly worse and unstable result.
  E.4 LEARNING CURVES AND ANALYSIS IN DBM
  We plot the learning curves of AdVIL and VCD in DBM, as shown in Fig. 8. AdVIL achieves a better result than VCD, which agrees with the quantitative results in Tab. 2. Note that we report the test NLL results in Tab. 2 according to the best validation performance.
  There are two types of biases introduced by using Q(h1, h2|v) = Q(h2|h1)Q(h1|v) in VCD. The first type of bias is introduced by using the approximate free energy in both the positive phase and the negative phase (See Eqn. (3) and Eqn. (4)). The second bias is introduced by the usage of Q(h1|v) to approximate P (h1|v) in the Gibbs sampling procedure to approximate the negative phase in Eqn. (4). The influence of the two types of biases on the negative phase is not clear, which can potentially explain the relatively inferior performance of VCD in DBM. In contrast, AdVIL approximates the negative phase by introducing one bias (i.e., the approximation error between q(v, h) and P (v, h)), whose effect on learning is theoretically characterized by Theorem 1. The above results and analysis essentially demonstrate the advantages of AdVIL over CD-based methods in DBM and the importance of developing black-box inference and learning algorithms for general MRFs.
  E.5 AIS RESULTS WITH STANDARD DEVIATIONS
  The AIS results in RBM and RBM with the means and standard deviations are shown in Tab. 6 and Tab. 7 respectively.