Newtonian Monte Carlo: single-site MCMC meets second-order gradient methods

Single-site Markov Chain Monte Carlo (MCMC) is a variant of MCMC in which a single coordinate in the state space is modified in each step. Structured relational models are a good candidate for this style of inference. In the single-site context, second order methods become feasible because the typical cubic costs associated with these methods is now restricted to the dimension of each coordinate. Our work, which we call Newtonian Monte Carlo (NMC), is a method to improve 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 to find a parameterized proposal density rather than the maxima. As a further improvement on existing first and second order methods, we show that random variables with constrained supports don't need to be transformed before taking a gradient step. 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.


Introduction
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(.|θ). The Metropolis-Hastings acceptance rule is then used to ac-Copyright c 2019, Association for the Advancement of Artificial Intelligence (www.aaai.org). All rights reserved. cept 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 Riemann manifold versions of MALA (MMALA) and HMC (RMHMC). They propose a Riemann 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.
Another 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 cases when the coordinate under consideration is a constrained random variable then we can propose from any parameterized density in the same constrained space by matching its curvature. 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.

Overview
This paper introduces the Newtonian Monte Carlo (NMC) technique for sampling from a target distribution via a proposal distribution that incorporates curvature around the current sample location. We wish to choose a proposal distribution that uses second order gradient information in order to closely match the target density. Whereas related MCMC techniques discussed in Section 1 may utilize second order gradient information, those techniques typically use it only to adjust the step size when simulating steps along the general direction of the target density's gradient.
Our proposed method involves matching the target density to a parameteric density that best explains the current state. We have a library of 2-parameter target densities F i , and simple inference rules such that, given the first and sec-ond order gradients, we can solve the following two equations: to determine α i and β i . For example, in the case of θ ∈ R k , we use either the multivariate Gaussian or the multivariate Cauchy. For the former, the update equation leads to the natural proposal, The update term in the mean of this multivariate Gaussian is precisely the update term of the Newton-Raphson Method (Whittaker and Robinson, 1967), which is where NMC gets its name from. However, if the negative Hessian inverse matrix is not positive definite, then the multivariate normal is not defined. In this case, we can instead use the Cauchy proposer as shown by Minka (2000). The full list of estimation methods are enumerated in Section 3. For example, for positive real values we use a Gamma proposer with parameters, and we don't need a log-transform to an unconstrained space. We rely on generic Tensor libraries such as PyTorch (Paszke et al., 2017) that make it easy to write statistical models and also automatically compute the gradients. This makes our approach easy to apply to models generically.
In the case of conjugate models, our estimation methods automatically recover the appropriate conditional posterior distribution, such as the ones used in BUGS (Spiegelhalter et al., 1996). However, even in cases of non-conjugacy, our proposal distributions pick out reasonable approximations to the conditional posterior of each variable.

Single Site Inference
An important observation related to our method is that we don't need to compute the joint Hessian of all the parameters in the latent space. Most statistical models with relational structure can be decomposed into multiple latent variables. This decomposition allows for single site MCMC methods that change the value of one variable at a time. In this case, we only need to compute the gradient and Hessian of the target density w.r.t. the variable being modified. Consider a model with N variables each drawn from R K . The full Hessian is of size (N K) 2 and has a cost of (N K) 3 to invert. On the other hand, a single site approach computes N Hessians each of size K 2 with a total cost of N K 3 to invert.

Estimation Rules
The estimation rules presented here are based on the work of Minka (2000).

Unconstrained spaces
We first consider distributions with the support R k .

Normal Distribution
The multivariate Normal distribution has the log-density: This leads to the natural estimation rule: In case the estimated Σ has a negative eigenvalue we set those negative eigenvalues to a very small positive number, and reconstruct Σ.

Cauchy Distribution
The multivariate Cauchy distribution has the log-density: , and Noting that the second term above is the outer product of the first gradient leads to the following estimation rules:

