Hybrid Monte Carlo algorithm for sampling rare events in space-time histories of stochastic fields

We introduce a variant of the Hybrid Monte Carlo (HMC) algorithm to address large-deviation statistics in stochastic hydrodynamics. Based on the path-integral approach to stochastic (partial) differential equations, our HMC algorithm samples space-time histories of the dynamical degrees of freedom under the influence of random noise. First, we validate and benchmark the HMC algorithm by reproducing multiscale properties of the one-dimensional Burgers equation driven by Gaussian and white-in-time noise. Second, we show how to implement an importance sampling protocol to significantly enhance, by orders of magnitudes, the probability to sample extreme and rare events, making it possible to estimate moments of field variables of extremely high order (up to 30 and more). By employing reweighting techniques, we map the biased configurations back to the original probability measure in order to probe their statistical importance. Finally, we show that by biasing the system towards very intense negative gradients, the HMC algorithm is able to explore the statistical fluctuations around instanton configurations. Our results will also be interesting and relevant in lattice gauge theory since they provide unique insights into reweighting techniques.

Figure 1

(a) Shock wave as a particular realization of the Burgers equation. Its two main features are the ramps and jumps of finite width. (b) PDF of velocity gradients of the Burgers equation for two different Reynolds numbers compared with a Gaussian distribution. The finite jumps in the shock wave profile are responsible for the heavy left tail, while the right tail is due to the ramps.

Figure 2

Schematic description of a single HMC iteration of trajectory length $\tau$. The dashes in the top arrow represent the $N_{\mathrm{\Delta }\tau }$ number of intermediate steps of size $\mathrm{\Delta }\tau$, so $\tau =N_{\mathrm{\Delta }\tau }\mathrm{\Delta }\tau$ holds. After the numerical integration of Eqs. (19), a configuration $v_{\tau }$ is proposed and is accepted with probability $p=min(1,e^{-\mathrm{\Delta }\mathcal{H}})$. If accepted, we use $v_{\tau }$ as the initial field for the next iteration; otherwise we restart from $v_0$. In both cases, we discard ${\pi }_{\tau }$ and resample the momenta according to the Gaussian distribution for the next HMC iteration. (Note that the conjugate momentum fields are not plotted.)

Figure 3

Kinetic energy as a function of time for three samples corresponding to three different viscosities for the (a) HMC and (b) DNS.

Figure 4

(a) Temporal evolution of the ensemble-average mean kinetic energy $\langle {\overline{\varepsilon }}_{\mathrm{kin}}\left(t\right)\rangle$. (b) Temporal evolution of the ensemble-average mean energy dissipation $\langle {\overline{\varepsilon }}_{\mathrm{diss}}\left(t\right)\rangle$. Note that $\langle {\overline{\varepsilon }}_{\mathrm{diss}}(t>t_s)\rangle =1$ because the energy injection is fixed at ${lim}_{k\to \infty }\langle {\overline{\varepsilon }}_{\mathrm{in}}\left(k\right)\rangle =1$. Closed gray symbols correspond to the DNS results (lines denote their interpolation) and open colored symbols correspond to the HMC results.

Figure 5

(a) Ensemble-average energy spectra $\langle \overline{E}\left(k\right)\rangle$. Closed gray symbols correspond to the DNS results (lines denote their interpolation) and open colored symbols correspond to the HMC results. (b) Probability distribution functions of velocity gradients. The gray dashed line corresponds to the DNS results, while open colored symbols correspond to the HMC results.

Figure 6

Results for the HMC using periodic boundary conditions in time. (a) Temporal evolution of the ensemble-average mean kinetic energy $\langle {\overline{\varepsilon }}_{\mathrm{kin}}\left(t\right)\rangle$. (b) Temporal evolution of the ensemble-average mean energy dissipation $\langle {\overline{\varepsilon }}_{\mathrm{diss}}\left(t\right)\rangle$. Closed gray symbols denote the DNS results (lines correspond to their interpolation) and open colored symbols denote the HMC results.

Figure 7

Sample configurations from the HMC simulation using the constraint $\mathrm{\Delta }{\mathcal{S}}_1$ with $c_1=1.9$, enforcing a negative velocity gradient maximization at $(x=0,t=t_f)$. The color spectrum corresponds to the intensity of the velocity field.

Figure 8

(a) Ensemble average of the velocity field $v(t=t_f)$ using the HMC with the action ${\mathcal{S}}^{\prime }=\mathcal{S}+\mathrm{\Delta }{\mathcal{S}}_1$ for different values of $c_1$. (b) PDF of velocity gradients (DNS versus the HMC). For the HMC we measure $P^{\prime }\left(w\right)$ only at the space-time point where we constrain the ensemble, i.e., at $(x=0,t=t_f)$.

