 # Hamiltonian Monte Carlo Acceleration Using Surrogate Functions with Random Bases

For big data analysis, high computational cost for Bayesian methods often limits their applications in practice. In recent years, there have been many attempts to improve computational efficiency of Bayesian inference. Here we propose an efficient and scalable computational technique for a state-of-the-art Markov Chain Monte Carlo (MCMC) methods, namely, Hamiltonian Monte Carlo (HMC). The key idea is to explore and exploit the structure and regularity in parameter space for the underlying probabilistic model to construct an effective approximation of its geometric properties. To this end, we build a surrogate function to approximate the target distribution using properly chosen random bases and an efficient optimization process. The resulting method provides a flexible, scalable, and efficient sampling algorithm, which converges to the correct target distribution. We show that by choosing the basis functions and optimization process differently, our method can be related to other approaches for the construction of surrogate functions such as generalized additive models or Gaussian process models. Experiments based on simulated and real data show that our approach leads to substantially more efficient sampling algorithms compared to existing state-of-the art methods.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Bayesian statistics has provided a principled and robust framework to create many important and powerful data analysis methods over the past several decades. Given a probabilistic model for the underlying mechanism of observed data, Bayesian methods properly quantify uncertainty and reveal the landscape or global structure of the parameter space. However, these methods tend to be computationally intensive since Bayesian inference usually requires the use of Markov Chain Monte Carlo (MCMC) algorithms to simulate samples from intractable distributions. Although the simple Metropolis algorithm  is often effective at exploring low-dimensional distributions, it can be very inefficient for complex, high-dimensional distributions: successive states may exhibit high autocorrelation, due to the random walk nature of the movement. As a result, the effective sample size tends to be quite low and the convergence to the true distribution is usually very slow. The celebrated Hamiltonian Monte Carlo (HMC) [3, 4]

reduces the random walk behavior of Metropolis by Hamiltonian dynamics, which uses gradient information to propose states that are distant from the current state, but nevertheless have a high probability of acceptance.

Although HMC explores the parameter space more efficiently than random walk Metropolis does, it does not fully exploit the structure (i.e., geometric properties) of parameter space  since dynamics are defined over Euclidean space. To address this issue, Girolami and Calderhead  proposed a new method, called Riemannian Manifold HMC (RMHMC), that uses the Riemannian geometry of the parameter space  to improve standard HMC’s efficiency by automatically adapting to local structures.

To make such geometrically motivated methods practical for big data analysis, one needs to combine them with efficient and scalable computational techniques. A common bottleneck for using such sampling algorithms for big data analysis is repetitive evaluations of functions, their derivatives, geometric and statistical quantities that involves the whole observed data and maybe a complicated model. A natural question is how to construct effective approximation of these quantities that provides a good balance between accuracy and computation cost. One common approach is subsampling (see, for example, [7, 8, 9, 10]), which restricts the computation to a subset of the observed data. This is based on the idea that big datasets contain a large amount of redundancy so the overall information can be retrieved from a small subset. However, in general applications, we cannot simply use random subsets for this purpose: the amount of information we lose as a result of random sampling leads to non-ignorable loss of accuracy, which in turn has a substantially negative impact on computational efficiency . Therefore, in practice, it is a challenge to find good criteria and strategies for an appropriate and effective subsampling.

Another approach is to exploit smoothness or regularity in parameter space, which is true for most statistical models. This way, one could find computationally cheaper surrogate functions to substitute the expensive target (or potential energy) functions [15, 16, 18, 25, 19, 20]. However, the usefulness of these methods is often limited to moderate dimensional problems because of the computational cost needed to achieve desired approximation accuracy.

In this work, our objective is to develop a faster alternative to the method of . To this end, we propose to use random nonlinear bases along with efficient learning algorithms to construct a surrogate functions that provides effective approximation of the probabilistic model based on the collective behavior of the large data. The randomized nonlinear basis functions combined with the computationally efficient learning process can incorporate correct criteria for an efficient implicit subsampling resulting in both flexible and scalable approximation [29, 30, 31, 32]. Because our method can be presented as a special case of shallow random networks implemented in HMC, we refer to it as random network surrogate function; however, we will show that our proposed method is related to (and can be extended to) other surrogate functions such as generalized additive models and Gaussian process models by constructing the surrogate functions using different bases and optimization process.

Our proposed method provides a natural framework to incorporate surrogate functions in the sampling algorithms such as HMC, and it can be easily extended to geometrically motivated methods such as Riemannian Manifold HMC. Further, for problems with a limited time budget, we propose an adaptive version of our method that substantially reduces the required number of training points. This way, the random bases surrogate function could be utilized earlier and its approximation accuracy could be improved adaptively as more training points become available. We show that theoretically the learning procedure for our surrogate function is asymptotically equivalent to potential matching, which is itself a novel distribution matching strategy similar to the score matching method discussed in [35, 20].

Finally, we should emphasize that out method is used to generate high quality proposals at low computational cost. However, when calculating the acceptance probability of these proposals, we use the original Hamiltonian (used in standard HMC) to ensure that the stationary distribution of the Markov chain will remain the correct target distribution.

Our paper is organized as follows. An overview of HMC and RMHMC is given in Section 2. Our random network surrogate HMC is explained in detail in Section 3. The adaptive model is developed in Section 4. Experimental results based on simulated and real data are presented in Section 5. Code for these examples is available at github.com/chengzhang-uci/RNSHMC. Finally, Section 6 is devoted to discussion and future work. As an interesting supplement, an overall description of the potential matching procedure is presented in Section B in the appendix.

## 2 Preliminaries

### 2.1 Hamiltonian Monte Carlo

In Bayesian Statistics, we are interested in sampling from the posterior distribution of the model parameters given the observed data, ,

 P(q|Y)∝exp(−U(q)), (1)

where the potential energy function is defined as

 U(q)=−N∑i=1logP(yi|q)−logP(q). (2)

The posterior distribution is almost always analytically intractable. Therefore, MCMC algorithms are typically used for sampling from the posterior distribution to perform statistical inference. As the number of parameters increases, however, simple methods such as random walk Metropolis  and Gibbs sampling  may require a long time to converge to the target distribution. Moreover, their explorations of parameter space become slow due to inefficient random walk proposal-generating mechanisms, especially when there exist strong dependencies between parameters in the target distribution. By inclusion of geometric information from the target distribution, HMC [3, 4] introduces a Hamiltonian dynamics system with auxiliary momentum variables to propose samples of in a Metropolis framework that explores the parameter space more efficiently compared to standard random walk proposals. More specifically, HMC generates proposals jointly for and using the following system of differential equations:

 dqidt =∂H∂pi (3) dpidt =−∂H∂qi (4)

