Year: 2018

Abstract:
In this paper we propose to view the acceptance rate of the Metropolis-Hastings algorithm as a universal objective for learning to sample from target distribution – given either as a set of samples or in the form of unnormalized density. This point of view unifies the goals of such approaches as Markov Chain Monte Carlo (MCMC), Generative Adversarial Networks (GANs), variational inference. To reveal the connection we derive the lower bound on the acceptance rate and treat it as the objective for learning explicit and implicit samplers. The form of the lower bound allows for doubly stochastic gradient optimization in case the target distribution factorizes (i.e. over data points). We empirically validate our approach on Bayesian inference for neural networks and generative models for images.

heading: 1 INTRODUCTION
text: Bayesian framework and deep learning have become more and more interrelated during recent years. Recently Bayesian deep neural networks were used for estimating uncertainty (Gal & Ghahramani, 2016), ensembling (Gal & Ghahramani, 2016) and model compression (Molchanov et al., 2017). On the other hand, deep neural networks may be used to improve approximate inference in Bayesian models (Kingma & Welling, 2014).
Learning modern Bayesian neural networks requires inference in the spaces with dimension up to several million by conditioning the weights of DNN on hundreds of thousands of objects. For such applications, one has to perform the approximate inference – predominantly by either sampling from the posterior with Markov Chain Monte Carlo (MCMC) methods or approximating the posterior with variational inference (VI) methods.
MCMC methods provide the unbiased (in the limit) estimate but require careful hyperparameter tuning especially for big datasets and high dimensional problems. The large dataset problem has been addressed for different MCMC algorithms: stochastic gradient Langevin dynamics (Welling & Teh, 2011), stochastic gradient Hamiltonian Monte Carlo (Chen et al., 2014), minibatch MetropolisHastings algorithms (Korattikara et al., 2014; Chen et al., 2016). One way to address the problem of high dimension is the design of a proposal distribution. For example, for the Metropolis-Hastings (MH) algorithm there exists a theoretical guideline for scaling the variance of a Gaussian proposal (Roberts et al., 1997; 2001). More complex proposal designs include adaptive updates of the proposal distribution during iterations of the MH algorithm (Holden et al., 2009; Giordani & Kohn, 2010). Another way to adapt the MH algorithm for high dimensions is combination of adaptive direction sampling and the multiple-try Metropolis algorithm as proposed in (Liu et al., 2000). Thorough overview of different extensions of the MH algorithm is presented in (Martino, 2018).
Variational inference is extremely scalable but provides a biased estimate of the target distribution. Using the doubly stochastic procedure (Titsias & Lázaro-Gredilla, 2014; Hoffman et al., 2013) VI can be applied to extremely large datasets and high dimensional spaces, such as a space of neural network weights (Kingma et al., 2015; Gal & Ghahramani, 2015; 2016). The bias introduced by variational approximation can be mitigated by using flexible approximations (Rezende & Mohamed, 2015) and resampling (Grover et al., 2018).
Generative Adversarial Networks (Goodfellow et al., 2014) (GANs) is a different approach to learn samplers. Under the framework of adversarial training different optimization problems could be solved efficiently (Arjovsky et al., 2017; Nowozin et al., 2016). The shared goal of ”learning to
sample” inspired the connection of GANs with VI (Mescheder et al., 2017) and MCMC (Song et al., 2017).
In this paper, we propose a novel perspective on learning to sample from a target distribution by optimizing parameters of either explicit or implicit probabilistic model. Our objective is inspired by the view on the acceptance rate of the Metropolis-Hastings algorithm as a quality measure of the sampler. We derive a lower bound on the acceptance rate and maximize it with respect to parameters of the sampler, treating the sampler as a proposal distribution in the Metropolis-Hastings scheme.
We consider two possible forms of the target distribution: unnormalized density (density-based setting) and a set of samples (sample-based setting). Each of these settings reveals a unifying property of the proposed perspective and the derived lower bound. In the density-based setting, the lower bound is the sum of forward and reverse KL-divergences between the true posterior and its approximation, connecting our approach to VI. In the sample-based setting, the lower bound admit a form of an adversarial game between the sampler and a discriminator, connecting our approach to GANs.
The closest work to ours is of Song et al. (2017). In contrast to their paper our approach (1) is free from hyperparameters; (2) is able to optimize the acceptance rate directly; (3) avoids minimax problem in the density based setting.
Our main contributions are as follows:
1. We introduce a novel perspective on learning to sample from the target distribution by treating the acceptance rate in the Metropolis-Hastings algorithm as a measure of sampler quality.
2. We derive the lower bound on the acceptance rate allowing for doubly stochastic optimization of the proposal distribution in case when the target distribution factorizes (i.e. over data points).
3. For sample-based and density-based forms of target distribution we show the connection of the proposed algorithm to variational inference and GANs.
The rest of the paper is organized as follows. In Section 2 we introduce the lower bound on the AR. Special forms of target distribution are addressed in Section 3. We validate our approach on the problems of approximate Bayesian inference in the space of high dimensional neural network weights and generative modeling in the space of images in Section 4. We discuss results and directions of the future work in Section 5.

heading: 2 ACCEPTANCE RATE FOR THE METROPOLIS-HASTINGS ALGORITHM