Figure 9

PDF of velocity gradients (DNS versus the HMC), i.e., $P^{\prime }\left(w\right)$, generated using the action ${\mathcal{S}}^{\prime }$ and measured only at the point that we constrain, i.e., at $(x=0,t=t_f)$, for different values of $c_i$ and $w_i$, $i=2,3$, using (a) $\mathrm{\Delta }{\mathcal{S}}_2$ and (b) $\mathrm{\Delta }{\mathcal{S}}_3$.

Figure 10

Ensemble-average mean kinetic energy of the HMC vs DNS using $\mathrm{\Delta }{\mathcal{S}}_1$ for different $c_1$ (a) before reweighting and (b) after reweighting. Notice that the error bars for the $c_1=1.9$ case after reweighting are pronounced as for this choice the fluctuations introduced by the reweighting factor become significantly large. (c) PDF of the nonreweighted mean kinetic energy ${\overline{\varepsilon }}_{\mathrm{kin}}\left(t\right)$ measured only for the last time slice $t=t_f$, in the case of the HMC, and using $\mathrm{\Delta }{\mathcal{S}}_1$ for different $c_1$.

Figure 11

PDF of the velocity gradients for the HMC (blue and red symbols) versus DNS (black line) (a) without rescaling and (b) rescaling the reweighted HMC histogram by dividing it by $\kappa =1.93$. For the HMC we use the constraint $\mathrm{\Delta }{\mathcal{S}}_1$ with $c_1=1.9$, and $P^{\prime }\left(w\right)$ is measured considering only the site $(x=0,t=t_f)$ on which the constraint is enforced. The inset shows the chosen interval $[-12,-6]$.

Figure 12

PDF of velocity gradients for the HMC (blue and red symbols) against DNS (black line) (a) without rescaling and (b) rescaling the reweighted HMC histogram by dividing it by $\kappa =2.63\times {10}^5$. We consider only the extracted histogram from the lattice point on which the constraint $\mathrm{\Delta }{\mathcal{S}}_2$ acted, i.e., $(x=0,t=t_f)$, in the case of the HMC, with $c_2=80$ and $w_2=18$. The inset shows the chosen interval $[-16,-11]$ for the rescaling.

Figure 13

PDF of velocity gradients $P\left(w\right)$ for the HMC and DNS. We consider here only the extracted histogram from the lattice point on which the constraint $\mathrm{\Delta }\mathcal{S}$ acted, i.e., $(x=0,t=t_f)$, in the case of the HMC. We show the effect of reweighting for different parameters of the constraints: (a) $\mathrm{\Delta }{\mathcal{S}}_1$, (b) $\mathrm{\Delta }{\mathcal{S}}_2$, (c) $\mathrm{\Delta }{\mathcal{S}}_3$, and (d) relative bin error $\delta P\left(w\right)/P\left(w\right)$, with the error being evaluated using Eq. (30). Regarding rescaling, for those $P\left(w\right)$ of which the overlap with the DNS was marginal or nonexistent, the rescaled $P\left(w\right)$ for $c_1=1.9$ was used.

Figure 14

(a)–(c) Velocity gradient PDF multiplied by a moment $w^q$. The $P\left(w\right)$ of the HMC is reweighted, rescaled with $\kappa$, and we consider only the lattice point at which the constraint acts. Here we used $\mathrm{\Delta }{\mathcal{S}}_1$ for $c_1=1.9$. (d)–(f) Computational time to the stabilized running average of the velocity gradient moment $\langle {\left({\partial }_{x}v\right)}^q\rangle$, divided with respect to the final stabilized value. Regarding DNS, any site belonging to the stationary regime is considered.

Figure 15

(a) Ensemble average of velocity configurations generated by the HMC using $\mathrm{\Delta }{\mathcal{S}}_1$ with $c_1=1.9$ compared to the classical instanton velocity-field profile generated for $\lambda =-1.148$ and $w=-24.23$. (b) PDF of the velocity gradients for the classical instanton (for a range of values of $\lambda$ and $w$), HMC simulation, and DNS.

Figure 16

(a) Ensemble average of velocity configurations generated by the HMC using $\mathrm{\Delta }{\mathcal{S}}_1$ with $c_1=1.9$. (b) Instanton velocity field for $\lambda =-1.148$ and $w=-24.23$. By averaging over the HMC velocity-field ensemble the spatiotemporal shape of the classical instanton is restored.