where the Hamiltonian function is defined as . The quadratic kinetic energy function

corresponds to the negative log-density of a zero-mean multivariate Gaussian distribution with the covariance

. Here,

is known as the mass matrix, which is often set to the identity matrix,

. Starting from the current state , the Hamiltonian dynamics system (3),(4) is simulated for steps using the leapfrog method, with a stepsize of . The proposed state, , which is at the end of the trajectory, is accepted with probability

. By simulating the Hamiltonian dynamics system together with the correction step, HMC generates samples from a joint distribution

 P(q,p)∝exp(−U(q)−12pTM−1p)

Notice that and are separated, the marginal distribution of then follows the target distribution. These steps are presented in Algorithm 1.

Following the dynamics of the assumed Hamiltonian system, HMC can generate distant proposals with high acceptance probability which allows an efficient exploration of parameter space.

### 2.2 Riemannian Manifold HMC

Although HMC explores the target distribution more efficiently than random walk Metropolis, it does not fully exploit the geometric structures of the underlying probabilistic model since a flat metric (i.e., ) is used. Using more geometrically motivated methods could substantially improves sampling algorithms’ efficiency. Recently, Girolami and Calderhead  proposed a new method, called Riemannian Manifold HMC (RMHMC), that exploits the Riemannian geometry of the target distribution to improve standard HMC’s efficiency by automatically adapting to local structures. To this end, instead of the identity mass matrix commonly used in standard HMC, they use a position-specific mass matrix . More specifically, they set to the Fisher information matrix, and define Hamiltonian as follows:

 H(q,p)=U(q)+12logdetG(q)+12pTG(q)−1p=ϕ(q)+12pTG(q)−1p (5)

where . Note that standard HMC is a special case of RMHMC with . Based on this dynamic, they propose the following HMC on Riemmanian manifold:

 ˙q=∇pH(q,p)=G(q)−1p˙p=−∇qH(q,p)=−∇qϕ(q)+12ν(q,p) (6)

where the

th element of the vector

is

 (ν(q,p))i=−pT∂i(G(q)−1)p=(G(q)−1p)T∂iG(q)G(q)−1p

with the shorthand notation for partial derivative.

The above dynamic is non-separable (it contains products of and ), and the resulting proposal generating mechanism based on the standard leapfrog method is neither time-reversible nor symplectic. Therefore, the standard leapfrog algorithm cannot be used for the above dynamic . Instead, we can use the Stömer-Verlet  method, known as generalized leapfrog ,

 p(t+1/2) = p(t)−ϵ2[∇qϕ(q(t))−12ν(q(t),p(t+1/2))] (7) q(t+1) = q(t)+ϵ2[G−1(q(t))+G−1(q(t+1))]p(t+1/2) (8) p(t+1) = p(t+1/2)−ϵ2[∇qϕ(q(t+1))−12ν(q(t+1),p(t+1/2))] (9)

The resulting map is 1) deterministic, 2) reversible, and 3) volume-preserving. However, it requires solving two computationally intensive implicit equations (Equations (7) and (8)) at each leapfrog step.

## 3 Random Network Surrogate HMC (RNS-HMC)

For HMC, the Hamiltonian dynamics contains the information from the target distribution through the potential energy and its gradient. For RMHMC, more geometric structure (i.e., the Fisher information) is included through the mass matrix for kinetic energy. It is the inclusion of such information in the Hamiltonian dynamics that allows HMC and RMHMC to improve upon random walk Metropolis. However, one common computational bottleneck for HMC and other Bayesian models for big data is repetitive evaluations of functions, their derivatives, geometric and statistical quantities. Typically, each evaluation involves the whole observed data. For example, one has to compute the potential and its gradient from the equation (2) for HMC and mass matrix , its inverse and the involved partial derivatives for RMHMC at every time step or M-H correction step. When is large, this can be extremely expensive to compute. In some problems, each evaluation may involve solving a computationally expensive problem. (See the inverse problem and Remark in Section 5.2.)

To alleviate this issue, in recent years several methods have been proposed to construct surrogate

Hamiltonians. For relatively low dimensional spaces, (sparse) grid based piecewise interpolative approximation using precomputing strategy was developed in

. Such grid based methods, however, are difficult to extend to high dimensional spaces due to the use of structured grids. Alternatively, we can use Gaussian process model, which are commonly used as surrogate models for emulating expensive-to-evaluate functions, to learn the target functions from early evaluations (training data) [16, 19, 25]. However, naive (but commonly used) implementations of Gaussian process models have high computation cost associated with inverting the covariance matrix, which grows cubically as the size of the training set increases. This is especially crucial in high dimensional spaces, where we need large training sets in order to achieve a reasonable level of accuracy. Recently, scalable Gaussian processes using induing point methods [23, 24] have been introduced to scale up GPs to larger datasets. While these methods have been quite succsessful in reducing computation cost, the tuning of inducing points could still be problematic in high dimensional spaces. (See a more detailed discussion in Section 3.3.)

The key in developing surrogate functions is to develop a method that can effectively capture the collective properties of large datasets with scalability, flexibility and efficiency. In this work, we propose to construct surrogate functions using proper random nonlinear bases and efficient optimization process on training data. More specifically, we present our method as a special case of shallow neural networks; although, we show that it is related to (and can be extended to) other surrogate functions such as generalized additive models and Gaussian process models. Random networks of nonlinear functions prove capable of approximating a rich class of functions arbitrarily well

[30, 32]. Using random nonlinear networks and algebraic learning algorithms can also be viewed as an effective implicit subsampling with desired criteria. The choice of hidden units (basis functions) and the fast learning process can be easily adapted to be problem specific and scalable. Unlike typical (naive) Gaussian process models, our random network scales linearly with the number of training points. In fact, a random nonlinear network can be considered as a standard regression model with randomly mapped features. For such shallow random networks, the computational cost for inference is cubic in the number of hidden nodes. Those differences in scaling allow us to explicitly trade off computational efficiency and approximation accuracy and construct more efficient surrogate in certain applications. As our empirical results suggest, with appropriate training data good approximation of smooth functions in high dimensional space can be achieved using a moderate and scalable number of hidden units. Therefore, our proposed method has the potential to scale up to large data sizes and provide effective and scalable surrogate Hamiltonians that balance accuracy and efficiency well.

