Numerical stability of force-gradient integrators and their Hessian-free variants in lattice QCD simulations

TL;DR

The paper conducts a linear stability analysis of force-gradient integrators and their Hessian-free variants in lattice QCD, using the harmonic oscillator as a testbed to derive stability thresholds and compare performance. It identifies minimum-error, self-adjoint integrators with favorable stability-accuracy trade-offs, and introduces the relative stability threshold $\mathrm{Eff}_{\mathtt{stab}}$ to contextualize stability against computational cost. Numerical tests in the Schwinger model and lattice QCD with heavy Wilson fermions confirm that Hessian-free force-gradient integrators can achieve comparable or better acceptance with reduced computational effort, while a sixth-order regime favors conventional splitting. The findings provide actionable guidance for selecting integrators in HMC simulations and motivate further development of problem-weighted efficiency measures and nested multilevel schemes.

Abstract

A comprehensive linear stability analysis of force-gradient integrators and their Hessian-free variants is carried out by investigating the harmonic oscillator as a test equation. The analysis reveals that the linear stability of conventional force-gradient integrators and their Hessian-free counterparts coincides. By performing detailed linear stability investigations for the entire family of self-adjoint integrators with up to eleven exponentials per time step, we detect promising integrator variants that are providing a good trade-off between accuracy and numerical stability. Special attention is given to the application of these promising integrator variants within the Hamiltonian Monte Carlo algorithm, particularly in the context of interacting field theories. Simulations for the two-dimensional Schwinger model are conducted to demonstrate that there are no significant differences in the stability domain of a force-gradient integrator and its Hessian-free counterpart. Lattice QCD simulations with two heavy Wilson fermions emphasize that Hessian-free force-gradient integrators with a larger stability threshold allow for a more efficient computational process compared to conventional splitting methods. Furthermore, detailed investigations of the stability threshold are performed by investigating Nf = 2 twisted-mass fermions and nested integrators, highlighting the reliability of the linear stability threshold for lattice QCD simulations.

Figures (7)

• Figure 1: The efficiency measures $\mathrm{Eff}_{\mathtt{err}}^{(4)}$ and $\mathrm{Eff}_{\mathtt{stab}}$ for the decomposition algorithms \ref{['eq:ABAXABA']} evaluated at different values of $b_1 \in (0,2]$.
• Figure 2: The efficiency measure $\mathrm{Eff}_{\mathtt{err}}^{(4)}$ for the decomposition algorithms \ref{['eq:BAXABAXAB']} evaluated at different values of $a_2 \in [0.15,0.3]$ and $b_2 \in [0.25,0.5]$. Moreover, the white contour lines indicate the values for the relative stability threshold $\mathrm{Eff}_{\mathtt{stab}}$.
• Figure 3: Visualization of $\mathrm{Eff}_{\mathtt{err}}^{(2)}(\lambda)$\ref{['eq:Eff_err_OMF2']} and $\mathrm{Eff}_{\mathtt{stab}}(\lambda)$\ref{['eq:Eff_stab_OMF2']} for the splitting methods BABAB \ref{['eq:BABAB']} and ABABA \ref{['eq:ABABA']}.
• Figure 4: Two-dimensional Schwinger model. Variance of $\Delta\mathcal{H}$ vs. the number of time steps per trajectory $N$ for the conventional force-gradient integrator BACAB \ref{['eq:BACAB']} ($\circ$) and its Hessian-free counterpart BADAB \ref{['eq:BADAB']} ($\triangle$). For any $N$, $1000$ trajectories of length $\tau = 4.0$ have been computed.
• Figure 5: Variance of $\Delta\mathcal{H}$ vs. the number of force evaluations per trajectory $(n_f + n_g)\cdot N$ for all non-dominated variants of self-adjoint Hessian-free force-gradient integrators (blue lines) and conventional splitting methods (red lines). For any integrator and step size, 100 trajectories have been computed.