heading: 2.1 PRELIMINARIES
text: In MH algorithm we need to sample from target distribution p(x) while we are only able to sample from proposal distribution q(x′ |x). One step of the MH algorithm can be described as follows.
1. sample proposal point x′ ∼ q(x′ |x), given previously accepted point x
2. accept
{ x′, if p(x
′)q(x | x′) p(x)q(x′ | x) > u, u ∼ Uniform[0, 1]
x, otherwise
If the proposal distribution q(x′ |x) does not depend on x, i.e. q(x′ |x) = q(x′), the algorithm is called independent MH algorithm.
The quality of the proposal distribution is measured by acceptance rate and mixing time. Mixing time defines the speed of convergence of the Markov chain to the stationary distribution. The acceptance rate of the MH algorithm is defined as
AR = Eξ min{1, ξ} = ∫ dxdx′p(x)q(x′ |x) min { 1, p(x′)q(x |x′) p(x)q(x′ |x) } , (1)
where
ξ = p(x′)q(x |x′) p(x)q(x′ |x) , x ∼ p(x), x′ ∼ q(x′ |x). (2)
In case of independent proposal distribution we show that the acceptance rate defines a semimetric in distribution space between p and q (see Appendix A.2).

heading: 2.2 OPTIMIZING THE LOWER BOUND ON ACCEPTANCE RATE
text: Although, we can maximize the acceptance rate of the MH algorithm (Eq. 1) directly w.r.t. parameters φ of the proposal distribution qφ(x′ |x), we propose to maximize the lower bound on the acceptance rate. As our experiments show (see Section 4) the optimization of the lower bound compares favorably to the direct optimization of the acceptance rate. To introduce this lower bound we first express the acceptance rate in terms of total variation distance.
Theorem 1 For random variable ξ = p(x ′)q(x | x′)
p(x)q(x′ | x) , x ∼ p(x), x ′ ∼ q(x′ |x)
Eξ min{1, ξ} = 1− 1
2 Eξ|ξ − 1| = 1− TV
( p(x′)q(x |x′) ∥∥∥∥p(x)q(x′ |x)), (3) where TV is the total variation distance.
The proof of Theorem 1 can be found in Appendix A.1. This reinterpretation in terms of total variation allows us to lower bound the acceptance rate through the Pinsker’s inequality
AR ≥ 1−
√ 1
2 · KL
( p(x′)q(x |x′) ∥∥∥∥p(x)q(x′ |x)). (4) The maximization of this lower bound can be equivalently formulated as the following optimization problem
min φ
KL ( p(x′)qφ(x |x′) ∥∥∥∥p(x)qφ(x′ |x)). (5) In the following sections, we show the benefits of this optimization problem in two different settings — when the target distribution is given in a form of unnormalized density and as a set of samples. In Appendix C.5 and C.1 we provide the empirical evidence that maximization of the proposed lower bound results in the maximization of the acceptance rate.

heading: 3 OPTIMIZATION OF PROPOSAL DISTRIBUTION
text: From now on we consider only optimization problem Eq. 5 but the proposed algorithms can be also used for the direct optimization of the acceptance rate (Eq. 1).
To estimate the loss function (Eq. 5) we need to evaluate the density ratio. In the density-based setting unnormalized density of the target distribution is given, so we suggest to use explicit proposal distribution to compute the density ratio explicitly. In the sample-based setting, however, we cannot compute the density ratio, so we propose to approximate it via adversarial training (Goodfellow et al., 2014). The brief summary of constraints for both settings is shown in Table 1.
The following subsections describe the algorithms in detail.

heading: 3.1 DENSITY-BASED SETTING
text: In the density-based setting, we assume the proposal to be an explicit probabilistic model, i.e. the model that we can sample from and evaluate its density at any point up to the normalization constant. We also assume that the proposal is reparameterisable (Kingma & Welling, 2014; Rezende et al., 2014; Gal, 2016).
During the optimization of the acceptance rate we might face a situation when proposal collapses to the delta-function. This problem usually arises when we use Markov chain proposal, for example, qφ(x′ |x) = N (x′ |x, σ). For such proposal we can obtain arbitrary high acceptance rate, making the σ small enough. However, this is not the case for the independent proposal distribution qφ(x
′ |x) = qφ(x′). In Appendix B.1 we provide more details and intuition on this property of acceptance rate maximization. We also provide empirical evidence in Section 4 that collapsing to the delta-function does not happen for the independent proposal.
In this paper, we consider two types of explicit proposals: simple parametric family (Section 4.2) and normalizing flows (Rezende & Mohamed, 2015; Dinh et al., 2016) (Section 4.1). Rich family of normalizing flows allows to learn expressive proposal and evaluate its density in any point of target distribution space. Moreover, an invertible model (such as normalizing flow) is a natural choice for the independent proposal due to its ergodicity. Indeed, choosing the arbitrary point in the target distribution space, we can obtain the corresponding point in the latent space using the inverse function. Since every point in the latent space has positive density, then every point in the target space also has positive density.
Considering qφ(x′) as the proposal, the objective of optimization problem 5 takes the form L(p, qφ) = KL ( p(x′)qφ(x) ∥∥∥∥p(x)qφ(x′)) = E x ∼ p(x) x′ ∼ qφ(x′) log p(x)qφ(x ′) p(x′)qφ(x) . (6) Explicit form of the proposal qφ(x′) and the target p(x) distributions allows us to obtain density ratios qφ(x)/qφ(x′) and p(x′)/p(x) for any points x, x′. But to estimate the loss in Eq. 6 we also need to obtain samples from the target distribution x ∼ p(x) during training. For this purpose, we use the current proposal qφ and run the independent MH algorithm. After obtaining samples from the target distribution it is possible to perform optimization step by taking stochastic gradients w.r.t. φ. Pseudo-code for the obtained procedure is shown in Algorithm 1.
Algorithm 1 Optimization of proposal distribution in density-based case Require: explicit probabilistic model qφ(x′) Require: density of target distribution p̂(x) ∝ p(x)
while φ not converged do sample {x′k}Kk=1 ∼ qφ(x′) sample {xk}Kk=1 ∼ p(x) using independent MH with current proposal qφ L(p, qφ) ' 1K ∑K k=1 log p(xk)qφ(x ′ k)
p(x′k)qφ(xk) . approximate loss with finite number of samples
φ← φ− α∇φL(p, qφ) . perform gradient descent step end while return optimal parameters φ
Algorithm 1 could also be employed for the direct optimization of the acceptance rate (Eq. 1). Now we apply this algorithm for Bayesian inference problem and show that during optimization of the lower bound we can use minibatches of data, while it is not the case for direct optimization of the acceptance rate. We consider Bayesian inference problem for discriminative model on dataset D = {(xi, yi)}Ni=1, where xi is the feature vector of ith object and yi is its label. For the discriminative model we know likelihood p(yi |xi, θ) and prior distribution p(θ). In order to obtain predictions for some object xi, we need to evaluate the predictive distribution p(yi |xi) = Ep(θ | D)p(yi |xi, θ). (7) To obtain samples from posterior distribution p(θ | D) we suggest to learn proposal distribution qφ(θ) and perform independent MH algorithm. Thus the objective 6 can be rewritten as
L ( p(θ | D), qφ(θ) ) = KL ( p(θ′ | D)qφ(θ) ∥∥∥∥p(θ | D)qφ(θ′)) (8)
Note that due to the usage of independent proposal, the minimized KL-divergence splits up into the sum of two KL-divergences.
KL ( p(θ′ | D)qφ(θ) ∥∥∥∥p(θ | D)qφ(θ′)) = KL(qφ(θ)∥∥∥∥p(θ | D))+ KL(p(θ′ | D)∥∥∥∥qφ(θ′)) (9) Minimization of the first KL-divergence corresponds to the variational inference procedure.
KL ( qφ(θ) ∥∥∥∥p(θ | D)) = −Eθ∼qφ(θ) N∑ i=1 log p(yi |xi, θ) + KL(qφ(θ)‖p(θ)) + log p(D) (10)
The second KL-divergence has the only term that depends on φ. Thus we obtain the following optimization problem
min φ
[ − Eθ∼qφ(θ) N∑ i=1 log p(yi |xi, θ) + KL(qφ(θ)‖p(θ))− Eθ∼p(θ | D) log qφ(θ) ] . (11)
The first summand here contains the sum over all objects in dataset D. We follow doubly stochastic variational inference and suggest to perform unbiased estimation of the gradient in problem 11 using only minibatches of data. Moreover, we can use recently proposed techniques (Korattikara et al., 2014; Chen et al., 2016) that perform the independent MH algorithm using only minibatches of data. Combination of these two techniques allows us to use only minibatches of data during iterations of algorithm 1. In the case of the direct optimization of the acceptance rate, straightforward usage of minibatches results in biased gradients. Indeed, for the direct optimization of the acceptance rate (Eq. 1) we have the product over the all training data inside min function.

heading: 3.2 SAMPLE-BASED SETTING
text: In the sample-based setting, we assume the proposal to be an implicit probabilistic model, i.e. the model that we can only sample from. As in the density-based setting, we assume that we are able to perform the reparameterization trick for the proposal.
In this subsection we consider only Markov chain proposal qφ(x′ |x), but everything can be applied to independent proposal qφ(x′) by simple substitution qφ(x′ |x) with qφ(x′). From now we will assume our proposal distribution to be a neural network that takes x as its input and outputs x′. Considering proposal distribution parameterized by a neural network allows us to easily exclude delta-function from the space of solutions. We avoid learning the identity mapping by using neural networks with the bottleneck and noisy layers. For the detailed description of the architectures see Appendix C.8.
The set of samples from the true distribution X ∼ p(x) allows for the Monte Carlo estimation of the loss
L(p, qφ) = E x ∼ p(x) x′ ∼ qφ(x′ |x) log p(x)qφ(x ′ |x) p(x′)qφ(x |x′) . (12)
To compute the density ratio p(x)qφ(x ′ | x)
p(x′)qφ(x | x′) we suggest to use well-known technique of density ratio estimation via training discriminator network. Denoting discriminator output as D(x, x′), we suggest the following optimization problem for the discriminator.
min D
[ −E x ∼ p(x)
x′ ∼ qφ(x′ |x) logD(x, x′)−E x ∼ p(x) x′ ∼ qφ(x′ |x) log(1−D(x′, x))
] (13)
Speaking informally, such discriminator takes two images as input and tries to figure out which image is sampled from true distribution and which one is generated by the one step of proposal distribution. It is easy to show that optimal discriminator in problem 13 will be
D(x, x′) = p(x)qφ(x ′ |x) p(x)qφ(x′ |x) + p(x′)qφ(x |x′) . (14)
Note that for optimal discriminator we have D(x, x′) = 1 − D(x′, x). In practice, we have no optimal discriminator and these values can differ significantly. Thus, we have four ways for density ratio estimation that may differ significantly.
p(x)qφ(x ′ |x)
p(x′)qφ(x |x′) ≈ D(x, x
′) 1−D(x, x′) ≈ 1−D(x ′, x) D(x′, x) ≈ 1−D(x ′, x) 1−D(x, x′) ≈ D(x, x ′) D(x′, x) (15)
To avoid the ambiguity we suggest to use the discriminator of a special structure. Let D̃(x, x′) be a convolutional neural network with scalar output. Then the output of discriminator D(x, x′) is defined as follows.
D(x, x′) = exp(D̃(x, x′))
exp(D̃(x, x′)) + exp(D̃(x′, x)) (16)
In other words, such discriminator can be described as the following procedure. For single neural network D̃(·, ·) we evaluate two outputs D̃(x, x′) and D̃(x′, x). Then we take softmax operation for these values. Summing up all the steps, we obtain algorithm 2.
Algorithm 2 Optimization of proposal distribution in sample-based case Require: implicit probabilistic model qφ(x′ |x) Require: large set of samples X ∼ p(x)
for n iterations do sample {xk}Kk=1 ∼ X sample {x′k}Kk=1 ∼ qφ(x′|x) train discriminator D by optimizing 13 L(p, qφ) ≈ 1K ∑K k=1 log D(xk,x ′ k)
1−D(xk,x′k) . approximate loss with finite number of samples
φ← φ− α∇φL(p, qφ) . perform gradient descent step end for return parameters φ
Algorithm 2 could also be employed for direct optimization of the acceptance rate (Eq. 1). But, in Appendix B.2 we provide an intuition for this setting that the direct optimization of the acceptance rate may struggle from vanishing gradients.

heading: 4 EXPERIMENTS
text: In this section, we provide experiments for both density-based and sample-based settings, showing the proposed procedure is applicable to high dimensional target distributions. Code for reproducing all of the experiments will be published with the camera-ready version of the paper.

heading: 4.1 SYNTHETIC DISTRIBUTIONS
text: To demonstrate performance of our approach we reproduce the experiment from (Song et al., 2017). For target distributions we use synthetic 2d distributions (see Appendix C.3 for densities): ring (a ring-shaped density), mog2 (a mixture of 2 Gaussians), mog6 (a mixture of 6 Gaussians), ring5 (a mixture of 5 distinct rings). We measure performance of learned samplers using Effective Sample Size (see Appendix C.4 for formulation). Since the unnormalized densities of target distributions are given, we can learn proposals as suggested in the density-based setting (Section 3.1).
To learn the independent proposal we use RealNVP model (Dinh et al., 2016) (see details in Appendix C.2) and compare the performance of proposals after optimization of different objectives: the acceptance rate (AR), our lower bound on the acceptance rate (ARLB), evidence lower bound that corresponds to the variational inference (VI). We also compare the performance of obtained independent proposals with the performance of Markov chain proposals: A-NICE-MC (Song et al., 2017), Hamiltonian Monte Carlo (HMC).
In Tables 2, 3 we see that our approach has comparable performance with A-NICE-MC (Song et al., 2017). However, comparison between A-NICE-MC and learning independent proposal is not the main subject of interest, since A-NICE-MC learns Markov chain proposal. On the one hand, Markov
chain proposal uses more information while generating a new sample, hence can learn more expressive stationary distribution, on the other hand, usage of previous sample increase autocorrelation between samples and reduces ESS. Thus, the main point of interest is the comparison of two independent proposals: one is learned by maximization of the acceptance rate (or its lower bound), and the second is learned by variational inference procedure, i.e. maximization of evidence lower bound. In Table 2 we see that both maximization of the acceptance rate and its lower bound outperform variational inference for all target distributions. Moreover, in Fig. 1 we show that variational inference fails to cover all the modes of mog6 in contrast to proposals learned via maximization of acceptance rate or its lower bound. Densities of learned proposals and histograms for all distributions are presented in Appendix C.6.
Table 3: Performance of learned proposals as measured by Effective Sample Size per second. All the numbers are rounded to integers. See description of the compared methods in the text.
Target A-NICE-MC AR (ours) ARLB (ours) VI ring 706417 275164 275488 227253 mog2 231253 267491 205552 101074 mog6 143877 94938 112034 3663 ring5 85186 95314 70635 48224
0 20 40 60 80
0
20
40
60
80
Variational Inference
0 20 40 60 80
Acceptance Rate
0 20 40 60 80
Acceptance Rate Lower Bound
0
150
300 450
Figure 1: 2d histrogram of samples from the MH algorithm with different proposals. From left to right proposals are learned by: variational inference, the acceptance rate maximization, the acceptance rate lower bound maximization.

heading: 4.2 DENSITY-BASED SETTING
text: In density-based setting, we consider Bayesian inference problem for the weights of a neural network. In our experiments we consider approximation of predictive distribution (Eq. 7) as our main goal. To estimate the goodness of the approximation we measure negative log-likelihood and accuracy on the test set.
In subsection 3.1 we show that lower bound on acceptance rate can be optimized more efficiently than acceptance rate due to the usage of minibatches. But other questions arise.
1. Does the proposed objective in Eq. 11 allow for better estimation of predictive distribution compared to the variational inference?
2. Does the application of the MH correction to the learned proposal distribution allow for better estimation of the predictive distribution (Eq. 7) than estimation via raw samples from the proposal?
To answer these questions we consider reduced LeNet-5 architecture (see Appendix C.7) for classification task on 20k images from MNIST dataset (for test data we use all of the MNIST test
set). Even after architecture reduction we still face a challenging task of learning a complex distribution in 8550-dimensional space. For the proposal distribution we use fully-factorized gaussian qφ(θ) = ∏d j=1N (θj |µj , σj) and standard normal distribution for prior p(θ) = ∏d j=1N (θj | 0, 1).
For variational inference, we train the model using different initialization and pick the model according to the best ELBO. For our procedure, we do the same and choose the model by the maximum value of the acceptance rate lower bound. In Algorithm 1 we propose to sample from the posterior distribution using the independent MH and the current proposal. It turns out in practice that it is better to use the currently learned proposal qφ(θ) = N (θ |µ,σ) as the initial state for random-walk MH algorithm. That is, we start with the mean µ as an initial point, and then use random-walk proposal q(θ′ | θ) = N (θ′ | θ,σ) with the variances σ of current independent proposal. This should be considered as a heuristic that improves the approximation of the loss function.
The optimization of the acceptance rate lower bound results in the better estimation of predictive distribution than the variational inference (see Fig. 2). Optimization of acceptance rate for the same number of epochs results in nearly 30% accuracy on the test set. That is why we do not report results for this procedure in Fig. 2.
To answer the second question we estimate predictive distribution in two ways. The first way is to perform 100 accept/reject steps of the independent MH algorithm with the learned proposal qφ(θ)
after each epoch, i.e. perform MH correction of the samples from the proposal. The second way is to take the same number of samples from qφ(θ) without MH correction. For both estimations of predictive distribution, we evaluate negative log-likelihood on the test set and compare them.
The MH correction of the learned proposal improves the estimation of predictive distribution for the variational inference (right plot of Fig. 3) but does not do so for the optimization of the acceptance rate lower bound (left plot of Fig. 3). This fact may be considered as an implicit evidence that our procedure learns the proposal distribution with higher acceptance rate.

heading: 4.3 SAMPLE-BASED SETTING
text: In the sample-based setting, we estimate density ratio using a discriminator. Hence we do not use the minibatching property (see subsection 3.1) of the obtained lower bound, and optimization problems for the acceptance rate and for the lower bound have the same efficiency in terms of using data. That is why our main goal in this setting is to compare the optimization of the acceptance rate and the optimization of the lower bound. Also, in this setting, we have Markov chain proposal that is interesting to compare with the independent proposal. Summing up, we formulate the following questions:
1. Does the optimization of the lower bound has any benefits compared to the direct optimization of the acceptance rate?
2. Do we have mixing issue while learning Markov chain proposal in practice? 3. Could we improve the visual quality of samples by applying the MH correction to the
learned proposal?
We use DCGAN architecture for the proposal and discriminator (see Appendix C.8) and apply our algorithm to MNIST dataset. We consider two optimization problems: direct optimization of the acceptance rate and its lower bound. We also consider two ways to obtain samples from the approximation of the target distribution — use raw samples from the learned proposal, or perform the MH algorithm, where we use the learned discriminator for density ratio estimation.
In case of the independent proposal, we show that the MH correction at evaluation step allows to improve visual quality of samples — figures 4(a) and 4(b) for the direct optimization of acceptance rate, figures 4(c) and 4(d) for the optimization of its lower bound. Note that in Algorithm 2 we do not apply the independent MH algorithm during training. Potentially, one can use the MH algorithm considering any generative model as a proposal distribution and learning a discriminator for density ratio estimation. Also, for this proposal, we demonstrate the negligible difference in visual quality of samples obtained by the direct optimization of acceptance rate (see Fig. 4(a)) and by the optimization of the lower bound (see Fig. 4(c)).
In the case of the Markov chain proposal, we show that the direct optimization of acceptance rate results in slow mixing (see Fig. 5(a)) — most of the time the proposal generates samples from one of the modes (digits) and rarely switches to another mode. When we perform the optimization of the lower bound the proposal switches between modes frequently (see Fig. 5(b)).
0 1 2 3 4 5 6 7 8 9
0
10
20
30
40
50
60
70
80
90
0 1 2 3 4 5 6 7 8 9
0
10
20
30
40
50
60
70
80
90
Figure 5: Samples from the chain obtained via the MH algorithm with the learned proposal and the learned discriminator for density ratio estimation. Fig. 5(a) corresponds to the direct optimization of the acceptance rate. Fig. 5(b) – to optimization of the lower bound on acceptance rate. Samples in the chain are obtained one by one from left to right from top to bottom.
To show that the learned proposal distribution has the Markov property rather than being totally independent, we show samples from the proposal conditioned on two different points in the dataset (see Fig. 6). The difference in samples from two these distributions (Fig. 6(a), 6(a)) reflects the dependence on the conditioning.
Additionally, in Appendix C.9 we present samples from the chain after 10000 accepted images and also samples from the chain that was initialized with noise.

heading: 5 DISCUSSION AND FUTURE WORK
text: This paper proposes to use the acceptance rate of the MH algorithm as the universal objective for learning to sample from some target distribution. We also propose the lower bound on the acceptance rate that should be preferred over the direct maximization of the acceptance rate in many cases. The proposed approach provides many ways of improvement by the combination with techniques from the recent developments in the field of MCMC, GANs, variational inference. For example
• The proposed loss function can be combined with the loss function from (Levy et al., 2017), thus allowing to learn the Markov chain proposal in the density-based setting.
• We can use stochastic Hamiltonian Monte Carlo (Chen et al., 2014) for the loss estimation in Algorithm 1.
• In sample-based setting one can use more advanced techniques of density ratio estimation.
Application of the MH algorithm to improve the quality of generative models also requires exhaustive further exploration and rigorous treatment.

heading: A ACCEPTANCE RATE OF THE MH ALGORITHM

heading: A.1 PROOF OF THEOREM 1
text: Remind that we have random variables ξ = p(x ′)q(x | x′)
p(x)q(x′ | x) , x ∼ p(x), x ′ ∼ q(x′ |x) and u ∼
Uniform[0, 1], and want to prove the following equalities.
Eξ min{1, ξ} = P{ξ > u} = 1− 1
2 Eξ|ξ − 1| (17)
Equality Eξ min{1, ξ} = P{ξ > u} is obvious. Eξ min{1, ξ} = ∫ ∞ 0 pξ(x) min{1, x}dx = ∫ x≥1 pξ(x)dx+ ∫ x<1 pξ(x)xdx (18)
P{ξ > u} = ∫ ∞ 0 dxpξ(x) ∫ x 0 [0 ≤ u ≤ 1]du = ∫ x≥1 pξ(x)dx+ ∫ x<1 pξ(x)xdx (19)
Equality P{ξ > u} = 1− 12Eξ|ξ − 1| can be proofed as follows.
P{ξ > u} = ∫ 1 0 du ∫ +∞ u pξ(x)dx = ∫ 1 0 (1− Fξ(u))du = (20)
= 1− [ uFξ(u) ∣∣∣∣1 0 − ∫ 1 0 upξ(u)du ] = 1− Fξ(1) + ∫ 1 0 upξ(u)du, (21)
where Fξ(u) is CDF of random variable ξ. Note that Fξ(0) = 0 since ξ ∈ (0,+∞]. Eq. 21 can be rewritten in two ways.
1− Fξ(1) + ∫ 1 0 upξ(u)du = 1 + ∫ 1 0 (u− 1)pξ(u)du = 1− ∫ 1 0 |u− 1|pξ(u)du (22)
To rewrite Eq. 21 in the second way we note that Eξ = 1. 1−Fξ(1)+ ∫ 1 0 upξ(u)du = ∫ +∞ 1 pξ(u)du+1− ∫ +∞ 1 upξ(u)du = 1− ∫ +∞ 1 |u−1|pξ(u)du (23)
Summing equations 22 and 23 results in the following formula
P{ξ > u} = 1− 1 2 Eξ|ξ − 1|. (24)
Using the form of ξ we can rewrite the acceptance rate as
1− 1 2 Eξ|ξ − 1| = 1− TV
( p(x′)q(x |x′) ∥∥∥∥p(x)q(x′ |x)). (25)

heading: A.2 ACCEPTANCE RATE OF INDEPENDENT MH DEFINES SEMIMETRIC IN DISTRIBUTION SPACE
text: In independent case we have ξ = p(x ′)q(x)
p(x)q(x′) , x ∼ p(x), x ′ ∼ q(x′) and we want to prove that
Eξ|ξ − 1| is semimetric (or pseudo-metric) in space of distributions. For this appendix, we denote D(p, q) = Eξ|ξ − 1|. The first two axioms for metric obviously holds
1. D(p, q) = 0 ⇐⇒ p = q 2. D(p, q) = D(q, p)
There is an example when triangle inequality does not hold. For distributions p = Uniform[0, 2/3], q = Uniform[1/3, 1], s = Uniform[0, 1]
D(p, s) +D(q, s) = 4
3 <
3 2 = D(p, q). (26)
But weaker inequality can be proved.
D(p, s) +D(q, s) = ∫ |p(x)s(y)− p(y)s(x)|dydx+ ∫ |q(x)s(y)− q(y)s(x)|dydx = (27)
= ∫ [ | p(x)s(y)q(z)︸ ︷︷ ︸
a
− p(y)s(x)q(z)︸ ︷︷ ︸ b |+ | q(x)s(y)p(z)︸ ︷︷ ︸ c − q(y)s(x)p(z)︸ ︷︷ ︸ d
| ] dxdydz (28)
D(p, s) +D(q, s) = ∫ |p(z)s(y)q(x)− p(y)s(z)q(x)|dxdydz+ (29)
+ ∫ |q(x)s(z)p(y)− q(z)s(x)p(y)|dxdydz ≥ ∫ ∣∣∣∣ q(x)s(y)p(z)︸ ︷︷ ︸ c − p(y)s(x)q(z)︸ ︷︷ ︸ b ∣∣∣∣dxdydz (30)
D(p, s) +D(q, s) = ∫ |p(z)s(x)q(y)− p(x)s(z)q(y)|dxdydz+ (31)
+ ∫ |q(y)s(z)p(x)− q(z)s(y)p(x)|dxdydz ≥ ∫ ∣∣∣∣ q(y)s(x)p(z)︸ ︷︷ ︸ d − p(x)s(y)q(z)︸ ︷︷ ︸ a ∣∣∣∣dxdydz (32) Summing up equations 28, 30 and 32 we obtain
3(D(p, s) +D(q, s)) ≥ ∫ dxdydz [ |a− b|+ |c− d|+ |c− b|+ |d− a| ] ≥ 2 ∫ dxdydz|d− b| =
(33)
= 2 ∫ dxdydzs(x) ∣∣∣∣q(y)p(z)− q(z)p(y)∣∣∣∣ = 2D(p, q) (34) D(p, s) +D(q, s) ≥ 2
3 D(p, q) (35)

heading: B OPTIMIZATION OF PROPOSAL DISTRIBUTION

heading: B.1 ON COLLAPSING TO THE DELTA-FUNCTION
text: Firstly, let’s consider the case of gaussian random-walk proposal q(x′ |x) = N (x′ |x, σ). The optimization problem for the acceptance rate takes the form
AR = ∫ dxdx′p(x)N (x′ |x, σ) min { 1, p(x′)
p(x)
} → max
σ . (36)
It is easy to see that we can obtain acceptance rate arbitrarly close to 1, taking σ small enough.
In the case of the independent proposal, we don’t have the collapsing to the delta-function problem. In our work, it is important to show non-collapsing during optimization of the lower bound, but the same hold for the direct optimization of the acceptance rate. To provide such intuition we consider one-dimensional case where we have some target distribution p(x) and independent proposal q(x) = N (x |µ, σ). Choosing σ small enough, we approximate sampling with the independent MH as sampling on some finite support x ∈ [µ − a, µ + a]. For this support, we approximate the target distribution with the uniform distribution (see Fig. 7).
For such approximation, optimization of lower bound takes the form
min q
[ KL(p(x)‖q(x)) + KL(q(x)‖p(x)) ] (37)
min σ
[ KL(Uniform[−a, a]‖N (x | 0, σ,−a, a)) + KL(N (x | 0, σ,−a, a)‖Uniform[−a, a]) ] (38)
Here N (x | 0, σ,−a, a) is truncated normal distribution. The first KL-divergence can be written as follows.
KL(Uniform[−a, a]‖N (x | 0, σ,−a, a)) = − 1 2a ∫ a −a dx logN (x | 0, σ,−a, a)− log 2a = (39)
= − 1 2a
[ − 2a log(σZ)− a log 2π − 1
2σ2 2a3 3
] − log 2a = (40)
= log σ + logZ + a2
6σ2 +
1 2 log 2π − log 2a (41)
Here Z is normalization constant of truncated log normal distribution and Z = Φ(a/σ)−Φ(−a/σ), where Φ(x) is CDF of standard normal distribution. The second KL-divergence is
KL(N (x | 0, σ,−a, a)‖Uniform[−a, a]) = (42)
= −1 2 log(2πe)− log σ − logZ + a√ 2πσZ exp
( − a 2
2σ2
) + log 2a (43)
Summing up two KL-divergencies and taking derivative w.r.t. σ we obtain ∂
∂σ
( KL(Uniform[−a, a]‖N (x | 0, σ,−a, a)) + KL(N (x | 0, σ,−a, a)‖Uniform[−a, a]) ) =
(44)
= − a 2
3σ3 + a3√ 2πσ4Z exp
( − a 2
2σ2
) +
a√ 2π exp
( − a 2
2σ2 )[ − 1 σ2Z − 1 σZ2 −2a σ2 √ 2π exp ( − a 2 2σ2 )] =
(45)
= 1
a
[ − a 3
3σ3 + a2√ 2πσ2Z exp
( − a 2
2σ2
)( a2
σ2 − 1 + 2a√ 2πσZ exp
( − a 2
2σ2
))] (46)
To show that the derivative of the lower bound w.r.t. σ is negative, we need to prove that the following inequality holds for positive x.
−1 3 x3+ x2√ 2π(Φ(x)− Φ(−x))
exp(−x2/2) ( x2−1+ 2x√
2π(Φ(x)− Φ(−x)) exp(−x2/2)
) < 0, x > 0
(47) Defining φ(x) = ∫ x 0 e−t 2/2dt and noting that 2φ(x) = √
2π(Φ(x) − Φ(−x)) we can rewrite inequality 47 as
1
φ(x) e−x 2/2
( x2 − 1 + 2xe −x2/2
φ(x)
) < 2x
3 , x > 0 (48)
By the fundamental theorem of calculus, we have
xe−x 2/2 = ∫ x 0 e−t 2/2(1− t2)dt (49)
Hence,
φ(x)− xe−x 2/2 = ∫ x 0 e−t 2/2t2dt ≥ e−x 2/2 ∫ x 0 t2dt = e−x 2/2x 3 3 (50)
Or equivalently,
φ(x) ≥ e−x 2/2x
3 + 3x
3 (51)
Using this inequality twice, we obtain
e−x 2/2
φ(x) ≤ 3 x(x2 + 3)
(52)
and
x2 − 1 + xe −x2/2
φ(x) ≤ x2 − 1 + 3 x2 + 3 = x2(2 + x2) x2 + 3 (53)
Thus, the target inequality can be verified by the verification of
3x(2 + x2)
(x2 + 3)2 ≤ 2x 3 . (54)
Thus, we show that partial derivative of our lower bound w.r.t. σ is negative. Using that knowledge we can improve our loss by taking a bigger value of σ. Hence, such proposal does not collapse to delta-function.
B.2 INTUITION FOR BETTER GRADIENTS IN SAMPLE-BASED SETTING
In this section, we provide an intuition for sample-based setting that the loss function for lower bound has better gradients than the loss function for acceptance rate. Firstly, we remind that in the sample-based setting we use a discriminator for density ratio estimation.
D(x, x′) = p(x)q(x′ |x)
p(x)q(x′ |x) + p(x′)q(x |x′) (55)
For this purpose we use the discriminator of special structure
D(x, x′) = exp(D̃(x, x′))
exp(D̃(x, x′)) + exp(D̃(x′, x)) =
1
1 + exp ( − (D̃(x, x′)− D̃(x′, x)) ) (56) We denote d(x, x′) = D̃(x, x′)− D̃(x′, x) and consider the case when the discriminator can easily distinguish fake pairs from valid pairs. So D(x, x′) is close to 1 and d(x, x′) 0 for x ∼ p(x) and x′ ∼ q(x′ |x). To evaluate gradients we consider Monte Carlo estimations of each loss and take gradients w.r.t. x′ in order to obtain gradients for parameters of proposal distribution. We do not introduce the reparameterization trick to simplify the notation but assume it to be performed. For the optimization of the acceptance rate we have∫
dxdx′p(x)q(x′ |x) ∣∣∣∣p(x′)q(x |x′)p(x)q(x′ |x) − 1 ∣∣∣∣ ' ∣∣∣∣p(x′)q(x |x′)p(x)q(x′ |x) − 1 ∣∣∣∣ (57)
LAR = ∣∣∣∣p(x′)q(x |x′)p(x)q(x′ |x) − 1 ∣∣∣∣ ≈ ∣∣∣∣1−D(x, x′)D(x, x′) − 1 ∣∣∣∣ (58) ∂LAR ∂x′ = 1 D2(x, x′) ∂D(x, x′) ∂x′ = exp(−d(x, x′))∂d(x, x ′) ∂x′ (59)
While for the optimization of the lower bound we have∫ dxdx′p(x)q(x′ |x) log ( p(x)q(x′ |x) p(x′)q(x |x′) ) ' log ( p(x)q(x′ |x) p(x′)q(x |x′) ) (60)
LLB = − log ( p(x′)q(x |x′) p(x)q(x′ |x) ) ≈ − log ( 1−D(x, x′) D(x, x′) ) (61)
∂LLB ∂x′ = 1 (1−D(x, x′))D(x, x′) ∂D(x, x′) ∂x′ = ∂d(x, x′) ∂x′ (62)
Now we compare Eq. 59 and Eq. 62. We see that in case of strong discriminator we have vanishing gradients in Eq. 59 due to exp(−d(x, x′)), while it is not the case for Eq. 62.

heading: C EXPERIMENTS

heading: C.1 TOY PROBLEM
text: This experiment shows that it is possible to optimize the acceptance rate, optimizing its lower bound. For the target distribution we consider bimodal Gaussian p(x) = 0.5 · N (x | − 2, 0.5) + 0.5 · N (x | 2, 0.7), for the independent proposal we consider unimodal gaussian q(x) = N (x |µ, σ). We perform stochastic gradient optimization using Algorithm 1 from the same initialization for both objectives (Fig. 8) and obtain approximately the same local maximums.

heading: C.2 ARCHITECTURE OF THE REALNVP PROPOSAL
text: For the proposal distribution we use similar architecture to the NICE proposal. The RealNVP model (Dinh et al., 2016) use the same strategy for evaluating the Jacobian as the NICE model does. Each coupling layer define the following function. Given a D dimensional input x and d < D, the output y is evaluated by the formula
y1:d = x1:d,
yd+1:D = xd+1:D exp(s(x1:d)) + t(x1:d),
where the functions s, t can be arbitrary complex, since the structure of the functions doesn’t influence the computation of the Jacobian.
For our proposal we use 4 coupling layers with s and t consist of two fully-connected layers with hidden dimension of 256.

heading: C.3 SYNTHETIC DISTRIBUTIONS DENSITIES
text: For synthetic distributions we consider the same distributions as in Song et al. (2017).
The analytic form of p(x) for ring is: p(x) ∝ exp(−U(x)), U(x) = ( √ x21 + x 2 2 − 2)2
0.32 (63)
The analytic form of p(x) for mog2 is:
p(x) = 1
2 N (x|µ1, σ1) +
1 2 N (x|µ2, σ2) (64)
where µ1 = [5, 0], µ2 = [−5, 0], σ1 = σ2 = [0.5, 0.5]. The analytic form of p(x) for mog6 is:
p(x) = 1
6 6∑ i=1 N (x|µi, σi) (65)
where µi = [sin iπ3 , cos iπ 3 ] and σi = [0.5, 0.5].
The analytic form of p(x) for ring5 is:
p(x) ∝ exp(−U(x)), U(x) = min(u1, u2, u3, u4, u5) (66) where ui = ( √ x21 + x 2 2 − i)2/0.04.

heading: C.4 EFFECTIVE SAMPLE SIZE FORMULATION
text: For the effective sample size formulation we follow Song et al. (2017).
Assume a target distribution p(x), and a Markov chain Monte Carlo (MCMC) sampler that produces a set of N correlated samples {xi}N1 from some distribution q({xi}N1 ) such that q(xi) = p(xi). Suppose we are estimating the mean of p(x) through sampling; we assume that increasing the number of samples will reduce the variance of that estimate.
Let V = Varq[ ∑N i=1 xi/N ] be the variance of the mean estimate through the MCMC samples. The effective sample size (ESS) of {xi}N1 , which we denote as M = ESS({xi}N1 ), is the number of independent samples from p(x) needed in order to achieve the same variance, i.e. Varp[ ∑M j=1 xj/M ] = V . A practical algorithm to compute the ESS given {xi}N1 is provided by:
ESS({xi}N1 ) = N 1 + 2 ∑N−1 s=1 (1− s N )ρs
(67)
where ρs denotes the autocorrelation under q of x at lag s. We compute the following empirical estimate ρ̂s for ρs:
ρ̂s = 1
σ̂2(N − s) N∑ n=s+1 (xn − µ̂)(xn−s − µ̂) (68)
where µ̂ and σ̂ are the empirical mean and variance obtained by an independent sampler.
Due to the noise in large lags s, we adopt the approach of Hoffman & Gelman (2014) where we truncate the sum over the autocorrelations when the autocorrelation goes below 0.05.

heading: C.5 OPTIMIZATION OF THE LOWER BOUND
text: In this section we provide the empirical evidence that maximization of the proposed lower bound on the acceptance rate (ARLB) results in maximization of the acceptance rate (AR). For that purpose we evaluate ARLB and AR at each iteration during the optimization of ARLB. After training we evaluate correlation coefficient between ARLB and logarithm of AR. The curves are shown in Fig. 9. Correlation coefficients for different distributions are: −0.914 (ring), −0.905 (mog2), −0.956 (mog6), −0.982 (ring5).

heading: C.6 LEARNED PROPOSALS
text: In this section we provide levelplots of learned proposals densities (see Fig. 10). We also provide 2d histrograms of samples from the MH algorithm using the corresponding proposals (see Fig. 11).

heading: C.7 ARCHITECTURE OF THE REDUCED LENET-5
text: class LeNet5(BayesNet): def __init__(self):
super(LeNet5, self).__init__() self.num_classes = 10 self.conv1 = layers.ConvFFG(1, 10, 5, padding=0) self.relu1 = nn.ReLU(True) self.pool1 = nn.MaxPool2d(2, padding=0) self.conv2 = layers.ConvFFG(10, 20, 5, padding=0) self.relu2 = nn.ReLU(True) self.pool2 = nn.MaxPool2d(2, padding=0) self.flatten = layers.ViewLayer([20*4*4]) self.dense1 = layers.LinearFFG(20*4*4, 10) self.relu3 = nn.ReLU() self.dense2 = layers.LinearFFG(10, 10)

heading: C.8 ARCHITECTURES OF NEURAL NETWORKS IN SAMPLE-BASED SETTING
text: In sample-based setting we use usual DCGAN architecture for independent proposal distribution
class Generator(layers.ModuleWrapper): def __init__(self):
super(Generator, self).__init__() self.fc = nn.Linear(100, 128*8*8) self.unflatten = layers.ViewLayer([128, 8, 8]) self.in1 = nn.InstanceNorm2d(128) self.us1 = nn.ConvTranspose2d(128, 128, 2, 2) self.conv1 = nn.Conv2d(128, 128, 3, stride=1, padding=1) self.in2 = nn.InstanceNorm2d(128, 0.8) self.lrelu1 = nn.LeakyReLU(0.2, inplace=True) self.us2 = nn.ConvTranspose2d(128, 128, 2, 2) self.conv2 = nn.Conv2d(128, 64, 3, stride=1, padding=1) self.in3 = nn.InstanceNorm2d(64, 0.8) self.lrelu2 = nn.LeakyReLU(0.2, inplace=True) self.conv3 = nn.Conv2d(64, 1, 3, stride=1, padding=1) self.tanh = nn.Tanh()
And a little be modified acrhitecture for Markov chain proposal distribution
class Generator(layers.ModuleWrapper): def __init__(self):
super(Generator, self).__init__()
self.d_conv1 = nn.Conv2d(1, 16, 5, stride=2, padding=2) self.d_lrelu1 = nn.LeakyReLU(0.2, inplace=True) self.d_do1 = nn.Dropout2d(0.5) self.d_conv2 = nn.Conv2d(16, 4, 5, stride=2, padding=2) self.d_in2 = nn.InstanceNorm2d(4, 0.8) self.d_lrelu2 = nn.LeakyReLU(0.2, inplace=True) self.d_do2 = nn.Dropout2d(0.5)
self.b_view = layers.ViewLayer([4*8*8]) self.b_fc = nn.Linear(4*8*8, 256) self.b_lrelu = nn.LeakyReLU(0.2, inplace=True) self.b_fc = nn.Linear(256, 128 * 8 * 8) self.b_do = layers.AdditiveNoise(0.5)
self.e_unflatten = layers.ViewLayer([128, 8, 8]) self.e_in1 = nn.InstanceNorm2d(128, 0.8) self.e_us1 = nn.ConvTranspose2d(128, 128, 2, 2) self.e_conv1 = nn.Conv2d(128, 128, 3, stride=1, padding=1) self.e_in2 = nn.InstanceNorm2d(128, 0.8) self.e_lrelu1 = nn.LeakyReLU(0.2, inplace=True) self.e_us2 = nn.ConvTranspose2d(128, 128, 2, 2) self.e_conv2 = nn.Conv2d(128, 64, 3, stride=1, padding=1) self.e_in3 = nn.InstanceNorm2d(64, 0.8) self.e_lrelu2 = nn.LeakyReLU(0.2, inplace=True) self.e_conv3 = nn.Conv2d(64, 1, 3, stride=1, padding=1) self.e_tanh = nn.Tanh()
For both proposals we use the proposed discriminator with the following architecture.
class Discriminator(nn.Module): def __init__(self):
super(Discriminator, self).__init__() self.conv1 = nn.Conv2d(2, 16, 3, 2, 1)
self.lrelu1 = nn.LeakyReLU(0.2, inplace=True) self.conv2 = nn.Conv2d(16, 32, 3, 2, 1) self.lrelu2 = nn.LeakyReLU(0.2, inplace=True) self.in2 = nn.InstanceNorm2d(32, 0.8) self.conv3 = nn.Conv2d(32, 64, 3, 2, 1) self.lrelu3 = nn.LeakyReLU(0.2, inplace=True) self.in3 = nn.InstanceNorm2d(64, 0.8) self.conv4 = nn.Conv2d(64, 128, 3, 2, 1) self.lrelu4 = nn.LeakyReLU(0.2, inplace=True) self.in4 = nn.InstanceNorm2d(128, 0.8) self.flatten = layers.ViewLayer([128*2*2]) self.fc = nn.Linear(128*2*2, 1)
def forward(self, x, y): xy = torch.cat([x, y], dim=1) for module in self.children():
xy = module(xy) yx = torch.cat([y, x], dim=1) for module in self.children():
yx = module(yx) return F.softmax(torch.cat([xy, yx], dim=1), dim=1)

heading: C.9 ADDITIONAL FIGURES FOR MARKOV CHAIN PROPOSALS IN SAMPLE-BASED SETTING
text: In this section, we show additional figures for Markov chain proposals. In Fig. 12 we show samples from the chain that was initialized by the noise. In Fig. 13 we show samples from the chain after 10000 accepted samples.

Review summaries:
One of the main contributions of the paper is showing how maximizing the acceptance rate in Metropolis Hastings (MH) translates in minimizing the symmetric KL divergence of the target and the proposal distribution. The result aligns with the intuition that if this variance is 0, detailed balanced is satisfied and, hence, we can always accept the proposal. Also Equation 11 nicely fits the intuition that a good (independent) proposal should minimize the KL divergence between proposal and target (as in VI) under the constraint that the proposal has full support compared to the target distribution, which is enforced by the last term. Theorem 1 and its proof are great.