We present Newtonian Monte Carlo (NMC), a method to improve Markov Chain Monte Carlo (MCMC) convergence by analyzing the first and second order gradients of the target density to determine a suitable proposal density at each point. Existing first order gradient-based methods suffer from the problem of determining an appropriate step size. Too small a step size and it will take a large number of steps to converge, while a very large step size will cause it to overshoot the high density region. NMC is similar to the Newton-Raphson update in optimization where the second order gradient is used to automatically scale the step size in each dimension. However, our objective is not to find a maxima but instead to find a parameterized density that can best match the local curvature of the target density.   This parameterized density is then used as a single-site Metropolis-Hastings proposal.

 As a further improvement on first order methods, we show that random variables with constrained supports don't need to be transformed before taking a gradient step. NMC directly matches constrained random variables to a proposal density with the same support thus keeping the curvature of the target density intact.

 We demonstrate the efficiency of NMC on a number of different domains. For statistical models where the prior is conjugate to the likelihood, our method recovers the posterior quite trivially in one step. However, we also show results on fairly large non-conjugate models, where NMC performs better than adaptive first order methods such as NUTS or other inexact scalable inference methods such as Stochastic Variational Inference or bootstrapping.
 Markov Chain Monte Carlo (MCMC) methods are often used to generate samples from an unnormalized probability density π(θ) that is easy to evaluate but hard to directly sample. Such densities arise quite often in Bayesian inference as the posterior of a generative model p(θ, Y ) conditioned on some observations Y = y, where π(θ) = p(θ, y). The typical setup is to select a proposal distribution q(.|θ) that proposes a move of the Markov chain to a new state θ * ∼ q(.|θ). This Metropolis-Hastings acceptance rule is then used to accept or reject this move with probability: min 1, π(θ * )q(θ|θ * ) π(θ)q(θ * |θ) . When θ ∈ R k , a common proposal density is the Gaussian distribution N (θ, 2 I k ) centered at θ with covariance 2 I k , where is the step size and I k is the identity matrix defined over R k,k . This proposal forms the basis of the so-called Random Walk MCMC (RWM) first proposed in Metropolis et al. (1953) . In cases where the target density π(θ) is differentiable, an improvement over the basic RWM method is to propose a new value in the direction of the gradient, as follows: This method is known as Metropolis Adjusted Langevin Algorithm (MALA), and arises from an Euler approximation of a Langevin diffusion process (Robert and Tweedie, 1996) . MALA has been shown to reduce the number of steps required for convergence to O(n 1/3 ) from O(n) for RWM (Roberts and Rosenthal, 1998 ). An alternate approach, which also uses the gradient, is to do an L-step Euler approximation of Hamiltonian dynamics known as Hamiltonian Monte Carlo (Neal, 1993) , although it was originally published under the name Hybrid Monte Carlo (Duane et al., 1987) . In HMC the number of steps, L, can be learned dynamically by the No-U-Turn Sampler (NUTS) algorithm (Hoffman and Gelman, 2014) . However, in all three of the above algorithms -RWM, MALA, and HMC -there is an open problem of selecting the optimal step size. Normally, the step size is adaptively learned by targeting a desired acceptance rate. This has the unfortunate effect of picking the same step size for all the dimensions of θ, which forces the step size to accomodate the dimension with the smallest variance as pointed out in Girolami and Calderhead (2011) . The same paper introduces alternate approaches, using Reimann manifold versions of MALA (MMALA) and HMC (RMHMC). They propose a Reimann manifold using the expected Fisher information matrix plus the negative Hessian of the log-prior as a metric tensor, −E y|θ ∂ 2 ∂θ 2 log{p(y, θ)} , and proceed to derive the Langevin diffusion equation and Hamiltonian dynamics in this manifold. The use of the above metric tensor does address the issue of differential scaling in each dimension. However, the method as presented requires analytic knowledge of the Fisher information matrix. This makes it difficult to design inference techniques in a generic way, and requires derivation on a per-model basis. A more practical approach involves using the negative Hessian of the log-probability as the metric tensor, ∂ 2 ∂θ 2 log{p(y, θ)}. However, this encounters the problem that this is not necessarily positive definite throughout the state space. An alternate approach for scaling the moves in each dimension is to use a preconditioning matrix M (Roberts and Stramer, 2002) in MALA, q(.|θ) = N θ + 2 M ∇ log{π(θ)}, 2 M , also known as the mass matrix in HMC and NUTS, but it's unclear how to compute this. An alternate approach is to approximately compute the Hessian (Zhang and Sutton, 2011) using ideas from quasi-Newton optimization methods such as L-BFGS (Nocedal and Wright, 2006) . This approach and its stochastic variant (Simsekli et al., 2016 ) use a fixed window of previous samples of size M to approximate the Hessian. However, this makes the chain an order M Markov chain, which introduces considerable complexity in designing the transition kernel in addition to introducing a new parameter M . The key observation in our work is that for single-site methods we only need to compute the Hessian of one coordinate at a time, and this is usually tractable. The other key observation is that we don't need to always make a Gaussian proposer using the Hessian. In some cases, other densities which are less concentrated such as Cauchy are more appropriate. In general, the Hessian can be used for the purpose of matching the curvature of any parameterized density that best approximates the conditional posterior. This approach of curvature-matching to an approximating density allows us to deal with constrained random variables without introducing a transformation such as in Stan (Carpenter et al., 2017) . In the rest of the paper, we will describe our approach to exploit the curvature of the target density, and show some results on multiple data sets. We have presented a novel MCMC method that uses the curvature of the target density to converge faster than existing state of the art methods, and without requiring any adaptive tuning. As next steps, we will fully integrate NMC into a production PPL and evaluate its performance across a wider spectrum of illustrative and real-world use cases.
Exploit curvature to make MCMC methods converge faster than state of the art.