### 3.1 Shallow Random Network Approximation

A typical shallow network architecture (i.e., a single-hidden layer feedforward scalar-output neural network) with

hidden units, a nonlinear activation function

, and a scalar (for simplicity) output for a given d-dimensional input is defined as

 z(q)=s∑i=1via(q;γi)+b (10)

where is the ith hidden node parameter, is the output weight for the ith hidden node, and is the output bias. Given a training data set

 T={(q(j),t(j))|q(j)∈Rd,t(j)∈R,j=1,…,N}

the neural network can be trained by finding the optimal model parameters to minimize the mean square error cost function,

 C(W|T)=1NN∑j=1∥z(q(j))−t(j)∥2 (11)

The most popular algorithm in machine learning to optimize (

11) is back-propagation . However, as a gradient descent-based iterative algorithm, back-propagation is usually quite slow and can be trapped at some local minimum since the cost function is nonlinear, and for most cases, non-convex. Motivated by the fact that randomization is computationally cheaper than optimization, alternative methods based on random nonlinear bases have been proposed [29, 31]

. These methods drastically decrease the computation cost while maintaining a reasonable level of approximation accuracy. The key feature of random networks is that they reduce the full optimization problem into standard linear regression by mapping the input data to a randomized feature space and then apply existing fast algebraic training methods (e.g., by minimizing squared error) to find the output weight. Given the design objective, algebraic training can achieve exact or approximate matching of the data at the training points. Compared to the gradient descent-based techniques, algebraic training methods can reduce computational complexity and provide better generalization properties. A typical algebraic approach for single-hidden layer feedforward random networks is extreme learning machine (ELM)

, which is summarized in Algorithm 2.

Using randomized nonlinear features, ELM estimates the output weight by finding the least-squares solution to the resulting linear equations system

. Note that presented this way, our method can also be regarded as a random version of Generalized Additive Model (GAM). In practice, people could add regularization to improve stability and generalizability.

### 3.2 Choice of Nonlinearity

There are many choices for nonlinear activation functions in random networks. Different types of activation functions can be used for different learning tasks. In this paper, we focus on random networks with two typical types of nonlinear nodes:

 a(q;γ)=a(w⋅q+d),w∈Rd,d∈R,γ={w,d}

where and are the weight vector and the bias of the hidden node.

• Radial basis functions (RBF) nodes:

 a(q;γ)=a(−∥q−c∥22ℓ2),c∈Rd,ℓ∈R+,γ={c,ℓ}

where and are the center and width of the hidden node.

Both random networks can approximate a rich class of functions arbitrarily well [30, 32]

. With randomly assigned input weights and biases composed linearly inside the nonlinear activation function, additive nodes form a set of basis functions, whose level sets are hyperplanes orientated by

and shifted by respectively. Random networks with additive nodes tend to reflect the global structure of the target function. On the other hand, RBF nodes are almost compactly supported (can be adjusted by the width ) rendering good local approximation for the corresponding random networks.

### 3.3 Connection to GPs and Sparse GPs

It is worth noting the connection between networks with RBF nodes and Gaussian processes models [21, 22]. Given a training data set

 T={(q(j),t(j))|q(j)∈Rd,t(j)∈R,j=1,…,N}

and using squared exponential covariance function

 K(q(j),q(j′))=σ2fexp(−∥q(j)−q(j′)∥22ℓ2),θ={σf,ℓ}

the standard GP regression with a Gaussian noise model has the following marginal likelihood

 p(t|Q,θ)=N(t|0,KN+σ2I)

where is the covariance matrix and is the noise parameter. Prediction on a new observation is made according to the conditional distribution

 p(t∗|q∗,T,θ)=N(t∗|kT∗(KN+σ2I)−1t,K∗∗−kT∗(KN+σ2I)−1k∗+σ2)

where and .

On the other hand, if we use as the th hidden node, , the output matrix becomes , and the output weight learned by algebraic approach to a regularized least square problem is

 ^v=argminv∥Hv−t∥2+σ2vTKNv=(KN+σ2I)−1t

Therefore, such a network provides the same prediction point estimate as the above full GP models. This way, a Gaussian process model can be interpreted as a self-organizing RBF network where new hidden nodes are added adaptively with new observations. This is also an alternative point of view to  where GP models were shown to be equivalent to single hidden-layer neural networks with infinite many hidden nodes.

Notice that the above GP model scales cubically with the number of data points which limits the application of GPs to relatively small datasets. To derive a GP model that is computationally tractable for large datasets, sparse GPs based on inducing point methods [23, 24] have been previously proposed. These sparse models introduce a set of inducing points and approximate the exact kernel by an approximation for fast computation. For example, the fully independent training conditional (FITC)  method uses the approximate kernel

 ~KFITC(q,q′)=kTqK−1Mkq′+δqq′(K(q,q′)−kTqK−1Mkq′)

where is the exact covariance matrix for the inducing points and are the exact covariance matrices between and the inducing points. Given the same training data set , the marginal likelihood is

 p(t|Q,θ)=N(t|0,KNMK−1MKMN+Λ+σ2I)

Here, is a diagonal matrix with that adjusts the diagonal elements to the ones from the exact covariance matrix . Simiarly, predictions can be made according to the conditional distribution

 p(t∗|q∗,T,θ)=N(t∗|kT∗Σ−1MKMN(Λ+σ2I)−1t,K∗∗−kT∗(K−1M−Σ−1M)k∗+σ2)

where . The computation cost is reduced to for learning and after that the predictive mean costs

per new observation. The hyperparameters

and inducing points can be tuned to maximize the marginal likelihood. However, in high dimension the tuning of becomes infeasible. Figure 1: Comparing different surrogate approximations with an increasing number of observations N=10,20,40 on target function y=x2/2. The observation points are nested samples from the standard normal distribution. For FITC and random networks, we choose 5 inducing points and 5 hidden neurons respectively. FITC and random networks are all run 100 times and averaged to reduce the random effects on the results.

