Variational Hybrid Monte Carlo for Efficient Multi-Modal Data Sampling

Abstract

The Hamiltonian Monte Carlo (HMC) sampling algorithm exploits Hamiltonian dynamics to construct efficient Markov Chain Monte Carlo (MCMC), which has become increasingly popular in machine learning and statistics. Since HMC uses the gradient information of the target distribution, it can explore the state space much more efficiently than random-walk proposals, but may suffer from high autocorrelation. In this paper, we propose Langevin Hamiltonian Monte Carlo (LHMC) to reduce the autocorrelation of the samples. Probabilistic inference involving multi-modal distributions is very difficult for dynamics-based MCMC samplers, which is easily trapped in the mode far away from other modes. To tackle this issue, we further propose a variational hybrid Monte Carlo (VHMC) which uses a variational distribution to explore the phase space and find new modes, and it is capable of sampling from multi-modal distributions effectively. A formal proof is provided that shows that the proposed method can converge to target distributions. Both synthetic and real datasets are used to evaluate its properties and performance. The experimental results verify the theory and show superior performance in multi-modal sampling.

1. Introduction

Generating samples from many distributions encountered in Bayesian inference and machine learning is difficult. Markov Chain Monte Carlo (MCMC) provides a robust framework to generate samples from complex target distributions. Through constructing specific Markov chains, the MCMC methods can converge to the correct target distribution with the chains evolving. Presently, MCMC plays an essential role in artificial intelligence applications and probability inference, especially for estimating the expectations of the target functions [1].

Sampling methods based on dynamics are one of the most popular MCMC methods. The most commonly used dynamics in MCMC are Langevin dynamics and Hamiltonian dynamics. Hamiltonian Monte Carlo (HMC) [2,3] has become one of the most popular MCMC algorithms in Bayesian inference and machine learning. Unlike previous MCMC algorithms [4], HMC takes advantage of the gradient information to explore the continuous probability density function (PDF), which makes HMC more efficient to converge to the target distribution. In particular, HMC transforms the PDF into the potential energy function and adds the kinetic energy function to simulate the motion of the particle in the particular phase space, and thus HMC can satisfy the ergodic property. In practice, HMC exploits the Hamiltonian equation to calculate the new state of the proposed points in the phase space. To maintain a detailed balance, the Metropolis–Hasting (MH) technique is adopted [5]. Gradient information helps to discover and explore the phase space more efficiently, and there is much further research on how to better leverage gradients for HMC [6,7].

Nevertheless, samplers based on dynamics still have some deficiencies. The traditional dynamics samplers [8] and their developments [7,9,10,11,12] have excellent performance in unimodal distributions. However, when facing multi-modal distributions, these algorithms may reveal some problems, especially when the modes are far away from each other. When the modes are close to each other, the momentum variable in dynamics samplers may offer chances for the sample to jump into different modes. When the modes are isolated, the momentum variable cannot jump out of the current mode, because the interval between the two modes has tremendous potential energy. Generally, objects tend to stay in low-energy places which are also referred as high-probability regions. Although we can enlarge the momentum variable to pass through high-potential-energy places, the momentum variable should be exponential-order large, which causes a rapid decrease in the performance of the samplers. To deal with the problem of multi-modal sampling, several studies have been developed [13,14,15,16]. The authors in [13] introduced temperature in the simulation, and simulated annealing is used to gradually reduce the temperature from a high initial value to the value at which we wish to sample. Parallel tempering [17] is one variant of the simulated annealing, which performs inter-distribution exchange with multiple Markov chains to make the sampler explore the state space more freely, without changing the resulting sample distribution. The authors in Sminchisescu and Welling [14] proposed a new dynamics sampler which is based on a darting algorithm [18]. Umbrella sampling [19] is another precise sampling method, which divides the reaction coordinate space into several windows with bias potential to calculate the unbiased free energy in each window. However, when the dimensions are high, these algorithms may have elaborate set-ups and large computational complexity, resulting in low efficiency. The authors in Lan et al. [15] used the natural gradient of the target distribution to establish paths between different modes, and thus samples can jump through the low-probability regions. This method may suffer from low effective sample size (ESS) [3], which means that the independence between two neighbor samples is low. The authors in Tripuraneni et al. [16] introduced the concept of magnetic fields. By means of constructing a dynamics system based on a magnetic field, this method can achieve great performance in multi-modal sampling. However, the setting of the magnetic field parameter is difficult, and this method may also suffer from low ESS in multi-modal sampling problems.

In this paper, we introduce a novel dynamics MCMC method called the variational hybrid Monte Carlo (VHMC). We first improve Hamiltonian dynamics through Langevin dynamics to reduce the autocorrelation of samples and accelerate the convergence of the dynamics sampler. Furthermore, we exploit the variational distribution [20] of the target distribution to help the dynamics sampler to find the new mode. A new Metropolis–Hasting criterion is proposed to satisfy the detailed balance condition [5]. We find that VHMC can overcome the distant multi-modal sampling problem since dynamics-based samplers can sample unimodal distributions well and variational distribution guides the dynamics-based sampler to jump between different modes. Finally, detailed proof is given to demonstrate that our algorithm can converge to the target distribution.

Both synthetic data and real data experiments are conducted to verify our theory. We sample points from seven different Gaussian mixture distributions whose dimensions range from 2 to 256. We apply our method to two-class classification exploiting Bayesian logistic regression [21] to test the performance of VHMC on real datasets. Evaluation indices such as maximum mean discrepancy [22] and autocorrelation are calculated to assess the quality of samples. Experiment results illustrate that the proposed method is capable of sampling from distant multi-modal distribution while obtaining better performance compared with other state-of-the-art methods [16,23].

The main contributions of this work can be summarized as follows. We propose a novel sampler called Langevin Hamiltonian Monte Carlo (LHMC), which achieves lower autocorrelation and faster convergence compared with the HMC sampler. Since we introduce random factors in LHMC, the total energy of the system changes during the simulation, so we design a new Metropolis–Hasting procedure to keep the detailed balance. In addition, to improve the poor performance of the LHMC sampler in multi-modal sampling, we propose a new method VHMC, which uses the variational distribution of the target distribution to guide the sampler to jump through different modes. We use Adam [24] to find the local points with the highest probability density value, and we use these points to construct a mixture of Gaussian as the variational distribution of the target distribution. Detailed proof is given to prove the correctness of our method.

