AUXILIARY VARIATIONAL MCMC
We introduce Auxiliary Variational MCMC, a novel framework for learning MCMC kernels that combines recent advances in variational inference with insights drawn from traditional auxiliary variable MCMC methods such as Hamiltonian Monte Carlo. Our framework exploits low dimensional structure in the target distribution in order to learn a more efficient MCMC sampler. The resulting sampler is able to suppress random walk behaviour and mix between modes efficiently, without the need to compute gradients of the target distribution. We test our sampler on a number of challenging distributions, where the underlying structure is known, and on the task of posterior sampling in Bayesian logistic regression. Code to reproduce all experiments is available at https://github.com/AVMCMC.
This paper proposes an auxiliary variable MCMC scheme involving variational inference for efficient MCMC. Given a target distribution p(x), the authors introduce an auxiliary variable a, and learn conditional distributions p(a|x) and q(a|x) by minimizing the KL divergence between p(x)p(a|x) and q(a)q(x|a), with q(a) something simple (the authors use Gaussian). A MH proposal step involves simulating x givea the current MCMC sample x (from p(a|x), taking a step in A-space, and then returning back to the X space (using q(x|a)).  The authors show how to calculate the acceptance probability. 
I think the idea is nice and useful (I'm surprised people haven't thought of this before), though I think the paper presents this in a less clear way (as an extension of ideas from Agakov and Barber's "Auxiliary variational method"). While this is correct and perhaps more general, in my mind it slightly obscures the main idea, as well as the strong ties with variational autoencoders: express a complex distribution as a (learnt) transformation of a simple distribution (this is the actual approach taken in the experiments). 
The motivation of the approach is that the nonlinear encoding network can transform the complex p(x) into a simpler q(a). 
For this reason, I think an important baseline is the independent MH sampler from equation 8 (I think this essentially uses a trained VAE generative model as a proposal distribution). The authors talk about how producing independent proposals can be sub-optimal, yet it seems to me that if the encoder and decoder neural networks are powerful enough, this should do a good job. I think excluding this baseline hurts the paper a bit.
The proof of correctness while correct is a bit unclear, can perhaps be simplified if you view the MCMC algorithm as operating on an augmented space (x,a,x') with stationary distribution p(x)q(a|x)q(x'|a) (writing writing q for \tilde(q)). This clearly has the right distribution over x. Each MCMC iteration starts with x and proceeds as follow:
  1) Given x, sample a and x' from q(a|x) and q(x'|a)
  2) Make a deterministic proposal on the augmented space to swap (x,x'). The acceptance probability is now equation 2.
  3) Discard a,x'.
In figure 4, the authors use HMC as an "improved MCMC algorithm", yet this is not an algorithm that deals with multimodality well. More useful would be to include some tempering algorithm like serial or parallel tempering.
While I like the idea, I unfortunately don't think the experiments are very convincing (and the authors barely discuss their results). Other than mixture of Gaussians, HMC (which involves no training) appears to be superior. With some tempering, I expect it to outperform the proposed method for the MoG case
Table 2 left: since HMC involves no training, does this mean that, taking training time into account, HMC is 5-6 orders of magnitude more efficient. L?ke I mentioned earlier, these results need more discussion. 
It would also help to provide absolute training and run times, so the reader can better understand whether the proposed method of ANICE is better.
Figure 3: why don't the authors also plot the histogram of values in the auxiliary space, p(a). It would be interesting to see how Gaussian this is (this is what variational inference is trying to achieve). Also, does Figure 3(a) mean that conditioned on x, p(a|x) is basically a delta function? This would suggest that the encoder is basically learning a deterministic transformation to a simpler low-dimensional space? There is some work in this direction in the statistics literature, e.g. 
"Variable transformation to obtain geometric ergodicity in the random-walk Metropolis algorithm"
The authors some refers to the distribution of a|x as q(a|x) sometimes (in section 2.1) and sometimes as p(a|x) which is a bit confusing.
Figure 2: the labels are wrong.