On the other hand, if we use the inducing points as the centers of a RBF network, the output matrix . Given the diagonal matrix , the output weight estimated by the algebraic approach to a weighted least square problem plus a regularization term is

 ^v=argminv∥D−12(Hv−t)∥2+vTKMv=Σ−1MKMN(Λ+σ2I)−1t

Therefore, the same predictive mean can be reproduced if we use the inducing points as centers and use the same hyperparameter configuration in our random network with RBF nodes and optimize with respect to the above cost function. Moreover, those hyperparameters in each basis (or hidden nodes) can also be trained (typically by gradient decent) if we abandon the use of random bases and simple linear regression.

Figure 1 compares random networks with different node types and related GP methods on fitting a simple function which corresponds to the negative log-density of the standard normal distributions. We used softplus function in additive nodes and exponential square kernels in RBF nodes and GP methods. As we can see from the graph, random networks generally perform better than GP methods when the number of observations is small. The randomness in the configuration of hidden nodes force networks to learn more globally. In contrast, GP models are more local and need more data to generalize well. By introducing sparsity, FITC tends to generalize better than full GP, especially on small datasizes. Since our goal is to fit negative log-posterior density function in (2) and as moves away from the high density domain, using softplus basis functions in random networks are more capable to capture this far field feature by striking a better balance between flexibility and generalization while being less demanding on the datasize. Also, the number of hidden neurons (bases) can be used to regularize the approximation to mitigate overfitting issue.

### 3.4 Surrogate Induced Hamiltonian Flow Figure 2: Comparing HMC and NNS-HMC based on a 2-dimensional banana-shaped distribution. The left panel shows the gradient fields (force map) for the original Hamiltonian flow (red) and the surrogate induced Hamiltonian flow (blue). The middle and right panel show the trajectories for HMC and NNS-HMC samplers. Both samplers start from the same point (red square) with same initial momentums. Blue points at the end of the trajectories are the proposals. The overall acceptance probability drops from 0.97 using HMC to 0.88 using NNS-HMC.

As mentioned in the previous sections, repetitive computation of Hamiltonian, its gradient and other quantities that involve all data set undermine the overall exploration efficiency of HMC. To alleviate this issue, we exploit the smoothness or regularity in parameter space, which is true for most statistical models. As discussed in [14, 15, 16], one can improve computational efficiency of HMC by approximating the energy function and using the resulting approximation to device a surrogate transition mechanism while still converging to the correct target distribution. To this end,  proposed to use pre-convergence samples (which are discarded during the burn-in period) to approximate the energy function using a Gaussian process model. Here, we define an alternative surrogate-induced Hamiltonian as follows:

 ~H(q,p)=~U(q)+12pTM−1p

where is the neural network surrogate function. now defines a surrogate-induced Hamiltonian flow, parametrized by a trajectory length , which is a map . Here, being the end-point of the trajectory governed by the following equations

 dqdt=∂~H∂p=M−1p,dpdt=−∂~H∂q=−∂~U∂q

When the original potential is computationally costly, simulating the surrogate induced Hamiltonian system provides a more efficient proposing mechanism for our HMC sampler. The introduced bias along with discretization error from the leap-frog integrator are all naturally corrected in the MH step where we use the original Hamiltonian in the computation of acceptance probability. As a result, the stationary distribution of the Markov chain will remain the correct target distribution (see Appendix A for a detailed proof). Note that by controlling the approximation quality of the surrogate function, we can maintain a relatively high acceptance probability. This is illustrated in Figure 2 for a two-dimensional banana-shaped distribution .

To construct such a surrogate, the early evaluations of the target function during the early iterations of MCMC will be used as the training set based on which we can train a shallow random network using fast algebraic approaches, such as ELM (Algorithm 2). The gradient of the scalar output (see 10) for a network with additive hidden nodes, for example, then can be computed as

 ∂z∂q=s∑i=1via′(wi⋅q+di)wi (12)

which costs only computations. To balance the efficiency in computation and flexibility in approximation, and to reduce the possibility of overfitting, the number of hidden nodes need to be small as long as a reasonable level of accuracy can be achieved. In practice, this can be done by monitoring the resulting acceptance rate using an initial chain. Figure 3: Comparing the efficiency of our random network surrogates and Gaussian process surrogates on a challenging 32dimensional Gaussian target whose covariance matrix has an eigenvector (1,1,…,1)Twith a corresponding eigenvalue of 1.0, and all other eigenvalues are 0.01. We set the step size to keep the acceptance probability around 70% for HMC and use the same step size in all surrogate methods. For FITC and random networks, the number of inducing points and hidden neurons are all set to be 1000 to allow reasonably accurate approximation. We ran each algorithm ten times and plot the medians and 80% error bars.

Following , we propose to run our method, henceforth called random network surrogate Hamiltonian Monte Carlo (RNS-HMC, see Algorithm 3), in two phases: exploration phase and exploitation phase. During the exploration phase, we initialize the training data set with an empty set or some samples from the prior distribution of parameters. We then run the standard HMC algorithm for some iterations and collect information from the new states (i.e., accepted proposals). When we have sufficiently explored the high density domain in parameter space and collected enough training data (during the burn-in period), a shallow random network is trained based on the collected training set to form a surrogate for the potential energy function. The surrogate function will be used to approximate the gradient information needed for HMC simulations later in the exploitation phase.

As an illustrative example, we compare the performance of different surrogate HMC methods on a challenging Gaussian target density in dimensions (A lower dimensional case was used in ). The target density has 31 confined directions and a main direction that is times wider, and all variables are correlated. Both the full GPs and FITC methods are implemented using GPML package . The results are presented in Figure 3. Compared to the full GPs, FITC and random networks (with additive and RBF nodes) all scales linearly with the number of observations. Both random network surrogates can start with fewer training data. We also compare the efficiency of the surrogate induced Hamiltonian flows in terms of time normalized mean effective sample sizes (ESS). The efficiency of FITC and random networks all increases as the number of observations increase until no more approximation gain can be obtained (see the acceptance probability in the middle panel). However, the efficiency of full GP begin to drop before the model reaches its full capacity. That is because its predictive complexity also grows with the number of observations, which in turn diminishes the overall efficiency. Overall, the random network with additive nodes outperform other methods based on these example.