The rest of this article is organized as follows. In Section 2, we review the preliminaries of our study, including Hamiltonian Monte Carlo and Langevin dynamics. Then, we introduce our LHMC sampler and show the objective function in Section 3. In Section 4, we propose the variational hybrid Monte Carlo, which aims to address the problem of multi-modal sampling. Experiments are reported in Section 5. Discussion and conclusions are summarized in Section 6 and Section 7.

2. Preliminaries

In this section, we introduce the basic methods exploited in our method. First, we introduce Hamiltonian Monte Carlo, a widely used MCMC sampler in Bayesian machine learning. Then, Langevin dynamics is introduced to improve the performance of the HMC sampler.

2.1. Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (HMC) [2,3] is one of the state-of-the-art Markov chain Monte Carlo algorithms. The use of gradient information of the target distribution makes HMC more efficient than the traditional Metropolis–Hasting (MH) algorithms, which employ the random-walk scheme to explore the state space.

HMC exploits Hamiltonian dynamics to calculate the new state, whose state space is composed of joint Gaussian momentum p and position $\theta$, where variable p is independent with variable $\theta$. Suppose we use a spring oscillator to describe the Hamiltonian dynamics. Then, we can obtain the following equation:

$$H(\theta ,p)=U\left(\theta \right)+K\left(p\right),$$

(1)

where $U\left(\theta \right)$ represents the potential energy of the ball at position $\theta$ at time t while $K\left(p\right)$ represents the kinetic energy of the ball at time t. $H(p,\theta )$ represents the total energy of the ball. To construct Hamiltonian dynamics, the derivatives of position $\theta$ and momentum p about time are required. The Hamiltonian equations are formed as follows:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\theta }{\mathrm{d}t}& =\frac{\partial H(p,\theta )}{\partial p}={\nabla }_{p}K\left(p\right)\hfill \\ \hfill \frac{\mathrm{d}p}{\mathrm{d}t}& =-\frac{\partial H(p,\theta )}{\partial \theta }=-{\nabla }_{\theta }U\left(\theta \right).\hfill \end{array}$$

(2)

In practice, it is difficult to obtain the exact solutions to these differential equations. HMC instead discretizes these equations by using non-zero time steps, which inevitably introduces some error. It is, nevertheless, necessary to use a discretization for which Liouville’s theorem holds exactly [8]. The common discretization method of HMC is leapfrog, which takes the following form:

$$\begin{array}{cc}\hfill p_{t+\frac{ϵ}{2}}& =p_t-\frac{ϵ}{2}{\nabla }_{\theta }U\left({\theta }_t\right)\hfill \\ \hfill {\theta }_{t+ϵ}& ={\theta }_t+ϵ{\nabla }_{p}K\left(p_{t+\frac{ϵ}{2}}\right)\hfill \\ \hfill p_{t+ϵ}& =p_{t+\frac{ϵ}{2}}-\frac{ϵ}{2}{\nabla }_{\theta }U\left({\theta }_{t+ϵ}\right),\hfill \end{array}$$

(3)

where $ϵ$ represents the step size. Leapfrog preserves the phase space volume and is also time-reversible. Through the discretization method, we can obtain the new state, and HMC methods then apply Metropolis–Hasting to the new state to decide whether to accept or reject the state, which takes the following form:

$$\begin{array}{c}\hfill \min \left[1,\exp \left(H({\theta }_{n-1},p_{n-1})-H(\widehat{\theta },\widehat{p})\right)\right],\end{array}$$

(4)

where $({\theta }_{(n-1)},p_{(n-1)})$ represents the last state and $(\widehat{\theta },\widehat{p})$ represents the newly proposed state. By means of controlling the leapfrog size L and small step $ϵ$, we can adjust the acceptance rate of the HMC sampler. Algorithm 1 gives the pseudo-code of HMC [3].

Suppose we need to sample from the distribution of $\theta$ given the observation data $\mathcal{D}$:

$$p\left(\theta \right|\mathcal{D})\propto e^{-U\left(\theta \right)},$$

(5)

where we have the form of potential energy:

$$U\left(\theta \right)\propto -\mathrm{l}np\left(\theta \right|\mathcal{D}).$$

(6)

Algorithm 1: Hamiltonian Monte Carlo [2]

Input: Step size $ϵ$, leapfrog size L, starting point ${\theta }^{\left(1\right)}$, sample number N Output: Samples ${\theta }_{(1:N)}$ for $n=1$ to N do     $p_{\left(n\right)}\sim \mathcal{N}(0,1)$     $({\theta }_0,p_0)=({\theta }_{\left(n\right)},p_{\left(n\right)})$     $p_0=p_0-\frac{ϵ}{2}{\nabla }_{\theta }U\left({\theta }_0\right)$     ${\theta }_0={\theta }_0+ϵ{\nabla }_{p}K\left(p_0\right)$     for $i=1$ to L do         $p_i=p_{i-1}-ϵ{\nabla }_{\theta }U\left({\theta }_{i-1}\right)$         ${\theta }_i={\theta }_{i-1}+ϵ{\nabla }_{p}K\left(p_i\right)$     end for     $p_L=p_L-\frac{ϵ}{2}{\nabla }_{\theta }U\left({\theta }_L\right)$     $(\widehat{\theta },\widehat{p})=({\theta }_L,p_L)$     $u\sim \mathrm{U}niform(0,1)$     $\alpha =\min \left[1,\exp \left((U\left({\theta }_{\left(n\right)}\right)+K\left(p_{\left(n\right)}\right))-(U\left(\widehat{\theta }\right)+K\left(\widehat{p}\right))\right)\right]$     if $\alpha >u$ then         $({\theta }_{(n+1)},p_{(n+1)})=(\widehat{\theta },\widehat{p})$     else         $({\theta }_{(n+1)},p_{(n+1)})=({\theta }_{\left(n\right)},p_{\left(n\right)})$     end if end for

    According to the Hamiltonian dynamics, through introducing a set of auxiliary momentum variables p, HMC sampler can sample the joint distribution $\pi (\theta ,p)$ defined as:

$$\pi (\theta ,p)\propto e^{-H(\theta ,p)}.$$

(7)

By using the leapfrog Equation (3), we obtain the new state of $\theta$ and p. Because the position variable $\theta$ and the momentum variable p are independent, sampling $\theta$ and p alternatively will not affect the results.

