Measuring the reliability of MCMC inference with bidirectional Monte Carlo
Markov chain Monte Carlo (MCMC) is one of the main workhorses of probabilistic inference, but it is notoriously hard to measure the quality of approximate posterior samples. This challenge is particularly salient in black box inference methods, which can hide details and obscure inference failures. In this work, we extend the recently introduced bidirectional Monte Carlo [GGA15] technique to evaluate MCMC-based posterior inference algorithms. By running annealed importance sampling (AIS) chains both from prior to posterior and vice versa on simulated data, we upper bound in expectation the symmetrized KL divergence between the true posterior distribution and the distribution of approximate samples. We integrate our method into two probabilistic programming languages, WebPPL [GS]  and Stan [CGHL+ p], and validate it on several models and datasets. As an example of how our method be used to guide the design of inference algorithms, we apply it to study the effectiveness of different model representations in WebPPL and Stan.
The authors propose a generic method for evaluation of MCMC transition operators based on a recently proposed method known as bidirectional Monte Carlo (BDMC). The BDMC method performs two annealed importance sampling (AIS) sweeps, which are essentially the same as those used in the Tempered Transitions (TT) method [Neal Stat & Comp 1996], but reinterprets the accumulated probabilities as a lower and upper bound on the marginal likelihood. This paper proposes to use the gap between these two bounds as a generic performance measure for MCMC transition operators. The presented technical contributions are implementations of BDMC for the Stan and WebPPL probabilistic programming systems, and a protocol for evaluation of quality of posterior inference.  
This paper has some strong points and some not so strong points. The main strong point is that using BDMC to assess convergence of MCMC operators is a beautifully simple idea, and easy to implement, which in my opinion means that this work is potentially high impact. This is particularly true in the context of probabilistic programming systems, which indeed are the envisioned use case here, and I think all such systems would do well to at least implement this method.   My main point of concern lie with the novelty and the depth of the technical contribution. The authors cite an arxiv submission on BDMC as existing work, but (I think wisely) choose to devote a relatively large amount of space to reiterating its description. Unfortunately this does mean that the main technical contributions presented in sections 3.1 and 3.2 are somewhat rushed, and it is unfortunately also here where the writing quality slips a bit. As a result some issues of clarity arise (see below).  Clarity issues aside, the technical contributions in this paper are somewhat modest. Relative to the existing work on BDMC, the main novelty appears to the realization that the difference between the lower and upper bound is itself provides an upper bound on the Jensen-Shannon divergence. This is indeed a very nice observation, though arguably not an entirely novel once.   This brings me a point about related work. The authors, somewhat curiously, do not cite Radford Neal's work on Tempered Transitions [Statistics and Computing (1996) 6, 353-366], which is highly relevant in this context. This citation is also missing from the manuscript on BDMC, which makes me wonder whether the authors are in fact aware of it. If so, this is unfortunate in the sense that there is a missed opportunity in pointing out connections between the two methods.  TT performs the same reverse "heating" and forward "annealing" sweep as BDMC, but with the interest of a "tempering" MCMC transtion move to improve mixing in multimodal distributions. In TT the final sample from the reverse sweep is used as the initial sample for the forward sweep, and the final sample from the forward sweep is then used as a proposal which is accepted or rejected according to the ratio of importance weights accumulated in both sweeps.   The acceptance ratio in TT is, as far a I can tell, precisely the ratio of the stochastic lower and upper bound. Section 3.3 of Neal's paper in fact presents an interpretation of this acceptance ratio in terms of two thermodynamic integration estimates that bound log Z_T - log Z_1 from above and below, noting that a TT sampler will accept with probability 1 when the two converge.   The analogy with TT also leads me to wonder why the authors do spend more attention on the choice of β_t values. In TT it is well known that too large steps in Δβ = β_{t+1}-β_t will lead to a poor acceptance ratio, wheras too small steps Δβ waste computation. The ratio of the stochastic lower and upper bound in BDMC will also depend on Δβ and this provides some very clear intuitions. If your two bounds are not close enough, you need to rerun BDMC with smaller Δβ. Similarly, I suspect that comparisons between MCMC transition operators are likely to be unreliable when the gap between the log upper and lower bounds is larger >1 nat. When this is the case, one should again lower Δβ.   Ιt hasn't quite become clear to me from reading the paper what the "auto" part in AutoBDMC is, but I think there is an opportunity here to come up with some schedule for iteratively decreasing Δβ until sufficiently reliable inference results are obtained. In this respect the experiments could be a little more thorough. It would be good to see a systematic study of how the lower and upper bound vary with the step size.   These criticisms aside, I think this is not a bad paper. As noted above, I think the method proposed by the authors is potentially quite useful. Depending on whether the author's response can convincingly argue that issues with clarity and discussion of related work can be addressed by camera ready, I could be inclined to argue for acceptance.     Clarity  - The definition of $\vec{v}$ as "all of the variables sampled in AIS before the final stage" is rather loose. The authors presumably mean ${z_{1:T-1}, \theta_{1:T-1}$, but this is not immediately unambiguous since there are two AIS sweeps in the algorithm. This is not helped by the use of inconsistent labels "rev" and "fwd" on the one hand and "for" and "back" on the other (is this difference significant, or a mere accident?)  - The authors do not define what estimator \hat{p}(y) they are taking the expectation of in equations (5-7) and (9-11), which makes it difficult to follow how these equations arise. Do the authors mean to write E[log \hat{R}_{rev}] and [log \hat{R}_{fwd}]?   - Equation (6), as written, makes no sense (the distributions inside the KL divergence are identical). Based on equations (4) and (10), I take it the right terms should read something like p(z, θ | y) q_back(v | z, θ, y)?  - As far as I can tell the authors never unambiguously define what they mean by AutoBDMC. This is presumably means (a) caculating the difference between the stochastic lower and upper bound in order to assess convergence and (b) using the fixed hyperparameter scheme described in section 3.2?   Experiments  - It is now clear how β_t is chosen. Looking at the code in the appendix I see `var step = 1/options.steps`. I assume that AutoBDMC is therefore parameterized by the number of steps T to go from β_T = 1 to β_1 = 1 / T? It would be helpful to state this somewhere. In figures 2 and 3, does "HMC / No-U-Turn steps" indicate that BDMC is performed with this number of steps?   - The acronym RAIS, which presumably stands for reverse AIS, is not defined in the caption, nor are the variants (u)c(R)AIS, which I am guessing refer to (un)collapsed variants?  - It would be helpful to increase the size of the axes and tick labels.   Typos  - Two closing parens in the D_KL terms are missing in the expression for D_JS - In the NUTS citation: Homan -> Hoffman