Our proposed method provides a natural framework to incorporate surrogate functions in HMC. Moreover, it can be easily extended to RMHMC. To this end, the Hessian matrix of the surrogate function can be used to construct a metric in parameter space and the third order derivatives can be used to simulate the corresponding modified Hamiltonian flow. We refer to this extended version of our method as RNS-RMHMC. Figure 4: Acceptance probability of surrogate induced Hamiltonian flow on simulated logistic regression models for different number of parameters, d, and hidden neurons, s. The step size is chosen to keep the acceptance probability around 70% for HMC. Left: Acceptance probability as a function of s (x-axis) and d (y-axis). Middle: Acceptance probability as a function of s for a fixed dimension d=32. Right: Acceptance probability as a function of d for a fixed s=2000.

Note that the approximation quality of the neural network surrogate function depends on several factors including the dimension of parameter space, , the number of hidden neurons, and the training size, . Here, we assume that is sufficiently large enough, and investigate the efficiency of RNS-HMC in terms of its acceptance probability for different values of and based on a standard logistic regression model with simulated data. Similar to the results presented in , Figure 4 shows the acceptance rate (over MCMC runs) as a function of and . For dimensions up to , RNS-HMC maintains a relatively high acceptance probability with a shallow random network trained in a few seconds on a laptop computer. Appendix B provides a theoretical justification for our method.

So far, we have assumed that the neural network model in our method is trained using a sufficiently large enough number of training points after waiting for an adequate number of iterations to allow the sampler explore the parameter space. This, however, could be very time consuming in practice: waiting for a long time to collect a large number of training points could undermine the benefits of using the surrogate Hamiltonian function.

Figure 5 shows the average acceptance probabilities (over MCMC chains) as a function of the number of training points, , and the number hidden neurons, , on a simulated logistic regression model for a fixed number of parameters, . While it takes around training data points to fulfill the network’s capability and reach a high acceptance comparable to HMC, only training points are enough to provide an acceptable surrogate Hamiltonian flow (around acceptance probability). Therefore, we can start training the neural network surrogate earlier and adapting it as more training points become available. Although adapting a Markov chain based on its history may undermine its ergodicity and consequently its convergence to the target distribution , we can enforce a vanishing adaption rate such that and and update the surrogate function with probability at iteration . By Theorem 5 of , the resulting algorithm is ergodic and converges to the right target distribution. Figure 5: Acceptance probability of the surrogate induced Hamiltonian flow based on a simulated logistic regression models with dimension d=32. Left: Acceptance probability as a function of the number of hidden neurons s (x-axis) and the number of training points N (y-axis). Middle: Acceptance probability as a function of N for a fixed s=1000 . Right: Acceptance probability as a function of s for a fixed N=1600 .

It is straightforward to adapt the neural network surrogate from the history of Markov chain. However, the estimator in Algorithm 2 costs computation and storage, where is the number of training data (e.g., the number of rows in the output matrix ). As increases, finding the right weight and bias for the output neuron becomes increasingly difficult. In , Grevillle shows that can be learned incrementally in real time as new data points become available. Based on Greville’s method, online and adaptive pseudoinverse solutions for updating ELM weights has been proposed in  which can be readily employed here to develop an adaptive version of RNS-HMC. To be more efficient, only the estimator is updated.

###### Proposition 1

Suppose the current output matrix is and the target vector is . At time , a new sample and the target (potential) are collected. Denote the output vector from the hidden layer at as . The adaptive updating formula for the empirical potential matching estimator is given by

 vk+1=vk+(tk+1−hTk+1vk)bk+1

where and auxiliary matrices are updated according to .

•  bk+1=Θkhk+11+hTk+1Θkhk+1,Φk+1=Φk,Θk+1=Θk−Θkhk+1bTk+1
•  bk+1=ck+1∥ck+1∥2,Φk+1=Φk−Φkhk+1bTk+1
 Θk+1=Θk−Θkhk+1bTk+1+(1+hTk+1Θkhk+1)bk+1bTk+1−bk+1hTk+1ΘTk

At time , the estimator takes a one-off computation and storage (only need to store and , not ). Starting at a previously computed solution , and two auxiliary matrices , this adaptive updating costs computation and storage, independent of the training data size . Further details are provided in Appendix C. We refer to this extended version of our method as Adaptive RNS-HMC (ARNS-HMC).

## 5 Experiments

In this section, we use several experiments based on logistic regression models and inverse problem for elliptic partial differential equation (PDE) to compare our proposed methods to standard HMC and RMHMC in terms of sampling efficiency defined as time-normalized effective sample size (ESS). Given

MCMC samples for each parameter, ESS = , where is the sum of monotone sample autocorrelations . We use the minimum ESS over all parameters normalized by the CPU time, (in seconds), as the overall measure of efficiency: . The corresponding step sizes and number of leapfrog steps for HMC and RMHMC are chosen to make them stable and efficient (e.g., reasonably high acceptance probability and fast mixing). The same settings are used for our methods. Note that while the acceptance rates are similar in the first two examples, they drop a little bit for the last two examples, which is mainly due to the constraints we imposed on our surrogate functions. To prevent non-ergodicity and ensure high ESS for both HMC and RNS-HMC, we follow the suggestion by  to uniformly sample from in each iteration. The number of hidden nodes in random network surrogates are not tuned too much and better results could be obtained by more careful tunings.

In what follows, we first compare different methods in terms of their time-normalized ESS after the burin-period. To this end, we collect samples after a reasonable large number of iterations () of burn-in to make sure the chains have reached their stationary states. For our methods, the accepted proposals during the burn-in period after a short warming-up session (the first 1000 iterations) are used as a training set for a shallow random network. Later, we show the advantages of our adaptive algorithm.

### 5.1 Logistic regression model

As our first example, we compare HMC, RMHMC, RNS-HMC, and RNS-RMHMC using a simulation study based on a logistic regression model with parameters and observations. The design matrix is and true parameters are uniformly sampled from , where . The binary responses

are sampled independently from Bernoulli distributions with probabilities

. We assume , and sample from the corresponding posterior distribution.

Notice that the potential energy function is now a convex function, the Hessian matrix is positive semi-definite everywhere. Therefore, we use the Hessian matrix of the surrogate as a local metric in RNS-RMHMC. For HMC and RNS-HMC, we set the step size and leapfrog steps . For RMHMC and RNS-RMHMC, we set the step size and leapfrog steps .