The Hamiltonian dynamics have three properties. First, it preserves total energy $H({\theta }_{\left(t\right)},p_{\left(t\right)})=H({\theta }_{\left(0\right)},p_{\left(0\right)})$, and thus the joint probability density has $p({\theta }_{\left(t\right)},p_{\left(t\right)})=p({\theta }_{\left(0\right)},p_{\left(0\right)})$. Second, it satisfies volume conservation. Last, it is time-reversible [25]. As a result, if the potential energy and the kinetic energy remain unchanged during the dynamics system, then the joint probability density of $\theta$ and p also remains unchanged.

Compared with the random-walk strategy, HMC methods explore the target distribution much more efficiently due to the use of gradient information. HMC can travel a long distance in the phase space $(\theta ,p)$, which enhances the acceptance rate. However, it is really difficult for HMC to travel across the low-probability region, in which the value of the gradient of the potential energy is very large. Enlarging the momentum variable p may be helpful to jump over these regions, but the efficiency of HMC may decrease significantly.

Recently, some new developments of HMC have been proposed to make the HMC sampler more flexible. For example, Riemann manifold HMC [12] exploits the Riemann geometry to tune the mass M, which tends to create a more efficient HMC sampler. The “No U-Turn” sampler [7] can tune the step size $ϵ$, leapfrog length L and the simulation steps automatically. We note that, in principle, these state-of-the-art HMC samplers can also be combined with our proposed method.

2.2. Langevin Dynamics

Langevin dynamics were first used to describe the diffusion process of molecular systems, which was then integrated into MCMC samplers [26,27,28,29]. Langevin dynamics are a system of Ito-type stochastic differential equations, which take the following form:

$$\begin{array}{c}\hfill \mathrm{d}\theta =-M^{-1}\nabla U\left(\theta \right)\mathrm{d}t+\sqrt{\frac{2{\beta }^{-1}}{M}}\mathrm{d}W,\end{array}$$

(8)

where W represents the stochastic Wiener process, and M represents the diagonal mass matrix. $U\left(\theta \right)$ represents the energy function and ${\beta }^{-1}=k_{B}T$, where $k_B$ is Boltzmann constant and T represents the temperature. Since solving (8) is difficult, the Euler–Maruyama discretization [30] is used to approximately solve the differential equation, which takes the following form:

$$\begin{array}{c}\hfill {\theta }_{n+1}={\theta }_n-\frac{{ϵ}^2}{2}{M}^{-1}\nabla U\left({\theta }_n\right)+\sigma {ϵ}^2\sqrt{\frac{2{\beta }^{-1}}{M}}{z}_n,\end{array}$$

(9)

where $z\sim \mathcal{N}\left(z\right|0,\mathbf{I})$, and $ϵ$ represents the integration step size. Please note that (9) only gives the solution to the overdamped Langevin dynamics [31], which means that the friction term has not been concerned. The Langevin dynamics with the friction term [32] is defined as:

$$\begin{array}{cc}\hfill \mathrm{d}\theta & =M^{-1}p\mathrm{d}t\hfill \\ \hfill \mathrm{d}p& =-\nabla U\left(\theta \right)\mathrm{d}t-\gamma p\mathrm{d}t+\sigma \sqrt{M}\mathrm{d}W,\hfill \end{array}$$

(10)

where $\gamma >0$ represents the friction factor and $\sigma =\sqrt{2\gamma {\beta }^{-1}}$. However, simulating (10) is very difficult, and thus a discretization method [31] was used to solve the above stochastic differential equation, which takes the following form:

$$\begin{array}{cc}\hfill p_{(n+1)/2}& =p_n-\frac{ϵ}{2}\nabla U\left({\theta }_n\right)\hfill \\ \hfill {\theta }_{(n+1)/2}& ={\theta }_n+\frac{ϵ}{2}{M}^{-1}{p}_{(n+1)/2}\hfill \\ \hfill {\widehat{p}}_{(n+1)/2}& =a_1{p}_{(n+1)/2}+a_2\sqrt{M}{z}_{n+1}\hfill \\ \hfill {\theta }_{n+1}& ={\theta }_{(n+1)/2}+\frac{ϵ}{2}{M}^{-1}{\widehat{p}}_{(n+1)/2}\hfill \\ \hfill p_{n+1}& ={\widehat{p}}_{(n+1)/2}-\frac{ϵ}{2}\nabla U\left({\theta }_{n+1}\right),\hfill \end{array}$$

(11)

where $z\sim \mathcal{N}\left(z\right|0,\mathbf{I})$, $ϵ$ represents the step size and $a_1={\mathrm{e}}^{-\gamma ϵ}$ and $a_2=\sqrt{{\beta }^{-1}\left(1-a_1^2\right)}$.

Compared with Hamiltonian dynamics, Langevin dynamics allow the exploration of the state space more freely, for Langevin dynamics concern more random factors. Please note that in (9), Langevin dynamics added Brownian movement, which provides more flexibility for exploring the state space.

3. Langevin Hamiltonian Monte Carlo

HMC exploits the Hamiltonian dynamics to propose the new sample. However, the HMC sampler may have large autocorrelation in that each new sample is obtained from the deterministic calculation of the last sample. Specifically, Equation (3) defines the process of calculating a new state $\left({\theta }_{t+1},p_{t+1}\right)$ through the old state $\left({\theta }_t,p_t\right)$. HMC can reduce autocorrelation by increasing the value of the leapfrog size. However, it will make the sampler inefficient.

To further improve the performance of the autocorrelation of the HMC sampler, we propose the Langevin Hamiltonian Monte Carlo (LHMC). The main idea of LHMC is to take advantage of Langevin dynamics to add randomness to the proposed state. In Langevin dynamics, we consider that the total energy consists of the potential energy, kinetic energy and internal energy, which takes the following:

$$\begin{array}{c}\hfill H=U\left(\theta \right)+K\left(p\right)+Q,\end{array}$$

(12)

where Q represents the internal energy. The random thermal motion consumes the internal energy which finally transforms into the kinetic energy, which is described as:

$$\begin{array}{cc}\hfill p_t& =a_1{p}_{t-1}+a_2\sqrt{M}z\hfill \\ \hfill \Delta E& =K\left(p_t\right)-K\left(p_{t-1}\right)\hfill \\ \hfill Q_t& =Q_{t-1}-\Delta E.\hfill \end{array}$$