Constrained Spaces
Half-Space We use the Gamma distribution for proposing values for variables which lie on any half-space constrained distribution, i.e. R + . The Gamma distribution has the logdensity: Which leads to the estimation rules: α = 1 − x 2 ∇ 2 log π(x), β = −x∇ 2 log π(x) − ∇ log π(x), Simplexes The K-simplexes refers to the set {x ∈ We use the Dirichlet distribution to propose random variables with this support. The log-density of the Dirichlet is given by, We consider the modified density, which includes the simplex constraint, Thus, Which leads to the following estimation rule,

Models
Neal's Funnel We first consider a toy model, which is considered difficult for MCMC methods. The model was first proposed in Neal and others (2003) and has been since called Neal's Funnel. The following equations define the joint density of the model, which is deceptively simple.
z ∼ N (0, 3) x ∼ N (0, e z 2 ) The difficulty for inference arises when we try to sample values of (x, z) for increasingly negative z. Since the scale of x varies exponentially with z it is hard to learn a good scale. Indeed Figure 1 shows that Stan, which uses NUTS, has a hard time sampling from this distribution as highlighted by the posterior marginal of z. On the other hand NMC, which effectively computes the scale dynamically at each point has no difficulty in generating good samples, as shown in Figure 2.
Bayesian Logistic Regression Next we consider the Bayesian Logistic Regression model that is commonly used in machine learning. The model is defined as follows: α ∼ N (0, 10, size = 1), β ∼ N (0, 2.5, size = K), X i ∼ N (0, 10, size = K) ∀i ∈ 1..N From this model we generate samples of α, β, X i and Y i for some given N and K. Half of the X, Y samples are given to the inference method to produce posterior samples of α and β. The held out values of X, Y are used to evaluate the predictive log-likelihood, which is averaged over the posterior samples. In this and the rest of the models, we draw exactly 1 000 samples using each method. For single-site methods a sample includes an update to each coordinate. Figure 3 shows results for N = 20 000 using a variety of methods. Table 1 shows the run times to produce 1 000 samples, plus the number of samples required to achieve convergence. We have defined samples to convergence to be the number of samples it takes for the predictive log-likelihood to stabilize to within 1% of the final value for the method. Figure 4 shows results for N = 2 000 000, where we only compare two of the methods since the other methods were too slow for such a large data set.
This model has a nice log concave posterior which is quite easy for MCMC inference and hence both JAGS and NMC converge in 1 sample. Stan, which uses NUTS, does take a bit longer to converge because it has to learn the optimum scale in each of the K dimensions of β. Bootstrappingbased approaches are often used for this model, but they do appear to incur additional cost of re-training the model for each bootstrap sample. Finally, we also used Stochastic Variational Inference as implemented in Pyro (Bingham et al., 2019), but it doesn't appear from this example that the loss of accuracy of using variational inference is worthwhile. In this model, NMC seems to be the fastest both in terms of time per samples and samples to convergence. Robust Regression Robust regression is a regression model in which an error distribution with a much wider tail than Gaussian such as the Student's t distribution is used to model data with outliers. We use the following model: As before we generate samples from the model of all the variables including N values of X i and Y i . Half of the generated samples are given to the inference algorithm and the other half are used for evaluating the posterior. This model is not log-concave because of the Student's t distribution, and as a result JAGS takes much longer to converge as shown in Figure 5 and and Table 2. We also ran an experiment for much larger N , Figure 6, where we left out JAGS because it was too slow. On this model, NMC converges much faster than Stan using merely 3 samples to converge for the larger data set, and with much faster runtimes as well. Annotation Model Our final model has a lot more relational structure than the regression models above. This is a slightly modified version of the model presented in Passonneau and Carpenter (2014) and Dawid and Skene (1979), and is designed to compute the true labels of items given noisy crowd-sourced labels. There are N items, K labelers, and each item could be one of C categories. Each item i is labeled by a set J i of labelers. Such that the size of J i is sampled randomly, and each labeler in J i is drawn uniformly without replacement from the set of all labelers. z i is the true label for item i and y ij is the label provided to item  i by labeler j. Each labeler l has a confusion matrix θ l such that θ lmn is the probability that an item with true class m is labeled n by l.
Here α m ∈ R + C . We set α mn = γ · ρ if m = n and α mn = γ · (1 − ρ) · 1 C−1 if m = n. Where γ is the concentration and ρ is the a-priori correctness of the labelers. In this model, Y il and J i are observed.
In our experiments, we fixed K = 100, C = 3, J loc = 2.5, γ = 10, and ρ = 0.5. As before we generated data for different sizes of N and used a random partition of the data for inference and evaluation. Since Stan doesn't support discrete variables such as z i above, we had to analytically integrate 1 the z's out of the model. For the purpose of evaluation, since Stan doesn't give us samples of z, we integrate over these variables to compute the predictive likelihood, which gives a disadvantage to methods such as JAGS and NMC where the samples depend on a specific value of z.
In this model each random variable has a conjugate conditional posterior, and since JAGS is designed to exploit conjugacy it really shines here. Unfortunately, the version of JAGS that we used kept crashing on larger data sets. The results for N = 10 000 are shown in Figure 7 and run times are in Table 3. NMC exploits the relational structure in this model to use single-site inference on each of the random variables such as θ lm individually rather than the entire θ jointly as in Stan. As such NMC is easily able to keep up with JAGS in terms of number of samples and is only a factor of 2,5 slower on the small data set. On the larger data set, Figure 8 and 9, NMC is nearly 7 times faster than Stan.

Conclusion
We have presented a new MCMC method that uses the first and second gradients of the target density for each coordinate in the latent state space to determine an appropriate proposal distribution. The method is shown to perform better than the existing state of the art NUTS implementation without requiring an adaptive phase or tuning of inference hyper-parameters.