To illustrate that our method indeed converges to the right target distribution, Figure 6 provides the one- and two-dimensional posterior marginals of some selected parameters obtained by HMC and RNS-HMC. Table 1 compares the performance of the four algorithms. As we can see, RNS-HMC has substantially improved the sampling efficiency in terms of time-normalized min(ESS). Figure 6: HMC vs RNS-HMC: Comparing one- and two-dimensional posterior marginals of β1,β11,β21,β31,β41 based on the logistic regression model with simulated data.

Next, we apply our method to two real datasets: Bank Marketing and the dataset . The Bank Marketing dataset (40197 observations and 24 features) is collected based on direct marketing campaigns of a Portuguese banking institution aiming at predicting if a client will subscribe to a term deposit . We set the step size and number of leapfrog steps for HMC and RNS-HMC; for RMHMC and RNS-RMHMC. The dataset (32561 features and 123 features) is complied from the UCI dataset  which has been used to determine whether a person makes over 50K a year. We reduce the number of features to by random projection (increasing the dimension to 100 results in a substantial drop in the acceptance probability). We set the step size and number of leapfrog steps for HMC and RNS-HMC;

for RMHMC and RNS-RMHMC. All datasets are normalized to have zero mean and unit standard deviation. The priors are the same as before. The results for the two data sets are summarized in Table 1. As before, both RNS-HMC and RNS-RMHMC significantly outperform their counterpart algorithms.

### 5.2 Elliptic PDE inverse problem

Another computationally intensive model is the elliptic PDE inverse problem discussed in [42, 43]. This classical inverse problem involves inference of the diffusion coefficient in an elliptic PDE which is usually used to model isothermal steady flow in porous media. Let be the unknown diffusion coefficient and be the pressure field, the forward model is governed by the elliptic PDE

 ∇x⋅(c(x,θ)∇xu(x,θ))=0, (13)

where is the spatial coordinate. The boundary conditions are

 u(x,θ)|x2=0=x1,u(x,θ)|x2=1=1−x1,∂u(x,θ)∂x1∣∣x1=0=0,∂u(x,θ)∂x1∣∣x1=1=0

In our numerical simulation, (13) is solved using standard continuous GFEM with bilinear basis functions on a uniform quadrilateral mesh.

A log-Gaussian process prior is used for with mean zero and an isotropic squared-exponential covariance kernel:

 C(x1,x2)=σ2exp(−∥x1−x2∥222l2)

for which we set the variance

and the length scale . Now, the diffusivity field can be easily parameterized with a Karhunen-Loeve (K-L) expansion:

 c(x,θ)≈exp(d∑i=1θi√λivi(x))

where and

are the eigenvalues and eigenfunctions of the integral operator defined by the kernel

, and the parameter are endowed with independent standard normal priors, , which are the targets of inference. In particular, we truncate the K-L expansion at modes and condition the corresponding mode weights on data. Data are generated by adding independent Gaussian noise to observations of the solution field on a uniform grid covering the unit square.

 yj=u(xj,θ)+ϵj,ϵj∼N(0,0.12),j=1,2,…,N

The number of leap frog steps and step sizes are set to be for both HMC and NNS-HMC. Note that the potential energy function is no longer convex; therefore, we can not construct a local metric from the Hessian matrix directly. However, the diagonal elements

 ∂2U∂θ2i=1σ2θ+N∑j=11σ2y(∂uj∂θi)2−N∑j=1ϵjσ2y∂2uj∂θ2i,σθ=0.5,σy=0.1,i=1,2,…,d

are highly likely to be positive in that the deterministic part (first two terms) is always positive and the noise part (last term) tends to cancel out. The diagonals of the Hessian matrix of surrogate therefore induce an effective local metric which can be used in RNS-RMHMC. A comparison of the results of all algorithms are presented in Table 1. As before, RNS-HMC provides a substantial improvement in the sampling efficiency. For the RMHMC methods, we set . As seen in the table, RMHMC is less efficient than HMC mainly due to the slow computation speed. However, RNS-RMHMC improves RMHMC substantially and outperforms HMC. Although the metric induced by the diagonals of the Hessian matrix of surrogate may not be as effective as Fisher information, it is much cheaper to compute and provide a good approximation.

Remark. In addition to the usual computational bottleneck as in previous examples, e.g., large amount of data, there is another challenge on top of that for this example due to the complicated forward model. Instead of a simple explicit probabilistic model that prescribes the likelihood of data given the parameter of interest, a PDE (13) is involved in the probabilistic model. The evaluation of geometrical and statistical quantities, therefore, involves solving a PDE similar to (13) in each iteration of HMC and RHMHC. This is a preventive factor in practice. Using our methods based on neural network surrogates provide a huge advantage. Numerical experiments show a gain of efficiency by more than 20 times. More improvement is expected as the amount of data increases. Figure 7: Median acceptance rate of ANNS-HMC along with the corresponding 90% interval (shaded area). The red line shows the average acceptance rate of standard HMC. Figure 8: Relative error of mean as a function of running time.

Next, using the above four examples we show that ARNS-HMC can start with far fewer training points and quickly reach the same level of performance as that of RNS-HMC. Figure 7 shows that as the number of training points (from initial MCMC iterations) increases, ARNS-HMC fully achieves the network’s capability and reaches a comparable acceptance rate to that of HMC.

We also compare ARNS-HMC to HMC and RNS-HMC in terms of the relative error of mean (REM) which is defined as , where means sample mean up to time . Figure 8 shows the results using the four examples discussed above. Note that before training the neural network models, both RNS-HMC and ARNS-HMC are simply standard HMC so the three algorithms have similar performance. As we can see, ARNS-HMC has the best overall performance: it tends to provide lower REM at early iterations. This could be useful if we have limited time budget to fit a model.

## 6 Discussion and Future Work

In this paper, we have proposed an efficient and scalable computational method for Bayesian inference by exploring and exploiting regularity of probability models in parameter space. Our method is based on training surrogate function of the potential energy after exploring the parameter space sufficiently well. For situations where it is not practical to wait for a thorough exploration of parameter space, we have proposed an adaptive version of our method that can start with fewer training points and can quickly reach its full potential.