(13)

We use the Metropolis–Hasting criterion to accept the samples, so the acceptance rate $\phi$ takes the following form:

$$\begin{array}{cc}\hfill \phi & =min\left(1,\exp \left((U_{t-1}+K_{t-1}+Q_{t-1})-(U_t+K_t+Q_t)\right)\right)\hfill \\ \hfill & =min\left(1,\exp \left((U_{t-1}+K_{t-1}+Q_{t-1})-(U_t+K_t+Q_{t-1}-\Delta E)\right)\right)\hfill \\ \hfill & =min\left(1,\exp \left((U_{t-1}+K_{t-1})-(U_t+K_t-\Delta E)\right)\right).\hfill \end{array}$$

(14)

The detailed algorithms of LHMC are illustrated in Algorithm 2. Given the target distribution, LHMC exploits Langevin dynamics and Hamiltonian dynamics to explore state space via the discretization of LHMC, which can be summarized as three substages shown in Algorithm 3. The first substage is Langevin dynamics, which takes form as (15)–(17). The second substage is Hamiltonian dynamics, which takes form as (3), and the last substage is also Langevin dynamics. Assume the initial state is $(\theta ,p)$; then, a half update of the Langevin dynamics can be written as:

$$\begin{array}{cc}\hfill p_{(n+1)/2}& =p_n-\frac{ϵ}{2}\nabla U\left({\theta }_n\right)\hfill \\ \hfill {\theta }_{(n+1)/2}& ={\theta }_n+\frac{ϵ}{2}{M}^{-1}{p}_{(n+1)/2}.\hfill \end{array}$$

(15)

The random thermal motion of molecules takes the following form:

$$\begin{array}{cc}\hfill {\widehat{p}}_{(n+1)/2}& =a_1{p}_{(n+1)/2}+a_2\sqrt{M}{z}_{n+1}.\hfill \end{array}$$

(16)

The other half update of the Langevin dynamics can be written as:

$$\begin{array}{cc}\hfill {\theta }_{n+1}& ={\theta }_{(n+1)/2}+\frac{ϵ}{2}{M}^{-1}\hfill \\ \hfill p_{n+1}& ={\widehat{p}}_{(n+1)/2}-\frac{ϵ}{2}\nabla U\left({\theta }_{n+1}\right).\hfill \end{array}$$

(17)

Algorithm 2: Langevin Hamiltonian Monte Carlo (LHMC)

Input: Step size $ϵ$, leapfrog size L, starting point ${\theta }_{\left(1\right)}$, sample number N Output: Samples ${\theta }_{(1:N)}$ for $n=1$ to Ndo      $\widehat{\theta },\widehat{p},\Delta E=\mathrm{DLHMC}({\theta }_n,p,ϵ,L)$ //$\phantom{(}$      # see in Algorithm 3$\phantom{(}$      $u\sim \mathrm{U}niform(0,1)$      $\alpha =\min \left(1,\exp \left(U({\theta }_{(n-1)}+K\left(p_{(n-1)}\right))-(U(\widehat{\theta })+K\left(\widehat{p}\right)-\Delta E)\right)\right)$      if $\alpha >u$ then          $({\theta }_{(n+1)},p_{(n+1)})=(\widehat{\theta },\widehat{p})$     else          $({\theta }_{(n+1)},p_{(n+1)})=({\theta }_{\left(n\right)},p_{\left(n\right)})$     end if end for

Algorithm 3: Discretization for Langevin Hamiltonian Monte Carlo (DLHMC)

Input: Step size $ϵ$, leapfrog size L, starting point $({\theta }_{\left(n\right)},p_{\left(n\right)})$ Output: ${\theta }_{(n+1)}$, $p_{(n+1)}$, $\Delta E$ $\Delta E=0$ Obtaining ${\theta }_{(n+1)/6},p_{(n+1)/6}$ through (15). $\widehat{p}=p_{(n+1)/6}$ Obtaining the thermal motion of molecules $p_{(n+1)/6}$ through (16). $\Delta E=\Delta E+K\left(p_{(n+1)/6}\right)-K\left(\widehat{p}\right)$ Obtaining ${\theta }_{(n+1)/3},p_{(n+1)/3}$ through (17). Obtaining ${\theta }_{2(n+1)/3},p_{2(n+1)/3}$ by simulating Hamiltonian dynamics through (3). Obtaining ${\theta }_{5(n+1)/6},p_{5(n+1)/6}$ through (15). $\widehat{p}=p_{5(n+1)/6}$ Obtaining the thermal motion of molecules $p_{5(n+1)/6}$ through (16). $\Delta E=\Delta E+K(p_{5(n+1)/6)}-K\left(\widehat{p}\right)$ Obtaining ${\theta }_{(n+1)},p_{(n+1)}$ through (17).

We demonstrate the performance of LHMC on a strongly correlated Gaussian with variances $[{10}^2,{10}^{-2}]$ rotated by $\frac{\pi }{4}$, which is an extreme circumstance of Brooks et al. [3]. Each method in Figure 1 is run for 10,000 iterations with 1000 burn-in samples. We experimented 100 times and calculated the mean and variance of autocorrelation and maximum mean discrepancy. We set $M=1.2\mathbf{I}$ where $\mathbf{I}$ is an identity matrix with the same dimension as the sample distribution, $\gamma =0.5$, $T=200$, $ϵ=0.05$, leapfrog size $L=40$. As Figure 1 illustrates, LHMC achieves lower autocorrelation and faster convergence rate compared with HMC and Metropolis-adjusted Langevin algorithm (MALA) [29].

4. Variational Hybrid Monte Carlo

LHMC method, as a dynamics-based MCMC sampler, cannot sample from multi-modal distributions well when the modes are far away from each other, especially in high dimensions [15,16]. In this section, we present a novel MCMC method called variational hybrid Monte Carlo (VHMC) to tackle the above issue and prove that it can target the correct distribution.

4.1. Variational Hybrid Monte Carlo

Recent studies [15,16] have shown that MCMC samplers based on dynamics are challenging when dealing with multi-modal distributions, since there exist low-probability regions between the isolated modes. Once the initial point is chosen, these methods could only sample from one mode closer to the initial point.

In this study, we are aiming to address the problem of multi-modal sampling. We propose a new concept referred to as guide points which are illustrated in Figure 2. Guide points are samples generated from the variational distribution. With the help of these guide points, VHMC can travel across the low-probability regions between two modes. The specific construction of the variational distribution is described as follows.