As an example, we used random networks and efficient learning algorithms to construct effective surrogate functions. These random bases surrogate functions provide good approximations of collective information of the full data set while striking a good balance between accuracy and computation cost for efficient computation. Random networks combined with the optimized learning process can provide flexibility, accuracy, and scalability. Note that in general the overall performance could be sensitive to the architecture of the random network. Our proposed random network surrogate method scales differently than GP emulators because of the specific constraints we imposed on its architecture. As our experimental results show, this approach could improve the performance of HMC in some applications.

In its current form, our method is more effective in problems with costly likelihood and a moderate number of parameters. In spite of improvements we have made to standard HMC, dealing with high dimensional and complex distributions still remains quite challenging. For multimodal distributions, for example, our method’s effectiveness largely depends on the quality of training samples. If these samples are collected from one mode only, the surrogate function will miss the remaining modes and the sampler might not be able to explore them (especially if they are isolated modes). A surrogate function based on Gaussian processes might have a better chance at finding these modes in the tails of the approximate distribution since it tends to go to zero gradually. To address this issue, we can utilize mode searching and mode exploring ideas such as those proposed by [48, 49]. For constrained target distributions, we can employ the method of  based on Spherical HMC.

For HMC, gradient of the potential function is an important driving force in the Hamiltonian dynamics. Although accurate approximation of a well sampled smooth function automatically leads to accurate approximation of its gradient, this is not the case when the sampling is not well distributed. For example, when dense and well sampled training data sets are difficult to obtain in very high dimensions, one can incorporate the gradient information in the training process. In future, we will study more effective way to utilize this information in the training process. As a common practice in adaptive MCMC methods, one may also learn the mass matrix adaptively together with the surrogate in ARNS-HMC.

#### Acknowledgments

This work is supported by NIH grant R01AI107034 and NSF grants DMS-1418422 and DMS-1622490

## References

•  N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. (1953) Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 21, pp. 1087–1092.
•  S. Geman and D. Geman. (1984) Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6: 721–741
•  S. Duane, A. D. Kennedy, B J. Pendleton, and D. Roweth. (1987) Hybrid Monte Carlo. Physics Letters B, 195, pp. 216 – 222.
•  R. M. Neal. (2011) MCMC using Hamiltonian dynamics, in Handbook of Markov Chain Monte Carlo, S. Brooks, A. Gelman, G. Jones, and X. L. Meng, eds., Chapman and Hall/CRC, pp. 113–162.
•  M. Girolami and B. Calderhead. (2011) Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society, Series B, (with discussion), 73, pp. 123–214.
•  S. Amari and H. Nagaoka. Methods of Information Geometry, volume 191 of Translations of Mathematical monographs, Oxford University Press, 2000.
•  M. Welling and Y. W. Teh. (2011) Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th International Conference on Machine Learning (ICML), pages 681–688.
•  M. D. Hoffman, D. M. Blei, F. R. Bach. (2010) Online Learning for Latent Dirichlet Allocation. In Neural Information Processing Systems (NIPS), pp. 856-864.
•  B. Shahbaba, S. Lan, W.O. Johnson, and R.M. Neal. (2014) Split Hamiltonian Monte Carlo. Statistics and Computing, 24, pp. 339–349.
•  T. Chen, E. B. Fox, and C. Guestrin. (2014) Stochastic gradient Hamiltonian Monte Carlo. Preprint.
•  M. J. Betancourt. (2015) The Fundamental Incompatibility of Hamiltonian Monte Carlo and Data Subsampling, preprint: arxiv:1502.01510.
•  L. Verlet. (1967) Computer ’experiments’ on classical fluids, I: Thermodynamical properties of Lennard-Jones molecules. it Phys. Rev, 159, pp. 98–103
•  B. Leimkuler and S. Reich. (2004) Simulating Hamiltonian Dynamics, Cambridge: Cambridge University Press.
•  R. Neal. Bayesian Learning in Neural Networks, Springer, 1996.
•  J. S. Liu. (2001) Monte Carlo Strategies in Scientific Computing, New York: Springer
•  C. E. Rasmussen. (2003) Gaussian processes to speed up hybrid monte carlo for expensive bayesian integrals. Bayesian Statistics, 7, pp. 651–659.
•  C. E. Rasmussen and H. Nickisch. (2010) Gaussian processes for machine learning (GPML) toolbox. The Journal of Machine Learning Research, 11, pp. 3011–3015.
•  C. Zhang, B. Shahbaba and H. Zhao. (2015) Precomputing Strategy for Hamiltonian Monte Carlo Method Based on Regularity in Parameter Space. Arxiv:1504.01418v1, Apr 2015. Unpublished manuscript.
•  S. Lan, T. Bui, M. Christie, M. Girolami. (2015) Emulation of Higher-Order Tensors in Manifold Monte Carlo Methods for Bayesian Inverse Problems, arXiv:1507.06244.
•  H. Strathmann, D. Sejdinovic, S. Livingstone, Z. Szabo, A. Gretton. (2015) Gradient-free Hamiltonian Monte Carlo with Efficient Kernel Exponential Families, arXiv:1506.02564.
•  C. E. Rasmussen. Evaluation of Gaussian Processes and Other Methods for Non-linear Regression, PhD thesis, University of Toronto, 1996.
•  R. Neal. (1999) Regression and classification using Gaussian process priors. Bayesian Statistics, 6, pp. 475–501
•  E. Snelson and Z. Ghahramani. (2006) Sparse Gaussian processes using pseudo-input. In Advances in neural information processing systems (NIPS), 18, Cambridge, MA. MIT Press.
•  J. Quinonero-Candela and C. E. Rasmussen. (2005) A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research (JMLR), 6, pp. 1939–1959
•  E. Meeds and M. Welling. (2014) GPS-ABC: Gaussian Process Surrogate Approximate Bayesian Computation, arXiv:1401.2838.
•  K. Hornik, M. Stinchcombe and H. White. (1989) Multilayer feedforward networks are universal approximators. Neural Networks, 2, pp. 349–358
•  D. E. Rumelhart, G. E. Hinton and R. J. Williams, (1986) Learning internal representations by error propagation. Parallel Distributed Processing: Explorations in the Microstructure of Cognition, I. Cambridge, MA: Bradford Books, pp. 318–362
•  S. Ferrari and R. F. Stengel. (2005) Smooth function approximation using neural networks. IEEE Trans. Neural Networks, 16(1), pp. 24–38
•  G. B. Huang, Q. Y. Zhu and C. K. Siew. (2006) Extreme learning machine: Theory and applications. Neurocomputing, 70(1-3), pp. 489–501
•  G. B. Huang, L. Chen and C. K. Siew. (2006) Universal approximation using incremental constructive feedforward networks with random hidden nodes. IEEE Trans. Neural Networks, 17(4), pp. 879–892
•  A. Rahimi and B. Recht. (2007) Random features for large-scale kernel machines. In Proceedings of the 21st conference on Advances in Neural Information Processing Systems (NIPS)
•  A. Rahimi and B. Recht. (2008) Uniform approximation of functions with random bases. Proc. 46th Ann. Allerton Conf. Commun., Contr. Comput.
•  P. L. Bartlett. (1998) The Sample complexity of pattern classification with neural networks: the size of the weights is more important than the size of the network. IEEE Trans. Inf. Theory, 44(2), pp. 525–536
•  C. J. Geyer. (1992) Practical Markov Chain Monte Carlo. Statistical Science, 7, pp. 473–483.
•  A. Hyvärinen. (2005) Estimation of non-normalized statistical models by score matching. JMLR, 6, pp. 695–709.
•  C. Andrieu and J.  Thoms. (2008) A tutorial on adaptive MCMC. Statistics and computing, 18(4), pp. 343–373.
•  G. O. Roberts and J. S. Rosenthal. (2007) Coupling and ergodicity of adaptive Markov chain Monte Carlo algorithm. J. Appl. Probab., 44(2), pp. 458–475.
•  T. Greville. (1960) Some applications of the pseudoinvese of a matrix. SIAM Rev., 2(1), pp. 15–22.
•  A. van Schaik and J. Tapson. (2015) Online and adaptive pseudoinverse solutions for ELM weights. Neurocomputing, 149, pp. 233–238.
•  T. Kohonen. (1988) Self-Organization and Associative Memory. Springer-Verlag
•  P. Kovanic. (1979) On the pseudo inverse of a sum of symmetric matrices with applications to estimation. Kybernetika, 15(5), pp. 341-348.
•  M.  Dashti and A. Stuart. (2011) Uncertainty Quantification and Weak Approximation of an Elliptic Inverse Problem. SIAM Journal on Scientific Computing, 49, pp. 2524–2542
•  P. R. Conard, Y. M. Marzouk, N. S. Pillai, and A. Smith. (2014) Asymptotically exact mcmc algorithms via local approximations of computationally intensive models, preprint arXiv:1402.1694v3, Nov 2014. Unpublished manuscript.
•  S. Moro, P. Cortez and P. Rita. (2014) A Data-Driven Approach to Predict the Success of Bank Telemarketing. Decision Support Systems, Elsevier, 62, pp. 22–31.
•  C. J. Lin, R. C. Weng and S. S. Keerthi. (2008) Trust region Newton method for large-scale logistic regression. Journal of Machine Learning Research, 9, pp. 627–650
•  K. Bache and M. Lichman. (2013) UCI machine learning repository. URL http:archive.ics.uci.edu/ml
•  W. R. Gilks, G. O. Roberts and S. K. Sahu. (1998) Adaptive Markov chain Monte Carlo through regeneration. Journal of the American Statistical Association, 93, pp. 1045–1054.
•  S. Ahn, Y. Chen, and M. Welling. Distributed and adaptive darting Monte Carlo through regenerations. In

Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AI Stat)