Suppose we want to sample from the distribution $p\left(\theta \right)$. To obtain the local points with the highest probability density value, K samples are generated from the solution space, and we then use these K initial points with Adam [24] to calculate the solutions which take the following form:

$$\begin{array}{cc}\hfill g_t& ={\nabla }_{{\theta }_k}\pi \left({\theta }_{k(t-1)}\right)\hfill \\ \hfill m_t& ={\beta }_1·m_{t-1}+(1-{\beta }_1)·g_t\hfill \\ \hfill v_t& ={\beta }_2·v_{t-1}+(1-{\beta }_2)·g_t^2\hfill \\ \hfill {\widehat{m}}_t& =\frac{m_t}{{(1-{\beta }_1)}^t}\hfill \\ \hfill {\widehat{v}}_t& =\frac{v_t}{{(1-{\beta }_2)}^t}\hfill \\ \hfill {\theta }_{k\left(t\right)}& ={\theta }_{k(t-1)}-\alpha ·\frac{{\widehat{m}}_t}{\sqrt{{\widehat{v}}_t}+\varphi },\hfill \end{array}$$

(18)

where ${\theta }_{k\left(t\right)}$ represents the solution of the k-th mode in the t-th iteration, $\pi$ represents the target distribution, $\alpha$ represents the step size and ${\beta }_1,{\beta }_2$ represent the exponential decay rates for the moment estimates. From these local points with the highest probability density value, we can obtain K modes. For each mode, we use these local points with the highest probability density value as the initial state of LHMC sampler to generate S samples. After that, we can obtain the sample sets $mod{e}_k=\left\{\theta \right|\theta \sim p_k\left(\theta \right)\},k=1,\dots ,K$, which are generated in the single mode distribution $p_k\left(\theta \right),k=1,\dots ,K$. We assume that each mode follows Gaussian distribution $q_k\left(\theta \right),k=1,\dots ,K$. To use $q_k\left(\theta \right)$ to approximate $p_k\left(\theta \right)$, we use the $\mathbb{KL}$ divergence to quantify the similarity of the two distributions:

$$\begin{array}{cc}\hfill \mathbb{KL}\left(p_k\left(\theta \right)\right|\left|q_k\left(\theta \right)\right)& =\int p_k\left(\theta \right)\ln \frac{p_k\left(\theta \right)}{q_k(\theta ;\mu ,\Sigma )}\mathrm{d}\theta ={\mathbb{E}}_{p_k\left(\theta \right)}\left[\ln \frac{p_k\left(\theta \right)}{q_k(\theta ;\mu ,\Sigma )}\right],\hfill \end{array}$$

(19)

$$\begin{array}{c}\hfill \min \left(\mathbb{KL}\left(p_k\left(\theta \right)\right|\left|q_k\left(\theta \right)\right)\right)=\min \left[-\int p_k\left(\theta \right)\ln q_k(\theta ;\mu ,\Sigma )\mathrm{d}\theta \right]=\max \left[{\mathbb{E}}_{p_k\left(\theta \right)}\left(\ln q_k(\theta ;\mu ,\Sigma )\right)\right].\end{array}$$

(20)

Since minimizing the integrator in (19) is equivalent to maximizing the likelihood, we translate the problem into solving the likelihood of the Gaussian mixture distribution, which takes the following form:

$$\begin{array}{c}\hfill q\left(\theta \right)=\frac{1}{K}\sum _{k=1}^K\mathcal{N}({\mu }_k,{\Sigma }_k),\end{array}$$

(21)

where ${\mu }_k$ and ${\sigma }_k$ represent the mean and variance of k-th mode, respectively. The whole process of generating the variational distribution (GVD) is demonstrated in Algorithm 4.

Since constructing the variational distribution is to find the modes of the target distribution, this process can be managed parallelly, so the total cost contains two aspects, and they are the cost of Adam and the LHMC algorithm. The effectiveness of Adam makes us find the local optimal solution rapidly. In addition, because the initial samples are given through the Adam algorithm, LHMC can explore the space rapidly without any burn-in samples, which reduces the computational cost.

As shown in Figure 2, we can find that the mean of the variational distribution with guide points is similar to the target distribution, but the variance is different, which demonstrates that our variational method can move the initial random Gaussian to the center of the corresponding mode. If the variance of the initial random Gaussian is too large, the guide points will be scattered without good alignment with the target distribution, which means that the variational distribution loses strong guidance for the sampler trapped in the low-probability regions. Therefore, empirically, we set the random standard Gaussian to initialize the ADAM optimizer to keep a good alignment between guide points and the high probability centers of multiple modes.

Algorithm 4: Generating the variational distribution (GVD)

Input: Population size K, number of samples for each population S, the target distribution $\pi$, step size $\alpha$, exponential decay rates for the moment estimates ${\beta }_1,{\beta }_2$, convergence error $ϵ$, $\varphi ={10}^{-7}$, step size $\gamma$, leapfrog size L Output: The variational distribution q of $\pi$ ${\theta }_k\sim \mathcal{N}(0,\mathbf{I}),k\in [1,K]$ for $k=1$ to K do      $m_0,v_0,t=0$     while ${\theta }_{k\left(t\right)}-{\theta }_{k(t-1)}>ϵ$ do          $t=t+1$          $g_t={\nabla }_{{\theta }_k}\pi \left({\theta }_{k(t-1)}\right)$          $m_t={\beta }_1·m_{t-1}+(1-{\beta }_1)·g_t$          $v_t={\beta }_2·v_{t-1}+(1-{\beta }_2)·g_t^2$          ${\widehat{m}}_t=\frac{m_t}{{(1-{\beta }_1)}^t}$          ${\widehat{v}}_t=\frac{v_t}{{(1-{\beta }_2)}^t}$          ${\theta }_{k\left(t\right)}={\theta }_{k(t-1)}-\alpha ·\frac{{\widehat{m}}_t}{\sqrt{{\widehat{v}}_t}+\varphi }$     end while end for for $k=1$ to K do     $mod{e}_k=\mathrm{LHMC}(\gamma ,L,{\theta }_{k\left(t\right)},S)$//$\phantom{(}$      # see in Algorithm 2$\phantom{(}$     Calculating the mean ${\mu }_k$ and variance ${\sigma }_k$ for $mod{e}_k$. end for $q=\frac{1}{K}{\sum }_{k=1}^K\mathcal{N}({\mu }_k,{\sigma }_k)$

We generate new samples with two strategies. First, a vast number of samples are generated through the LHMC sampler with probability $\phi$. When LHMC rejects the sample with probability $1-\phi$, we then generate a sample with probability $\frac{\pi \left(x\right)}{cq\left(x\right)}$ from the variational distribution, whose derivation is demonstrated in Appendix A. We finally accept the samples generated from the variational distribution through the following equation:

$$\begin{array}{c}\hfill \min \left[1,\frac{r(x,x_r)}{1-\phi }\right],\end{array}$$

(22)

where $r(.)$ represents the rejection probability of the current state jumping into other states. The overall algorithm flow of VHMC is given in Algorithm 5 and the theoretical convergence analysis of VHMC is provided to demonstrate that VHMC could converge to the target distribution in Appendix B.

In the LHMC sampler, we found that in the Metropolis–Hasting procedure, the sampler rejects the proposed sample with probability $1-\phi$, where $\phi$ is defined in (14). This is an interesting phenomenon in that if the proposed sample is rejected, the traditional MCMC samplers would put the sample back to the last position. In our study, if one sample is rejected in the MH step, we generate a new sample from $q\left(\theta \right)$. We accept the samples generated from the variational distribution through the newly proposed acceptance rate, and thus the detailed balance holds.

Algorithm 5: Variational Hybrid Monte Carlo (VHMC)

Input: Step size $ϵ$, leapfrog size L, starting point ${\theta }_{\left(1\right)}$, sample number N, target distribution $\pi$, population size k, number of samples for each population s Output: Samples ${\theta }_{(1:N)}$ $q=\mathrm{GVD}(k,s,\pi )$//$\phantom{(}$      # see in Algorithm 4$\phantom{(}$ for $n=1$ to N do         $u\sim \mathrm{Uniform}(0,1)$         $\widehat{\theta },\widehat{p},\Delta E=\mathrm{DLHMC}({\theta }_n,p,ϵ,L)$ //$\phantom{(}$      # see in Algorithm 3$\phantom{(}$         $u\sim \mathrm{Uniform}(0,1)$         $\phi =\min \left(1,\exp \left(U({\theta }_{\left(n\right)}+K\left(p_{\left(n\right)}\right))-(U(\widehat{\theta })+K\left(\widehat{p}\right)-\Delta E)\right)\right)$         if $\phi >u$ then             ${\theta }_{(n+1)}=\widehat{\theta }$         else             ${\theta }^{*}\sim q\left(\theta \right)$             $h\sim \mathrm{Uniform}(0,1)$             while $h>\frac{\pi \left({\theta }^{*}\right)}{cq\left({\theta }^{*}\right)}$ do                ${\theta }^{*}\sim q\left(\theta \right)$             end while             $u\sim \mathrm{Uniform}(0,1)$             if $u<min(1,\frac{1-r({\theta }^{*},{\theta }_r^{*})}{1-\phi })$ then                ${\theta }_{(n+1)}={\theta }^{*}$             else                ${\theta }_{(n+1)}={\theta }_{\left(n\right)}$             end if      end if end for

4.2. Deficiency of Parallel HMC

Although we can run N HMC samplers in parallel to approximately sample from a multi-modal distribution, in high dimensions, this kind of method is inaccurate. In other words, the probability of each mode may be the same, which cannot reflect the actual distribution. Figure 3 shows that parallel HMC cannot sample from the actual distribution. The problem with parallel HMC is that the gradient direction cannot determine the probability of each mode.

5. Experiments

In this section, we investigate the performance of VHMC on multi-modal distributions and real datasets and compare our method with state-of-art algorithms. All our experiments are conducted on a standard computer with a 4.0 GHz Intel core i7 CPU. First, we introduce the performance index that will be used in the following parts.

Effective sample size. The variance of a Monte Carlo sampler is determined by its effective sample size (ESS) [3], which is defined as:

$$\mathrm{ESS}=N/(1+2\times \sum _{s=1}^{\infty }\rho \left(s\right)),$$

(23)

where N represents the number of all the samples and $\rho \left(s\right)$ represents the s—step autocorrelation where autocorrelation is an index that considers the correlation between two samples. Let X be a set of samples, and t be the number of iterations (t is an integer). Then, $X_t$ is the sample at time t of X. The autocorrelation between time s and t is defined as:

$$R(s,t)=\frac{E\left[(X_t-{\mu }_t)(X_s-{\mu }_s)\right]}{{\sigma }_t{\sigma }_s},$$

(24)

where E is the expected value operator. The correlation between two nearby samples can be measured with autocorrelation. The lower the value of autocorrelation, the more independent the samples.

Maximum mean discrepancy. The difference between samples drawn from two distributions can be measured as maximum mean discrepancy (MMD) [22], which takes the following form:

$$\begin{array}{cc}\hfill MM{D}^2[X,Y]=& \frac{1}{M^2}\sum _{i,j=1}^{M}k(x_i,x_j)-\frac{2}{MN}\sum _{i,j=1}^{M,N}k(x_i,y_j)+\frac{1}{N^2}\sum _{i,j=1}^{N}k(y_i,y_j),\hfill \end{array}$$

(25)

where M represents the sample number in X, N represents the sample number in Y and k represents the kernel function. By calculating the MMD value, we can analyze the convergence rate of the proposed methods.

Relative error of mean. This is a summary of the errors in approximating the expectation of variables across all dimensions [33], which is computed as:

$$RE{M}_t=\frac{{\sum }_{i=1}^d\mid \overline{{\theta }_i^t}-{\theta }_i^{*}\mid }{{\sum }_i\mid {\theta }_i^{*}\mid },$$

(26)

where $\overline{{\theta }_i^t}$ is the average of the i’th variable at time t, ${\theta }_i^{*}$ is the actual mean value, d denotes the dimension of sampling distribution and the denominator ${\sum }_i\mid {\theta }_i^{*}\mid$ represents the sum of $\mid {\theta }_i^{*}\mid$ on the true distribution.

5.1. Mixture of Isotropic Gaussians