, 2013.
•  S. Lan, J. Streets and B. Shahbaba. (2014) Wormhole Hamiltonian Monte Carlo. In Twenty-Eighth AAAI Conference on Artificial Intelligence
•  S. Lan, B. Zhou, and B. Shahbaba. Spherical Hamiltonian Monte Carlo for constrained target distributions. In Proceedings of the 31th International Conference on Machine Learning (ICML), (2014).

## Appendix A Convergence to the correct target distribution

In order to prove that the equilibrium distribution remains the same, it suffices to show that the detailed balance condition still holds. Denote . In the Metropolis-Hasting step, we use the original Hamiltonian to compute the acceptance probability

 α(θ,θ′)=min(1,exp[−H(θ′)+H(θ)])

therefore,

 α(θ,θ′)P(dθ) =α(θ,θ′)exp[−H(θ)]dθ θ=~ϕ−1t(θ′)=min(exp[−H(θ)],exp[−H(θ′)])∣∣∣dθdθ′∣∣∣dθ′ =α(θ′,θ)exp[−H(θ′)]dθ′ =α(θ′,θ)P(dθ′)

since due to the volume conservation property of the surrogate induced Hamiltonian flow . Now that we showed the detailed balance condition is satisfied, along with the reversibility of the surrogate induced Hamiltonian flow, the modified Markov chain will converge to the correct target distribution.

## Appendix B Potential Matching

In the paper, training data collected from the history of Markov chain are used without a detailed explanation. Here, we will analyze the asymptotical behavior of surrogate induced distribution and try to explain why the history of the Markov chain is a reasonable choice for training data. Recall that we find our neural network surrogate function by minimizing the mean square error (11). Similarly to , it turns out that minimizing (11) is asymptotically equivalent to minimizing a new distance between the surrogate induced distribution and the underlying target distribution, independent of their corresponding normalizing constants.

Suppose we know the density of the underlying intractable target distribution up to a constant

 P(q|Y)=1Zexp(−U(q))

where is the unknown normalizing constant. Our goal is to approximate this distribution using a parametrized density model, also known up to a constant,

 Q(q,θ)=1Z(θ)exp(−V(q,θ))

Ignoring the multiplicative constant, the corresponding potential energy functions are and respectively. The straightforward square distance between the two potentials will not be a well-defined measure between the two distributions distributions because of the unknown normalizing constants. Therefore, we use the following distance instead:

 K(θ) =mind∫∥V(q,θ)−U(q)−d∥2P(q|Y)dq =∫∥V(q,θ)−U(q)∥2P(q|Y)dq−[Eq(V(θ)−U)]2=Varq(V(θ)−U) (14)

Because of its similarity to score matching , we refer to the approximation method based on this new distance as potential matching; the corresponding potential matching estimator of is given by

 ^θ=argminθK(θ)

It is easy to verify that , so is a well-defined squared distance. Exact evaluation of (B) is analytically intractable. In practice, given samples from the target distribution