We conduct our first experiment on two multi-modal distributions where we consider two simple 2D Gaussian mixtures whose distributions are analytically available. First, we consider a Gaussian mixture distribution whose modes are close to each other, and then we consider a Gaussian mixture whose modes are isolated and far away from each other. The distributions are given as follows: $p\left(\mathbf{x}\right)=\frac{1}{2}\mathcal{N}(\mathbf{x};\mu ,\Sigma )+\frac{1}{2}(\mathbf{x};-\mu ,\Sigma )$ for ${\sigma }_x^2={\sigma }_y^2=1$, ${\rho }_{xy}=0$, $\mathbf{x}=(x,y)\in {\mathbb{R}}^2$ and $\mu =(2.5,-2.5)$ (modes are close to each other) or $\mu =(6.5,-6.5)$ (modes are far away from each other). The experiment setting is the same with Tripuraneni et al. [16]. The purpose of the experiments is to sample points that are independently identically distributed in the multi-modal distribution correctly.

In this experiment, we compare MHMC, HMC, MGHMC [23] and parallel tempering (SAMCMC) [17] against VHMC. First, we compare the sample result of these methods intuitively. Then, averaged autocorrelation and MMD are used to compare the performance of each method further. Each method is run for 10,000 iterations with 1000 burn-in samples. The number of leapfrog steps is uniformly drawn from $(100-l,100+l)$ with $l=20$, which is suggested by Livingstone et al. [34]. We set step size as $ϵ=0.05$, friction coefficient as $\gamma =0.1$ and the initiate position as $\theta =(0,0)$. The authors in Tripuraneni et al. [16] indicated that the multi-modal problem is a challenge for HMC samplers. However, we find that HMC samplers can sample points from the multi-modal distribution, especially when the modes are close to each other.

Figure 4 clearly shows that when $\mu =(2.5,-2.5)$, three methods can sample the multi-modal distribution. Nevertheless, there is some difference between them. MHMC may sample from this mixture on a Gaussian distribution, but it updates state and position only according to random-walk proposals, which hardly jumps to the other sampling mode from a far place. The HMC sampler changes its sampling mode more frequently due to the gradient information of the target distribution. Obviously, combining the help of guide points, Brownian movement and gradient information VHMC changes its mode much more frequently than MHMC and HMC. From the result, we can also conclude that when the modes are close to each other, MHMC and HMC may sample this multi-modal distribution. However, when the modes are isolated and far away from each other with larger $\mu$, for instance, $\mu =(6.5,-6.5)$, both MHMC and HMC cannot sample from the target distribution, shown in the second row of Figure 4. This is because only random-walk proposals and gradient information cannot make them directly travel across the large low-probability regions. Nevertheless, our method still performs well by taking advantage of the samples generated from variational distribution to explore the phase space in different modes freely.

For the multiple modes far away from each other, HMC hardly changes its mode, so it converges to the target distribution slowly. On the contrary, VHMC changes the sampling mode very frequently, thus it converges to the target distribution quickly. To compare the convergence rate and the independence of the samples with state-of-the-art sampling methods, we exploit MMD and autocorrelation to describe the performance when sampling the Gaussian mixture.

MMD between exact samples generated from the target density and generated samples is used to describe the convergence performance of the samplers. We use a quadratic kernel [35] $k(x,x^{{}^{\prime }})={(1+\langle x,x^{{}^{\prime }}\rangle )}^2$ where $\langle ·\rangle$ denotes dot product, averaged over 100 runs of the Markov chains. Figure 5 demonstrates that our method achieves the best performance in convergence rate and autocorrelation after the burn-in period. Since our method converges to the target distribution quickly, we furthermore narrow the number of the first 500 samples. In general, comparing Figure 4 and Figure 5, we conclude that the convergence rate is inversely proportional to the MMD and autocorrelation for HMCM, HMC and VHMC.

To test the performance of the proposed method on high-dimensional multi-modal distribution, we conduct our experiments on 2 to 128 dimensions. The target distribution is given as:

$$\begin{array}{c}\hfill p\left(\theta \right)=\frac{1}{{\left(2\pi \right)}^{\frac{n}{2}}}\left(0.7·\exp \left[\frac{-{(x-{\mu }_0)}^{\top }(x-{\mu }_0)}{2}\right]+0.3·\exp \left[\frac{-{(x-{\mu }_1)}^{\top }(x-{\mu }_1)}{2}\right]\right),\end{array}$$

(27)

where ${\mu }_0=(a_1,\dots ,a_n),{\mu }_1=(b_1,\dots b_n),a_i=-1,b_i=1$ and n equals dimensions. Figure 6 shows that the proposed method has lower REM in high dimensions, which indicates that VHMC can sample from the high-dimensional distant multi-modal distributions.

5.2. Mixture of Heterogeneous Gaussians

In the first experiment, we have already discussed the Gaussian mixture when the variance of the modes is the same. In practice, real data distributions often have different variances and probability of modes. To demonstrate strong stability, we construct two new Gaussian mixtures with different variances and probability of modes. The first one is given as follows: $p\left(\theta \right)={\pi }_1\mathcal{N}(\theta ;{\mu }_1,{\sigma }_1)+{\pi }_2\mathcal{N}(\theta ;{\mu }_2,{\sigma }_2)+{\pi }_3\mathcal{N}(\theta ;{\mu }_3,{\sigma }_3).$ We set ${\pi }_1=0.1$, ${\pi }_2=0.8$, ${\pi }_3=0.1$, ${\sigma }_1=1$, ${\sigma }_2=3$, ${\sigma }_3=2$. The second one takes the following form: $p\left(\theta \right)=0.5\mathcal{N}(\theta ;{\mu }_1,{\sigma }_1)+0.5\mathcal{N}(\theta ;{\mu }_2,{\sigma }_2)$. Here, we set ${\sigma }_x^2=0.01$, ${\sigma }_y^2=1$, ${\rho }_{xy}=0.0$. Similar to the previous experiment, our method runs 10,000 iterations with 1000 burn-in samples. Figure 7 shows that VHMC has strong stability. Even when the variance becomes tiny, our method still shows advanced performance. From the second column of Figure 7, we can also observe that the HMC sampler may sample multi-modal distribution especially when the HMC sampler has chances to jump out of one mode. Although the distance between the left mode and the middle mode is the same as the distance between the middle mode and right mode, the different variances force the HMC sampler to sample from the left two modes.

5.3. Bayesian Logistic Regression

Logistic regression (LR) [36] is a traditional method for classification. We optimize the parameters by maximizing the logistic likelihood function. By exploiting the parameters, we can predict the class of the data.

To verify the performance on real datasets, we apply the proposed method to Bayesian logistic regression (BLR) [21] and our method is compared with logistic regression (LR), variational Bayesian logistic regression (VBLR) and HMC.

The likelihood function of a two-class classification problem can be defined as:

$$p\left(\mathrm{t}\right|w)=\prod _{n=1}^N{y}_n^{t_n}{[1-y_n]}^{1-t_n},$$

(28)

where $t_n\in \{0,1\}$ and $\mathrm{t}={(t_1,\dots ,t_N)}^{\top }$ and $y_n=p\left({\mathcal{C}}_1\right|{\varphi }_n)=\sigma \left(w^{\top }\varphi \right)$. $t_n$ represents the label of the data and $y_n$ represents the predicted value. We obtain the class of the data by means of integrating the logistic function on the posterior distribution.

We evaluate our methods on eight real-world datasets from the UCI repository [37]: Pima Indian (PI), Haberman (HA), Mammographic (MA), Blood (BL), Cryotherapy (CR), Immunotherapy (IM), Indian (IN), Diabetic (DI) using Bayesian logistic regression. The eight datasets are normalized to have zero mean value and unit variance. We give the Gaussian distribution $\mathcal{N}(0,100\mathbf{I})$ as the prior distribution of the parameters.

In each experiment, we run 10,000 iterations with 2000 burn-in samples. We draw leapfrog steps from a uniform distribution $\mathrm{Uniform}[80,120]$. We set step size $ϵ=0.00045$ and mass matrix $m=3\mathbf{I}$. We use the uniform distribution as the variational distribution. We run each dataset 100 times to calculate the mean and the standard deviation.

Results in terms of the accurate rate of prediction and area under the ROC curve (AUC) [38] are summarized in

Table 1. Classification accuracies for VBLR, LR, HMC and VHMC on eight data sets.

DATA	LR	VBLR	HMC	VHMC
HA	67.9 ± 0.6	67.7 ± 0.5	67.7 ± 0.4	68.2 ± 0.5
PI	82.5 ± 0.2	82.7 ± 0.2	82.4 ± 0.3	83.1 ± 0.2
MA	89.8 ± 0.1	89.9 ± 0.1	89.9 ± 0.2	89.9 ± 0.2
BL	75.1 ± 0.3	75.3 ± 0.3	71.2 ± 0.6	75.4 ± 0.3
CR	95.4 ± 0.1	95.6 ± 0.1	92.1 ± 0.2	95.7 ± 0.1
IM	77.5 ± 0.4	77.5 ± 0.4	77.4 ± 0.4	77.6 ± 0.4
IN	75.5 ± 0.2	75.8 ± 0.2	75.3 ± 0.3	75.9 ± 0.3
DI	82.6 ± 0.2	82.5 ± 0.2	81.8 ± 0.3	82.5 ± 0.2
RANK	2.75	2.125	3.5	1.25
PVALUE	0.012	0.041	0.056	/

and

Table 2. Area under the ROC curve (AUC) for VBLR, LR, HMC and VHMC on eight data sets.

DATA	LR	VBLR	HMC	VHMC
HA	73.8 ± 0.4	74.6 ± 0.4	74.6 ± 0.2	74.6 ± 0.4
PI	77.0 ± 0.2	77.1 ± 0.2	77.2 ± 0.3	77.6 ± 0.2
MA	82.6 ± 0.2	82.7 ± 0.1	82.5 ± 0.2	82.8 ± 0.1
BL	77.1 ± 0.2	77.1 ± 0.2	74.4 ± 0.5	77.2 ± 0.2
CR	86.9 ± 0.2	87.2 ± 0.2	84.2 ± 0.3	87.3 ± 0.2
IM	84.7 ± 0.2	84.8 ± 0.2	83.8 ± 0.2	84.9 ± 0.2
IN	72.7 ± 0.2	72.8 ± 0.2	71.3 ± 0.3	72.9 ± 0.2
DI	74.4 ± 0.2	74.6 ± 0.2	74.0 ± 0.2	74.7 ± 0.2
RANK	3.125	2.125	3.375	1.125
PVALUE	0.004	0.036	0.018	/

. The results show that in these eight datasets, VHMC achieves better performance in classification accuracy rate than in AUC. Compared with other methods, VHMC outperforms HMC and provides a similar performance to VBLR, which indicates that the method proposed in this paper can sample actual posterior distribution.

6. Discussion

The most innovative part of our method is that when LHMC refuses to sample, we use the sample generated by the variational distribution to replace the previous sample with a new MH step, thus accelerating the sampling from different modes with a higher acceptance rate. The convergence rate and autocorrelation of LHMC, VHMC and other typical MCMC sampling methods have been shown in Figure 1 and Figure 5. The effectiveness and stability of our method have been demonstrated on both isotropic and heterogeneous Gaussian mixture distributions in Figure 4 and Figure 7. We also compare VHMC with VBLR on the real datasets, where similar performance indicates that our VHMC can successfully sample actual posterior distribution. The object of our method is to approximate target distribution in the probability model, which may not be directly used in other research fields, such as molecular simulation with a more complex system. However, our work can still provide some inspiration for other researchers. It is likely to expand the exploration space and accelerate the target convergence by integrating the variational distribution samples into the MCMC simulation process. For the generation of variational mixture distribution, the traditional machine-learning method uses the EM (Expectation Maximization) algorithm with a slow convergence rate. Here, we adopt the gradient descent method combined with LHMC to estimate the parameters of the variational distribution and demonstrate its effectiveness and stability in different datasets. Of course, more methods are worth further exploration for the effectiveness and robustness of variational distribution in other fields.

7. Conclusions

In this study, we propose LHMC to further accelerate the convergence rate and reduce the autocorrelation of samples. We then present VHMC to exploit the information from the variational distribution of the target distribution to make effective distant multi-modal sampling available. A formal theoretical analysis is provided, which demonstrated that VHMC could converge to the target distribution. Our findings are supported by synthetic and real data experiments, which showed that VHMC brings multiple benefits, such as providing superior performance in multi-modal sampling and lower autocorrelation. In the future, we plan to introduce equipotential transformation to the LHMC sampler. If an equipotential sample can be obtained, the autocorrelation of the samples can be significantly reduced.